Actual source code: Generator.hh

  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) {
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, false);
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:               mesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
339:             }
340:           }
341:           if (interpolate) {
342:             for(int e = 0; e < out.numberofedges; e++) {
343:               if (out.edgemarkerlist[e]) {
344:                 const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
345:                 const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
346:                 const Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(vertexA, vertexB, 1);

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

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

379:         in.numberofpoints = vertices->size();
380:         if (in.numberofpoints > 0) {
381:           const typename Mesh::label_sequence::iterator vEnd = vertices->end();

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

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

398:         in.numberofcorners   = 3;
399:         in.numberoftriangles = faces->size();
400:         in.trianglearealist  = (double *) maxVolumes;
401:         if (in.numberoftriangles > 0) {
402:           const typename Mesh::label_sequence::iterator fEnd = faces->end();

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

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

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

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

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

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

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

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

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

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

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

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

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

534:                 refMesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
535:               }
536:             }
537:           }
538:         }

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

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

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

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

587:         in.numberofpoints = vertices->size();
588:         if (in.numberofpoints > 0) {
589:           const typename Mesh::label_sequence::iterator vEnd = vertices->end();

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

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

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

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

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

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

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

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

665:               if (edgeMarker) {
666:                 sieve->cone(*e_iter, pV);
667:                 const size_t      n    = pV.getSize();
668:                 const point_type *cone = pV.getPoints();

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

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

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

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

735:         if (refMesh->commRank() == 0) {
736:           for(int v = 0; v < out.numberofpoints; v++) {
737:             if (out.pointmarkerlist[v]) {
738:               refMesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
739:             }
740:           }
741:           if (interpolate) {
742:             for(int e = 0; e < out.numberofedges; e++) {
743:               if (out.edgemarkerlist[e]) {
744:                 const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
745:                 const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
746:                 const Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(vertexA, vertexB, 1);

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

777:         for(int c = 0; c < numCells; c++) {
778:           serialMaxVolumes[c] = maxVolume;
779:         }
780:         const Obj<Mesh> refMesh = refineMeshV(mesh, serialMaxVolumes, interpolate);
781:         delete [] serialMaxVolumes;
782:         return refMesh;
783:       };
784:       static Obj<Mesh> refineMeshLocal(const Obj<Mesh>& mesh, const double maxVolumes[], const bool interpolate = false) {
785:         const int                    dim     = mesh->getDimension();
786:         const Obj<Mesh>              refMesh = new Mesh(mesh->comm(), dim, mesh->debug());
787:         const Obj<typename Mesh::sieve_type>& sieve   = mesh->getSieve();
788:         struct triangulateio in;
789:         struct triangulateio out;
790:         PetscErrorCode       ierr;

792:         Generator<Mesh>::initInput(&in);
793:         Generator<Mesh>::initOutput(&out);
794:         const Obj<typename Mesh::label_sequence>&    vertices    = mesh->depthStratum(0);
795:         const Obj<typename Mesh::label_type>&        markers     = mesh->getLabel("marker");
796:         const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
797:         const Obj<typename Mesh::numbering_type>&    vNumbering  = mesh->getFactory()->getLocalNumbering(mesh, 0);

799:         in.numberofpoints = vertices->size();
800:         if (in.numberofpoints > 0) {
801:           PetscMalloc(in.numberofpoints * dim * sizeof(double), &in.pointlist);CHKERRXX(ierr);
802:           PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);CHKERRXX(ierr);
803:           for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
804:             const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
805:             const int                                  idx   = vNumbering->getIndex(*v_iter);

807:             for(int d = 0; d < dim; d++) {
808:               in.pointlist[idx*dim + d] = array[d];
809:             }
810:             in.pointmarkerlist[idx] = mesh->getValue(markers, *v_iter);
811:           }
812:         }
813:         const Obj<typename Mesh::label_sequence>& faces      = mesh->heightStratum(0);
814:         const Obj<typename Mesh::numbering_type>& fNumbering = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());

816:         in.numberofcorners   = 3;
817:         in.numberoftriangles = faces->size();
818:         in.trianglearealist  = (double *) maxVolumes;
819:         if (in.numberoftriangles > 0) {
820:           PetscMalloc(in.numberoftriangles*in.numberofcorners * sizeof(int), &in.trianglelist);CHKERRXX(ierr);
821:           if (mesh->depth() == 1) {
822:             for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != faces->end(); ++f_iter) {
823:               const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = sieve->cone(*f_iter);
824:               const int                                          idx  = fNumbering->getIndex(*f_iter);
825:               int                                                v    = 0;

827:               for(typename Mesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
828:                 in.trianglelist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*c_iter);
829:               }
830:             }
831:           } else if (mesh->depth() == 2) {
832:             for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != faces->end(); ++f_iter) {
833:               typedef ALE::SieveAlg<Mesh> sieve_alg_type;
834:               const Obj<typename sieve_alg_type::coneArray>& cone = sieve_alg_type::nCone(mesh, *f_iter, 2);
835:               const int                             idx  = fNumbering->getIndex(*f_iter);
836:               int                                   v    = 0;

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

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

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

866:                 for(typename Mesh::sieve_type::traits::coneSequence::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
867:                   in.segmentlist[s*2 + (p++)] = vNumbering->getIndex(*v_iter);
868:                 }
869:                 in.segmentmarkerlist[s++] = edgeMarker;
870:               }
871:             }
872:           }
873:         }

875:         in.numberofholes = 0;
876:         if (in.numberofholes > 0) {
877:           PetscMalloc(in.numberofholes * dim * sizeof(int), &in.holelist);CHKERRXX(ierr);
878:         }
879:         std::string args("pqezQra");

881:         triangulate((char *) args.c_str(), &in, &out, NULL);
882:         if (in.pointlist)         {PetscFree(in.pointlist);CHKERRXX(ierr);}
883:         if (in.pointmarkerlist)   {PetscFree(in.pointmarkerlist);CHKERRXX(ierr);}
884:         if (in.segmentlist)       {PetscFree(in.segmentlist);CHKERRXX(ierr);}
885:         if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);CHKERRXX(ierr);}
886:         if (in.trianglelist)      {PetscFree(in.trianglelist);CHKERRXX(ierr);}
887:         const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
888:         int     numCorners  = 3;
889:         int     numCells    = out.numberoftriangles;
890:         int    *cells       = out.trianglelist;
891:         int     numVertices = out.numberofpoints;
892:         double *coords      = out.pointlist;

894:         ALE::SieveBuilder<Mesh>::buildTopologyMultiple(newSieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, refMesh->getArrowSection("orientation"));
895:         refMesh->setSieve(newSieve);
896:         refMesh->stratify();
897:         ALE::SieveBuilder<Mesh>::buildCoordinatesMultiple(refMesh, dim, coords);
898:         const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");

900:         for(int v = 0; v < out.numberofpoints; v++) {
901:           if (out.pointmarkerlist[v]) {
902:             refMesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
903:           }
904:         }
905:         if (interpolate) {
906:           for(int e = 0; e < out.numberofedges; e++) {
907:             if (out.edgemarkerlist[e]) {
908:               const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
909:               const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
910:               const Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(vertexA, vertexB, 1);

912:               refMesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
913:             }
914:           }
915:         }
916:         Generator<Mesh>::finiOutput(&out);
917:         return refMesh;
918:       };
919:       static Obj<Mesh> refineMeshLocal(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false) {
920:         const int numLocalFaces   = mesh->heightStratum(0)->size();
921:         double   *localMaxVolumes = new double[numLocalFaces];

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

963:         const Obj<typename Mesh::label_sequence>&    vertices    = boundary->depthStratum(0);
964:         const Obj<typename Mesh::numbering_type>&    vNumbering  = boundary->getFactory()->getLocalNumbering(boundary, 0);
965:         const Obj<typename Mesh::real_section_type>& coordinates = boundary->getRealSection("coordinates");
966:         const Obj<typename Mesh::label_type>&        markers     = boundary->getLabel("marker");


969:         in.numberofpoints = vertices->size();
970:         if (in.numberofpoints > 0) {
971:           in.pointlist       = new double[in.numberofpoints*dim];
972:           in.pointmarkerlist = new int[in.numberofpoints];
973:           for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
974:             const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
975:             const int                                  idx   = vNumbering->getIndex(*v_iter);

977:             for(int d = 0; d < dim; d++) {
978:               in.pointlist[idx*dim + d] = array[d];
979:             }
980:             in.pointmarkerlist[idx] = boundary->getValue(markers, *v_iter);
981:           }
982:         }

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

989:           in.numberoffacets = facets->size();
990:           if (in.numberoffacets > 0) {
991:             in.facetlist       = new tetgenio::facet[in.numberoffacets];
992:             in.facetmarkerlist = new int[in.numberoffacets];
993:             for(typename Mesh::label_sequence::iterator f_iter = facets->begin(); f_iter != facets->end(); ++f_iter) {
994:               const Obj<typename sieve_alg_type::coneArray>& cone = sieve_alg_type::nCone(boundary, *f_iter, boundary->depth());
995:               const int                             idx  = fNumbering->getIndex(*f_iter);

997:               in.facetlist[idx].numberofpolygons = 1;
998:               in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
999:               in.facetlist[idx].numberofholes    = 0;
1000:               in.facetlist[idx].holelist         = NULL;

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

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

1010:                 poly->vertexlist[c++] = vIdx;
1011:               }
1012:               in.facetmarkerlist[idx] = boundary->getValue(markers, *f_iter);
1013:             }
1014:           }
1015:         }else {
1016:           createConvexHull = true;
1017:         }

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

1041:           //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.
1042:           ::tetrahedralize((char *) args.c_str(), &in, &out);
1043:         }
1044:         const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
1045:         int     numCorners  = 4;
1046:         int     numCells    = out.numberoftetrahedra;
1047:         int    *cells       = out.tetrahedronlist;
1048:         int     numVertices = out.numberofpoints;
1049:         double *coords      = out.pointlist;

1051:         if (!interpolate) {
1052:           for(int c = 0; c < numCells; ++c) {
1053:             int tmp = cells[c*4+0];
1054:             cells[c*4+0] = cells[c*4+1];
1055:             cells[c*4+1] = tmp;
1056:           }
1057:         }
1058:         ALE::SieveBuilder<Mesh>::buildTopology(newSieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, mesh->getArrowSection("orientation"));
1059:         mesh->setSieve(newSieve);
1060:         mesh->stratify();
1061:         ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, dim, coords);
1062:         const Obj<typename Mesh::label_type>& newMarkers = mesh->createLabel("marker");

1064:         for(int v = 0; v < out.numberofpoints; v++) {
1065:           if (out.pointmarkerlist[v]) {
1066:             mesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
1067:           }
1068:         }
1069:         if (interpolate) {
1070:           if (out.edgemarkerlist) {
1071:             for(int e = 0; e < out.numberofedges; e++) {
1072:               if (out.edgemarkerlist[e]) {
1073:                 typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
1074:                 typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
1075:                 Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(endpointA, endpointB, 1);

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

1109:                 for(typename Mesh::coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
1110:                   mesh->setValue(newMarkers, *cl_iter, faceMarker);
1111:                 }
1112:               }
1113:             }
1114:           }
1115:         }
1116:         return mesh;
1117:       };
1120:       static Obj<Mesh> generateMeshV(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false) {
1121:         typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
1122:         typedef typename Mesh::real_section_type::value_type real;
1123:         const int                             dim   = 3;
1124:         Obj<Mesh>                             mesh  = new Mesh(boundary->comm(), dim, boundary->debug());
1125:         const Obj<typename Mesh::sieve_type>& sieve = boundary->getSieve();
1126:         const PetscMPIInt                     rank  = mesh->commRank();
1127:         bool                                  createConvexHull = false;
1128:         PetscErrorCode                        ierr;
1129:         ::tetgenio in;
1130:         ::tetgenio out;

1132:         const Obj<typename Mesh::label_sequence>&    vertices    = boundary->depthStratum(0);
1133:         const Obj<typename Mesh::label_type>&        markers     = boundary->getLabel("marker");
1134:         const Obj<typename Mesh::real_section_type>& coordinates = boundary->getRealSection("coordinates");
1135:         const Obj<typename Mesh::numbering_type>&    vNumbering  = boundary->getFactory()->getLocalNumbering(boundary, 0);

1137:         in.numberofpoints = vertices->size();
1138:         if (in.numberofpoints > 0) {
1139:           const typename Mesh::label_sequence::iterator vEnd = vertices->end();

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

1147:             for(int d = 0; d < dim; ++d) {
1148:               in.pointlist[idx*dim + d] = array[d];
1149:             }
1150:             in.pointmarkerlist[idx] = boundary->getValue(markers, *v_iter);
1151:           }
1152:         }

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

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

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

1170:               in.facetlist[idx].numberofpolygons = 1;
1171:               in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
1172:               in.facetlist[idx].numberofholes    = 0;
1173:               in.facetlist[idx].holelist         = NULL;

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

1179:               poly->numberofvertices = n;
1180:               poly->vertexlist       = new int[poly->numberofvertices];
1181:               for(size_t c = 0; c < n; ++c) {
1182:                 poly->vertexlist[c] = vNumbering->getIndex(cone[c]);
1183:               }
1184:               in.facetmarkerlist[idx] = boundary->getValue(markers, *f_iter);
1185:               ncV.clear();
1186:             }
1187:           }
1188:         } else {
1189:           createConvexHull = true;
1190:         }
1191:         const typename Mesh::holes_type& holes = mesh->getHoles();

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

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

1265:         for(int v = 0; v < out.numberofpoints; v++) {
1266:           if (out.pointmarkerlist[v]) {
1267:             mesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
1268:           }
1269:         }
1270:         if (interpolate) {
1271:           if (out.edgemarkerlist) {
1272:             for(int e = 0; e < out.numberofedges; e++) {
1273:               if (out.edgemarkerlist[e]) {
1274:                 typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
1275:                 typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
1276:                 Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(endpointA, endpointB, 1);

1278:                 mesh->setValue(newMarkers, *edge->begin(), out.edgemarkerlist[e]);
1279:               }
1280:             }
1281:           }
1282:           if (out.trifacemarkerlist) {
1283:             // Work around TetGen bug for raw tetrahedralization
1284:             //   The boundary faces are 0,1,4,5,8,9,11,12,13,15,16,17
1285: //             for(int f = 0; f < out.numberoftrifaces; f++) {
1286: //               if (out.trifacemarkerlist[f]) {
1287: //                 out.trifacemarkerlist[f] = 0;
1288: //               } else {
1289: //                 out.trifacemarkerlist[f] = 1;
1290: //               }
1291: //             }
1292:             for(int f = 0; f < out.numberoftrifaces; f++) {
1293:               if (out.trifacemarkerlist[f]) {
1294:                 typename Mesh::point_type cornerA(out.trifacelist[f*3+0]+out.numberoftetrahedra);
1295:                 typename Mesh::point_type cornerB(out.trifacelist[f*3+1]+out.numberoftetrahedra);
1296:                 typename Mesh::point_type cornerC(out.trifacelist[f*3+2]+out.numberoftetrahedra);
1297:                 Obj<typename Mesh::sieve_type::supportSet> corners = typename Mesh::sieve_type::supportSet();
1298:                 Obj<typename Mesh::sieve_type::supportSet> edges   = typename Mesh::sieve_type::supportSet();
1299:                 corners->clear();corners->insert(cornerA);corners->insert(cornerB);
1300:                 edges->insert(*newS->nJoin1(corners)->begin());
1301:                 corners->clear();corners->insert(cornerB);corners->insert(cornerC);
1302:                 edges->insert(*newS->nJoin1(corners)->begin());
1303:                 corners->clear();corners->insert(cornerC);corners->insert(cornerA);
1304:                 edges->insert(*newS->nJoin1(corners)->begin());
1305:                 ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(30);
1306:                 const typename Mesh::point_type           face       = *newS->nJoin1(edges)->begin();
1307:                 const int                                 faceMarker = out.trifacemarkerlist[f];

1309:                 ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*newSieve, face, pV);
1310:                 const size_t                     n    = pV.getSize();
1311:                 const typename Mesh::point_type *cone = pV.getPoints();

1313:                 for(size_t c = 0; c < n; ++c) {
1314:                   mesh->setValue(newMarkers, cone[c], faceMarker);
1315:                 }
1316:                 pV.clear();
1317:               }
1318:             }
1319:           }
1320:         }
1321:         mesh->copyHoles(boundary);
1322:         return mesh;
1323:       };
1324:     };
1325:     template<typename Mesh>
1326:     class Refiner {
1327:     public:
1328:       static Obj<Mesh> refineMesh(const Obj<Mesh>& serialMesh, const double maxVolumes[], const bool interpolate = false) {
1329:         typedef ALE::SieveAlg<Mesh> sieve_alg_type;
1330:         const int       dim     = serialMesh->getDimension();
1331:         const int       depth   = serialMesh->depth();
1332:         const Obj<Mesh> refMesh = new Mesh(serialMesh->comm(), dim, serialMesh->debug());
1333:         ::tetgenio      in;
1334:         ::tetgenio      out;

1336:         const Obj<typename Mesh::label_sequence>&    vertices    = serialMesh->depthStratum(0);
1337:         const Obj<typename Mesh::label_type>&        markers     = serialMesh->getLabel("marker");
1338:         const Obj<typename Mesh::real_section_type>& coordinates = serialMesh->getRealSection("coordinates");
1339:         const Obj<typename Mesh::numbering_type>&    vNumbering  = serialMesh->getFactory()->getLocalNumbering(serialMesh, 0);

1341:         in.numberofpoints = vertices->size();
1342:         if (in.numberofpoints > 0) {
1343:           in.pointlist       = new double[in.numberofpoints*dim];
1344:           in.pointmarkerlist = new int[in.numberofpoints];
1345:           for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
1346:             const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
1347:             const int                                  idx   = vNumbering->getIndex(*v_iter);

1349:             for(int d = 0; d < dim; d++) {
1350:               in.pointlist[idx*dim + d] = array[d];
1351:             }
1352:             in.pointmarkerlist[idx] = serialMesh->getValue(markers, *v_iter);
1353:           }
1354:         }
1355:         const Obj<typename Mesh::label_sequence>& cells      = serialMesh->heightStratum(0);
1356:         const Obj<typename Mesh::numbering_type>& cNumbering = serialMesh->getFactory()->getLocalNumbering(serialMesh, depth);

1358:         in.numberofcorners       = 4;
1359:         in.numberoftetrahedra    = cells->size();
1360:         in.tetrahedronvolumelist = (double *) maxVolumes;
1361:         if (in.numberoftetrahedra > 0) {
1362:           in.tetrahedronlist     = new int[in.numberoftetrahedra*in.numberofcorners];
1363:           for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
1364:             typedef ALE::SieveAlg<Mesh> sieve_alg_type;
1365:             const Obj<typename sieve_alg_type::coneArray>& cone = sieve_alg_type::nCone(serialMesh, *c_iter, depth);
1366:             const int                             idx  = cNumbering->getIndex(*c_iter);
1367:             int                                   v    = 0;

1369:             for(typename Mesh::sieve_type::coneArray::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
1370:               in.tetrahedronlist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*v_iter);
1371:             }
1372:           }
1373:         }
1374:         if (serialMesh->depth() == 3) {
1375:           const Obj<typename Mesh::label_sequence>& boundary = serialMesh->getLabelStratum("marker", 1);

1377:           in.numberoftrifaces = 0;
1378:           for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != boundary->end(); ++b_iter) {
1379:             if (serialMesh->height(*b_iter) == 1) {
1380:               in.numberoftrifaces++;
1381:             }
1382:           }
1383:           if (in.numberoftrifaces > 0) {
1384:             int f = 0;

1386:             in.trifacelist       = new int[in.numberoftrifaces*3];
1387:             in.trifacemarkerlist = new int[in.numberoftrifaces];
1388:             for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != boundary->end(); ++b_iter) {
1389:               if (serialMesh->height(*b_iter) == 1) {
1390:                 const Obj<typename Mesh::coneArray>& cone = sieve_alg_type::nCone(serialMesh, *b_iter, 2);
1391:                 int                         p    = 0;

1393:                 for(typename Mesh::coneArray::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
1394:                   in.trifacelist[f*3 + (p++)] = vNumbering->getIndex(*v_iter);
1395:                 }
1396:                 in.trifacemarkerlist[f++] = serialMesh->getValue(markers, *b_iter);
1397:               }
1398:             }
1399:           }
1400:         }

1402:         in.numberofholes = 0;
1403:         if (serialMesh->commRank() == 0) {
1404:           std::string args("qezQra");

1406:           ::tetrahedralize((char *) args.c_str(), &in, &out);
1407:         }
1408:         in.tetrahedronvolumelist = NULL;
1409:         const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(refMesh->comm(), refMesh->debug());
1410:         int     numCorners  = 4;
1411:         int     numCells    = out.numberoftetrahedra;
1412:         int    *newCells       = out.tetrahedronlist;
1413:         int     numVertices = out.numberofpoints;
1414:         double *coords      = out.pointlist;

1416:         if (!interpolate) {
1417:           for(int c = 0; c < numCells; ++c) {
1418:             int tmp = newCells[c*4+0];
1419:             newCells[c*4+0] = newCells[c*4+1];
1420:             newCells[c*4+1] = tmp;
1421:           }
1422:         }
1423:         ALE::SieveBuilder<Mesh>::buildTopology(newSieve, dim, numCells, newCells, numVertices, interpolate, numCorners, -1, refMesh->getArrowSection("orientation"));
1424:         refMesh->setSieve(newSieve);
1425:         refMesh->stratify();
1426:         ALE::SieveBuilder<Mesh>::buildCoordinates(refMesh, dim, coords);
1427:         const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");


1430:         for(int v = 0; v < out.numberofpoints; v++) {
1431:           if (out.pointmarkerlist[v]) {
1432:             refMesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
1433:           }
1434:         }
1435:         if (interpolate) {
1436:           if (out.edgemarkerlist) {
1437:             for(int e = 0; e < out.numberofedges; e++) {
1438:               if (out.edgemarkerlist[e]) {
1439:                 typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
1440:                 typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
1441:                 Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(endpointA, endpointB, 1);

1443:                 refMesh->setValue(newMarkers, *edge->begin(), out.edgemarkerlist[e]);
1444:               }
1445:             }
1446:           }
1447:           if (out.trifacemarkerlist) {
1448:             for(int f = 0; f < out.numberoftrifaces; f++) {
1449:               if (out.trifacemarkerlist[f]) {
1450:                 typename Mesh::point_type cornerA(out.trifacelist[f*3+0]+out.numberoftetrahedra);
1451:                 typename Mesh::point_type cornerB(out.trifacelist[f*3+1]+out.numberoftetrahedra);
1452:                 typename Mesh::point_type cornerC(out.trifacelist[f*3+2]+out.numberoftetrahedra);
1453:                 Obj<typename Mesh::sieve_type::supportSet> corners = typename Mesh::sieve_type::supportSet();
1454:                 Obj<typename Mesh::sieve_type::supportSet> edges   = typename Mesh::sieve_type::supportSet();
1455:                 corners->clear();corners->insert(cornerA);corners->insert(cornerB);
1456:                 edges->insert(*newSieve->nJoin1(corners)->begin());
1457:                 corners->clear();corners->insert(cornerB);corners->insert(cornerC);
1458:                 edges->insert(*newSieve->nJoin1(corners)->begin());
1459:                 corners->clear();corners->insert(cornerC);corners->insert(cornerA);
1460:                 edges->insert(*newSieve->nJoin1(corners)->begin());
1461:                 const typename Mesh::point_type          face       = *newSieve->nJoin1(edges)->begin();
1462:                 const int                       faceMarker = out.trifacemarkerlist[f];
1463:                 const Obj<typename Mesh::coneArray>      closure    = sieve_alg_type::closure(refMesh, face);
1464:                 const typename Mesh::coneArray::iterator end        = closure->end();

1466:                 for(typename Mesh::coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
1467:                   refMesh->setValue(newMarkers, *cl_iter, faceMarker);
1468:                 }
1469:               }
1470:             }
1471:           }
1472:         }
1473:         if (refMesh->commSize() > 1) {
1474:           return ALE::Distribution<Mesh>::distributeMesh(refMesh);
1475:         }
1476:         return refMesh;
1477:       };
1478:       static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false) {
1479:         Obj<Mesh>                          serialMesh       = ALE::Distribution<Mesh>::unifyMesh(mesh);
1480:         const Obj<typename Mesh::real_section_type> serialMaxVolumes = ALE::Distribution<Mesh>::distributeSection(maxVolumes, serialMesh, serialMesh->getDistSendOverlap(), serialMesh->getDistRecvOverlap());

1482:         return refineMesh(serialMesh, serialMaxVolumes->restrictSpace(), interpolate);
1483:       };
1484:       static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false) {
1485:         Obj<Mesh> serialMesh;
1486:         if (mesh->commSize() > 1) {
1487:           serialMesh = ALE::Distribution<Mesh>::unifyMesh(mesh);
1488:         } else {
1489:           serialMesh = mesh;
1490:         }
1491:         const int numCells         = serialMesh->heightStratum(0)->size();
1492:         double   *serialMaxVolumes = new double[numCells];

1494:         for(int c = 0; c < numCells; c++) {
1495:           serialMaxVolumes[c] = maxVolume;
1496:         }
1497:         const Obj<Mesh> refMesh = refineMesh(serialMesh, serialMaxVolumes, interpolate);
1498:         delete [] serialMaxVolumes;
1499:         return refMesh;
1500:       };
1501:       static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolumes[], const bool interpolate = false, const bool forceSerial = false) {
1502:         typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
1503:         typedef typename Mesh::real_section_type::value_type real;
1504:         typedef typename Mesh::point_type point_type;
1505:         const int                             dim     = mesh->getDimension();
1506:         const int                             depth   = mesh->depth();
1507:         const Obj<Mesh>                       refMesh = new Mesh(mesh->comm(), dim, mesh->debug());
1508:         const Obj<typename Mesh::sieve_type>& sieve   = mesh->getSieve();
1509:         PetscErrorCode                        ierr;
1510:         ::tetgenio in;
1511:         ::tetgenio out;

1513:         const Obj<typename Mesh::label_sequence>&    vertices    = mesh->depthStratum(0);
1514:         const Obj<typename Mesh::label_type>&        markers     = mesh->getLabel("marker");
1515:         const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
1516:         const Obj<typename Mesh::numbering_type>&    vNumbering  = mesh->getFactory()->getLocalNumbering(mesh, 0);

1518:         in.numberofpoints = vertices->size();
1519:         if (in.numberofpoints > 0) {
1520:           const typename Mesh::label_sequence::iterator vEnd = vertices->end();

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

1528:             for(int d = 0; d < dim; d++) {
1529:               in.pointlist[idx*dim + d] = array[d];
1530:             }
1531:             in.pointmarkerlist[idx] = mesh->getValue(markers, *v_iter);
1532:           }
1533:         }
1534:         const Obj<typename Mesh::label_sequence>& cells      = mesh->heightStratum(0);
1535:         const Obj<typename Mesh::numbering_type>& cNumbering = mesh->getFactory()->getLocalNumbering(mesh, depth);

1537:         in.numberofcorners       = 4;
1538:         in.numberoftetrahedra    = cells->size();
1539:         in.tetrahedronvolumelist = (double *) maxVolumes;
1540:         if (in.numberoftetrahedra > 0) {
1541:           in.tetrahedronlist     = new int[in.numberoftetrahedra*in.numberofcorners];
1542:           if (mesh->depth() == 1) {
1543:             ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(4);
1544:             const typename Mesh::label_sequence::iterator cEnd = cells->end();

1546:             for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
1547:               sieve->cone(*c_iter, pV);
1548:               const int         idx  = cNumbering->getIndex(*c_iter);
1549:               const size_t      n    = pV.getSize();
1550:               const point_type *cone = pV.getPoints();

1552:               assert(n == 4);
1553:               for(int v = 0; v < 4; ++v) {
1554:                 in.tetrahedronlist[idx*in.numberofcorners + v] = vNumbering->getIndex(cone[v]);
1555:               }
1556:               pV.clear();
1557:             }
1558:           } else if (mesh->depth() == 3) {
1559:             // Need extra space due to early error checking
1560:             ALE::ISieveVisitor::NConeRetriever<typename Mesh::sieve_type> ncV(*sieve, 5);
1561:             const typename Mesh::label_sequence::iterator cEnd = cells->end();

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

1569:               assert(n == 4);
1570:               for(int v = 0; v < 4; ++v) {
1571:                 in.tetrahedronlist[idx*in.numberofcorners + v] = vNumbering->getIndex(cone[v]);
1572:               }
1573:               ncV.clear();
1574:             }
1575:           } else {
1576:             throw ALE::Exception("Invalid sieve: Cannot gives sieves of arbitrary depth to TetGen");
1577:           }
1578:         }
1579:         if (depth == 3) {
1580:           const Obj<typename Mesh::label_sequence>&     boundary = mesh->getLabelStratum("marker", 1);
1581:           const typename Mesh::label_sequence::iterator bEnd     = boundary->end();

1583:           in.numberoftrifaces = 0;
1584:           for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != bEnd; ++b_iter) {
1585:             if (mesh->height(*b_iter) == 1) {
1586:               in.numberoftrifaces++;
1587:             }
1588:           }
1589:           if (in.numberoftrifaces > 0) {
1590:             ALE::ISieveVisitor::NConeRetriever<typename Mesh::sieve_type> ncV(*sieve, 5);
1591:             int f = 0;

1593:             in.trifacelist       = new int[in.numberoftrifaces*3];
1594:             in.trifacemarkerlist = new int[in.numberoftrifaces];
1595:             for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != bEnd; ++b_iter) {
1596:               if (mesh->height(*b_iter) == 1) {
1597:                 ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*sieve, *b_iter, ncV);
1598:                 const size_t      n    = ncV.getSize();
1599:                 const point_type *cone = ncV.getPoints();

1601:                 for(size_t p = 0; p < n; ++p) {
1602:                   in.trifacelist[f*3 + p] = vNumbering->getIndex(cone[p]);
1603:                 }
1604:                 in.trifacemarkerlist[f++] = mesh->getValue(markers, *b_iter);
1605:                 ncV.clear();
1606:               }
1607:             }
1608:           }
1609:         }
1610:         const typename Mesh::holes_type& holes = mesh->getHoles();

1612:         in.numberofholes = holes.size();
1613:         if (in.numberofholes > 0) {
1614:           PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);CHKERRXX(ierr);
1615:           for(int h = 0; h < in.numberofholes; ++h) {
1616:             for(int d = 0; d < dim; ++d) {
1617:               in.holelist[h*dim+d] = holes[h][d];
1618:             }
1619:           }
1620:         }
1621:         if (mesh->commRank() == 0) {
1622:           std::string args("qezQra");

1624:           ::tetrahedralize((char *) args.c_str(), &in, &out);
1625:         }
1626:         in.tetrahedronvolumelist = NULL;
1627:         const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(refMesh->comm(), refMesh->debug());
1628:         const Obj<FlexMesh>                  m        = new FlexMesh(mesh->comm(), dim, mesh->debug());
1629:         const Obj<FlexMesh::sieve_type>      newS     = new FlexMesh::sieve_type(m->comm(), m->debug());
1630:         int     numCorners  = 4;
1631:         int     numCells    = out.numberoftetrahedra;
1632:         int    *newCells    = out.tetrahedronlist;
1633:         int     numVertices = out.numberofpoints;
1634:         double *coords      = out.pointlist;
1635:         real   *coordsR;

1637:         if (!interpolate) {
1638:           for(int c = 0; c < numCells; ++c) {
1639:             int tmp = newCells[c*4+0];
1640:             newCells[c*4+0] = newCells[c*4+1];
1641:             newCells[c*4+1] = tmp;
1642:           }
1643:         }
1644:         ALE::SieveBuilder<FlexMesh>::buildTopology(newS, dim, numCells, newCells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
1645:         m->setSieve(newS);
1646:         m->stratify();
1647:         refMesh->setSieve(newSieve);
1648:         std::map<typename Mesh::point_type,typename Mesh::point_type> renumbering;
1649:         ALE::ISieveConverter::convertSieve(*newS, *newSieve, renumbering, false);
1650:         refMesh->stratify();
1651:         ALE::ISieveConverter::convertOrientation(*newS, *newSieve, renumbering, m->getArrowSection("orientation").ptr());
1652:         {
1653:           if (sizeof(double) == sizeof(real)) {
1654:             coordsR = (real *) coords;
1655:           } else {
1656:             coordsR = new real[numVertices*dim];
1657:             for(int i = 0; i < numVertices*dim; ++i) coordsR[i] = coords[i];
1658:           }
1659:         }
1660:         ALE::SieveBuilder<Mesh>::buildCoordinates(refMesh, dim, coordsR);
1661:         {
1662:           if (sizeof(double) != sizeof(real)) {
1663:             delete [] coordsR;
1664:           }
1665:         }
1666:         const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");

1668:         for(int v = 0; v < out.numberofpoints; v++) {
1669:           if (out.pointmarkerlist[v]) {
1670:             refMesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
1671:           }
1672:         }
1673:         if (interpolate) {
1674:           if (out.edgemarkerlist) {
1675:             for(int e = 0; e < out.numberofedges; e++) {
1676:               if (out.edgemarkerlist[e]) {
1677:                 typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
1678:                 typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
1679:                 Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(endpointA, endpointB, 1);

1681:                 refMesh->setValue(newMarkers, *edge->begin(), out.edgemarkerlist[e]);
1682:               }
1683:             }
1684:           }
1685:           if (out.trifacemarkerlist) {
1686:             for(int f = 0; f < out.numberoftrifaces; f++) {
1687:               if (out.trifacemarkerlist[f]) {
1688:                 typename Mesh::point_type cornerA(out.trifacelist[f*3+0]+out.numberoftetrahedra);
1689:                 typename Mesh::point_type cornerB(out.trifacelist[f*3+1]+out.numberoftetrahedra);
1690:                 typename Mesh::point_type cornerC(out.trifacelist[f*3+2]+out.numberoftetrahedra);
1691:                 Obj<typename Mesh::sieve_type::supportSet> corners = typename Mesh::sieve_type::supportSet();
1692:                 Obj<typename Mesh::sieve_type::supportSet> edges   = typename Mesh::sieve_type::supportSet();
1693:                 corners->clear();corners->insert(cornerA);corners->insert(cornerB);
1694:                 edges->insert(*newS->nJoin1(corners)->begin());
1695:                 corners->clear();corners->insert(cornerB);corners->insert(cornerC);
1696:                 edges->insert(*newS->nJoin1(corners)->begin());
1697:                 corners->clear();corners->insert(cornerC);corners->insert(cornerA);
1698:                 edges->insert(*newS->nJoin1(corners)->begin());
1699:                 ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(30);
1700:                 const typename Mesh::point_type          face       = *newS->nJoin1(edges)->begin();
1701:                 const int                                faceMarker = out.trifacemarkerlist[f];

1703:                 ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*newSieve, face, pV);
1704:                 const size_t                     n    = pV.getSize();
1705:                 const typename Mesh::point_type *cone = pV.getPoints();

1707:                 for(size_t c = 0; c < n; ++c) {
1708:                   refMesh->setValue(newMarkers, cone[c], faceMarker);
1709:                 }
1710:                 pV.clear();
1711:               }
1712:             }
1713:           }
1714:         }
1715: #if 0
1716:         if (refMesh->commSize() > 1) {
1717:           return ALE::Distribution<Mesh>::distributeMesh(refMesh);
1718:         }
1719: #endif
1720:         return refMesh;
1721:       };
1722:       static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false) {
1723:         throw ALE::Exception("Not yet implemented");
1724:       };
1725:       static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false) {
1726:         const int numCells         = mesh->heightStratum(0)->size();
1727:         double   *serialMaxVolumes = new double[numCells];

1729:         for(int c = 0; c < numCells; c++) {
1730:           serialMaxVolumes[c] = maxVolume;
1731:         }
1732:         const Obj<Mesh> refMesh = refineMeshV(mesh, serialMaxVolumes, interpolate);
1733:         delete [] serialMaxVolumes;
1734:         return refMesh;
1735:       };
1736:     };
1737:   };
1738: #endif
1739:   template<typename Mesh>
1740:   class Generator {
1741:   public:
1742:     static Obj<Mesh> generateMesh(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false) {
1743:       int dim = boundary->getDimension();

1745:       if (dim == 1) {
1746: #ifdef PETSC_HAVE_TRIANGLE
1747:         return ALE::Triangle::Generator<Mesh>::generateMesh(boundary, interpolate, constrained);
1748: #else
1749:         throw ALE::Exception("Mesh generation currently requires Triangle to be installed. Use --download-triangle during configure.");
1750: #endif
1751:       } else if (dim == 2) {
1752: #ifdef PETSC_HAVE_TETGEN
1753:         return ALE::TetGen::Generator<Mesh>::generateMesh(boundary, interpolate, constrained);
1754: #else
1755:         throw ALE::Exception("Mesh generation currently requires TetGen to be installed. Use --download-tetgen during configure.");
1756: #endif
1757:       }
1758:       return NULL;
1759:     };
1760:     static Obj<Mesh> generateMeshV(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false) {
1761:       int dim = boundary->getDimension();

1763:       if (dim == 1) {
1764: #ifdef PETSC_HAVE_TRIANGLE
1765:         return ALE::Triangle::Generator<Mesh>::generateMeshV(boundary, interpolate, constrained);
1766: #else
1767:         throw ALE::Exception("Mesh generation currently requires Triangle to be installed. Use --download-triangle during configure.");
1768: #endif
1769:       } else if (dim == 2) {
1770: #ifdef PETSC_HAVE_TETGEN
1771:         return ALE::TetGen::Generator<Mesh>::generateMeshV(boundary, interpolate, constrained);
1772: #else
1773:         throw ALE::Exception("Mesh generation currently requires TetGen to be installed. Use --download-tetgen during configure.");
1774: #endif
1775:       }
1776:       return NULL;
1777:     };
1778:     static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false) {
1779:       int dim = mesh->getDimension();

1781:       if (dim == 2) {
1782: #ifdef PETSC_HAVE_TRIANGLE
1783:         return ALE::Triangle::Refiner<Mesh>::refineMesh(mesh, maxVolumes, interpolate);
1784: #else
1785:         throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
1786: #endif
1787:       } else if (dim == 3) {
1788: #ifdef PETSC_HAVE_TETGEN
1789:         return ALE::TetGen::Refiner<Mesh>::refineMesh(mesh, maxVolumes, interpolate);
1790: #else
1791:         throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
1792: #endif
1793:       }
1794:       return NULL;
1795:     };
1796:     static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false, const bool forceSerial = false) {
1797:       int dim = mesh->getDimension();

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

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

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

1853:       if (dim == 2) {
1854: #ifdef PETSC_HAVE_TRIANGLE
1855:         return ALE::Triangle::Refiner<Mesh>::refineMeshLocal(mesh, maxVolume, interpolate);
1856: #else
1857:         throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
1858: #endif
1859:       } else if (dim == 3) {
1860: #ifdef PETSC_HAVE_TETGEN
1861:         return ALE::TetGen::Refiner<Mesh>::refineMesh(mesh, maxVolume, interpolate);
1862: #else
1863:         throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
1864: #endif
1865:       }
1866:       return NULL;
1867:     };
1868:   };
1869: }

1871: #endif