xref: /petsc/src/dm/impls/plex/plexgenerate.c (revision 1118d4bc1f37eaf681057353470764a7fc73fc53) 
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 static 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   struct triangulateio in;
193   struct triangulateio out;
194   DMLabel              label;
195   PetscInt             vStart, vEnd, v, eStart, eEnd, e;
196   PetscMPIInt          rank;
197   PetscErrorCode       ierr;
198 
199   PetscFunctionBegin;
200   ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr);
201   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
202   ierr = InitInput_Triangle(&in);CHKERRQ(ierr);
203   ierr = InitOutput_Triangle(&out);CHKERRQ(ierr);
204   ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr);
205   ierr = DMGetLabel(boundary, labelName, &label);CHKERRQ(ierr);
206 
207   in.numberofpoints = vEnd - vStart;
208   if (in.numberofpoints > 0) {
209     PetscSection coordSection;
210     Vec          coordinates;
211     PetscScalar *array;
212 
213     ierr = PetscMalloc1(in.numberofpoints*dim, &in.pointlist);CHKERRQ(ierr);
214     ierr = PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);CHKERRQ(ierr);
215     ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr);
216     ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr);
217     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
218     for (v = vStart; v < vEnd; ++v) {
219       const PetscInt idx = v - vStart;
220       PetscInt       off, d;
221 
222       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
223       for (d = 0; d < dim; ++d) {
224         in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
225       }
226       if (label) {ierr = DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);CHKERRQ(ierr);}
227     }
228     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
229   }
230   ierr  = DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);CHKERRQ(ierr);
231   in.numberofsegments = eEnd - eStart;
232   if (in.numberofsegments > 0) {
233     ierr = PetscMalloc1(in.numberofsegments*2, &in.segmentlist);CHKERRQ(ierr);
234     ierr = PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist);CHKERRQ(ierr);
235     for (e = eStart; e < eEnd; ++e) {
236       const PetscInt  idx = e - eStart;
237       const PetscInt *cone;
238 
239       ierr = DMPlexGetCone(boundary, e, &cone);CHKERRQ(ierr);
240 
241       in.segmentlist[idx*2+0] = cone[0] - vStart;
242       in.segmentlist[idx*2+1] = cone[1] - vStart;
243 
244       if (label) {ierr = DMLabelGetValue(label, e, &in.segmentmarkerlist[idx]);CHKERRQ(ierr);}
245     }
246   }
247 #if 0 /* Do not currently support holes */
248   PetscReal *holeCoords;
249   PetscInt   h, d;
250 
251   ierr = DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);CHKERRQ(ierr);
252   if (in.numberofholes > 0) {
253     ierr = PetscMalloc1(in.numberofholes*dim, &in.holelist);CHKERRQ(ierr);
254     for (h = 0; h < in.numberofholes; ++h) {
255       for (d = 0; d < dim; ++d) {
256         in.holelist[h*dim+d] = holeCoords[h*dim+d];
257       }
258     }
259   }
260 #endif
261   if (!rank) {
262     char args[32];
263 
264     /* Take away 'Q' for verbose output */
265     ierr = PetscStrcpy(args, "pqezQ");CHKERRQ(ierr);
266     if (createConvexHull)   {ierr = PetscStrcat(args, "c");CHKERRQ(ierr);}
267     if (constrained)        {ierr = PetscStrcpy(args, "zepDQ");CHKERRQ(ierr);}
268     if (mesh->triangleOpts) {triangulate(mesh->triangleOpts, &in, &out, NULL);}
269     else                    {triangulate(args, &in, &out, NULL);}
270   }
271   ierr = PetscFree(in.pointlist);CHKERRQ(ierr);
272   ierr = PetscFree(in.pointmarkerlist);CHKERRQ(ierr);
273   ierr = PetscFree(in.segmentlist);CHKERRQ(ierr);
274   ierr = PetscFree(in.segmentmarkerlist);CHKERRQ(ierr);
275   ierr = PetscFree(in.holelist);CHKERRQ(ierr);
276 
277   {
278     DMLabel        glabel      = NULL;
279     const PetscInt numCorners  = 3;
280     const PetscInt numCells    = out.numberoftriangles;
281     const PetscInt numVertices = out.numberofpoints;
282     const int     *cells      = out.trianglelist;
283     const double  *meshCoords = out.pointlist;
284 
285     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr);
286     if (label) {ierr = DMCreateLabel(*dm, labelName); ierr = DMGetLabel(*dm, labelName, &glabel);}
287     /* Set labels */
288     for (v = 0; v < numVertices; ++v) {
289       if (out.pointmarkerlist[v]) {
290         if (glabel) {ierr = DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);}
291       }
292     }
293     if (interpolate) {
294       for (e = 0; e < out.numberofedges; e++) {
295         if (out.edgemarkerlist[e]) {
296           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
297           const PetscInt *edges;
298           PetscInt        numEdges;
299 
300           ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
301           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
302           if (glabel) {ierr = DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);}
303           ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
304         }
305       }
306     }
307     ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
308   }
309 #if 0 /* Do not currently support holes */
310   ierr = DMPlexCopyHoles(*dm, boundary);CHKERRQ(ierr);
311 #endif
312   ierr = FiniOutput_Triangle(&out);CHKERRQ(ierr);
313   PetscFunctionReturn(0);
314 }
315 
316 static PetscErrorCode DMPlexRefine_Triangle(DM dm, double *maxVolumes, DM *dmRefined)
317 {
318   MPI_Comm             comm;
319   PetscInt             dim       = 2;
320   const char          *labelName = "marker";
321   struct triangulateio in;
322   struct triangulateio out;
323   DMLabel              label;
324   PetscInt             vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
325   PetscMPIInt          rank;
326   PetscErrorCode       ierr;
327 
328   PetscFunctionBegin;
329   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
330   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
331   ierr = InitInput_Triangle(&in);CHKERRQ(ierr);
332   ierr = InitOutput_Triangle(&out);CHKERRQ(ierr);
333   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
334   ierr = MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
335   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
336   ierr = DMGetLabel(dm, labelName, &label);CHKERRQ(ierr);
337 
338   in.numberofpoints = vEnd - vStart;
339   if (in.numberofpoints > 0) {
340     PetscSection coordSection;
341     Vec          coordinates;
342     PetscScalar *array;
343 
344     ierr = PetscMalloc1(in.numberofpoints*dim, &in.pointlist);CHKERRQ(ierr);
345     ierr = PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);CHKERRQ(ierr);
346     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
347     ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
348     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
349     for (v = vStart; v < vEnd; ++v) {
350       const PetscInt idx = v - vStart;
351       PetscInt       off, d;
352 
353       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
354       for (d = 0; d < dim; ++d) {
355         in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
356       }
357       if (label) {ierr = DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);CHKERRQ(ierr);}
358     }
359     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
360   }
361   ierr  = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
362 
363   in.numberofcorners   = 3;
364   in.numberoftriangles = cEnd - cStart;
365 
366   in.trianglearealist  = (double*) maxVolumes;
367   if (in.numberoftriangles > 0) {
368     ierr = PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);CHKERRQ(ierr);
369     for (c = cStart; c < cEnd; ++c) {
370       const PetscInt idx      = c - cStart;
371       PetscInt      *closure = NULL;
372       PetscInt       closureSize;
373 
374       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
375       if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize);
376       for (v = 0; v < 3; ++v) {
377         in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart;
378       }
379       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
380     }
381   }
382   /* TODO: Segment markers are missing on input */
383 #if 0 /* Do not currently support holes */
384   PetscReal *holeCoords;
385   PetscInt   h, d;
386 
387   ierr = DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);CHKERRQ(ierr);
388   if (in.numberofholes > 0) {
389     ierr = PetscMalloc1(in.numberofholes*dim, &in.holelist);CHKERRQ(ierr);
390     for (h = 0; h < in.numberofholes; ++h) {
391       for (d = 0; d < dim; ++d) {
392         in.holelist[h*dim+d] = holeCoords[h*dim+d];
393       }
394     }
395   }
396 #endif
397   if (!rank) {
398     char args[32];
399 
400     /* Take away 'Q' for verbose output */
401     ierr = PetscStrcpy(args, "pqezQra");CHKERRQ(ierr);
402     triangulate(args, &in, &out, NULL);
403   }
404   ierr = PetscFree(in.pointlist);CHKERRQ(ierr);
405   ierr = PetscFree(in.pointmarkerlist);CHKERRQ(ierr);
406   ierr = PetscFree(in.segmentlist);CHKERRQ(ierr);
407   ierr = PetscFree(in.segmentmarkerlist);CHKERRQ(ierr);
408   ierr = PetscFree(in.trianglelist);CHKERRQ(ierr);
409 
410   {
411     DMLabel        rlabel      = NULL;
412     const PetscInt numCorners  = 3;
413     const PetscInt numCells    = out.numberoftriangles;
414     const PetscInt numVertices = out.numberofpoints;
415     const int     *cells      = out.trianglelist;
416     const double  *meshCoords = out.pointlist;
417     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
418 
419     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr);
420     if (label) {ierr = DMCreateLabel(*dmRefined, labelName); ierr = DMGetLabel(*dmRefined, labelName, &rlabel);}
421     /* Set labels */
422     for (v = 0; v < numVertices; ++v) {
423       if (out.pointmarkerlist[v]) {
424         if (rlabel) {ierr = DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);}
425       }
426     }
427     if (interpolate) {
428       PetscInt e;
429 
430       for (e = 0; e < out.numberofedges; e++) {
431         if (out.edgemarkerlist[e]) {
432           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
433           const PetscInt *edges;
434           PetscInt        numEdges;
435 
436           ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
437           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
438           if (rlabel) {ierr = DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);}
439           ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
440         }
441       }
442     }
443     ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr);
444   }
445 #if 0 /* Do not currently support holes */
446   ierr = DMPlexCopyHoles(*dm, boundary);CHKERRQ(ierr);
447 #endif
448   ierr = FiniOutput_Triangle(&out);CHKERRQ(ierr);
449   PetscFunctionReturn(0);
450 }
451 #endif /* PETSC_HAVE_TRIANGLE */
452 
453 #if defined(PETSC_HAVE_TETGEN)
454 #include <tetgen.h>
455 static PetscErrorCode DMPlexGenerate_Tetgen(DM boundary, PetscBool interpolate, DM *dm)
456 {
457   MPI_Comm       comm;
458   DM_Plex       *mesh      = (DM_Plex *) boundary->data;
459   const PetscInt dim       = 3;
460   const char    *labelName = "marker";
461   ::tetgenio     in;
462   ::tetgenio     out;
463   DMLabel        label;
464   PetscInt       vStart, vEnd, v, fStart, fEnd, f;
465   PetscMPIInt    rank;
466   PetscErrorCode ierr;
467 
468   PetscFunctionBegin;
469   ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr);
470   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
471   ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr);
472   ierr = DMGetLabel(boundary, labelName, &label);CHKERRQ(ierr);
473 
474   in.numberofpoints = vEnd - vStart;
475   if (in.numberofpoints > 0) {
476     PetscSection coordSection;
477     Vec          coordinates;
478     PetscScalar *array;
479 
480     in.pointlist       = new double[in.numberofpoints*dim];
481     in.pointmarkerlist = new int[in.numberofpoints];
482 
483     ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr);
484     ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr);
485     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
486     for (v = vStart; v < vEnd; ++v) {
487       const PetscInt idx = v - vStart;
488       PetscInt       off, d;
489 
490       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
491       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
492       if (label) {
493         PetscInt val;
494 
495         ierr = DMLabelGetValue(label, v, &val);CHKERRQ(ierr);
496         in.pointmarkerlist[idx] = (int) val;
497       }
498     }
499     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
500   }
501   ierr  = DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);CHKERRQ(ierr);
502 
503   in.numberoffacets = fEnd - fStart;
504   if (in.numberoffacets > 0) {
505     in.facetlist       = new tetgenio::facet[in.numberoffacets];
506     in.facetmarkerlist = new int[in.numberoffacets];
507     for (f = fStart; f < fEnd; ++f) {
508       const PetscInt idx     = f - fStart;
509       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v;
510 
511       in.facetlist[idx].numberofpolygons = 1;
512       in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
513       in.facetlist[idx].numberofholes    = 0;
514       in.facetlist[idx].holelist         = NULL;
515 
516       ierr = DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
517       for (p = 0; p < numPoints*2; p += 2) {
518         const PetscInt point = points[p];
519         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
520       }
521 
522       tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
523       poly->numberofvertices = numVertices;
524       poly->vertexlist       = new int[poly->numberofvertices];
525       for (v = 0; v < numVertices; ++v) {
526         const PetscInt vIdx = points[v] - vStart;
527         poly->vertexlist[v] = vIdx;
528       }
529       if (label) {
530         PetscInt val;
531 
532         ierr = DMLabelGetValue(label, f, &val);CHKERRQ(ierr);
533         in.facetmarkerlist[idx] = (int) val;
534       }
535       ierr = DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
536     }
537   }
538   if (!rank) {
539     char args[32];
540 
541     /* Take away 'Q' for verbose output */
542     ierr = PetscStrcpy(args, "pqezQ");CHKERRQ(ierr);
543     if (mesh->tetgenOpts) {::tetrahedralize(mesh->tetgenOpts, &in, &out);}
544     else                  {::tetrahedralize(args, &in, &out);}
545   }
546   {
547     DMLabel        glabel      = NULL;
548     const PetscInt numCorners  = 4;
549     const PetscInt numCells    = out.numberoftetrahedra;
550     const PetscInt numVertices = out.numberofpoints;
551     const double   *meshCoords = out.pointlist;
552     int            *cells      = out.tetrahedronlist;
553 
554     ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr);
555     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr);
556     if (label) {ierr = DMCreateLabel(*dm, labelName); ierr = DMGetLabel(*dm, labelName, &glabel);}
557     /* Set labels */
558     for (v = 0; v < numVertices; ++v) {
559       if (out.pointmarkerlist[v]) {
560         if (glabel) {ierr = DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);}
561       }
562     }
563     if (interpolate) {
564 #if 0
565       PetscInt e;
566 
567       /* This check is never actually executed for ctetgen (which never returns edgemarkers) and seems to be broken for
568        * tetgen */
569       for (e = 0; e < out.numberofedges; e++) {
570         if (out.edgemarkerlist[e]) {
571           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
572           const PetscInt *edges;
573           PetscInt        numEdges;
574 
575           ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
576           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
577           if (glabel) {ierr = DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);}
578           ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
579         }
580       }
581 #endif
582       for (f = 0; f < out.numberoftrifaces; f++) {
583         if (out.trifacemarkerlist[f]) {
584           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
585           const PetscInt *faces;
586           PetscInt        numFaces;
587 
588           ierr = DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
589           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
590           if (glabel) {ierr = DMLabelSetValue(glabel, faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);}
591           ierr = DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
592         }
593       }
594     }
595     ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
596   }
597   PetscFunctionReturn(0);
598 }
599 
600 static PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
601 {
602   MPI_Comm       comm;
603   const PetscInt dim       = 3;
604   const char    *labelName = "marker";
605   ::tetgenio     in;
606   ::tetgenio     out;
607   DMLabel        label;
608   PetscInt       vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
609   PetscMPIInt    rank;
610   PetscErrorCode ierr;
611 
612   PetscFunctionBegin;
613   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
614   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
615   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
616   ierr = MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
617   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
618   ierr = DMGetLabel(dm, labelName, &label);CHKERRQ(ierr);
619 
620   in.numberofpoints = vEnd - vStart;
621   if (in.numberofpoints > 0) {
622     PetscSection coordSection;
623     Vec          coordinates;
624     PetscScalar *array;
625 
626     in.pointlist       = new double[in.numberofpoints*dim];
627     in.pointmarkerlist = new int[in.numberofpoints];
628 
629     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
630     ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
631     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
632     for (v = vStart; v < vEnd; ++v) {
633       const PetscInt idx = v - vStart;
634       PetscInt       off, d;
635 
636       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
637       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
638       if (label) {
639         PetscInt val;
640 
641         ierr = DMLabelGetValue(label, v, &val);CHKERRQ(ierr);
642         in.pointmarkerlist[idx] = (int) val;
643       }
644     }
645     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
646   }
647   ierr  = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
648 
649   in.numberofcorners       = 4;
650   in.numberoftetrahedra    = cEnd - cStart;
651   in.tetrahedronvolumelist = (double*) maxVolumes;
652   if (in.numberoftetrahedra > 0) {
653     in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
654     for (c = cStart; c < cEnd; ++c) {
655       const PetscInt idx      = c - cStart;
656       PetscInt      *closure = NULL;
657       PetscInt       closureSize;
658 
659       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
660       if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
661       for (v = 0; v < 4; ++v) {
662         in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
663       }
664       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
665     }
666   }
667   /* TODO: Put in boundary faces with markers */
668   if (!rank) {
669     char args[32];
670 
671 #if 1
672     /* Take away 'Q' for verbose output */
673     ierr = PetscStrcpy(args, "qezQra");CHKERRQ(ierr);
674 #else
675     ierr = PetscStrcpy(args, "qezraVVVV");CHKERRQ(ierr);
676 #endif
677     ::tetrahedralize(args, &in, &out);
678   }
679   in.tetrahedronvolumelist = NULL;
680 
681   {
682     DMLabel        rlabel      = NULL;
683     const PetscInt numCorners  = 4;
684     const PetscInt numCells    = out.numberoftetrahedra;
685     const PetscInt numVertices = out.numberofpoints;
686     const double   *meshCoords = out.pointlist;
687     int            *cells      = out.tetrahedronlist;
688 
689     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
690 
691     ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr);
692     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr);
693     if (label) {ierr = DMCreateLabel(*dmRefined, labelName); ierr = DMGetLabel(*dmRefined, labelName, &rlabel);}
694     /* Set labels */
695     for (v = 0; v < numVertices; ++v) {
696       if (out.pointmarkerlist[v]) {
697         if (rlabel) {ierr = DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);}
698       }
699     }
700     if (interpolate) {
701       PetscInt f;
702 #if 0
703       PetscInt e;
704 
705       for (e = 0; e < out.numberofedges; e++) {
706         if (out.edgemarkerlist[e]) {
707           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
708           const PetscInt *edges;
709           PetscInt        numEdges;
710 
711           ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
712           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
713           if (rlabel) {ierr = DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);}
714           ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
715         }
716       }
717 #endif
718       for (f = 0; f < out.numberoftrifaces; f++) {
719         if (out.trifacemarkerlist[f]) {
720           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
721           const PetscInt *faces;
722           PetscInt        numFaces;
723 
724           ierr = DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
725           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
726           if (rlabel) {ierr = DMLabelSetValue(rlabel, faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);}
727           ierr = DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
728         }
729       }
730     }
731     ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr);
732   }
733   PetscFunctionReturn(0);
734 }
735 #endif /* PETSC_HAVE_TETGEN */
736 
737 #if defined(PETSC_HAVE_CTETGEN)
738 #include <ctetgen.h>
739 
740 static PetscErrorCode DMPlexGenerate_CTetgen(DM boundary, PetscBool interpolate, DM *dm)
741 {
742   MPI_Comm       comm;
743   const PetscInt dim       = 3;
744   const char    *labelName = "marker";
745   PLC           *in, *out;
746   DMLabel        label;
747   PetscInt       verbose = 0, vStart, vEnd, v, fStart, fEnd, f;
748   PetscMPIInt    rank;
749   PetscErrorCode ierr;
750 
751   PetscFunctionBegin;
752   ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr);
753   ierr = PetscOptionsGetInt(NULL,((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);CHKERRQ(ierr);
754   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
755   ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr);
756   ierr = DMGetLabel(boundary, labelName, &label);CHKERRQ(ierr);
757   ierr = PLCCreate(&in);CHKERRQ(ierr);
758   ierr = PLCCreate(&out);CHKERRQ(ierr);
759 
760   in->numberofpoints = vEnd - vStart;
761   if (in->numberofpoints > 0) {
762     PetscSection coordSection;
763     Vec          coordinates;
764     PetscScalar *array;
765 
766     ierr = PetscMalloc1(in->numberofpoints*dim, &in->pointlist);CHKERRQ(ierr);
767     ierr = PetscMalloc1(in->numberofpoints,       &in->pointmarkerlist);CHKERRQ(ierr);
768     ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr);
769     ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr);
770     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
771     for (v = vStart; v < vEnd; ++v) {
772       const PetscInt idx = v - vStart;
773       PetscInt       off, d, m = -1;
774 
775       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
776       for (d = 0; d < dim; ++d) {
777         in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
778       }
779       if (label) {ierr = DMLabelGetValue(label, v, &m);CHKERRQ(ierr);}
780 
781       in->pointmarkerlist[idx] = (int) m;
782     }
783     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
784   }
785   ierr  = DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);CHKERRQ(ierr);
786 
787   in->numberoffacets = fEnd - fStart;
788   if (in->numberoffacets > 0) {
789     ierr = PetscMalloc1(in->numberoffacets, &in->facetlist);CHKERRQ(ierr);
790     ierr = PetscMalloc1(in->numberoffacets,   &in->facetmarkerlist);CHKERRQ(ierr);
791     for (f = fStart; f < fEnd; ++f) {
792       const PetscInt idx     = f - fStart;
793       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v, m = -1;
794       polygon       *poly;
795 
796       in->facetlist[idx].numberofpolygons = 1;
797 
798       ierr = PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist);CHKERRQ(ierr);
799 
800       in->facetlist[idx].numberofholes    = 0;
801       in->facetlist[idx].holelist         = NULL;
802 
803       ierr = DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
804       for (p = 0; p < numPoints*2; p += 2) {
805         const PetscInt point = points[p];
806         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
807       }
808 
809       poly                   = in->facetlist[idx].polygonlist;
810       poly->numberofvertices = numVertices;
811       ierr                   = PetscMalloc1(poly->numberofvertices, &poly->vertexlist);CHKERRQ(ierr);
812       for (v = 0; v < numVertices; ++v) {
813         const PetscInt vIdx = points[v] - vStart;
814         poly->vertexlist[v] = vIdx;
815       }
816       if (label) {ierr = DMLabelGetValue(label, f, &m);CHKERRQ(ierr);}
817       in->facetmarkerlist[idx] = (int) m;
818       ierr = DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
819     }
820   }
821   if (!rank) {
822     TetGenOpts t;
823 
824     ierr        = TetGenOptsInitialize(&t);CHKERRQ(ierr);
825     t.in        = boundary; /* Should go away */
826     t.plc       = 1;
827     t.quality   = 1;
828     t.edgesout  = 1;
829     t.zeroindex = 1;
830     t.quiet     = 1;
831     t.verbose   = verbose;
832     ierr        = TetGenCheckOpts(&t);CHKERRQ(ierr);
833     ierr        = TetGenTetrahedralize(&t, in, out);CHKERRQ(ierr);
834   }
835   {
836     DMLabel        glabel      = NULL;
837     const PetscInt numCorners  = 4;
838     const PetscInt numCells    = out->numberoftetrahedra;
839     const PetscInt numVertices = out->numberofpoints;
840     double         *meshCoords;
841     int            *cells      = out->tetrahedronlist;
842 
843     if (sizeof (PetscReal) == sizeof (double)) {
844       meshCoords = (double *) out->pointlist;
845     }
846     else {
847       PetscInt i;
848 
849       ierr = PetscMalloc1(3 * numVertices,&meshCoords);CHKERRQ(ierr);
850       for (i = 0; i < 3 * numVertices; i++) {
851         meshCoords[i] = (PetscReal) out->pointlist[i];
852       }
853     }
854 
855     ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr);
856     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr);
857     if (sizeof (PetscReal) != sizeof (double)) {
858       ierr = PetscFree(meshCoords);CHKERRQ(ierr);
859     }
860     if (label) {ierr = DMCreateLabel(*dm, labelName); ierr = DMGetLabel(*dm, labelName, &glabel);}
861     /* Set labels */
862     for (v = 0; v < numVertices; ++v) {
863       if (out->pointmarkerlist[v]) {
864         if (glabel) {ierr = DMLabelSetValue(glabel, v+numCells, out->pointmarkerlist[v]);CHKERRQ(ierr);}
865       }
866     }
867     if (interpolate) {
868       PetscInt e;
869 
870       for (e = 0; e < out->numberofedges; e++) {
871         if (out->edgemarkerlist[e]) {
872           const PetscInt  vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
873           const PetscInt *edges;
874           PetscInt        numEdges;
875 
876           ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
877           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
878           if (glabel) {ierr = DMLabelSetValue(glabel, edges[0], out->edgemarkerlist[e]);CHKERRQ(ierr);}
879           ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
880         }
881       }
882       for (f = 0; f < out->numberoftrifaces; f++) {
883         if (out->trifacemarkerlist[f]) {
884           const PetscInt  vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
885           const PetscInt *faces;
886           PetscInt        numFaces;
887 
888           ierr = DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
889           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
890           if (glabel) {ierr = DMLabelSetValue(glabel, faces[0], out->trifacemarkerlist[f]);CHKERRQ(ierr);}
891           ierr = DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
892         }
893       }
894     }
895     ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
896   }
897 
898   ierr = PLCDestroy(&in);CHKERRQ(ierr);
899   ierr = PLCDestroy(&out);CHKERRQ(ierr);
900   PetscFunctionReturn(0);
901 }
902 
903 static PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
904 {
905   MPI_Comm       comm;
906   const PetscInt dim       = 3;
907   const char    *labelName = "marker";
908   PLC           *in, *out;
909   DMLabel        label;
910   PetscInt       verbose = 0, vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
911   PetscMPIInt    rank;
912   PetscErrorCode ierr;
913 
914   PetscFunctionBegin;
915   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
916   ierr = PetscOptionsGetInt(NULL,((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);CHKERRQ(ierr);
917   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
918   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
919   ierr = MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
920   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
921   ierr = DMGetLabel(dm, labelName, &label);CHKERRQ(ierr);
922   ierr = PLCCreate(&in);CHKERRQ(ierr);
923   ierr = PLCCreate(&out);CHKERRQ(ierr);
924 
925   in->numberofpoints = vEnd - vStart;
926   if (in->numberofpoints > 0) {
927     PetscSection coordSection;
928     Vec          coordinates;
929     PetscScalar *array;
930 
931     ierr = PetscMalloc1(in->numberofpoints*dim, &in->pointlist);CHKERRQ(ierr);
932     ierr = PetscMalloc1(in->numberofpoints,       &in->pointmarkerlist);CHKERRQ(ierr);
933     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
934     ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
935     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
936     for (v = vStart; v < vEnd; ++v) {
937       const PetscInt idx = v - vStart;
938       PetscInt       off, d, m = -1;
939 
940       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
941       for (d = 0; d < dim; ++d) {
942         in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
943       }
944       if (label) {ierr = DMLabelGetValue(label, v, &m);CHKERRQ(ierr);}
945 
946       in->pointmarkerlist[idx] = (int) m;
947     }
948     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
949   }
950   ierr  = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
951 
952   in->numberofcorners       = 4;
953   in->numberoftetrahedra    = cEnd - cStart;
954   in->tetrahedronvolumelist = maxVolumes;
955   if (in->numberoftetrahedra > 0) {
956     ierr = PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);CHKERRQ(ierr);
957     for (c = cStart; c < cEnd; ++c) {
958       const PetscInt idx      = c - cStart;
959       PetscInt      *closure = NULL;
960       PetscInt       closureSize;
961 
962       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
963       if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
964       for (v = 0; v < 4; ++v) {
965         in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
966       }
967       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
968     }
969   }
970   if (!rank) {
971     TetGenOpts t;
972 
973     ierr = TetGenOptsInitialize(&t);CHKERRQ(ierr);
974 
975     t.in        = dm; /* Should go away */
976     t.refine    = 1;
977     t.varvolume = 1;
978     t.quality   = 1;
979     t.edgesout  = 1;
980     t.zeroindex = 1;
981     t.quiet     = 1;
982     t.verbose   = verbose; /* Change this */
983 
984     ierr = TetGenCheckOpts(&t);CHKERRQ(ierr);
985     ierr = TetGenTetrahedralize(&t, in, out);CHKERRQ(ierr);
986   }
987   in->tetrahedronvolumelist = NULL;
988   {
989     DMLabel        rlabel      = NULL;
990     const PetscInt numCorners  = 4;
991     const PetscInt numCells    = out->numberoftetrahedra;
992     const PetscInt numVertices = out->numberofpoints;
993     double         *meshCoords;
994     int            *cells      = out->tetrahedronlist;
995     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
996 
997     if (sizeof (PetscReal) == sizeof (double)) {
998       meshCoords = (double *) out->pointlist;
999     }
1000     else {
1001       PetscInt i;
1002 
1003       ierr = PetscMalloc1(3 * numVertices,&meshCoords);CHKERRQ(ierr);
1004       for (i = 0; i < 3 * numVertices; i++) {
1005         meshCoords[i] = (PetscReal) out->pointlist[i];
1006       }
1007     }
1008 
1009     ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr);
1010     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr);
1011     if (sizeof (PetscReal) != sizeof (double)) {
1012       ierr = PetscFree(meshCoords);CHKERRQ(ierr);
1013     }
1014     if (label) {ierr = DMCreateLabel(*dmRefined, labelName); ierr = DMGetLabel(*dmRefined, labelName, &rlabel);}
1015     /* Set labels */
1016     for (v = 0; v < numVertices; ++v) {
1017       if (out->pointmarkerlist[v]) {
1018         if (rlabel) {ierr = DMLabelSetValue(rlabel, v+numCells, out->pointmarkerlist[v]);CHKERRQ(ierr);}
1019       }
1020     }
1021     if (interpolate) {
1022       PetscInt e, f;
1023 
1024       for (e = 0; e < out->numberofedges; e++) {
1025         if (out->edgemarkerlist[e]) {
1026           const PetscInt  vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
1027           const PetscInt *edges;
1028           PetscInt        numEdges;
1029 
1030           ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
1031           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
1032           if (rlabel) {ierr = DMLabelSetValue(rlabel, edges[0], out->edgemarkerlist[e]);CHKERRQ(ierr);}
1033           ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
1034         }
1035       }
1036       for (f = 0; f < out->numberoftrifaces; f++) {
1037         if (out->trifacemarkerlist[f]) {
1038           const PetscInt  vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
1039           const PetscInt *faces;
1040           PetscInt        numFaces;
1041 
1042           ierr = DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
1043           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
1044           if (rlabel) {ierr = DMLabelSetValue(rlabel, faces[0], out->trifacemarkerlist[f]);CHKERRQ(ierr);}
1045           ierr = DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
1046         }
1047       }
1048     }
1049     ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr);
1050   }
1051   ierr = PLCDestroy(&in);CHKERRQ(ierr);
1052   ierr = PLCDestroy(&out);CHKERRQ(ierr);
1053   PetscFunctionReturn(0);
1054 }
1055 #endif /* PETSC_HAVE_CTETGEN */
1056 
1057 /*@C
1058   DMPlexGenerate - Generates a mesh.
1059 
1060   Not Collective
1061 
1062   Input Parameters:
1063 + boundary - The DMPlex boundary object
1064 . name - The mesh generation package name
1065 - interpolate - Flag to create intermediate mesh elements
1066 
1067   Output Parameter:
1068 . mesh - The DMPlex object
1069 
1070   Level: intermediate
1071 
1072 .keywords: mesh, elements
1073 .seealso: DMPlexCreate(), DMRefine()
1074 @*/
1075 PetscErrorCode DMPlexGenerate(DM boundary, const char name[], PetscBool interpolate, DM *mesh)
1076 {
1077   PetscInt       dim;
1078   char           genname[1024];
1079   PetscBool      isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;
1080   PetscErrorCode ierr;
1081 
1082   PetscFunctionBegin;
1083   PetscValidHeaderSpecific(boundary, DM_CLASSID, 1);
1084   PetscValidLogicalCollectiveBool(boundary, interpolate, 2);
1085   ierr = DMGetDimension(boundary, &dim);CHKERRQ(ierr);
1086   ierr = PetscOptionsGetString(((PetscObject) boundary)->options,((PetscObject) boundary)->prefix, "-dm_plex_generator", genname, 1024, &flg);CHKERRQ(ierr);
1087   if (flg) name = genname;
1088   if (name) {
1089     ierr = PetscStrcmp(name, "triangle", &isTriangle);CHKERRQ(ierr);
1090     ierr = PetscStrcmp(name, "tetgen",   &isTetgen);CHKERRQ(ierr);
1091     ierr = PetscStrcmp(name, "ctetgen",  &isCTetgen);CHKERRQ(ierr);
1092   }
1093   switch (dim) {
1094   case 1:
1095     if (!name || isTriangle) {
1096 #if defined(PETSC_HAVE_TRIANGLE)
1097       ierr = DMPlexGenerate_Triangle(boundary, interpolate, mesh);CHKERRQ(ierr);
1098 #else
1099       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation needs external package support.\nPlease reconfigure with --download-triangle.");
1100 #endif
1101     } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
1102     break;
1103   case 2:
1104     if (!name) {
1105 #if defined(PETSC_HAVE_CTETGEN)
1106       ierr = DMPlexGenerate_CTetgen(boundary, interpolate, mesh);CHKERRQ(ierr);
1107 #elif defined(PETSC_HAVE_TETGEN)
1108       ierr = DMPlexGenerate_Tetgen(boundary, interpolate, mesh);CHKERRQ(ierr);
1109 #else
1110       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "External package CTetgen or Tetgen needed.\nPlease reconfigure with '--download-ctetgen' or '--with-clanguage=cxx --download-tetgen'.");
1111 #endif
1112     } else if (isCTetgen) {
1113 #if defined(PETSC_HAVE_CTETGEN)
1114       ierr = DMPlexGenerate_CTetgen(boundary, interpolate, mesh);CHKERRQ(ierr);
1115 #else
1116       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
1117 #endif
1118     } else if (isTetgen) {
1119 #if defined(PETSC_HAVE_TETGEN)
1120       ierr = DMPlexGenerate_Tetgen(boundary, interpolate, mesh);CHKERRQ(ierr);
1121 #else
1122       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --download-tetgen.");
1123 #endif
1124     } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
1125     break;
1126   default:
1127     SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation for a dimension %d boundary is not supported.", dim);
1128   }
1129   PetscFunctionReturn(0);
1130 }
1131 
1132 static PetscErrorCode DMRefine_Plex_Label(DM dm, DMLabel adaptLabel, PetscInt cStart, PetscInt cEnd, PetscReal maxVolumes[])
1133 {
1134   PetscInt       dim, c;
1135   PetscReal      refRatio;
1136   PetscErrorCode ierr;
1137 
1138   PetscFunctionBegin;
1139   ierr     = DMGetDimension(dm,&dim);CHKERRQ(ierr);
1140   refRatio = (PetscReal) ((PetscInt) 1 << dim);
1141   for (c = cStart; c < cEnd; c++) {
1142     PetscReal vol;
1143     PetscInt  i, closureSize = 0;
1144     PetscInt  *closure = NULL;
1145     PetscBool anyRefine  = PETSC_FALSE;
1146     PetscBool anyCoarsen = PETSC_FALSE;
1147     PetscBool anyKeep    = PETSC_FALSE;
1148 
1149     ierr = DMPlexComputeCellGeometryFVM(dm,c,&vol,NULL,NULL);CHKERRQ(ierr);
1150     maxVolumes[c - cStart] = vol;
1151     ierr = DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
1152     for (i = 0; i < closureSize; i++) {
1153       PetscInt point = closure[2 * i], refFlag;
1154 
1155       ierr = DMLabelGetValue(adaptLabel,point,&refFlag);CHKERRQ(ierr);
1156       switch (refFlag) {
1157       case DM_ADAPT_REFINE:
1158         anyRefine = PETSC_TRUE;
1159         break;
1160       case DM_ADAPT_COARSEN:
1161         anyCoarsen = PETSC_TRUE;
1162         break;
1163       case DM_ADAPT_KEEP:
1164         anyKeep = PETSC_TRUE;
1165         break;
1166       case DM_ADAPT_DETERMINE:
1167         break;
1168       default:
1169         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"DMPlex does not support refinement flag %D\n",refFlag);
1170         break;
1171       }
1172       if (anyRefine) break;
1173     }
1174     ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
1175     if (anyRefine) {
1176       maxVolumes[c - cStart] = vol / refRatio;
1177     } else if (anyKeep) {
1178       maxVolumes[c - cStart] = vol;
1179     } else if (anyCoarsen) {
1180       maxVolumes[c - cStart] = vol * refRatio;
1181     }
1182   }
1183   PetscFunctionReturn(0);
1184 }
1185 
1186 static PetscErrorCode DMRefine_Plex_Internal(DM dm, MPI_Comm comm, DMLabel adaptLabel, DM *dmRefined)
1187 {
1188   PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *);
1189   PetscReal        refinementLimit;
1190   PetscInt         dim, cStart, cEnd;
1191   char             genname[1024], *name = NULL;
1192   PetscBool        isUniform, isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg, localized;
1193   PetscErrorCode   ierr;
1194 
1195   PetscFunctionBegin;
1196   ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr);
1197   ierr = DMPlexGetRefinementUniform(dm, &isUniform);CHKERRQ(ierr);
1198   if (isUniform) {
1199     CellRefiner cellRefiner;
1200 
1201     ierr = DMPlexGetCellRefiner_Internal(dm, &cellRefiner);CHKERRQ(ierr);
1202     ierr = DMPlexRefineUniform_Internal(dm, cellRefiner, dmRefined);CHKERRQ(ierr);
1203     ierr = DMCopyBoundary(dm, *dmRefined);CHKERRQ(ierr);
1204     if (localized) {
1205       ierr = DMLocalizeCoordinates(*dmRefined);CHKERRQ(ierr);
1206     }
1207     PetscFunctionReturn(0);
1208   }
1209   ierr = DMPlexGetRefinementLimit(dm, &refinementLimit);CHKERRQ(ierr);
1210   ierr = DMPlexGetRefinementFunction(dm, &refinementFunc);CHKERRQ(ierr);
1211   if (refinementLimit == 0.0 && !refinementFunc) PetscFunctionReturn(0);
1212   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1213   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1214   ierr = PetscOptionsGetString(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_generator", genname, 1024, &flg);CHKERRQ(ierr);
1215   if (flg) name = genname;
1216   if (name) {
1217     ierr = PetscStrcmp(name, "triangle", &isTriangle);CHKERRQ(ierr);
1218     ierr = PetscStrcmp(name, "tetgen",   &isTetgen);CHKERRQ(ierr);
1219     ierr = PetscStrcmp(name, "ctetgen",  &isCTetgen);CHKERRQ(ierr);
1220   }
1221   switch (dim) {
1222   case 2:
1223     if (!name || isTriangle) {
1224 #if defined(PETSC_HAVE_TRIANGLE)
1225       PetscReal *maxVolumes;
1226       PetscInt  c;
1227 
1228       ierr = PetscMalloc1(cEnd - cStart, &maxVolumes);CHKERRQ(ierr);
1229       if (adaptLabel) {
1230         ierr = DMRefine_Plex_Label(dm,adaptLabel,cStart,cEnd,maxVolumes);CHKERRQ(ierr);
1231       } else if (refinementFunc) {
1232         for (c = cStart; c < cEnd; ++c) {
1233           PetscReal vol, centroid[3];
1234           PetscReal maxVol;
1235 
1236           ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);CHKERRQ(ierr);
1237           ierr = (*refinementFunc)(centroid, &maxVol);CHKERRQ(ierr);
1238           maxVolumes[c - cStart] = (double) maxVol;
1239         }
1240       } else {
1241         for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1242       }
1243 #if !defined(PETSC_USE_REAL_DOUBLE)
1244       {
1245         double *mvols;
1246         ierr = PetscMalloc1(cEnd - cStart,&mvols);CHKERRQ(ierr);
1247         for (c = 0; c < cEnd-cStart; ++c) mvols[c] = (double)maxVolumes[c];
1248         ierr = DMPlexRefine_Triangle(dm, mvols, dmRefined);CHKERRQ(ierr);
1249         ierr = PetscFree(mvols);CHKERRQ(ierr);
1250       }
1251 #else
1252       ierr = DMPlexRefine_Triangle(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
1253 #endif
1254       ierr = PetscFree(maxVolumes);CHKERRQ(ierr);
1255 #else
1256       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement needs external package support.\nPlease reconfigure with --download-triangle.");
1257 #endif
1258     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
1259     break;
1260   case 3:
1261     if (!name || isCTetgen || isTetgen) {
1262       PetscReal *maxVolumes;
1263       PetscInt   c;
1264 
1265       ierr = PetscMalloc1(cEnd - cStart, &maxVolumes);CHKERRQ(ierr);
1266       if (adaptLabel) {
1267         ierr = DMRefine_Plex_Label(dm,adaptLabel,cStart,cEnd,maxVolumes);CHKERRQ(ierr);
1268       } else if (refinementFunc) {
1269         for (c = cStart; c < cEnd; ++c) {
1270           PetscReal vol, centroid[3];
1271 
1272           ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);CHKERRQ(ierr);
1273           ierr = (*refinementFunc)(centroid, &maxVolumes[c-cStart]);CHKERRQ(ierr);
1274         }
1275       } else {
1276         for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1277       }
1278       if (!name) {
1279 #if defined(PETSC_HAVE_CTETGEN)
1280         ierr = DMPlexRefine_CTetgen(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
1281 #elif defined(PETSC_HAVE_TETGEN)
1282         ierr = DMPlexRefine_Tetgen(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
1283 #else
1284         SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "External package CTetgen or Tetgen needed.\nPlease reconfigure with '--download-ctetgen' or '--with-clanguage=cxx --download-tetgen'.");
1285 #endif
1286       } else if (isCTetgen) {
1287 #if defined(PETSC_HAVE_CTETGEN)
1288         ierr = DMPlexRefine_CTetgen(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
1289 #else
1290         SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
1291 #endif
1292       } else {
1293 #if defined(PETSC_HAVE_TETGEN)
1294         ierr = DMPlexRefine_Tetgen(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
1295 #else
1296         SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --download-tetgen.");
1297 #endif
1298       }
1299       ierr = PetscFree(maxVolumes);CHKERRQ(ierr);
1300     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
1301     break;
1302   default:
1303     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %d is not supported.", dim);
1304   }
1305   ierr = DMCopyBoundary(dm, *dmRefined);CHKERRQ(ierr);
1306   if (localized) {
1307     ierr = DMLocalizeCoordinates(*dmRefined);CHKERRQ(ierr);
1308   }
1309   PetscFunctionReturn(0);
1310 }
1311 
1312 PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined)
1313 {
1314   PetscErrorCode ierr;
1315 
1316   PetscFunctionBegin;
1317   ierr = DMRefine_Plex_Internal(dm,comm,NULL,dmRefined);CHKERRQ(ierr);
1318   PetscFunctionReturn(0);
1319 }
1320 
1321 PetscErrorCode DMAdaptLabel_Plex(DM dm, DMLabel adaptLabel, DM *dmAdapted)
1322 {
1323   PetscInt       defFlag, minFlag, maxFlag, numFlags, i;
1324   const PetscInt *flags;
1325   IS             flagIS;
1326   PetscErrorCode ierr;
1327 
1328   PetscFunctionBegin;
1329   ierr = DMLabelGetDefaultValue(adaptLabel,&defFlag);CHKERRQ(ierr);
1330   minFlag = defFlag;
1331   maxFlag = defFlag;
1332   ierr = DMLabelGetValueIS(adaptLabel,&flagIS);CHKERRQ(ierr);
1333   ierr = ISGetLocalSize(flagIS,&numFlags);CHKERRQ(ierr);
1334   ierr = ISGetIndices(flagIS,&flags);CHKERRQ(ierr);
1335   for (i = 0; i < numFlags; i++) {
1336     PetscInt flag = flags[i];
1337 
1338     minFlag = PetscMin(minFlag,flag);
1339     maxFlag = PetscMax(maxFlag,flag);
1340   }
1341   ierr = ISRestoreIndices(flagIS,&flags);CHKERRQ(ierr);
1342   ierr = ISDestroy(&flagIS);CHKERRQ(ierr);
1343   {
1344     PetscInt minMaxFlag[2], minMaxFlagGlobal[2];
1345 
1346     minMaxFlag[0] = minFlag;
1347     minMaxFlag[1] = -maxFlag;
1348     ierr = MPI_Allreduce(minMaxFlag,minMaxFlagGlobal,2,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
1349     minFlag = minMaxFlagGlobal[0];
1350     maxFlag = -minMaxFlagGlobal[1];
1351   }
1352   if (minFlag == maxFlag) {
1353     switch (minFlag) {
1354     case DM_ADAPT_DETERMINE:
1355       *dmAdapted = NULL;
1356       break;
1357     case DM_ADAPT_REFINE:
1358       ierr = DMPlexSetRefinementUniform(dm,PETSC_TRUE);CHKERRQ(ierr);
1359       ierr = DMRefine(dm,MPI_COMM_NULL,dmAdapted);CHKERRQ(ierr);
1360       break;
1361     case DM_ADAPT_COARSEN:
1362       ierr = DMCoarsen(dm,MPI_COMM_NULL,dmAdapted);CHKERRQ(ierr);
1363       break;
1364     default:
1365       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"DMPlex does not support refinement flag %D\n",minFlag);
1366       break;
1367     }
1368   } else {
1369     ierr = DMPlexSetRefinementUniform(dm,PETSC_FALSE);CHKERRQ(ierr);
1370     ierr = DMRefine_Plex_Internal(dm,MPI_COMM_NULL,adaptLabel,dmAdapted);CHKERRQ(ierr);
1371   }
1372   PetscFunctionReturn(0);
1373 }
1374 
1375 PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[])
1376 {
1377   DM             cdm = dm;
1378   PetscInt       r;
1379   PetscBool      isUniform, localized;
1380   PetscErrorCode ierr;
1381 
1382   PetscFunctionBegin;
1383   ierr = DMPlexGetRefinementUniform(dm, &isUniform);CHKERRQ(ierr);
1384   ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr);
1385   if (!isUniform) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Non-uniform refinement is incompatible with the hierarchy");
1386   for (r = 0; r < nlevels; ++r) {
1387     CellRefiner cellRefiner;
1388 
1389     ierr = DMPlexGetCellRefiner_Internal(cdm, &cellRefiner);CHKERRQ(ierr);
1390     ierr = DMPlexRefineUniform_Internal(cdm, cellRefiner, &dmRefined[r]);CHKERRQ(ierr);
1391     ierr = DMCopyBoundary(cdm, dmRefined[r]);CHKERRQ(ierr);
1392     if (localized) {
1393       ierr = DMLocalizeCoordinates(dmRefined[r]);CHKERRQ(ierr);
1394     }
1395     ierr = DMSetCoarseDM(dmRefined[r], cdm);CHKERRQ(ierr);
1396     ierr = DMPlexSetRegularRefinement(dmRefined[r], PETSC_TRUE);CHKERRQ(ierr);
1397     cdm  = dmRefined[r];
1398   }
1399   PetscFunctionReturn(0);
1400 }
1401