Actual source code: plexgenerate.c
petsc-3.8.4 2018-03-24
1: #include <petsc/private/dmpleximpl.h>
3: PetscErrorCode DMPlexInvertCell_Internal(PetscInt dim, PetscInt numCorners, PetscInt cone[])
4: {
5: int tmpc;
8: if (dim != 3) return(0);
9: switch (numCorners) {
10: case 4:
11: tmpc = cone[0];
12: cone[0] = cone[1];
13: cone[1] = tmpc;
14: break;
15: case 8:
16: tmpc = cone[1];
17: cone[1] = cone[3];
18: cone[3] = tmpc;
19: break;
20: default: break;
21: }
22: return(0);
23: }
25: /*@C
26: DMPlexInvertCell - This flips tetrahedron and hexahedron orientation since Plex stores them internally with outward normals. Other cells are left untouched.
28: Input Parameters:
29: + numCorners - The number of vertices in a cell
30: - cone - The incoming cone
32: Output Parameter:
33: . cone - The inverted cone (in-place)
35: Level: developer
37: .seealso: DMPlexGenerate()
38: @*/
39: PetscErrorCode DMPlexInvertCell(PetscInt dim, PetscInt numCorners, int cone[])
40: {
41: int tmpc;
44: if (dim != 3) return(0);
45: switch (numCorners) {
46: case 4:
47: tmpc = cone[0];
48: cone[0] = cone[1];
49: cone[1] = tmpc;
50: break;
51: case 8:
52: tmpc = cone[1];
53: cone[1] = cone[3];
54: cone[3] = tmpc;
55: break;
56: default: break;
57: }
58: return(0);
59: }
61: /* This is to fix the tetrahedron orientation from TetGen */
62: PETSC_UNUSED static PetscErrorCode DMPlexInvertCells_Internal(PetscInt dim, PetscInt numCells, PetscInt numCorners, int cells[])
63: {
64: PetscInt bound = numCells*numCorners, coff;
68: for (coff = 0; coff < bound; coff += numCorners) {
69: DMPlexInvertCell(dim, numCorners, &cells[coff]);
70: }
71: return(0);
72: }
74: /*@C
75: DMPlexTriangleSetOptions - Set the options used for the Triangle mesh generator
77: Not Collective
79: Inputs Parameters:
80: + dm - The DMPlex object
81: - opts - The command line options
83: Level: developer
85: .keywords: mesh, points
86: .seealso: DMPlexTetgenSetOptions(), DMPlexGenerate()
87: @*/
88: PetscErrorCode DMPlexTriangleSetOptions(DM dm, const char *opts)
89: {
90: DM_Plex *mesh = (DM_Plex*) dm->data;
96: PetscFree(mesh->triangleOpts);
97: PetscStrallocpy(opts, &mesh->triangleOpts);
98: return(0);
99: }
101: /*@C
102: DMPlexTetgenSetOptions - Set the options used for the Tetgen mesh generator
104: Not Collective
106: Inputs Parameters:
107: + dm - The DMPlex object
108: - opts - The command line options
110: Level: developer
112: .keywords: mesh, points
113: .seealso: DMPlexTriangleSetOptions(), DMPlexGenerate()
114: @*/
115: PetscErrorCode DMPlexTetgenSetOptions(DM dm, const char *opts)
116: {
117: DM_Plex *mesh = (DM_Plex*) dm->data;
123: PetscFree(mesh->tetgenOpts);
124: PetscStrallocpy(opts, &mesh->tetgenOpts);
125: return(0);
126: }
128: #if defined(PETSC_HAVE_TRIANGLE)
129: #include <triangle.h>
131: static PetscErrorCode InitInput_Triangle(struct triangulateio *inputCtx)
132: {
134: inputCtx->numberofpoints = 0;
135: inputCtx->numberofpointattributes = 0;
136: inputCtx->pointlist = NULL;
137: inputCtx->pointattributelist = NULL;
138: inputCtx->pointmarkerlist = NULL;
139: inputCtx->numberofsegments = 0;
140: inputCtx->segmentlist = NULL;
141: inputCtx->segmentmarkerlist = NULL;
142: inputCtx->numberoftriangleattributes = 0;
143: inputCtx->trianglelist = NULL;
144: inputCtx->numberofholes = 0;
145: inputCtx->holelist = NULL;
146: inputCtx->numberofregions = 0;
147: inputCtx->regionlist = NULL;
148: return(0);
149: }
151: static PetscErrorCode InitOutput_Triangle(struct triangulateio *outputCtx)
152: {
154: outputCtx->numberofpoints = 0;
155: outputCtx->pointlist = NULL;
156: outputCtx->pointattributelist = NULL;
157: outputCtx->pointmarkerlist = NULL;
158: outputCtx->numberoftriangles = 0;
159: outputCtx->trianglelist = NULL;
160: outputCtx->triangleattributelist = NULL;
161: outputCtx->neighborlist = NULL;
162: outputCtx->segmentlist = NULL;
163: outputCtx->segmentmarkerlist = NULL;
164: outputCtx->numberofedges = 0;
165: outputCtx->edgelist = NULL;
166: outputCtx->edgemarkerlist = NULL;
167: return(0);
168: }
170: static PetscErrorCode FiniOutput_Triangle(struct triangulateio *outputCtx)
171: {
173: free(outputCtx->pointlist);
174: free(outputCtx->pointmarkerlist);
175: free(outputCtx->segmentlist);
176: free(outputCtx->segmentmarkerlist);
177: free(outputCtx->edgelist);
178: free(outputCtx->edgemarkerlist);
179: free(outputCtx->trianglelist);
180: free(outputCtx->neighborlist);
181: return(0);
182: }
184: PetscErrorCode DMPlexGenerate_Triangle(DM boundary, PetscBool interpolate, DM *dm)
185: {
186: MPI_Comm comm;
187: DM_Plex *mesh = (DM_Plex *) boundary->data;
188: PetscInt dim = 2;
189: const PetscBool createConvexHull = PETSC_FALSE;
190: const PetscBool constrained = PETSC_FALSE;
191: const char *labelName = "marker";
192: const char *labelName2 = "Face Sets";
193: struct triangulateio in;
194: struct triangulateio out;
195: DMLabel label, label2;
196: PetscInt vStart, vEnd, v, eStart, eEnd, e;
197: PetscMPIInt rank;
198: PetscErrorCode ierr;
201: PetscObjectGetComm((PetscObject)boundary,&comm);
202: MPI_Comm_rank(comm, &rank);
203: InitInput_Triangle(&in);
204: InitOutput_Triangle(&out);
205: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
206: DMGetLabel(boundary, labelName, &label);
207: DMGetLabel(boundary, labelName2, &label2);
209: in.numberofpoints = vEnd - vStart;
210: if (in.numberofpoints > 0) {
211: PetscSection coordSection;
212: Vec coordinates;
213: PetscScalar *array;
215: PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
216: PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
217: DMGetCoordinatesLocal(boundary, &coordinates);
218: DMGetCoordinateSection(boundary, &coordSection);
219: VecGetArray(coordinates, &array);
220: for (v = vStart; v < vEnd; ++v) {
221: const PetscInt idx = v - vStart;
222: PetscInt off, d;
224: PetscSectionGetOffset(coordSection, v, &off);
225: for (d = 0; d < dim; ++d) {
226: in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
227: }
228: if (label) {DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);}
229: }
230: VecRestoreArray(coordinates, &array);
231: }
232: DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);
233: in.numberofsegments = eEnd - eStart;
234: if (in.numberofsegments > 0) {
235: PetscMalloc1(in.numberofsegments*2, &in.segmentlist);
236: PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist);
237: for (e = eStart; e < eEnd; ++e) {
238: const PetscInt idx = e - eStart;
239: const PetscInt *cone;
241: DMPlexGetCone(boundary, e, &cone);
243: in.segmentlist[idx*2+0] = cone[0] - vStart;
244: in.segmentlist[idx*2+1] = cone[1] - vStart;
246: if (label) {DMLabelGetValue(label, e, &in.segmentmarkerlist[idx]);}
247: }
248: }
249: #if 0 /* Do not currently support holes */
250: PetscReal *holeCoords;
251: PetscInt h, d;
253: DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
254: if (in.numberofholes > 0) {
255: PetscMalloc1(in.numberofholes*dim, &in.holelist);
256: for (h = 0; h < in.numberofholes; ++h) {
257: for (d = 0; d < dim; ++d) {
258: in.holelist[h*dim+d] = holeCoords[h*dim+d];
259: }
260: }
261: }
262: #endif
263: if (!rank) {
264: char args[32];
266: /* Take away 'Q' for verbose output */
267: PetscStrcpy(args, "pqezQ");
268: if (createConvexHull) {PetscStrcat(args, "c");}
269: if (constrained) {PetscStrcpy(args, "zepDQ");}
270: if (mesh->triangleOpts) {triangulate(mesh->triangleOpts, &in, &out, NULL);}
271: else {triangulate(args, &in, &out, NULL);}
272: }
273: PetscFree(in.pointlist);
274: PetscFree(in.pointmarkerlist);
275: PetscFree(in.segmentlist);
276: PetscFree(in.segmentmarkerlist);
277: PetscFree(in.holelist);
279: {
280: DMLabel glabel = NULL;
281: DMLabel glabel2 = NULL;
282: const PetscInt numCorners = 3;
283: const PetscInt numCells = out.numberoftriangles;
284: const PetscInt numVertices = out.numberofpoints;
285: const int *cells = out.trianglelist;
286: const double *meshCoords = out.pointlist;
288: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
289: if (label) {DMCreateLabel(*dm, labelName); DMGetLabel(*dm, labelName, &glabel);}
290: if (label2) {DMCreateLabel(*dm, labelName2); DMGetLabel(*dm, labelName2, &glabel2);}
291: /* Set labels */
292: for (v = 0; v < numVertices; ++v) {
293: if (out.pointmarkerlist[v]) {
294: if (glabel) {DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);}
295: }
296: }
297: if (interpolate) {
298: for (e = 0; e < out.numberofedges; e++) {
299: if (out.edgemarkerlist[e]) {
300: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
301: const PetscInt *edges;
302: PetscInt numEdges;
304: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
305: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
306: if (glabel) {DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);}
307: if (glabel2) {DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]);}
308: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
309: }
310: }
311: }
312: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
313: }
314: #if 0 /* Do not currently support holes */
315: DMPlexCopyHoles(*dm, boundary);
316: #endif
317: FiniOutput_Triangle(&out);
318: return(0);
319: }
321: PetscErrorCode DMPlexRefine_Triangle(DM dm, double *maxVolumes, DM *dmRefined)
322: {
323: MPI_Comm comm;
324: PetscInt dim = 2;
325: const char *labelName = "marker";
326: struct triangulateio in;
327: struct triangulateio out;
328: DMLabel label;
329: PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
330: PetscMPIInt rank;
331: PetscErrorCode ierr;
334: PetscObjectGetComm((PetscObject)dm,&comm);
335: MPI_Comm_rank(comm, &rank);
336: InitInput_Triangle(&in);
337: InitOutput_Triangle(&out);
338: DMPlexGetDepth(dm, &depth);
339: MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
340: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
341: DMGetLabel(dm, labelName, &label);
343: in.numberofpoints = vEnd - vStart;
344: if (in.numberofpoints > 0) {
345: PetscSection coordSection;
346: Vec coordinates;
347: PetscScalar *array;
349: PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
350: PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
351: DMGetCoordinatesLocal(dm, &coordinates);
352: DMGetCoordinateSection(dm, &coordSection);
353: VecGetArray(coordinates, &array);
354: for (v = vStart; v < vEnd; ++v) {
355: const PetscInt idx = v - vStart;
356: PetscInt off, d;
358: PetscSectionGetOffset(coordSection, v, &off);
359: for (d = 0; d < dim; ++d) {
360: in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
361: }
362: if (label) {DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);}
363: }
364: VecRestoreArray(coordinates, &array);
365: }
366: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
368: in.numberofcorners = 3;
369: in.numberoftriangles = cEnd - cStart;
371: in.trianglearealist = (double*) maxVolumes;
372: if (in.numberoftriangles > 0) {
373: PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);
374: for (c = cStart; c < cEnd; ++c) {
375: const PetscInt idx = c - cStart;
376: PetscInt *closure = NULL;
377: PetscInt closureSize;
379: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
380: if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize);
381: for (v = 0; v < 3; ++v) {
382: in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart;
383: }
384: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
385: }
386: }
387: /* TODO: Segment markers are missing on input */
388: #if 0 /* Do not currently support holes */
389: PetscReal *holeCoords;
390: PetscInt h, d;
392: DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
393: if (in.numberofholes > 0) {
394: PetscMalloc1(in.numberofholes*dim, &in.holelist);
395: for (h = 0; h < in.numberofholes; ++h) {
396: for (d = 0; d < dim; ++d) {
397: in.holelist[h*dim+d] = holeCoords[h*dim+d];
398: }
399: }
400: }
401: #endif
402: if (!rank) {
403: char args[32];
405: /* Take away 'Q' for verbose output */
406: PetscStrcpy(args, "pqezQra");
407: triangulate(args, &in, &out, NULL);
408: }
409: PetscFree(in.pointlist);
410: PetscFree(in.pointmarkerlist);
411: PetscFree(in.segmentlist);
412: PetscFree(in.segmentmarkerlist);
413: PetscFree(in.trianglelist);
415: {
416: DMLabel rlabel = NULL;
417: const PetscInt numCorners = 3;
418: const PetscInt numCells = out.numberoftriangles;
419: const PetscInt numVertices = out.numberofpoints;
420: const int *cells = out.trianglelist;
421: const double *meshCoords = out.pointlist;
422: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
424: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
425: if (label) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*dmRefined, labelName, &rlabel);}
426: /* Set labels */
427: for (v = 0; v < numVertices; ++v) {
428: if (out.pointmarkerlist[v]) {
429: if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);}
430: }
431: }
432: if (interpolate) {
433: PetscInt e;
435: for (e = 0; e < out.numberofedges; e++) {
436: if (out.edgemarkerlist[e]) {
437: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
438: const PetscInt *edges;
439: PetscInt numEdges;
441: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
442: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
443: if (rlabel) {DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);}
444: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
445: }
446: }
447: }
448: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
449: }
450: #if 0 /* Do not currently support holes */
451: DMPlexCopyHoles(*dm, boundary);
452: #endif
453: FiniOutput_Triangle(&out);
454: return(0);
455: }
456: #endif /* PETSC_HAVE_TRIANGLE */
458: #if defined(PETSC_HAVE_TETGEN)
459: #include <tetgen.h>
460: PetscErrorCode DMPlexGenerate_Tetgen(DM boundary, PetscBool interpolate, DM *dm)
461: {
462: MPI_Comm comm;
463: DM_Plex *mesh = (DM_Plex *) boundary->data;
464: const PetscInt dim = 3;
465: const char *labelName = "marker";
466: ::tetgenio in;
467: ::tetgenio out;
468: DMLabel label;
469: PetscInt vStart, vEnd, v, fStart, fEnd, f;
470: PetscMPIInt rank;
474: PetscObjectGetComm((PetscObject)boundary,&comm);
475: MPI_Comm_rank(comm, &rank);
476: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
477: DMGetLabel(boundary, labelName, &label);
479: in.numberofpoints = vEnd - vStart;
480: if (in.numberofpoints > 0) {
481: PetscSection coordSection;
482: Vec coordinates;
483: PetscScalar *array;
485: in.pointlist = new double[in.numberofpoints*dim];
486: in.pointmarkerlist = new int[in.numberofpoints];
488: DMGetCoordinatesLocal(boundary, &coordinates);
489: DMGetCoordinateSection(boundary, &coordSection);
490: VecGetArray(coordinates, &array);
491: for (v = vStart; v < vEnd; ++v) {
492: const PetscInt idx = v - vStart;
493: PetscInt off, d;
495: PetscSectionGetOffset(coordSection, v, &off);
496: for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
497: if (label) {
498: PetscInt val;
500: DMLabelGetValue(label, v, &val);
501: in.pointmarkerlist[idx] = (int) val;
502: }
503: }
504: VecRestoreArray(coordinates, &array);
505: }
506: DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
508: in.numberoffacets = fEnd - fStart;
509: if (in.numberoffacets > 0) {
510: in.facetlist = new tetgenio::facet[in.numberoffacets];
511: in.facetmarkerlist = new int[in.numberoffacets];
512: for (f = fStart; f < fEnd; ++f) {
513: const PetscInt idx = f - fStart;
514: PetscInt *points = NULL, numPoints, p, numVertices = 0, v;
516: in.facetlist[idx].numberofpolygons = 1;
517: in.facetlist[idx].polygonlist = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
518: in.facetlist[idx].numberofholes = 0;
519: in.facetlist[idx].holelist = NULL;
521: DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
522: for (p = 0; p < numPoints*2; p += 2) {
523: const PetscInt point = points[p];
524: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
525: }
527: tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
528: poly->numberofvertices = numVertices;
529: poly->vertexlist = new int[poly->numberofvertices];
530: for (v = 0; v < numVertices; ++v) {
531: const PetscInt vIdx = points[v] - vStart;
532: poly->vertexlist[v] = vIdx;
533: }
534: if (label) {
535: PetscInt val;
537: DMLabelGetValue(label, f, &val);
538: in.facetmarkerlist[idx] = (int) val;
539: }
540: DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
541: }
542: }
543: if (!rank) {
544: char args[32];
546: /* Take away 'Q' for verbose output */
547: PetscStrcpy(args, "pqezQ");
548: if (mesh->tetgenOpts) {::tetrahedralize(mesh->tetgenOpts, &in, &out);}
549: else {::tetrahedralize(args, &in, &out);}
550: }
551: {
552: DMLabel glabel = NULL;
553: const PetscInt numCorners = 4;
554: const PetscInt numCells = out.numberoftetrahedra;
555: const PetscInt numVertices = out.numberofpoints;
556: const double *meshCoords = out.pointlist;
557: int *cells = out.tetrahedronlist;
559: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
560: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
561: if (label) {DMCreateLabel(*dm, labelName); DMGetLabel(*dm, labelName, &glabel);}
562: /* Set labels */
563: for (v = 0; v < numVertices; ++v) {
564: if (out.pointmarkerlist[v]) {
565: if (glabel) {DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);}
566: }
567: }
568: if (interpolate) {
569: #if 0
570: PetscInt e;
572: /* This check is never actually executed for ctetgen (which never returns edgemarkers) and seems to be broken for
573: * tetgen */
574: for (e = 0; e < out.numberofedges; e++) {
575: if (out.edgemarkerlist[e]) {
576: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
577: const PetscInt *edges;
578: PetscInt numEdges;
580: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
581: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
582: if (glabel) {DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);}
583: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
584: }
585: }
586: #endif
587: for (f = 0; f < out.numberoftrifaces; f++) {
588: if (out.trifacemarkerlist[f]) {
589: const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
590: const PetscInt *faces;
591: PetscInt numFaces;
593: DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
594: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
595: if (glabel) {DMLabelSetValue(glabel, faces[0], out.trifacemarkerlist[f]);}
596: DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
597: }
598: }
599: }
600: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
601: }
602: return(0);
603: }
605: PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
606: {
607: MPI_Comm comm;
608: const PetscInt dim = 3;
609: const char *labelName = "marker";
610: ::tetgenio in;
611: ::tetgenio out;
612: DMLabel label;
613: PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
614: PetscMPIInt rank;
618: PetscObjectGetComm((PetscObject)dm,&comm);
619: MPI_Comm_rank(comm, &rank);
620: DMPlexGetDepth(dm, &depth);
621: MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
622: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
623: DMGetLabel(dm, labelName, &label);
625: in.numberofpoints = vEnd - vStart;
626: if (in.numberofpoints > 0) {
627: PetscSection coordSection;
628: Vec coordinates;
629: PetscScalar *array;
631: in.pointlist = new double[in.numberofpoints*dim];
632: in.pointmarkerlist = new int[in.numberofpoints];
634: DMGetCoordinatesLocal(dm, &coordinates);
635: DMGetCoordinateSection(dm, &coordSection);
636: VecGetArray(coordinates, &array);
637: for (v = vStart; v < vEnd; ++v) {
638: const PetscInt idx = v - vStart;
639: PetscInt off, d;
641: PetscSectionGetOffset(coordSection, v, &off);
642: for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
643: if (label) {
644: PetscInt val;
645:
646: DMLabelGetValue(label, v, &val);
647: in.pointmarkerlist[idx] = (int) val;
648: }
649: }
650: VecRestoreArray(coordinates, &array);
651: }
652: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
654: in.numberofcorners = 4;
655: in.numberoftetrahedra = cEnd - cStart;
656: in.tetrahedronvolumelist = (double*) maxVolumes;
657: if (in.numberoftetrahedra > 0) {
658: in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
659: for (c = cStart; c < cEnd; ++c) {
660: const PetscInt idx = c - cStart;
661: PetscInt *closure = NULL;
662: PetscInt closureSize;
664: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
665: if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
666: for (v = 0; v < 4; ++v) {
667: in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
668: }
669: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
670: }
671: }
672: /* TODO: Put in boundary faces with markers */
673: if (!rank) {
674: char args[32];
676: #if 1
677: /* Take away 'Q' for verbose output */
678: PetscStrcpy(args, "qezQra");
679: #else
680: PetscStrcpy(args, "qezraVVVV");
681: #endif
682: ::tetrahedralize(args, &in, &out);
683: }
684: in.tetrahedronvolumelist = NULL;
686: {
687: DMLabel rlabel = NULL;
688: const PetscInt numCorners = 4;
689: const PetscInt numCells = out.numberoftetrahedra;
690: const PetscInt numVertices = out.numberofpoints;
691: const double *meshCoords = out.pointlist;
692: int *cells = out.tetrahedronlist;
694: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
696: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
697: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
698: if (label) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*dmRefined, labelName, &rlabel);}
699: /* Set labels */
700: for (v = 0; v < numVertices; ++v) {
701: if (out.pointmarkerlist[v]) {
702: if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);}
703: }
704: }
705: if (interpolate) {
706: PetscInt f;
707: #if 0
708: PetscInt e;
710: for (e = 0; e < out.numberofedges; e++) {
711: if (out.edgemarkerlist[e]) {
712: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
713: const PetscInt *edges;
714: PetscInt numEdges;
716: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
717: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
718: if (rlabel) {DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);}
719: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
720: }
721: }
722: #endif
723: for (f = 0; f < out.numberoftrifaces; f++) {
724: if (out.trifacemarkerlist[f]) {
725: const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
726: const PetscInt *faces;
727: PetscInt numFaces;
729: DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
730: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
731: if (rlabel) {DMLabelSetValue(rlabel, faces[0], out.trifacemarkerlist[f]);}
732: DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
733: }
734: }
735: }
736: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
737: }
738: return(0);
739: }
740: #endif /* PETSC_HAVE_TETGEN */
742: #if defined(PETSC_HAVE_CTETGEN)
743: #include <ctetgen.h>
745: PetscErrorCode DMPlexGenerate_CTetgen(DM boundary, PetscBool interpolate, DM *dm)
746: {
747: MPI_Comm comm;
748: const PetscInt dim = 3;
749: const char *labelName = "marker";
750: PLC *in, *out;
751: DMLabel label;
752: PetscInt verbose = 0, vStart, vEnd, v, fStart, fEnd, f;
753: PetscMPIInt rank;
757: PetscObjectGetComm((PetscObject)boundary,&comm);
758: PetscOptionsGetInt(NULL,((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);
759: MPI_Comm_rank(comm, &rank);
760: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
761: DMGetLabel(boundary, labelName, &label);
762: PLCCreate(&in);
763: PLCCreate(&out);
765: in->numberofpoints = vEnd - vStart;
766: if (in->numberofpoints > 0) {
767: PetscSection coordSection;
768: Vec coordinates;
769: PetscScalar *array;
771: PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
772: PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);
773: DMGetCoordinatesLocal(boundary, &coordinates);
774: DMGetCoordinateSection(boundary, &coordSection);
775: VecGetArray(coordinates, &array);
776: for (v = vStart; v < vEnd; ++v) {
777: const PetscInt idx = v - vStart;
778: PetscInt off, d, m = -1;
780: PetscSectionGetOffset(coordSection, v, &off);
781: for (d = 0; d < dim; ++d) {
782: in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
783: }
784: if (label) {DMLabelGetValue(label, v, &m);}
786: in->pointmarkerlist[idx] = (int) m;
787: }
788: VecRestoreArray(coordinates, &array);
789: }
790: DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
792: in->numberoffacets = fEnd - fStart;
793: if (in->numberoffacets > 0) {
794: PetscMalloc1(in->numberoffacets, &in->facetlist);
795: PetscMalloc1(in->numberoffacets, &in->facetmarkerlist);
796: for (f = fStart; f < fEnd; ++f) {
797: const PetscInt idx = f - fStart;
798: PetscInt *points = NULL, numPoints, p, numVertices = 0, v, m = -1;
799: polygon *poly;
801: in->facetlist[idx].numberofpolygons = 1;
803: PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist);
805: in->facetlist[idx].numberofholes = 0;
806: in->facetlist[idx].holelist = NULL;
808: DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
809: for (p = 0; p < numPoints*2; p += 2) {
810: const PetscInt point = points[p];
811: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
812: }
814: poly = in->facetlist[idx].polygonlist;
815: poly->numberofvertices = numVertices;
816: PetscMalloc1(poly->numberofvertices, &poly->vertexlist);
817: for (v = 0; v < numVertices; ++v) {
818: const PetscInt vIdx = points[v] - vStart;
819: poly->vertexlist[v] = vIdx;
820: }
821: if (label) {DMLabelGetValue(label, f, &m);}
822: in->facetmarkerlist[idx] = (int) m;
823: DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
824: }
825: }
826: if (!rank) {
827: TetGenOpts t;
829: TetGenOptsInitialize(&t);
830: t.in = boundary; /* Should go away */
831: t.plc = 1;
832: t.quality = 1;
833: t.edgesout = 1;
834: t.zeroindex = 1;
835: t.quiet = 1;
836: t.verbose = verbose;
837: TetGenCheckOpts(&t);
838: TetGenTetrahedralize(&t, in, out);
839: }
840: {
841: DMLabel glabel = NULL;
842: const PetscInt numCorners = 4;
843: const PetscInt numCells = out->numberoftetrahedra;
844: const PetscInt numVertices = out->numberofpoints;
845: double *meshCoords;
846: int *cells = out->tetrahedronlist;
848: if (sizeof (PetscReal) == sizeof (double)) {
849: meshCoords = (double *) out->pointlist;
850: }
851: else {
852: PetscInt i;
854: PetscMalloc1(3 * numVertices,&meshCoords);
855: for (i = 0; i < 3 * numVertices; i++) {
856: meshCoords[i] = (PetscReal) out->pointlist[i];
857: }
858: }
860: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
861: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
862: if (sizeof (PetscReal) != sizeof (double)) {
863: PetscFree(meshCoords);
864: }
865: if (label) {DMCreateLabel(*dm, labelName); DMGetLabel(*dm, labelName, &glabel);}
866: /* Set labels */
867: for (v = 0; v < numVertices; ++v) {
868: if (out->pointmarkerlist[v]) {
869: if (glabel) {DMLabelSetValue(glabel, v+numCells, out->pointmarkerlist[v]);}
870: }
871: }
872: if (interpolate) {
873: PetscInt e;
875: for (e = 0; e < out->numberofedges; e++) {
876: if (out->edgemarkerlist[e]) {
877: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
878: const PetscInt *edges;
879: PetscInt numEdges;
881: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
882: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
883: if (glabel) {DMLabelSetValue(glabel, edges[0], out->edgemarkerlist[e]);}
884: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
885: }
886: }
887: for (f = 0; f < out->numberoftrifaces; f++) {
888: if (out->trifacemarkerlist[f]) {
889: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
890: const PetscInt *faces;
891: PetscInt numFaces;
893: DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
894: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
895: if (glabel) {DMLabelSetValue(glabel, faces[0], out->trifacemarkerlist[f]);}
896: DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
897: }
898: }
899: }
900: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
901: }
903: PLCDestroy(&in);
904: PLCDestroy(&out);
905: return(0);
906: }
908: PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
909: {
910: MPI_Comm comm;
911: const PetscInt dim = 3;
912: const char *labelName = "marker";
913: PLC *in, *out;
914: DMLabel label;
915: PetscInt verbose = 0, vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
916: PetscMPIInt rank;
920: PetscObjectGetComm((PetscObject)dm,&comm);
921: PetscOptionsGetInt(NULL,((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);
922: MPI_Comm_rank(comm, &rank);
923: DMPlexGetDepth(dm, &depth);
924: MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
925: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
926: DMGetLabel(dm, labelName, &label);
927: PLCCreate(&in);
928: PLCCreate(&out);
930: in->numberofpoints = vEnd - vStart;
931: if (in->numberofpoints > 0) {
932: PetscSection coordSection;
933: Vec coordinates;
934: PetscScalar *array;
936: PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
937: PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);
938: DMGetCoordinatesLocal(dm, &coordinates);
939: DMGetCoordinateSection(dm, &coordSection);
940: VecGetArray(coordinates, &array);
941: for (v = vStart; v < vEnd; ++v) {
942: const PetscInt idx = v - vStart;
943: PetscInt off, d, m = -1;
945: PetscSectionGetOffset(coordSection, v, &off);
946: for (d = 0; d < dim; ++d) {
947: in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
948: }
949: if (label) {DMLabelGetValue(label, v, &m);}
951: in->pointmarkerlist[idx] = (int) m;
952: }
953: VecRestoreArray(coordinates, &array);
954: }
955: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
957: in->numberofcorners = 4;
958: in->numberoftetrahedra = cEnd - cStart;
959: in->tetrahedronvolumelist = maxVolumes;
960: if (in->numberoftetrahedra > 0) {
961: PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);
962: for (c = cStart; c < cEnd; ++c) {
963: const PetscInt idx = c - cStart;
964: PetscInt *closure = NULL;
965: PetscInt closureSize;
967: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
968: if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
969: for (v = 0; v < 4; ++v) {
970: in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
971: }
972: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
973: }
974: }
975: if (!rank) {
976: TetGenOpts t;
978: TetGenOptsInitialize(&t);
980: t.in = dm; /* Should go away */
981: t.refine = 1;
982: t.varvolume = 1;
983: t.quality = 1;
984: t.edgesout = 1;
985: t.zeroindex = 1;
986: t.quiet = 1;
987: t.verbose = verbose; /* Change this */
989: TetGenCheckOpts(&t);
990: TetGenTetrahedralize(&t, in, out);
991: }
992: in->tetrahedronvolumelist = NULL;
993: {
994: DMLabel rlabel = NULL;
995: const PetscInt numCorners = 4;
996: const PetscInt numCells = out->numberoftetrahedra;
997: const PetscInt numVertices = out->numberofpoints;
998: double *meshCoords;
999: int *cells = out->tetrahedronlist;
1000: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
1002: if (sizeof (PetscReal) == sizeof (double)) {
1003: meshCoords = (double *) out->pointlist;
1004: }
1005: else {
1006: PetscInt i;
1008: PetscMalloc1(3 * numVertices,&meshCoords);
1009: for (i = 0; i < 3 * numVertices; i++) {
1010: meshCoords[i] = (PetscReal) out->pointlist[i];
1011: }
1012: }
1014: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
1015: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
1016: if (sizeof (PetscReal) != sizeof (double)) {
1017: PetscFree(meshCoords);
1018: }
1019: if (label) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*dmRefined, labelName, &rlabel);}
1020: /* Set labels */
1021: for (v = 0; v < numVertices; ++v) {
1022: if (out->pointmarkerlist[v]) {
1023: if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out->pointmarkerlist[v]);}
1024: }
1025: }
1026: if (interpolate) {
1027: PetscInt e, f;
1029: for (e = 0; e < out->numberofedges; e++) {
1030: if (out->edgemarkerlist[e]) {
1031: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
1032: const PetscInt *edges;
1033: PetscInt numEdges;
1035: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
1036: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
1037: if (rlabel) {DMLabelSetValue(rlabel, edges[0], out->edgemarkerlist[e]);}
1038: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
1039: }
1040: }
1041: for (f = 0; f < out->numberoftrifaces; f++) {
1042: if (out->trifacemarkerlist[f]) {
1043: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
1044: const PetscInt *faces;
1045: PetscInt numFaces;
1047: DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
1048: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
1049: if (rlabel) {DMLabelSetValue(rlabel, faces[0], out->trifacemarkerlist[f]);}
1050: DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
1051: }
1052: }
1053: }
1054: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
1055: }
1056: PLCDestroy(&in);
1057: PLCDestroy(&out);
1058: return(0);
1059: }
1060: #endif /* PETSC_HAVE_CTETGEN */
1062: /*@C
1063: DMPlexGenerate - Generates a mesh.
1065: Not Collective
1067: Input Parameters:
1068: + boundary - The DMPlex boundary object
1069: . name - The mesh generation package name
1070: - interpolate - Flag to create intermediate mesh elements
1072: Output Parameter:
1073: . mesh - The DMPlex object
1075: Level: intermediate
1077: .keywords: mesh, elements
1078: .seealso: DMPlexCreate(), DMRefine()
1079: @*/
1080: PetscErrorCode DMPlexGenerate(DM boundary, const char name[], PetscBool interpolate, DM *mesh)
1081: {
1082: PetscInt dim;
1083: char genname[1024];
1084: PetscBool isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;
1090: DMGetDimension(boundary, &dim);
1091: PetscOptionsGetString(((PetscObject) boundary)->options,((PetscObject) boundary)->prefix, "-dm_plex_generator", genname, 1024, &flg);
1092: if (flg) name = genname;
1093: if (name) {
1094: PetscStrcmp(name, "triangle", &isTriangle);
1095: PetscStrcmp(name, "tetgen", &isTetgen);
1096: PetscStrcmp(name, "ctetgen", &isCTetgen);
1097: }
1098: switch (dim) {
1099: case 1:
1100: if (!name || isTriangle) {
1101: #if defined(PETSC_HAVE_TRIANGLE)
1102: DMPlexGenerate_Triangle(boundary, interpolate, mesh);
1103: #else
1104: SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation needs external package support.\nPlease reconfigure with --download-triangle.");
1105: #endif
1106: } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
1107: break;
1108: case 2:
1109: if (!name) {
1110: #if defined(PETSC_HAVE_CTETGEN)
1111: DMPlexGenerate_CTetgen(boundary, interpolate, mesh);
1112: #elif defined(PETSC_HAVE_TETGEN)
1113: DMPlexGenerate_Tetgen(boundary, interpolate, mesh);
1114: #else
1115: SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "External package CTetgen or Tetgen needed.\nPlease reconfigure with '--download-ctetgen' or '--with-clanguage=cxx --download-tetgen'.");
1116: #endif
1117: } else if (isCTetgen) {
1118: #if defined(PETSC_HAVE_CTETGEN)
1119: DMPlexGenerate_CTetgen(boundary, interpolate, mesh);
1120: #else
1121: SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
1122: #endif
1123: } else if (isTetgen) {
1124: #if defined(PETSC_HAVE_TETGEN)
1125: DMPlexGenerate_Tetgen(boundary, interpolate, mesh);
1126: #else
1127: SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --download-tetgen.");
1128: #endif
1129: } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
1130: break;
1131: default:
1132: SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation for a dimension %d boundary is not supported.", dim);
1133: }
1134: return(0);
1135: }