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;
76: PetscFunctionBegin;
77: PetscCall(PetscObjectGetComm((PetscObject)boundary, &comm));
78: PetscCallMPI(MPI_Comm_rank(comm, &rank));
79: PetscCall(InitInput_Triangle(&in));
80: PetscCall(InitOutput_Triangle(&out));
81: PetscCall(DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd));
82: PetscCall(DMGetLabel(boundary, labelName, &label));
83: PetscCall(DMGetLabel(boundary, labelName2, &label2));
85: in.numberofpoints = vEnd - vStart;
86: if (in.numberofpoints > 0) {
87: PetscSection coordSection;
88: Vec coordinates;
89: PetscScalar *array;
91: PetscCall(PetscMalloc1(in.numberofpoints * dim, &in.pointlist));
92: PetscCall(PetscMalloc1(in.numberofpoints, &in.pointmarkerlist));
93: PetscCall(DMGetCoordinatesLocal(boundary, &coordinates));
94: PetscCall(DMGetCoordinateSection(boundary, &coordSection));
95: PetscCall(VecGetArray(coordinates, &array));
96: for (v = vStart; v < vEnd; ++v) {
97: const PetscInt idx = v - vStart;
98: PetscInt val, off, d;
100: PetscCall(PetscSectionGetOffset(coordSection, v, &off));
101: for (d = 0; d < dim; ++d) in.pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
102: if (label) {
103: PetscCall(DMLabelGetValue(label, v, &val));
104: in.pointmarkerlist[idx] = val;
105: }
106: }
107: PetscCall(VecRestoreArray(coordinates, &array));
108: }
109: PetscCall(DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd));
110: in.numberofsegments = eEnd - eStart;
111: if (in.numberofsegments > 0) {
112: PetscCall(PetscMalloc1(in.numberofsegments * 2, &in.segmentlist));
113: PetscCall(PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist));
114: for (e = eStart; e < eEnd; ++e) {
115: const PetscInt idx = e - eStart;
116: const PetscInt *cone;
117: PetscInt val;
119: PetscCall(DMPlexGetCone(boundary, e, &cone));
121: in.segmentlist[idx * 2 + 0] = cone[0] - vStart;
122: in.segmentlist[idx * 2 + 1] = cone[1] - vStart;
124: if (label) {
125: PetscCall(DMLabelGetValue(label, e, &val));
126: in.segmentmarkerlist[idx] = val;
127: }
128: }
129: }
130: #if 0 /* Do not currently support holes */
131: PetscReal *holeCoords;
132: PetscInt h, d;
134: PetscCall(DMPlexGetHoles(boundary, &in.numberofholes, &holeCords));
135: if (in.numberofholes > 0) {
136: PetscCall(PetscMalloc1(in.numberofholes*dim, &in.holelist));
137: for (h = 0; h < in.numberofholes; ++h) {
138: for (d = 0; d < dim; ++d) {
139: in.holelist[h*dim+d] = holeCoords[h*dim+d];
140: }
141: }
142: }
143: #endif
144: if (rank == 0) {
145: char args[32];
147: /* Take away 'Q' for verbose output */
148: PetscCall(PetscStrncpy(args, "pqezQ", sizeof(args)));
149: if (createConvexHull) PetscCall(PetscStrlcat(args, "c", sizeof(args)));
150: if (constrained) PetscCall(PetscStrncpy(args, "zepDQ", sizeof(args)));
151: if (mesh->triangleOpts) {
152: triangulate(mesh->triangleOpts, &in, &out, NULL);
153: } else {
154: triangulate(args, &in, &out, NULL);
155: }
156: }
157: PetscCall(PetscFree(in.pointlist));
158: PetscCall(PetscFree(in.pointmarkerlist));
159: PetscCall(PetscFree(in.segmentlist));
160: PetscCall(PetscFree(in.segmentmarkerlist));
161: PetscCall(PetscFree(in.holelist));
163: {
164: DMLabel glabel = NULL;
165: DMLabel glabel2 = NULL;
166: const PetscInt numCorners = 3;
167: const PetscInt numCells = out.numberoftriangles;
168: const PetscInt numVertices = out.numberofpoints;
169: PetscInt *cells;
170: PetscReal *meshCoords;
172: if (sizeof(PetscReal) == sizeof(out.pointlist[0])) {
173: meshCoords = (PetscReal *)out.pointlist;
174: } else {
175: PetscInt i;
177: PetscCall(PetscMalloc1(dim * numVertices, &meshCoords));
178: for (i = 0; i < dim * numVertices; i++) meshCoords[i] = (PetscReal)out.pointlist[i];
179: }
180: if (sizeof(PetscInt) == sizeof(out.trianglelist[0])) {
181: cells = (PetscInt *)out.trianglelist;
182: } else {
183: PetscInt i;
185: PetscCall(PetscMalloc1(numCells * numCorners, &cells));
186: for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out.trianglelist[i];
187: }
188: PetscCall(DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm));
189: if (sizeof(PetscReal) != sizeof(out.pointlist[0])) PetscCall(PetscFree(meshCoords));
190: if (sizeof(PetscInt) != sizeof(out.trianglelist[0])) PetscCall(PetscFree(cells));
191: if (label) {
192: PetscCall(DMCreateLabel(*dm, labelName));
193: PetscCall(DMGetLabel(*dm, labelName, &glabel));
194: }
195: if (label2) {
196: PetscCall(DMCreateLabel(*dm, labelName2));
197: PetscCall(DMGetLabel(*dm, labelName2, &glabel2));
198: }
199: /* Set labels */
200: for (v = 0; v < numVertices; ++v) {
201: if (out.pointmarkerlist[v]) {
202: if (glabel) PetscCall(DMLabelSetValue(glabel, v + numCells, out.pointmarkerlist[v]));
203: }
204: }
205: if (interpolate) {
206: for (e = 0; e < out.numberofedges; e++) {
207: if (out.edgemarkerlist[e]) {
208: const PetscInt vertices[2] = {out.edgelist[e * 2 + 0] + numCells, out.edgelist[e * 2 + 1] + numCells};
209: const PetscInt *edges;
210: PetscInt numEdges;
212: PetscCall(DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges));
213: PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %" PetscInt_FMT, numEdges);
214: if (glabel) PetscCall(DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]));
215: if (glabel2) PetscCall(DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]));
216: PetscCall(DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges));
217: }
218: }
219: }
220: PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_FALSE));
221: }
222: #if 0 /* Do not currently support holes */
223: PetscCall(DMPlexCopyHoles(*dm, boundary));
224: #endif
225: PetscCall(FiniOutput_Triangle(&out));
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: PETSC_EXTERN PetscErrorCode DMPlexRefine_Triangle(DM dm, PetscReal *inmaxVolumes, DM *dmRefined)
230: {
231: MPI_Comm comm;
232: PetscInt dim = 2;
233: const char *labelName = "marker";
234: struct triangulateio in;
235: struct triangulateio out;
236: DMLabel label;
237: PetscInt vStart, vEnd, v, gcStart, cStart, cEnd, c, depth, depthGlobal;
238: PetscMPIInt rank;
239: double *maxVolumes;
241: PetscFunctionBegin;
242: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
243: PetscCallMPI(MPI_Comm_rank(comm, &rank));
244: PetscCall(InitInput_Triangle(&in));
245: PetscCall(InitOutput_Triangle(&out));
246: PetscCall(DMPlexGetDepth(dm, &depth));
247: PetscCall(MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm));
248: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
249: PetscCall(DMGetLabel(dm, labelName, &label));
251: in.numberofpoints = vEnd - vStart;
252: if (in.numberofpoints > 0) {
253: PetscSection coordSection;
254: Vec coordinates;
255: PetscScalar *array;
257: PetscCall(PetscMalloc1(in.numberofpoints * dim, &in.pointlist));
258: PetscCall(PetscMalloc1(in.numberofpoints, &in.pointmarkerlist));
259: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
260: PetscCall(DMGetCoordinateSection(dm, &coordSection));
261: PetscCall(VecGetArray(coordinates, &array));
262: for (v = vStart; v < vEnd; ++v) {
263: const PetscInt idx = v - vStart;
264: PetscInt off, d, val;
266: PetscCall(PetscSectionGetOffset(coordSection, v, &off));
267: for (d = 0; d < dim; ++d) in.pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
268: if (label) {
269: PetscCall(DMLabelGetValue(label, v, &val));
270: in.pointmarkerlist[idx] = val;
271: }
272: }
273: PetscCall(VecRestoreArray(coordinates, &array));
274: }
275: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
276: PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &gcStart, NULL));
277: if (gcStart >= 0) cEnd = gcStart;
279: in.numberofcorners = 3;
280: in.numberoftriangles = cEnd - cStart;
282: #if !defined(PETSC_USE_REAL_DOUBLE)
283: PetscCall(PetscMalloc1(cEnd - cStart, &maxVolumes));
284: for (c = 0; c < cEnd - cStart; ++c) maxVolumes[c] = (double)inmaxVolumes[c];
285: #else
286: maxVolumes = inmaxVolumes;
287: #endif
289: in.trianglearealist = (double *)maxVolumes;
290: if (in.numberoftriangles > 0) {
291: PetscCall(PetscMalloc1(in.numberoftriangles * in.numberofcorners, &in.trianglelist));
292: for (c = cStart; c < cEnd; ++c) {
293: const PetscInt idx = c - cStart;
294: PetscInt *closure = NULL;
295: PetscInt closureSize;
297: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
298: PetscCheck(!(closureSize != 4) || !(closureSize != 7), comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %" PetscInt_FMT " vertices in closure", closureSize);
299: for (v = 0; v < 3; ++v) in.trianglelist[idx * in.numberofcorners + v] = closure[(v + closureSize - 3) * 2] - vStart;
300: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
301: }
302: }
303: /* TODO: Segment markers are missing on input */
304: #if 0 /* Do not currently support holes */
305: PetscReal *holeCoords;
306: PetscInt h, d;
308: PetscCall(DMPlexGetHoles(boundary, &in.numberofholes, &holeCords));
309: if (in.numberofholes > 0) {
310: PetscCall(PetscMalloc1(in.numberofholes*dim, &in.holelist));
311: for (h = 0; h < in.numberofholes; ++h) {
312: for (d = 0; d < dim; ++d) {
313: in.holelist[h*dim+d] = holeCoords[h*dim+d];
314: }
315: }
316: }
317: #endif
318: if (rank == 0) {
319: char args[32];
321: /* Take away 'Q' for verbose output */
322: PetscCall(PetscStrncpy(args, "pqezQra", sizeof(args)));
323: triangulate(args, &in, &out, NULL);
324: }
325: PetscCall(PetscFree(in.pointlist));
326: PetscCall(PetscFree(in.pointmarkerlist));
327: PetscCall(PetscFree(in.segmentlist));
328: PetscCall(PetscFree(in.segmentmarkerlist));
329: PetscCall(PetscFree(in.trianglelist));
331: {
332: DMLabel rlabel = NULL;
333: const PetscInt numCorners = 3;
334: const PetscInt numCells = out.numberoftriangles;
335: const PetscInt numVertices = out.numberofpoints;
336: PetscInt *cells;
337: PetscReal *meshCoords;
338: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
340: if (sizeof(PetscReal) == sizeof(out.pointlist[0])) {
341: meshCoords = (PetscReal *)out.pointlist;
342: } else {
343: PetscInt i;
345: PetscCall(PetscMalloc1(dim * numVertices, &meshCoords));
346: for (i = 0; i < dim * numVertices; i++) meshCoords[i] = (PetscReal)out.pointlist[i];
347: }
348: if (sizeof(PetscInt) == sizeof(out.trianglelist[0])) {
349: cells = (PetscInt *)out.trianglelist;
350: } else {
351: PetscInt i;
353: PetscCall(PetscMalloc1(numCells * numCorners, &cells));
354: for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out.trianglelist[i];
355: }
357: PetscCall(DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined));
358: if (label) {
359: PetscCall(DMCreateLabel(*dmRefined, labelName));
360: PetscCall(DMGetLabel(*dmRefined, labelName, &rlabel));
361: }
362: if (sizeof(PetscReal) != sizeof(out.pointlist[0])) PetscCall(PetscFree(meshCoords));
363: if (sizeof(PetscInt) != sizeof(out.trianglelist[0])) PetscCall(PetscFree(cells));
364: /* Set labels */
365: for (v = 0; v < numVertices; ++v) {
366: if (out.pointmarkerlist[v]) {
367: if (rlabel) PetscCall(DMLabelSetValue(rlabel, v + numCells, out.pointmarkerlist[v]));
368: }
369: }
370: if (interpolate) {
371: PetscInt e;
373: for (e = 0; e < out.numberofedges; e++) {
374: if (out.edgemarkerlist[e]) {
375: const PetscInt vertices[2] = {out.edgelist[e * 2 + 0] + numCells, out.edgelist[e * 2 + 1] + numCells};
376: const PetscInt *edges;
377: PetscInt numEdges;
379: PetscCall(DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges));
380: PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %" PetscInt_FMT, numEdges);
381: if (rlabel) PetscCall(DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]));
382: PetscCall(DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges));
383: }
384: }
385: }
386: PetscCall(DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE));
387: }
388: #if 0 /* Do not currently support holes */
389: PetscCall(DMPlexCopyHoles(*dm, boundary));
390: #endif
391: PetscCall(FiniOutput_Triangle(&out));
392: #if !defined(PETSC_USE_REAL_DOUBLE)
393: PetscCall(PetscFree(maxVolumes));
394: #endif
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }