xref: /petsc/src/dm/impls/plex/generators/triangle/trigenerate.c (revision 534a8f05a7a8aff70dd8cfd53d9cd834400a8dbf) 
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(0);
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(0);
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(0);
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   PetscErrorCode       ierr;
76 
77   PetscFunctionBegin;
78   ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr);
79   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
80   ierr = InitInput_Triangle(&in);CHKERRQ(ierr);
81   ierr = InitOutput_Triangle(&out);CHKERRQ(ierr);
82   ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr);
83   ierr = DMGetLabel(boundary, labelName,  &label);CHKERRQ(ierr);
84   ierr = DMGetLabel(boundary, labelName2, &label2);CHKERRQ(ierr);
85 
86   in.numberofpoints = vEnd - vStart;
87   if (in.numberofpoints > 0) {
88     PetscSection coordSection;
89     Vec          coordinates;
90     PetscScalar *array;
91 
92     ierr = PetscMalloc1(in.numberofpoints*dim, &in.pointlist);CHKERRQ(ierr);
93     ierr = PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);CHKERRQ(ierr);
94     ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr);
95     ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr);
96     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
97     for (v = vStart; v < vEnd; ++v) {
98       const PetscInt idx = v - vStart;
99       PetscInt       val, off, d;
100 
101       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
102       for (d = 0; d < dim; ++d) {
103         in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
104       }
105       if (label) {
106         ierr = DMLabelGetValue(label, v, &val);CHKERRQ(ierr);
107         in.pointmarkerlist[idx] = val;
108       }
109     }
110     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
111   }
112   ierr  = DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);CHKERRQ(ierr);
113   in.numberofsegments = eEnd - eStart;
114   if (in.numberofsegments > 0) {
115     ierr = PetscMalloc1(in.numberofsegments*2, &in.segmentlist);CHKERRQ(ierr);
116     ierr = PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist);CHKERRQ(ierr);
117     for (e = eStart; e < eEnd; ++e) {
118       const PetscInt  idx = e - eStart;
119       const PetscInt *cone;
120       PetscInt        val;
121 
122       ierr = DMPlexGetCone(boundary, e, &cone);CHKERRQ(ierr);
123 
124       in.segmentlist[idx*2+0] = cone[0] - vStart;
125       in.segmentlist[idx*2+1] = cone[1] - vStart;
126 
127       if (label) {
128         ierr = DMLabelGetValue(label, e, &val);CHKERRQ(ierr);
129         in.segmentmarkerlist[idx] = val;
130       }
131     }
132   }
133 #if 0 /* Do not currently support holes */
134   PetscReal *holeCoords;
135   PetscInt   h, d;
136 
137   ierr = DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);CHKERRQ(ierr);
138   if (in.numberofholes > 0) {
139     ierr = PetscMalloc1(in.numberofholes*dim, &in.holelist);CHKERRQ(ierr);
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];
149 
150     /* Take away 'Q' for verbose output */
151     ierr = PetscStrcpy(args, "pqezQ");CHKERRQ(ierr);
152     if (createConvexHull)   {ierr = PetscStrcat(args, "c");CHKERRQ(ierr);}
153     if (constrained)        {ierr = PetscStrcpy(args, "zepDQ");CHKERRQ(ierr);}
154     if (mesh->triangleOpts) {triangulate(mesh->triangleOpts, &in, &out, NULL);}
155     else                    {triangulate(args, &in, &out, NULL);}
156   }
157   ierr = PetscFree(in.pointlist);CHKERRQ(ierr);
158   ierr = PetscFree(in.pointmarkerlist);CHKERRQ(ierr);
159   ierr = PetscFree(in.segmentlist);CHKERRQ(ierr);
160   ierr = PetscFree(in.segmentmarkerlist);CHKERRQ(ierr);
161   ierr = PetscFree(in.holelist);CHKERRQ(ierr);
162 
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;
171 
172     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr);
173     if (label)  {ierr = DMCreateLabel(*dm, labelName);  ierr = DMGetLabel(*dm, labelName,  &glabel);}
174     if (label2) {ierr = DMCreateLabel(*dm, labelName2); ierr = DMGetLabel(*dm, labelName2, &glabel2);}
175     /* Set labels */
176     for (v = 0; v < numVertices; ++v) {
177       if (out.pointmarkerlist[v]) {
178         if (glabel) {ierr = DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);}
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;
187 
188           ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
189           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
190           if (glabel)  {ierr = DMLabelSetValue(glabel,  edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);}
191           if (glabel2) {ierr = DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);}
192           ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
193         }
194       }
195     }
196     ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
197   }
198 #if 0 /* Do not currently support holes */
199   ierr = DMPlexCopyHoles(*dm, boundary);CHKERRQ(ierr);
200 #endif
201   ierr = FiniOutput_Triangle(&out);CHKERRQ(ierr);
202   PetscFunctionReturn(0);
203 }
204 
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;
217 
218   PetscFunctionBegin;
219   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
220   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
221   ierr = InitInput_Triangle(&in);CHKERRQ(ierr);
222   ierr = InitOutput_Triangle(&out);CHKERRQ(ierr);
223   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
224   ierr = MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
225   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
226   ierr = DMGetLabel(dm, labelName, &label);CHKERRQ(ierr);
227 
228   in.numberofpoints = vEnd - vStart;
229   if (in.numberofpoints > 0) {
230     PetscSection coordSection;
231     Vec          coordinates;
232     PetscScalar *array;
233 
234     ierr = PetscMalloc1(in.numberofpoints*dim, &in.pointlist);CHKERRQ(ierr);
235     ierr = PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);CHKERRQ(ierr);
236     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
237     ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
238     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
239     for (v = vStart; v < vEnd; ++v) {
240       const PetscInt idx = v - vStart;
241       PetscInt       off, d, val;
242 
243       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
244       for (d = 0; d < dim; ++d) {
245         in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
246       }
247       if (label) {
248         ierr = DMLabelGetValue(label, v, &val);CHKERRQ(ierr);
249         in.pointmarkerlist[idx] = val;
250       }
251     }
252     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
253   }
254   ierr  = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
255 
256   in.numberofcorners   = 3;
257   in.numberoftriangles = cEnd - cStart;
258 
259 #if !defined(PETSC_USE_REAL_DOUBLE)
260   ierr = PetscMalloc1(cEnd - cStart,&maxVolumes);CHKERRQ(ierr);
261   for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = (double)inmaxVolumes[c];
262 #else
263   maxVolumes = inmaxVolumes;
264 #endif
265 
266   in.trianglearealist  = (double*) maxVolumes;
267   if (in.numberoftriangles > 0) {
268     ierr = PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);CHKERRQ(ierr);
269     for (c = cStart; c < cEnd; ++c) {
270       const PetscInt idx      = c - cStart;
271       PetscInt      *closure = NULL;
272       PetscInt       closureSize;
273 
274       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
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       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
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;
286 
287   ierr = DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);CHKERRQ(ierr);
288   if (in.numberofholes > 0) {
289     ierr = PetscMalloc1(in.numberofholes*dim, &in.holelist);CHKERRQ(ierr);
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];
299 
300     /* Take away 'Q' for verbose output */
301     ierr = PetscStrcpy(args, "pqezQra");CHKERRQ(ierr);
302     triangulate(args, &in, &out, NULL);
303   }
304   ierr = PetscFree(in.pointlist);CHKERRQ(ierr);
305   ierr = PetscFree(in.pointmarkerlist);CHKERRQ(ierr);
306   ierr = PetscFree(in.segmentlist);CHKERRQ(ierr);
307   ierr = PetscFree(in.segmentmarkerlist);CHKERRQ(ierr);
308   ierr = PetscFree(in.trianglelist);CHKERRQ(ierr);
309 
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;
318 
319     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr);
320     if (label) {ierr = DMCreateLabel(*dmRefined, labelName); ierr = DMGetLabel(*dmRefined, labelName, &rlabel);}
321     /* Set labels */
322     for (v = 0; v < numVertices; ++v) {
323       if (out.pointmarkerlist[v]) {
324         if (rlabel) {ierr = DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);}
325       }
326     }
327     if (interpolate) {
328       PetscInt e;
329 
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;
335 
336           ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
337           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
338           if (rlabel) {ierr = DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);}
339           ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
340         }
341       }
342     }
343     ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr);
344   }
345 #if 0 /* Do not currently support holes */
346   ierr = DMPlexCopyHoles(*dm, boundary);CHKERRQ(ierr);
347 #endif
348   ierr = FiniOutput_Triangle(&out);CHKERRQ(ierr);
349 #if !defined(PETSC_USE_REAL_DOUBLE)
350   ierr = PetscFree(maxVolumes);CHKERRQ(ierr);
351 #endif
352   PetscFunctionReturn(0);
353 }
354