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