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