Actual source code: trigenerate.c
1: #include <petsc/private/dmpleximpl.h>
3: #if !defined(ANSI_DECLARATORS)
4: #define ANSI_DECLARATORS
5: #endif
6: #include <triangle.h>
8: static PetscErrorCode InitInput_Triangle(struct triangulateio *inputCtx)
9: {
10: PetscFunctionBegin;
11: inputCtx->numberofpoints = 0;
12: inputCtx->numberofpointattributes = 0;
13: inputCtx->pointlist = NULL;
14: inputCtx->pointattributelist = NULL;
15: inputCtx->pointmarkerlist = NULL;
16: inputCtx->numberofsegments = 0;
17: inputCtx->segmentlist = NULL;
18: inputCtx->segmentmarkerlist = NULL;
19: inputCtx->numberoftriangleattributes = 0;
20: inputCtx->trianglelist = NULL;
21: inputCtx->numberofholes = 0;
22: inputCtx->holelist = NULL;
23: inputCtx->numberofregions = 0;
24: inputCtx->regionlist = NULL;
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: static PetscErrorCode InitOutput_Triangle(struct triangulateio *outputCtx)
29: {
30: PetscFunctionBegin;
31: outputCtx->numberofpoints = 0;
32: outputCtx->pointlist = NULL;
33: outputCtx->pointattributelist = NULL;
34: outputCtx->pointmarkerlist = NULL;
35: outputCtx->numberoftriangles = 0;
36: outputCtx->trianglelist = NULL;
37: outputCtx->triangleattributelist = NULL;
38: outputCtx->neighborlist = NULL;
39: outputCtx->segmentlist = NULL;
40: outputCtx->segmentmarkerlist = NULL;
41: outputCtx->numberofedges = 0;
42: outputCtx->edgelist = NULL;
43: outputCtx->edgemarkerlist = NULL;
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: static PetscErrorCode FiniOutput_Triangle(struct triangulateio *outputCtx)
48: {
49: PetscFunctionBegin;
50: free(outputCtx->pointlist);
51: free(outputCtx->pointmarkerlist);
52: free(outputCtx->segmentlist);
53: free(outputCtx->segmentmarkerlist);
54: free(outputCtx->edgelist);
55: free(outputCtx->edgemarkerlist);
56: free(outputCtx->trianglelist);
57: free(outputCtx->neighborlist);
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: PETSC_EXTERN PetscErrorCode DMPlexGenerate_Triangle(DM boundary, PetscBool interpolate, DM *dm)
62: {
63: MPI_Comm comm;
64: DM_Plex *mesh = (DM_Plex *)boundary->data;
65: PetscInt dim = 2;
66: const PetscBool createConvexHull = PETSC_FALSE;
67: const PetscBool constrained = PETSC_FALSE;
68: const char *labelName = "marker";
69: const char *labelName2 = "Face Sets";
70: struct triangulateio in;
71: struct triangulateio out;
72: DMLabel label, label2;
73: PetscInt vStart, vEnd, v, eStart, eEnd, e;
74: PetscMPIInt rank;
75: PetscBool flg;
76: char opts[64];
78: PetscFunctionBegin;
79: PetscCall(PetscObjectGetComm((PetscObject)boundary, &comm));
80: PetscCallMPI(MPI_Comm_rank(comm, &rank));
81: PetscCall(InitInput_Triangle(&in));
82: PetscCall(InitOutput_Triangle(&out));
83: PetscCall(DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd));
84: PetscCall(DMGetLabel(boundary, labelName, &label));
85: PetscCall(DMGetLabel(boundary, labelName2, &label2));
86: PetscCall(PetscOptionsGetString(((PetscObject)boundary)->options, ((PetscObject)boundary)->prefix, "-dm_plex_generate_triangle_opts", opts, sizeof(opts), &flg));
87: if (flg) PetscCall(DMPlexTriangleSetOptions(boundary, opts));
89: PetscCall(PetscCIntCast(vEnd - vStart, &in.numberofpoints));
90: if (in.numberofpoints > 0) {
91: PetscSection coordSection;
92: Vec coordinates;
93: PetscScalar *array;
95: PetscCall(PetscMalloc1(in.numberofpoints * dim, &in.pointlist));
96: PetscCall(PetscMalloc1(in.numberofpoints, &in.pointmarkerlist));
97: PetscCall(DMGetCoordinatesLocal(boundary, &coordinates));
98: PetscCall(DMGetCoordinateSection(boundary, &coordSection));
99: PetscCall(VecGetArray(coordinates, &array));
100: for (v = vStart; v < vEnd; ++v) {
101: const PetscInt idx = v - vStart;
102: PetscInt val, off, d;
104: PetscCall(PetscSectionGetOffset(coordSection, v, &off));
105: for (d = 0; d < dim; ++d) in.pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
106: if (label) {
107: PetscCall(DMLabelGetValue(label, v, &val));
108: PetscCall(PetscCIntCast(val, &in.pointmarkerlist[idx]));
109: }
110: }
111: PetscCall(VecRestoreArray(coordinates, &array));
112: }
113: PetscCall(DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd));
114: PetscCall(PetscCIntCast(eEnd - eStart, &in.numberofsegments));
115: if (in.numberofsegments > 0) {
116: PetscCall(PetscMalloc1(in.numberofsegments * 2, &in.segmentlist));
117: PetscCall(PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist));
118: for (e = eStart; e < eEnd; ++e) {
119: const PetscInt idx = e - eStart;
120: const PetscInt *cone;
121: PetscInt val;
123: PetscCall(DMPlexGetCone(boundary, e, &cone));
125: PetscCall(PetscCIntCast(cone[0] - vStart, &in.segmentlist[idx * 2 + 0]));
126: PetscCall(PetscCIntCast(cone[1] - vStart, &in.segmentlist[idx * 2 + 1]));
128: if (label) {
129: PetscCall(DMLabelGetValue(label, e, &val));
130: PetscCall(PetscCIntCast(val, &in.segmentmarkerlist[idx]));
131: }
132: }
133: }
134: #if 0 /* Do not currently support holes */
135: PetscReal *holeCoords;
136: PetscInt h, d;
138: PetscCall(DMPlexGetHoles(boundary, &in.numberofholes, &holeCords));
139: if (in.numberofholes > 0) {
140: PetscCall(PetscMalloc1(in.numberofholes*dim, &in.holelist));
141: for (h = 0; h < in.numberofholes; ++h) {
142: for (d = 0; d < dim; ++d) {
143: in.holelist[h*dim+d] = holeCoords[h*dim+d];
144: }
145: }
146: }
147: #endif
148: if (rank == 0) {
149: char args[32];
151: /* Take away 'Q' for verbose output */
152: PetscCall(PetscStrncpy(args, "pqezQ", sizeof(args)));
153: if (createConvexHull) PetscCall(PetscStrlcat(args, "c", sizeof(args)));
154: if (constrained) PetscCall(PetscStrncpy(args, "zepDQ", sizeof(args)));
155: if (mesh->triangleOpts) {
156: triangulate(mesh->triangleOpts, &in, &out, NULL);
157: } else {
158: triangulate(args, &in, &out, NULL);
159: }
160: }
161: PetscCall(PetscFree(in.pointlist));
162: PetscCall(PetscFree(in.pointmarkerlist));
163: PetscCall(PetscFree(in.segmentlist));
164: PetscCall(PetscFree(in.segmentmarkerlist));
165: PetscCall(PetscFree(in.holelist));
167: {
168: DMLabel glabel = NULL;
169: DMLabel glabel2 = NULL;
170: const PetscInt numCorners = 3;
171: const PetscInt numCells = out.numberoftriangles;
172: const PetscInt numVertices = out.numberofpoints;
173: PetscInt *cells;
174: PetscReal *meshCoords;
176: if (sizeof(PetscReal) == sizeof(out.pointlist[0])) {
177: meshCoords = (PetscReal *)out.pointlist;
178: } else {
179: PetscInt i;
181: PetscCall(PetscMalloc1(dim * numVertices, &meshCoords));
182: for (i = 0; i < dim * numVertices; i++) meshCoords[i] = (PetscReal)out.pointlist[i];
183: }
184: if (sizeof(PetscInt) == sizeof(out.trianglelist[0])) {
185: cells = (PetscInt *)out.trianglelist;
186: } else {
187: PetscInt i;
189: PetscCall(PetscMalloc1(numCells * numCorners, &cells));
190: for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out.trianglelist[i];
191: }
192: PetscCall(DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm));
193: if (sizeof(PetscReal) != sizeof(out.pointlist[0])) PetscCall(PetscFree(meshCoords));
194: if (sizeof(PetscInt) != sizeof(out.trianglelist[0])) PetscCall(PetscFree(cells));
195: if (label) {
196: PetscCall(DMCreateLabel(*dm, labelName));
197: PetscCall(DMGetLabel(*dm, labelName, &glabel));
198: }
199: if (label2) {
200: PetscCall(DMCreateLabel(*dm, labelName2));
201: PetscCall(DMGetLabel(*dm, labelName2, &glabel2));
202: }
203: /* Set labels */
204: for (v = 0; v < numVertices; ++v) {
205: if (out.pointmarkerlist[v]) {
206: if (glabel) PetscCall(DMLabelSetValue(glabel, v + numCells, out.pointmarkerlist[v]));
207: }
208: }
209: if (interpolate) {
210: for (e = 0; e < out.numberofedges; e++) {
211: if (out.edgemarkerlist[e]) {
212: const PetscInt vertices[2] = {out.edgelist[e * 2 + 0] + numCells, out.edgelist[e * 2 + 1] + numCells};
213: const PetscInt *edges;
214: PetscInt numEdges;
216: PetscCall(DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges));
217: PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %" PetscInt_FMT, numEdges);
218: if (glabel) PetscCall(DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]));
219: if (glabel2) PetscCall(DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]));
220: PetscCall(DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges));
221: }
222: }
223: }
224: PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_FALSE));
225: }
226: #if 0 /* Do not currently support holes */
227: PetscCall(DMPlexCopyHoles(*dm, boundary));
228: #endif
229: PetscCall(FiniOutput_Triangle(&out));
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: PETSC_EXTERN PetscErrorCode DMPlexRefine_Triangle(DM dm, PetscReal *inmaxVolumes, DM *dmRefined)
234: {
235: MPI_Comm comm;
236: PetscInt dim = 2;
237: const char *labelName = "marker";
238: struct triangulateio in;
239: struct triangulateio out;
240: DMLabel label;
241: PetscInt vStart, vEnd, v, gcStart, cStart, cEnd, c, depth, depthGlobal;
242: PetscMPIInt rank;
243: double *maxVolumes;
245: PetscFunctionBegin;
246: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
247: PetscCallMPI(MPI_Comm_rank(comm, &rank));
248: PetscCall(InitInput_Triangle(&in));
249: PetscCall(InitOutput_Triangle(&out));
250: PetscCall(DMPlexGetDepth(dm, &depth));
251: PetscCallMPI(MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm));
252: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
253: PetscCall(DMGetLabel(dm, labelName, &label));
255: PetscCall(PetscCIntCast(vEnd - vStart, &in.numberofpoints));
256: if (in.numberofpoints > 0) {
257: PetscSection coordSection;
258: Vec coordinates;
259: PetscScalar *array;
261: PetscCall(PetscMalloc1(in.numberofpoints * dim, &in.pointlist));
262: PetscCall(PetscMalloc1(in.numberofpoints, &in.pointmarkerlist));
263: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
264: PetscCall(DMGetCoordinateSection(dm, &coordSection));
265: PetscCall(VecGetArray(coordinates, &array));
266: for (v = vStart; v < vEnd; ++v) {
267: const PetscInt idx = v - vStart;
268: PetscInt off, d, val;
270: PetscCall(PetscSectionGetOffset(coordSection, v, &off));
271: for (d = 0; d < dim; ++d) in.pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
272: if (label) {
273: PetscCall(DMLabelGetValue(label, v, &val));
274: PetscCall(PetscCIntCast(val, &in.pointmarkerlist[idx]));
275: }
276: }
277: PetscCall(VecRestoreArray(coordinates, &array));
278: }
279: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
280: PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &gcStart, NULL));
281: if (gcStart >= 0) cEnd = gcStart;
283: in.numberofcorners = 3;
284: PetscCall(PetscCIntCast(cEnd - cStart, &in.numberoftriangles));
286: #if !defined(PETSC_USE_REAL_DOUBLE)
287: PetscCall(PetscMalloc1(cEnd - cStart, &maxVolumes));
288: for (c = 0; c < cEnd - cStart; ++c) maxVolumes[c] = (double)inmaxVolumes[c];
289: #else
290: maxVolumes = inmaxVolumes;
291: #endif
293: in.trianglearealist = (double *)maxVolumes;
294: if (in.numberoftriangles > 0) {
295: PetscCall(PetscMalloc1(in.numberoftriangles * in.numberofcorners, &in.trianglelist));
296: for (c = cStart; c < cEnd; ++c) {
297: const PetscInt idx = c - cStart;
298: PetscInt *closure = NULL;
299: PetscInt closureSize;
301: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
302: PetscCheck(!(closureSize != 4) || !(closureSize != 7), comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %" PetscInt_FMT " vertices in closure", closureSize);
303: for (v = 0; v < 3; ++v) PetscCall(PetscCIntCast(closure[(v + closureSize - 3) * 2] - vStart, &in.trianglelist[idx * in.numberofcorners + v]));
304: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
305: }
306: }
307: /* TODO: Segment markers are missing on input */
308: #if 0 /* Do not currently support holes */
309: PetscReal *holeCoords;
310: PetscInt h, d;
312: PetscCall(DMPlexGetHoles(boundary, &in.numberofholes, &holeCords));
313: if (in.numberofholes > 0) {
314: PetscCall(PetscMalloc1(in.numberofholes*dim, &in.holelist));
315: for (h = 0; h < in.numberofholes; ++h) {
316: for (d = 0; d < dim; ++d) {
317: in.holelist[h*dim+d] = holeCoords[h*dim+d];
318: }
319: }
320: }
321: #endif
322: if (rank == 0) {
323: char args[32];
325: /* Take away 'Q' for verbose output */
326: PetscCall(PetscStrncpy(args, "pqezQra", sizeof(args)));
327: triangulate(args, &in, &out, NULL);
328: }
329: PetscCall(PetscFree(in.pointlist));
330: PetscCall(PetscFree(in.pointmarkerlist));
331: PetscCall(PetscFree(in.segmentlist));
332: PetscCall(PetscFree(in.segmentmarkerlist));
333: PetscCall(PetscFree(in.trianglelist));
335: {
336: DMLabel rlabel = NULL;
337: const PetscInt numCorners = 3;
338: const PetscInt numCells = out.numberoftriangles;
339: const PetscInt numVertices = out.numberofpoints;
340: PetscInt *cells;
341: PetscReal *meshCoords;
342: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
344: if (sizeof(PetscReal) == sizeof(out.pointlist[0])) {
345: meshCoords = (PetscReal *)out.pointlist;
346: } else {
347: PetscInt i;
349: PetscCall(PetscMalloc1(dim * numVertices, &meshCoords));
350: for (i = 0; i < dim * numVertices; i++) meshCoords[i] = (PetscReal)out.pointlist[i];
351: }
352: if (sizeof(PetscInt) == sizeof(out.trianglelist[0])) {
353: cells = (PetscInt *)out.trianglelist;
354: } else {
355: PetscInt i;
357: PetscCall(PetscMalloc1(numCells * numCorners, &cells));
358: for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out.trianglelist[i];
359: }
361: PetscCall(DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined));
362: if (label) {
363: PetscCall(DMCreateLabel(*dmRefined, labelName));
364: PetscCall(DMGetLabel(*dmRefined, labelName, &rlabel));
365: }
366: if (sizeof(PetscReal) != sizeof(out.pointlist[0])) PetscCall(PetscFree(meshCoords));
367: if (sizeof(PetscInt) != sizeof(out.trianglelist[0])) PetscCall(PetscFree(cells));
368: /* Set labels */
369: for (v = 0; v < numVertices; ++v) {
370: if (out.pointmarkerlist[v]) {
371: if (rlabel) PetscCall(DMLabelSetValue(rlabel, v + numCells, out.pointmarkerlist[v]));
372: }
373: }
374: if (interpolate) {
375: PetscInt e;
377: for (e = 0; e < out.numberofedges; e++) {
378: if (out.edgemarkerlist[e]) {
379: const PetscInt vertices[2] = {out.edgelist[e * 2 + 0] + numCells, out.edgelist[e * 2 + 1] + numCells};
380: const PetscInt *edges;
381: PetscInt numEdges;
383: PetscCall(DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges));
384: PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %" PetscInt_FMT, numEdges);
385: if (rlabel) PetscCall(DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]));
386: PetscCall(DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges));
387: }
388: }
389: }
390: PetscCall(DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE));
391: }
392: #if 0 /* Do not currently support holes */
393: PetscCall(DMPlexCopyHoles(*dm, boundary));
394: #endif
395: PetscCall(FiniOutput_Triangle(&out));
396: #if !defined(PETSC_USE_REAL_DOUBLE)
397: PetscCall(PetscFree(maxVolumes));
398: #endif
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }