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

      1: #ifndef included_ALE_Generator_hh
      2: #define included_ALE_Generator_hh

      4: #ifndef  included_ALE_Distribution_hh
      5: #include <sieve/Distribution.hh>
      6: #endif

      8: #ifdef PETSC_HAVE_TRIANGLE
      9: #include <triangle.h>
     10: #endif
     11: #ifdef PETSC_HAVE_TETGEN
     12: #include <tetgen.h>
     13: #endif

     15: namespace ALE {
     16: #ifdef PETSC_HAVE_TRIANGLE
     17:   namespace Triangle {
     18:     template<typename Mesh>
     19:     class Generator {
     20:       class SegmentVisitor {
     21:       protected:
     22:         const int dim;
     23:         int *segmentlist;
     24:         typename Mesh::numbering_type& vNumbering;
     25:         int idx, v;
     26:       public:
     27:         SegmentVisitor(const int dim, int segmentlist[], typename Mesh::numbering_type& vNumbering) : dim(dim), segmentlist(segmentlist), vNumbering(vNumbering), idx(0), v(0) {};
     28:         ~SegmentVisitor() {};
     29:       public:
     30:         template<typename Point>
     31:         void visitPoint(const Point& point) {
     32:           this->segmentlist[this->idx*dim + (this->v++)] = this->vNumbering.getIndex(point);
     33:         }
     34:         template<typename Arrow>
     35:         void visitArrow(const Arrow& arrow) {}
     36:         void setIndex(const int idx) {this->idx = idx; this->v = 0;};
     37:       };
     38:     public:
     39:       static void initInput(struct triangulateio *inputCtx) {
     40:         inputCtx->numberofpoints = 0;
     41:         inputCtx->numberofpointattributes = 0;
     42:         inputCtx->pointlist = NULL;
     43:         inputCtx->pointattributelist = NULL;
     44:         inputCtx->pointmarkerlist = NULL;
     45:         inputCtx->numberofsegments = 0;
     46:         inputCtx->segmentlist = NULL;
     47:         inputCtx->segmentmarkerlist = NULL;
     48:         inputCtx->numberoftriangleattributes = 0;
     49:         inputCtx->trianglelist = NULL;
     50:         inputCtx->numberofholes = 0;
     51:         inputCtx->holelist = NULL;
     52:         inputCtx->numberofregions = 0;
     53:         inputCtx->regionlist = NULL;
     54:       };
     55:       static void initOutput(struct triangulateio *outputCtx) {
     56:         outputCtx->numberofpoints = 0;
     57:         outputCtx->pointlist = NULL;
     58:         outputCtx->pointattributelist = NULL;
     59:         outputCtx->pointmarkerlist = NULL;
     60:         outputCtx->numberoftriangles = 0;
     61:         outputCtx->trianglelist = NULL;
     62:         outputCtx->triangleattributelist = NULL;
     63:         outputCtx->neighborlist = NULL;
     64:         outputCtx->segmentlist = NULL;
     65:         outputCtx->segmentmarkerlist = NULL;
     66:         outputCtx->edgelist = NULL;
     67:         outputCtx->edgemarkerlist = NULL;
     68:       };
     69:       static void finiOutput(struct triangulateio *outputCtx) {
     70:         free(outputCtx->pointmarkerlist);
     71:         free(outputCtx->edgelist);
     72:         free(outputCtx->edgemarkerlist);
     73:         free(outputCtx->trianglelist);
     74:         free(outputCtx->neighborlist);
     75:       };
     78:       static Obj<Mesh> generateMesh(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false) {
     79:         int                          dim   = 2;
     80:         Obj<Mesh>                    mesh  = new Mesh(boundary->comm(), dim, boundary->debug());
     81:         const Obj<typename Mesh::sieve_type>& sieve = boundary->getSieve();
     82:         const bool                   createConvexHull = false;
     83:         struct triangulateio in;
     84:         struct triangulateio out;

     87:         initInput(&in);
     88:         initOutput(&out);
     89:         const Obj<typename Mesh::label_sequence>&    vertices    = boundary->depthStratum(0);
     90:         const Obj<typename Mesh::label_type>&        markers     = boundary->getLabel("marker");
     91:         const Obj<typename Mesh::real_section_type>& coordinates = boundary->getRealSection("coordinates");
     92:         const Obj<typename Mesh::numbering_type>&    vNumbering  = boundary->getFactory()->getLocalNumbering(boundary, 0);

     94:         in.numberofpoints = vertices->size();
     95:         if (in.numberofpoints > 0) {
     96:           const typename Mesh::label_sequence::iterator vEnd = vertices->end();

     98:           PetscMalloc(in.numberofpoints * dim * sizeof(double), &in.pointlist);CHKERRXX(ierr);
     99:           PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);CHKERRXX(ierr);
    100:           for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
    101:             const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
    102:             const int                                  idx   = vNumbering->getIndex(*v_iter);

    104:             for(int d = 0; d < dim; d++) {
    105:               in.pointlist[idx*dim + d] = array[d];
    106:             }
    107:             in.pointmarkerlist[idx] = boundary->getValue(markers, *v_iter);
    108:           }
    109:         }
    110:         const Obj<typename Mesh::label_sequence>& edges      = boundary->depthStratum(1);
    111:         const Obj<typename Mesh::numbering_type>& eNumbering = boundary->getFactory()->getLocalNumbering(boundary, 1);

    113:         in.numberofsegments = edges->size();
    114:         if (in.numberofsegments > 0) {
    115:           const typename Mesh::label_sequence::iterator eEnd = edges->end();

    117:           PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);CHKERRXX(ierr);
    118:           PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);CHKERRXX(ierr);
    119:           for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != eEnd; ++e_iter) {
    120:             const Obj<typename Mesh::sieve_type::traits::coneSequence>&     cone = sieve->cone(*e_iter);
    121:             const typename Mesh::sieve_type::traits::coneSequence::iterator cEnd = cone->end();
    122:             const int                                                       idx  = eNumbering->getIndex(*e_iter);
    123:             int                                                             v    = 0;

    125:             for(typename Mesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
    126:               in.segmentlist[idx*dim + (v++)] = vNumbering->getIndex(*c_iter);
    127:             }
    128:             in.segmentmarkerlist[idx] = boundary->getValue(markers, *e_iter);
    129:           }
    130:         }
    131:         const typename Mesh::holes_type& holes = boundary->getHoles();

    133:         in.numberofholes = holes.size();
    134:         if (in.numberofholes > 0) {
    135:           PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);CHKERRXX(ierr);
    136:           for(int h = 0; h < in.numberofholes; ++h) {
    137:             for(int d = 0; d < dim; ++d) {
    138:               in.holelist[h*dim+d] = holes[h][d];
    139:             }
    140:           }
    141:         }
    142:         if (mesh->commRank() == 0) {
    143:           std::string args("pqezQ");

    145:           if (createConvexHull) {
    146:             args += "c";
    147:           }
    148:           if (constrained) {
    149:             args = "zepDQ";
    150:           }
    151:           triangulate((char *) args.c_str(), &in, &out, NULL);
    152:         }

    154:         if (in.pointlist)         {PetscFree(in.pointlist);CHKERRXX(ierr);}
    155:         if (in.pointmarkerlist)   {PetscFree(in.pointmarkerlist);CHKERRXX(ierr);}
    156:         if (in.segmentlist)       {PetscFree(in.segmentlist);CHKERRXX(ierr);}
    157:         if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);CHKERRXX(ierr);}
    158:         if (in.holelist)          {PetscFree(in.holelist);CHKERRXX(ierr);}
    159:         const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
    160:         int     numCorners  = 3;
    161:         int     numCells    = out.numberoftriangles;
    162:         int    *cells       = out.trianglelist;
    163:         int     numVertices = out.numberofpoints;
    164:         double *coords      = out.pointlist;

    166:         ALE::SieveBuilder<Mesh>::buildTopology(newSieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, mesh->getArrowSection("orientation"));
    167:         mesh->setSieve(newSieve);
    168:         mesh->stratify();
    169:         ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, dim, coords);
    170:         const Obj<typename Mesh::label_type>& newMarkers = mesh->createLabel("marker");

    172:         if (mesh->commRank() == 0) {
    173:           for(int v = 0; v < out.numberofpoints; v++) {
    174:             if (out.pointmarkerlist[v]) {
    175:               mesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
    176:             }
    177:           }
    178:           if (interpolate) {
    179:             for(int e = 0; e < out.numberofedges; e++) {
    180:               if (out.edgemarkerlist[e]) {
    181:                 const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
    182:                 const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
    183:                 const Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(vertexA, vertexB, 1);

    185:                 mesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
    186:               }
    187:             }
    188:           }
    189:         }
    190:         mesh->copyHoles(boundary);
    191:         finiOutput(&out);
    192:         return mesh;
    193:       };
    196:       static Obj<Mesh> generateMeshV(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false, const bool renumber = false) {
    197:         typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
    198:         typedef typename Mesh::real_section_type::value_type real;
    199:         int                                   dim   = 2;
    200:         const Obj<Mesh>                       mesh  = new Mesh(boundary->comm(), dim, boundary->debug());
    201:         const Obj<typename Mesh::sieve_type>& sieve = boundary->getSieve();
    202:         const bool                            createConvexHull = false;
    203:         struct triangulateio in;
    204:         struct triangulateio out;
    205:         PetscErrorCode       ierr;

    207:         initInput(&in);
    208:         initOutput(&out);
    209:         const Obj<typename Mesh::label_sequence>&    vertices    = boundary->depthStratum(0);
    210:         const Obj<typename Mesh::label_type>&        markers     = boundary->getLabel("marker");
    211:         const Obj<typename Mesh::real_section_type>& coordinates = boundary->getRealSection("coordinates");
    212:         const Obj<typename Mesh::numbering_type>&    vNumbering  = boundary->getFactory()->getLocalNumbering(boundary, 0);

    214:         in.numberofpoints = vertices->size();
    215:         if (in.numberofpoints > 0) {
    216:           const typename Mesh::label_sequence::iterator vEnd = vertices->end();

    218:           PetscMalloc(in.numberofpoints * dim * sizeof(double), &in.pointlist);CHKERRXX(ierr);
    219:           PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);CHKERRXX(ierr);
    220:           for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
    221:             const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
    222:             const int                                           idx   = vNumbering->getIndex(*v_iter);

    224:             for(int d = 0; d < dim; ++d) {
    225:               in.pointlist[idx*dim + d] = array[d];
    226:             }
    227:             in.pointmarkerlist[idx] = boundary->getValue(markers, *v_iter);
    228:           }
    229:         }
    230:         const Obj<typename Mesh::label_sequence>& edges      = boundary->depthStratum(1);
    231:         const Obj<typename Mesh::numbering_type>& eNumbering = boundary->getFactory()->getLocalNumbering(boundary, 1);

    233:         in.numberofsegments = edges->size();
    234:         if (in.numberofsegments > 0) {
    235:           const typename Mesh::label_sequence::iterator eEnd = edges->end();

    237:           PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);CHKERRXX(ierr);
    238:           PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);CHKERRXX(ierr);
    239:           SegmentVisitor sV(dim, in.segmentlist, *vNumbering);
    240:           for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != eEnd; ++e_iter) {
    241:             const int idx = eNumbering->getIndex(*e_iter);

    243:             sV.setIndex(idx);
    244:             sieve->cone(*e_iter, sV);
    245:             in.segmentmarkerlist[idx] = boundary->getValue(markers, *e_iter);
    246:           }
    247:         }
    248:         const typename Mesh::holes_type& holes = boundary->getHoles();

    250:         in.numberofholes = holes.size();
    251:         if (in.numberofholes > 0) {
    252:           PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);CHKERRXX(ierr);
    253:           for(int h = 0; h < in.numberofholes; ++h) {
    254:             for(int d = 0; d < dim; ++d) {
    255:               in.holelist[h*dim+d] = holes[h][d];
    256:             }
    257:           }
    258:         }
    259:         if (mesh->commRank() == 0) {
    260:           std::string args("pqezQ");

    262:           if (createConvexHull) {
    263:             args += "c";
    264:           }
    265:           if (constrained) {
    266:             args = "zepDQ";
    267:           }
    268:           triangulate((char *) args.c_str(), &in, &out, NULL);
    269:         }
    270:         if (in.pointlist)         {PetscFree(in.pointlist);CHKERRXX(ierr);}
    271:         if (in.pointmarkerlist)   {PetscFree(in.pointmarkerlist);CHKERRXX(ierr);}
    272:         if (in.segmentlist)       {PetscFree(in.segmentlist);CHKERRXX(ierr);}
    273:         if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);CHKERRXX(ierr);}
    274:         if (in.holelist)          {PetscFree(in.holelist);CHKERRXX(ierr);}
    275:         const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
    276:         const Obj<FlexMesh>                  m        = new FlexMesh(boundary->comm(), dim, boundary->debug());
    277:         const Obj<FlexMesh::sieve_type>      newS     = new FlexMesh::sieve_type(m->comm(), m->debug());
    278:         int     numCorners  = 3;
    279:         int     numCells    = out.numberoftriangles;
    280:         int    *cells       = out.trianglelist;
    281:         int     numVertices = out.numberofpoints;
    282:         double *coords      = out.pointlist;
    283:         real   *coordsR;

    285:         ALE::SieveBuilder<FlexMesh>::buildTopology(newS, dim, numCells, cells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
    286:         m->setSieve(newS);
    287:         m->stratify();
    288:         mesh->setSieve(newSieve);
    289:         std::map<typename Mesh::point_type,typename Mesh::point_type> renumbering;
    290:         ALE::ISieveConverter::convertSieve(*newS, *newSieve, renumbering, renumber);
    291:         mesh->stratify();
    292:         ALE::ISieveConverter::convertOrientation(*newS, *newSieve, renumbering,
    293:         m->getArrowSection("orientation").ptr());
    294:         {
    295:           if (sizeof(double) == sizeof(real)) {
    296:             coordsR = (real *) coords;
    297:           } else {
    298:             coordsR = new real[numVertices*dim];
    299:             for(int i = 0; i < numVertices*dim; ++i) coordsR[i] = coords[i];
    300:           }
    301:         }
    302:         ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, dim, coordsR);
    303:         {
    304:           if (sizeof(double) != sizeof(real)) {
    305:             delete [] coordsR;
    306:           }
    307:         }
    308:         const Obj<typename Mesh::label_type>& newMarkers = mesh->createLabel("marker");

    310:         if (mesh->commRank() == 0) {
    311: #ifdef IMESH_NEW_LABELS
    312:           int size = 0;

    314:           newMarkers->setChart(mesh->getSieve()->getChart());
    315:           for(int v = 0; v < out.numberofpoints; v++) {
    316:             if (out.pointmarkerlist[v]) {
    317:               newMarkers->setConeSize(v+out.numberoftriangles, 1);
    318:               size++;
    319:             }
    320:           }
    321:           if (interpolate) {
    322:             for(int e = 0; e < out.numberofedges; e++) {
    323:               if (out.edgemarkerlist[e]) {
    324:                 const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
    325:                 const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
    326:                 const Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(vertexA, vertexB, 1);

    328:                 newMarkers->setConeSize(*edge->begin(), 1);
    329:                 size++;
    330:               }
    331:             }
    332:           }
    333:           newMarkers->setSupportSize(0, size);
    334:           newMarkers->allocate();
    335: #endif
    336:           for(int v = 0; v < out.numberofpoints; v++) {
    337:             if (out.pointmarkerlist[v]) {
    338:               if (renumber) {
    339:                 mesh->setValue(newMarkers, renumbering[v+out.numberoftriangles], out.pointmarkerlist[v]);
    340:               } else {
    341:                 mesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
    342:               }
    343:             }
    344:           }
    345:           if (interpolate) {
    346:             for(int e = 0; e < out.numberofedges; e++) {
    347:               if (out.edgemarkerlist[e]) {
    348:                 const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
    349:                 const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
    350:                 const Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(vertexA, vertexB, 1);

    352:                 mesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
    353:               }
    354:             }
    355:           }
    356: #ifdef IMESH_NEW_LABELS
    357:           newMarkers->recalculateLabel();
    358: #endif
    359:         }
    360:         mesh->copyHoles(boundary);
    361:         finiOutput(&out);
    362:         return mesh;
    363:       };
    364:     };
    365:     template<typename Mesh>
    366:     class Refiner {
    367:     public:
    368:       static Obj<Mesh> refineMesh(const Obj<Mesh>& serialMesh, const double maxVolumes[], const bool interpolate = false, const bool forceSerial = false) {
    369:         const int                    dim         = serialMesh->getDimension();
    370:         const Obj<Mesh>              refMesh     = new Mesh(serialMesh->comm(), dim, serialMesh->debug());
    371:         const Obj<typename Mesh::sieve_type>& serialSieve = serialMesh->getSieve();
    372:         struct triangulateio in;
    373:         struct triangulateio out;
    374:         PetscErrorCode       ierr;

    376:         Generator<Mesh>::initInput(&in);
    377:         Generator<Mesh>::initOutput(&out);
    378:         const Obj<typename Mesh::label_sequence>&    vertices    = serialMesh->depthStratum(0);
    379:         const Obj<typename Mesh::label_type>&        markers     = serialMesh->getLabel("marker");
    380:         const Obj<typename Mesh::real_section_type>& coordinates = serialMesh->getRealSection("coordinates");
    381:         const Obj<typename Mesh::numbering_type>&    vNumbering  = serialMesh->getFactory()->getLocalNumbering(serialMesh, 0);

    383:         in.numberofpoints = vertices->size();
    384:         if (in.numberofpoints > 0) {
    385:           const typename Mesh::label_sequence::iterator vEnd = vertices->end();

    387:           PetscMalloc(in.numberofpoints * dim * sizeof(double), &in.pointlist);CHKERRXX(ierr);
    388:           PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);CHKERRXX(ierr);
    389:           for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
    390:             const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
    391:             const int                                  idx   = vNumbering->getIndex(*v_iter);

    393:             for(int d = 0; d < dim; d++) {
    394:               in.pointlist[idx*dim + d] = array[d];
    395:             }
    396:             in.pointmarkerlist[idx] = serialMesh->getValue(markers, *v_iter);
    397:           }
    398:         }
    399:         const Obj<typename Mesh::label_sequence>& faces      = serialMesh->heightStratum(0);
    400:         const Obj<typename Mesh::numbering_type>& fNumbering = serialMesh->getFactory()->getLocalNumbering(serialMesh, serialMesh->depth());

    402:         in.numberofcorners   = 3;
    403:         in.numberoftriangles = faces->size();
    404:         in.trianglearealist  = (double *) maxVolumes;
    405:         if (in.numberoftriangles > 0) {
    406:           const typename Mesh::label_sequence::iterator fEnd = faces->end();

    408:           PetscMalloc(in.numberoftriangles*in.numberofcorners * sizeof(int), &in.trianglelist);CHKERRXX(ierr);
    409:           if (serialMesh->depth() == 1) {
    410:             for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
    411:               const Obj<typename Mesh::sieve_type::traits::coneSequence>&     cone = serialSieve->cone(*f_iter);
    412:               const typename Mesh::sieve_type::traits::coneSequence::iterator cEnd = cone->end();
    413:               const int                                                       idx  = fNumbering->getIndex(*f_iter);
    414:               int                                                             v    = 0;

    416:               for(typename Mesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
    417:                 in.trianglelist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*c_iter);
    418:               }
    419:             }
    420:           } else if (serialMesh->depth() == 2) {
    421:             for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
    422:               typedef ALE::SieveAlg<Mesh> sieve_alg_type;
    423:               const Obj<typename sieve_alg_type::coneArray>&       cone = sieve_alg_type::nCone(serialMesh, *f_iter, 2);
    424:               const typename Mesh::sieve_type::coneArray::iterator cEnd = cone->end();
    425:               const int                                            idx  = fNumbering->getIndex(*f_iter);
    426:               int                                                  v    = 0;

    428:               for(typename Mesh::sieve_type::coneArray::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
    429:                 in.trianglelist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*c_iter);
    430:               }
    431:             }
    432:           } else {
    433:             throw ALE::Exception("Invalid sieve: Cannot gives sieves of arbitrary depth to Triangle");
    434:           }
    435:         }
    436:         if (serialMesh->depth() == 2) {
    437:           const Obj<typename Mesh::label_sequence>&           edges    = serialMesh->depthStratum(1);
    438: #define NEW_LABEL
    439: #ifdef NEW_LABEL
    440:           for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
    441:             if (serialMesh->getValue(markers, *e_iter)) {
    442:               in.numberofsegments++;
    443:             }
    444:           }
    445:           //std::cout << "Number of segments: " << in.numberofsegments << std::endl;
    446:           if (in.numberofsegments > 0) {
    447:             int s = 0;

    449:             PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);CHKERRXX(ierr);
    450:             PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);CHKERRXX(ierr);
    451:             for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
    452:               const int edgeMarker = serialMesh->getValue(markers, *e_iter);

    454:               if (edgeMarker) {
    455:                 const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = serialSieve->cone(*e_iter);
    456:                 int                                                p    = 0;

    458:                 for(typename Mesh::sieve_type::traits::coneSequence::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
    459:                   in.segmentlist[s*2 + (p++)] = vNumbering->getIndex(*v_iter);
    460:                 }
    461:                 in.segmentmarkerlist[s++] = edgeMarker;
    462:               }
    463:             }
    464:           }
    465: #else
    466:           const Obj<typename Mesh::label_type::baseSequence>& boundary = markers->base();

    468:           in.numberofsegments = 0;
    469:           for(typename Mesh::label_type::baseSequence::iterator b_iter = boundary->begin(); b_iter != boundary->end(); ++b_iter) {
    470:             for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
    471:               if (*b_iter == *e_iter) {
    472:                 in.numberofsegments++;
    473:               }
    474:             }
    475:           }
    476:           if (in.numberofsegments > 0) {
    477:             int s = 0;

    479:             PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);CHKERRXX(ierr);
    480:             PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);CHKERRXX(ierr);
    481:             for(typename Mesh::label_type::baseSequence::iterator b_iter = boundary->begin(); b_iter != boundary->end(); ++b_iter) {
    482:               for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
    483:                 if (*b_iter == *e_iter) {
    484:                   const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = serialSieve->cone(*e_iter);
    485:                   int                                                p    = 0;

    487:                   for(typename Mesh::sieve_type::traits::coneSequence::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
    488:                     in.segmentlist[s*2 + (p++)] = vNumbering->getIndex(*v_iter);
    489:                   }
    490:                   in.segmentmarkerlist[s++] = serialMesh->getValue(markers, *e_iter);
    491:                 }
    492:               }
    493:             }
    494:           }
    495: #endif
    496:         }

    498:         in.numberofholes = 0;
    499:         if (in.numberofholes > 0) {
    500:           PetscMalloc(in.numberofholes * dim * sizeof(int), &in.holelist);CHKERRXX(ierr);
    501:         }
    502:         if (serialMesh->commRank() == 0) {
    503:           std::string args("pqezQra");

    505:           triangulate((char *) args.c_str(), &in, &out, NULL);
    506:         }
    507:         if (in.pointlist)         {PetscFree(in.pointlist);CHKERRXX(ierr);}
    508:         if (in.pointmarkerlist)   {PetscFree(in.pointmarkerlist);CHKERRXX(ierr);}
    509:         if (in.segmentlist)       {PetscFree(in.segmentlist);CHKERRXX(ierr);}
    510:         if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);CHKERRXX(ierr);}
    511:         if (in.trianglelist)      {PetscFree(in.trianglelist);CHKERRXX(ierr);}
    512:         const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(serialMesh->comm(), serialMesh->debug());
    513:         int     numCorners  = 3;
    514:         int     numCells    = out.numberoftriangles;
    515:         int    *cells       = out.trianglelist;
    516:         int     numVertices = out.numberofpoints;
    517:         double *coords      = out.pointlist;

    519:         ALE::SieveBuilder<Mesh>::buildTopology(newSieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, refMesh->getArrowSection("orientation"));
    520:         refMesh->setSieve(newSieve);
    521:         refMesh->stratify();
    522:         ALE::SieveBuilder<Mesh>::buildCoordinates(refMesh, dim, coords);
    523:         const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");

    525:         if (refMesh->commRank() == 0) {
    526:           for(int v = 0; v < out.numberofpoints; v++) {
    527:             if (out.pointmarkerlist[v]) {
    528:               refMesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
    529:             }
    530:           }
    531:           if (interpolate) {
    532:             for(int e = 0; e < out.numberofedges; e++) {
    533:               if (out.edgemarkerlist[e]) {
    534:                 const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
    535:                 const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
    536:                 const Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(vertexA, vertexB, 1);

    538:                 refMesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
    539:               }
    540:             }
    541:           }
    542:         }

    544:         Generator<Mesh>::finiOutput(&out);
    545:         if ((refMesh->commSize() > 1) && (!forceSerial)) {
    546:           return ALE::Distribution<Mesh>::distributeMesh(refMesh);
    547:         }
    548:         return refMesh;
    549:       };
    550:       static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false) {
    551:         Obj<Mesh>                          serialMesh       = ALE::Distribution<Mesh>::unifyMesh(mesh);
    552:         const Obj<typename Mesh::real_section_type> serialMaxVolumes = ALE::Distribution<Mesh>::distributeSection(maxVolumes, serialMesh, serialMesh->getDistSendOverlap(), serialMesh->getDistRecvOverlap());

    554:         return refineMesh(serialMesh, serialMaxVolumes->restrictSpace(), interpolate);
    555:       };
    556:       static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false, const bool forceSerial = false) {
    557:         Obj<Mesh> serialMesh;
    558:         if ((mesh->commSize() > 1) && (!forceSerial)) {
    559:           serialMesh = ALE::Distribution<Mesh>::unifyMesh(mesh);
    560:         } else {
    561:           serialMesh = mesh;
    562:         }
    563:         const int numFaces         = serialMesh->heightStratum(0)->size();
    564:         double   *serialMaxVolumes = new double[numFaces];

    566:         for(int f = 0; f < numFaces; f++) {
    567:           serialMaxVolumes[f] = maxVolume;
    568:         }
    569:         const Obj<Mesh> refMesh = refineMesh(serialMesh, serialMaxVolumes, interpolate, forceSerial);
    570:         delete [] serialMaxVolumes;
    571:         return refMesh;
    572:       };
    573:       static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolumes[], const bool interpolate = false, const bool forceSerial = false, const bool renumber = false) {
    574:         typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
    575:         typedef typename Mesh::real_section_type::value_type real;
    576:         typedef typename Mesh::point_type point_type;
    577:         const int                             dim     = mesh->getDimension();
    578:         const Obj<Mesh>                       refMesh = new Mesh(mesh->comm(), dim, mesh->debug());
    579:         const Obj<typename Mesh::sieve_type>& sieve   = mesh->getSieve();
    580:         struct triangulateio in;
    581:         struct triangulateio out;
    582:         PetscErrorCode       ierr;

    584:         Generator<Mesh>::initInput(&in);
    585:         Generator<Mesh>::initOutput(&out);
    586:         const Obj<typename Mesh::label_sequence>&    vertices    = mesh->depthStratum(0);
    587:         const Obj<typename Mesh::label_type>&        markers     = mesh->getLabel("marker");
    588:         const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
    589:         const Obj<typename Mesh::numbering_type>&    vNumbering  = mesh->getFactory()->getLocalNumbering(mesh, 0);

    591:         in.numberofpoints = vertices->size();
    592:         if (in.numberofpoints > 0) {
    593:           const typename Mesh::label_sequence::iterator vEnd = vertices->end();

    595:           PetscMalloc(in.numberofpoints * dim * sizeof(double), &in.pointlist);CHKERRXX(ierr);
    596:           PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);CHKERRXX(ierr);
    597:           for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
    598:             const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
    599:             const int                                           idx   = vNumbering->getIndex(*v_iter);

    601:             for(int d = 0; d < dim; d++) {
    602:               in.pointlist[idx*dim + d] = array[d];
    603:             }
    604:             in.pointmarkerlist[idx] = mesh->getValue(markers, *v_iter);
    605:           }
    606:         }
    607:         const Obj<typename Mesh::label_sequence>& faces      = mesh->heightStratum(0);
    608:         const Obj<typename Mesh::numbering_type>& fNumbering = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());

    610:         in.numberofcorners   = 3;
    611:         in.numberoftriangles = faces->size();
    612:         in.trianglearealist  = (double *) maxVolumes;
    613:         if (in.numberoftriangles > 0) {
    614:           PetscMalloc(in.numberoftriangles*in.numberofcorners * sizeof(int), &in.trianglelist);CHKERRXX(ierr);
    615:           if (mesh->depth() == 1) {
    616:             ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(3);
    617:             const typename Mesh::label_sequence::iterator fEnd = faces->end();

    619:             for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
    620:               sieve->cone(*f_iter, pV);
    621:               const int         idx  = fNumbering->getIndex(*f_iter);
    622:               const size_t      n    = pV.getSize();
    623:               const point_type *cone = pV.getPoints();

    625:               assert(n == 3);
    626:               for(int v = 0; v < 3; ++v) {
    627:                 in.trianglelist[idx*in.numberofcorners + v] = vNumbering->getIndex(cone[v]);
    628:               }
    629:               pV.clear();
    630:             }
    631:           } else if (mesh->depth() == 2) {
    632:             // Need extra space due to early error checking
    633:             ALE::ISieveVisitor::NConeRetriever<typename Mesh::sieve_type> ncV(*sieve, 4);
    634:             const typename Mesh::label_sequence::iterator fEnd = faces->end();

    636:             for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
    637:               ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*sieve, *f_iter, ncV);
    638:               const int         idx  = fNumbering->getIndex(*f_iter);
    639:               const size_t      n    = ncV.getSize();
    640:               const point_type *cone = ncV.getPoints();

    642:               assert(n == 3);
    643:               for(int v = 0; v < 3; ++v) {
    644:                 in.trianglelist[idx*in.numberofcorners + v] = vNumbering->getIndex(cone[v]);
    645:               }
    646:               ncV.clear();
    647:             }
    648:           } else {
    649:             throw ALE::Exception("Invalid sieve: Cannot gives sieves of arbitrary depth to Triangle");
    650:           }
    651:         }
    652:         if (mesh->depth() == 2) {
    653:           const Obj<typename Mesh::label_sequence>& edges = mesh->depthStratum(1);
    654:           for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
    655:             if (mesh->getValue(markers, *e_iter)) {
    656:               in.numberofsegments++;
    657:             }
    658:           }
    659:           //std::cout << "Number of segments: " << in.numberofsegments << std::endl;
    660:           if (in.numberofsegments > 0) {
    661:             ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(2);
    662:             int s = 0;

    664:             PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);CHKERRXX(ierr);
    665:             PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);CHKERRXX(ierr);
    666:             for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
    667:               const int edgeMarker = mesh->getValue(markers, *e_iter);

    669:               if (edgeMarker) {
    670:                 sieve->cone(*e_iter, pV);
    671:                 const size_t      n    = pV.getSize();
    672:                 const point_type *cone = pV.getPoints();

    674:                 assert(n == 2);
    675:                 for(int v = 0; v < 2; ++v) {
    676:                   in.segmentlist[s*2 + v] = vNumbering->getIndex(cone[v]);
    677:                 }
    678:                 in.segmentmarkerlist[s++] = edgeMarker;
    679:                 pV.clear();
    680:               }
    681:             }
    682:           }
    683:         }
    684:         const typename Mesh::holes_type& holes = mesh->getHoles();

    686:         in.numberofholes = holes.size();
    687:         if (in.numberofholes > 0) {
    688:           PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);CHKERRXX(ierr);
    689:           for(int h = 0; h < in.numberofholes; ++h) {
    690:             for(int d = 0; d < dim; ++d) {
    691:               in.holelist[h*dim+d] = holes[h][d];
    692:             }
    693:           }
    694:         }
    695:         if (mesh->commRank() == 0) {
    696:           std::string args("pqezQra");

    698:           triangulate((char *) args.c_str(), &in, &out, NULL);
    699:         }
    700:         if (in.pointlist)         {PetscFree(in.pointlist);CHKERRXX(ierr);}
    701:         if (in.pointmarkerlist)   {PetscFree(in.pointmarkerlist);CHKERRXX(ierr);}
    702:         if (in.segmentlist)       {PetscFree(in.segmentlist);CHKERRXX(ierr);}
    703:         if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);CHKERRXX(ierr);}
    704:         if (in.trianglelist)      {PetscFree(in.trianglelist);CHKERRXX(ierr);}
    705:         const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
    706:         const Obj<FlexMesh>                  m        = new FlexMesh(mesh->comm(), dim, mesh->debug());
    707:         const Obj<FlexMesh::sieve_type>      newS     = new FlexMesh::sieve_type(m->comm(), m->debug());
    708:         int     numCorners  = 3;
    709:         int     numCells    = out.numberoftriangles;
    710:         int    *cells       = out.trianglelist;
    711:         int     numVertices = out.numberofpoints;
    712:         double *coords      = out.pointlist;
    713:         real   *coordsR;

    715:         ALE::SieveBuilder<FlexMesh>::buildTopology(newS, dim, numCells, cells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
    716:         m->setSieve(newS);
    717:         m->stratify();
    718:         refMesh->setSieve(newSieve);
    719:         std::map<typename Mesh::point_type,typename Mesh::point_type> renumbering;
    720:         ALE::ISieveConverter::convertSieve(*newS, *newSieve, renumbering, renumber);
    721:         refMesh->stratify();
    722:         ALE::ISieveConverter::convertOrientation(*newS, *newSieve, renumbering, m->getArrowSection("orientation").ptr());
    723:         {
    724:           if (sizeof(double) == sizeof(real)) {
    725:             coordsR = (real *) coords;
    726:           } else {
    727:             coordsR = new real[numVertices*dim];
    728:             for(int i = 0; i < numVertices*dim; ++i) coordsR[i] = coords[i];
    729:           }
    730:         }
    731:         ALE::SieveBuilder<Mesh>::buildCoordinates(refMesh, dim, coordsR);
    732:         {
    733:           if (sizeof(double) != sizeof(real)) {
    734:             delete [] coordsR;
    735:           }
    736:         }
    737:         const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");

    739:         if (refMesh->commRank() == 0) {
    740:           for(int v = 0; v < out.numberofpoints; v++) {
    741:             if (out.pointmarkerlist[v]) {
    742:               if (renumber) {
    743:                 refMesh->setValue(newMarkers, renumbering[v+out.numberoftriangles], out.pointmarkerlist[v]);
    744:               } else {
    745:                 refMesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
    746:               }
    747:             }
    748:           }
    749:           if (interpolate) {
    750:             for(int e = 0; e < out.numberofedges; e++) {
    751:               if (out.edgemarkerlist[e]) {
    752:                 const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
    753:                 const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
    754:                 const Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(vertexA, vertexB, 1);

    756:                 refMesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
    757:               }
    758:             }
    759:           }
    760:         }
    761:         refMesh->copyHoles(mesh);
    762:         Generator<Mesh>::finiOutput(&out);
    763: #if 0
    764:         if ((refMesh->commSize() > 1) && (!forceSerial)) {
    765:           return ALE::Distribution<Mesh>::distributeMesh(refMesh);
    766:         }
    767: #endif
    768:         return refMesh;
    769:       };
    770:       static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false, const bool renumber = false) {
    771:         throw ALE::Exception("Not yet implemented");
    772:       };
    773:       static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false, const bool forceSerial = false, const bool renumber = false) {
    774: #if 0
    775:         Obj<Mesh> serialMesh;
    776:         if (mesh->commSize() > 1) {
    777:           serialMesh = ALE::Distribution<Mesh>::unifyMesh(mesh);
    778:         } else {
    779:           serialMesh = mesh;
    780:         }
    781: #endif
    782:         const int numCells         = mesh->heightStratum(0)->size();
    783:         double   *serialMaxVolumes = new double[numCells];

    785:         for(int c = 0; c < numCells; c++) {
    786:           serialMaxVolumes[c] = maxVolume;
    787:         }
    788:         const Obj<Mesh> refMesh = refineMeshV(mesh, serialMaxVolumes, interpolate, forceSerial, renumber);
    789:         delete [] serialMaxVolumes;
    790:         return refMesh;
    791:       };
    792:       static Obj<Mesh> refineMeshLocal(const Obj<Mesh>& mesh, const double maxVolumes[], const bool interpolate = false) {
    793:         const int                    dim     = mesh->getDimension();
    794:         const Obj<Mesh>              refMesh = new Mesh(mesh->comm(), dim, mesh->debug());
    795:         const Obj<typename Mesh::sieve_type>& sieve   = mesh->getSieve();
    796:         struct triangulateio in;
    797:         struct triangulateio out;
    798:         PetscErrorCode       ierr;

    800:         Generator<Mesh>::initInput(&in);
    801:         Generator<Mesh>::initOutput(&out);
    802:         const Obj<typename Mesh::label_sequence>&    vertices    = mesh->depthStratum(0);
    803:         const Obj<typename Mesh::label_type>&        markers     = mesh->getLabel("marker");
    804:         const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
    805:         const Obj<typename Mesh::numbering_type>&    vNumbering  = mesh->getFactory()->getLocalNumbering(mesh, 0);

    807:         in.numberofpoints = vertices->size();
    808:         if (in.numberofpoints > 0) {
    809:           PetscMalloc(in.numberofpoints * dim * sizeof(double), &in.pointlist);CHKERRXX(ierr);
    810:           PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);CHKERRXX(ierr);
    811:           for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
    812:             const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
    813:             const int                                  idx   = vNumbering->getIndex(*v_iter);

    815:             for(int d = 0; d < dim; d++) {
    816:               in.pointlist[idx*dim + d] = array[d];
    817:             }
    818:             in.pointmarkerlist[idx] = mesh->getValue(markers, *v_iter);
    819:           }
    820:         }
    821:         const Obj<typename Mesh::label_sequence>& faces      = mesh->heightStratum(0);
    822:         const Obj<typename Mesh::numbering_type>& fNumbering = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());

    824:         in.numberofcorners   = 3;
    825:         in.numberoftriangles = faces->size();
    826:         in.trianglearealist  = (double *) maxVolumes;
    827:         if (in.numberoftriangles > 0) {
    828:           PetscMalloc(in.numberoftriangles*in.numberofcorners * sizeof(int), &in.trianglelist);CHKERRXX(ierr);
    829:           if (mesh->depth() == 1) {
    830:             for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != faces->end(); ++f_iter) {
    831:               const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = sieve->cone(*f_iter);
    832:               const int                                          idx  = fNumbering->getIndex(*f_iter);
    833:               int                                                v    = 0;

    835:               for(typename Mesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
    836:                 in.trianglelist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*c_iter);
    837:               }
    838:             }
    839:           } else if (mesh->depth() == 2) {
    840:             for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != faces->end(); ++f_iter) {
    841:               typedef ALE::SieveAlg<Mesh> sieve_alg_type;
    842:               const Obj<typename sieve_alg_type::coneArray>& cone = sieve_alg_type::nCone(mesh, *f_iter, 2);
    843:               const int                             idx  = fNumbering->getIndex(*f_iter);
    844:               int                                   v    = 0;

    846:               for(typename Mesh::sieve_type::coneArray::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
    847:                 in.trianglelist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*c_iter);
    848:               }
    849:             }
    850:           } else {
    851:             throw ALE::Exception("Invalid sieve: Cannot gives sieves of arbitrary depth to Triangle");
    852:           }
    853:         }
    854:         if (mesh->depth() == 2) {
    855:           const Obj<typename Mesh::label_sequence>& edges = mesh->depthStratum(1);
    856:           for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
    857:             if (mesh->getValue(markers, *e_iter)) {
    858:               in.numberofsegments++;
    859:             }
    860:           }
    861:           //std::cout << "Number of segments: " << in.numberofsegments << std::endl;
    862:           if (in.numberofsegments > 0) {
    863:             int s = 0;

    865:             PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);CHKERRXX(ierr);
    866:             PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);CHKERRXX(ierr);
    867:             for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
    868:               const int edgeMarker = mesh->getValue(markers, *e_iter);

    870:               if (edgeMarker) {
    871:                 const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = sieve->cone(*e_iter);
    872:                 int                                                p    = 0;

    874:                 for(typename Mesh::sieve_type::traits::coneSequence::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
    875:                   in.segmentlist[s*2 + (p++)] = vNumbering->getIndex(*v_iter);
    876:                 }
    877:                 in.segmentmarkerlist[s++] = edgeMarker;
    878:               }
    879:             }
    880:           }
    881:         }

    883:         in.numberofholes = 0;
    884:         if (in.numberofholes > 0) {
    885:           PetscMalloc(in.numberofholes * dim * sizeof(int), &in.holelist);CHKERRXX(ierr);
    886:         }
    887:         std::string args("pqezQra");

    889:         triangulate((char *) args.c_str(), &in, &out, NULL);
    890:         if (in.pointlist)         {PetscFree(in.pointlist);CHKERRXX(ierr);}
    891:         if (in.pointmarkerlist)   {PetscFree(in.pointmarkerlist);CHKERRXX(ierr);}
    892:         if (in.segmentlist)       {PetscFree(in.segmentlist);CHKERRXX(ierr);}
    893:         if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);CHKERRXX(ierr);}
    894:         if (in.trianglelist)      {PetscFree(in.trianglelist);CHKERRXX(ierr);}
    895:         const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
    896:         int     numCorners  = 3;
    897:         int     numCells    = out.numberoftriangles;
    898:         int    *cells       = out.trianglelist;
    899:         int     numVertices = out.numberofpoints;
    900:         double *coords      = out.pointlist;

    902:         ALE::SieveBuilder<Mesh>::buildTopologyMultiple(newSieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, refMesh->getArrowSection("orientation"));
    903:         refMesh->setSieve(newSieve);
    904:         refMesh->stratify();
    905:         ALE::SieveBuilder<Mesh>::buildCoordinatesMultiple(refMesh, dim, coords);
    906:         const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");

    908:         for(int v = 0; v < out.numberofpoints; v++) {
    909:           if (out.pointmarkerlist[v]) {
    910:             refMesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
    911:           }
    912:         }
    913:         if (interpolate) {
    914:           for(int e = 0; e < out.numberofedges; e++) {
    915:             if (out.edgemarkerlist[e]) {
    916:               const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
    917:               const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
    918:               const Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(vertexA, vertexB, 1);

    920:               refMesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
    921:             }
    922:           }
    923:         }
    924:         Generator<Mesh>::finiOutput(&out);
    925:         return refMesh;
    926:       };
    927:       static Obj<Mesh> refineMeshLocal(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false) {
    928:         const int numLocalFaces   = mesh->heightStratum(0)->size();
    929:         double   *localMaxVolumes = new double[numLocalFaces];

    931:         for(int f = 0; f < numLocalFaces; f++) {
    932:           localMaxVolumes[f] = maxVolume;
    933:         }
    934:         const Obj<Mesh> refMesh = refineMeshLocal(mesh, localMaxVolumes, interpolate);
    935:         const Obj<typename Mesh::sieve_type> refSieve = refMesh->getSieve();
    936:         delete [] localMaxVolumes;
    937: #if 0
    938:         typedef typename ALE::New::Completion<Mesh, typename Mesh::sieve_type::point_type> sieveCompletion;
    939:         // This is where we enforce consistency over the overlap
    940:         //   We need somehow to update the overlap to account for the new stuff
    941:         //
    942:         //   1) Since we are refining only, the vertices are invariant
    943:         //   2) We need to make another label for the interprocess boundaries so
    944:         //      that Triangle will respect them
    945:         //   3) We then throw all that label into the new overlap
    946:         //
    947:         // Alternative: Figure out explicitly which segments were refined, and then
    948:         //   communicated the refinement over the old overlap. Use this info to locally
    949:         //   construct the new overlap and flip to get a decent mesh
    950:         sieveCompletion::scatterCones(refSieve, refSieve, reMesh->getDistSendOverlap(), refMesh->getDistRecvOverlap(), refMesh);
    951: #endif
    952:         return refMesh;
    953:       };
    954:     };
    955:   };
    956: #endif
    957: #ifdef PETSC_HAVE_TETGEN
    958:   namespace TetGen {
    959:     template<typename Mesh>
    960:     class Generator {
    961:     public:
    962:       static Obj<Mesh> generateMesh(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false) {
    963:         typedef ALE::SieveAlg<Mesh> sieve_alg_type;
    964:         const int         dim   = 3;
    965:         Obj<Mesh>         mesh  = new Mesh(boundary->comm(), dim, boundary->debug());
    966:         const PetscMPIInt rank  = mesh->commRank();
    967:         bool              createConvexHull = false;
    968:         ::tetgenio        in;
    969:         ::tetgenio        out;

    971:         const Obj<typename Mesh::label_sequence>&    vertices    = boundary->depthStratum(0);
    972:         const Obj<typename Mesh::numbering_type>&    vNumbering  = boundary->getFactory()->getLocalNumbering(boundary, 0);
    973:         const Obj<typename Mesh::real_section_type>& coordinates = boundary->getRealSection("coordinates");
    974:         const Obj<typename Mesh::label_type>&        markers     = boundary->getLabel("marker");

    977:         in.numberofpoints = vertices->size();
    978:         if (in.numberofpoints > 0) {
    979:           in.pointlist       = new double[in.numberofpoints*dim];
    980:           in.pointmarkerlist = new int[in.numberofpoints];
    981:           for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
    982:             const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
    983:             const int                                  idx   = vNumbering->getIndex(*v_iter);

    985:             for(int d = 0; d < dim; d++) {
    986:               in.pointlist[idx*dim + d] = array[d];
    987:             }
    988:             in.pointmarkerlist[idx] = boundary->getValue(markers, *v_iter);
    989:           }
    990:         }

    992:         if (boundary->depth() != 0) {  //our boundary mesh COULD be just a pointset; in which case depth = height = 0;
    993:           const Obj<typename Mesh::label_sequence>& facets     = boundary->depthStratum(boundary->depth());
    994:           //PetscPrintf(boundary->comm(), "%d facets on the boundary\n", facets->size());
    995:           const Obj<typename Mesh::numbering_type>& fNumbering = boundary->getFactory()->getLocalNumbering(boundary, boundary->depth());

    997:           in.numberoffacets = facets->size();
    998:           if (in.numberoffacets > 0) {
    999:             in.facetlist       = new tetgenio::facet[in.numberoffacets];
    1000:             in.facetmarkerlist = new int[in.numberoffacets];
    1001:             for(typename Mesh::label_sequence::iterator f_iter = facets->begin(); f_iter != facets->end(); ++f_iter) {
    1002:               const Obj<typename sieve_alg_type::coneArray>& cone = sieve_alg_type::nCone(boundary, *f_iter, boundary->depth());
    1003:               const int                             idx  = fNumbering->getIndex(*f_iter);

    1005:               in.facetlist[idx].numberofpolygons = 1;
    1006:               in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
    1007:               in.facetlist[idx].numberofholes    = 0;
    1008:               in.facetlist[idx].holelist         = NULL;

    1010:               tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
    1011:               int                c    = 0;

    1013:               poly->numberofvertices = cone->size();
    1014:               poly->vertexlist       = new int[poly->numberofvertices];
    1015:               for(typename sieve_alg_type::coneArray::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
    1016:                 const int vIdx = vNumbering->getIndex(*c_iter);

    1018:                 poly->vertexlist[c++] = vIdx;
    1019:               }
    1020:               in.facetmarkerlist[idx] = boundary->getValue(markers, *f_iter);
    1021:             }
    1022:           }
    1023:         }else {
    1024:           createConvexHull = true;
    1025:         }

    1027:         in.numberofholes = 0;
    1028:         if (rank == 0) {
    1029:           // Normal operation
    1030:           std::string args("pqezQ");
    1031:           //constrained operation
    1032:           if (constrained) {
    1033:             args = "pezQ";
    1034:             if (createConvexHull) {
    1035:               args = "ezQ";
    1036:               //PetscPrintf(boundary->comm(), "createConvexHull\n");
    1037:             }
    1038:           }
    1039:           // Just make tetrahedrons
    1040: //           std::string args("efzV");
    1041:           // Adds a center point
    1042: //           std::string args("pqezQi");
    1043: //           in.numberofaddpoints = 1;
    1044: //           in.addpointlist      = new double[in.numberofaddpoints*dim];
    1045: //           in.addpointlist[0]   = 0.5;
    1046: //           in.addpointlist[1]   = 0.5;
    1047: //           in.addpointlist[2]   = 0.5;

    1049:           //if (createConvexHull) args += "c";  NOT SURE, but this was basically unused before, and the convex hull should be filled in if "p" isn't included.
    1050:           ::tetrahedralize((char *) args.c_str(), &in, &out);
    1051:         }
    1052:         const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
    1053:         int     numCorners  = 4;
    1054:         int     numCells    = out.numberoftetrahedra;
    1055:         int    *cells       = out.tetrahedronlist;
    1056:         int     numVertices = out.numberofpoints;
    1057:         double *coords      = out.pointlist;

    1059:         if (!interpolate) {
    1060:           for(int c = 0; c < numCells; ++c) {
    1061:             int tmp = cells[c*4+0];
    1062:             cells[c*4+0] = cells[c*4+1];
    1063:             cells[c*4+1] = tmp;
    1064:           }
    1065:         }
    1066:         ALE::SieveBuilder<Mesh>::buildTopology(newSieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, mesh->getArrowSection("orientation"));
    1067:         mesh->setSieve(newSieve);
    1068:         mesh->stratify();
    1069:         ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, dim, coords);
    1070:         const Obj<typename Mesh::label_type>& newMarkers = mesh->createLabel("marker");

    1072:         for(int v = 0; v < out.numberofpoints; v++) {
    1073:           if (out.pointmarkerlist[v]) {
    1074:             mesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
    1075:           }
    1076:         }
    1077:         if (interpolate) {
    1078:           if (out.edgemarkerlist) {
    1079:             for(int e = 0; e < out.numberofedges; e++) {
    1080:               if (out.edgemarkerlist[e]) {
    1081:                 typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
    1082:                 typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
    1083:                 Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(endpointA, endpointB, 1);

    1085:                 mesh->setValue(newMarkers, *edge->begin(), out.edgemarkerlist[e]);
    1086:               }
    1087:             }
    1088:           }
    1089:           if (out.trifacemarkerlist) {
    1090:             // Work around TetGen bug for raw tetrahedralization
    1091:             //   The boundary faces are 0,1,4,5,8,9,11,12,13,15,16,17
    1092: //             for(int f = 0; f < out.numberoftrifaces; f++) {
    1093: //               if (out.trifacemarkerlist[f]) {
    1094: //                 out.trifacemarkerlist[f] = 0;
    1095: //               } else {
    1096: //                 out.trifacemarkerlist[f] = 1;
    1097: //               }
    1098: //             }
    1099:             for(int f = 0; f < out.numberoftrifaces; f++) {
    1100:               if (out.trifacemarkerlist[f]) {
    1101:                 typename Mesh::point_type cornerA(out.trifacelist[f*3+0]+out.numberoftetrahedra);
    1102:                 typename Mesh::point_type cornerB(out.trifacelist[f*3+1]+out.numberoftetrahedra);
    1103:                 typename Mesh::point_type cornerC(out.trifacelist[f*3+2]+out.numberoftetrahedra);
    1104:                 Obj<typename Mesh::sieve_type::supportSet> corners = typename Mesh::sieve_type::supportSet();
    1105:                 Obj<typename Mesh::sieve_type::supportSet> edges   = typename Mesh::sieve_type::supportSet();
    1106:                 corners->clear();corners->insert(cornerA);corners->insert(cornerB);
    1107:                 edges->insert(*newSieve->nJoin1(corners)->begin());
    1108:                 corners->clear();corners->insert(cornerB);corners->insert(cornerC);
    1109:                 edges->insert(*newSieve->nJoin1(corners)->begin());
    1110:                 corners->clear();corners->insert(cornerC);corners->insert(cornerA);
    1111:                 edges->insert(*newSieve->nJoin1(corners)->begin());
    1112:                 const typename Mesh::point_type          face       = *newSieve->nJoin1(edges)->begin();
    1113:                 const int                       faceMarker = out.trifacemarkerlist[f];
    1114:                 const Obj<typename Mesh::coneArray>      closure    = sieve_alg_type::closure(mesh, face);
    1115:                 const typename Mesh::coneArray::iterator end        = closure->end();

    1117:                 for(typename Mesh::coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
    1118:                   mesh->setValue(newMarkers, *cl_iter, faceMarker);
    1119:                 }
    1120:               }
    1121:             }
    1122:           }
    1123:         }
    1124:         return mesh;
    1125:       };
    1128:       static Obj<Mesh> generateMeshV(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false, const bool renumber = false) {
    1129:         typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
    1130:         typedef typename Mesh::real_section_type::value_type real;
    1131:         const int                             dim   = 3;
    1132:         Obj<Mesh>                             mesh  = new Mesh(boundary->comm(), dim, boundary->debug());
    1133:         const Obj<typename Mesh::sieve_type>& sieve = boundary->getSieve();
    1134:         const PetscMPIInt                     rank  = mesh->commRank();
    1135:         bool                                  createConvexHull = false;
    1136:         PetscErrorCode                        ierr;
    1137:         ::tetgenio in;
    1138:         ::tetgenio out;

    1140:         const Obj<typename Mesh::label_sequence>&    vertices    = boundary->depthStratum(0);
    1141:         const Obj<typename Mesh::label_type>&        markers     = boundary->getLabel("marker");
    1142:         const Obj<typename Mesh::real_section_type>& coordinates = boundary->getRealSection("coordinates");
    1143:         const Obj<typename Mesh::numbering_type>&    vNumbering  = boundary->getFactory()->getLocalNumbering(boundary, 0);

    1145:         in.numberofpoints = vertices->size();
    1146:         if (in.numberofpoints > 0) {
    1147:           const typename Mesh::label_sequence::iterator vEnd = vertices->end();

    1149:           in.pointlist       = new double[in.numberofpoints*dim];
    1150:           in.pointmarkerlist = new int[in.numberofpoints];
    1151:           for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
    1152:             const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
    1153:             const int                                           idx   = vNumbering->getIndex(*v_iter);

    1155:             for(int d = 0; d < dim; ++d) {
    1156:               in.pointlist[idx*dim + d] = array[d];
    1157:             }
    1158:             in.pointmarkerlist[idx] = boundary->getValue(markers, *v_iter);
    1159:           }
    1160:         }

    1162:         // Our boundary mesh COULD be just a pointset; in which case depth = height = 0;
    1163:         if (boundary->depth() != 0) {
    1164:           const Obj<typename Mesh::label_sequence>& facets     = boundary->depthStratum(boundary->depth());
    1165:           const Obj<typename Mesh::numbering_type>& fNumbering = boundary->getFactory()->getLocalNumbering(boundary, boundary->depth());

    1167:           in.numberoffacets = facets->size();
    1168:           if (in.numberoffacets > 0) {
    1169:             ALE::ISieveVisitor::NConeRetriever<typename Mesh::sieve_type> ncV(*sieve, 5);
    1170:             const typename Mesh::label_sequence::iterator fEnd = facets->end();

    1172:             in.facetlist       = new tetgenio::facet[in.numberoffacets];
    1173:             in.facetmarkerlist = new int[in.numberoffacets];
    1174:             for(typename Mesh::label_sequence::iterator f_iter = facets->begin(); f_iter != fEnd; ++f_iter) {
    1175:               ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*sieve, *f_iter, ncV);
    1176:               const int idx  = fNumbering->getIndex(*f_iter);

    1178:               in.facetlist[idx].numberofpolygons = 1;
    1179:               in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
    1180:               in.facetlist[idx].numberofholes    = 0;
    1181:               in.facetlist[idx].holelist         = NULL;

    1183:               tetgenio::polygon               *poly = in.facetlist[idx].polygonlist;
    1184:               const size_t                     n    = ncV.getSize();
    1185:               const typename Mesh::point_type *cone = ncV.getPoints();

    1187:               poly->numberofvertices = n;
    1188:               poly->vertexlist       = new int[poly->numberofvertices];
    1189:               for(size_t c = 0; c < n; ++c) {
    1190:                 poly->vertexlist[c] = vNumbering->getIndex(cone[c]);
    1191:               }
    1192:               in.facetmarkerlist[idx] = boundary->getValue(markers, *f_iter);
    1193:               ncV.clear();
    1194:             }
    1195:           }
    1196:         } else {
    1197:           createConvexHull = true;
    1198:         }
    1199:         const typename Mesh::holes_type& holes = mesh->getHoles();

    1201:         in.numberofholes = holes.size();
    1202:         if (in.numberofholes > 0) {
    1203:           PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);CHKERRXX(ierr);
    1204:           for(int h = 0; h < in.numberofholes; ++h) {
    1205:             for(int d = 0; d < dim; ++d) {
    1206:               in.holelist[h*dim+d] = holes[h][d];
    1207:             }
    1208:           }
    1209:         }
    1210:         if (rank == 0) {
    1211:           // Normal operation
    1212:           std::string args("pqezQ");
    1213:           // Constrained operation
    1214:           if (constrained) {
    1215:             args = "pezQ";
    1216:             if (createConvexHull) {
    1217:               args = "ezQ";
    1218:             }
    1219:           }
    1220:           // Just make tetrahedrons
    1221: //           std::string args("efzV");
    1222:           // Adds a center point
    1223: //           std::string args("pqezQi");
    1224: //           in.numberofaddpoints = 1;
    1225: //           in.addpointlist      = new double[in.numberofaddpoints*dim];
    1226: //           in.addpointlist[0]   = 0.5;
    1227: //           in.addpointlist[1]   = 0.5;
    1228: //           in.addpointlist[2]   = 0.5;
    1229:           ::tetrahedralize((char *) args.c_str(), &in, &out);
    1230:         }
    1231:         const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
    1232:         const Obj<FlexMesh>                  m        = new FlexMesh(mesh->comm(), dim, mesh->debug());
    1233:         const Obj<FlexMesh::sieve_type>      newS     = new FlexMesh::sieve_type(m->comm(), m->debug());
    1234:         int     numCorners  = 4;
    1235:         int     numCells    = out.numberoftetrahedra;
    1236:         int    *cells       = out.tetrahedronlist;
    1237:         int     numVertices = out.numberofpoints;
    1238:         double *coords      = out.pointlist;
    1239:         real   *coordsR;

    1241:         if (!interpolate) {
    1242:           // TetGen reports tetrahedra with the opposite orientation from what we expect
    1243:           for(int c = 0; c < numCells; ++c) {
    1244:             int tmp = cells[c*4+0];
    1245:             cells[c*4+0] = cells[c*4+1];
    1246:             cells[c*4+1] = tmp;
    1247:           }
    1248:         }
    1249:         ALE::SieveBuilder<FlexMesh>::buildTopology(newS, dim, numCells, cells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
    1250:         m->setSieve(newS);
    1251:         m->stratify();
    1252:         mesh->setSieve(newSieve);
    1253:         std::map<typename Mesh::point_type,typename Mesh::point_type> renumbering;
    1254:         ALE::ISieveConverter::convertSieve(*newS, *newSieve, renumbering, renumber);
    1255:         mesh->stratify();
    1256:         ALE::ISieveConverter::convertOrientation(*newS, *newSieve, renumbering, m->getArrowSection("orientation").ptr());
    1257:         {
    1258:           if (sizeof(double) == sizeof(real)) {
    1259:             coordsR = (real *) coords;
    1260:           } else {
    1261:             coordsR = new real[numVertices*dim];
    1262:             for(int i = 0; i < numVertices*dim; ++i) coordsR[i] = coords[i];
    1263:           }
    1264:         }
    1265:         ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, dim, coordsR);
    1266:         {
    1267:           if (sizeof(double) != sizeof(real)) {
    1268:             delete [] coordsR;
    1269:           }
    1270:         }
    1271:         const Obj<typename Mesh::label_type>& newMarkers = mesh->createLabel("marker");

    1273:         for(int v = 0; v < out.numberofpoints; v++) {
    1274:           if (out.pointmarkerlist[v]) {
    1275:             if (renumber) {
    1276:               mesh->setValue(newMarkers, renumbering[v+out.numberoftetrahedra], out.pointmarkerlist[v]);
    1277:             } else {
    1278:               mesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
    1279:             }
    1280:           }
    1281:         }
    1282:         if (interpolate) {
    1283:           // This does not work anymore (edgemarkerlist is always empty). I tried -ee and it gave bogus results
    1284:           if (out.edgemarkerlist) {
    1285:             for(int e = 0; e < out.numberofedges; e++) {
    1286:               if (out.edgemarkerlist[e]) {
    1287:                 typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
    1288:                 typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
    1289:                 Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(endpointA, endpointB, 1);

    1291:                 if (renumber) {
    1292:                   mesh->setValue(newMarkers, renumbering[*edge->begin()], out.edgemarkerlist[e]);
    1293:                 } else {
    1294:                   mesh->setValue(newMarkers, *edge->begin(), out.edgemarkerlist[e]);
    1295:                 }
    1296:               }
    1297:             }
    1298:           }
    1299:           if (out.trifacemarkerlist) {
    1300:             // Work around TetGen bug for raw tetrahedralization
    1301:             //   The boundary faces are 0,1,4,5,8,9,11,12,13,15,16,17
    1302: //             for(int f = 0; f < out.numberoftrifaces; f++) {
    1303: //               if (out.trifacemarkerlist[f]) {
    1304: //                 out.trifacemarkerlist[f] = 0;
    1305: //               } else {
    1306: //                 out.trifacemarkerlist[f] = 1;
    1307: //               }
    1308: //             }
    1309:             for(int f = 0; f < out.numberoftrifaces; f++) {
    1310:               if (out.trifacemarkerlist[f]) {
    1311:                 typename Mesh::point_type cornerA = out.trifacelist[f*3+0]+out.numberoftetrahedra;
    1312:                 typename Mesh::point_type cornerB = out.trifacelist[f*3+1]+out.numberoftetrahedra;
    1313:                 typename Mesh::point_type cornerC = out.trifacelist[f*3+2]+out.numberoftetrahedra;
    1314:                 Obj<typename Mesh::sieve_type::supportSet> corners = typename Mesh::sieve_type::supportSet();
    1315:                 Obj<typename Mesh::sieve_type::supportSet> edges   = typename Mesh::sieve_type::supportSet();
    1316:                 corners->clear();corners->insert(cornerA);corners->insert(cornerB);
    1317:                 typename Mesh::point_type edgeA = *newS->nJoin1(corners)->begin();
    1318:                 corners->clear();corners->insert(cornerB);corners->insert(cornerC);
    1319:                 typename Mesh::point_type edgeB = *newS->nJoin1(corners)->begin();
    1320:                 corners->clear();corners->insert(cornerC);corners->insert(cornerA);
    1321:                 typename Mesh::point_type edgeC = *newS->nJoin1(corners)->begin();
    1322:                 edges->insert(edgeA); edges->insert(edgeB); edges->insert(edgeC);
    1323:                 ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(30);
    1324:                 typename Mesh::point_type face       = *newS->nJoin1(edges)->begin();
    1325:                 const int                 faceMarker = out.trifacemarkerlist[f];

    1327:                 if (renumber) {face = renumbering[face];}
    1328:                 ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*newSieve, face, pV);
    1329:                 const size_t                     n    = pV.getSize();
    1330:                 const typename Mesh::point_type *cone = pV.getPoints();

    1332:                 for(size_t c = 0; c < n; ++c) {
    1333:                   mesh->setValue(newMarkers, cone[c], faceMarker);
    1334:                 }
    1335:                 pV.clear();
    1336:               }
    1337:             }
    1338:           }
    1339:         }
    1340:         mesh->copyHoles(boundary);
    1341:         return mesh;
    1342:       };
    1343:     };
    1344:     template<typename Mesh>
    1345:     class Refiner {
    1346:     public:
    1347:       static Obj<Mesh> refineMesh(const Obj<Mesh>& serialMesh, const double maxVolumes[], const bool interpolate = false) {
    1348:         typedef ALE::SieveAlg<Mesh> sieve_alg_type;
    1349:         const int       dim     = serialMesh->getDimension();
    1350:         const int       depth   = serialMesh->depth();
    1351:         const Obj<Mesh> refMesh = new Mesh(serialMesh->comm(), dim, serialMesh->debug());
    1352:         ::tetgenio      in;
    1353:         ::tetgenio      out;

    1355:         const Obj<typename Mesh::label_sequence>&    vertices    = serialMesh->depthStratum(0);
    1356:         const Obj<typename Mesh::label_type>&        markers     = serialMesh->getLabel("marker");
    1357:         const Obj<typename Mesh::real_section_type>& coordinates = serialMesh->getRealSection("coordinates");
    1358:         const Obj<typename Mesh::numbering_type>&    vNumbering  = serialMesh->getFactory()->getLocalNumbering(serialMesh, 0);

    1360:         in.numberofpoints = vertices->size();
    1361:         if (in.numberofpoints > 0) {
    1362:           in.pointlist       = new double[in.numberofpoints*dim];
    1363:           in.pointmarkerlist = new int[in.numberofpoints];
    1364:           for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
    1365:             const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
    1366:             const int                                  idx   = vNumbering->getIndex(*v_iter);

    1368:             for(int d = 0; d < dim; d++) {
    1369:               in.pointlist[idx*dim + d] = array[d];
    1370:             }
    1371:             in.pointmarkerlist[idx] = serialMesh->getValue(markers, *v_iter);
    1372:           }
    1373:         }
    1374:         const Obj<typename Mesh::label_sequence>& cells      = serialMesh->heightStratum(0);
    1375:         const Obj<typename Mesh::numbering_type>& cNumbering = serialMesh->getFactory()->getLocalNumbering(serialMesh, depth);

    1377:         in.numberofcorners       = 4;
    1378:         in.numberoftetrahedra    = cells->size();
    1379:         in.tetrahedronvolumelist = (double *) maxVolumes;
    1380:         if (in.numberoftetrahedra > 0) {
    1381:           in.tetrahedronlist     = new int[in.numberoftetrahedra*in.numberofcorners];
    1382:           for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
    1383:             typedef ALE::SieveAlg<Mesh> sieve_alg_type;
    1384:             const Obj<typename sieve_alg_type::coneArray>& cone = sieve_alg_type::nCone(serialMesh, *c_iter, depth);
    1385:             const int                             idx  = cNumbering->getIndex(*c_iter);
    1386:             int                                   v    = 0;

    1388:             for(typename Mesh::sieve_type::coneArray::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
    1389:               in.tetrahedronlist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*v_iter);
    1390:             }
    1391:           }
    1392:         }
    1393:         if (serialMesh->depth() == 3) {
    1394:           const Obj<typename Mesh::label_sequence>& boundary = serialMesh->getLabelStratum("marker", 1);

    1396:           in.numberoftrifaces = 0;
    1397:           for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != boundary->end(); ++b_iter) {
    1398:             if (serialMesh->height(*b_iter) == 1) {
    1399:               in.numberoftrifaces++;
    1400:             }
    1401:           }
    1402:           if (in.numberoftrifaces > 0) {
    1403:             int f = 0;

    1405:             in.trifacelist       = new int[in.numberoftrifaces*3];
    1406:             in.trifacemarkerlist = new int[in.numberoftrifaces];
    1407:             for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != boundary->end(); ++b_iter) {
    1408:               if (serialMesh->height(*b_iter) == 1) {
    1409:                 const Obj<typename Mesh::coneArray>& cone = sieve_alg_type::nCone(serialMesh, *b_iter, 2);
    1410:                 int                         p    = 0;

    1412:                 for(typename Mesh::coneArray::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
    1413:                   in.trifacelist[f*3 + (p++)] = vNumbering->getIndex(*v_iter);
    1414:                 }
    1415:                 in.trifacemarkerlist[f++] = serialMesh->getValue(markers, *b_iter);
    1416:               }
    1417:             }
    1418:           }
    1419:         }

    1421:         in.numberofholes = 0;
    1422:         if (serialMesh->commRank() == 0) {
    1423:           std::string args("qezQra");

    1425:           ::tetrahedralize((char *) args.c_str(), &in, &out);
    1426:         }
    1427:         in.tetrahedronvolumelist = NULL;
    1428:         const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(refMesh->comm(), refMesh->debug());
    1429:         int     numCorners  = 4;
    1430:         int     numCells    = out.numberoftetrahedra;
    1431:         int    *newCells       = out.tetrahedronlist;
    1432:         int     numVertices = out.numberofpoints;
    1433:         double *coords      = out.pointlist;

    1435:         if (!interpolate) {
    1436:           for(int c = 0; c < numCells; ++c) {
    1437:             int tmp = newCells[c*4+0];
    1438:             newCells[c*4+0] = newCells[c*4+1];
    1439:             newCells[c*4+1] = tmp;
    1440:           }
    1441:         }
    1442:         ALE::SieveBuilder<Mesh>::buildTopology(newSieve, dim, numCells, newCells, numVertices, interpolate, numCorners, -1, refMesh->getArrowSection("orientation"));
    1443:         refMesh->setSieve(newSieve);
    1444:         refMesh->stratify();
    1445:         ALE::SieveBuilder<Mesh>::buildCoordinates(refMesh, dim, coords);
    1446:         const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");

    1449:         for(int v = 0; v < out.numberofpoints; v++) {
    1450:           if (out.pointmarkerlist[v]) {
    1451:             refMesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
    1452:           }
    1453:         }
    1454:         if (interpolate) {
    1455:           if (out.edgemarkerlist) {
    1456:             for(int e = 0; e < out.numberofedges; e++) {
    1457:               if (out.edgemarkerlist[e]) {
    1458:                 typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
    1459:                 typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
    1460:                 Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(endpointA, endpointB, 1);

    1462:                 refMesh->setValue(newMarkers, *edge->begin(), out.edgemarkerlist[e]);
    1463:               }
    1464:             }
    1465:           }
    1466:           if (out.trifacemarkerlist) {
    1467:             for(int f = 0; f < out.numberoftrifaces; f++) {
    1468:               if (out.trifacemarkerlist[f]) {
    1469:                 typename Mesh::point_type cornerA(out.trifacelist[f*3+0]+out.numberoftetrahedra);
    1470:                 typename Mesh::point_type cornerB(out.trifacelist[f*3+1]+out.numberoftetrahedra);
    1471:                 typename Mesh::point_type cornerC(out.trifacelist[f*3+2]+out.numberoftetrahedra);
    1472:                 Obj<typename Mesh::sieve_type::supportSet> corners = typename Mesh::sieve_type::supportSet();
    1473:                 Obj<typename Mesh::sieve_type::supportSet> edges   = typename Mesh::sieve_type::supportSet();
    1474:                 corners->clear();corners->insert(cornerA);corners->insert(cornerB);
    1475:                 edges->insert(*newSieve->nJoin1(corners)->begin());
    1476:                 corners->clear();corners->insert(cornerB);corners->insert(cornerC);
    1477:                 edges->insert(*newSieve->nJoin1(corners)->begin());
    1478:                 corners->clear();corners->insert(cornerC);corners->insert(cornerA);
    1479:                 edges->insert(*newSieve->nJoin1(corners)->begin());
    1480:                 const typename Mesh::point_type          face       = *newSieve->nJoin1(edges)->begin();
    1481:                 const int                       faceMarker = out.trifacemarkerlist[f];
    1482:                 const Obj<typename Mesh::coneArray>      closure    = sieve_alg_type::closure(refMesh, face);
    1483:                 const typename Mesh::coneArray::iterator end        = closure->end();

    1485:                 for(typename Mesh::coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
    1486:                   refMesh->setValue(newMarkers, *cl_iter, faceMarker);
    1487:                 }
    1488:               }
    1489:             }
    1490:           }
    1491:         }
    1492:         if (refMesh->commSize() > 1) {
    1493:           return ALE::Distribution<Mesh>::distributeMesh(refMesh);
    1494:         }
    1495:         return refMesh;
    1496:       };
    1497:       static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false) {
    1498:         Obj<Mesh>                          serialMesh       = ALE::Distribution<Mesh>::unifyMesh(mesh);
    1499:         const Obj<typename Mesh::real_section_type> serialMaxVolumes = ALE::Distribution<Mesh>::distributeSection(maxVolumes, serialMesh, serialMesh->getDistSendOverlap(), serialMesh->getDistRecvOverlap());

    1501:         return refineMesh(serialMesh, serialMaxVolumes->restrictSpace(), interpolate);
    1502:       };
    1503:       static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false) {
    1504:         Obj<Mesh> serialMesh;
    1505:         if (mesh->commSize() > 1) {
    1506:           serialMesh = ALE::Distribution<Mesh>::unifyMesh(mesh);
    1507:         } else {
    1508:           serialMesh = mesh;
    1509:         }
    1510:         const int numCells         = serialMesh->heightStratum(0)->size();
    1511:         double   *serialMaxVolumes = new double[numCells];

    1513:         for(int c = 0; c < numCells; c++) {
    1514:           serialMaxVolumes[c] = maxVolume;
    1515:         }
    1516:         const Obj<Mesh> refMesh = refineMesh(serialMesh, serialMaxVolumes, interpolate);
    1517:         delete [] serialMaxVolumes;
    1518:         return refMesh;
    1519:       };
    1520:       static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolumes[], const bool interpolate = false, const bool forceSerial = false, const bool renumber = false) {
    1521:         typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
    1522:         typedef typename Mesh::real_section_type::value_type real;
    1523:         typedef typename Mesh::point_type point_type;
    1524:         const int                             dim     = mesh->getDimension();
    1525:         const int                             depth   = mesh->depth();
    1526:         const Obj<Mesh>                       refMesh = new Mesh(mesh->comm(), dim, mesh->debug());
    1527:         const Obj<typename Mesh::sieve_type>& sieve   = mesh->getSieve();
    1528:         PetscErrorCode                        ierr;
    1529:         ::tetgenio in;
    1530:         ::tetgenio out;

    1532:         const Obj<typename Mesh::label_sequence>&    vertices    = mesh->depthStratum(0);
    1533:         const Obj<typename Mesh::label_type>&        markers     = mesh->getLabel("marker");
    1534:         const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
    1535:         const Obj<typename Mesh::numbering_type>&    vNumbering  = mesh->getFactory()->getLocalNumbering(mesh, 0);

    1537:         in.numberofpoints = vertices->size();
    1538:         if (in.numberofpoints > 0) {
    1539:           const typename Mesh::label_sequence::iterator vEnd = vertices->end();

    1541:           in.pointlist       = new double[in.numberofpoints*dim];
    1542:           in.pointmarkerlist = new int[in.numberofpoints];
    1543:           for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
    1544:             const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
    1545:             const int                                           idx   = vNumbering->getIndex(*v_iter);

    1547:             for(int d = 0; d < dim; d++) {
    1548:               in.pointlist[idx*dim + d] = array[d];
    1549:             }
    1550:             in.pointmarkerlist[idx] = mesh->getValue(markers, *v_iter);
    1551:           }
    1552:         }
    1553:         const Obj<typename Mesh::label_sequence>& cells      = mesh->heightStratum(0);
    1554:         const Obj<typename Mesh::numbering_type>& cNumbering = mesh->getFactory()->getLocalNumbering(mesh, depth);

    1556:         in.numberofcorners       = 4;
    1557:         in.numberoftetrahedra    = cells->size();
    1558:         in.tetrahedronvolumelist = (double *) maxVolumes;
    1559:         if (in.numberoftetrahedra > 0) {
    1560:           in.tetrahedronlist     = new int[in.numberoftetrahedra*in.numberofcorners];
    1561:           if (mesh->depth() == 1) {
    1562:             ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(4);
    1563:             const typename Mesh::label_sequence::iterator cEnd = cells->end();

    1565:             for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
    1566:               sieve->cone(*c_iter, pV);
    1567:               const int         idx  = cNumbering->getIndex(*c_iter);
    1568:               const size_t      n    = pV.getSize();
    1569:               const point_type *cone = pV.getPoints();

    1571:               assert(n == 4);
    1572:               for(int v = 0; v < 4; ++v) {
    1573:                 in.tetrahedronlist[idx*in.numberofcorners + v] = vNumbering->getIndex(cone[v]);
    1574:               }
    1575:               pV.clear();
    1576:             }
    1577:           } else if (mesh->depth() == 3) {
    1578:             // Need extra space due to early error checking
    1579:             ALE::ISieveVisitor::NConeRetriever<typename Mesh::sieve_type> ncV(*sieve, 5);
    1580:             const typename Mesh::label_sequence::iterator cEnd = cells->end();

    1582:             for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
    1583:               ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
    1584:               const int         idx  = cNumbering->getIndex(*c_iter);
    1585:               const size_t      n    = ncV.getSize();
    1586:               const point_type *cone = ncV.getPoints();

    1588:               assert(n == 4);
    1589:               for(int v = 0; v < 4; ++v) {
    1590:                 in.tetrahedronlist[idx*in.numberofcorners + v] = vNumbering->getIndex(cone[v]);
    1591:               }
    1592:               ncV.clear();
    1593:             }
    1594:           } else {
    1595:             throw ALE::Exception("Invalid sieve: Cannot gives sieves of arbitrary depth to TetGen");
    1596:           }
    1597:         }
    1598:         if (depth == 3) {
    1599:           const Obj<typename Mesh::label_sequence>&     boundary = mesh->getLabelStratum("marker", 1);
    1600:           const typename Mesh::label_sequence::iterator bEnd     = boundary->end();

    1602:           in.numberoftrifaces = 0;
    1603:           for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != bEnd; ++b_iter) {
    1604:             if (mesh->height(*b_iter) == 1) {
    1605:               in.numberoftrifaces++;
    1606:             }
    1607:           }
    1608:           if (in.numberoftrifaces > 0) {
    1609:             ALE::ISieveVisitor::NConeRetriever<typename Mesh::sieve_type> ncV(*sieve, 5);
    1610:             int f = 0;

    1612:             in.trifacelist       = new int[in.numberoftrifaces*3];
    1613:             in.trifacemarkerlist = new int[in.numberoftrifaces];
    1614:             for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != bEnd; ++b_iter) {
    1615:               if (mesh->height(*b_iter) == 1) {
    1616:                 ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*sieve, *b_iter, ncV);
    1617:                 const size_t      n    = ncV.getSize();
    1618:                 const point_type *cone = ncV.getPoints();

    1620:                 for(size_t p = 0; p < n; ++p) {
    1621:                   in.trifacelist[f*3 + p] = vNumbering->getIndex(cone[p]);
    1622:                 }
    1623:                 in.trifacemarkerlist[f++] = mesh->getValue(markers, *b_iter);
    1624:                 ncV.clear();
    1625:               }
    1626:             }
    1627:           }
    1628:         }
    1629:         const typename Mesh::holes_type& holes = mesh->getHoles();

    1631:         in.numberofholes = holes.size();
    1632:         if (in.numberofholes > 0) {
    1633:           PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);CHKERRXX(ierr);
    1634:           for(int h = 0; h < in.numberofholes; ++h) {
    1635:             for(int d = 0; d < dim; ++d) {
    1636:               in.holelist[h*dim+d] = holes[h][d];
    1637:             }
    1638:           }
    1639:         }
    1640:         if (mesh->commRank() == 0) {
    1641:           std::string args("qezQra");

    1643:           ::tetrahedralize((char *) args.c_str(), &in, &out);
    1644:         }
    1645:         in.tetrahedronvolumelist = NULL;
    1646:         const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(refMesh->comm(), refMesh->debug());
    1647:         const Obj<FlexMesh>                  m        = new FlexMesh(mesh->comm(), dim, mesh->debug());
    1648:         const Obj<FlexMesh::sieve_type>      newS     = new FlexMesh::sieve_type(m->comm(), m->debug());
    1649:         int     numCorners  = 4;
    1650:         int     numCells    = out.numberoftetrahedra;
    1651:         int    *newCells    = out.tetrahedronlist;
    1652:         int     numVertices = out.numberofpoints;
    1653:         double *coords      = out.pointlist;
    1654:         real   *coordsR;

    1656:         if (!interpolate) {
    1657:           for(int c = 0; c < numCells; ++c) {
    1658:             int tmp = newCells[c*4+0];
    1659:             newCells[c*4+0] = newCells[c*4+1];
    1660:             newCells[c*4+1] = tmp;
    1661:           }
    1662:         }
    1663:         ALE::SieveBuilder<FlexMesh>::buildTopology(newS, dim, numCells, newCells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
    1664:         m->setSieve(newS);
    1665:         m->stratify();
    1666:         refMesh->setSieve(newSieve);
    1667:         std::map<typename Mesh::point_type,typename Mesh::point_type> renumbering;
    1668:         ALE::ISieveConverter::convertSieve(*newS, *newSieve, renumbering, renumber);
    1669:         refMesh->stratify();
    1670:         ALE::ISieveConverter::convertOrientation(*newS, *newSieve, renumbering, m->getArrowSection("orientation").ptr());
    1671:         {
    1672:           if (sizeof(double) == sizeof(real)) {
    1673:             coordsR = (real *) coords;
    1674:           } else {
    1675:             coordsR = new real[numVertices*dim];
    1676:             for(int i = 0; i < numVertices*dim; ++i) coordsR[i] = coords[i];
    1677:           }
    1678:         }
    1679:         ALE::SieveBuilder<Mesh>::buildCoordinates(refMesh, dim, coordsR);
    1680:         {
    1681:           if (sizeof(double) != sizeof(real)) {
    1682:             delete [] coordsR;
    1683:           }
    1684:         }
    1685:         const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");

    1687:         for(int v = 0; v < out.numberofpoints; v++) {
    1688:           if (out.pointmarkerlist[v]) {
    1689:             if (renumber) {
    1690:               refMesh->setValue(newMarkers, renumbering[v+out.numberoftetrahedra], out.pointmarkerlist[v]);
    1691:             } else {
    1692:               refMesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
    1693:             }
    1694:           }
    1695:         }
    1696:         if (interpolate) {
    1697:           // This does not work anymore (edgemarkerlist is always empty). I tried -ee and it gave bogus results
    1698:           if (out.edgemarkerlist) {
    1699:             for(int e = 0; e < out.numberofedges; e++) {
    1700:               if (out.edgemarkerlist[e]) {
    1701:                 typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
    1702:                 typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
    1703:                 Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(endpointA, endpointB, 1);

    1705:                 if (renumber) {
    1706:                   refMesh->setValue(newMarkers, renumbering[*edge->begin()], out.edgemarkerlist[e]);
    1707:                 } else {
    1708:                   refMesh->setValue(newMarkers, *edge->begin(), out.edgemarkerlist[e]);
    1709:                 }
    1710:               }
    1711:             }
    1712:           }
    1713:           if (out.trifacemarkerlist) {
    1714:             for(int f = 0; f < out.numberoftrifaces; f++) {
    1715:               if (out.trifacemarkerlist[f]) {
    1716:                 typename Mesh::point_type cornerA(out.trifacelist[f*3+0]+out.numberoftetrahedra);
    1717:                 typename Mesh::point_type cornerB(out.trifacelist[f*3+1]+out.numberoftetrahedra);
    1718:                 typename Mesh::point_type cornerC(out.trifacelist[f*3+2]+out.numberoftetrahedra);
    1719:                 Obj<typename Mesh::sieve_type::supportSet> corners = typename Mesh::sieve_type::supportSet();
    1720:                 Obj<typename Mesh::sieve_type::supportSet> edges   = typename Mesh::sieve_type::supportSet();
    1721:                 corners->clear();corners->insert(cornerA);corners->insert(cornerB);
    1722:                 edges->insert(*newS->nJoin1(corners)->begin());
    1723:                 corners->clear();corners->insert(cornerB);corners->insert(cornerC);
    1724:                 edges->insert(*newS->nJoin1(corners)->begin());
    1725:                 corners->clear();corners->insert(cornerC);corners->insert(cornerA);
    1726:                 edges->insert(*newS->nJoin1(corners)->begin());
    1727:                 ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(30);
    1728:                 typename Mesh::point_type face       = *newS->nJoin1(edges)->begin();
    1729:                 const int                 faceMarker = out.trifacemarkerlist[f];

    1731:                 if (renumber) {face = renumbering[face];}
    1732:                 ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*newSieve, face, pV);
    1733:                 const size_t                     n    = pV.getSize();
    1734:                 const typename Mesh::point_type *cone = pV.getPoints();

    1736:                 for(size_t c = 0; c < n; ++c) {
    1737:                   refMesh->setValue(newMarkers, cone[c], faceMarker);
    1738:                 }
    1739:                 pV.clear();
    1740:               }
    1741:             }
    1742:           }
    1743:         }
    1744:         return refMesh;
    1745:       };
    1746:       static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false, const bool forceSerial = false, const bool renumber = false) {
    1747:         throw ALE::Exception("Not yet implemented");
    1748:       };
    1749:       static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false, const bool forceSerial = false, const bool renumber = false) {
    1750:         const int numCells         = mesh->heightStratum(0)->size();
    1751:         double   *serialMaxVolumes = new double[numCells];

    1753:         for(int c = 0; c < numCells; c++) {
    1754:           serialMaxVolumes[c] = maxVolume;
    1755:         }
    1756:         const Obj<Mesh> refMesh = refineMeshV(mesh, serialMaxVolumes, interpolate, forceSerial, renumber);
    1757:         delete [] serialMaxVolumes;
    1758:         return refMesh;
    1759:       };
    1760:     };
    1761:   };
    1762: #endif
    1763:   template<typename Mesh>
    1764:   class Generator {
    1765:   public:
    1766:     static Obj<Mesh> generateMesh(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false) {
    1767:       int dim = boundary->getDimension();

    1769:       if (dim == 1) {
    1770: #ifdef PETSC_HAVE_TRIANGLE
    1771:         return ALE::Triangle::Generator<Mesh>::generateMesh(boundary, interpolate, constrained);
    1772: #else
    1773:         throw ALE::Exception("Mesh generation currently requires Triangle to be installed. Use --download-triangle during configure.");
    1774: #endif
    1775:       } else if (dim == 2) {
    1776: #ifdef PETSC_HAVE_TETGEN
    1777:         return ALE::TetGen::Generator<Mesh>::generateMesh(boundary, interpolate, constrained);
    1778: #else
    1779:         throw ALE::Exception("Mesh generation currently requires TetGen to be installed. Use --download-tetgen during configure.");
    1780: #endif
    1781:       }
    1782:       return NULL;
    1783:     };
    1784:     static Obj<Mesh> generateMeshV(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false, const bool renumber = false) {
    1785:       int dim = boundary->getDimension();

    1787:       if (dim == 1) {
    1788: #ifdef PETSC_HAVE_TRIANGLE
    1789:         return ALE::Triangle::Generator<Mesh>::generateMeshV(boundary, interpolate, constrained, renumber);
    1790: #else
    1791:         throw ALE::Exception("Mesh generation currently requires Triangle to be installed. Use --download-triangle during configure.");
    1792: #endif
    1793:       } else if (dim == 2) {
    1794: #ifdef PETSC_HAVE_TETGEN
    1795:         return ALE::TetGen::Generator<Mesh>::generateMeshV(boundary, interpolate, constrained, renumber);
    1796: #else
    1797:         throw ALE::Exception("Mesh generation currently requires TetGen to be installed. Use --download-tetgen during configure.");
    1798: #endif
    1799:       }
    1800:       return NULL;
    1801:     };
    1802:     static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false) {
    1803:       int dim = mesh->getDimension();

    1805:       if (dim == 2) {
    1806: #ifdef PETSC_HAVE_TRIANGLE
    1807:         return ALE::Triangle::Refiner<Mesh>::refineMesh(mesh, maxVolumes, interpolate);
    1808: #else
    1809:         throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
    1810: #endif
    1811:       } else if (dim == 3) {
    1812: #ifdef PETSC_HAVE_TETGEN
    1813:         return ALE::TetGen::Refiner<Mesh>::refineMesh(mesh, maxVolumes, interpolate);
    1814: #else
    1815:         throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
    1816: #endif
    1817:       }
    1818:       return NULL;
    1819:     };
    1820:     static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false, const bool forceSerial = false) {
    1821:       int dim = mesh->getDimension();

    1823:       if (dim == 2) {
    1824: #ifdef PETSC_HAVE_TRIANGLE
    1825:         return ALE::Triangle::Refiner<Mesh>::refineMesh(mesh, maxVolume, interpolate, forceSerial);
    1826: #else
    1827:         throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
    1828: #endif
    1829:       } else if (dim == 3) {
    1830: #ifdef PETSC_HAVE_TETGEN
    1831:         return ALE::TetGen::Refiner<Mesh>::refineMesh(mesh, maxVolume, interpolate);
    1832: #else
    1833:         throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
    1834: #endif
    1835:       }
    1836:       return NULL;
    1837:     };
    1838:     static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false, const bool renumber = false) {
    1839:       int dim = mesh->getDimension();

    1841:       if (dim == 2) {
    1842: #ifdef PETSC_HAVE_TRIANGLE
    1843:         return ALE::Triangle::Refiner<Mesh>::refineMeshV(mesh, maxVolumes, interpolate, renumber);
    1844: #else
    1845:         throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
    1846: #endif
    1847:       } else if (dim == 3) {
    1848: #ifdef PETSC_HAVE_TETGEN
    1849:         return ALE::TetGen::Refiner<Mesh>::refineMeshV(mesh, maxVolumes, interpolate, renumber);
    1850: #else
    1851:         throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
    1852: #endif
    1853:       }
    1854:       return NULL;
    1855:     };
    1856:     static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false, const bool forceSerial = false, const bool renumber = false) {
    1857:       int dim = mesh->getDimension();

    1859:       if (dim == 2) {
    1860: #ifdef PETSC_HAVE_TRIANGLE
    1861:         return ALE::Triangle::Refiner<Mesh>::refineMeshV(mesh, maxVolume, interpolate, forceSerial, renumber);
    1862: #else
    1863:         throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
    1864: #endif
    1865:       } else if (dim == 3) {
    1866: #ifdef PETSC_HAVE_TETGEN
    1867:         return ALE::TetGen::Refiner<Mesh>::refineMeshV(mesh, maxVolume, interpolate, forceSerial, renumber);
    1868: #else
    1869:         throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
    1870: #endif
    1871:       }
    1872:       return NULL;
    1873:     };
    1874:     static Obj<Mesh> refineMeshLocal(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false) {
    1875:       int dim = mesh->getDimension();

    1877:       if (dim == 2) {
    1878: #ifdef PETSC_HAVE_TRIANGLE
    1879:         return ALE::Triangle::Refiner<Mesh>::refineMeshLocal(mesh, maxVolume, interpolate);
    1880: #else
    1881:         throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
    1882: #endif
    1883:       } else if (dim == 3) {
    1884: #ifdef PETSC_HAVE_TETGEN
    1885:         return ALE::TetGen::Refiner<Mesh>::refineMesh(mesh, maxVolume, interpolate);
    1886: #else
    1887:         throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
    1888: #endif
    1889:       }
    1890:       return NULL;
    1891:     };
    1892:   };
    1893: }

    1895: #endif