Actual source code: trigenerate.c
petsc-3.11.4 2019-09-28
1: #include <petsc/private/dmpleximpl.h>
3: #include <triangle.h>
5: static PetscErrorCode InitInput_Triangle(struct triangulateio *inputCtx)
6: {
8: inputCtx->numberofpoints = 0;
9: inputCtx->numberofpointattributes = 0;
10: inputCtx->pointlist = NULL;
11: inputCtx->pointattributelist = NULL;
12: inputCtx->pointmarkerlist = NULL;
13: inputCtx->numberofsegments = 0;
14: inputCtx->segmentlist = NULL;
15: inputCtx->segmentmarkerlist = NULL;
16: inputCtx->numberoftriangleattributes = 0;
17: inputCtx->trianglelist = NULL;
18: inputCtx->numberofholes = 0;
19: inputCtx->holelist = NULL;
20: inputCtx->numberofregions = 0;
21: inputCtx->regionlist = NULL;
22: return(0);
23: }
25: static PetscErrorCode InitOutput_Triangle(struct triangulateio *outputCtx)
26: {
28: outputCtx->numberofpoints = 0;
29: outputCtx->pointlist = NULL;
30: outputCtx->pointattributelist = NULL;
31: outputCtx->pointmarkerlist = NULL;
32: outputCtx->numberoftriangles = 0;
33: outputCtx->trianglelist = NULL;
34: outputCtx->triangleattributelist = NULL;
35: outputCtx->neighborlist = NULL;
36: outputCtx->segmentlist = NULL;
37: outputCtx->segmentmarkerlist = NULL;
38: outputCtx->numberofedges = 0;
39: outputCtx->edgelist = NULL;
40: outputCtx->edgemarkerlist = NULL;
41: return(0);
42: }
44: static PetscErrorCode FiniOutput_Triangle(struct triangulateio *outputCtx)
45: {
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;
72: PetscErrorCode ierr;
75: PetscObjectGetComm((PetscObject)boundary,&comm);
76: MPI_Comm_rank(comm, &rank);
77: InitInput_Triangle(&in);
78: InitOutput_Triangle(&out);
79: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
80: DMGetLabel(boundary, labelName, &label);
81: DMGetLabel(boundary, labelName2, &label2);
83: in.numberofpoints = vEnd - vStart;
84: if (in.numberofpoints > 0) {
85: PetscSection coordSection;
86: Vec coordinates;
87: PetscScalar *array;
89: PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
90: PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
91: DMGetCoordinatesLocal(boundary, &coordinates);
92: DMGetCoordinateSection(boundary, &coordSection);
93: VecGetArray(coordinates, &array);
94: for (v = vStart; v < vEnd; ++v) {
95: const PetscInt idx = v - vStart;
96: PetscInt val, off, d;
98: PetscSectionGetOffset(coordSection, v, &off);
99: for (d = 0; d < dim; ++d) {
100: in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
101: }
102: if (label) {
103: DMLabelGetValue(label, v, &val);
104: in.pointmarkerlist[idx] = val;
105: }
106: }
107: VecRestoreArray(coordinates, &array);
108: }
109: DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);
110: in.numberofsegments = eEnd - eStart;
111: if (in.numberofsegments > 0) {
112: PetscMalloc1(in.numberofsegments*2, &in.segmentlist);
113: 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: 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: 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: DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
135: if (in.numberofholes > 0) {
136: 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) {
145: char args[32];
147: /* Take away 'Q' for verbose output */
148: PetscStrcpy(args, "pqezQ");
149: if (createConvexHull) {PetscStrcat(args, "c");}
150: if (constrained) {PetscStrcpy(args, "zepDQ");}
151: if (mesh->triangleOpts) {triangulate(mesh->triangleOpts, &in, &out, NULL);}
152: else {triangulate(args, &in, &out, NULL);}
153: }
154: PetscFree(in.pointlist);
155: PetscFree(in.pointmarkerlist);
156: PetscFree(in.segmentlist);
157: PetscFree(in.segmentmarkerlist);
158: PetscFree(in.holelist);
160: {
161: DMLabel glabel = NULL;
162: DMLabel glabel2 = NULL;
163: const PetscInt numCorners = 3;
164: const PetscInt numCells = out.numberoftriangles;
165: const PetscInt numVertices = out.numberofpoints;
166: const int *cells = out.trianglelist;
167: const double *meshCoords = out.pointlist;
169: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
170: if (label) {DMCreateLabel(*dm, labelName); DMGetLabel(*dm, labelName, &glabel);}
171: if (label2) {DMCreateLabel(*dm, labelName2); DMGetLabel(*dm, labelName2, &glabel2);}
172: /* Set labels */
173: for (v = 0; v < numVertices; ++v) {
174: if (out.pointmarkerlist[v]) {
175: if (glabel) {DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);}
176: }
177: }
178: if (interpolate) {
179: for (e = 0; e < out.numberofedges; e++) {
180: if (out.edgemarkerlist[e]) {
181: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
182: const PetscInt *edges;
183: PetscInt numEdges;
185: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
186: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
187: if (glabel) {DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);}
188: if (glabel2) {DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]);}
189: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
190: }
191: }
192: }
193: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
194: }
195: #if 0 /* Do not currently support holes */
196: DMPlexCopyHoles(*dm, boundary);
197: #endif
198: FiniOutput_Triangle(&out);
199: return(0);
200: }
202: PETSC_EXTERN PetscErrorCode DMPlexRefine_Triangle(DM dm, PetscReal *inmaxVolumes, DM *dmRefined)
203: {
204: MPI_Comm comm;
205: PetscInt dim = 2;
206: const char *labelName = "marker";
207: struct triangulateio in;
208: struct triangulateio out;
209: DMLabel label;
210: PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
211: PetscMPIInt rank;
212: PetscErrorCode ierr;
213: double *maxVolumes;
216: PetscObjectGetComm((PetscObject)dm,&comm);
217: MPI_Comm_rank(comm, &rank);
218: InitInput_Triangle(&in);
219: InitOutput_Triangle(&out);
220: DMPlexGetDepth(dm, &depth);
221: MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
222: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
223: DMGetLabel(dm, labelName, &label);
225: in.numberofpoints = vEnd - vStart;
226: if (in.numberofpoints > 0) {
227: PetscSection coordSection;
228: Vec coordinates;
229: PetscScalar *array;
231: PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
232: PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
233: DMGetCoordinatesLocal(dm, &coordinates);
234: DMGetCoordinateSection(dm, &coordSection);
235: VecGetArray(coordinates, &array);
236: for (v = vStart; v < vEnd; ++v) {
237: const PetscInt idx = v - vStart;
238: PetscInt off, d, val;
240: PetscSectionGetOffset(coordSection, v, &off);
241: for (d = 0; d < dim; ++d) {
242: in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
243: }
244: if (label) {
245: DMLabelGetValue(label, v, &val);
246: in.pointmarkerlist[idx] = val;
247: }
248: }
249: VecRestoreArray(coordinates, &array);
250: }
251: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
253: in.numberofcorners = 3;
254: in.numberoftriangles = cEnd - cStart;
256: #if !defined(PETSC_USE_REAL_DOUBLE)
257: PetscMalloc1(cEnd - cStart,&maxVolumes);
258: for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = (double)inmaxVolumes[c];
259: #else
260: maxVolumes = inmaxVolumes;
261: #endif
263: in.trianglearealist = (double*) maxVolumes;
264: if (in.numberoftriangles > 0) {
265: PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);
266: for (c = cStart; c < cEnd; ++c) {
267: const PetscInt idx = c - cStart;
268: PetscInt *closure = NULL;
269: PetscInt closureSize;
271: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
272: if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize);
273: for (v = 0; v < 3; ++v) {
274: in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart;
275: }
276: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
277: }
278: }
279: /* TODO: Segment markers are missing on input */
280: #if 0 /* Do not currently support holes */
281: PetscReal *holeCoords;
282: PetscInt h, d;
284: DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
285: if (in.numberofholes > 0) {
286: PetscMalloc1(in.numberofholes*dim, &in.holelist);
287: for (h = 0; h < in.numberofholes; ++h) {
288: for (d = 0; d < dim; ++d) {
289: in.holelist[h*dim+d] = holeCoords[h*dim+d];
290: }
291: }
292: }
293: #endif
294: if (!rank) {
295: char args[32];
297: /* Take away 'Q' for verbose output */
298: PetscStrcpy(args, "pqezQra");
299: triangulate(args, &in, &out, NULL);
300: }
301: PetscFree(in.pointlist);
302: PetscFree(in.pointmarkerlist);
303: PetscFree(in.segmentlist);
304: PetscFree(in.segmentmarkerlist);
305: PetscFree(in.trianglelist);
307: {
308: DMLabel rlabel = NULL;
309: const PetscInt numCorners = 3;
310: const PetscInt numCells = out.numberoftriangles;
311: const PetscInt numVertices = out.numberofpoints;
312: const int *cells = out.trianglelist;
313: const double *meshCoords = out.pointlist;
314: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
316: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
317: if (label) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*dmRefined, labelName, &rlabel);}
318: /* Set labels */
319: for (v = 0; v < numVertices; ++v) {
320: if (out.pointmarkerlist[v]) {
321: if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);}
322: }
323: }
324: if (interpolate) {
325: PetscInt e;
327: for (e = 0; e < out.numberofedges; e++) {
328: if (out.edgemarkerlist[e]) {
329: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
330: const PetscInt *edges;
331: PetscInt numEdges;
333: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
334: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
335: if (rlabel) {DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);}
336: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
337: }
338: }
339: }
340: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
341: }
342: #if 0 /* Do not currently support holes */
343: DMPlexCopyHoles(*dm, boundary);
344: #endif
345: FiniOutput_Triangle(&out);
346: #if !defined(PETSC_USE_REAL_DOUBLE)
347: PetscFree(maxVolumes);
348: #endif
349: return(0);
350: }