xref: /petsc/src/dm/impls/plex/plexgenerate.c (revision 276c5506312304d9c58828fa177e3f5c76c17b23) 
1 #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 
3 PetscErrorCode DMPlexInvertCell_Internal(PetscInt dim, PetscInt numCorners, PetscInt cone[])
4 {
5   int tmpc;
6 
7   PetscFunctionBegin;
8   if (dim != 3) PetscFunctionReturn(0);
9   switch (numCorners) {
10   case 4:
11     tmpc    = cone[0];
12     cone[0] = cone[1];
13     cone[1] = tmpc;
14     break;
15   case 8:
16     tmpc    = cone[1];
17     cone[1] = cone[3];
18     cone[3] = tmpc;
19     break;
20   default: break;
21   }
22   PetscFunctionReturn(0);
23 }
24 
25 /*@C
26   DMPlexInvertCell - This flips tetrahedron and hexahedron orientation since Plex stores them internally with outward normals. Other cells are left untouched.
27 
28   Input Parameters:
29 + numCorners - The number of vertices in a cell
30 - cone - The incoming cone
31 
32   Output Parameter:
33 . cone - The inverted cone (in-place)
34 
35   Level: developer
36 
37 .seealso: DMPlexGenerate()
38 @*/
39 PetscErrorCode DMPlexInvertCell(PetscInt dim, PetscInt numCorners, int cone[])
40 {
41   int tmpc;
42 
43   PetscFunctionBegin;
44   if (dim != 3) PetscFunctionReturn(0);
45   switch (numCorners) {
46   case 4:
47     tmpc    = cone[0];
48     cone[0] = cone[1];
49     cone[1] = tmpc;
50     break;
51   case 8:
52     tmpc    = cone[1];
53     cone[1] = cone[3];
54     cone[3] = tmpc;
55     break;
56   default: break;
57   }
58   PetscFunctionReturn(0);
59 }
60 
61 /* This is to fix the tetrahedron orientation from TetGen */
62 PETSC_UNUSED static PetscErrorCode DMPlexInvertCells_Internal(PetscInt dim, PetscInt numCells, PetscInt numCorners, int cells[])
63 {
64   PetscInt       bound = numCells*numCorners, coff;
65   PetscErrorCode ierr;
66 
67   PetscFunctionBegin;
68   for (coff = 0; coff < bound; coff += numCorners) {
69     ierr = DMPlexInvertCell(dim, numCorners, &cells[coff]);CHKERRQ(ierr);
70   }
71   PetscFunctionReturn(0);
72 }
73 
74 /*@C
75   DMPlexTriangleSetOptions - Set the options used for the Triangle mesh generator
76 
77   Not Collective
78 
79   Inputs Parameters:
80 + dm - The DMPlex object
81 - opts - The command line options
82 
83   Level: developer
84 
85 .keywords: mesh, points
86 .seealso: DMPlexTetgenSetOptions(), DMPlexGenerate()
87 @*/
88 PetscErrorCode DMPlexTriangleSetOptions(DM dm, const char *opts)
89 {
90   DM_Plex       *mesh = (DM_Plex*) dm->data;
91   PetscErrorCode ierr;
92 
93   PetscFunctionBegin;
94   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
95   PetscValidPointer(opts, 2);
96   ierr = PetscFree(mesh->triangleOpts);CHKERRQ(ierr);
97   ierr = PetscStrallocpy(opts, &mesh->triangleOpts);CHKERRQ(ierr);
98   PetscFunctionReturn(0);
99 }
100 
101 /*@C
102   DMPlexTetgenSetOptions - Set the options used for the Tetgen mesh generator
103 
104   Not Collective
105 
106   Inputs Parameters:
107 + dm - The DMPlex object
108 - opts - The command line options
109 
110   Level: developer
111 
112 .keywords: mesh, points
113 .seealso: DMPlexTriangleSetOptions(), DMPlexGenerate()
114 @*/
115 PetscErrorCode DMPlexTetgenSetOptions(DM dm, const char *opts)
116 {
117   DM_Plex       *mesh = (DM_Plex*) dm->data;
118   PetscErrorCode ierr;
119 
120   PetscFunctionBegin;
121   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
122   PetscValidPointer(opts, 2);
123   ierr = PetscFree(mesh->tetgenOpts);CHKERRQ(ierr);
124   ierr = PetscStrallocpy(opts, &mesh->tetgenOpts);CHKERRQ(ierr);
125   PetscFunctionReturn(0);
126 }
127 
128 #if defined(PETSC_HAVE_TRIANGLE)
129 #include <triangle.h>
130 
131 static PetscErrorCode InitInput_Triangle(struct triangulateio *inputCtx)
132 {
133   PetscFunctionBegin;
134   inputCtx->numberofpoints             = 0;
135   inputCtx->numberofpointattributes    = 0;
136   inputCtx->pointlist                  = NULL;
137   inputCtx->pointattributelist         = NULL;
138   inputCtx->pointmarkerlist            = NULL;
139   inputCtx->numberofsegments           = 0;
140   inputCtx->segmentlist                = NULL;
141   inputCtx->segmentmarkerlist          = NULL;
142   inputCtx->numberoftriangleattributes = 0;
143   inputCtx->trianglelist               = NULL;
144   inputCtx->numberofholes              = 0;
145   inputCtx->holelist                   = NULL;
146   inputCtx->numberofregions            = 0;
147   inputCtx->regionlist                 = NULL;
148   PetscFunctionReturn(0);
149 }
150 
151 static PetscErrorCode InitOutput_Triangle(struct triangulateio *outputCtx)
152 {
153   PetscFunctionBegin;
154   outputCtx->numberofpoints        = 0;
155   outputCtx->pointlist             = NULL;
156   outputCtx->pointattributelist    = NULL;
157   outputCtx->pointmarkerlist       = NULL;
158   outputCtx->numberoftriangles     = 0;
159   outputCtx->trianglelist          = NULL;
160   outputCtx->triangleattributelist = NULL;
161   outputCtx->neighborlist          = NULL;
162   outputCtx->segmentlist           = NULL;
163   outputCtx->segmentmarkerlist     = NULL;
164   outputCtx->numberofedges         = 0;
165   outputCtx->edgelist              = NULL;
166   outputCtx->edgemarkerlist        = NULL;
167   PetscFunctionReturn(0);
168 }
169 
170 static PetscErrorCode FiniOutput_Triangle(struct triangulateio *outputCtx)
171 {
172   PetscFunctionBegin;
173   free(outputCtx->pointlist);
174   free(outputCtx->pointmarkerlist);
175   free(outputCtx->segmentlist);
176   free(outputCtx->segmentmarkerlist);
177   free(outputCtx->edgelist);
178   free(outputCtx->edgemarkerlist);
179   free(outputCtx->trianglelist);
180   free(outputCtx->neighborlist);
181   PetscFunctionReturn(0);
182 }
183 
184 PetscErrorCode DMPlexGenerate_Triangle(DM boundary, PetscBool interpolate, DM *dm)
185 {
186   MPI_Comm             comm;
187   DM_Plex             *mesh             = (DM_Plex *) boundary->data;
188   PetscInt             dim              = 2;
189   const PetscBool      createConvexHull = PETSC_FALSE;
190   const PetscBool      constrained      = PETSC_FALSE;
191   const char          *labelName        = "marker";
192   const char          *labelName2       = "Face Sets";
193   struct triangulateio in;
194   struct triangulateio out;
195   DMLabel              label, label2;
196   PetscInt             vStart, vEnd, v, eStart, eEnd, e;
197   PetscMPIInt          rank;
198   PetscErrorCode       ierr;
199 
200   PetscFunctionBegin;
201   ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr);
202   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
203   ierr = InitInput_Triangle(&in);CHKERRQ(ierr);
204   ierr = InitOutput_Triangle(&out);CHKERRQ(ierr);
205   ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr);
206   ierr = DMGetLabel(boundary, labelName,  &label);CHKERRQ(ierr);
207   ierr = DMGetLabel(boundary, labelName2, &label2);CHKERRQ(ierr);
208 
209   in.numberofpoints = vEnd - vStart;
210   if (in.numberofpoints > 0) {
211     PetscSection coordSection;
212     Vec          coordinates;
213     PetscScalar *array;
214 
215     ierr = PetscMalloc1(in.numberofpoints*dim, &in.pointlist);CHKERRQ(ierr);
216     ierr = PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);CHKERRQ(ierr);
217     ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr);
218     ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr);
219     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
220     for (v = vStart; v < vEnd; ++v) {
221       const PetscInt idx = v - vStart;
222       PetscInt       off, d;
223 
224       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
225       for (d = 0; d < dim; ++d) {
226         in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
227       }
228       if (label) {ierr = DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);CHKERRQ(ierr);}
229     }
230     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
231   }
232   ierr  = DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);CHKERRQ(ierr);
233   in.numberofsegments = eEnd - eStart;
234   if (in.numberofsegments > 0) {
235     ierr = PetscMalloc1(in.numberofsegments*2, &in.segmentlist);CHKERRQ(ierr);
236     ierr = PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist);CHKERRQ(ierr);
237     for (e = eStart; e < eEnd; ++e) {
238       const PetscInt  idx = e - eStart;
239       const PetscInt *cone;
240 
241       ierr = DMPlexGetCone(boundary, e, &cone);CHKERRQ(ierr);
242 
243       in.segmentlist[idx*2+0] = cone[0] - vStart;
244       in.segmentlist[idx*2+1] = cone[1] - vStart;
245 
246       if (label) {ierr = DMLabelGetValue(label, e, &in.segmentmarkerlist[idx]);CHKERRQ(ierr);}
247     }
248   }
249 #if 0 /* Do not currently support holes */
250   PetscReal *holeCoords;
251   PetscInt   h, d;
252 
253   ierr = DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);CHKERRQ(ierr);
254   if (in.numberofholes > 0) {
255     ierr = PetscMalloc1(in.numberofholes*dim, &in.holelist);CHKERRQ(ierr);
256     for (h = 0; h < in.numberofholes; ++h) {
257       for (d = 0; d < dim; ++d) {
258         in.holelist[h*dim+d] = holeCoords[h*dim+d];
259       }
260     }
261   }
262 #endif
263   if (!rank) {
264     char args[32];
265 
266     /* Take away 'Q' for verbose output */
267     ierr = PetscStrcpy(args, "pqezQ");CHKERRQ(ierr);
268     if (createConvexHull)   {ierr = PetscStrcat(args, "c");CHKERRQ(ierr);}
269     if (constrained)        {ierr = PetscStrcpy(args, "zepDQ");CHKERRQ(ierr);}
270     if (mesh->triangleOpts) {triangulate(mesh->triangleOpts, &in, &out, NULL);}
271     else                    {triangulate(args, &in, &out, NULL);}
272   }
273   ierr = PetscFree(in.pointlist);CHKERRQ(ierr);
274   ierr = PetscFree(in.pointmarkerlist);CHKERRQ(ierr);
275   ierr = PetscFree(in.segmentlist);CHKERRQ(ierr);
276   ierr = PetscFree(in.segmentmarkerlist);CHKERRQ(ierr);
277   ierr = PetscFree(in.holelist);CHKERRQ(ierr);
278 
279   {
280     DMLabel        glabel      = NULL;
281     DMLabel        glabel2     = NULL;
282     const PetscInt numCorners  = 3;
283     const PetscInt numCells    = out.numberoftriangles;
284     const PetscInt numVertices = out.numberofpoints;
285     const int     *cells      = out.trianglelist;
286     const double  *meshCoords = out.pointlist;
287 
288     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr);
289     if (label)  {ierr = DMCreateLabel(*dm, labelName);  ierr = DMGetLabel(*dm, labelName,  &glabel);}
290     if (label2) {ierr = DMCreateLabel(*dm, labelName2); ierr = DMGetLabel(*dm, labelName2, &glabel2);}
291     /* Set labels */
292     for (v = 0; v < numVertices; ++v) {
293       if (out.pointmarkerlist[v]) {
294         if (glabel) {ierr = DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);}
295       }
296     }
297     if (interpolate) {
298       for (e = 0; e < out.numberofedges; e++) {
299         if (out.edgemarkerlist[e]) {
300           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
301           const PetscInt *edges;
302           PetscInt        numEdges;
303 
304           ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
305           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
306           if (glabel)  {ierr = DMLabelSetValue(glabel,  edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);}
307           if (glabel2) {ierr = DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);}
308           ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
309         }
310       }
311     }
312     ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
313   }
314 #if 0 /* Do not currently support holes */
315   ierr = DMPlexCopyHoles(*dm, boundary);CHKERRQ(ierr);
316 #endif
317   ierr = FiniOutput_Triangle(&out);CHKERRQ(ierr);
318   PetscFunctionReturn(0);
319 }
320 
321 PetscErrorCode DMPlexRefine_Triangle(DM dm, double *maxVolumes, DM *dmRefined)
322 {
323   MPI_Comm             comm;
324   PetscInt             dim       = 2;
325   const char          *labelName = "marker";
326   struct triangulateio in;
327   struct triangulateio out;
328   DMLabel              label;
329   PetscInt             vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
330   PetscMPIInt          rank;
331   PetscErrorCode       ierr;
332 
333   PetscFunctionBegin;
334   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
335   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
336   ierr = InitInput_Triangle(&in);CHKERRQ(ierr);
337   ierr = InitOutput_Triangle(&out);CHKERRQ(ierr);
338   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
339   ierr = MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
340   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
341   ierr = DMGetLabel(dm, labelName, &label);CHKERRQ(ierr);
342 
343   in.numberofpoints = vEnd - vStart;
344   if (in.numberofpoints > 0) {
345     PetscSection coordSection;
346     Vec          coordinates;
347     PetscScalar *array;
348 
349     ierr = PetscMalloc1(in.numberofpoints*dim, &in.pointlist);CHKERRQ(ierr);
350     ierr = PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);CHKERRQ(ierr);
351     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
352     ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
353     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
354     for (v = vStart; v < vEnd; ++v) {
355       const PetscInt idx = v - vStart;
356       PetscInt       off, d;
357 
358       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
359       for (d = 0; d < dim; ++d) {
360         in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
361       }
362       if (label) {ierr = DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);CHKERRQ(ierr);}
363     }
364     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
365   }
366   ierr  = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
367 
368   in.numberofcorners   = 3;
369   in.numberoftriangles = cEnd - cStart;
370 
371   in.trianglearealist  = (double*) maxVolumes;
372   if (in.numberoftriangles > 0) {
373     ierr = PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);CHKERRQ(ierr);
374     for (c = cStart; c < cEnd; ++c) {
375       const PetscInt idx      = c - cStart;
376       PetscInt      *closure = NULL;
377       PetscInt       closureSize;
378 
379       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
380       if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize);
381       for (v = 0; v < 3; ++v) {
382         in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart;
383       }
384       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
385     }
386   }
387   /* TODO: Segment markers are missing on input */
388 #if 0 /* Do not currently support holes */
389   PetscReal *holeCoords;
390   PetscInt   h, d;
391 
392   ierr = DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);CHKERRQ(ierr);
393   if (in.numberofholes > 0) {
394     ierr = PetscMalloc1(in.numberofholes*dim, &in.holelist);CHKERRQ(ierr);
395     for (h = 0; h < in.numberofholes; ++h) {
396       for (d = 0; d < dim; ++d) {
397         in.holelist[h*dim+d] = holeCoords[h*dim+d];
398       }
399     }
400   }
401 #endif
402   if (!rank) {
403     char args[32];
404 
405     /* Take away 'Q' for verbose output */
406     ierr = PetscStrcpy(args, "pqezQra");CHKERRQ(ierr);
407     triangulate(args, &in, &out, NULL);
408   }
409   ierr = PetscFree(in.pointlist);CHKERRQ(ierr);
410   ierr = PetscFree(in.pointmarkerlist);CHKERRQ(ierr);
411   ierr = PetscFree(in.segmentlist);CHKERRQ(ierr);
412   ierr = PetscFree(in.segmentmarkerlist);CHKERRQ(ierr);
413   ierr = PetscFree(in.trianglelist);CHKERRQ(ierr);
414 
415   {
416     DMLabel        rlabel      = NULL;
417     const PetscInt numCorners  = 3;
418     const PetscInt numCells    = out.numberoftriangles;
419     const PetscInt numVertices = out.numberofpoints;
420     const int     *cells      = out.trianglelist;
421     const double  *meshCoords = out.pointlist;
422     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
423 
424     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr);
425     if (label) {ierr = DMCreateLabel(*dmRefined, labelName); ierr = DMGetLabel(*dmRefined, labelName, &rlabel);}
426     /* Set labels */
427     for (v = 0; v < numVertices; ++v) {
428       if (out.pointmarkerlist[v]) {
429         if (rlabel) {ierr = DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);}
430       }
431     }
432     if (interpolate) {
433       PetscInt e;
434 
435       for (e = 0; e < out.numberofedges; e++) {
436         if (out.edgemarkerlist[e]) {
437           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
438           const PetscInt *edges;
439           PetscInt        numEdges;
440 
441           ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
442           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
443           if (rlabel) {ierr = DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);}
444           ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
445         }
446       }
447     }
448     ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr);
449   }
450 #if 0 /* Do not currently support holes */
451   ierr = DMPlexCopyHoles(*dm, boundary);CHKERRQ(ierr);
452 #endif
453   ierr = FiniOutput_Triangle(&out);CHKERRQ(ierr);
454   PetscFunctionReturn(0);
455 }
456 #endif /* PETSC_HAVE_TRIANGLE */
457 
458 #if defined(PETSC_HAVE_TETGEN)
459 #include <tetgen.h>
460 PetscErrorCode DMPlexGenerate_Tetgen(DM boundary, PetscBool interpolate, DM *dm)
461 {
462   MPI_Comm       comm;
463   DM_Plex       *mesh      = (DM_Plex *) boundary->data;
464   const PetscInt dim       = 3;
465   const char    *labelName = "marker";
466   ::tetgenio     in;
467   ::tetgenio     out;
468   DMLabel        label;
469   PetscInt       vStart, vEnd, v, fStart, fEnd, f;
470   PetscMPIInt    rank;
471   PetscErrorCode ierr;
472 
473   PetscFunctionBegin;
474   ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr);
475   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
476   ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr);
477   ierr = DMGetLabel(boundary, labelName, &label);CHKERRQ(ierr);
478 
479   in.numberofpoints = vEnd - vStart;
480   if (in.numberofpoints > 0) {
481     PetscSection coordSection;
482     Vec          coordinates;
483     PetscScalar *array;
484 
485     in.pointlist       = new double[in.numberofpoints*dim];
486     in.pointmarkerlist = new int[in.numberofpoints];
487 
488     ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr);
489     ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr);
490     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
491     for (v = vStart; v < vEnd; ++v) {
492       const PetscInt idx = v - vStart;
493       PetscInt       off, d;
494 
495       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
496       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
497       if (label) {
498         PetscInt val;
499 
500         ierr = DMLabelGetValue(label, v, &val);CHKERRQ(ierr);
501         in.pointmarkerlist[idx] = (int) val;
502       }
503     }
504     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
505   }
506   ierr  = DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);CHKERRQ(ierr);
507 
508   in.numberoffacets = fEnd - fStart;
509   if (in.numberoffacets > 0) {
510     in.facetlist       = new tetgenio::facet[in.numberoffacets];
511     in.facetmarkerlist = new int[in.numberoffacets];
512     for (f = fStart; f < fEnd; ++f) {
513       const PetscInt idx     = f - fStart;
514       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v;
515 
516       in.facetlist[idx].numberofpolygons = 1;
517       in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
518       in.facetlist[idx].numberofholes    = 0;
519       in.facetlist[idx].holelist         = NULL;
520 
521       ierr = DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
522       for (p = 0; p < numPoints*2; p += 2) {
523         const PetscInt point = points[p];
524         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
525       }
526 
527       tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
528       poly->numberofvertices = numVertices;
529       poly->vertexlist       = new int[poly->numberofvertices];
530       for (v = 0; v < numVertices; ++v) {
531         const PetscInt vIdx = points[v] - vStart;
532         poly->vertexlist[v] = vIdx;
533       }
534       if (label) {
535         PetscInt val;
536 
537         ierr = DMLabelGetValue(label, f, &val);CHKERRQ(ierr);
538         in.facetmarkerlist[idx] = (int) val;
539       }
540       ierr = DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
541     }
542   }
543   if (!rank) {
544     char args[32];
545 
546     /* Take away 'Q' for verbose output */
547     ierr = PetscStrcpy(args, "pqezQ");CHKERRQ(ierr);
548     if (mesh->tetgenOpts) {::tetrahedralize(mesh->tetgenOpts, &in, &out);}
549     else                  {::tetrahedralize(args, &in, &out);}
550   }
551   {
552     DMLabel        glabel      = NULL;
553     const PetscInt numCorners  = 4;
554     const PetscInt numCells    = out.numberoftetrahedra;
555     const PetscInt numVertices = out.numberofpoints;
556     const double   *meshCoords = out.pointlist;
557     int            *cells      = out.tetrahedronlist;
558 
559     ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr);
560     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr);
561     if (label) {ierr = DMCreateLabel(*dm, labelName); ierr = DMGetLabel(*dm, labelName, &glabel);}
562     /* Set labels */
563     for (v = 0; v < numVertices; ++v) {
564       if (out.pointmarkerlist[v]) {
565         if (glabel) {ierr = DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);}
566       }
567     }
568     if (interpolate) {
569 #if 0
570       PetscInt e;
571 
572       /* This check is never actually executed for ctetgen (which never returns edgemarkers) and seems to be broken for
573        * tetgen */
574       for (e = 0; e < out.numberofedges; e++) {
575         if (out.edgemarkerlist[e]) {
576           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
577           const PetscInt *edges;
578           PetscInt        numEdges;
579 
580           ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
581           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
582           if (glabel) {ierr = DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);}
583           ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
584         }
585       }
586 #endif
587       for (f = 0; f < out.numberoftrifaces; f++) {
588         if (out.trifacemarkerlist[f]) {
589           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
590           const PetscInt *faces;
591           PetscInt        numFaces;
592 
593           ierr = DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
594           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
595           if (glabel) {ierr = DMLabelSetValue(glabel, faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);}
596           ierr = DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
597         }
598       }
599     }
600     ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
601   }
602   PetscFunctionReturn(0);
603 }
604 
605 PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
606 {
607   MPI_Comm       comm;
608   const PetscInt dim       = 3;
609   const char    *labelName = "marker";
610   ::tetgenio     in;
611   ::tetgenio     out;
612   DMLabel        label;
613   PetscInt       vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
614   PetscMPIInt    rank;
615   PetscErrorCode ierr;
616 
617   PetscFunctionBegin;
618   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
619   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
620   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
621   ierr = MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
622   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
623   ierr = DMGetLabel(dm, labelName, &label);CHKERRQ(ierr);
624 
625   in.numberofpoints = vEnd - vStart;
626   if (in.numberofpoints > 0) {
627     PetscSection coordSection;
628     Vec          coordinates;
629     PetscScalar *array;
630 
631     in.pointlist       = new double[in.numberofpoints*dim];
632     in.pointmarkerlist = new int[in.numberofpoints];
633 
634     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
635     ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
636     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
637     for (v = vStart; v < vEnd; ++v) {
638       const PetscInt idx = v - vStart;
639       PetscInt       off, d;
640 
641       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
642       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
643       if (label) {
644         PetscInt val;
645 
646         ierr = DMLabelGetValue(label, v, &val);CHKERRQ(ierr);
647         in.pointmarkerlist[idx] = (int) val;
648       }
649     }
650     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
651   }
652   ierr  = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
653 
654   in.numberofcorners       = 4;
655   in.numberoftetrahedra    = cEnd - cStart;
656   in.tetrahedronvolumelist = (double*) maxVolumes;
657   if (in.numberoftetrahedra > 0) {
658     in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
659     for (c = cStart; c < cEnd; ++c) {
660       const PetscInt idx      = c - cStart;
661       PetscInt      *closure = NULL;
662       PetscInt       closureSize;
663 
664       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
665       if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
666       for (v = 0; v < 4; ++v) {
667         in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
668       }
669       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
670     }
671   }
672   /* TODO: Put in boundary faces with markers */
673   if (!rank) {
674     char args[32];
675 
676 #if 1
677     /* Take away 'Q' for verbose output */
678     ierr = PetscStrcpy(args, "qezQra");CHKERRQ(ierr);
679 #else
680     ierr = PetscStrcpy(args, "qezraVVVV");CHKERRQ(ierr);
681 #endif
682     ::tetrahedralize(args, &in, &out);
683   }
684   in.tetrahedronvolumelist = NULL;
685 
686   {
687     DMLabel        rlabel      = NULL;
688     const PetscInt numCorners  = 4;
689     const PetscInt numCells    = out.numberoftetrahedra;
690     const PetscInt numVertices = out.numberofpoints;
691     const double   *meshCoords = out.pointlist;
692     int            *cells      = out.tetrahedronlist;
693 
694     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
695 
696     ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr);
697     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr);
698     if (label) {ierr = DMCreateLabel(*dmRefined, labelName); ierr = DMGetLabel(*dmRefined, labelName, &rlabel);}
699     /* Set labels */
700     for (v = 0; v < numVertices; ++v) {
701       if (out.pointmarkerlist[v]) {
702         if (rlabel) {ierr = DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);}
703       }
704     }
705     if (interpolate) {
706       PetscInt f;
707 #if 0
708       PetscInt e;
709 
710       for (e = 0; e < out.numberofedges; e++) {
711         if (out.edgemarkerlist[e]) {
712           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
713           const PetscInt *edges;
714           PetscInt        numEdges;
715 
716           ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
717           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
718           if (rlabel) {ierr = DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);}
719           ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
720         }
721       }
722 #endif
723       for (f = 0; f < out.numberoftrifaces; f++) {
724         if (out.trifacemarkerlist[f]) {
725           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
726           const PetscInt *faces;
727           PetscInt        numFaces;
728 
729           ierr = DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
730           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
731           if (rlabel) {ierr = DMLabelSetValue(rlabel, faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);}
732           ierr = DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
733         }
734       }
735     }
736     ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr);
737   }
738   PetscFunctionReturn(0);
739 }
740 #endif /* PETSC_HAVE_TETGEN */
741 
742 #if defined(PETSC_HAVE_CTETGEN)
743 #include <ctetgen.h>
744 
745 PetscErrorCode DMPlexGenerate_CTetgen(DM boundary, PetscBool interpolate, DM *dm)
746 {
747   MPI_Comm       comm;
748   const PetscInt dim       = 3;
749   const char    *labelName = "marker";
750   PLC           *in, *out;
751   DMLabel        label;
752   PetscInt       verbose = 0, vStart, vEnd, v, fStart, fEnd, f;
753   PetscMPIInt    rank;
754   PetscErrorCode ierr;
755 
756   PetscFunctionBegin;
757   ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr);
758   ierr = PetscOptionsGetInt(NULL,((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);CHKERRQ(ierr);
759   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
760   ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr);
761   ierr = DMGetLabel(boundary, labelName, &label);CHKERRQ(ierr);
762   ierr = PLCCreate(&in);CHKERRQ(ierr);
763   ierr = PLCCreate(&out);CHKERRQ(ierr);
764 
765   in->numberofpoints = vEnd - vStart;
766   if (in->numberofpoints > 0) {
767     PetscSection coordSection;
768     Vec          coordinates;
769     PetscScalar *array;
770 
771     ierr = PetscMalloc1(in->numberofpoints*dim, &in->pointlist);CHKERRQ(ierr);
772     ierr = PetscMalloc1(in->numberofpoints,       &in->pointmarkerlist);CHKERRQ(ierr);
773     ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr);
774     ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr);
775     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
776     for (v = vStart; v < vEnd; ++v) {
777       const PetscInt idx = v - vStart;
778       PetscInt       off, d, m = -1;
779 
780       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
781       for (d = 0; d < dim; ++d) {
782         in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
783       }
784       if (label) {ierr = DMLabelGetValue(label, v, &m);CHKERRQ(ierr);}
785 
786       in->pointmarkerlist[idx] = (int) m;
787     }
788     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
789   }
790   ierr  = DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);CHKERRQ(ierr);
791 
792   in->numberoffacets = fEnd - fStart;
793   if (in->numberoffacets > 0) {
794     ierr = PetscMalloc1(in->numberoffacets, &in->facetlist);CHKERRQ(ierr);
795     ierr = PetscMalloc1(in->numberoffacets,   &in->facetmarkerlist);CHKERRQ(ierr);
796     for (f = fStart; f < fEnd; ++f) {
797       const PetscInt idx     = f - fStart;
798       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v, m = -1;
799       polygon       *poly;
800 
801       in->facetlist[idx].numberofpolygons = 1;
802 
803       ierr = PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist);CHKERRQ(ierr);
804 
805       in->facetlist[idx].numberofholes    = 0;
806       in->facetlist[idx].holelist         = NULL;
807 
808       ierr = DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
809       for (p = 0; p < numPoints*2; p += 2) {
810         const PetscInt point = points[p];
811         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
812       }
813 
814       poly                   = in->facetlist[idx].polygonlist;
815       poly->numberofvertices = numVertices;
816       ierr                   = PetscMalloc1(poly->numberofvertices, &poly->vertexlist);CHKERRQ(ierr);
817       for (v = 0; v < numVertices; ++v) {
818         const PetscInt vIdx = points[v] - vStart;
819         poly->vertexlist[v] = vIdx;
820       }
821       if (label) {ierr = DMLabelGetValue(label, f, &m);CHKERRQ(ierr);}
822       in->facetmarkerlist[idx] = (int) m;
823       ierr = DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
824     }
825   }
826   if (!rank) {
827     TetGenOpts t;
828 
829     ierr        = TetGenOptsInitialize(&t);CHKERRQ(ierr);
830     t.in        = boundary; /* Should go away */
831     t.plc       = 1;
832     t.quality   = 1;
833     t.edgesout  = 1;
834     t.zeroindex = 1;
835     t.quiet     = 1;
836     t.verbose   = verbose;
837     ierr        = TetGenCheckOpts(&t);CHKERRQ(ierr);
838     ierr        = TetGenTetrahedralize(&t, in, out);CHKERRQ(ierr);
839   }
840   {
841     DMLabel        glabel      = NULL;
842     const PetscInt numCorners  = 4;
843     const PetscInt numCells    = out->numberoftetrahedra;
844     const PetscInt numVertices = out->numberofpoints;
845     double         *meshCoords;
846     int            *cells      = out->tetrahedronlist;
847 
848     if (sizeof (PetscReal) == sizeof (double)) {
849       meshCoords = (double *) out->pointlist;
850     }
851     else {
852       PetscInt i;
853 
854       ierr = PetscMalloc1(3 * numVertices,&meshCoords);CHKERRQ(ierr);
855       for (i = 0; i < 3 * numVertices; i++) {
856         meshCoords[i] = (PetscReal) out->pointlist[i];
857       }
858     }
859 
860     ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr);
861     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr);
862     if (sizeof (PetscReal) != sizeof (double)) {
863       ierr = PetscFree(meshCoords);CHKERRQ(ierr);
864     }
865     if (label) {ierr = DMCreateLabel(*dm, labelName); ierr = DMGetLabel(*dm, labelName, &glabel);}
866     /* Set labels */
867     for (v = 0; v < numVertices; ++v) {
868       if (out->pointmarkerlist[v]) {
869         if (glabel) {ierr = DMLabelSetValue(glabel, v+numCells, out->pointmarkerlist[v]);CHKERRQ(ierr);}
870       }
871     }
872     if (interpolate) {
873       PetscInt e;
874 
875       for (e = 0; e < out->numberofedges; e++) {
876         if (out->edgemarkerlist[e]) {
877           const PetscInt  vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
878           const PetscInt *edges;
879           PetscInt        numEdges;
880 
881           ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
882           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
883           if (glabel) {ierr = DMLabelSetValue(glabel, edges[0], out->edgemarkerlist[e]);CHKERRQ(ierr);}
884           ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
885         }
886       }
887       for (f = 0; f < out->numberoftrifaces; f++) {
888         if (out->trifacemarkerlist[f]) {
889           const PetscInt  vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
890           const PetscInt *faces;
891           PetscInt        numFaces;
892 
893           ierr = DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
894           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
895           if (glabel) {ierr = DMLabelSetValue(glabel, faces[0], out->trifacemarkerlist[f]);CHKERRQ(ierr);}
896           ierr = DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
897         }
898       }
899     }
900     ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
901   }
902 
903   ierr = PLCDestroy(&in);CHKERRQ(ierr);
904   ierr = PLCDestroy(&out);CHKERRQ(ierr);
905   PetscFunctionReturn(0);
906 }
907 
908 PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
909 {
910   MPI_Comm       comm;
911   const PetscInt dim       = 3;
912   const char    *labelName = "marker";
913   PLC           *in, *out;
914   DMLabel        label;
915   PetscInt       verbose = 0, vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
916   PetscMPIInt    rank;
917   PetscErrorCode ierr;
918 
919   PetscFunctionBegin;
920   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
921   ierr = PetscOptionsGetInt(NULL,((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);CHKERRQ(ierr);
922   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
923   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
924   ierr = MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
925   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
926   ierr = DMGetLabel(dm, labelName, &label);CHKERRQ(ierr);
927   ierr = PLCCreate(&in);CHKERRQ(ierr);
928   ierr = PLCCreate(&out);CHKERRQ(ierr);
929 
930   in->numberofpoints = vEnd - vStart;
931   if (in->numberofpoints > 0) {
932     PetscSection coordSection;
933     Vec          coordinates;
934     PetscScalar *array;
935 
936     ierr = PetscMalloc1(in->numberofpoints*dim, &in->pointlist);CHKERRQ(ierr);
937     ierr = PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);CHKERRQ(ierr);
938     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
939     ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
940     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
941     for (v = vStart; v < vEnd; ++v) {
942       const PetscInt idx = v - vStart;
943       PetscInt       off, d, m = -1;
944 
945       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
946       for (d = 0; d < dim; ++d) {
947         in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
948       }
949       if (label) {ierr = DMLabelGetValue(label, v, &m);CHKERRQ(ierr);}
950 
951       in->pointmarkerlist[idx] = (int) m;
952     }
953     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
954   }
955   ierr  = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
956 
957   in->numberofcorners       = 4;
958   in->numberoftetrahedra    = cEnd - cStart;
959   in->tetrahedronvolumelist = maxVolumes;
960   if (in->numberoftetrahedra > 0) {
961     ierr = PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);CHKERRQ(ierr);
962     for (c = cStart; c < cEnd; ++c) {
963       const PetscInt idx      = c - cStart;
964       PetscInt      *closure = NULL;
965       PetscInt       closureSize;
966 
967       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
968       if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
969       for (v = 0; v < 4; ++v) {
970         in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
971       }
972       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
973     }
974   }
975   if (!rank) {
976     TetGenOpts t;
977 
978     ierr = TetGenOptsInitialize(&t);CHKERRQ(ierr);
979 
980     t.in        = dm; /* Should go away */
981     t.refine    = 1;
982     t.varvolume = 1;
983     t.quality   = 1;
984     t.edgesout  = 1;
985     t.zeroindex = 1;
986     t.quiet     = 1;
987     t.verbose   = verbose; /* Change this */
988 
989     ierr = TetGenCheckOpts(&t);CHKERRQ(ierr);
990     ierr = TetGenTetrahedralize(&t, in, out);CHKERRQ(ierr);
991   }
992   in->tetrahedronvolumelist = NULL;
993   {
994     DMLabel        rlabel      = NULL;
995     const PetscInt numCorners  = 4;
996     const PetscInt numCells    = out->numberoftetrahedra;
997     const PetscInt numVertices = out->numberofpoints;
998     double         *meshCoords;
999     int            *cells      = out->tetrahedronlist;
1000     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
1001 
1002     if (sizeof (PetscReal) == sizeof (double)) {
1003       meshCoords = (double *) out->pointlist;
1004     }
1005     else {
1006       PetscInt i;
1007 
1008       ierr = PetscMalloc1(3 * numVertices,&meshCoords);CHKERRQ(ierr);
1009       for (i = 0; i < 3 * numVertices; i++) {
1010         meshCoords[i] = (PetscReal) out->pointlist[i];
1011       }
1012     }
1013 
1014     ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr);
1015     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr);
1016     if (sizeof (PetscReal) != sizeof (double)) {
1017       ierr = PetscFree(meshCoords);CHKERRQ(ierr);
1018     }
1019     if (label) {ierr = DMCreateLabel(*dmRefined, labelName); ierr = DMGetLabel(*dmRefined, labelName, &rlabel);}
1020     /* Set labels */
1021     for (v = 0; v < numVertices; ++v) {
1022       if (out->pointmarkerlist[v]) {
1023         if (rlabel) {ierr = DMLabelSetValue(rlabel, v+numCells, out->pointmarkerlist[v]);CHKERRQ(ierr);}
1024       }
1025     }
1026     if (interpolate) {
1027       PetscInt e, f;
1028 
1029       for (e = 0; e < out->numberofedges; e++) {
1030         if (out->edgemarkerlist[e]) {
1031           const PetscInt  vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
1032           const PetscInt *edges;
1033           PetscInt        numEdges;
1034 
1035           ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
1036           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
1037           if (rlabel) {ierr = DMLabelSetValue(rlabel, edges[0], out->edgemarkerlist[e]);CHKERRQ(ierr);}
1038           ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
1039         }
1040       }
1041       for (f = 0; f < out->numberoftrifaces; f++) {
1042         if (out->trifacemarkerlist[f]) {
1043           const PetscInt  vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
1044           const PetscInt *faces;
1045           PetscInt        numFaces;
1046 
1047           ierr = DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
1048           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
1049           if (rlabel) {ierr = DMLabelSetValue(rlabel, faces[0], out->trifacemarkerlist[f]);CHKERRQ(ierr);}
1050           ierr = DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
1051         }
1052       }
1053     }
1054     ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr);
1055   }
1056   ierr = PLCDestroy(&in);CHKERRQ(ierr);
1057   ierr = PLCDestroy(&out);CHKERRQ(ierr);
1058   PetscFunctionReturn(0);
1059 }
1060 #endif /* PETSC_HAVE_CTETGEN */
1061 
1062 /*@C
1063   DMPlexGenerate - Generates a mesh.
1064 
1065   Not Collective
1066 
1067   Input Parameters:
1068 + boundary - The DMPlex boundary object
1069 . name - The mesh generation package name
1070 - interpolate - Flag to create intermediate mesh elements
1071 
1072   Output Parameter:
1073 . mesh - The DMPlex object
1074 
1075   Level: intermediate
1076 
1077 .keywords: mesh, elements
1078 .seealso: DMPlexCreate(), DMRefine()
1079 @*/
1080 PetscErrorCode DMPlexGenerate(DM boundary, const char name[], PetscBool interpolate, DM *mesh)
1081 {
1082   PetscInt       dim;
1083   char           genname[1024];
1084   PetscBool      isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;
1085   PetscErrorCode ierr;
1086 
1087   PetscFunctionBegin;
1088   PetscValidHeaderSpecific(boundary, DM_CLASSID, 1);
1089   PetscValidLogicalCollectiveBool(boundary, interpolate, 2);
1090   ierr = DMGetDimension(boundary, &dim);CHKERRQ(ierr);
1091   ierr = PetscOptionsGetString(((PetscObject) boundary)->options,((PetscObject) boundary)->prefix, "-dm_plex_generator", genname, 1024, &flg);CHKERRQ(ierr);
1092   if (flg) name = genname;
1093   if (name) {
1094     ierr = PetscStrcmp(name, "triangle", &isTriangle);CHKERRQ(ierr);
1095     ierr = PetscStrcmp(name, "tetgen",   &isTetgen);CHKERRQ(ierr);
1096     ierr = PetscStrcmp(name, "ctetgen",  &isCTetgen);CHKERRQ(ierr);
1097   }
1098   switch (dim) {
1099   case 1:
1100     if (!name || isTriangle) {
1101 #if defined(PETSC_HAVE_TRIANGLE)
1102       ierr = DMPlexGenerate_Triangle(boundary, interpolate, mesh);CHKERRQ(ierr);
1103 #else
1104       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation needs external package support.\nPlease reconfigure with --download-triangle.");
1105 #endif
1106     } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
1107     break;
1108   case 2:
1109     if (!name) {
1110 #if defined(PETSC_HAVE_CTETGEN)
1111       ierr = DMPlexGenerate_CTetgen(boundary, interpolate, mesh);CHKERRQ(ierr);
1112 #elif defined(PETSC_HAVE_TETGEN)
1113       ierr = DMPlexGenerate_Tetgen(boundary, interpolate, mesh);CHKERRQ(ierr);
1114 #else
1115       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "External package CTetgen or Tetgen needed.\nPlease reconfigure with '--download-ctetgen' or '--with-clanguage=cxx --download-tetgen'.");
1116 #endif
1117     } else if (isCTetgen) {
1118 #if defined(PETSC_HAVE_CTETGEN)
1119       ierr = DMPlexGenerate_CTetgen(boundary, interpolate, mesh);CHKERRQ(ierr);
1120 #else
1121       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
1122 #endif
1123     } else if (isTetgen) {
1124 #if defined(PETSC_HAVE_TETGEN)
1125       ierr = DMPlexGenerate_Tetgen(boundary, interpolate, mesh);CHKERRQ(ierr);
1126 #else
1127       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --download-tetgen.");
1128 #endif
1129     } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
1130     break;
1131   default:
1132     SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation for a dimension %d boundary is not supported.", dim);
1133   }
1134   PetscFunctionReturn(0);
1135 }
1136