xref: /petsc/src/dm/impls/plex/generators/triangle/trigenerate.c (revision bfe80ac4a46d58cb7760074b25f5e81b2f541d8a) 
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 
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 
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 
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 
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) {
143         in.holelist[h*dim+d] = holeCoords[h*dim+d];
144       }
145     }
146   }
147 #endif
148   if (rank == 0) {
149     char args[32];
150 
151     /* Take away 'Q' for verbose output */
152     PetscCall(PetscStrncpy(args, "pqezQ", sizeof(args)));
153     if (createConvexHull) PetscCall(PetscStrlcat(args, "c", sizeof(args)));
154     if (constrained) PetscCall(PetscStrncpy(args, "zepDQ", sizeof(args)));
155     if (mesh->triangleOpts) {
156       triangulate(mesh->triangleOpts, &in, &out, NULL);
157     } else {
158       triangulate(args, &in, &out, NULL);
159     }
160   }
161   PetscCall(PetscFree(in.pointlist));
162   PetscCall(PetscFree(in.pointmarkerlist));
163   PetscCall(PetscFree(in.segmentlist));
164   PetscCall(PetscFree(in.segmentmarkerlist));
165   PetscCall(PetscFree(in.holelist));
166 
167   {
168     DMLabel        glabel      = NULL;
169     DMLabel        glabel2     = NULL;
170     const PetscInt numCorners  = 3;
171     const PetscInt numCells    = out.numberoftriangles;
172     const PetscInt numVertices = out.numberofpoints;
173     PetscInt      *cells;
174     PetscReal     *meshCoords;
175 
176     if (sizeof(PetscReal) == sizeof(out.pointlist[0])) {
177       meshCoords = (PetscReal *)out.pointlist;
178     } else {
179       PetscInt i;
180 
181       PetscCall(PetscMalloc1(dim * numVertices, &meshCoords));
182       for (i = 0; i < dim * numVertices; i++) meshCoords[i] = (PetscReal)out.pointlist[i];
183     }
184     if (sizeof(PetscInt) == sizeof(out.trianglelist[0])) {
185       cells = (PetscInt *)out.trianglelist;
186     } else {
187       PetscInt i;
188 
189       PetscCall(PetscMalloc1(numCells * numCorners, &cells));
190       for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out.trianglelist[i];
191     }
192     PetscCall(DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm));
193     if (sizeof(PetscReal) != sizeof(out.pointlist[0])) PetscCall(PetscFree(meshCoords));
194     if (sizeof(PetscInt) != sizeof(out.trianglelist[0])) PetscCall(PetscFree(cells));
195     if (label) {
196       PetscCall(DMCreateLabel(*dm, labelName));
197       PetscCall(DMGetLabel(*dm, labelName, &glabel));
198     }
199     if (label2) {
200       PetscCall(DMCreateLabel(*dm, labelName2));
201       PetscCall(DMGetLabel(*dm, labelName2, &glabel2));
202     }
203     /* Set labels */
204     for (v = 0; v < numVertices; ++v) {
205       if (out.pointmarkerlist[v]) {
206         if (glabel) PetscCall(DMLabelSetValue(glabel, v + numCells, out.pointmarkerlist[v]));
207       }
208     }
209     if (interpolate) {
210       for (e = 0; e < out.numberofedges; e++) {
211         if (out.edgemarkerlist[e]) {
212           const PetscInt  vertices[2] = {out.edgelist[e * 2 + 0] + numCells, out.edgelist[e * 2 + 1] + numCells};
213           const PetscInt *edges;
214           PetscInt        numEdges;
215 
216           PetscCall(DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges));
217           PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %" PetscInt_FMT, numEdges);
218           if (glabel) PetscCall(DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]));
219           if (glabel2) PetscCall(DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]));
220           PetscCall(DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges));
221         }
222       }
223     }
224     PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_FALSE));
225   }
226 #if 0 /* Do not currently support holes */
227   PetscCall(DMPlexCopyHoles(*dm, boundary));
228 #endif
229   PetscCall(FiniOutput_Triangle(&out));
230   PetscFunctionReturn(PETSC_SUCCESS);
231 }
232 
233 PETSC_EXTERN PetscErrorCode DMPlexRefine_Triangle(DM dm, PetscReal *inmaxVolumes, DM *dmRefined)
234 {
235   MPI_Comm             comm;
236   PetscInt             dim       = 2;
237   const char          *labelName = "marker";
238   struct triangulateio in;
239   struct triangulateio out;
240   DMLabel              label;
241   PetscInt             vStart, vEnd, v, gcStart, cStart, cEnd, c, depth, depthGlobal;
242   PetscMPIInt          rank;
243   double              *maxVolumes;
244 
245   PetscFunctionBegin;
246   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
247   PetscCallMPI(MPI_Comm_rank(comm, &rank));
248   PetscCall(InitInput_Triangle(&in));
249   PetscCall(InitOutput_Triangle(&out));
250   PetscCall(DMPlexGetDepth(dm, &depth));
251   PetscCallMPI(MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm));
252   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
253   PetscCall(DMGetLabel(dm, labelName, &label));
254 
255   PetscCall(PetscCIntCast(vEnd - vStart, &in.numberofpoints));
256   if (in.numberofpoints > 0) {
257     PetscSection coordSection;
258     Vec          coordinates;
259     PetscScalar *array;
260 
261     PetscCall(PetscMalloc1(in.numberofpoints * dim, &in.pointlist));
262     PetscCall(PetscMalloc1(in.numberofpoints, &in.pointmarkerlist));
263     PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
264     PetscCall(DMGetCoordinateSection(dm, &coordSection));
265     PetscCall(VecGetArray(coordinates, &array));
266     for (v = vStart; v < vEnd; ++v) {
267       const PetscInt idx = v - vStart;
268       PetscInt       off, d, val;
269 
270       PetscCall(PetscSectionGetOffset(coordSection, v, &off));
271       for (d = 0; d < dim; ++d) in.pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
272       if (label) {
273         PetscCall(DMLabelGetValue(label, v, &val));
274         PetscCall(PetscCIntCast(val, &in.pointmarkerlist[idx]));
275       }
276     }
277     PetscCall(VecRestoreArray(coordinates, &array));
278   }
279   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
280   PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &gcStart, NULL));
281   if (gcStart >= 0) cEnd = gcStart;
282 
283   in.numberofcorners = 3;
284   PetscCall(PetscCIntCast(cEnd - cStart, &in.numberoftriangles));
285 
286 #if !defined(PETSC_USE_REAL_DOUBLE)
287   PetscCall(PetscMalloc1(cEnd - cStart, &maxVolumes));
288   for (c = 0; c < cEnd - cStart; ++c) maxVolumes[c] = (double)inmaxVolumes[c];
289 #else
290   maxVolumes = inmaxVolumes;
291 #endif
292 
293   in.trianglearealist = (double *)maxVolumes;
294   if (in.numberoftriangles > 0) {
295     PetscCall(PetscMalloc1(in.numberoftriangles * in.numberofcorners, &in.trianglelist));
296     for (c = cStart; c < cEnd; ++c) {
297       const PetscInt idx     = c - cStart;
298       PetscInt      *closure = NULL;
299       PetscInt       closureSize;
300 
301       PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
302       PetscCheck(!(closureSize != 4) || !(closureSize != 7), comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %" PetscInt_FMT " vertices in closure", closureSize);
303       for (v = 0; v < 3; ++v) PetscCall(PetscCIntCast(closure[(v + closureSize - 3) * 2] - vStart, &in.trianglelist[idx * in.numberofcorners + v]));
304       PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
305     }
306   }
307   /* TODO: Segment markers are missing on input */
308 #if 0 /* Do not currently support holes */
309   PetscReal *holeCoords;
310   PetscInt   h, d;
311 
312   PetscCall(DMPlexGetHoles(boundary, &in.numberofholes, &holeCords));
313   if (in.numberofholes > 0) {
314     PetscCall(PetscMalloc1(in.numberofholes*dim, &in.holelist));
315     for (h = 0; h < in.numberofholes; ++h) {
316       for (d = 0; d < dim; ++d) {
317         in.holelist[h*dim+d] = holeCoords[h*dim+d];
318       }
319     }
320   }
321 #endif
322   if (rank == 0) {
323     char args[32];
324 
325     /* Take away 'Q' for verbose output */
326     PetscCall(PetscStrncpy(args, "pqezQra", sizeof(args)));
327     triangulate(args, &in, &out, NULL);
328   }
329   PetscCall(PetscFree(in.pointlist));
330   PetscCall(PetscFree(in.pointmarkerlist));
331   PetscCall(PetscFree(in.segmentlist));
332   PetscCall(PetscFree(in.segmentmarkerlist));
333   PetscCall(PetscFree(in.trianglelist));
334 
335   {
336     DMLabel        rlabel      = NULL;
337     const PetscInt numCorners  = 3;
338     const PetscInt numCells    = out.numberoftriangles;
339     const PetscInt numVertices = out.numberofpoints;
340     PetscInt      *cells;
341     PetscReal     *meshCoords;
342     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
343 
344     if (sizeof(PetscReal) == sizeof(out.pointlist[0])) {
345       meshCoords = (PetscReal *)out.pointlist;
346     } else {
347       PetscInt i;
348 
349       PetscCall(PetscMalloc1(dim * numVertices, &meshCoords));
350       for (i = 0; i < dim * numVertices; i++) meshCoords[i] = (PetscReal)out.pointlist[i];
351     }
352     if (sizeof(PetscInt) == sizeof(out.trianglelist[0])) {
353       cells = (PetscInt *)out.trianglelist;
354     } else {
355       PetscInt i;
356 
357       PetscCall(PetscMalloc1(numCells * numCorners, &cells));
358       for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out.trianglelist[i];
359     }
360 
361     PetscCall(DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined));
362     if (label) {
363       PetscCall(DMCreateLabel(*dmRefined, labelName));
364       PetscCall(DMGetLabel(*dmRefined, labelName, &rlabel));
365     }
366     if (sizeof(PetscReal) != sizeof(out.pointlist[0])) PetscCall(PetscFree(meshCoords));
367     if (sizeof(PetscInt) != sizeof(out.trianglelist[0])) PetscCall(PetscFree(cells));
368     /* Set labels */
369     for (v = 0; v < numVertices; ++v) {
370       if (out.pointmarkerlist[v]) {
371         if (rlabel) PetscCall(DMLabelSetValue(rlabel, v + numCells, out.pointmarkerlist[v]));
372       }
373     }
374     if (interpolate) {
375       PetscInt e;
376 
377       for (e = 0; e < out.numberofedges; e++) {
378         if (out.edgemarkerlist[e]) {
379           const PetscInt  vertices[2] = {out.edgelist[e * 2 + 0] + numCells, out.edgelist[e * 2 + 1] + numCells};
380           const PetscInt *edges;
381           PetscInt        numEdges;
382 
383           PetscCall(DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges));
384           PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %" PetscInt_FMT, numEdges);
385           if (rlabel) PetscCall(DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]));
386           PetscCall(DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges));
387         }
388       }
389     }
390     PetscCall(DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE));
391   }
392 #if 0 /* Do not currently support holes */
393   PetscCall(DMPlexCopyHoles(*dm, boundary));
394 #endif
395   PetscCall(FiniOutput_Triangle(&out));
396 #if !defined(PETSC_USE_REAL_DOUBLE)
397   PetscCall(PetscFree(maxVolumes));
398 #endif
399   PetscFunctionReturn(PETSC_SUCCESS);
400 }
401