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: inputCtx->numberofpoints = 0;
11: inputCtx->numberofpointattributes = 0;
12: inputCtx->pointlist = NULL;
13: inputCtx->pointattributelist = NULL;
14: inputCtx->pointmarkerlist = NULL;
15: inputCtx->numberofsegments = 0;
16: inputCtx->segmentlist = NULL;
17: inputCtx->segmentmarkerlist = NULL;
18: inputCtx->numberoftriangleattributes = 0;
19: inputCtx->trianglelist = NULL;
20: inputCtx->numberofholes = 0;
21: inputCtx->holelist = NULL;
22: inputCtx->numberofregions = 0;
23: inputCtx->regionlist = NULL;
24: return 0;
25: }
27: static PetscErrorCode InitOutput_Triangle(struct triangulateio *outputCtx)
28: {
29: outputCtx->numberofpoints = 0;
30: outputCtx->pointlist = NULL;
31: outputCtx->pointattributelist = NULL;
32: outputCtx->pointmarkerlist = NULL;
33: outputCtx->numberoftriangles = 0;
34: outputCtx->trianglelist = NULL;
35: outputCtx->triangleattributelist = NULL;
36: outputCtx->neighborlist = NULL;
37: outputCtx->segmentlist = NULL;
38: outputCtx->segmentmarkerlist = NULL;
39: outputCtx->numberofedges = 0;
40: outputCtx->edgelist = NULL;
41: outputCtx->edgemarkerlist = NULL;
42: return 0;
43: }
45: static PetscErrorCode FiniOutput_Triangle(struct triangulateio *outputCtx)
46: {
47: free(outputCtx->pointlist);
48: free(outputCtx->pointmarkerlist);
49: free(outputCtx->segmentlist);
50: free(outputCtx->segmentmarkerlist);
51: free(outputCtx->edgelist);
52: free(outputCtx->edgemarkerlist);
53: free(outputCtx->trianglelist);
54: free(outputCtx->neighborlist);
55: return 0;
56: }
58: PETSC_EXTERN PetscErrorCode DMPlexGenerate_Triangle(DM boundary, PetscBool interpolate, DM *dm)
59: {
60: MPI_Comm comm;
61: DM_Plex *mesh = (DM_Plex *)boundary->data;
62: PetscInt dim = 2;
63: const PetscBool createConvexHull = PETSC_FALSE;
64: const PetscBool constrained = PETSC_FALSE;
65: const char *labelName = "marker";
66: const char *labelName2 = "Face Sets";
67: struct triangulateio in;
68: struct triangulateio out;
69: DMLabel label, label2;
70: PetscInt vStart, vEnd, v, eStart, eEnd, e;
71: PetscMPIInt rank;
73: PetscObjectGetComm((PetscObject)boundary, &comm);
74: MPI_Comm_rank(comm, &rank);
75: InitInput_Triangle(&in);
76: InitOutput_Triangle(&out);
77: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
78: DMGetLabel(boundary, labelName, &label);
79: DMGetLabel(boundary, labelName2, &label2);
81: in.numberofpoints = vEnd - vStart;
82: if (in.numberofpoints > 0) {
83: PetscSection coordSection;
84: Vec coordinates;
85: PetscScalar *array;
87: PetscMalloc1(in.numberofpoints * dim, &in.pointlist);
88: PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
89: DMGetCoordinatesLocal(boundary, &coordinates);
90: DMGetCoordinateSection(boundary, &coordSection);
91: VecGetArray(coordinates, &array);
92: for (v = vStart; v < vEnd; ++v) {
93: const PetscInt idx = v - vStart;
94: PetscInt val, off, d;
96: PetscSectionGetOffset(coordSection, v, &off);
97: for (d = 0; d < dim; ++d) in.pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
98: if (label) {
99: DMLabelGetValue(label, v, &val);
100: in.pointmarkerlist[idx] = val;
101: }
102: }
103: VecRestoreArray(coordinates, &array);
104: }
105: DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);
106: in.numberofsegments = eEnd - eStart;
107: if (in.numberofsegments > 0) {
108: PetscMalloc1(in.numberofsegments * 2, &in.segmentlist);
109: PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist);
110: for (e = eStart; e < eEnd; ++e) {
111: const PetscInt idx = e - eStart;
112: const PetscInt *cone;
113: PetscInt val;
115: DMPlexGetCone(boundary, e, &cone);
117: in.segmentlist[idx * 2 + 0] = cone[0] - vStart;
118: in.segmentlist[idx * 2 + 1] = cone[1] - vStart;
120: if (label) {
121: DMLabelGetValue(label, e, &val);
122: in.segmentmarkerlist[idx] = val;
123: }
124: }
125: }
126: #if 0 /* Do not currently support holes */
127: PetscReal *holeCoords;
128: PetscInt h, d;
130: DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
131: if (in.numberofholes > 0) {
132: PetscMalloc1(in.numberofholes*dim, &in.holelist);
133: for (h = 0; h < in.numberofholes; ++h) {
134: for (d = 0; d < dim; ++d) {
135: in.holelist[h*dim+d] = holeCoords[h*dim+d];
136: }
137: }
138: }
139: #endif
140: if (rank == 0) {
141: char args[32];
143: /* Take away 'Q' for verbose output */
144: PetscStrcpy(args, "pqezQ");
145: if (createConvexHull) PetscStrcat(args, "c");
146: if (constrained) PetscStrcpy(args, "zepDQ");
147: if (mesh->triangleOpts) {
148: triangulate(mesh->triangleOpts, &in, &out, NULL);
149: } else {
150: triangulate(args, &in, &out, NULL);
151: }
152: }
153: PetscFree(in.pointlist);
154: PetscFree(in.pointmarkerlist);
155: PetscFree(in.segmentlist);
156: PetscFree(in.segmentmarkerlist);
157: PetscFree(in.holelist);
159: {
160: DMLabel glabel = NULL;
161: DMLabel glabel2 = NULL;
162: const PetscInt numCorners = 3;
163: const PetscInt numCells = out.numberoftriangles;
164: const PetscInt numVertices = out.numberofpoints;
165: PetscInt *cells;
166: PetscReal *meshCoords;
168: if (sizeof(PetscReal) == sizeof(out.pointlist[0])) {
169: meshCoords = (PetscReal *)out.pointlist;
170: } else {
171: PetscInt i;
173: PetscMalloc1(dim * numVertices, &meshCoords);
174: for (i = 0; i < dim * numVertices; i++) meshCoords[i] = (PetscReal)out.pointlist[i];
175: }
176: if (sizeof(PetscInt) == sizeof(out.trianglelist[0])) {
177: cells = (PetscInt *)out.trianglelist;
178: } else {
179: PetscInt i;
181: PetscMalloc1(numCells * numCorners, &cells);
182: for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out.trianglelist[i];
183: }
184: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
185: if (sizeof(PetscReal) != sizeof(out.pointlist[0])) PetscFree(meshCoords);
186: if (sizeof(PetscInt) != sizeof(out.trianglelist[0])) PetscFree(cells);
187: if (label) {
188: DMCreateLabel(*dm, labelName);
189: DMGetLabel(*dm, labelName, &glabel);
190: }
191: if (label2) {
192: DMCreateLabel(*dm, labelName2);
193: DMGetLabel(*dm, labelName2, &glabel2);
194: }
195: /* Set labels */
196: for (v = 0; v < numVertices; ++v) {
197: if (out.pointmarkerlist[v]) {
198: if (glabel) DMLabelSetValue(glabel, v + numCells, out.pointmarkerlist[v]);
199: }
200: }
201: if (interpolate) {
202: for (e = 0; e < out.numberofedges; e++) {
203: if (out.edgemarkerlist[e]) {
204: const PetscInt vertices[2] = {out.edgelist[e * 2 + 0] + numCells, out.edgelist[e * 2 + 1] + numCells};
205: const PetscInt *edges;
206: PetscInt numEdges;
208: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
210: if (glabel) DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);
211: if (glabel2) DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]);
212: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
213: }
214: }
215: }
216: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
217: }
218: #if 0 /* Do not currently support holes */
219: DMPlexCopyHoles(*dm, boundary);
220: #endif
221: FiniOutput_Triangle(&out);
222: return 0;
223: }
225: PETSC_EXTERN PetscErrorCode DMPlexRefine_Triangle(DM dm, PetscReal *inmaxVolumes, DM *dmRefined)
226: {
227: MPI_Comm comm;
228: PetscInt dim = 2;
229: const char *labelName = "marker";
230: struct triangulateio in;
231: struct triangulateio out;
232: DMLabel label;
233: PetscInt vStart, vEnd, v, gcStart, cStart, cEnd, c, depth, depthGlobal;
234: PetscMPIInt rank;
235: double *maxVolumes;
237: PetscObjectGetComm((PetscObject)dm, &comm);
238: MPI_Comm_rank(comm, &rank);
239: InitInput_Triangle(&in);
240: InitOutput_Triangle(&out);
241: DMPlexGetDepth(dm, &depth);
242: MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
243: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
244: DMGetLabel(dm, labelName, &label);
246: in.numberofpoints = vEnd - vStart;
247: if (in.numberofpoints > 0) {
248: PetscSection coordSection;
249: Vec coordinates;
250: PetscScalar *array;
252: PetscMalloc1(in.numberofpoints * dim, &in.pointlist);
253: PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
254: DMGetCoordinatesLocal(dm, &coordinates);
255: DMGetCoordinateSection(dm, &coordSection);
256: VecGetArray(coordinates, &array);
257: for (v = vStart; v < vEnd; ++v) {
258: const PetscInt idx = v - vStart;
259: PetscInt off, d, val;
261: PetscSectionGetOffset(coordSection, v, &off);
262: for (d = 0; d < dim; ++d) in.pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
263: if (label) {
264: DMLabelGetValue(label, v, &val);
265: in.pointmarkerlist[idx] = val;
266: }
267: }
268: VecRestoreArray(coordinates, &array);
269: }
270: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
271: DMPlexGetGhostCellStratum(dm, &gcStart, NULL);
272: if (gcStart >= 0) cEnd = gcStart;
274: in.numberofcorners = 3;
275: in.numberoftriangles = cEnd - cStart;
277: #if !defined(PETSC_USE_REAL_DOUBLE)
278: PetscMalloc1(cEnd - cStart, &maxVolumes);
279: for (c = 0; c < cEnd - cStart; ++c) maxVolumes[c] = (double)inmaxVolumes[c];
280: #else
281: maxVolumes = inmaxVolumes;
282: #endif
284: in.trianglearealist = (double *)maxVolumes;
285: if (in.numberoftriangles > 0) {
286: PetscMalloc1(in.numberoftriangles * in.numberofcorners, &in.trianglelist);
287: for (c = cStart; c < cEnd; ++c) {
288: const PetscInt idx = c - cStart;
289: PetscInt *closure = NULL;
290: PetscInt closureSize;
292: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
294: for (v = 0; v < 3; ++v) in.trianglelist[idx * in.numberofcorners + v] = closure[(v + closureSize - 3) * 2] - vStart;
295: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
296: }
297: }
298: /* TODO: Segment markers are missing on input */
299: #if 0 /* Do not currently support holes */
300: PetscReal *holeCoords;
301: PetscInt h, d;
303: DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
304: if (in.numberofholes > 0) {
305: PetscMalloc1(in.numberofholes*dim, &in.holelist);
306: for (h = 0; h < in.numberofholes; ++h) {
307: for (d = 0; d < dim; ++d) {
308: in.holelist[h*dim+d] = holeCoords[h*dim+d];
309: }
310: }
311: }
312: #endif
313: if (rank == 0) {
314: char args[32];
316: /* Take away 'Q' for verbose output */
317: PetscStrcpy(args, "pqezQra");
318: triangulate(args, &in, &out, NULL);
319: }
320: PetscFree(in.pointlist);
321: PetscFree(in.pointmarkerlist);
322: PetscFree(in.segmentlist);
323: PetscFree(in.segmentmarkerlist);
324: PetscFree(in.trianglelist);
326: {
327: DMLabel rlabel = NULL;
328: const PetscInt numCorners = 3;
329: const PetscInt numCells = out.numberoftriangles;
330: const PetscInt numVertices = out.numberofpoints;
331: PetscInt *cells;
332: PetscReal *meshCoords;
333: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
335: if (sizeof(PetscReal) == sizeof(out.pointlist[0])) {
336: meshCoords = (PetscReal *)out.pointlist;
337: } else {
338: PetscInt i;
340: PetscMalloc1(dim * numVertices, &meshCoords);
341: for (i = 0; i < dim * numVertices; i++) meshCoords[i] = (PetscReal)out.pointlist[i];
342: }
343: if (sizeof(PetscInt) == sizeof(out.trianglelist[0])) {
344: cells = (PetscInt *)out.trianglelist;
345: } else {
346: PetscInt i;
348: PetscMalloc1(numCells * numCorners, &cells);
349: for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out.trianglelist[i];
350: }
352: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
353: if (label) {
354: DMCreateLabel(*dmRefined, labelName);
355: DMGetLabel(*dmRefined, labelName, &rlabel);
356: }
357: if (sizeof(PetscReal) != sizeof(out.pointlist[0])) PetscFree(meshCoords);
358: if (sizeof(PetscInt) != sizeof(out.trianglelist[0])) PetscFree(cells);
359: /* Set labels */
360: for (v = 0; v < numVertices; ++v) {
361: if (out.pointmarkerlist[v]) {
362: if (rlabel) DMLabelSetValue(rlabel, v + numCells, out.pointmarkerlist[v]);
363: }
364: }
365: if (interpolate) {
366: PetscInt e;
368: for (e = 0; e < out.numberofedges; e++) {
369: if (out.edgemarkerlist[e]) {
370: const PetscInt vertices[2] = {out.edgelist[e * 2 + 0] + numCells, out.edgelist[e * 2 + 1] + numCells};
371: const PetscInt *edges;
372: PetscInt numEdges;
374: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
376: if (rlabel) DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);
377: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
378: }
379: }
380: }
381: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
382: }
383: #if 0 /* Do not currently support holes */
384: DMPlexCopyHoles(*dm, boundary);
385: #endif
386: FiniOutput_Triangle(&out);
387: #if !defined(PETSC_USE_REAL_DOUBLE)
388: PetscFree(maxVolumes);
389: #endif
390: return 0;
391: }