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 <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;
85: PetscErrorCode ierr;
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);
99: PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);
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);
118: PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);
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;
124:
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);
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);}
155: if (in.pointmarkerlist) {PetscFree(in.pointmarkerlist);}
156: if (in.segmentlist) {PetscFree(in.segmentlist);}
157: if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);}
158: if (in.holelist) {PetscFree(in.holelist);}
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: int dim = 2;
198: const Obj<Mesh> mesh = new Mesh(boundary->comm(), dim, boundary->debug());
199: const Obj<typename Mesh::sieve_type>& sieve = boundary->getSieve();
200: const bool createConvexHull = false;
201: struct triangulateio in;
202: struct triangulateio out;
203: PetscErrorCode ierr;
205: initInput(&in);
206: initOutput(&out);
207: const Obj<typename Mesh::label_sequence>& vertices = boundary->depthStratum(0);
208: const Obj<typename Mesh::label_type>& markers = boundary->getLabel("marker");
209: const Obj<typename Mesh::real_section_type>& coordinates = boundary->getRealSection("coordinates");
210: const Obj<typename Mesh::numbering_type>& vNumbering = boundary->getFactory()->getLocalNumbering(boundary, 0);
212: in.numberofpoints = vertices->size();
213: if (in.numberofpoints > 0) {
214: const typename Mesh::label_sequence::iterator vEnd = vertices->end();
216: PetscMalloc(in.numberofpoints * dim * sizeof(double), &in.pointlist);
217: PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);
218: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
219: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
220: const int idx = vNumbering->getIndex(*v_iter);
222: for(int d = 0; d < dim; ++d) {
223: in.pointlist[idx*dim + d] = array[d];
224: }
225: in.pointmarkerlist[idx] = boundary->getValue(markers, *v_iter);
226: }
227: }
228: const Obj<typename Mesh::label_sequence>& edges = boundary->depthStratum(1);
229: const Obj<typename Mesh::numbering_type>& eNumbering = boundary->getFactory()->getLocalNumbering(boundary, 1);
231: in.numberofsegments = edges->size();
232: if (in.numberofsegments > 0) {
233: const typename Mesh::label_sequence::iterator eEnd = edges->end();
235: PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);
236: PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);
237: SegmentVisitor sV(dim, in.segmentlist, *vNumbering);
238: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != eEnd; ++e_iter) {
239: const int idx = eNumbering->getIndex(*e_iter);
241: sV.setIndex(idx);
242: sieve->cone(*e_iter, sV);
243: in.segmentmarkerlist[idx] = boundary->getValue(markers, *e_iter);
244: }
245: }
246: const typename Mesh::holes_type& holes = boundary->getHoles();
248: in.numberofholes = holes.size();
249: if (in.numberofholes > 0) {
250: PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);
251: for(int h = 0; h < in.numberofholes; ++h) {
252: for(int d = 0; d < dim; ++d) {
253: in.holelist[h*dim+d] = holes[h][d];
254: }
255: }
256: }
257: if (mesh->commRank() == 0) {
258: std::string args("pqezQ");
260: if (createConvexHull) {
261: args += "c";
262: }
263: if (constrained) {
264: args = "zepDQ";
265: }
266: triangulate((char *) args.c_str(), &in, &out, NULL);
267: }
268: if (in.pointlist) {PetscFree(in.pointlist);}
269: if (in.pointmarkerlist) {PetscFree(in.pointmarkerlist);}
270: if (in.segmentlist) {PetscFree(in.segmentlist);}
271: if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);}
272: if (in.holelist) {PetscFree(in.holelist);}
273: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
274: const Obj<ALE::Mesh> m = new ALE::Mesh(boundary->comm(), dim, boundary->debug());
275: const Obj<ALE::Mesh::sieve_type> newS = new ALE::Mesh::sieve_type(m->comm(), m->debug());
276: int numCorners = 3;
277: int numCells = out.numberoftriangles;
278: int *cells = out.trianglelist;
279: int numVertices = out.numberofpoints;
280: double *coords = out.pointlist;
282: ALE::SieveBuilder<ALE::Mesh>::buildTopology(newS, dim, numCells, cells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
283: m->setSieve(newS);
284: m->stratify();
285: mesh->setSieve(newSieve);
286: std::map<typename Mesh::point_type,typename Mesh::point_type> renumbering;
287: ALE::ISieveConverter::convertSieve(*newS, *newSieve, renumbering, false);
288: mesh->stratify();
289: ALE::ISieveConverter::convertOrientation(*newS, *newSieve, renumbering, m->getArrowSection("orientation").ptr());
290: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, dim, coords);
291: const Obj<typename Mesh::label_type>& newMarkers = mesh->createLabel("marker");
293: if (mesh->commRank() == 0) {
294: #ifdef IMESH_NEW_LABELS
295: int size = 0;
297: newMarkers->setChart(mesh->getSieve()->getChart());
298: for(int v = 0; v < out.numberofpoints; v++) {
299: if (out.pointmarkerlist[v]) {
300: newMarkers->setConeSize(v+out.numberoftriangles, 1);
301: size++;
302: }
303: }
304: if (interpolate) {
305: for(int e = 0; e < out.numberofedges; e++) {
306: if (out.edgemarkerlist[e]) {
307: const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
308: const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
309: const Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(vertexA, vertexB, 1);
311: newMarkers->setConeSize(*edge->begin(), 1);
312: size++;
313: }
314: }
315: }
316: newMarkers->setSupportSize(0, size);
317: newMarkers->allocate();
318: #endif
319: for(int v = 0; v < out.numberofpoints; v++) {
320: if (out.pointmarkerlist[v]) {
321: mesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
322: }
323: }
324: if (interpolate) {
325: for(int e = 0; e < out.numberofedges; e++) {
326: if (out.edgemarkerlist[e]) {
327: const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
328: const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
329: const Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(vertexA, vertexB, 1);
331: mesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
332: }
333: }
334: }
335: #ifdef IMESH_NEW_LABELS
336: newMarkers->recalculateLabel();
337: #endif
338: }
339: mesh->copyHoles(boundary);
340: finiOutput(&out);
341: return mesh;
342: };
343: };
344: template<typename Mesh>
345: class Refiner {
346: public:
347: static Obj<Mesh> refineMesh(const Obj<Mesh>& serialMesh, const double maxVolumes[], const bool interpolate = false, const bool forceSerial = false) {
348: const int dim = serialMesh->getDimension();
349: const Obj<Mesh> refMesh = new Mesh(serialMesh->comm(), dim, serialMesh->debug());
350: const Obj<typename Mesh::sieve_type>& serialSieve = serialMesh->getSieve();
351: struct triangulateio in;
352: struct triangulateio out;
353: PetscErrorCode ierr;
355: Generator<Mesh>::initInput(&in);
356: Generator<Mesh>::initOutput(&out);
357: const Obj<typename Mesh::label_sequence>& vertices = serialMesh->depthStratum(0);
358: const Obj<typename Mesh::label_type>& markers = serialMesh->getLabel("marker");
359: const Obj<typename Mesh::real_section_type>& coordinates = serialMesh->getRealSection("coordinates");
360: const Obj<typename Mesh::numbering_type>& vNumbering = serialMesh->getFactory()->getLocalNumbering(serialMesh, 0);
362: in.numberofpoints = vertices->size();
363: if (in.numberofpoints > 0) {
364: const typename Mesh::label_sequence::iterator vEnd = vertices->end();
366: PetscMalloc(in.numberofpoints * dim * sizeof(double), &in.pointlist);
367: PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);
368: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
369: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
370: const int idx = vNumbering->getIndex(*v_iter);
372: for(int d = 0; d < dim; d++) {
373: in.pointlist[idx*dim + d] = array[d];
374: }
375: in.pointmarkerlist[idx] = serialMesh->getValue(markers, *v_iter);
376: }
377: }
378: const Obj<typename Mesh::label_sequence>& faces = serialMesh->heightStratum(0);
379: const Obj<typename Mesh::numbering_type>& fNumbering = serialMesh->getFactory()->getLocalNumbering(serialMesh, serialMesh->depth());
381: in.numberofcorners = 3;
382: in.numberoftriangles = faces->size();
383: in.trianglearealist = (double *) maxVolumes;
384: if (in.numberoftriangles > 0) {
385: const typename Mesh::label_sequence::iterator fEnd = faces->end();
387: PetscMalloc(in.numberoftriangles*in.numberofcorners * sizeof(int), &in.trianglelist);
388: if (serialMesh->depth() == 1) {
389: for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
390: const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = serialSieve->cone(*f_iter);
391: const typename Mesh::sieve_type::traits::coneSequence::iterator cEnd = cone->end();
392: const int idx = fNumbering->getIndex(*f_iter);
393: int v = 0;
395: for(typename Mesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
396: in.trianglelist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*c_iter);
397: }
398: }
399: } else if (serialMesh->depth() == 2) {
400: for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
401: typedef ALE::SieveAlg<Mesh> sieve_alg_type;
402: const Obj<typename sieve_alg_type::coneArray>& cone = sieve_alg_type::nCone(serialMesh, *f_iter, 2);
403: const typename Mesh::sieve_type::coneArray::iterator cEnd = cone->end();
404: const int idx = fNumbering->getIndex(*f_iter);
405: int v = 0;
407: for(typename Mesh::sieve_type::coneArray::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
408: in.trianglelist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*c_iter);
409: }
410: }
411: } else {
412: throw ALE::Exception("Invalid sieve: Cannot gives sieves of arbitrary depth to Triangle");
413: }
414: }
415: if (serialMesh->depth() == 2) {
416: const Obj<typename Mesh::label_sequence>& edges = serialMesh->depthStratum(1);
417: #define NEW_LABEL
418: #ifdef NEW_LABEL
419: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
420: if (serialMesh->getValue(markers, *e_iter)) {
421: in.numberofsegments++;
422: }
423: }
424: std::cout << "Number of segments: " << in.numberofsegments << std::endl;
425: if (in.numberofsegments > 0) {
426: int s = 0;
428: PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);
429: PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);
430: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
431: const int edgeMarker = serialMesh->getValue(markers, *e_iter);
433: if (edgeMarker) {
434: const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = serialSieve->cone(*e_iter);
435: int p = 0;
437: for(typename Mesh::sieve_type::traits::coneSequence::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
438: in.segmentlist[s*2 + (p++)] = vNumbering->getIndex(*v_iter);
439: }
440: in.segmentmarkerlist[s++] = edgeMarker;
441: }
442: }
443: }
444: #else
445: const Obj<typename Mesh::label_type::baseSequence>& boundary = markers->base();
447: in.numberofsegments = 0;
448: for(typename Mesh::label_type::baseSequence::iterator b_iter = boundary->begin(); b_iter != boundary->end(); ++b_iter) {
449: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
450: if (*b_iter == *e_iter) {
451: in.numberofsegments++;
452: }
453: }
454: }
455: if (in.numberofsegments > 0) {
456: int s = 0;
458: PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);
459: PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);
460: for(typename Mesh::label_type::baseSequence::iterator b_iter = boundary->begin(); b_iter != boundary->end(); ++b_iter) {
461: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
462: if (*b_iter == *e_iter) {
463: const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = serialSieve->cone(*e_iter);
464: int p = 0;
466: for(typename Mesh::sieve_type::traits::coneSequence::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
467: in.segmentlist[s*2 + (p++)] = vNumbering->getIndex(*v_iter);
468: }
469: in.segmentmarkerlist[s++] = serialMesh->getValue(markers, *e_iter);
470: }
471: }
472: }
473: }
474: #endif
475: }
477: in.numberofholes = 0;
478: if (in.numberofholes > 0) {
479: PetscMalloc(in.numberofholes * dim * sizeof(int), &in.holelist);
480: }
481: if (serialMesh->commRank() == 0) {
482: std::string args("pqezQra");
484: triangulate((char *) args.c_str(), &in, &out, NULL);
485: }
486: if (in.pointlist) {PetscFree(in.pointlist);}
487: if (in.pointmarkerlist) {PetscFree(in.pointmarkerlist);}
488: if (in.segmentlist) {PetscFree(in.segmentlist);}
489: if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);}
490: if (in.trianglelist) {PetscFree(in.trianglelist);}
491: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(serialMesh->comm(), serialMesh->debug());
492: int numCorners = 3;
493: int numCells = out.numberoftriangles;
494: int *cells = out.trianglelist;
495: int numVertices = out.numberofpoints;
496: double *coords = out.pointlist;
498: ALE::SieveBuilder<Mesh>::buildTopology(newSieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, refMesh->getArrowSection("orientation"));
499: refMesh->setSieve(newSieve);
500: refMesh->stratify();
501: ALE::SieveBuilder<Mesh>::buildCoordinates(refMesh, dim, coords);
502: const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");
504: if (refMesh->commRank() == 0) {
505: for(int v = 0; v < out.numberofpoints; v++) {
506: if (out.pointmarkerlist[v]) {
507: refMesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
508: }
509: }
510: if (interpolate) {
511: for(int e = 0; e < out.numberofedges; e++) {
512: if (out.edgemarkerlist[e]) {
513: const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
514: const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
515: const Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(vertexA, vertexB, 1);
517: refMesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
518: }
519: }
520: }
521: }
523: Generator<Mesh>::finiOutput(&out);
524: if ((refMesh->commSize() > 1) && (!forceSerial)) {
525: return ALE::Distribution<Mesh>::distributeMesh(refMesh);
526: }
527: return refMesh;
528: };
529: static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false) {
530: Obj<Mesh> serialMesh = ALE::Distribution<Mesh>::unifyMesh(mesh);
531: const Obj<typename Mesh::real_section_type> serialMaxVolumes = ALE::Distribution<Mesh>::distributeSection(maxVolumes, serialMesh, serialMesh->getDistSendOverlap(), serialMesh->getDistRecvOverlap());
533: return refineMesh(serialMesh, serialMaxVolumes->restrictSpace(), interpolate);
534: };
535: static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false, const bool forceSerial = false) {
536: Obj<Mesh> serialMesh;
537: if ((mesh->commSize() > 1) && (!forceSerial)) {
538: serialMesh = ALE::Distribution<Mesh>::unifyMesh(mesh);
539: } else {
540: serialMesh = mesh;
541: }
542: const int numFaces = serialMesh->heightStratum(0)->size();
543: double *serialMaxVolumes = new double[numFaces];
545: for(int f = 0; f < numFaces; f++) {
546: serialMaxVolumes[f] = maxVolume;
547: }
548: const Obj<Mesh> refMesh = refineMesh(serialMesh, serialMaxVolumes, interpolate, forceSerial);
549: delete [] serialMaxVolumes;
550: return refMesh;
551: };
552: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolumes[], const bool interpolate = false, const bool forceSerial = false) {
553: typedef typename Mesh::point_type point_type;
554: const int dim = mesh->getDimension();
555: const Obj<Mesh> refMesh = new Mesh(mesh->comm(), dim, mesh->debug());
556: const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
557: struct triangulateio in;
558: struct triangulateio out;
559: PetscErrorCode ierr;
561: Generator<Mesh>::initInput(&in);
562: Generator<Mesh>::initOutput(&out);
563: const Obj<typename Mesh::label_sequence>& vertices = mesh->depthStratum(0);
564: const Obj<typename Mesh::label_type>& markers = mesh->getLabel("marker");
565: const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
566: const Obj<typename Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
568: in.numberofpoints = vertices->size();
569: if (in.numberofpoints > 0) {
570: const typename Mesh::label_sequence::iterator vEnd = vertices->end();
572: PetscMalloc(in.numberofpoints * dim * sizeof(double), &in.pointlist);
573: PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);
574: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
575: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
576: const int idx = vNumbering->getIndex(*v_iter);
578: for(int d = 0; d < dim; d++) {
579: in.pointlist[idx*dim + d] = array[d];
580: }
581: in.pointmarkerlist[idx] = mesh->getValue(markers, *v_iter);
582: }
583: }
584: const Obj<typename Mesh::label_sequence>& faces = mesh->heightStratum(0);
585: const Obj<typename Mesh::numbering_type>& fNumbering = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());
587: in.numberofcorners = 3;
588: in.numberoftriangles = faces->size();
589: in.trianglearealist = (double *) maxVolumes;
590: if (in.numberoftriangles > 0) {
591: PetscMalloc(in.numberoftriangles*in.numberofcorners * sizeof(int), &in.trianglelist);
592: if (mesh->depth() == 1) {
593: ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(3);
594: const typename Mesh::label_sequence::iterator fEnd = faces->end();
596: for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
597: sieve->cone(*f_iter, pV);
598: const int idx = fNumbering->getIndex(*f_iter);
599: const size_t n = pV.getSize();
600: const point_type *cone = pV.getPoints();
602: assert(n == 3);
603: for(int v = 0; v < 3; ++v) {
604: in.trianglelist[idx*in.numberofcorners + v] = vNumbering->getIndex(cone[v]);
605: }
606: pV.clear();
607: }
608: } else if (mesh->depth() == 2) {
609: // Need extra space due to early error checking
610: ALE::ISieveVisitor::NConeRetriever<typename Mesh::sieve_type> ncV(*sieve, 4);
611: const typename Mesh::label_sequence::iterator fEnd = faces->end();
613: for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
614: ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*sieve, *f_iter, ncV);
615: const int idx = fNumbering->getIndex(*f_iter);
616: const size_t n = ncV.getSize();
617: const point_type *cone = ncV.getPoints();
619: assert(n == 3);
620: for(int v = 0; v < 3; ++v) {
621: in.trianglelist[idx*in.numberofcorners + v] = vNumbering->getIndex(cone[v]);
622: }
623: ncV.clear();
624: }
625: } else {
626: throw ALE::Exception("Invalid sieve: Cannot gives sieves of arbitrary depth to Triangle");
627: }
628: }
629: if (mesh->depth() == 2) {
630: const Obj<typename Mesh::label_sequence>& edges = mesh->depthStratum(1);
631: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
632: if (mesh->getValue(markers, *e_iter)) {
633: in.numberofsegments++;
634: }
635: }
636: std::cout << "Number of segments: " << in.numberofsegments << std::endl;
637: if (in.numberofsegments > 0) {
638: ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(2);
639: int s = 0;
641: PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);
642: PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);
643: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
644: const int edgeMarker = mesh->getValue(markers, *e_iter);
646: if (edgeMarker) {
647: sieve->cone(*e_iter, pV);
648: const size_t n = pV.getSize();
649: const point_type *cone = pV.getPoints();
651: assert(n == 2);
652: for(int v = 0; v < 2; ++v) {
653: in.segmentlist[s*2 + v] = vNumbering->getIndex(cone[v]);
654: }
655: in.segmentmarkerlist[s++] = edgeMarker;
656: pV.clear();
657: }
658: }
659: }
660: }
661: const typename Mesh::holes_type& holes = mesh->getHoles();
663: in.numberofholes = holes.size();
664: if (in.numberofholes > 0) {
665: PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);
666: for(int h = 0; h < in.numberofholes; ++h) {
667: for(int d = 0; d < dim; ++d) {
668: in.holelist[h*dim+d] = holes[h][d];
669: }
670: }
671: }
672: if (mesh->commRank() == 0) {
673: std::string args("pqezQra");
675: triangulate((char *) args.c_str(), &in, &out, NULL);
676: }
677: if (in.pointlist) {PetscFree(in.pointlist);}
678: if (in.pointmarkerlist) {PetscFree(in.pointmarkerlist);}
679: if (in.segmentlist) {PetscFree(in.segmentlist);}
680: if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);}
681: if (in.trianglelist) {PetscFree(in.trianglelist);}
682: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
683: const Obj<ALE::Mesh> m = new ALE::Mesh(mesh->comm(), dim, mesh->debug());
684: const Obj<ALE::Mesh::sieve_type> newS = new ALE::Mesh::sieve_type(m->comm(), m->debug());
685: int numCorners = 3;
686: int numCells = out.numberoftriangles;
687: int *cells = out.trianglelist;
688: int numVertices = out.numberofpoints;
689: double *coords = out.pointlist;
691: ALE::SieveBuilder<ALE::Mesh>::buildTopology(newS, dim, numCells, cells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
692: m->setSieve(newS);
693: m->stratify();
694: refMesh->setSieve(newSieve);
695: std::map<typename Mesh::point_type,typename Mesh::point_type> renumbering;
696: ALE::ISieveConverter::convertSieve(*newS, *newSieve, renumbering, false);
697: refMesh->stratify();
698: ALE::ISieveConverter::convertOrientation(*newS, *newSieve, renumbering, m->getArrowSection("orientation").ptr());
699: ALE::SieveBuilder<Mesh>::buildCoordinates(refMesh, dim, coords);
700: const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");
702: if (refMesh->commRank() == 0) {
703: for(int v = 0; v < out.numberofpoints; v++) {
704: if (out.pointmarkerlist[v]) {
705: refMesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
706: }
707: }
708: if (interpolate) {
709: for(int e = 0; e < out.numberofedges; e++) {
710: if (out.edgemarkerlist[e]) {
711: const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
712: const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
713: const Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(vertexA, vertexB, 1);
715: refMesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
716: }
717: }
718: }
719: }
720: refMesh->copyHoles(mesh);
721: Generator<Mesh>::finiOutput(&out);
722: #if 0
723: if ((refMesh->commSize() > 1) && (!forceSerial)) {
724: return ALE::Distribution<Mesh>::distributeMesh(refMesh);
725: }
726: #endif
727: return refMesh;
728: };
729: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false) {
730: throw ALE::Exception("Not yet implemented");
731: };
732: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false, const bool forceSerial = false) {
733: #if 0
734: Obj<Mesh> serialMesh;
735: if (mesh->commSize() > 1) {
736: serialMesh = ALE::Distribution<Mesh>::unifyMesh(mesh);
737: } else {
738: serialMesh = mesh;
739: }
740: #endif
741: const int numCells = mesh->heightStratum(0)->size();
742: double *serialMaxVolumes = new double[numCells];
744: for(int c = 0; c < numCells; c++) {
745: serialMaxVolumes[c] = maxVolume;
746: }
747: const Obj<Mesh> refMesh = refineMeshV(mesh, serialMaxVolumes, interpolate);
748: delete [] serialMaxVolumes;
749: return refMesh;
750: };
751: static Obj<Mesh> refineMeshLocal(const Obj<Mesh>& mesh, const double maxVolumes[], const bool interpolate = false) {
752: const int dim = mesh->getDimension();
753: const Obj<Mesh> refMesh = new Mesh(mesh->comm(), dim, mesh->debug());
754: const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
755: struct triangulateio in;
756: struct triangulateio out;
757: PetscErrorCode ierr;
759: Generator<Mesh>::initInput(&in);
760: Generator<Mesh>::initOutput(&out);
761: const Obj<typename Mesh::label_sequence>& vertices = mesh->depthStratum(0);
762: const Obj<typename Mesh::label_type>& markers = mesh->getLabel("marker");
763: const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
764: const Obj<typename Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
766: in.numberofpoints = vertices->size();
767: if (in.numberofpoints > 0) {
768: PetscMalloc(in.numberofpoints * dim * sizeof(double), &in.pointlist);
769: PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);
770: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
771: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
772: const int idx = vNumbering->getIndex(*v_iter);
774: for(int d = 0; d < dim; d++) {
775: in.pointlist[idx*dim + d] = array[d];
776: }
777: in.pointmarkerlist[idx] = mesh->getValue(markers, *v_iter);
778: }
779: }
780: const Obj<typename Mesh::label_sequence>& faces = mesh->heightStratum(0);
781: const Obj<typename Mesh::numbering_type>& fNumbering = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());
783: in.numberofcorners = 3;
784: in.numberoftriangles = faces->size();
785: in.trianglearealist = (double *) maxVolumes;
786: if (in.numberoftriangles > 0) {
787: PetscMalloc(in.numberoftriangles*in.numberofcorners * sizeof(int), &in.trianglelist);
788: if (mesh->depth() == 1) {
789: for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != faces->end(); ++f_iter) {
790: const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = sieve->cone(*f_iter);
791: const int idx = fNumbering->getIndex(*f_iter);
792: int v = 0;
794: for(typename Mesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
795: in.trianglelist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*c_iter);
796: }
797: }
798: } else if (mesh->depth() == 2) {
799: for(typename Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != faces->end(); ++f_iter) {
800: typedef ALE::SieveAlg<Mesh> sieve_alg_type;
801: const Obj<typename sieve_alg_type::coneArray>& cone = sieve_alg_type::nCone(mesh, *f_iter, 2);
802: const int idx = fNumbering->getIndex(*f_iter);
803: int v = 0;
805: for(typename Mesh::sieve_type::coneArray::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
806: in.trianglelist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*c_iter);
807: }
808: }
809: } else {
810: throw ALE::Exception("Invalid sieve: Cannot gives sieves of arbitrary depth to Triangle");
811: }
812: }
813: if (mesh->depth() == 2) {
814: const Obj<typename Mesh::label_sequence>& edges = mesh->depthStratum(1);
815: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
816: if (mesh->getValue(markers, *e_iter)) {
817: in.numberofsegments++;
818: }
819: }
820: std::cout << "Number of segments: " << in.numberofsegments << std::endl;
821: if (in.numberofsegments > 0) {
822: int s = 0;
824: PetscMalloc(in.numberofsegments * 2 * sizeof(int), &in.segmentlist);
825: PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);
826: for(typename Mesh::label_sequence::iterator e_iter = edges->begin(); e_iter != edges->end(); ++e_iter) {
827: const int edgeMarker = mesh->getValue(markers, *e_iter);
829: if (edgeMarker) {
830: const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = sieve->cone(*e_iter);
831: int p = 0;
833: for(typename Mesh::sieve_type::traits::coneSequence::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
834: in.segmentlist[s*2 + (p++)] = vNumbering->getIndex(*v_iter);
835: }
836: in.segmentmarkerlist[s++] = edgeMarker;
837: }
838: }
839: }
840: }
842: in.numberofholes = 0;
843: if (in.numberofholes > 0) {
844: PetscMalloc(in.numberofholes * dim * sizeof(int), &in.holelist);
845: }
846: std::string args("pqezQra");
848: triangulate((char *) args.c_str(), &in, &out, NULL);
849: if (in.pointlist) {PetscFree(in.pointlist);}
850: if (in.pointmarkerlist) {PetscFree(in.pointmarkerlist);}
851: if (in.segmentlist) {PetscFree(in.segmentlist);}
852: if (in.segmentmarkerlist) {PetscFree(in.segmentmarkerlist);}
853: if (in.trianglelist) {PetscFree(in.trianglelist);}
854: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
855: int numCorners = 3;
856: int numCells = out.numberoftriangles;
857: int *cells = out.trianglelist;
858: int numVertices = out.numberofpoints;
859: double *coords = out.pointlist;
861: ALE::SieveBuilder<Mesh>::buildTopologyMultiple(newSieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, refMesh->getArrowSection("orientation"));
862: refMesh->setSieve(newSieve);
863: refMesh->stratify();
864: ALE::SieveBuilder<Mesh>::buildCoordinatesMultiple(refMesh, dim, coords);
865: const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");
867: for(int v = 0; v < out.numberofpoints; v++) {
868: if (out.pointmarkerlist[v]) {
869: refMesh->setValue(newMarkers, v+out.numberoftriangles, out.pointmarkerlist[v]);
870: }
871: }
872: if (interpolate) {
873: for(int e = 0; e < out.numberofedges; e++) {
874: if (out.edgemarkerlist[e]) {
875: const typename Mesh::point_type vertexA(out.edgelist[e*2+0]+out.numberoftriangles);
876: const typename Mesh::point_type vertexB(out.edgelist[e*2+1]+out.numberoftriangles);
877: const Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(vertexA, vertexB, 1);
879: refMesh->setValue(newMarkers, *(edge->begin()), out.edgemarkerlist[e]);
880: }
881: }
882: }
883: Generator<Mesh>::finiOutput(&out);
884: return refMesh;
885: };
886: static Obj<Mesh> refineMeshLocal(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false) {
887: const int numLocalFaces = mesh->heightStratum(0)->size();
888: double *localMaxVolumes = new double[numLocalFaces];
890: for(int f = 0; f < numLocalFaces; f++) {
891: localMaxVolumes[f] = maxVolume;
892: }
893: const Obj<Mesh> refMesh = refineMeshLocal(mesh, localMaxVolumes, interpolate);
894: const Obj<typename Mesh::sieve_type> refSieve = refMesh->getSieve();
895: delete [] localMaxVolumes;
896: #if 0
897: typedef typename ALE::New::Completion<Mesh, typename Mesh::sieve_type::point_type> sieveCompletion;
898: // This is where we enforce consistency over the overlap
899: // We need somehow to update the overlap to account for the new stuff
900: //
901: // 1) Since we are refining only, the vertices are invariant
902: // 2) We need to make another label for the interprocess boundaries so
903: // that Triangle will respect them
904: // 3) We then throw all that label into the new overlap
905: //
906: // Alternative: Figure out explicitly which segments were refined, and then
907: // communicated the refinement over the old overlap. Use this info to locally
908: // construct the new overlap and flip to get a decent mesh
909: sieveCompletion::scatterCones(refSieve, refSieve, reMesh->getDistSendOverlap(), refMesh->getDistRecvOverlap(), refMesh);
910: #endif
911: return refMesh;
912: };
913: };
914: };
915: #endif
916: #ifdef PETSC_HAVE_TETGEN
917: namespace TetGen {
918: template<typename Mesh>
919: class Generator {
920: public:
921: static Obj<Mesh> generateMesh(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false) {
922: typedef ALE::SieveAlg<Mesh> sieve_alg_type;
923: const int dim = 3;
924: Obj<Mesh> mesh = new Mesh(boundary->comm(), dim, boundary->debug());
925: const PetscMPIInt rank = mesh->commRank();
926: bool createConvexHull = false;
927: ::tetgenio in;
928: ::tetgenio out;
930: const Obj<typename Mesh::label_sequence>& vertices = boundary->depthStratum(0);
931: const Obj<typename Mesh::numbering_type>& vNumbering = boundary->getFactory()->getLocalNumbering(boundary, 0);
932: const Obj<typename Mesh::real_section_type>& coordinates = boundary->getRealSection("coordinates");
933: const Obj<typename Mesh::label_type>& markers = boundary->getLabel("marker");
936: in.numberofpoints = vertices->size();
937: if (in.numberofpoints > 0) {
938: in.pointlist = new double[in.numberofpoints*dim];
939: in.pointmarkerlist = new int[in.numberofpoints];
940: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
941: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
942: const int idx = vNumbering->getIndex(*v_iter);
944: for(int d = 0; d < dim; d++) {
945: in.pointlist[idx*dim + d] = array[d];
946: }
947: in.pointmarkerlist[idx] = boundary->getValue(markers, *v_iter);
948: }
949: }
950:
951: if (boundary->depth() != 0) { //our boundary mesh COULD be just a pointset; in which case depth = height = 0;
952: const Obj<typename Mesh::label_sequence>& facets = boundary->depthStratum(boundary->depth());
953: //PetscPrintf(boundary->comm(), "%d facets on the boundary\n", facets->size());
954: const Obj<typename Mesh::numbering_type>& fNumbering = boundary->getFactory()->getLocalNumbering(boundary, boundary->depth());
956: in.numberoffacets = facets->size();
957: if (in.numberoffacets > 0) {
958: in.facetlist = new tetgenio::facet[in.numberoffacets];
959: in.facetmarkerlist = new int[in.numberoffacets];
960: for(typename Mesh::label_sequence::iterator f_iter = facets->begin(); f_iter != facets->end(); ++f_iter) {
961: const Obj<typename sieve_alg_type::coneArray>& cone = sieve_alg_type::nCone(boundary, *f_iter, boundary->depth());
962: const int idx = fNumbering->getIndex(*f_iter);
963:
964: in.facetlist[idx].numberofpolygons = 1;
965: in.facetlist[idx].polygonlist = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
966: in.facetlist[idx].numberofholes = 0;
967: in.facetlist[idx].holelist = NULL;
968:
969: tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
970: int c = 0;
972: poly->numberofvertices = cone->size();
973: poly->vertexlist = new int[poly->numberofvertices];
974: for(typename sieve_alg_type::coneArray::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
975: const int vIdx = vNumbering->getIndex(*c_iter);
976:
977: poly->vertexlist[c++] = vIdx;
978: }
979: in.facetmarkerlist[idx] = boundary->getValue(markers, *f_iter);
980: }
981: }
982: }else {
983: createConvexHull = true;
984: }
986: in.numberofholes = 0;
987: if (rank == 0) {
988: // Normal operation
989: std::string args("pqezQ");
990: //constrained operation
991: if (constrained) {
992: args = "pezQ";
993: if (createConvexHull) {
994: args = "ezQ";
995: //PetscPrintf(boundary->comm(), "createConvexHull\n");
996: }
997: }
998: // Just make tetrahedrons
999: // std::string args("efzV");
1000: // Adds a center point
1001: // std::string args("pqezQi");
1002: // in.numberofaddpoints = 1;
1003: // in.addpointlist = new double[in.numberofaddpoints*dim];
1004: // in.addpointlist[0] = 0.5;
1005: // in.addpointlist[1] = 0.5;
1006: // in.addpointlist[2] = 0.5;
1008: //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.
1009: ::tetrahedralize((char *) args.c_str(), &in, &out);
1010: }
1011: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
1012: int numCorners = 4;
1013: int numCells = out.numberoftetrahedra;
1014: int *cells = out.tetrahedronlist;
1015: int numVertices = out.numberofpoints;
1016: double *coords = out.pointlist;
1018: if (!interpolate) {
1019: for(int c = 0; c < numCells; ++c) {
1020: int tmp = cells[c*4+0];
1021: cells[c*4+0] = cells[c*4+1];
1022: cells[c*4+1] = tmp;
1023: }
1024: }
1025: ALE::SieveBuilder<Mesh>::buildTopology(newSieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, mesh->getArrowSection("orientation"));
1026: mesh->setSieve(newSieve);
1027: mesh->stratify();
1028: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, dim, coords);
1029: const Obj<typename Mesh::label_type>& newMarkers = mesh->createLabel("marker");
1030:
1031: for(int v = 0; v < out.numberofpoints; v++) {
1032: if (out.pointmarkerlist[v]) {
1033: mesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
1034: }
1035: }
1036: if (interpolate) {
1037: if (out.edgemarkerlist) {
1038: for(int e = 0; e < out.numberofedges; e++) {
1039: if (out.edgemarkerlist[e]) {
1040: typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
1041: typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
1042: Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(endpointA, endpointB, 1);
1044: mesh->setValue(newMarkers, *edge->begin(), out.edgemarkerlist[e]);
1045: }
1046: }
1047: }
1048: if (out.trifacemarkerlist) {
1049: // Work around TetGen bug for raw tetrahedralization
1050: // The boundary faces are 0,1,4,5,8,9,11,12,13,15,16,17
1051: // for(int f = 0; f < out.numberoftrifaces; f++) {
1052: // if (out.trifacemarkerlist[f]) {
1053: // out.trifacemarkerlist[f] = 0;
1054: // } else {
1055: // out.trifacemarkerlist[f] = 1;
1056: // }
1057: // }
1058: for(int f = 0; f < out.numberoftrifaces; f++) {
1059: if (out.trifacemarkerlist[f]) {
1060: typename Mesh::point_type cornerA(out.trifacelist[f*3+0]+out.numberoftetrahedra);
1061: typename Mesh::point_type cornerB(out.trifacelist[f*3+1]+out.numberoftetrahedra);
1062: typename Mesh::point_type cornerC(out.trifacelist[f*3+2]+out.numberoftetrahedra);
1063: Obj<typename Mesh::sieve_type::supportSet> corners = typename Mesh::sieve_type::supportSet();
1064: Obj<typename Mesh::sieve_type::supportSet> edges = typename Mesh::sieve_type::supportSet();
1065: corners->clear();corners->insert(cornerA);corners->insert(cornerB);
1066: edges->insert(*newSieve->nJoin1(corners)->begin());
1067: corners->clear();corners->insert(cornerB);corners->insert(cornerC);
1068: edges->insert(*newSieve->nJoin1(corners)->begin());
1069: corners->clear();corners->insert(cornerC);corners->insert(cornerA);
1070: edges->insert(*newSieve->nJoin1(corners)->begin());
1071: const typename Mesh::point_type face = *newSieve->nJoin1(edges)->begin();
1072: const int faceMarker = out.trifacemarkerlist[f];
1073: const Obj<typename Mesh::coneArray> closure = sieve_alg_type::closure(mesh, face);
1074: const typename Mesh::coneArray::iterator end = closure->end();
1076: for(typename Mesh::coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
1077: mesh->setValue(newMarkers, *cl_iter, faceMarker);
1078: }
1079: }
1080: }
1081: }
1082: }
1083: return mesh;
1084: };
1087: static Obj<Mesh> generateMeshV(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false) {
1088: const int dim = 3;
1089: Obj<Mesh> mesh = new Mesh(boundary->comm(), dim, boundary->debug());
1090: const Obj<typename Mesh::sieve_type>& sieve = boundary->getSieve();
1091: const PetscMPIInt rank = mesh->commRank();
1092: bool createConvexHull = false;
1093: PetscErrorCode ierr;
1094: ::tetgenio in;
1095: ::tetgenio out;
1097: const Obj<typename Mesh::label_sequence>& vertices = boundary->depthStratum(0);
1098: const Obj<typename Mesh::label_type>& markers = boundary->getLabel("marker");
1099: const Obj<typename Mesh::real_section_type>& coordinates = boundary->getRealSection("coordinates");
1100: const Obj<typename Mesh::numbering_type>& vNumbering = boundary->getFactory()->getLocalNumbering(boundary, 0);
1102: in.numberofpoints = vertices->size();
1103: if (in.numberofpoints > 0) {
1104: const typename Mesh::label_sequence::iterator vEnd = vertices->end();
1106: in.pointlist = new double[in.numberofpoints*dim];
1107: in.pointmarkerlist = new int[in.numberofpoints];
1108: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
1109: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
1110: const int idx = vNumbering->getIndex(*v_iter);
1112: for(int d = 0; d < dim; ++d) {
1113: in.pointlist[idx*dim + d] = array[d];
1114: }
1115: in.pointmarkerlist[idx] = boundary->getValue(markers, *v_iter);
1116: }
1117: }
1118:
1119: // Our boundary mesh COULD be just a pointset; in which case depth = height = 0;
1120: if (boundary->depth() != 0) {
1121: const Obj<typename Mesh::label_sequence>& facets = boundary->depthStratum(boundary->depth());
1122: const Obj<typename Mesh::numbering_type>& fNumbering = boundary->getFactory()->getLocalNumbering(boundary, boundary->depth());
1124: in.numberoffacets = facets->size();
1125: if (in.numberoffacets > 0) {
1126: ALE::ISieveVisitor::NConeRetriever<typename Mesh::sieve_type> ncV(*sieve, 5);
1127: const typename Mesh::label_sequence::iterator fEnd = facets->end();
1129: in.facetlist = new tetgenio::facet[in.numberoffacets];
1130: in.facetmarkerlist = new int[in.numberoffacets];
1131: for(typename Mesh::label_sequence::iterator f_iter = facets->begin(); f_iter != fEnd; ++f_iter) {
1132: ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*sieve, *f_iter, ncV);
1133: const int idx = fNumbering->getIndex(*f_iter);
1134:
1135: in.facetlist[idx].numberofpolygons = 1;
1136: in.facetlist[idx].polygonlist = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
1137: in.facetlist[idx].numberofholes = 0;
1138: in.facetlist[idx].holelist = NULL;
1139:
1140: tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
1141: const size_t n = ncV.getSize();
1142: const typename Mesh::point_type *cone = ncV.getPoints();
1144: poly->numberofvertices = n;
1145: poly->vertexlist = new int[poly->numberofvertices];
1146: for(size_t c = 0; c < n; ++c) {
1147: poly->vertexlist[c] = vNumbering->getIndex(cone[c]);
1148: }
1149: in.facetmarkerlist[idx] = boundary->getValue(markers, *f_iter);
1150: ncV.clear();
1151: }
1152: }
1153: } else {
1154: createConvexHull = true;
1155: }
1156: const typename Mesh::holes_type& holes = mesh->getHoles();
1158: in.numberofholes = holes.size();
1159: if (in.numberofholes > 0) {
1160: PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);
1161: for(int h = 0; h < in.numberofholes; ++h) {
1162: for(int d = 0; d < dim; ++d) {
1163: in.holelist[h*dim+d] = holes[h][d];
1164: }
1165: }
1166: }
1167: if (rank == 0) {
1168: // Normal operation
1169: std::string args("pqezQ");
1170: // Constrained operation
1171: if (constrained) {
1172: args = "pezQ";
1173: if (createConvexHull) {
1174: args = "ezQ";
1175: }
1176: }
1177: // Just make tetrahedrons
1178: // std::string args("efzV");
1179: // Adds a center point
1180: // std::string args("pqezQi");
1181: // in.numberofaddpoints = 1;
1182: // in.addpointlist = new double[in.numberofaddpoints*dim];
1183: // in.addpointlist[0] = 0.5;
1184: // in.addpointlist[1] = 0.5;
1185: // in.addpointlist[2] = 0.5;
1186: ::tetrahedralize((char *) args.c_str(), &in, &out);
1187: }
1188: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
1189: const Obj<ALE::Mesh> m = new ALE::Mesh(mesh->comm(), dim, mesh->debug());
1190: const Obj<ALE::Mesh::sieve_type> newS = new ALE::Mesh::sieve_type(m->comm(), m->debug());
1191: int numCorners = 4;
1192: int numCells = out.numberoftetrahedra;
1193: int *cells = out.tetrahedronlist;
1194: int numVertices = out.numberofpoints;
1195: double *coords = out.pointlist;
1197: if (!interpolate) {
1198: // TetGen reports tetrahedra with the opposite orientation from what we expect
1199: for(int c = 0; c < numCells; ++c) {
1200: int tmp = cells[c*4+0];
1201: cells[c*4+0] = cells[c*4+1];
1202: cells[c*4+1] = tmp;
1203: }
1204: }
1205: ALE::SieveBuilder<ALE::Mesh>::buildTopology(newS, dim, numCells, cells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
1206: m->setSieve(newS);
1207: m->stratify();
1208: mesh->setSieve(newSieve);
1209: std::map<typename Mesh::point_type,typename Mesh::point_type> renumbering;
1210: ALE::ISieveConverter::convertSieve(*newS, *newSieve, renumbering, false);
1211: mesh->stratify();
1212: ALE::ISieveConverter::convertOrientation(*newS, *newSieve, renumbering, m->getArrowSection("orientation").ptr());
1213: ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, dim, coords);
1214: const Obj<typename Mesh::label_type>& newMarkers = mesh->createLabel("marker");
1215:
1216: for(int v = 0; v < out.numberofpoints; v++) {
1217: if (out.pointmarkerlist[v]) {
1218: mesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
1219: }
1220: }
1221: if (interpolate) {
1222: if (out.edgemarkerlist) {
1223: for(int e = 0; e < out.numberofedges; e++) {
1224: if (out.edgemarkerlist[e]) {
1225: typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
1226: typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
1227: Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(endpointA, endpointB, 1);
1229: mesh->setValue(newMarkers, *edge->begin(), out.edgemarkerlist[e]);
1230: }
1231: }
1232: }
1233: if (out.trifacemarkerlist) {
1234: // Work around TetGen bug for raw tetrahedralization
1235: // The boundary faces are 0,1,4,5,8,9,11,12,13,15,16,17
1236: // for(int f = 0; f < out.numberoftrifaces; f++) {
1237: // if (out.trifacemarkerlist[f]) {
1238: // out.trifacemarkerlist[f] = 0;
1239: // } else {
1240: // out.trifacemarkerlist[f] = 1;
1241: // }
1242: // }
1243: for(int f = 0; f < out.numberoftrifaces; f++) {
1244: if (out.trifacemarkerlist[f]) {
1245: typename Mesh::point_type cornerA(out.trifacelist[f*3+0]+out.numberoftetrahedra);
1246: typename Mesh::point_type cornerB(out.trifacelist[f*3+1]+out.numberoftetrahedra);
1247: typename Mesh::point_type cornerC(out.trifacelist[f*3+2]+out.numberoftetrahedra);
1248: Obj<typename Mesh::sieve_type::supportSet> corners = typename Mesh::sieve_type::supportSet();
1249: Obj<typename Mesh::sieve_type::supportSet> edges = typename Mesh::sieve_type::supportSet();
1250: corners->clear();corners->insert(cornerA);corners->insert(cornerB);
1251: edges->insert(*newS->nJoin1(corners)->begin());
1252: corners->clear();corners->insert(cornerB);corners->insert(cornerC);
1253: edges->insert(*newS->nJoin1(corners)->begin());
1254: corners->clear();corners->insert(cornerC);corners->insert(cornerA);
1255: edges->insert(*newS->nJoin1(corners)->begin());
1256: ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(30);
1257: const typename Mesh::point_type face = *newS->nJoin1(edges)->begin();
1258: const int faceMarker = out.trifacemarkerlist[f];
1260: ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*newSieve, face, pV);
1261: const size_t n = pV.getSize();
1262: const typename Mesh::point_type *cone = pV.getPoints();
1264: for(size_t c = 0; c < n; ++c) {
1265: mesh->setValue(newMarkers, cone[c], faceMarker);
1266: }
1267: pV.clear();
1268: }
1269: }
1270: }
1271: }
1272: mesh->copyHoles(boundary);
1273: return mesh;
1274: };
1275: };
1276: template<typename Mesh>
1277: class Refiner {
1278: public:
1279: static Obj<Mesh> refineMesh(const Obj<Mesh>& serialMesh, const double maxVolumes[], const bool interpolate = false) {
1280: typedef ALE::SieveAlg<Mesh> sieve_alg_type;
1281: const int dim = serialMesh->getDimension();
1282: const int depth = serialMesh->depth();
1283: const Obj<Mesh> refMesh = new Mesh(serialMesh->comm(), dim, serialMesh->debug());
1284: ::tetgenio in;
1285: ::tetgenio out;
1287: const Obj<typename Mesh::label_sequence>& vertices = serialMesh->depthStratum(0);
1288: const Obj<typename Mesh::label_type>& markers = serialMesh->getLabel("marker");
1289: const Obj<typename Mesh::real_section_type>& coordinates = serialMesh->getRealSection("coordinates");
1290: const Obj<typename Mesh::numbering_type>& vNumbering = serialMesh->getFactory()->getLocalNumbering(serialMesh, 0);
1292: in.numberofpoints = vertices->size();
1293: if (in.numberofpoints > 0) {
1294: in.pointlist = new double[in.numberofpoints*dim];
1295: in.pointmarkerlist = new int[in.numberofpoints];
1296: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
1297: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
1298: const int idx = vNumbering->getIndex(*v_iter);
1300: for(int d = 0; d < dim; d++) {
1301: in.pointlist[idx*dim + d] = array[d];
1302: }
1303: in.pointmarkerlist[idx] = serialMesh->getValue(markers, *v_iter);
1304: }
1305: }
1306: const Obj<typename Mesh::label_sequence>& cells = serialMesh->heightStratum(0);
1307: const Obj<typename Mesh::numbering_type>& cNumbering = serialMesh->getFactory()->getLocalNumbering(serialMesh, depth);
1309: in.numberofcorners = 4;
1310: in.numberoftetrahedra = cells->size();
1311: in.tetrahedronvolumelist = (double *) maxVolumes;
1312: if (in.numberoftetrahedra > 0) {
1313: in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
1314: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
1315: typedef ALE::SieveAlg<Mesh> sieve_alg_type;
1316: const Obj<typename sieve_alg_type::coneArray>& cone = sieve_alg_type::nCone(serialMesh, *c_iter, depth);
1317: const int idx = cNumbering->getIndex(*c_iter);
1318: int v = 0;
1320: for(typename Mesh::sieve_type::coneArray::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
1321: in.tetrahedronlist[idx*in.numberofcorners + v++] = vNumbering->getIndex(*v_iter);
1322: }
1323: }
1324: }
1325: if (serialMesh->depth() == 3) {
1326: const Obj<typename Mesh::label_sequence>& boundary = serialMesh->getLabelStratum("marker", 1);
1328: in.numberoftrifaces = 0;
1329: for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != boundary->end(); ++b_iter) {
1330: if (serialMesh->height(*b_iter) == 1) {
1331: in.numberoftrifaces++;
1332: }
1333: }
1334: if (in.numberoftrifaces > 0) {
1335: int f = 0;
1337: in.trifacelist = new int[in.numberoftrifaces*3];
1338: in.trifacemarkerlist = new int[in.numberoftrifaces];
1339: for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != boundary->end(); ++b_iter) {
1340: if (serialMesh->height(*b_iter) == 1) {
1341: const Obj<typename Mesh::coneArray>& cone = sieve_alg_type::nCone(serialMesh, *b_iter, 2);
1342: int p = 0;
1344: for(typename Mesh::coneArray::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
1345: in.trifacelist[f*3 + (p++)] = vNumbering->getIndex(*v_iter);
1346: }
1347: in.trifacemarkerlist[f++] = serialMesh->getValue(markers, *b_iter);
1348: }
1349: }
1350: }
1351: }
1353: in.numberofholes = 0;
1354: if (serialMesh->commRank() == 0) {
1355: std::string args("qezQra");
1357: ::tetrahedralize((char *) args.c_str(), &in, &out);
1358: }
1359: in.tetrahedronvolumelist = NULL;
1360: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(refMesh->comm(), refMesh->debug());
1361: int numCorners = 4;
1362: int numCells = out.numberoftetrahedra;
1363: int *newCells = out.tetrahedronlist;
1364: int numVertices = out.numberofpoints;
1365: double *coords = out.pointlist;
1367: if (!interpolate) {
1368: for(int c = 0; c < numCells; ++c) {
1369: int tmp = newCells[c*4+0];
1370: newCells[c*4+0] = newCells[c*4+1];
1371: newCells[c*4+1] = tmp;
1372: }
1373: }
1374: ALE::SieveBuilder<Mesh>::buildTopology(newSieve, dim, numCells, newCells, numVertices, interpolate, numCorners, -1, refMesh->getArrowSection("orientation"));
1375: refMesh->setSieve(newSieve);
1376: refMesh->stratify();
1377: ALE::SieveBuilder<Mesh>::buildCoordinates(refMesh, dim, coords);
1378: const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");
1381: for(int v = 0; v < out.numberofpoints; v++) {
1382: if (out.pointmarkerlist[v]) {
1383: refMesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
1384: }
1385: }
1386: if (interpolate) {
1387: if (out.edgemarkerlist) {
1388: for(int e = 0; e < out.numberofedges; e++) {
1389: if (out.edgemarkerlist[e]) {
1390: typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
1391: typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
1392: Obj<typename Mesh::sieve_type::supportSet> edge = newSieve->nJoin(endpointA, endpointB, 1);
1394: refMesh->setValue(newMarkers, *edge->begin(), out.edgemarkerlist[e]);
1395: }
1396: }
1397: }
1398: if (out.trifacemarkerlist) {
1399: for(int f = 0; f < out.numberoftrifaces; f++) {
1400: if (out.trifacemarkerlist[f]) {
1401: typename Mesh::point_type cornerA(out.trifacelist[f*3+0]+out.numberoftetrahedra);
1402: typename Mesh::point_type cornerB(out.trifacelist[f*3+1]+out.numberoftetrahedra);
1403: typename Mesh::point_type cornerC(out.trifacelist[f*3+2]+out.numberoftetrahedra);
1404: Obj<typename Mesh::sieve_type::supportSet> corners = typename Mesh::sieve_type::supportSet();
1405: Obj<typename Mesh::sieve_type::supportSet> edges = typename Mesh::sieve_type::supportSet();
1406: corners->clear();corners->insert(cornerA);corners->insert(cornerB);
1407: edges->insert(*newSieve->nJoin1(corners)->begin());
1408: corners->clear();corners->insert(cornerB);corners->insert(cornerC);
1409: edges->insert(*newSieve->nJoin1(corners)->begin());
1410: corners->clear();corners->insert(cornerC);corners->insert(cornerA);
1411: edges->insert(*newSieve->nJoin1(corners)->begin());
1412: const typename Mesh::point_type face = *newSieve->nJoin1(edges)->begin();
1413: const int faceMarker = out.trifacemarkerlist[f];
1414: const Obj<typename Mesh::coneArray> closure = sieve_alg_type::closure(refMesh, face);
1415: const typename Mesh::coneArray::iterator end = closure->end();
1417: for(typename Mesh::coneArray::iterator cl_iter = closure->begin(); cl_iter != end; ++cl_iter) {
1418: refMesh->setValue(newMarkers, *cl_iter, faceMarker);
1419: }
1420: }
1421: }
1422: }
1423: }
1424: if (refMesh->commSize() > 1) {
1425: return ALE::Distribution<Mesh>::distributeMesh(refMesh);
1426: }
1427: return refMesh;
1428: };
1429: static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false) {
1430: Obj<Mesh> serialMesh = ALE::Distribution<Mesh>::unifyMesh(mesh);
1431: const Obj<typename Mesh::real_section_type> serialMaxVolumes = ALE::Distribution<Mesh>::distributeSection(maxVolumes, serialMesh, serialMesh->getDistSendOverlap(), serialMesh->getDistRecvOverlap());
1433: return refineMesh(serialMesh, serialMaxVolumes->restrictSpace(), interpolate);
1434: };
1435: static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false) {
1436: Obj<Mesh> serialMesh;
1437: if (mesh->commSize() > 1) {
1438: serialMesh = ALE::Distribution<Mesh>::unifyMesh(mesh);
1439: } else {
1440: serialMesh = mesh;
1441: }
1442: const int numCells = serialMesh->heightStratum(0)->size();
1443: double *serialMaxVolumes = new double[numCells];
1445: for(int c = 0; c < numCells; c++) {
1446: serialMaxVolumes[c] = maxVolume;
1447: }
1448: const Obj<Mesh> refMesh = refineMesh(serialMesh, serialMaxVolumes, interpolate);
1449: delete [] serialMaxVolumes;
1450: return refMesh;
1451: };
1452: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolumes[], const bool interpolate = false, const bool forceSerial = false) {
1453: typedef typename Mesh::point_type point_type;
1454: const int dim = mesh->getDimension();
1455: const int depth = mesh->depth();
1456: const Obj<Mesh> refMesh = new Mesh(mesh->comm(), dim, mesh->debug());
1457: const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
1458: PetscErrorCode ierr;
1459: ::tetgenio in;
1460: ::tetgenio out;
1462: const Obj<typename Mesh::label_sequence>& vertices = mesh->depthStratum(0);
1463: const Obj<typename Mesh::label_type>& markers = mesh->getLabel("marker");
1464: const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
1465: const Obj<typename Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
1467: in.numberofpoints = vertices->size();
1468: if (in.numberofpoints > 0) {
1469: const typename Mesh::label_sequence::iterator vEnd = vertices->end();
1471: in.pointlist = new double[in.numberofpoints*dim];
1472: in.pointmarkerlist = new int[in.numberofpoints];
1473: for(typename Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
1474: const typename Mesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
1475: const int idx = vNumbering->getIndex(*v_iter);
1477: for(int d = 0; d < dim; d++) {
1478: in.pointlist[idx*dim + d] = array[d];
1479: }
1480: in.pointmarkerlist[idx] = mesh->getValue(markers, *v_iter);
1481: }
1482: }
1483: const Obj<typename Mesh::label_sequence>& cells = mesh->heightStratum(0);
1484: const Obj<typename Mesh::numbering_type>& cNumbering = mesh->getFactory()->getLocalNumbering(mesh, depth);
1486: in.numberofcorners = 4;
1487: in.numberoftetrahedra = cells->size();
1488: in.tetrahedronvolumelist = (double *) maxVolumes;
1489: if (in.numberoftetrahedra > 0) {
1490: in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
1491: if (mesh->depth() == 1) {
1492: ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(4);
1493: const typename Mesh::label_sequence::iterator cEnd = cells->end();
1495: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
1496: sieve->cone(*c_iter, pV);
1497: const int idx = cNumbering->getIndex(*c_iter);
1498: const size_t n = pV.getSize();
1499: const point_type *cone = pV.getPoints();
1501: assert(n == 4);
1502: for(int v = 0; v < 4; ++v) {
1503: in.tetrahedronlist[idx*in.numberofcorners + v] = vNumbering->getIndex(cone[v]);
1504: }
1505: pV.clear();
1506: }
1507: } else if (mesh->depth() == 3) {
1508: // Need extra space due to early error checking
1509: ALE::ISieveVisitor::NConeRetriever<typename Mesh::sieve_type> ncV(*sieve, 5);
1510: const typename Mesh::label_sequence::iterator cEnd = cells->end();
1512: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
1513: ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
1514: const int idx = cNumbering->getIndex(*c_iter);
1515: const size_t n = ncV.getSize();
1516: const point_type *cone = ncV.getPoints();
1518: assert(n == 4);
1519: for(int v = 0; v < 4; ++v) {
1520: in.tetrahedronlist[idx*in.numberofcorners + v] = vNumbering->getIndex(cone[v]);
1521: }
1522: ncV.clear();
1523: }
1524: } else {
1525: throw ALE::Exception("Invalid sieve: Cannot gives sieves of arbitrary depth to TetGen");
1526: }
1527: }
1528: if (depth == 3) {
1529: const Obj<typename Mesh::label_sequence>& boundary = mesh->getLabelStratum("marker", 1);
1530: const typename Mesh::label_sequence::iterator bEnd = boundary->end();
1532: in.numberoftrifaces = 0;
1533: for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != bEnd; ++b_iter) {
1534: if (mesh->height(*b_iter) == 1) {
1535: in.numberoftrifaces++;
1536: }
1537: }
1538: if (in.numberoftrifaces > 0) {
1539: ALE::ISieveVisitor::NConeRetriever<typename Mesh::sieve_type> ncV(*sieve, 5);
1540: int f = 0;
1542: in.trifacelist = new int[in.numberoftrifaces*3];
1543: in.trifacemarkerlist = new int[in.numberoftrifaces];
1544: for(typename Mesh::label_sequence::iterator b_iter = boundary->begin(); b_iter != bEnd; ++b_iter) {
1545: if (mesh->height(*b_iter) == 1) {
1546: ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*sieve, *b_iter, ncV);
1547: const size_t n = ncV.getSize();
1548: const point_type *cone = ncV.getPoints();
1550: for(size_t p = 0; p < n; ++p) {
1551: in.trifacelist[f*3 + p] = vNumbering->getIndex(cone[p]);
1552: }
1553: in.trifacemarkerlist[f++] = mesh->getValue(markers, *b_iter);
1554: ncV.clear();
1555: }
1556: }
1557: }
1558: }
1559: const typename Mesh::holes_type& holes = mesh->getHoles();
1561: in.numberofholes = holes.size();
1562: if (in.numberofholes > 0) {
1563: PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);
1564: for(int h = 0; h < in.numberofholes; ++h) {
1565: for(int d = 0; d < dim; ++d) {
1566: in.holelist[h*dim+d] = holes[h][d];
1567: }
1568: }
1569: }
1570: if (mesh->commRank() == 0) {
1571: std::string args("qezQra");
1573: ::tetrahedralize((char *) args.c_str(), &in, &out);
1574: }
1575: in.tetrahedronvolumelist = NULL;
1576: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(refMesh->comm(), refMesh->debug());
1577: const Obj<ALE::Mesh> m = new ALE::Mesh(mesh->comm(), dim, mesh->debug());
1578: const Obj<ALE::Mesh::sieve_type> newS = new ALE::Mesh::sieve_type(m->comm(), m->debug());
1579: int numCorners = 4;
1580: int numCells = out.numberoftetrahedra;
1581: int *newCells = out.tetrahedronlist;
1582: int numVertices = out.numberofpoints;
1583: double *coords = out.pointlist;
1585: if (!interpolate) {
1586: for(int c = 0; c < numCells; ++c) {
1587: int tmp = newCells[c*4+0];
1588: newCells[c*4+0] = newCells[c*4+1];
1589: newCells[c*4+1] = tmp;
1590: }
1591: }
1592: ALE::SieveBuilder<ALE::Mesh>::buildTopology(newS, dim, numCells, newCells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
1593: m->setSieve(newS);
1594: m->stratify();
1595: refMesh->setSieve(newSieve);
1596: std::map<typename Mesh::point_type,typename Mesh::point_type> renumbering;
1597: ALE::ISieveConverter::convertSieve(*newS, *newSieve, renumbering, false);
1598: refMesh->stratify();
1599: ALE::ISieveConverter::convertOrientation(*newS, *newSieve, renumbering, m->getArrowSection("orientation").ptr());
1600: ALE::SieveBuilder<Mesh>::buildCoordinates(refMesh, dim, coords);
1601: const Obj<typename Mesh::label_type>& newMarkers = refMesh->createLabel("marker");
1603: for(int v = 0; v < out.numberofpoints; v++) {
1604: if (out.pointmarkerlist[v]) {
1605: refMesh->setValue(newMarkers, v+out.numberoftetrahedra, out.pointmarkerlist[v]);
1606: }
1607: }
1608: if (interpolate) {
1609: if (out.edgemarkerlist) {
1610: for(int e = 0; e < out.numberofedges; e++) {
1611: if (out.edgemarkerlist[e]) {
1612: typename Mesh::point_type endpointA(out.edgelist[e*2+0]+out.numberoftetrahedra);
1613: typename Mesh::point_type endpointB(out.edgelist[e*2+1]+out.numberoftetrahedra);
1614: Obj<typename Mesh::sieve_type::supportSet> edge = newS->nJoin(endpointA, endpointB, 1);
1616: refMesh->setValue(newMarkers, *edge->begin(), out.edgemarkerlist[e]);
1617: }
1618: }
1619: }
1620: if (out.trifacemarkerlist) {
1621: for(int f = 0; f < out.numberoftrifaces; f++) {
1622: if (out.trifacemarkerlist[f]) {
1623: typename Mesh::point_type cornerA(out.trifacelist[f*3+0]+out.numberoftetrahedra);
1624: typename Mesh::point_type cornerB(out.trifacelist[f*3+1]+out.numberoftetrahedra);
1625: typename Mesh::point_type cornerC(out.trifacelist[f*3+2]+out.numberoftetrahedra);
1626: Obj<typename Mesh::sieve_type::supportSet> corners = typename Mesh::sieve_type::supportSet();
1627: Obj<typename Mesh::sieve_type::supportSet> edges = typename Mesh::sieve_type::supportSet();
1628: corners->clear();corners->insert(cornerA);corners->insert(cornerB);
1629: edges->insert(*newS->nJoin1(corners)->begin());
1630: corners->clear();corners->insert(cornerB);corners->insert(cornerC);
1631: edges->insert(*newS->nJoin1(corners)->begin());
1632: corners->clear();corners->insert(cornerC);corners->insert(cornerA);
1633: edges->insert(*newS->nJoin1(corners)->begin());
1634: ALE::ISieveVisitor::PointRetriever<typename Mesh::sieve_type> pV(30);
1635: const typename Mesh::point_type face = *newS->nJoin1(edges)->begin();
1636: const int faceMarker = out.trifacemarkerlist[f];
1638: ALE::ISieveTraversal<typename Mesh::sieve_type>::orientedClosure(*newSieve, face, pV);
1639: const size_t n = pV.getSize();
1640: const typename Mesh::point_type *cone = pV.getPoints();
1642: for(size_t c = 0; c < n; ++c) {
1643: refMesh->setValue(newMarkers, cone[c], faceMarker);
1644: }
1645: pV.clear();
1646: }
1647: }
1648: }
1649: }
1650: #if 0
1651: if (refMesh->commSize() > 1) {
1652: return ALE::Distribution<Mesh>::distributeMesh(refMesh);
1653: }
1654: #endif
1655: return refMesh;
1656: };
1657: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false) {
1658: throw ALE::Exception("Not yet implemented");
1659: };
1660: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false) {
1661: const int numCells = mesh->heightStratum(0)->size();
1662: double *serialMaxVolumes = new double[numCells];
1664: for(int c = 0; c < numCells; c++) {
1665: serialMaxVolumes[c] = maxVolume;
1666: }
1667: const Obj<Mesh> refMesh = refineMeshV(mesh, serialMaxVolumes, interpolate);
1668: delete [] serialMaxVolumes;
1669: return refMesh;
1670: };
1671: };
1672: };
1673: #endif
1674: template<typename Mesh>
1675: class Generator {
1676: public:
1677: static Obj<Mesh> generateMesh(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false) {
1678: int dim = boundary->getDimension();
1680: if (dim == 1) {
1681: #ifdef PETSC_HAVE_TRIANGLE
1682: return ALE::Triangle::Generator<Mesh>::generateMesh(boundary, interpolate, constrained);
1683: #else
1684: throw ALE::Exception("Mesh generation currently requires Triangle to be installed. Use --download-triangle during configure.");
1685: #endif
1686: } else if (dim == 2) {
1687: #ifdef PETSC_HAVE_TETGEN
1688: return ALE::TetGen::Generator<Mesh>::generateMesh(boundary, interpolate, constrained);
1689: #else
1690: throw ALE::Exception("Mesh generation currently requires TetGen to be installed. Use --download-tetgen during configure.");
1691: #endif
1692: }
1693: return NULL;
1694: };
1695: static Obj<Mesh> generateMeshV(const Obj<Mesh>& boundary, const bool interpolate = false, const bool constrained = false) {
1696: int dim = boundary->getDimension();
1698: if (dim == 1) {
1699: #ifdef PETSC_HAVE_TRIANGLE
1700: return ALE::Triangle::Generator<Mesh>::generateMeshV(boundary, interpolate, constrained);
1701: #else
1702: throw ALE::Exception("Mesh generation currently requires Triangle to be installed. Use --download-triangle during configure.");
1703: #endif
1704: } else if (dim == 2) {
1705: #ifdef PETSC_HAVE_TETGEN
1706: return ALE::TetGen::Generator<Mesh>::generateMeshV(boundary, interpolate, constrained);
1707: #else
1708: throw ALE::Exception("Mesh generation currently requires TetGen to be installed. Use --download-tetgen during configure.");
1709: #endif
1710: }
1711: return NULL;
1712: };
1713: static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false) {
1714: int dim = mesh->getDimension();
1716: if (dim == 2) {
1717: #ifdef PETSC_HAVE_TRIANGLE
1718: return ALE::Triangle::Refiner<Mesh>::refineMesh(mesh, maxVolumes, interpolate);
1719: #else
1720: throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
1721: #endif
1722: } else if (dim == 3) {
1723: #ifdef PETSC_HAVE_TETGEN
1724: return ALE::TetGen::Refiner<Mesh>::refineMesh(mesh, maxVolumes, interpolate);
1725: #else
1726: throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
1727: #endif
1728: }
1729: return NULL;
1730: };
1731: static Obj<Mesh> refineMesh(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false, const bool forceSerial = false) {
1732: int dim = mesh->getDimension();
1734: if (dim == 2) {
1735: #ifdef PETSC_HAVE_TRIANGLE
1736: return ALE::Triangle::Refiner<Mesh>::refineMesh(mesh, maxVolume, interpolate, forceSerial);
1737: #else
1738: throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
1739: #endif
1740: } else if (dim == 3) {
1741: #ifdef PETSC_HAVE_TETGEN
1742: return ALE::TetGen::Refiner<Mesh>::refineMesh(mesh, maxVolume, interpolate);
1743: #else
1744: throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
1745: #endif
1746: }
1747: return NULL;
1748: };
1749: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const Obj<typename Mesh::real_section_type>& maxVolumes, const bool interpolate = false) {
1750: int dim = mesh->getDimension();
1752: if (dim == 2) {
1753: #ifdef PETSC_HAVE_TRIANGLE
1754: return ALE::Triangle::Refiner<Mesh>::refineMeshV(mesh, maxVolumes, interpolate);
1755: #else
1756: throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
1757: #endif
1758: } else if (dim == 3) {
1759: #ifdef PETSC_HAVE_TETGEN
1760: return ALE::TetGen::Refiner<Mesh>::refineMeshV(mesh, maxVolumes, interpolate);
1761: #else
1762: throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
1763: #endif
1764: }
1765: return NULL;
1766: };
1767: static Obj<Mesh> refineMeshV(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false, const bool forceSerial = false) {
1768: int dim = mesh->getDimension();
1770: if (dim == 2) {
1771: #ifdef PETSC_HAVE_TRIANGLE
1772: return ALE::Triangle::Refiner<Mesh>::refineMeshV(mesh, maxVolume, interpolate, forceSerial);
1773: #else
1774: throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
1775: #endif
1776: } else if (dim == 3) {
1777: #ifdef PETSC_HAVE_TETGEN
1778: return ALE::TetGen::Refiner<Mesh>::refineMeshV(mesh, maxVolume, interpolate);
1779: #else
1780: throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
1781: #endif
1782: }
1783: return NULL;
1784: };
1785: static Obj<Mesh> refineMeshLocal(const Obj<Mesh>& mesh, const double maxVolume, const bool interpolate = false) {
1786: int dim = mesh->getDimension();
1788: if (dim == 2) {
1789: #ifdef PETSC_HAVE_TRIANGLE
1790: return ALE::Triangle::Refiner<Mesh>::refineMeshLocal(mesh, maxVolume, interpolate);
1791: #else
1792: throw ALE::Exception("Mesh refinement currently requires Triangle to be installed. Use --download-triangle during configure.");
1793: #endif
1794: } else if (dim == 3) {
1795: #ifdef PETSC_HAVE_TETGEN
1796: return ALE::TetGen::Refiner<Mesh>::refineMesh(mesh, maxVolume, interpolate);
1797: #else
1798: throw ALE::Exception("Mesh refinement currently requires TetGen to be installed. Use --download-tetgen during configure.");
1799: #endif
1800: }
1801: return NULL;
1802: };
1803: };
1804: }
1806: #endif