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) {
98: in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
99: }
100: if (label) {
101: DMLabelGetValue(label, v, &val);
102: in.pointmarkerlist[idx] = val;
103: }
104: }
105: VecRestoreArray(coordinates, &array);
106: }
107: DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);
108: in.numberofsegments = eEnd - eStart;
109: if (in.numberofsegments > 0) {
110: PetscMalloc1(in.numberofsegments*2, &in.segmentlist);
111: PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist);
112: for (e = eStart; e < eEnd; ++e) {
113: const PetscInt idx = e - eStart;
114: const PetscInt *cone;
115: PetscInt val;
117: DMPlexGetCone(boundary, e, &cone);
119: in.segmentlist[idx*2+0] = cone[0] - vStart;
120: in.segmentlist[idx*2+1] = cone[1] - vStart;
122: if (label) {
123: DMLabelGetValue(label, e, &val);
124: in.segmentmarkerlist[idx] = val;
125: }
126: }
127: }
128: #if 0 /* Do not currently support holes */
129: PetscReal *holeCoords;
130: PetscInt h, d;
132: DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
133: if (in.numberofholes > 0) {
134: PetscMalloc1(in.numberofholes*dim, &in.holelist);
135: for (h = 0; h < in.numberofholes; ++h) {
136: for (d = 0; d < dim; ++d) {
137: in.holelist[h*dim+d] = holeCoords[h*dim+d];
138: }
139: }
140: }
141: #endif
142: if (rank == 0) {
143: char args[32];
145: /* Take away 'Q' for verbose output */
146: PetscStrcpy(args, "pqezQ");
147: if (createConvexHull) PetscStrcat(args, "c");
148: if (constrained) PetscStrcpy(args, "zepDQ");
149: if (mesh->triangleOpts) {triangulate(mesh->triangleOpts, &in, &out, NULL);}
150: else {triangulate(args, &in, &out, NULL);}
151: }
152: PetscFree(in.pointlist);
153: PetscFree(in.pointmarkerlist);
154: PetscFree(in.segmentlist);
155: PetscFree(in.segmentmarkerlist);
156: PetscFree(in.holelist);
158: {
159: DMLabel glabel = NULL;
160: DMLabel glabel2 = NULL;
161: const PetscInt numCorners = 3;
162: const PetscInt numCells = out.numberoftriangles;
163: const PetscInt numVertices = out.numberofpoints;
164: PetscInt *cells;
165: PetscReal *meshCoords;
167: if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
168: meshCoords = (PetscReal *) out.pointlist;
169: } else {
170: PetscInt i;
172: PetscMalloc1(dim * numVertices,&meshCoords);
173: for (i = 0; i < dim * numVertices; i++) {
174: meshCoords[i] = (PetscReal) out.pointlist[i];
175: }
176: }
177: if (sizeof (PetscInt) == sizeof (out.trianglelist[0])) {
178: cells = (PetscInt *) out.trianglelist;
179: } else {
180: PetscInt i;
182: PetscMalloc1(numCells * numCorners, &cells);
183: for (i = 0; i < numCells * numCorners; i++) {
184: cells[i] = (PetscInt) out.trianglelist[i];
185: }
186: }
187: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
188: if (sizeof (PetscReal) != sizeof (out.pointlist[0])) {
189: PetscFree(meshCoords);
190: }
191: if (sizeof (PetscInt) != sizeof (out.trianglelist[0])) {
192: PetscFree(cells);
193: }
194: if (label) {
195: DMCreateLabel(*dm, labelName);
196: DMGetLabel(*dm, labelName, &glabel);
197: }
198: if (label2) {
199: DMCreateLabel(*dm, labelName2);
200: DMGetLabel(*dm, labelName2, &glabel2);
201: }
202: /* Set labels */
203: for (v = 0; v < numVertices; ++v) {
204: if (out.pointmarkerlist[v]) {
205: if (glabel) DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);
206: }
207: }
208: if (interpolate) {
209: for (e = 0; e < out.numberofedges; e++) {
210: if (out.edgemarkerlist[e]) {
211: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
212: const PetscInt *edges;
213: PetscInt numEdges;
215: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
217: if (glabel) DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);
218: if (glabel2) DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]);
219: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
220: }
221: }
222: }
223: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
224: }
225: #if 0 /* Do not currently support holes */
226: DMPlexCopyHoles(*dm, boundary);
227: #endif
228: FiniOutput_Triangle(&out);
229: return 0;
230: }
232: PETSC_EXTERN PetscErrorCode DMPlexRefine_Triangle(DM dm, PetscReal *inmaxVolumes, DM *dmRefined)
233: {
234: MPI_Comm comm;
235: PetscInt dim = 2;
236: const char *labelName = "marker";
237: struct triangulateio in;
238: struct triangulateio out;
239: DMLabel label;
240: PetscInt vStart, vEnd, v, gcStart, cStart, cEnd, c, depth, depthGlobal;
241: PetscMPIInt rank;
242: double *maxVolumes;
244: PetscObjectGetComm((PetscObject)dm,&comm);
245: MPI_Comm_rank(comm, &rank);
246: InitInput_Triangle(&in);
247: InitOutput_Triangle(&out);
248: DMPlexGetDepth(dm, &depth);
249: MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
250: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
251: DMGetLabel(dm, labelName, &label);
253: in.numberofpoints = vEnd - vStart;
254: if (in.numberofpoints > 0) {
255: PetscSection coordSection;
256: Vec coordinates;
257: PetscScalar *array;
259: PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
260: PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
261: DMGetCoordinatesLocal(dm, &coordinates);
262: DMGetCoordinateSection(dm, &coordSection);
263: VecGetArray(coordinates, &array);
264: for (v = vStart; v < vEnd; ++v) {
265: const PetscInt idx = v - vStart;
266: PetscInt off, d, val;
268: PetscSectionGetOffset(coordSection, v, &off);
269: for (d = 0; d < dim; ++d) {
270: in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
271: }
272: if (label) {
273: DMLabelGetValue(label, v, &val);
274: in.pointmarkerlist[idx] = val;
275: }
276: }
277: VecRestoreArray(coordinates, &array);
278: }
279: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
280: DMPlexGetGhostCellStratum(dm, &gcStart, NULL);
281: if (gcStart >= 0) cEnd = gcStart;
283: in.numberofcorners = 3;
284: in.numberoftriangles = cEnd - cStart;
286: #if !defined(PETSC_USE_REAL_DOUBLE)
287: 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: 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: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
303: for (v = 0; v < 3; ++v) {
304: in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart;
305: }
306: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
307: }
308: }
309: /* TODO: Segment markers are missing on input */
310: #if 0 /* Do not currently support holes */
311: PetscReal *holeCoords;
312: PetscInt h, d;
314: DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
315: if (in.numberofholes > 0) {
316: PetscMalloc1(in.numberofholes*dim, &in.holelist);
317: for (h = 0; h < in.numberofholes; ++h) {
318: for (d = 0; d < dim; ++d) {
319: in.holelist[h*dim+d] = holeCoords[h*dim+d];
320: }
321: }
322: }
323: #endif
324: if (rank == 0) {
325: char args[32];
327: /* Take away 'Q' for verbose output */
328: PetscStrcpy(args, "pqezQra");
329: triangulate(args, &in, &out, NULL);
330: }
331: PetscFree(in.pointlist);
332: PetscFree(in.pointmarkerlist);
333: PetscFree(in.segmentlist);
334: PetscFree(in.segmentmarkerlist);
335: PetscFree(in.trianglelist);
337: {
338: DMLabel rlabel = NULL;
339: const PetscInt numCorners = 3;
340: const PetscInt numCells = out.numberoftriangles;
341: const PetscInt numVertices = out.numberofpoints;
342: PetscInt *cells;
343: PetscReal *meshCoords;
344: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
346: if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
347: meshCoords = (PetscReal *) out.pointlist;
348: } else {
349: PetscInt i;
351: PetscMalloc1(dim * numVertices,&meshCoords);
352: for (i = 0; i < dim * numVertices; i++) {
353: meshCoords[i] = (PetscReal) out.pointlist[i];
354: }
355: }
356: if (sizeof (PetscInt) == sizeof (out.trianglelist[0])) {
357: cells = (PetscInt *) out.trianglelist;
358: } else {
359: PetscInt i;
361: PetscMalloc1(numCells * numCorners, &cells);
362: for (i = 0; i < numCells * numCorners; i++) {
363: cells[i] = (PetscInt) out.trianglelist[i];
364: }
365: }
367: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
368: if (label) {
369: DMCreateLabel(*dmRefined, labelName);
370: DMGetLabel(*dmRefined, labelName, &rlabel);
371: }
372: if (sizeof (PetscReal) != sizeof (out.pointlist[0])) {
373: PetscFree(meshCoords);
374: }
375: if (sizeof (PetscInt) != sizeof (out.trianglelist[0])) {
376: PetscFree(cells);
377: }
378: /* Set labels */
379: for (v = 0; v < numVertices; ++v) {
380: if (out.pointmarkerlist[v]) {
381: if (rlabel) DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);
382: }
383: }
384: if (interpolate) {
385: PetscInt e;
387: for (e = 0; e < out.numberofedges; e++) {
388: if (out.edgemarkerlist[e]) {
389: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
390: const PetscInt *edges;
391: PetscInt numEdges;
393: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
395: if (rlabel) DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);
396: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
397: }
398: }
399: }
400: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
401: }
402: #if 0 /* Do not currently support holes */
403: DMPlexCopyHoles(*dm, boundary);
404: #endif
405: FiniOutput_Triangle(&out);
406: #if !defined(PETSC_USE_REAL_DOUBLE)
407: PetscFree(maxVolumes);
408: #endif
409: return 0;
410: }