Actual source code: Generator.hh
petsc-3.4.5 2014-06-29
1: #ifndef included_ALE_Generator_hh
2: #define included_ALE_Generator_hh
4: #ifndef included_ALE_Distribution_hh
5: #include <sieve/Distribution.hh>
6: #endif
8: #ifdef PETSC_HAVE_TRIANGLE
9: #include <triangle.h>
10: #endif
11: #ifdef PETSC_HAVE_TETGEN
12: #include <tetgen.h>
13: #endif
15: namespace ALE {
16: #ifdef PETSC_HAVE_TRIANGLE
17: namespace Triangle {
18: template<typename Mesh>
19: class Generator {
20: class SegmentVisitor {
21: protected:
22: const int dim;
23: int *segmentlist;
24: typename Mesh::numbering_type& vNumbering;
25: int idx, v;
26: public:
27: SegmentVisitor(const int dim, int segmentlist[], typename Mesh::numbering_type& vNumbering) : dim(dim), segmentlist(segmentlist), vNumbering(vNumbering), idx(0), v(0) {};
28: ~SegmentVisitor() {};
29: public:
30: template<typename Point>
31: void visitPoint(const Point& point) {
32: this->segmentlist[this->idx*dim + (this->v++)] = this->vNumbering.getIndex(point);
33: }
34: template<typename Arrow>
35: void visitArrow(const Arrow& arrow) {}
36: void setIndex(const int idx) {this->idx = idx; this->v = 0;};
37: };
38: public:
39: static void initInput(struct triangulateio *inputCtx) {
40: inputCtx->numberofpoints = 0;
41: inputCtx->numberofpointattributes = 0;
42: inputCtx->pointlist = NULL;
43: inputCtx->pointattributelist = NULL;
44: inputCtx->pointmarkerlist = NULL;
45: inputCtx->numberofsegments = 0;
46: inputCtx->segmentlist = NULL;
47: inputCtx->segmentmarkerlist = NULL;
48: inputCtx->numberoftriangleattributes = 0;
49: inputCtx->trianglelist = NULL;
50: inputCtx->numberofholes = 0;
51: inputCtx->holelist = NULL;
52: inputCtx->numberofregions = 0;
53: inputCtx->regionlist = NULL;
54: };
55: static void initOutput(struct triangulateio *outputCtx) {
56: outputCtx->numberofpoints = 0;
57: outputCtx->pointlist = NULL;
58: outputCtx->pointattributelist = NULL;
59: outputCtx->pointmarkerlist = NULL;
60: outputCtx->numberoftriangles = 0;
61: outputCtx->trianglelist = NULL;
62: outputCtx->triangleattributelist = NULL;
63: outputCtx->neighborlist = NULL;
64: outputCtx->segmentlist = NULL;
65: outputCtx->segmentmarkerlist = NULL;
66: outputCtx->edgelist = NULL;
67: outputCtx->edgemarkerlist = NULL;
68: };
69: static void finiOutput(struct triangulateio *outputCtx) {
70: free(outputCtx->pointmarkerlist);
71: free(outputCtx->edgelist);
72: free(outputCtx->edgemarkerlist);
73: free(outputCtx->trianglelist);
74: free(outputCtx->neighborlist);
75: };
78: static Obj<Mesh> generateMesh(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false) {
79: int dim = 2;
80: Obj<Mesh> mesh = new Mesh(boundary->comm(), dim, boundary->debug());
81: const Obj<typename Mesh::sieve_type>& sieve = boundary->getSieve();
82: const bool createConvexHull = false;
83: struct triangulateio in;
84: struct triangulateio out;
87: initInput(&in);
88: initOutput(&out);
89: const Obj<typename Mesh::label_sequence>& vertices = boundary->depthStratum(0);
90: const Obj<typename Mesh::label_type>& markers = boundary->getLabel("marker");
91: const Obj<typename Mesh::real_section_type>& coordinates = boundary->getRealSection("coordinates");
92: const Obj<typename Mesh::numbering_type>& vNumbering = boundary->getFactory()->getLocalNumbering(boundary, 0);
94: in.numberofpoints = vertices->size();
95: if (in.numberofpoints > 0) {
96: const typename Mesh::label_sequence::iterator vEnd = vertices->end();
98: PetscMalloc(in.numberofpoints * dim * sizeof(double), &in.pointlist);CHKERRXX(ierr);
99: PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);CHKERRXX(ierr);
100: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
101: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
102: const int idx = vNumbering->getIndex(*v_iter);
104: for(int d = 0; d < dim; d++) {
105: in.pointlist[idx*dim + d] = array[d];
106: }
107: in.pointmarkerlist[idx] = boundary->getValue(markers, *v_iter);
108: }
109: }
110: const Obj<typename Mesh::label_sequence>& edges = boundary->depthStratum(1);
111: const Obj<typename Mesh::numbering_type>& eNumbering = boundary->getFactory()->getLocalNumbering(boundary, 1);
113: in.numberofsegments = edges->size();
114: if (in.numberofsegments > 0) {
115: const typename Mesh::label_sequence::iterator eEnd = edges->end();
117: PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);CHKERRXX(ierr);
118: PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);CHKERRXX(ierr);
119: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != eEnd; ++e_iter) {
120: const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = sieve->cone(*e_iter);
121: const typename Mesh::sieve_type::traits::coneSequence::iterator cEnd = cone->end();
122: const int idx = eNumbering->getIndex(*e_iter);
123: int v = 0;
125: for(typename Mesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
126: in.segmentlist[idx*dim + (v++)] = vNumbering->getIndex(*c_iter);
127: }
128: in.segmentmarkerlist[idx] = boundary->getValue(markers, *e_iter);
129: }
130: }
131: const typename Mesh::holes_type& holes = boundary->getHoles();
133: in.numberofholes = holes.size();
134: if (in.numberofholes > 0) {
135: PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);CHKERRXX(ierr);
136: for(int h = 0; h < in.numberofholes; ++h) {
137: for(int d = 0; d < dim; ++d) {
138: in.holelist[h*dim+d] = holes[h][d];
139: }
140: }
141: }
142: if (mesh->commRank() == 0) {
143: std::string args("pqezQ");
145: if (createConvexHull) {
146: args += "c";
147: }
148: if (constrained) {
149: args = "zepDQ";
150: }
151: triangulate((char *) args.c_str(), &in, &out, NULL);
152: }
154: if (in.pointlist) {PetscFree(in.pointlist);CHKERRXX(ierr);}
155: if (in.pointmarkerlist) {PetscFree(in.pointmarkerlist);CHKERRXX(ierr);}
156: if (in.segmentlist) {PetscFree(in.segmentlist);CHKERRXX(ierr);}
157: if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);CHKERRXX(ierr);}
158: if (in.holelist) {PetscFree(in.holelist);CHKERRXX(ierr);}
159: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
160: int numCorners = 3;
161: int numCells = out.numberoftriangles;
162: int *cells = out.trianglelist;
163: int numVertices = out.numberofpoints;
164: double *coords = out.pointlist;
166: ALE::SieveBuilder<Mesh>::buildTopology(newSieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, mesh->getArrowSection("orientation"));
167: mesh->setSieve(newSieve);
168: mesh->stratify();
169: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, dim, coords);
170: const Obj<typename Mesh::label_type>& newMarkers = mesh->createLabel("marker");
172: if (mesh->commRank() == 0) {
173: for(int v = 0; v < out.numberofpoints; v++) {
174: if (out.pointmarkerlist[v]) {
175: mesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
176: }
177: }
178: if (interpolate) {
179: for(int e = 0; e < out.numberofedges; e++) {
180: if (out.edgemarkerlist[e]) {
181: const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
182: const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
183: const Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(vertexA, vertexB, 1);
185: mesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
186: }
187: }
188: }
189: }
190: mesh->copyHoles(boundary);
191: finiOutput(&out);
192: return mesh;
193: };
196: static Obj<Mesh> generateMeshV(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false, const bool renumber = false) {
197: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
198: typedef typename Mesh::real_section_type::value_type real;
199: int dim = 2;
200: const Obj<Mesh> mesh = new Mesh(boundary->comm(), dim, boundary->debug());
201: const Obj<typename Mesh::sieve_type>& sieve = boundary->getSieve();
202: const bool createConvexHull = false;
203: struct triangulateio in;
204: struct triangulateio out;
205: PetscErrorCode ierr;
207: initInput(&in);
208: initOutput(&out);
209: const Obj<typename Mesh::label_sequence>& vertices = boundary->depthStratum(0);
210: const Obj<typename Mesh::label_type>& markers = boundary->getLabel("marker");
211: const Obj<typename Mesh::real_section_type>& coordinates = boundary->getRealSection("coordinates");
212: const Obj<typename Mesh::numbering_type>& vNumbering = boundary->getFactory()->getLocalNumbering(boundary, 0);
214: in.numberofpoints = vertices->size();
215: if (in.numberofpoints > 0) {
216: const typename Mesh::label_sequence::iterator vEnd = vertices->end();
218: PetscMalloc(in.numberofpoints * dim * sizeof(double), &in.pointlist);CHKERRXX(ierr);
219: PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);CHKERRXX(ierr);
220: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
221: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
222: const int idx = vNumbering->getIndex(*v_iter);
224: for(int d = 0; d < dim; ++d) {
225: in.pointlist[idx*dim + d] = array[d];
226: }
227: in.pointmarkerlist[idx] = boundary->getValue(markers, *v_iter);
228: }
229: }
230: const Obj<typename Mesh::label_sequence>& edges = boundary->depthStratum(1);
231: const Obj<typename Mesh::numbering_type>& eNumbering = boundary->getFactory()->getLocalNumbering(boundary, 1);
233: in.numberofsegments = edges->size();
234: if (in.numberofsegments > 0) {
235: const typename Mesh::label_sequence::iterator eEnd = edges->end();
237: PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);CHKERRXX(ierr);
238: PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);CHKERRXX(ierr);
239: SegmentVisitor sV(dim, in.segmentlist, *vNumbering);
240: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != eEnd; ++e_iter) {
241: const int idx = eNumbering->getIndex(*e_iter);
243: sV.setIndex(idx);
244: sieve->cone(*e_iter, sV);
245: in.segmentmarkerlist[idx] = boundary->getValue(markers, *e_iter);
246: }
247: }
248: const typename Mesh::holes_type& holes = boundary->getHoles();
250: in.numberofholes = holes.size();
251: if (in.numberofholes > 0) {
252: PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);CHKERRXX(ierr);
253: for(int h = 0; h < in.numberofholes; ++h) {
254: for(int d = 0; d < dim; ++d) {
255: in.holelist[h*dim+d] = holes[h][d];
256: }
257: }
258: }
259: if (mesh->commRank() == 0) {
260: std::string args("pqezQ");
262: if (createConvexHull) {
263: args += "c";
264: }
265: if (constrained) {
266: args = "zepDQ";
267: }
268: triangulate((char *) args.c_str(), &in, &out, NULL);
269: }
270: if (in.pointlist) {PetscFree(in.pointlist);CHKERRXX(ierr);}
271: if (in.pointmarkerlist) {PetscFree(in.pointmarkerlist);CHKERRXX(ierr);}
272: if (in.segmentlist) {PetscFree(in.segmentlist);CHKERRXX(ierr);}
273: if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);CHKERRXX(ierr);}
274: if (in.holelist) {PetscFree(in.holelist);CHKERRXX(ierr);}
275: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
276: const Obj<FlexMesh> m = new FlexMesh(boundary->comm(), dim, boundary->debug());
277: const Obj<FlexMesh::sieve_type> newS = new FlexMesh::sieve_type(m->comm(), m->debug());
278: int numCorners = 3;
279: int numCells = out.numberoftriangles;
280: int *cells = out.trianglelist;
281: int numVertices = out.numberofpoints;
282: double *coords = out.pointlist;
283: real *coordsR;
285: ALE::SieveBuilder<FlexMesh>::buildTopology(newS, dim, numCells, cells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
286: m->setSieve(newS);
287: m->stratify();
288: mesh->setSieve(newSieve);
289: std::map<typename Mesh::point_type,typename Mesh::point_type> renumbering;
290: ALE::ISieveConverter::convertSieve(*newS, *newSieve, renumbering, renumber);
291: mesh->stratify();
292: ALE::ISieveConverter::convertOrientation(*newS, *newSieve, renumbering,
293: m->getArrowSection("orientation").ptr());
294: {
295: if (sizeof(double) == sizeof(real)) {
296: coordsR = (real *) coords;
297: } else {
298: coordsR = new real[numVertices*dim];
299: for(int i = 0; i < numVertices*dim; ++i) coordsR[i] = coords[i];
300: }
301: }
302: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, dim, coordsR);
303: {
304: if (sizeof(double) != sizeof(real)) {
305: delete [] coordsR;
306: }
307: }
308: const Obj<typename Mesh::label_type>& newMarkers = mesh->createLabel("marker");
310: if (mesh->commRank() == 0) {
311: #ifdef IMESH_NEW_LABELS
312: int size = 0;
314: newMarkers->setChart(mesh->getSieve()->getChart());
315: for(int v = 0; v < out.numberofpoints; v++) {
316: if (out.pointmarkerlist[v]) {
317: newMarkers->setConeSize(v+out.numberoftriangles, 1);
318: size++;
319: }
320: }
321: if (interpolate) {
322: for(int e = 0; e < out.numberofedges; e++) {
323: if (out.edgemarkerlist[e]) {
324: const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
325: const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
326: const Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(vertexA, vertexB, 1);
328: newMarkers->setConeSize(*edge->begin(), 1);
329: size++;
330: }
331: }
332: }
333: newMarkers->setSupportSize(0, size);
334: newMarkers->allocate();
335: #endif
336: for(int v = 0; v < out.numberofpoints; v++) {
337: if (out.pointmarkerlist[v]) {
338: if (renumber) {
339: mesh->setValue(newMarkers, renumbering[v+out.numberoftriangles], out.pointmarkerlist[v]);
340: } else {
341: mesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
342: }
343: }
344: }
345: if (interpolate) {
346: for(int e = 0; e < out.numberofedges; e++) {
347: if (out.edgemarkerlist[e]) {
348: const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
349: const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
350: const Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(vertexA, vertexB, 1);
352: mesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
353: }
354: }
355: }
356: #ifdef IMESH_NEW_LABELS
357: newMarkers->recalculateLabel();
358: #endif
359: }
360: mesh->copyHoles(boundary);
361: finiOutput(&out);
362: return mesh;
363: };
364: };
365: template<typename Mesh>
366: class Refiner {
367: public:
368: static Obj<Mesh> refineMesh(const Obj<Mesh>& serialMesh, const double maxVolumes[], const bool interpolate = false, const bool forceSerial = false) {
369: const int dim = serialMesh->getDimension();
370: const Obj<Mesh> refMesh = new Mesh(serialMesh->comm(), dim, serialMesh->debug());
371: const Obj<typename Mesh::sieve_type>& serialSieve = serialMesh->getSieve();
372: struct triangulateio in;
373: struct triangulateio out;
374: PetscErrorCode ierr;
376: Generator<Mesh>::initInput(&in);
377: Generator<Mesh>::initOutput(&out);
378: const Obj<typename Mesh::label_sequence>& vertices = serialMesh->depthStratum(0);
379: const Obj<typename Mesh::label_type>& markers = serialMesh->getLabel("marker");
380: const Obj<typename Mesh::real_section_type>& coordinates = serialMesh->getRealSection("coordinates");
381: const Obj<typename Mesh::numbering_type>& vNumbering = serialMesh->getFactory()->getLocalNumbering(serialMesh, 0);
383: in.numberofpoints = vertices->size();
384: if (in.numberofpoints > 0) {
385: const typename Mesh::label_sequence::iterator vEnd = vertices->end();
387: PetscMalloc(in.numberofpoints * dim * sizeof(double), &in.pointlist);CHKERRXX(ierr);
388: PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);CHKERRXX(ierr);
389: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
390: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
391: const int idx = vNumbering->getIndex(*v_iter);
393: for(int d = 0; d < dim; d++) {
394: in.pointlist[idx*dim + d] = array[d];
395: }
396: in.pointmarkerlist[idx] = serialMesh->getValue(markers, *v_iter);
397: }
398: }
399: const Obj<typename Mesh::label_sequence>& faces = serialMesh->heightStratum(0);
400: const Obj<typename Mesh::numbering_type>& fNumbering = serialMesh->getFactory()->getLocalNumbering(serialMesh, serialMesh->depth());
402: in.numberofcorners = 3;
403: in.numberoftriangles = faces->size();
404: in.trianglearealist = (double *) maxVolumes;
405: if (in.numberoftriangles > 0) {
406: const typename Mesh::label_sequence::iterator fEnd = faces->end();
408: PetscMalloc(in.numberoftriangles*in.numberofcorners * sizeof(int), &in.trianglelist);CHKERRXX(ierr);
409: if (serialMesh->depth() == 1) {
410: for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
411: const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = serialSieve->cone(*f_iter);
412: const typename Mesh::sieve_type::traits::coneSequence::iterator cEnd = cone->end();
413: const int idx = fNumbering->getIndex(*f_iter);
414: int v = 0;
416: for(typename Mesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
417: in.trianglelist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*c_iter);
418: }
419: }
420: } else if (serialMesh->depth() == 2) {
421: for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
422: typedef ALE::SieveAlg<Mesh> sieve_alg_type;
423: const Obj<typename sieve_alg_type::coneArray>& cone = sieve_alg_type::nCone(serialMesh, *f_iter, 2);
424: const typename Mesh::sieve_type::coneArray::iterator cEnd = cone->end();
425: const int idx = fNumbering->getIndex(*f_iter);
426: int v = 0;
428: for(typename Mesh::sieve_type::coneArray::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
429: in.trianglelist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*c_iter);
430: }
431: }
432: } else {
433: throw ALE::Exception("Invalid sieve: Cannot gives sieves of arbitrary depth to Triangle");
434: }
435: }
436: if (serialMesh->depth() == 2) {
437: const Obj<typename Mesh::label_sequence>& edges = serialMesh->depthStratum(1);
438: #define NEW_LABEL
439: #ifdef NEW_LABEL
440: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
441: if (serialMesh->getValue(markers, *e_iter)) {
442: in.numberofsegments++;
443: }
444: }
445: //std::cout << "Number of segments: " << in.numberofsegments << std::endl;
446: if (in.numberofsegments > 0) {
447: int s = 0;
449: PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);CHKERRXX(ierr);
450: PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);CHKERRXX(ierr);
451: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
452: const int edgeMarker = serialMesh->getValue(markers, *e_iter);
454: if (edgeMarker) {
455: const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = serialSieve->cone(*e_iter);
456: int p = 0;
458: for(typename Mesh::sieve_type::traits::coneSequence::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
459: in.segmentlist[s*2 + (p++)] = vNumbering->getIndex(*v_iter);
460: }
461: in.segmentmarkerlist[s++] = edgeMarker;
462: }
463: }
464: }
465: #else
466: const Obj<typename Mesh::label_type::baseSequence>& boundary = markers->base();
468: in.numberofsegments = 0;
469: for(typename Mesh::label_type::baseSequence::iterator b_iter = boundary->begin(); b_iter != boundary->end(); ++b_iter) {
470: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
471: if (*b_iter == *e_iter) {
472: in.numberofsegments++;
473: }
474: }
475: }
476: if (in.numberofsegments > 0) {
477: int s = 0;
479: PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);CHKERRXX(ierr);
480: PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);CHKERRXX(ierr);
481: for(typename Mesh::label_type::baseSequence::iterator b_iter = boundary->begin(); b_iter != boundary->end(); ++b_iter) {
482: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
483: if (*b_iter == *e_iter) {
484: const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = serialSieve->cone(*e_iter);
485: int p = 0;
487: for(typename Mesh::sieve_type::traits::coneSequence::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
488: in.segmentlist[s*2 + (p++)] = vNumbering->getIndex(*v_iter);
489: }
490: in.segmentmarkerlist[s++] = serialMesh->getValue(markers, *e_iter);
491: }
492: }
493: }
494: }
495: #endif
496: }
498: in.numberofholes = 0;
499: if (in.numberofholes > 0) {
500: PetscMalloc(in.numberofholes * dim * sizeof(int), &in.holelist);CHKERRXX(ierr);
501: }
502: if (serialMesh->commRank() == 0) {
503: std::string args("pqezQra");
505: triangulate((char *) args.c_str(), &in, &out, NULL);
506: }
507: if (in.pointlist) {PetscFree(in.pointlist);CHKERRXX(ierr);}
508: if (in.pointmarkerlist) {PetscFree(in.pointmarkerlist);CHKERRXX(ierr);}
509: if (in.segmentlist) {PetscFree(in.segmentlist);CHKERRXX(ierr);}
510: if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);CHKERRXX(ierr);}
511: if (in.trianglelist) {PetscFree(in.trianglelist);CHKERRXX(ierr);}
512: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(serialMesh->comm(), serialMesh->debug());
513: int numCorners = 3;
514: int numCells = out.numberoftriangles;
515: int *cells = out.trianglelist;
516: int numVertices = out.numberofpoints;
517: double *coords = out.pointlist;
519: ALE::SieveBuilder<Mesh>::buildTopology(newSieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, refMesh->getArrowSection("orientation"));
520: refMesh->setSieve(newSieve);
521: refMesh->stratify();
522: ALE::SieveBuilder<Mesh>::buildCoordinates(refMesh, dim, coords);
523: const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");
525: if (refMesh->commRank() == 0) {
526: for(int v = 0; v < out.numberofpoints; v++) {
527: if (out.pointmarkerlist[v]) {
528: refMesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
529: }
530: }
531: if (interpolate) {
532: for(int e = 0; e < out.numberofedges; e++) {
533: if (out.edgemarkerlist[e]) {
534: const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
535: const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
536: const Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(vertexA, vertexB, 1);
538: refMesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
539: }
540: }
541: }
542: }
544: Generator<Mesh>::finiOutput(&out);
545: if ((refMesh->commSize() > 1) && (!forceSerial)) {
546: return ALE::Distribution<Mesh>::distributeMesh(refMesh);
547: }
548: return refMesh;
549: };
550: static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false) {
551: Obj<Mesh> serialMesh = ALE::Distribution<Mesh>::unifyMesh(mesh);
552: const Obj<typename Mesh::real_section_type> serialMaxVolumes = ALE::Distribution<Mesh>::distributeSection(maxVolumes, serialMesh, serialMesh->getDistSendOverlap(), serialMesh->getDistRecvOverlap());
554: return refineMesh(serialMesh, serialMaxVolumes->restrictSpace(), interpolate);
555: };
556: static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false, const bool forceSerial = false) {
557: Obj<Mesh> serialMesh;
558: if ((mesh->commSize() > 1) && (!forceSerial)) {
559: serialMesh = ALE::Distribution<Mesh>::unifyMesh(mesh);
560: } else {
561: serialMesh = mesh;
562: }
563: const int numFaces = serialMesh->heightStratum(0)->size();
564: double *serialMaxVolumes = new double[numFaces];
566: for(int f = 0; f < numFaces; f++) {
567: serialMaxVolumes[f] = maxVolume;
568: }
569: const Obj<Mesh> refMesh = refineMesh(serialMesh, serialMaxVolumes, interpolate, forceSerial);
570: delete [] serialMaxVolumes;
571: return refMesh;
572: };
573: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolumes[], const bool interpolate = false, const bool forceSerial = false, const bool renumber = false) {
574: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
575: typedef typename Mesh::real_section_type::value_type real;
576: typedef typename Mesh::point_type point_type;
577: const int dim = mesh->getDimension();
578: const Obj<Mesh> refMesh = new Mesh(mesh->comm(), dim, mesh->debug());
579: const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
580: struct triangulateio in;
581: struct triangulateio out;
582: PetscErrorCode ierr;
584: Generator<Mesh>::initInput(&in);
585: Generator<Mesh>::initOutput(&out);
586: const Obj<typename Mesh::label_sequence>& vertices = mesh->depthStratum(0);
587: const Obj<typename Mesh::label_type>& markers = mesh->getLabel("marker");
588: const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
589: const Obj<typename Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
591: in.numberofpoints = vertices->size();
592: if (in.numberofpoints > 0) {
593: const typename Mesh::label_sequence::iterator vEnd = vertices->end();
595: PetscMalloc(in.numberofpoints * dim * sizeof(double), &in.pointlist);CHKERRXX(ierr);
596: PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);CHKERRXX(ierr);
597: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
598: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
599: const int idx = vNumbering->getIndex(*v_iter);
601: for(int d = 0; d < dim; d++) {
602: in.pointlist[idx*dim + d] = array[d];
603: }
604: in.pointmarkerlist[idx] = mesh->getValue(markers, *v_iter);
605: }
606: }
607: const Obj<typename Mesh::label_sequence>& faces = mesh->heightStratum(0);
608: const Obj<typename Mesh::numbering_type>& fNumbering = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());
610: in.numberofcorners = 3;
611: in.numberoftriangles = faces->size();
612: in.trianglearealist = (double *) maxVolumes;
613: if (in.numberoftriangles > 0) {
614: PetscMalloc(in.numberoftriangles*in.numberofcorners * sizeof(int), &in.trianglelist);CHKERRXX(ierr);
615: if (mesh->depth() == 1) {
616: ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(3);
617: const typename Mesh::label_sequence::iterator fEnd = faces->end();
619: for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
620: sieve->cone(*f_iter, pV);
621: const int idx = fNumbering->getIndex(*f_iter);
622: const size_t n = pV.getSize();
623: const point_type *cone = pV.getPoints();
625: assert(n == 3);
626: for(int v = 0; v < 3; ++v) {
627: in.trianglelist[idx*in.numberofcorners + v] = vNumbering->getIndex(cone[v]);
628: }
629: pV.clear();
630: }
631: } else if (mesh->depth() == 2) {
632: // Need extra space due to early error checking
633: ALE::ISieveVisitor::NConeRetriever<typename Mesh::sieve_type> ncV(*sieve, 4);
634: const typename Mesh::label_sequence::iterator fEnd = faces->end();
636: for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
637: ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*sieve, *f_iter, ncV);
638: const int idx = fNumbering->getIndex(*f_iter);
639: const size_t n = ncV.getSize();
640: const point_type *cone = ncV.getPoints();
642: assert(n == 3);
643: for(int v = 0; v < 3; ++v) {
644: in.trianglelist[idx*in.numberofcorners + v] = vNumbering->getIndex(cone[v]);
645: }
646: ncV.clear();
647: }
648: } else {
649: throw ALE::Exception("Invalid sieve: Cannot gives sieves of arbitrary depth to Triangle");
650: }
651: }
652: if (mesh->depth() == 2) {
653: const Obj<typename Mesh::label_sequence>& edges = mesh->depthStratum(1);
654: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
655: if (mesh->getValue(markers, *e_iter)) {
656: in.numberofsegments++;
657: }
658: }
659: //std::cout << "Number of segments: " << in.numberofsegments << std::endl;
660: if (in.numberofsegments > 0) {
661: ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(2);
662: int s = 0;
664: PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);CHKERRXX(ierr);
665: PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);CHKERRXX(ierr);
666: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
667: const int edgeMarker = mesh->getValue(markers, *e_iter);
669: if (edgeMarker) {
670: sieve->cone(*e_iter, pV);
671: const size_t n = pV.getSize();
672: const point_type *cone = pV.getPoints();
674: assert(n == 2);
675: for(int v = 0; v < 2; ++v) {
676: in.segmentlist[s*2 + v] = vNumbering->getIndex(cone[v]);
677: }
678: in.segmentmarkerlist[s++] = edgeMarker;
679: pV.clear();
680: }
681: }
682: }
683: }
684: const typename Mesh::holes_type& holes = mesh->getHoles();
686: in.numberofholes = holes.size();
687: if (in.numberofholes > 0) {
688: PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);CHKERRXX(ierr);
689: for(int h = 0; h < in.numberofholes; ++h) {
690: for(int d = 0; d < dim; ++d) {
691: in.holelist[h*dim+d] = holes[h][d];
692: }
693: }
694: }
695: if (mesh->commRank() == 0) {
696: std::string args("pqezQra");
698: triangulate((char *) args.c_str(), &in, &out, NULL);
699: }
700: if (in.pointlist) {PetscFree(in.pointlist);CHKERRXX(ierr);}
701: if (in.pointmarkerlist) {PetscFree(in.pointmarkerlist);CHKERRXX(ierr);}
702: if (in.segmentlist) {PetscFree(in.segmentlist);CHKERRXX(ierr);}
703: if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);CHKERRXX(ierr);}
704: if (in.trianglelist) {PetscFree(in.trianglelist);CHKERRXX(ierr);}
705: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
706: const Obj<FlexMesh> m = new FlexMesh(mesh->comm(), dim, mesh->debug());
707: const Obj<FlexMesh::sieve_type> newS = new FlexMesh::sieve_type(m->comm(), m->debug());
708: int numCorners = 3;
709: int numCells = out.numberoftriangles;
710: int *cells = out.trianglelist;
711: int numVertices = out.numberofpoints;
712: double *coords = out.pointlist;
713: real *coordsR;
715: ALE::SieveBuilder<FlexMesh>::buildTopology(newS, dim, numCells, cells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
716: m->setSieve(newS);
717: m->stratify();
718: refMesh->setSieve(newSieve);
719: std::map<typename Mesh::point_type,typename Mesh::point_type> renumbering;
720: ALE::ISieveConverter::convertSieve(*newS, *newSieve, renumbering, renumber);
721: refMesh->stratify();
722: ALE::ISieveConverter::convertOrientation(*newS, *newSieve, renumbering, m->getArrowSection("orientation").ptr());
723: {
724: if (sizeof(double) == sizeof(real)) {
725: coordsR = (real *) coords;
726: } else {
727: coordsR = new real[numVertices*dim];
728: for(int i = 0; i < numVertices*dim; ++i) coordsR[i] = coords[i];
729: }
730: }
731: ALE::SieveBuilder<Mesh>::buildCoordinates(refMesh, dim, coordsR);
732: {
733: if (sizeof(double) != sizeof(real)) {
734: delete [] coordsR;
735: }
736: }
737: const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");
739: if (refMesh->commRank() == 0) {
740: for(int v = 0; v < out.numberofpoints; v++) {
741: if (out.pointmarkerlist[v]) {
742: if (renumber) {
743: refMesh->setValue(newMarkers, renumbering[v+out.numberoftriangles], out.pointmarkerlist[v]);
744: } else {
745: refMesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
746: }
747: }
748: }
749: if (interpolate) {
750: for(int e = 0; e < out.numberofedges; e++) {
751: if (out.edgemarkerlist[e]) {
752: const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
753: const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
754: const Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(vertexA, vertexB, 1);
756: refMesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
757: }
758: }
759: }
760: }
761: refMesh->copyHoles(mesh);
762: Generator<Mesh>::finiOutput(&out);
763: #if 0
764: if ((refMesh->commSize() > 1) && (!forceSerial)) {
765: return ALE::Distribution<Mesh>::distributeMesh(refMesh);
766: }
767: #endif
768: return refMesh;
769: };
770: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false, const bool renumber = false) {
771: throw ALE::Exception("Not yet implemented");
772: };
773: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false, const bool forceSerial = false, const bool renumber = false) {
774: #if 0
775: Obj<Mesh> serialMesh;
776: if (mesh->commSize() > 1) {
777: serialMesh = ALE::Distribution<Mesh>::unifyMesh(mesh);
778: } else {
779: serialMesh = mesh;
780: }
781: #endif
782: const int numCells = mesh->heightStratum(0)->size();
783: double *serialMaxVolumes = new double[numCells];
785: for(int c = 0; c < numCells; c++) {
786: serialMaxVolumes[c] = maxVolume;
787: }
788: const Obj<Mesh> refMesh = refineMeshV(mesh, serialMaxVolumes, interpolate, forceSerial, renumber);
789: delete [] serialMaxVolumes;
790: return refMesh;
791: };
792: static Obj<Mesh> refineMeshLocal(const Obj<Mesh>& mesh, const double maxVolumes[], const bool interpolate = false) {
793: const int dim = mesh->getDimension();
794: const Obj<Mesh> refMesh = new Mesh(mesh->comm(), dim, mesh->debug());
795: const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
796: struct triangulateio in;
797: struct triangulateio out;
798: PetscErrorCode ierr;
800: Generator<Mesh>::initInput(&in);
801: Generator<Mesh>::initOutput(&out);
802: const Obj<typename Mesh::label_sequence>& vertices = mesh->depthStratum(0);
803: const Obj<typename Mesh::label_type>& markers = mesh->getLabel("marker");
804: const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
805: const Obj<typename Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
807: in.numberofpoints = vertices->size();
808: if (in.numberofpoints > 0) {
809: PetscMalloc(in.numberofpoints * dim * sizeof(double), &in.pointlist);CHKERRXX(ierr);
810: PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);CHKERRXX(ierr);
811: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
812: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
813: const int idx = vNumbering->getIndex(*v_iter);
815: for(int d = 0; d < dim; d++) {
816: in.pointlist[idx*dim + d] = array[d];
817: }
818: in.pointmarkerlist[idx] = mesh->getValue(markers, *v_iter);
819: }
820: }
821: const Obj<typename Mesh::label_sequence>& faces = mesh->heightStratum(0);
822: const Obj<typename Mesh::numbering_type>& fNumbering = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());
824: in.numberofcorners = 3;
825: in.numberoftriangles = faces->size();
826: in.trianglearealist = (double *) maxVolumes;
827: if (in.numberoftriangles > 0) {
828: PetscMalloc(in.numberoftriangles*in.numberofcorners * sizeof(int), &in.trianglelist);CHKERRXX(ierr);
829: if (mesh->depth() == 1) {
830: for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != faces->end(); ++f_iter) {
831: const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = sieve->cone(*f_iter);
832: const int idx = fNumbering->getIndex(*f_iter);
833: int v = 0;
835: for(typename Mesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
836: in.trianglelist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*c_iter);
837: }
838: }
839: } else if (mesh->depth() == 2) {
840: for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != faces->end(); ++f_iter) {
841: typedef ALE::SieveAlg<Mesh> sieve_alg_type;
842: const Obj<typename sieve_alg_type::coneArray>& cone = sieve_alg_type::nCone(mesh, *f_iter, 2);
843: const int idx = fNumbering->getIndex(*f_iter);
844: int v = 0;
846: for(typename Mesh::sieve_type::coneArray::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
847: in.trianglelist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*c_iter);
848: }
849: }
850: } else {
851: throw ALE::Exception("Invalid sieve: Cannot gives sieves of arbitrary depth to Triangle");
852: }
853: }
854: if (mesh->depth() == 2) {
855: const Obj<typename Mesh::label_sequence>& edges = mesh->depthStratum(1);
856: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
857: if (mesh->getValue(markers, *e_iter)) {
858: in.numberofsegments++;
859: }
860: }
861: //std::cout << "Number of segments: " << in.numberofsegments << std::endl;
862: if (in.numberofsegments > 0) {
863: int s = 0;
865: PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);CHKERRXX(ierr);
866: PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);CHKERRXX(ierr);
867: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
868: const int edgeMarker = mesh->getValue(markers, *e_iter);
870: if (edgeMarker) {
871: const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = sieve->cone(*e_iter);
872: int p = 0;
874: for(typename Mesh::sieve_type::traits::coneSequence::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
875: in.segmentlist[s*2 + (p++)] = vNumbering->getIndex(*v_iter);
876: }
877: in.segmentmarkerlist[s++] = edgeMarker;
878: }
879: }
880: }
881: }
883: in.numberofholes = 0;
884: if (in.numberofholes > 0) {
885: PetscMalloc(in.numberofholes * dim * sizeof(int), &in.holelist);CHKERRXX(ierr);
886: }
887: std::string args("pqezQra");
889: triangulate((char *) args.c_str(), &in, &out, NULL);
890: if (in.pointlist) {PetscFree(in.pointlist);CHKERRXX(ierr);}
891: if (in.pointmarkerlist) {PetscFree(in.pointmarkerlist);CHKERRXX(ierr);}
892: if (in.segmentlist) {PetscFree(in.segmentlist);CHKERRXX(ierr);}
893: if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);CHKERRXX(ierr);}
894: if (in.trianglelist) {PetscFree(in.trianglelist);CHKERRXX(ierr);}
895: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
896: int numCorners = 3;
897: int numCells = out.numberoftriangles;
898: int *cells = out.trianglelist;
899: int numVertices = out.numberofpoints;
900: double *coords = out.pointlist;
902: ALE::SieveBuilder<Mesh>::buildTopologyMultiple(newSieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, refMesh->getArrowSection("orientation"));
903: refMesh->setSieve(newSieve);
904: refMesh->stratify();
905: ALE::SieveBuilder<Mesh>::buildCoordinatesMultiple(refMesh, dim, coords);
906: const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");
908: for(int v = 0; v < out.numberofpoints; v++) {
909: if (out.pointmarkerlist[v]) {
910: refMesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
911: }
912: }
913: if (interpolate) {
914: for(int e = 0; e < out.numberofedges; e++) {
915: if (out.edgemarkerlist[e]) {
916: const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
917: const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
918: const Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(vertexA, vertexB, 1);
920: refMesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
921: }
922: }
923: }
924: Generator<Mesh>::finiOutput(&out);
925: return refMesh;
926: };
927: static Obj<Mesh> refineMeshLocal(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false) {
928: const int numLocalFaces = mesh->heightStratum(0)->size();
929: double *localMaxVolumes = new double[numLocalFaces];
931: for(int f = 0; f < numLocalFaces; f++) {
932: localMaxVolumes[f] = maxVolume;
933: }
934: const Obj<Mesh> refMesh = refineMeshLocal(mesh, localMaxVolumes, interpolate);
935: const Obj<typename Mesh::sieve_type> refSieve = refMesh->getSieve();
936: delete [] localMaxVolumes;
937: #if 0
938: typedef typename ALE::New::Completion<Mesh, typename Mesh::sieve_type::point_type> sieveCompletion;
939: // This is where we enforce consistency over the overlap
940: // We need somehow to update the overlap to account for the new stuff
941: //
942: // 1) Since we are refining only, the vertices are invariant
943: // 2) We need to make another label for the interprocess boundaries so
944: // that Triangle will respect them
945: // 3) We then throw all that label into the new overlap
946: //
947: // Alternative: Figure out explicitly which segments were refined, and then
948: // communicated the refinement over the old overlap. Use this info to locally
949: // construct the new overlap and flip to get a decent mesh
950: sieveCompletion::scatterCones(refSieve, refSieve, reMesh->getDistSendOverlap(), refMesh->getDistRecvOverlap(), refMesh);
951: #endif
952: return refMesh;
953: };
954: };
955: };
956: #endif
957: #ifdef PETSC_HAVE_TETGEN
958: namespace TetGen {
959: template<typename Mesh>
960: class Generator {
961: public:
962: static Obj<Mesh> generateMesh(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false) {
963: typedef ALE::SieveAlg<Mesh> sieve_alg_type;
964: const int dim = 3;
965: Obj<Mesh> mesh = new Mesh(boundary->comm(), dim, boundary->debug());
966: const PetscMPIInt rank = mesh->commRank();
967: bool createConvexHull = false;
968: ::tetgenio in;
969: ::tetgenio out;
971: const Obj<typename Mesh::label_sequence>& vertices = boundary->depthStratum(0);
972: const Obj<typename Mesh::numbering_type>& vNumbering = boundary->getFactory()->getLocalNumbering(boundary, 0);
973: const Obj<typename Mesh::real_section_type>& coordinates = boundary->getRealSection("coordinates");
974: const Obj<typename Mesh::label_type>& markers = boundary->getLabel("marker");
977: in.numberofpoints = vertices->size();
978: if (in.numberofpoints > 0) {
979: in.pointlist = new double[in.numberofpoints*dim];
980: in.pointmarkerlist = new int[in.numberofpoints];
981: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
982: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
983: const int idx = vNumbering->getIndex(*v_iter);
985: for(int d = 0; d < dim; d++) {
986: in.pointlist[idx*dim + d] = array[d];
987: }
988: in.pointmarkerlist[idx] = boundary->getValue(markers, *v_iter);
989: }
990: }
992: if (boundary->depth() != 0) { //our boundary mesh COULD be just a pointset; in which case depth = height = 0;
993: const Obj<typename Mesh::label_sequence>& facets = boundary->depthStratum(boundary->depth());
994: //PetscPrintf(boundary->comm(), "%d facets on the boundary\n", facets->size());
995: const Obj<typename Mesh::numbering_type>& fNumbering = boundary->getFactory()->getLocalNumbering(boundary, boundary->depth());
997: in.numberoffacets = facets->size();
998: if (in.numberoffacets > 0) {
999: in.facetlist = new tetgenio::facet[in.numberoffacets];
1000: in.facetmarkerlist = new int[in.numberoffacets];
1001: for(typename Mesh::label_sequence::iterator f_iter = facets->begin(); f_iter != facets->end(); ++f_iter) {
1002: const Obj<typename sieve_alg_type::coneArray>& cone = sieve_alg_type::nCone(boundary, *f_iter, boundary->depth());
1003: const int idx = fNumbering->getIndex(*f_iter);
1005: in.facetlist[idx].numberofpolygons = 1;
1006: in.facetlist[idx].polygonlist = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
1007: in.facetlist[idx].numberofholes = 0;
1008: in.facetlist[idx].holelist = NULL;
1010: tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
1011: int c = 0;
1013: poly->numberofvertices = cone->size();
1014: poly->vertexlist = new int[poly->numberofvertices];
1015: for(typename sieve_alg_type::coneArray::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
1016: const int vIdx = vNumbering->getIndex(*c_iter);
1018: poly->vertexlist[c++] = vIdx;
1019: }
1020: in.facetmarkerlist[idx] = boundary->getValue(markers, *f_iter);
1021: }
1022: }
1023: }else {
1024: createConvexHull = true;
1025: }
1027: in.numberofholes = 0;
1028: if (rank == 0) {
1029: // Normal operation
1030: std::string args("pqezQ");
1031: //constrained operation
1032: if (constrained) {
1033: args = "pezQ";
1034: if (createConvexHull) {
1035: args = "ezQ";
1036: //PetscPrintf(boundary->comm(), "createConvexHull\n");
1037: }
1038: }
1039: // Just make tetrahedrons
1040: // std::string args("efzV");
1041: // Adds a center point
1042: // std::string args("pqezQi");
1043: // in.numberofaddpoints = 1;
1044: // in.addpointlist = new double[in.numberofaddpoints*dim];
1045: // in.addpointlist[0] = 0.5;
1046: // in.addpointlist[1] = 0.5;
1047: // in.addpointlist[2] = 0.5;
1049: //if (createConvexHull) args += "c"; NOT SURE, but this was basically unused before, and the convex hull should be filled in if "p" isn't included.
1050: ::tetrahedralize((char *) args.c_str(), &in, &out);
1051: }
1052: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
1053: int numCorners = 4;
1054: int numCells = out.numberoftetrahedra;
1055: int *cells = out.tetrahedronlist;
1056: int numVertices = out.numberofpoints;
1057: double *coords = out.pointlist;
1059: if (!interpolate) {
1060: for(int c = 0; c < numCells; ++c) {
1061: int tmp = cells[c*4+0];
1062: cells[c*4+0] = cells[c*4+1];
1063: cells[c*4+1] = tmp;
1064: }
1065: }
1066: ALE::SieveBuilder<Mesh>::buildTopology(newSieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, mesh->getArrowSection("orientation"));
1067: mesh->setSieve(newSieve);
1068: mesh->stratify();
1069: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, dim, coords);
1070: const Obj<typename Mesh::label_type>& newMarkers = mesh->createLabel("marker");
1072: for(int v = 0; v < out.numberofpoints; v++) {
1073: if (out.pointmarkerlist[v]) {
1074: mesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
1075: }
1076: }
1077: if (interpolate) {
1078: if (out.edgemarkerlist) {
1079: for(int e = 0; e < out.numberofedges; e++) {
1080: if (out.edgemarkerlist[e]) {
1081: typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
1082: typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
1083: Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(endpointA, endpointB, 1);
1085: mesh->setValue(newMarkers, *edge->begin(), out.edgemarkerlist[e]);
1086: }
1087: }
1088: }
1089: if (out.trifacemarkerlist) {
1090: // Work around TetGen bug for raw tetrahedralization
1091: // The boundary faces are 0,1,4,5,8,9,11,12,13,15,16,17
1092: // for(int f = 0; f < out.numberoftrifaces; f++) {
1093: // if (out.trifacemarkerlist[f]) {
1094: // out.trifacemarkerlist[f] = 0;
1095: // } else {
1096: // out.trifacemarkerlist[f] = 1;
1097: // }
1098: // }
1099: for(int f = 0; f < out.numberoftrifaces; f++) {
1100: if (out.trifacemarkerlist[f]) {
1101: typename Mesh::point_type cornerA(out.trifacelist[f*3+0]+out.numberoftetrahedra);
1102: typename Mesh::point_type cornerB(out.trifacelist[f*3+1]+out.numberoftetrahedra);
1103: typename Mesh::point_type cornerC(out.trifacelist[f*3+2]+out.numberoftetrahedra);
1104: Obj<typename Mesh::sieve_type::supportSet> corners = typename Mesh::sieve_type::supportSet();
1105: Obj<typename Mesh::sieve_type::supportSet> edges = typename Mesh::sieve_type::supportSet();
1106: corners->clear();corners->insert(cornerA);corners->insert(cornerB);
1107: edges->insert(*newSieve->nJoin1(corners)->begin());
1108: corners->clear();corners->insert(cornerB);corners->insert(cornerC);
1109: edges->insert(*newSieve->nJoin1(corners)->begin());
1110: corners->clear();corners->insert(cornerC);corners->insert(cornerA);
1111: edges->insert(*newSieve->nJoin1(corners)->begin());
1112: const typename Mesh::point_type face = *newSieve->nJoin1(edges)->begin();
1113: const int faceMarker = out.trifacemarkerlist[f];
1114: const Obj<typename Mesh::coneArray> closure = sieve_alg_type::closure(mesh, face);
1115: const typename Mesh::coneArray::iterator end = closure->end();
1117: for(typename Mesh::coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
1118: mesh->setValue(newMarkers, *cl_iter, faceMarker);
1119: }
1120: }
1121: }
1122: }
1123: }
1124: return mesh;
1125: };
1128: static Obj<Mesh> generateMeshV(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false, const bool renumber = false) {
1129: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
1130: typedef typename Mesh::real_section_type::value_type real;
1131: const int dim = 3;
1132: Obj<Mesh> mesh = new Mesh(boundary->comm(), dim, boundary->debug());
1133: const Obj<typename Mesh::sieve_type>& sieve = boundary->getSieve();
1134: const PetscMPIInt rank = mesh->commRank();
1135: bool createConvexHull = false;
1136: PetscErrorCode ierr;
1137: ::tetgenio in;
1138: ::tetgenio out;
1140: const Obj<typename Mesh::label_sequence>& vertices = boundary->depthStratum(0);
1141: const Obj<typename Mesh::label_type>& markers = boundary->getLabel("marker");
1142: const Obj<typename Mesh::real_section_type>& coordinates = boundary->getRealSection("coordinates");
1143: const Obj<typename Mesh::numbering_type>& vNumbering = boundary->getFactory()->getLocalNumbering(boundary, 0);
1145: in.numberofpoints = vertices->size();
1146: if (in.numberofpoints > 0) {
1147: const typename Mesh::label_sequence::iterator vEnd = vertices->end();
1149: in.pointlist = new double[in.numberofpoints*dim];
1150: in.pointmarkerlist = new int[in.numberofpoints];
1151: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
1152: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
1153: const int idx = vNumbering->getIndex(*v_iter);
1155: for(int d = 0; d < dim; ++d) {
1156: in.pointlist[idx*dim + d] = array[d];
1157: }
1158: in.pointmarkerlist[idx] = boundary->getValue(markers, *v_iter);
1159: }
1160: }
1162: // Our boundary mesh COULD be just a pointset; in which case depth = height = 0;
1163: if (boundary->depth() != 0) {
1164: const Obj<typename Mesh::label_sequence>& facets = boundary->depthStratum(boundary->depth());
1165: const Obj<typename Mesh::numbering_type>& fNumbering = boundary->getFactory()->getLocalNumbering(boundary, boundary->depth());
1167: in.numberoffacets = facets->size();
1168: if (in.numberoffacets > 0) {
1169: ALE::ISieveVisitor::NConeRetriever<typename Mesh::sieve_type> ncV(*sieve, 5);
1170: const typename Mesh::label_sequence::iterator fEnd = facets->end();
1172: in.facetlist = new tetgenio::facet[in.numberoffacets];
1173: in.facetmarkerlist = new int[in.numberoffacets];
1174: for(typename Mesh::label_sequence::iterator f_iter = facets->begin(); f_iter != fEnd; ++f_iter) {
1175: ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*sieve, *f_iter, ncV);
1176: const int idx = fNumbering->getIndex(*f_iter);
1178: in.facetlist[idx].numberofpolygons = 1;
1179: in.facetlist[idx].polygonlist = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
1180: in.facetlist[idx].numberofholes = 0;
1181: in.facetlist[idx].holelist = NULL;
1183: tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
1184: const size_t n = ncV.getSize();
1185: const typename Mesh::point_type *cone = ncV.getPoints();
1187: poly->numberofvertices = n;
1188: poly->vertexlist = new int[poly->numberofvertices];
1189: for(size_t c = 0; c < n; ++c) {
1190: poly->vertexlist[c] = vNumbering->getIndex(cone[c]);
1191: }
1192: in.facetmarkerlist[idx] = boundary->getValue(markers, *f_iter);
1193: ncV.clear();
1194: }
1195: }
1196: } else {
1197: createConvexHull = true;
1198: }
1199: const typename Mesh::holes_type& holes = mesh->getHoles();
1201: in.numberofholes = holes.size();
1202: if (in.numberofholes > 0) {
1203: PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);CHKERRXX(ierr);
1204: for(int h = 0; h < in.numberofholes; ++h) {
1205: for(int d = 0; d < dim; ++d) {
1206: in.holelist[h*dim+d] = holes[h][d];
1207: }
1208: }
1209: }
1210: if (rank == 0) {
1211: // Normal operation
1212: std::string args("pqezQ");
1213: // Constrained operation
1214: if (constrained) {
1215: args = "pezQ";
1216: if (createConvexHull) {
1217: args = "ezQ";
1218: }
1219: }
1220: // Just make tetrahedrons
1221: // std::string args("efzV");
1222: // Adds a center point
1223: // std::string args("pqezQi");
1224: // in.numberofaddpoints = 1;
1225: // in.addpointlist = new double[in.numberofaddpoints*dim];
1226: // in.addpointlist[0] = 0.5;
1227: // in.addpointlist[1] = 0.5;
1228: // in.addpointlist[2] = 0.5;
1229: ::tetrahedralize((char *) args.c_str(), &in, &out);
1230: }
1231: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
1232: const Obj<FlexMesh> m = new FlexMesh(mesh->comm(), dim, mesh->debug());
1233: const Obj<FlexMesh::sieve_type> newS = new FlexMesh::sieve_type(m->comm(), m->debug());
1234: int numCorners = 4;
1235: int numCells = out.numberoftetrahedra;
1236: int *cells = out.tetrahedronlist;
1237: int numVertices = out.numberofpoints;
1238: double *coords = out.pointlist;
1239: real *coordsR;
1241: if (!interpolate) {
1242: // TetGen reports tetrahedra with the opposite orientation from what we expect
1243: for(int c = 0; c < numCells; ++c) {
1244: int tmp = cells[c*4+0];
1245: cells[c*4+0] = cells[c*4+1];
1246: cells[c*4+1] = tmp;
1247: }
1248: }
1249: ALE::SieveBuilder<FlexMesh>::buildTopology(newS, dim, numCells, cells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
1250: m->setSieve(newS);
1251: m->stratify();
1252: mesh->setSieve(newSieve);
1253: std::map<typename Mesh::point_type,typename Mesh::point_type> renumbering;
1254: ALE::ISieveConverter::convertSieve(*newS, *newSieve, renumbering, renumber);
1255: mesh->stratify();
1256: ALE::ISieveConverter::convertOrientation(*newS, *newSieve, renumbering, m->getArrowSection("orientation").ptr());
1257: {
1258: if (sizeof(double) == sizeof(real)) {
1259: coordsR = (real *) coords;
1260: } else {
1261: coordsR = new real[numVertices*dim];
1262: for(int i = 0; i < numVertices*dim; ++i) coordsR[i] = coords[i];
1263: }
1264: }
1265: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, dim, coordsR);
1266: {
1267: if (sizeof(double) != sizeof(real)) {
1268: delete [] coordsR;
1269: }
1270: }
1271: const Obj<typename Mesh::label_type>& newMarkers = mesh->createLabel("marker");
1273: for(int v = 0; v < out.numberofpoints; v++) {
1274: if (out.pointmarkerlist[v]) {
1275: if (renumber) {
1276: mesh->setValue(newMarkers, renumbering[v+out.numberoftetrahedra], out.pointmarkerlist[v]);
1277: } else {
1278: mesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
1279: }
1280: }
1281: }
1282: if (interpolate) {
1283: // This does not work anymore (edgemarkerlist is always empty). I tried -ee and it gave bogus results
1284: if (out.edgemarkerlist) {
1285: for(int e = 0; e < out.numberofedges; e++) {
1286: if (out.edgemarkerlist[e]) {
1287: typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
1288: typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
1289: Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(endpointA, endpointB, 1);
1291: if (renumber) {
1292: mesh->setValue(newMarkers, renumbering[*edge->begin()], out.edgemarkerlist[e]);
1293: } else {
1294: mesh->setValue(newMarkers, *edge->begin(), out.edgemarkerlist[e]);
1295: }
1296: }
1297: }
1298: }
1299: if (out.trifacemarkerlist) {
1300: // Work around TetGen bug for raw tetrahedralization
1301: // The boundary faces are 0,1,4,5,8,9,11,12,13,15,16,17
1302: // for(int f = 0; f < out.numberoftrifaces; f++) {
1303: // if (out.trifacemarkerlist[f]) {
1304: // out.trifacemarkerlist[f] = 0;
1305: // } else {
1306: // out.trifacemarkerlist[f] = 1;
1307: // }
1308: // }
1309: for(int f = 0; f < out.numberoftrifaces; f++) {
1310: if (out.trifacemarkerlist[f]) {
1311: typename Mesh::point_type cornerA = out.trifacelist[f*3+0]+out.numberoftetrahedra;
1312: typename Mesh::point_type cornerB = out.trifacelist[f*3+1]+out.numberoftetrahedra;
1313: typename Mesh::point_type cornerC = out.trifacelist[f*3+2]+out.numberoftetrahedra;
1314: Obj<typename Mesh::sieve_type::supportSet> corners = typename Mesh::sieve_type::supportSet();
1315: Obj<typename Mesh::sieve_type::supportSet> edges = typename Mesh::sieve_type::supportSet();
1316: corners->clear();corners->insert(cornerA);corners->insert(cornerB);
1317: typename Mesh::point_type edgeA = *newS->nJoin1(corners)->begin();
1318: corners->clear();corners->insert(cornerB);corners->insert(cornerC);
1319: typename Mesh::point_type edgeB = *newS->nJoin1(corners)->begin();
1320: corners->clear();corners->insert(cornerC);corners->insert(cornerA);
1321: typename Mesh::point_type edgeC = *newS->nJoin1(corners)->begin();
1322: edges->insert(edgeA); edges->insert(edgeB); edges->insert(edgeC);
1323: ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(30);
1324: typename Mesh::point_type face = *newS->nJoin1(edges)->begin();
1325: const int faceMarker = out.trifacemarkerlist[f];
1327: if (renumber) {face = renumbering[face];}
1328: ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*newSieve, face, pV);
1329: const size_t n = pV.getSize();
1330: const typename Mesh::point_type *cone = pV.getPoints();
1332: for(size_t c = 0; c < n; ++c) {
1333: mesh->setValue(newMarkers, cone[c], faceMarker);
1334: }
1335: pV.clear();
1336: }
1337: }
1338: }
1339: }
1340: mesh->copyHoles(boundary);
1341: return mesh;
1342: };
1343: };
1344: template<typename Mesh>
1345: class Refiner {
1346: public:
1347: static Obj<Mesh> refineMesh(const Obj<Mesh>& serialMesh, const double maxVolumes[], const bool interpolate = false) {
1348: typedef ALE::SieveAlg<Mesh> sieve_alg_type;
1349: const int dim = serialMesh->getDimension();
1350: const int depth = serialMesh->depth();
1351: const Obj<Mesh> refMesh = new Mesh(serialMesh->comm(), dim, serialMesh->debug());
1352: ::tetgenio in;
1353: ::tetgenio out;
1355: const Obj<typename Mesh::label_sequence>& vertices = serialMesh->depthStratum(0);
1356: const Obj<typename Mesh::label_type>& markers = serialMesh->getLabel("marker");
1357: const Obj<typename Mesh::real_section_type>& coordinates = serialMesh->getRealSection("coordinates");
1358: const Obj<typename Mesh::numbering_type>& vNumbering = serialMesh->getFactory()->getLocalNumbering(serialMesh, 0);
1360: in.numberofpoints = vertices->size();
1361: if (in.numberofpoints > 0) {
1362: in.pointlist = new double[in.numberofpoints*dim];
1363: in.pointmarkerlist = new int[in.numberofpoints];
1364: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
1365: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
1366: const int idx = vNumbering->getIndex(*v_iter);
1368: for(int d = 0; d < dim; d++) {
1369: in.pointlist[idx*dim + d] = array[d];
1370: }
1371: in.pointmarkerlist[idx] = serialMesh->getValue(markers, *v_iter);
1372: }
1373: }
1374: const Obj<typename Mesh::label_sequence>& cells = serialMesh->heightStratum(0);
1375: const Obj<typename Mesh::numbering_type>& cNumbering = serialMesh->getFactory()->getLocalNumbering(serialMesh, depth);
1377: in.numberofcorners = 4;
1378: in.numberoftetrahedra = cells->size();
1379: in.tetrahedronvolumelist = (double *) maxVolumes;
1380: if (in.numberoftetrahedra > 0) {
1381: in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
1382: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
1383: typedef ALE::SieveAlg<Mesh> sieve_alg_type;
1384: const Obj<typename sieve_alg_type::coneArray>& cone = sieve_alg_type::nCone(serialMesh, *c_iter, depth);
1385: const int idx = cNumbering->getIndex(*c_iter);
1386: int v = 0;
1388: for(typename Mesh::sieve_type::coneArray::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
1389: in.tetrahedronlist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*v_iter);
1390: }
1391: }
1392: }
1393: if (serialMesh->depth() == 3) {
1394: const Obj<typename Mesh::label_sequence>& boundary = serialMesh->getLabelStratum("marker", 1);
1396: in.numberoftrifaces = 0;
1397: for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != boundary->end(); ++b_iter) {
1398: if (serialMesh->height(*b_iter) == 1) {
1399: in.numberoftrifaces++;
1400: }
1401: }
1402: if (in.numberoftrifaces > 0) {
1403: int f = 0;
1405: in.trifacelist = new int[in.numberoftrifaces*3];
1406: in.trifacemarkerlist = new int[in.numberoftrifaces];
1407: for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != boundary->end(); ++b_iter) {
1408: if (serialMesh->height(*b_iter) == 1) {
1409: const Obj<typename Mesh::coneArray>& cone = sieve_alg_type::nCone(serialMesh, *b_iter, 2);
1410: int p = 0;
1412: for(typename Mesh::coneArray::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
1413: in.trifacelist[f*3 + (p++)] = vNumbering->getIndex(*v_iter);
1414: }
1415: in.trifacemarkerlist[f++] = serialMesh->getValue(markers, *b_iter);
1416: }
1417: }
1418: }
1419: }
1421: in.numberofholes = 0;
1422: if (serialMesh->commRank() == 0) {
1423: std::string args("qezQra");
1425: ::tetrahedralize((char *) args.c_str(), &in, &out);
1426: }
1427: in.tetrahedronvolumelist = NULL;
1428: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(refMesh->comm(), refMesh->debug());
1429: int numCorners = 4;
1430: int numCells = out.numberoftetrahedra;
1431: int *newCells = out.tetrahedronlist;
1432: int numVertices = out.numberofpoints;
1433: double *coords = out.pointlist;
1435: if (!interpolate) {
1436: for(int c = 0; c < numCells; ++c) {
1437: int tmp = newCells[c*4+0];
1438: newCells[c*4+0] = newCells[c*4+1];
1439: newCells[c*4+1] = tmp;
1440: }
1441: }
1442: ALE::SieveBuilder<Mesh>::buildTopology(newSieve, dim, numCells, newCells, numVertices, interpolate, numCorners, -1, refMesh->getArrowSection("orientation"));
1443: refMesh->setSieve(newSieve);
1444: refMesh->stratify();
1445: ALE::SieveBuilder<Mesh>::buildCoordinates(refMesh, dim, coords);
1446: const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");
1449: for(int v = 0; v < out.numberofpoints; v++) {
1450: if (out.pointmarkerlist[v]) {
1451: refMesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
1452: }
1453: }
1454: if (interpolate) {
1455: if (out.edgemarkerlist) {
1456: for(int e = 0; e < out.numberofedges; e++) {
1457: if (out.edgemarkerlist[e]) {
1458: typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
1459: typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
1460: Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(endpointA, endpointB, 1);
1462: refMesh->setValue(newMarkers, *edge->begin(), out.edgemarkerlist[e]);
1463: }
1464: }
1465: }
1466: if (out.trifacemarkerlist) {
1467: for(int f = 0; f < out.numberoftrifaces; f++) {
1468: if (out.trifacemarkerlist[f]) {
1469: typename Mesh::point_type cornerA(out.trifacelist[f*3+0]+out.numberoftetrahedra);
1470: typename Mesh::point_type cornerB(out.trifacelist[f*3+1]+out.numberoftetrahedra);
1471: typename Mesh::point_type cornerC(out.trifacelist[f*3+2]+out.numberoftetrahedra);
1472: Obj<typename Mesh::sieve_type::supportSet> corners = typename Mesh::sieve_type::supportSet();
1473: Obj<typename Mesh::sieve_type::supportSet> edges = typename Mesh::sieve_type::supportSet();
1474: corners->clear();corners->insert(cornerA);corners->insert(cornerB);
1475: edges->insert(*newSieve->nJoin1(corners)->begin());
1476: corners->clear();corners->insert(cornerB);corners->insert(cornerC);
1477: edges->insert(*newSieve->nJoin1(corners)->begin());
1478: corners->clear();corners->insert(cornerC);corners->insert(cornerA);
1479: edges->insert(*newSieve->nJoin1(corners)->begin());
1480: const typename Mesh::point_type face = *newSieve->nJoin1(edges)->begin();
1481: const int faceMarker = out.trifacemarkerlist[f];
1482: const Obj<typename Mesh::coneArray> closure = sieve_alg_type::closure(refMesh, face);
1483: const typename Mesh::coneArray::iterator end = closure->end();
1485: for(typename Mesh::coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
1486: refMesh->setValue(newMarkers, *cl_iter, faceMarker);
1487: }
1488: }
1489: }
1490: }
1491: }
1492: if (refMesh->commSize() > 1) {
1493: return ALE::Distribution<Mesh>::distributeMesh(refMesh);
1494: }
1495: return refMesh;
1496: };
1497: static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false) {
1498: Obj<Mesh> serialMesh = ALE::Distribution<Mesh>::unifyMesh(mesh);
1499: const Obj<typename Mesh::real_section_type> serialMaxVolumes = ALE::Distribution<Mesh>::distributeSection(maxVolumes, serialMesh, serialMesh->getDistSendOverlap(), serialMesh->getDistRecvOverlap());
1501: return refineMesh(serialMesh, serialMaxVolumes->restrictSpace(), interpolate);
1502: };
1503: static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false) {
1504: Obj<Mesh> serialMesh;
1505: if (mesh->commSize() > 1) {
1506: serialMesh = ALE::Distribution<Mesh>::unifyMesh(mesh);
1507: } else {
1508: serialMesh = mesh;
1509: }
1510: const int numCells = serialMesh->heightStratum(0)->size();
1511: double *serialMaxVolumes = new double[numCells];
1513: for(int c = 0; c < numCells; c++) {
1514: serialMaxVolumes[c] = maxVolume;
1515: }
1516: const Obj<Mesh> refMesh = refineMesh(serialMesh, serialMaxVolumes, interpolate);
1517: delete [] serialMaxVolumes;
1518: return refMesh;
1519: };
1520: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolumes[], const bool interpolate = false, const bool forceSerial = false, const bool renumber = false) {
1521: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
1522: typedef typename Mesh::real_section_type::value_type real;
1523: typedef typename Mesh::point_type point_type;
1524: const int dim = mesh->getDimension();
1525: const int depth = mesh->depth();
1526: const Obj<Mesh> refMesh = new Mesh(mesh->comm(), dim, mesh->debug());
1527: const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
1528: PetscErrorCode ierr;
1529: ::tetgenio in;
1530: ::tetgenio out;
1532: const Obj<typename Mesh::label_sequence>& vertices = mesh->depthStratum(0);
1533: const Obj<typename Mesh::label_type>& markers = mesh->getLabel("marker");
1534: const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
1535: const Obj<typename Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
1537: in.numberofpoints = vertices->size();
1538: if (in.numberofpoints > 0) {
1539: const typename Mesh::label_sequence::iterator vEnd = vertices->end();
1541: in.pointlist = new double[in.numberofpoints*dim];
1542: in.pointmarkerlist = new int[in.numberofpoints];
1543: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
1544: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
1545: const int idx = vNumbering->getIndex(*v_iter);
1547: for(int d = 0; d < dim; d++) {
1548: in.pointlist[idx*dim + d] = array[d];
1549: }
1550: in.pointmarkerlist[idx] = mesh->getValue(markers, *v_iter);
1551: }
1552: }
1553: const Obj<typename Mesh::label_sequence>& cells = mesh->heightStratum(0);
1554: const Obj<typename Mesh::numbering_type>& cNumbering = mesh->getFactory()->getLocalNumbering(mesh, depth);
1556: in.numberofcorners = 4;
1557: in.numberoftetrahedra = cells->size();
1558: in.tetrahedronvolumelist = (double *) maxVolumes;
1559: if (in.numberoftetrahedra > 0) {
1560: in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
1561: if (mesh->depth() == 1) {
1562: ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(4);
1563: const typename Mesh::label_sequence::iterator cEnd = cells->end();
1565: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
1566: sieve->cone(*c_iter, pV);
1567: const int idx = cNumbering->getIndex(*c_iter);
1568: const size_t n = pV.getSize();
1569: const point_type *cone = pV.getPoints();
1571: assert(n == 4);
1572: for(int v = 0; v < 4; ++v) {
1573: in.tetrahedronlist[idx*in.numberofcorners + v] = vNumbering->getIndex(cone[v]);
1574: }
1575: pV.clear();
1576: }
1577: } else if (mesh->depth() == 3) {
1578: // Need extra space due to early error checking
1579: ALE::ISieveVisitor::NConeRetriever<typename Mesh::sieve_type> ncV(*sieve, 5);
1580: const typename Mesh::label_sequence::iterator cEnd = cells->end();
1582: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
1583: ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
1584: const int idx = cNumbering->getIndex(*c_iter);
1585: const size_t n = ncV.getSize();
1586: const point_type *cone = ncV.getPoints();
1588: assert(n == 4);
1589: for(int v = 0; v < 4; ++v) {
1590: in.tetrahedronlist[idx*in.numberofcorners + v] = vNumbering->getIndex(cone[v]);
1591: }
1592: ncV.clear();
1593: }
1594: } else {
1595: throw ALE::Exception("Invalid sieve: Cannot gives sieves of arbitrary depth to TetGen");
1596: }
1597: }
1598: if (depth == 3) {
1599: const Obj<typename Mesh::label_sequence>& boundary = mesh->getLabelStratum("marker", 1);
1600: const typename Mesh::label_sequence::iterator bEnd = boundary->end();
1602: in.numberoftrifaces = 0;
1603: for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != bEnd; ++b_iter) {
1604: if (mesh->height(*b_iter) == 1) {
1605: in.numberoftrifaces++;
1606: }
1607: }
1608: if (in.numberoftrifaces > 0) {
1609: ALE::ISieveVisitor::NConeRetriever<typename Mesh::sieve_type> ncV(*sieve, 5);
1610: int f = 0;
1612: in.trifacelist = new int[in.numberoftrifaces*3];
1613: in.trifacemarkerlist = new int[in.numberoftrifaces];
1614: for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != bEnd; ++b_iter) {
1615: if (mesh->height(*b_iter) == 1) {
1616: ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*sieve, *b_iter, ncV);
1617: const size_t n = ncV.getSize();
1618: const point_type *cone = ncV.getPoints();
1620: for(size_t p = 0; p < n; ++p) {
1621: in.trifacelist[f*3 + p] = vNumbering->getIndex(cone[p]);
1622: }
1623: in.trifacemarkerlist[f++] = mesh->getValue(markers, *b_iter);
1624: ncV.clear();
1625: }
1626: }
1627: }
1628: }
1629: const typename Mesh::holes_type& holes = mesh->getHoles();
1631: in.numberofholes = holes.size();
1632: if (in.numberofholes > 0) {
1633: PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);CHKERRXX(ierr);
1634: for(int h = 0; h < in.numberofholes; ++h) {
1635: for(int d = 0; d < dim; ++d) {
1636: in.holelist[h*dim+d] = holes[h][d];
1637: }
1638: }
1639: }
1640: if (mesh->commRank() == 0) {
1641: std::string args("qezQra");
1643: ::tetrahedralize((char *) args.c_str(), &in, &out);
1644: }
1645: in.tetrahedronvolumelist = NULL;
1646: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(refMesh->comm(), refMesh->debug());
1647: const Obj<FlexMesh> m = new FlexMesh(mesh->comm(), dim, mesh->debug());
1648: const Obj<FlexMesh::sieve_type> newS = new FlexMesh::sieve_type(m->comm(), m->debug());
1649: int numCorners = 4;
1650: int numCells = out.numberoftetrahedra;
1651: int *newCells = out.tetrahedronlist;
1652: int numVertices = out.numberofpoints;
1653: double *coords = out.pointlist;
1654: real *coordsR;
1656: if (!interpolate) {
1657: for(int c = 0; c < numCells; ++c) {
1658: int tmp = newCells[c*4+0];
1659: newCells[c*4+0] = newCells[c*4+1];
1660: newCells[c*4+1] = tmp;
1661: }
1662: }
1663: ALE::SieveBuilder<FlexMesh>::buildTopology(newS, dim, numCells, newCells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
1664: m->setSieve(newS);
1665: m->stratify();
1666: refMesh->setSieve(newSieve);
1667: std::map<typename Mesh::point_type,typename Mesh::point_type> renumbering;
1668: ALE::ISieveConverter::convertSieve(*newS, *newSieve, renumbering, renumber);
1669: refMesh->stratify();
1670: ALE::ISieveConverter::convertOrientation(*newS, *newSieve, renumbering, m->getArrowSection("orientation").ptr());
1671: {
1672: if (sizeof(double) == sizeof(real)) {
1673: coordsR = (real *) coords;
1674: } else {
1675: coordsR = new real[numVertices*dim];
1676: for(int i = 0; i < numVertices*dim; ++i) coordsR[i] = coords[i];
1677: }
1678: }
1679: ALE::SieveBuilder<Mesh>::buildCoordinates(refMesh, dim, coordsR);
1680: {
1681: if (sizeof(double) != sizeof(real)) {
1682: delete [] coordsR;
1683: }
1684: }
1685: const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");
1687: for(int v = 0; v < out.numberofpoints; v++) {
1688: if (out.pointmarkerlist[v]) {
1689: if (renumber) {
1690: refMesh->setValue(newMarkers, renumbering[v+out.numberoftetrahedra], out.pointmarkerlist[v]);
1691: } else {
1692: refMesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
1693: }
1694: }
1695: }
1696: if (interpolate) {
1697: // This does not work anymore (edgemarkerlist is always empty). I tried -ee and it gave bogus results
1698: if (out.edgemarkerlist) {
1699: for(int e = 0; e < out.numberofedges; e++) {
1700: if (out.edgemarkerlist[e]) {
1701: typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
1702: typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
1703: Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(endpointA, endpointB, 1);
1705: if (renumber) {
1706: refMesh->setValue(newMarkers, renumbering[*edge->begin()], out.edgemarkerlist[e]);
1707: } else {
1708: refMesh->setValue(newMarkers, *edge->begin(), out.edgemarkerlist[e]);
1709: }
1710: }
1711: }
1712: }
1713: if (out.trifacemarkerlist) {
1714: for(int f = 0; f < out.numberoftrifaces; f++) {
1715: if (out.trifacemarkerlist[f]) {
1716: typename Mesh::point_type cornerA(out.trifacelist[f*3+0]+out.numberoftetrahedra);
1717: typename Mesh::point_type cornerB(out.trifacelist[f*3+1]+out.numberoftetrahedra);
1718: typename Mesh::point_type cornerC(out.trifacelist[f*3+2]+out.numberoftetrahedra);
1719: Obj<typename Mesh::sieve_type::supportSet> corners = typename Mesh::sieve_type::supportSet();
1720: Obj<typename Mesh::sieve_type::supportSet> edges = typename Mesh::sieve_type::supportSet();
1721: corners->clear();corners->insert(cornerA);corners->insert(cornerB);
1722: edges->insert(*newS->nJoin1(corners)->begin());
1723: corners->clear();corners->insert(cornerB);corners->insert(cornerC);
1724: edges->insert(*newS->nJoin1(corners)->begin());
1725: corners->clear();corners->insert(cornerC);corners->insert(cornerA);
1726: edges->insert(*newS->nJoin1(corners)->begin());
1727: ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(30);
1728: typename Mesh::point_type face = *newS->nJoin1(edges)->begin();
1729: const int faceMarker = out.trifacemarkerlist[f];
1731: if (renumber) {face = renumbering[face];}
1732: ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*newSieve, face, pV);
1733: const size_t n = pV.getSize();
1734: const typename Mesh::point_type *cone = pV.getPoints();
1736: for(size_t c = 0; c < n; ++c) {
1737: refMesh->setValue(newMarkers, cone[c], faceMarker);
1738: }
1739: pV.clear();
1740: }
1741: }
1742: }
1743: }
1744: return refMesh;
1745: };
1746: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false, const bool forceSerial = false, const bool renumber = false) {
1747: throw ALE::Exception("Not yet implemented");
1748: };
1749: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false, const bool forceSerial = false, const bool renumber = false) {
1750: const int numCells = mesh->heightStratum(0)->size();
1751: double *serialMaxVolumes = new double[numCells];
1753: for(int c = 0; c < numCells; c++) {
1754: serialMaxVolumes[c] = maxVolume;
1755: }
1756: const Obj<Mesh> refMesh = refineMeshV(mesh, serialMaxVolumes, interpolate, forceSerial, renumber);
1757: delete [] serialMaxVolumes;
1758: return refMesh;
1759: };
1760: };
1761: };
1762: #endif
1763: template<typename Mesh>
1764: class Generator {
1765: public:
1766: static Obj<Mesh> generateMesh(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false) {
1767: int dim = boundary->getDimension();
1769: if (dim == 1) {
1770: #ifdef PETSC_HAVE_TRIANGLE
1771: return ALE::Triangle::Generator<Mesh>::generateMesh(boundary, interpolate, constrained);
1772: #else
1773: throw ALE::Exception("Mesh generation currently requires Triangle to be installed. Use --download-triangle during configure.");
1774: #endif
1775: } else if (dim == 2) {
1776: #ifdef PETSC_HAVE_TETGEN
1777: return ALE::TetGen::Generator<Mesh>::generateMesh(boundary, interpolate, constrained);
1778: #else
1779: throw ALE::Exception("Mesh generation currently requires TetGen to be installed. Use --download-tetgen during configure.");
1780: #endif
1781: }
1782: return NULL;
1783: };
1784: static Obj<Mesh> generateMeshV(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false, const bool renumber = false) {
1785: int dim = boundary->getDimension();
1787: if (dim == 1) {
1788: #ifdef PETSC_HAVE_TRIANGLE
1789: return ALE::Triangle::Generator<Mesh>::generateMeshV(boundary, interpolate, constrained, renumber);
1790: #else
1791: throw ALE::Exception("Mesh generation currently requires Triangle to be installed. Use --download-triangle during configure.");
1792: #endif
1793: } else if (dim == 2) {
1794: #ifdef PETSC_HAVE_TETGEN
1795: return ALE::TetGen::Generator<Mesh>::generateMeshV(boundary, interpolate, constrained, renumber);
1796: #else
1797: throw ALE::Exception("Mesh generation currently requires TetGen to be installed. Use --download-tetgen during configure.");
1798: #endif
1799: }
1800: return NULL;
1801: };
1802: static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false) {
1803: int dim = mesh->getDimension();
1805: if (dim == 2) {
1806: #ifdef PETSC_HAVE_TRIANGLE
1807: return ALE::Triangle::Refiner<Mesh>::refineMesh(mesh, maxVolumes, interpolate);
1808: #else
1809: throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
1810: #endif
1811: } else if (dim == 3) {
1812: #ifdef PETSC_HAVE_TETGEN
1813: return ALE::TetGen::Refiner<Mesh>::refineMesh(mesh, maxVolumes, interpolate);
1814: #else
1815: throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
1816: #endif
1817: }
1818: return NULL;
1819: };
1820: static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false, const bool forceSerial = false) {
1821: int dim = mesh->getDimension();
1823: if (dim == 2) {
1824: #ifdef PETSC_HAVE_TRIANGLE
1825: return ALE::Triangle::Refiner<Mesh>::refineMesh(mesh, maxVolume, interpolate, forceSerial);
1826: #else
1827: throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
1828: #endif
1829: } else if (dim == 3) {
1830: #ifdef PETSC_HAVE_TETGEN
1831: return ALE::TetGen::Refiner<Mesh>::refineMesh(mesh, maxVolume, interpolate);
1832: #else
1833: throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
1834: #endif
1835: }
1836: return NULL;
1837: };
1838: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false, const bool renumber = false) {
1839: int dim = mesh->getDimension();
1841: if (dim == 2) {
1842: #ifdef PETSC_HAVE_TRIANGLE
1843: return ALE::Triangle::Refiner<Mesh>::refineMeshV(mesh, maxVolumes, interpolate, renumber);
1844: #else
1845: throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
1846: #endif
1847: } else if (dim == 3) {
1848: #ifdef PETSC_HAVE_TETGEN
1849: return ALE::TetGen::Refiner<Mesh>::refineMeshV(mesh, maxVolumes, interpolate, renumber);
1850: #else
1851: throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
1852: #endif
1853: }
1854: return NULL;
1855: };
1856: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false, const bool forceSerial = false, const bool renumber = false) {
1857: int dim = mesh->getDimension();
1859: if (dim == 2) {
1860: #ifdef PETSC_HAVE_TRIANGLE
1861: return ALE::Triangle::Refiner<Mesh>::refineMeshV(mesh, maxVolume, interpolate, forceSerial, renumber);
1862: #else
1863: throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
1864: #endif
1865: } else if (dim == 3) {
1866: #ifdef PETSC_HAVE_TETGEN
1867: return ALE::TetGen::Refiner<Mesh>::refineMeshV(mesh, maxVolume, interpolate, forceSerial, renumber);
1868: #else
1869: throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
1870: #endif
1871: }
1872: return NULL;
1873: };
1874: static Obj<Mesh> refineMeshLocal(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false) {
1875: int dim = mesh->getDimension();
1877: if (dim == 2) {
1878: #ifdef PETSC_HAVE_TRIANGLE
1879: return ALE::Triangle::Refiner<Mesh>::refineMeshLocal(mesh, maxVolume, interpolate);
1880: #else
1881: throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
1882: #endif
1883: } else if (dim == 3) {
1884: #ifdef PETSC_HAVE_TETGEN
1885: return ALE::TetGen::Refiner<Mesh>::refineMesh(mesh, maxVolume, interpolate);
1886: #else
1887: throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
1888: #endif
1889: }
1890: return NULL;
1891: };
1892: };
1893: }
1895: #endif