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