Actual source code: trigenerate.c
petsc-3.13.6 2020-09-29
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: const int *cells = out.trianglelist;
170: const double *meshCoords = out.pointlist;
172: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
173: if (label) {DMCreateLabel(*dm, labelName); DMGetLabel(*dm, labelName, &glabel);}
174: if (label2) {DMCreateLabel(*dm, labelName2); DMGetLabel(*dm, labelName2, &glabel2);}
175: /* Set labels */
176: for (v = 0; v < numVertices; ++v) {
177: if (out.pointmarkerlist[v]) {
178: if (glabel) {DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);}
179: }
180: }
181: if (interpolate) {
182: for (e = 0; e < out.numberofedges; e++) {
183: if (out.edgemarkerlist[e]) {
184: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
185: const PetscInt *edges;
186: PetscInt numEdges;
188: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
189: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
190: if (glabel) {DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);}
191: if (glabel2) {DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]);}
192: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
193: }
194: }
195: }
196: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
197: }
198: #if 0 /* Do not currently support holes */
199: DMPlexCopyHoles(*dm, boundary);
200: #endif
201: FiniOutput_Triangle(&out);
202: return(0);
203: }
205: PETSC_EXTERN PetscErrorCode DMPlexRefine_Triangle(DM dm, PetscReal *inmaxVolumes, DM *dmRefined)
206: {
207: MPI_Comm comm;
208: PetscInt dim = 2;
209: const char *labelName = "marker";
210: struct triangulateio in;
211: struct triangulateio out;
212: DMLabel label;
213: PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
214: PetscMPIInt rank;
215: PetscErrorCode ierr;
216: double *maxVolumes;
219: PetscObjectGetComm((PetscObject)dm,&comm);
220: MPI_Comm_rank(comm, &rank);
221: InitInput_Triangle(&in);
222: InitOutput_Triangle(&out);
223: DMPlexGetDepth(dm, &depth);
224: MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
225: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
226: DMGetLabel(dm, labelName, &label);
228: in.numberofpoints = vEnd - vStart;
229: if (in.numberofpoints > 0) {
230: PetscSection coordSection;
231: Vec coordinates;
232: PetscScalar *array;
234: PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
235: PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
236: DMGetCoordinatesLocal(dm, &coordinates);
237: DMGetCoordinateSection(dm, &coordSection);
238: VecGetArray(coordinates, &array);
239: for (v = vStart; v < vEnd; ++v) {
240: const PetscInt idx = v - vStart;
241: PetscInt off, d, val;
243: PetscSectionGetOffset(coordSection, v, &off);
244: for (d = 0; d < dim; ++d) {
245: in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
246: }
247: if (label) {
248: DMLabelGetValue(label, v, &val);
249: in.pointmarkerlist[idx] = val;
250: }
251: }
252: VecRestoreArray(coordinates, &array);
253: }
254: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
256: in.numberofcorners = 3;
257: in.numberoftriangles = cEnd - cStart;
259: #if !defined(PETSC_USE_REAL_DOUBLE)
260: PetscMalloc1(cEnd - cStart,&maxVolumes);
261: for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = (double)inmaxVolumes[c];
262: #else
263: maxVolumes = inmaxVolumes;
264: #endif
266: in.trianglearealist = (double*) maxVolumes;
267: if (in.numberoftriangles > 0) {
268: PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);
269: for (c = cStart; c < cEnd; ++c) {
270: const PetscInt idx = c - cStart;
271: PetscInt *closure = NULL;
272: PetscInt closureSize;
274: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
275: if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize);
276: for (v = 0; v < 3; ++v) {
277: in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart;
278: }
279: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
280: }
281: }
282: /* TODO: Segment markers are missing on input */
283: #if 0 /* Do not currently support holes */
284: PetscReal *holeCoords;
285: PetscInt h, d;
287: DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
288: if (in.numberofholes > 0) {
289: PetscMalloc1(in.numberofholes*dim, &in.holelist);
290: for (h = 0; h < in.numberofholes; ++h) {
291: for (d = 0; d < dim; ++d) {
292: in.holelist[h*dim+d] = holeCoords[h*dim+d];
293: }
294: }
295: }
296: #endif
297: if (!rank) {
298: char args[32];
300: /* Take away 'Q' for verbose output */
301: PetscStrcpy(args, "pqezQra");
302: triangulate(args, &in, &out, NULL);
303: }
304: PetscFree(in.pointlist);
305: PetscFree(in.pointmarkerlist);
306: PetscFree(in.segmentlist);
307: PetscFree(in.segmentmarkerlist);
308: PetscFree(in.trianglelist);
310: {
311: DMLabel rlabel = NULL;
312: const PetscInt numCorners = 3;
313: const PetscInt numCells = out.numberoftriangles;
314: const PetscInt numVertices = out.numberofpoints;
315: const int *cells = out.trianglelist;
316: const double *meshCoords = out.pointlist;
317: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
319: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
320: if (label) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*dmRefined, labelName, &rlabel);}
321: /* Set labels */
322: for (v = 0; v < numVertices; ++v) {
323: if (out.pointmarkerlist[v]) {
324: if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);}
325: }
326: }
327: if (interpolate) {
328: PetscInt e;
330: for (e = 0; e < out.numberofedges; e++) {
331: if (out.edgemarkerlist[e]) {
332: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
333: const PetscInt *edges;
334: PetscInt numEdges;
336: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
337: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
338: if (rlabel) {DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);}
339: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
340: }
341: }
342: }
343: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
344: }
345: #if 0 /* Do not currently support holes */
346: DMPlexCopyHoles(*dm, boundary);
347: #endif
348: FiniOutput_Triangle(&out);
349: #if !defined(PETSC_USE_REAL_DOUBLE)
350: PetscFree(maxVolumes);
351: #endif
352: return(0);
353: }