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