Actual source code: trigenerate.c
petsc-3.9.4 2018-09-11
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 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) {DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);}
103: }
104: VecRestoreArray(coordinates, &array);
105: }
106: DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);
107: in.numberofsegments = eEnd - eStart;
108: if (in.numberofsegments > 0) {
109: PetscMalloc1(in.numberofsegments*2, &in.segmentlist);
110: PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist);
111: for (e = eStart; e < eEnd; ++e) {
112: const PetscInt idx = e - eStart;
113: const PetscInt *cone;
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) {DMLabelGetValue(label, e, &in.segmentmarkerlist[idx]);}
121: }
122: }
123: #if 0 /* Do not currently support holes */
124: PetscReal *holeCoords;
125: PetscInt h, d;
127: DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
128: if (in.numberofholes > 0) {
129: PetscMalloc1(in.numberofholes*dim, &in.holelist);
130: for (h = 0; h < in.numberofholes; ++h) {
131: for (d = 0; d < dim; ++d) {
132: in.holelist[h*dim+d] = holeCoords[h*dim+d];
133: }
134: }
135: }
136: #endif
137: if (!rank) {
138: char args[32];
140: /* Take away 'Q' for verbose output */
141: PetscStrcpy(args, "pqezQ");
142: if (createConvexHull) {PetscStrcat(args, "c");}
143: if (constrained) {PetscStrcpy(args, "zepDQ");}
144: if (mesh->triangleOpts) {triangulate(mesh->triangleOpts, &in, &out, NULL);}
145: else {triangulate(args, &in, &out, NULL);}
146: }
147: PetscFree(in.pointlist);
148: PetscFree(in.pointmarkerlist);
149: PetscFree(in.segmentlist);
150: PetscFree(in.segmentmarkerlist);
151: PetscFree(in.holelist);
153: {
154: DMLabel glabel = NULL;
155: DMLabel glabel2 = NULL;
156: const PetscInt numCorners = 3;
157: const PetscInt numCells = out.numberoftriangles;
158: const PetscInt numVertices = out.numberofpoints;
159: const int *cells = out.trianglelist;
160: const double *meshCoords = out.pointlist;
162: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
163: if (label) {DMCreateLabel(*dm, labelName); DMGetLabel(*dm, labelName, &glabel);}
164: if (label2) {DMCreateLabel(*dm, labelName2); DMGetLabel(*dm, labelName2, &glabel2);}
165: /* Set labels */
166: for (v = 0; v < numVertices; ++v) {
167: if (out.pointmarkerlist[v]) {
168: if (glabel) {DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);}
169: }
170: }
171: if (interpolate) {
172: for (e = 0; e < out.numberofedges; e++) {
173: if (out.edgemarkerlist[e]) {
174: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
175: const PetscInt *edges;
176: PetscInt numEdges;
178: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
179: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
180: if (glabel) {DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);}
181: if (glabel2) {DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]);}
182: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
183: }
184: }
185: }
186: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
187: }
188: #if 0 /* Do not currently support holes */
189: DMPlexCopyHoles(*dm, boundary);
190: #endif
191: FiniOutput_Triangle(&out);
192: return(0);
193: }
195: PETSC_EXTERN PetscErrorCode DMPlexRefine_Triangle(DM dm, PetscReal *inmaxVolumes, DM *dmRefined)
196: {
197: MPI_Comm comm;
198: PetscInt dim = 2;
199: const char *labelName = "marker";
200: struct triangulateio in;
201: struct triangulateio out;
202: DMLabel label;
203: PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
204: PetscMPIInt rank;
205: PetscErrorCode ierr;
206: double *maxVolumes;
209: PetscObjectGetComm((PetscObject)dm,&comm);
210: MPI_Comm_rank(comm, &rank);
211: InitInput_Triangle(&in);
212: InitOutput_Triangle(&out);
213: DMPlexGetDepth(dm, &depth);
214: MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
215: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
216: DMGetLabel(dm, labelName, &label);
218: in.numberofpoints = vEnd - vStart;
219: if (in.numberofpoints > 0) {
220: PetscSection coordSection;
221: Vec coordinates;
222: PetscScalar *array;
224: PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
225: PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
226: DMGetCoordinatesLocal(dm, &coordinates);
227: DMGetCoordinateSection(dm, &coordSection);
228: VecGetArray(coordinates, &array);
229: for (v = vStart; v < vEnd; ++v) {
230: const PetscInt idx = v - vStart;
231: PetscInt off, d;
233: PetscSectionGetOffset(coordSection, v, &off);
234: for (d = 0; d < dim; ++d) {
235: in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
236: }
237: if (label) {DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);}
238: }
239: VecRestoreArray(coordinates, &array);
240: }
241: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
243: in.numberofcorners = 3;
244: in.numberoftriangles = cEnd - cStart;
246: #if !defined(PETSC_USE_REAL_DOUBLE)
247: PetscMalloc1(cEnd - cStart,&maxVolumes);
248: for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = (double)inmaxVolumes[c];
249: #else
250: maxVolumes = inmaxVolumes;
251: #endif
253: in.trianglearealist = (double*) maxVolumes;
254: if (in.numberoftriangles > 0) {
255: PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);
256: for (c = cStart; c < cEnd; ++c) {
257: const PetscInt idx = c - cStart;
258: PetscInt *closure = NULL;
259: PetscInt closureSize;
261: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
262: if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize);
263: for (v = 0; v < 3; ++v) {
264: in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart;
265: }
266: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
267: }
268: }
269: /* TODO: Segment markers are missing on input */
270: #if 0 /* Do not currently support holes */
271: PetscReal *holeCoords;
272: PetscInt h, d;
274: DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
275: if (in.numberofholes > 0) {
276: PetscMalloc1(in.numberofholes*dim, &in.holelist);
277: for (h = 0; h < in.numberofholes; ++h) {
278: for (d = 0; d < dim; ++d) {
279: in.holelist[h*dim+d] = holeCoords[h*dim+d];
280: }
281: }
282: }
283: #endif
284: if (!rank) {
285: char args[32];
287: /* Take away 'Q' for verbose output */
288: PetscStrcpy(args, "pqezQra");
289: triangulate(args, &in, &out, NULL);
290: }
291: PetscFree(in.pointlist);
292: PetscFree(in.pointmarkerlist);
293: PetscFree(in.segmentlist);
294: PetscFree(in.segmentmarkerlist);
295: PetscFree(in.trianglelist);
297: {
298: DMLabel rlabel = NULL;
299: const PetscInt numCorners = 3;
300: const PetscInt numCells = out.numberoftriangles;
301: const PetscInt numVertices = out.numberofpoints;
302: const int *cells = out.trianglelist;
303: const double *meshCoords = out.pointlist;
304: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
306: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
307: if (label) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*dmRefined, labelName, &rlabel);}
308: /* Set labels */
309: for (v = 0; v < numVertices; ++v) {
310: if (out.pointmarkerlist[v]) {
311: if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);}
312: }
313: }
314: if (interpolate) {
315: PetscInt e;
317: for (e = 0; e < out.numberofedges; e++) {
318: if (out.edgemarkerlist[e]) {
319: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
320: const PetscInt *edges;
321: PetscInt numEdges;
323: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
324: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
325: if (rlabel) {DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);}
326: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
327: }
328: }
329: }
330: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
331: }
332: #if 0 /* Do not currently support holes */
333: DMPlexCopyHoles(*dm, boundary);
334: #endif
335: FiniOutput_Triangle(&out);
336: #if !defined(PETSC_USE_REAL_DOUBLE)
337: PetscFree(maxVolumes);
338: #endif
339: return(0);
340: }