xref: /petsc/src/dm/impls/plex/plexgenerate.c (revision 4446f3b4cb0d5758d26ec359716e4cff36668235) 
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   struct triangulateio in;
210   struct triangulateio out;
211   PetscInt             vStart, vEnd, v, eStart, eEnd, e;
212   PetscMPIInt          rank;
213   PetscErrorCode       ierr;
214 
215   PetscFunctionBegin;
216   ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr);
217   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
218   ierr = InitInput_Triangle(&in);CHKERRQ(ierr);
219   ierr = InitOutput_Triangle(&out);CHKERRQ(ierr);
220   ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr);
221 
222   in.numberofpoints = vEnd - vStart;
223   if (in.numberofpoints > 0) {
224     PetscSection coordSection;
225     Vec          coordinates;
226     PetscScalar *array;
227 
228     ierr = PetscMalloc1(in.numberofpoints*dim, &in.pointlist);CHKERRQ(ierr);
229     ierr = PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);CHKERRQ(ierr);
230     ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr);
231     ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr);
232     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
233     for (v = vStart; v < vEnd; ++v) {
234       const PetscInt idx = v - vStart;
235       PetscInt       off, d;
236 
237       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
238       for (d = 0; d < dim; ++d) {
239         in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
240       }
241       ierr = DMPlexGetLabelValue(boundary, "marker", v, &in.pointmarkerlist[idx]);CHKERRQ(ierr);
242     }
243     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
244   }
245   ierr  = DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);CHKERRQ(ierr);
246   in.numberofsegments = eEnd - eStart;
247   if (in.numberofsegments > 0) {
248     ierr = PetscMalloc1(in.numberofsegments*2, &in.segmentlist);CHKERRQ(ierr);
249     ierr = PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist);CHKERRQ(ierr);
250     for (e = eStart; e < eEnd; ++e) {
251       const PetscInt  idx = e - eStart;
252       const PetscInt *cone;
253 
254       ierr = DMPlexGetCone(boundary, e, &cone);CHKERRQ(ierr);
255 
256       in.segmentlist[idx*2+0] = cone[0] - vStart;
257       in.segmentlist[idx*2+1] = cone[1] - vStart;
258 
259       ierr = DMPlexGetLabelValue(boundary, "marker", e, &in.segmentmarkerlist[idx]);CHKERRQ(ierr);
260     }
261   }
262 #if 0 /* Do not currently support holes */
263   PetscReal *holeCoords;
264   PetscInt   h, d;
265 
266   ierr = DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);CHKERRQ(ierr);
267   if (in.numberofholes > 0) {
268     ierr = PetscMalloc1(in.numberofholes*dim, &in.holelist);CHKERRQ(ierr);
269     for (h = 0; h < in.numberofholes; ++h) {
270       for (d = 0; d < dim; ++d) {
271         in.holelist[h*dim+d] = holeCoords[h*dim+d];
272       }
273     }
274   }
275 #endif
276   if (!rank) {
277     char args[32];
278 
279     /* Take away 'Q' for verbose output */
280     ierr = PetscStrcpy(args, "pqezQ");CHKERRQ(ierr);
281     if (createConvexHull)   {ierr = PetscStrcat(args, "c");CHKERRQ(ierr);}
282     if (constrained)        {ierr = PetscStrcpy(args, "zepDQ");CHKERRQ(ierr);}
283     if (mesh->triangleOpts) {triangulate(mesh->triangleOpts, &in, &out, NULL);}
284     else                    {triangulate(args, &in, &out, NULL);}
285   }
286   ierr = PetscFree(in.pointlist);CHKERRQ(ierr);
287   ierr = PetscFree(in.pointmarkerlist);CHKERRQ(ierr);
288   ierr = PetscFree(in.segmentlist);CHKERRQ(ierr);
289   ierr = PetscFree(in.segmentmarkerlist);CHKERRQ(ierr);
290   ierr = PetscFree(in.holelist);CHKERRQ(ierr);
291 
292   {
293     const PetscInt numCorners  = 3;
294     const PetscInt numCells    = out.numberoftriangles;
295     const PetscInt numVertices = out.numberofpoints;
296     const int     *cells      = out.trianglelist;
297     const double  *meshCoords = out.pointlist;
298 
299     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr);
300     /* Set labels */
301     for (v = 0; v < numVertices; ++v) {
302       if (out.pointmarkerlist[v]) {
303         ierr = DMPlexSetLabelValue(*dm, "marker", v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);
304       }
305     }
306     if (interpolate) {
307       for (e = 0; e < out.numberofedges; e++) {
308         if (out.edgemarkerlist[e]) {
309           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
310           const PetscInt *edges;
311           PetscInt        numEdges;
312 
313           ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
314           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
315           ierr = DMPlexSetLabelValue(*dm, "marker", edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);
316           ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
317         }
318       }
319     }
320     ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
321   }
322 #if 0 /* Do not currently support holes */
323   ierr = DMPlexCopyHoles(*dm, boundary);CHKERRQ(ierr);
324 #endif
325   ierr = FiniOutput_Triangle(&out);CHKERRQ(ierr);
326   PetscFunctionReturn(0);
327 }
328 
329 #undef __FUNCT__
330 #define __FUNCT__ "DMPlexRefine_Triangle"
331 PetscErrorCode DMPlexRefine_Triangle(DM dm, double *maxVolumes, DM *dmRefined)
332 {
333   MPI_Comm             comm;
334   PetscInt             dim  = 2;
335   struct triangulateio in;
336   struct triangulateio out;
337   PetscInt             vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
338   PetscMPIInt          rank;
339   PetscErrorCode       ierr;
340 
341   PetscFunctionBegin;
342   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
343   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
344   ierr = InitInput_Triangle(&in);CHKERRQ(ierr);
345   ierr = InitOutput_Triangle(&out);CHKERRQ(ierr);
346   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
347   ierr = MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
348   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
349 
350   in.numberofpoints = vEnd - vStart;
351   if (in.numberofpoints > 0) {
352     PetscSection coordSection;
353     Vec          coordinates;
354     PetscScalar *array;
355 
356     ierr = PetscMalloc1(in.numberofpoints*dim, &in.pointlist);CHKERRQ(ierr);
357     ierr = PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);CHKERRQ(ierr);
358     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
359     ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
360     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
361     for (v = vStart; v < vEnd; ++v) {
362       const PetscInt idx = v - vStart;
363       PetscInt       off, d;
364 
365       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
366       for (d = 0; d < dim; ++d) {
367         in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
368       }
369       ierr = DMPlexGetLabelValue(dm, "marker", v, &in.pointmarkerlist[idx]);CHKERRQ(ierr);
370     }
371     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
372   }
373   ierr  = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
374 
375   in.numberofcorners   = 3;
376   in.numberoftriangles = cEnd - cStart;
377 
378   in.trianglearealist  = (double*) maxVolumes;
379   if (in.numberoftriangles > 0) {
380     ierr = PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);CHKERRQ(ierr);
381     for (c = cStart; c < cEnd; ++c) {
382       const PetscInt idx      = c - cStart;
383       PetscInt      *closure = NULL;
384       PetscInt       closureSize;
385 
386       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
387       if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize);
388       for (v = 0; v < 3; ++v) {
389         in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart;
390       }
391       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
392     }
393   }
394   /* TODO: Segment markers are missing on input */
395 #if 0 /* Do not currently support holes */
396   PetscReal *holeCoords;
397   PetscInt   h, d;
398 
399   ierr = DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);CHKERRQ(ierr);
400   if (in.numberofholes > 0) {
401     ierr = PetscMalloc1(in.numberofholes*dim, &in.holelist);CHKERRQ(ierr);
402     for (h = 0; h < in.numberofholes; ++h) {
403       for (d = 0; d < dim; ++d) {
404         in.holelist[h*dim+d] = holeCoords[h*dim+d];
405       }
406     }
407   }
408 #endif
409   if (!rank) {
410     char args[32];
411 
412     /* Take away 'Q' for verbose output */
413     ierr = PetscStrcpy(args, "pqezQra");CHKERRQ(ierr);
414     triangulate(args, &in, &out, NULL);
415   }
416   ierr = PetscFree(in.pointlist);CHKERRQ(ierr);
417   ierr = PetscFree(in.pointmarkerlist);CHKERRQ(ierr);
418   ierr = PetscFree(in.segmentlist);CHKERRQ(ierr);
419   ierr = PetscFree(in.segmentmarkerlist);CHKERRQ(ierr);
420   ierr = PetscFree(in.trianglelist);CHKERRQ(ierr);
421 
422   {
423     const PetscInt numCorners  = 3;
424     const PetscInt numCells    = out.numberoftriangles;
425     const PetscInt numVertices = out.numberofpoints;
426     const int     *cells      = out.trianglelist;
427     const double  *meshCoords = out.pointlist;
428     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
429 
430     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr);
431     /* Set labels */
432     for (v = 0; v < numVertices; ++v) {
433       if (out.pointmarkerlist[v]) {
434         ierr = DMPlexSetLabelValue(*dmRefined, "marker", v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);
435       }
436     }
437     if (interpolate) {
438       PetscInt e;
439 
440       for (e = 0; e < out.numberofedges; e++) {
441         if (out.edgemarkerlist[e]) {
442           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
443           const PetscInt *edges;
444           PetscInt        numEdges;
445 
446           ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
447           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
448           ierr = DMPlexSetLabelValue(*dmRefined, "marker", edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);
449           ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
450         }
451       }
452     }
453     ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr);
454   }
455 #if 0 /* Do not currently support holes */
456   ierr = DMPlexCopyHoles(*dm, boundary);CHKERRQ(ierr);
457 #endif
458   ierr = FiniOutput_Triangle(&out);CHKERRQ(ierr);
459   PetscFunctionReturn(0);
460 }
461 #endif /* PETSC_HAVE_TRIANGLE */
462 
463 #if defined(PETSC_HAVE_TETGEN)
464 #include <tetgen.h>
465 #undef __FUNCT__
466 #define __FUNCT__ "DMPlexGenerate_Tetgen"
467 PetscErrorCode DMPlexGenerate_Tetgen(DM boundary, PetscBool interpolate, DM *dm)
468 {
469   MPI_Comm       comm;
470   DM_Plex       *mesh = (DM_Plex *) boundary->data;
471   const PetscInt dim  = 3;
472   ::tetgenio     in;
473   ::tetgenio     out;
474   PetscInt       vStart, vEnd, v, fStart, fEnd, f;
475   PetscMPIInt    rank;
476   PetscErrorCode ierr;
477 
478   PetscFunctionBegin;
479   ierr              = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr);
480   ierr              = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
481   ierr              = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr);
482   in.numberofpoints = vEnd - vStart;
483   if (in.numberofpoints > 0) {
484     PetscSection coordSection;
485     Vec          coordinates;
486     PetscScalar *array;
487 
488     in.pointlist       = new double[in.numberofpoints*dim];
489     in.pointmarkerlist = new int[in.numberofpoints];
490 
491     ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr);
492     ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr);
493     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
494     for (v = vStart; v < vEnd; ++v) {
495       const PetscInt idx = v - vStart;
496       PetscInt       off, d;
497 
498       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
499       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
500       ierr = DMPlexGetLabelValue(boundary, "marker", v, &in.pointmarkerlist[idx]);CHKERRQ(ierr);
501     }
502     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
503   }
504   ierr  = DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);CHKERRQ(ierr);
505 
506   in.numberoffacets = fEnd - fStart;
507   if (in.numberoffacets > 0) {
508     in.facetlist       = new tetgenio::facet[in.numberoffacets];
509     in.facetmarkerlist = new int[in.numberoffacets];
510     for (f = fStart; f < fEnd; ++f) {
511       const PetscInt idx     = f - fStart;
512       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v;
513 
514       in.facetlist[idx].numberofpolygons = 1;
515       in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
516       in.facetlist[idx].numberofholes    = 0;
517       in.facetlist[idx].holelist         = NULL;
518 
519       ierr = DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
520       for (p = 0; p < numPoints*2; p += 2) {
521         const PetscInt point = points[p];
522         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
523       }
524 
525       tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
526       poly->numberofvertices = numVertices;
527       poly->vertexlist       = new int[poly->numberofvertices];
528       for (v = 0; v < numVertices; ++v) {
529         const PetscInt vIdx = points[v] - vStart;
530         poly->vertexlist[v] = vIdx;
531       }
532       ierr = DMPlexGetLabelValue(boundary, "marker", f, &in.facetmarkerlist[idx]);CHKERRQ(ierr);
533       ierr = DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
534     }
535   }
536   if (!rank) {
537     char args[32];
538 
539     /* Take away 'Q' for verbose output */
540     ierr = PetscStrcpy(args, "pqezQ");CHKERRQ(ierr);
541     if (mesh->tetgenOpts) {::tetrahedralize(mesh->tetgenOpts, &in, &out);}
542     else                  {::tetrahedralize(args, &in, &out);}
543   }
544   {
545     const PetscInt numCorners  = 4;
546     const PetscInt numCells    = out.numberoftetrahedra;
547     const PetscInt numVertices = out.numberofpoints;
548     const double   *meshCoords = out.pointlist;
549     int            *cells      = out.tetrahedronlist;
550 
551     ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr);
552     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr);
553     /* Set labels */
554     for (v = 0; v < numVertices; ++v) {
555       if (out.pointmarkerlist[v]) {
556         ierr = DMPlexSetLabelValue(*dm, "marker", v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);
557       }
558     }
559     if (interpolate) {
560       PetscInt e;
561 
562       for (e = 0; e < out.numberofedges; e++) {
563         if (out.edgemarkerlist[e]) {
564           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
565           const PetscInt *edges;
566           PetscInt        numEdges;
567 
568           ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
569           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
570           ierr = DMPlexSetLabelValue(*dm, "marker", edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);
571           ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
572         }
573       }
574       for (f = 0; f < out.numberoftrifaces; f++) {
575         if (out.trifacemarkerlist[f]) {
576           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
577           const PetscInt *faces;
578           PetscInt        numFaces;
579 
580           ierr = DMPlexGetJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
581           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
582           ierr = DMPlexSetLabelValue(*dm, "marker", faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);
583           ierr = DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
584         }
585       }
586     }
587     ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
588   }
589   PetscFunctionReturn(0);
590 }
591 
592 #undef __FUNCT__
593 #define __FUNCT__ "DMPlexRefine_Tetgen"
594 PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
595 {
596   MPI_Comm       comm;
597   const PetscInt dim  = 3;
598   ::tetgenio     in;
599   ::tetgenio     out;
600   PetscInt       vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
601   PetscMPIInt    rank;
602   PetscErrorCode ierr;
603 
604   PetscFunctionBegin;
605   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
606   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
607   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
608   ierr = MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
609   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
610 
611   in.numberofpoints = vEnd - vStart;
612   if (in.numberofpoints > 0) {
613     PetscSection coordSection;
614     Vec          coordinates;
615     PetscScalar *array;
616 
617     in.pointlist       = new double[in.numberofpoints*dim];
618     in.pointmarkerlist = new int[in.numberofpoints];
619 
620     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
621     ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
622     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
623     for (v = vStart; v < vEnd; ++v) {
624       const PetscInt idx = v - vStart;
625       PetscInt       off, d;
626 
627       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
628       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
629       ierr = DMPlexGetLabelValue(dm, "marker", v, &in.pointmarkerlist[idx]);CHKERRQ(ierr);
630     }
631     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
632   }
633   ierr  = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
634 
635   in.numberofcorners       = 4;
636   in.numberoftetrahedra    = cEnd - cStart;
637   in.tetrahedronvolumelist = (double*) maxVolumes;
638   if (in.numberoftetrahedra > 0) {
639     in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
640     for (c = cStart; c < cEnd; ++c) {
641       const PetscInt idx      = c - cStart;
642       PetscInt      *closure = NULL;
643       PetscInt       closureSize;
644 
645       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
646       if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
647       for (v = 0; v < 4; ++v) {
648         in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
649       }
650       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
651     }
652   }
653   /* TODO: Put in boundary faces with markers */
654   if (!rank) {
655     char args[32];
656 
657     /* Take away 'Q' for verbose output */
658     /*ierr = PetscStrcpy(args, "qezQra");CHKERRQ(ierr); */
659     ierr = PetscStrcpy(args, "qezraVVVV");CHKERRQ(ierr);
660     ::tetrahedralize(args, &in, &out);
661   }
662   in.tetrahedronvolumelist = NULL;
663 
664   {
665     const PetscInt numCorners  = 4;
666     const PetscInt numCells    = out.numberoftetrahedra;
667     const PetscInt numVertices = out.numberofpoints;
668     const double   *meshCoords = out.pointlist;
669     int            *cells      = out.tetrahedronlist;
670 
671     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
672 
673     ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr);
674     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr);
675     /* Set labels */
676     for (v = 0; v < numVertices; ++v) {
677       if (out.pointmarkerlist[v]) {
678         ierr = DMPlexSetLabelValue(*dmRefined, "marker", v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);
679       }
680     }
681     if (interpolate) {
682       PetscInt e, f;
683 
684       for (e = 0; e < out.numberofedges; e++) {
685         if (out.edgemarkerlist[e]) {
686           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
687           const PetscInt *edges;
688           PetscInt        numEdges;
689 
690           ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
691           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
692           ierr = DMPlexSetLabelValue(*dmRefined, "marker", edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);
693           ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
694         }
695       }
696       for (f = 0; f < out.numberoftrifaces; f++) {
697         if (out.trifacemarkerlist[f]) {
698           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
699           const PetscInt *faces;
700           PetscInt        numFaces;
701 
702           ierr = DMPlexGetJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
703           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
704           ierr = DMPlexSetLabelValue(*dmRefined, "marker", faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);
705           ierr = DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
706         }
707       }
708     }
709     ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr);
710   }
711   PetscFunctionReturn(0);
712 }
713 #endif /* PETSC_HAVE_TETGEN */
714 
715 #if defined(PETSC_HAVE_CTETGEN)
716 #include <ctetgen.h>
717 
718 #undef __FUNCT__
719 #define __FUNCT__ "DMPlexGenerate_CTetgen"
720 PetscErrorCode DMPlexGenerate_CTetgen(DM boundary, PetscBool interpolate, DM *dm)
721 {
722   MPI_Comm       comm;
723   const PetscInt dim  = 3;
724   PLC           *in, *out;
725   PetscInt       verbose = 0, vStart, vEnd, v, fStart, fEnd, f;
726   PetscMPIInt    rank;
727   PetscErrorCode ierr;
728 
729   PetscFunctionBegin;
730   ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr);
731   ierr = PetscOptionsGetInt(((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);CHKERRQ(ierr);
732   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
733   ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr);
734   ierr = PLCCreate(&in);CHKERRQ(ierr);
735   ierr = PLCCreate(&out);CHKERRQ(ierr);
736 
737   in->numberofpoints = vEnd - vStart;
738   if (in->numberofpoints > 0) {
739     PetscSection coordSection;
740     Vec          coordinates;
741     PetscScalar *array;
742 
743     ierr = PetscMalloc1(in->numberofpoints*dim, &in->pointlist);CHKERRQ(ierr);
744     ierr = PetscMalloc1(in->numberofpoints,       &in->pointmarkerlist);CHKERRQ(ierr);
745     ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr);
746     ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr);
747     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
748     for (v = vStart; v < vEnd; ++v) {
749       const PetscInt idx = v - vStart;
750       PetscInt       off, d, m;
751 
752       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
753       for (d = 0; d < dim; ++d) {
754         in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
755       }
756       ierr = DMPlexGetLabelValue(boundary, "marker", v, &m);CHKERRQ(ierr);
757 
758       in->pointmarkerlist[idx] = (int) m;
759     }
760     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
761   }
762   ierr  = DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);CHKERRQ(ierr);
763 
764   in->numberoffacets = fEnd - fStart;
765   if (in->numberoffacets > 0) {
766     ierr = PetscMalloc1(in->numberoffacets, &in->facetlist);CHKERRQ(ierr);
767     ierr = PetscMalloc1(in->numberoffacets,   &in->facetmarkerlist);CHKERRQ(ierr);
768     for (f = fStart; f < fEnd; ++f) {
769       const PetscInt idx     = f - fStart;
770       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v, m;
771       polygon       *poly;
772 
773       in->facetlist[idx].numberofpolygons = 1;
774 
775       ierr = PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist);CHKERRQ(ierr);
776 
777       in->facetlist[idx].numberofholes    = 0;
778       in->facetlist[idx].holelist         = NULL;
779 
780       ierr = DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
781       for (p = 0; p < numPoints*2; p += 2) {
782         const PetscInt point = points[p];
783         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
784       }
785 
786       poly                   = in->facetlist[idx].polygonlist;
787       poly->numberofvertices = numVertices;
788       ierr                   = PetscMalloc1(poly->numberofvertices, &poly->vertexlist);CHKERRQ(ierr);
789       for (v = 0; v < numVertices; ++v) {
790         const PetscInt vIdx = points[v] - vStart;
791         poly->vertexlist[v] = vIdx;
792       }
793       ierr                     = DMPlexGetLabelValue(boundary, "marker", f, &m);CHKERRQ(ierr);
794       in->facetmarkerlist[idx] = (int) m;
795       ierr                     = DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
796     }
797   }
798   if (!rank) {
799     TetGenOpts t;
800 
801     ierr        = TetGenOptsInitialize(&t);CHKERRQ(ierr);
802     t.in        = boundary; /* Should go away */
803     t.plc       = 1;
804     t.quality   = 1;
805     t.edgesout  = 1;
806     t.zeroindex = 1;
807     t.quiet     = 1;
808     t.verbose   = verbose;
809     ierr        = TetGenCheckOpts(&t);CHKERRQ(ierr);
810     ierr        = TetGenTetrahedralize(&t, in, out);CHKERRQ(ierr);
811   }
812   {
813     const PetscInt numCorners  = 4;
814     const PetscInt numCells    = out->numberoftetrahedra;
815     const PetscInt numVertices = out->numberofpoints;
816     const double   *meshCoords = out->pointlist;
817     int            *cells      = out->tetrahedronlist;
818 
819     ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr);
820     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr);
821     /* Set labels */
822     for (v = 0; v < numVertices; ++v) {
823       if (out->pointmarkerlist[v]) {
824         ierr = DMPlexSetLabelValue(*dm, "marker", v+numCells, out->pointmarkerlist[v]);CHKERRQ(ierr);
825       }
826     }
827     if (interpolate) {
828       PetscInt e;
829 
830       for (e = 0; e < out->numberofedges; e++) {
831         if (out->edgemarkerlist[e]) {
832           const PetscInt  vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
833           const PetscInt *edges;
834           PetscInt        numEdges;
835 
836           ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
837           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
838           ierr = DMPlexSetLabelValue(*dm, "marker", edges[0], out->edgemarkerlist[e]);CHKERRQ(ierr);
839           ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
840         }
841       }
842       for (f = 0; f < out->numberoftrifaces; f++) {
843         if (out->trifacemarkerlist[f]) {
844           const PetscInt  vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
845           const PetscInt *faces;
846           PetscInt        numFaces;
847 
848           ierr = DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
849           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
850           ierr = DMPlexSetLabelValue(*dm, "marker", faces[0], out->trifacemarkerlist[f]);CHKERRQ(ierr);
851           ierr = DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
852         }
853       }
854     }
855     ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
856   }
857 
858   ierr = PLCDestroy(&in);CHKERRQ(ierr);
859   ierr = PLCDestroy(&out);CHKERRQ(ierr);
860   PetscFunctionReturn(0);
861 }
862 
863 #undef __FUNCT__
864 #define __FUNCT__ "DMPlexRefine_CTetgen"
865 PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
866 {
867   MPI_Comm       comm;
868   const PetscInt dim  = 3;
869   PLC           *in, *out;
870   PetscInt       verbose = 0, vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
871   PetscMPIInt    rank;
872   PetscErrorCode ierr;
873 
874   PetscFunctionBegin;
875   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
876   ierr = PetscOptionsGetInt(((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);CHKERRQ(ierr);
877   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
878   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
879   ierr = MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
880   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
881   ierr = PLCCreate(&in);CHKERRQ(ierr);
882   ierr = PLCCreate(&out);CHKERRQ(ierr);
883 
884   in->numberofpoints = vEnd - vStart;
885   if (in->numberofpoints > 0) {
886     PetscSection coordSection;
887     Vec          coordinates;
888     PetscScalar *array;
889 
890     ierr = PetscMalloc1(in->numberofpoints*dim, &in->pointlist);CHKERRQ(ierr);
891     ierr = PetscMalloc1(in->numberofpoints,       &in->pointmarkerlist);CHKERRQ(ierr);
892     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
893     ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
894     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
895     for (v = vStart; v < vEnd; ++v) {
896       const PetscInt idx = v - vStart;
897       PetscInt       off, d, m;
898 
899       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
900       for (d = 0; d < dim; ++d) {
901         in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
902       }
903       ierr = DMPlexGetLabelValue(dm, "marker", v, &m);CHKERRQ(ierr);
904 
905       in->pointmarkerlist[idx] = (int) m;
906     }
907     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
908   }
909   ierr  = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
910 
911   in->numberofcorners       = 4;
912   in->numberoftetrahedra    = cEnd - cStart;
913   in->tetrahedronvolumelist = maxVolumes;
914   if (in->numberoftetrahedra > 0) {
915     ierr = PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);CHKERRQ(ierr);
916     for (c = cStart; c < cEnd; ++c) {
917       const PetscInt idx      = c - cStart;
918       PetscInt      *closure = NULL;
919       PetscInt       closureSize;
920 
921       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
922       if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
923       for (v = 0; v < 4; ++v) {
924         in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
925       }
926       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
927     }
928   }
929   if (!rank) {
930     TetGenOpts t;
931 
932     ierr = TetGenOptsInitialize(&t);CHKERRQ(ierr);
933 
934     t.in        = dm; /* Should go away */
935     t.refine    = 1;
936     t.varvolume = 1;
937     t.quality   = 1;
938     t.edgesout  = 1;
939     t.zeroindex = 1;
940     t.quiet     = 1;
941     t.verbose   = verbose; /* Change this */
942 
943     ierr = TetGenCheckOpts(&t);CHKERRQ(ierr);
944     ierr = TetGenTetrahedralize(&t, in, out);CHKERRQ(ierr);
945   }
946   {
947     const PetscInt numCorners  = 4;
948     const PetscInt numCells    = out->numberoftetrahedra;
949     const PetscInt numVertices = out->numberofpoints;
950     const double   *meshCoords = out->pointlist;
951     int            *cells      = out->tetrahedronlist;
952     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
953 
954     ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr);
955     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr);
956     /* Set labels */
957     for (v = 0; v < numVertices; ++v) {
958       if (out->pointmarkerlist[v]) {
959         ierr = DMPlexSetLabelValue(*dmRefined, "marker", v+numCells, out->pointmarkerlist[v]);CHKERRQ(ierr);
960       }
961     }
962     if (interpolate) {
963       PetscInt e, f;
964 
965       for (e = 0; e < out->numberofedges; e++) {
966         if (out->edgemarkerlist[e]) {
967           const PetscInt  vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
968           const PetscInt *edges;
969           PetscInt        numEdges;
970 
971           ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
972           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
973           ierr = DMPlexSetLabelValue(*dmRefined, "marker", edges[0], out->edgemarkerlist[e]);CHKERRQ(ierr);
974           ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
975         }
976       }
977       for (f = 0; f < out->numberoftrifaces; f++) {
978         if (out->trifacemarkerlist[f]) {
979           const PetscInt  vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
980           const PetscInt *faces;
981           PetscInt        numFaces;
982 
983           ierr = DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
984           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
985           ierr = DMPlexSetLabelValue(*dmRefined, "marker", faces[0], out->trifacemarkerlist[f]);CHKERRQ(ierr);
986           ierr = DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
987         }
988       }
989     }
990     ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr);
991   }
992   ierr = PLCDestroy(&in);CHKERRQ(ierr);
993   ierr = PLCDestroy(&out);CHKERRQ(ierr);
994   PetscFunctionReturn(0);
995 }
996 #endif /* PETSC_HAVE_CTETGEN */
997 
998 #undef __FUNCT__
999 #define __FUNCT__ "DMPlexGenerate"
1000 /*@C
1001   DMPlexGenerate - Generates a mesh.
1002 
1003   Not Collective
1004 
1005   Input Parameters:
1006 + boundary - The DMPlex boundary object
1007 . name - The mesh generation package name
1008 - interpolate - Flag to create intermediate mesh elements
1009 
1010   Output Parameter:
1011 . mesh - The DMPlex object
1012 
1013   Level: intermediate
1014 
1015 .keywords: mesh, elements
1016 .seealso: DMPlexCreate(), DMRefine()
1017 @*/
1018 PetscErrorCode DMPlexGenerate(DM boundary, const char name[], PetscBool interpolate, DM *mesh)
1019 {
1020   PetscInt       dim;
1021   char           genname[1024];
1022   PetscBool      isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;
1023   PetscErrorCode ierr;
1024 
1025   PetscFunctionBegin;
1026   PetscValidHeaderSpecific(boundary, DM_CLASSID, 1);
1027   PetscValidLogicalCollectiveBool(boundary, interpolate, 2);
1028   ierr = DMGetDimension(boundary, &dim);CHKERRQ(ierr);
1029   ierr = PetscOptionsGetString(((PetscObject) boundary)->prefix, "-dm_plex_generator", genname, 1024, &flg);CHKERRQ(ierr);
1030   if (flg) name = genname;
1031   if (name) {
1032     ierr = PetscStrcmp(name, "triangle", &isTriangle);CHKERRQ(ierr);
1033     ierr = PetscStrcmp(name, "tetgen",   &isTetgen);CHKERRQ(ierr);
1034     ierr = PetscStrcmp(name, "ctetgen",  &isCTetgen);CHKERRQ(ierr);
1035   }
1036   switch (dim) {
1037   case 1:
1038     if (!name || isTriangle) {
1039 #if defined(PETSC_HAVE_TRIANGLE)
1040       ierr = DMPlexGenerate_Triangle(boundary, interpolate, mesh);CHKERRQ(ierr);
1041 #else
1042       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation needs external package support.\nPlease reconfigure with --download-triangle.");
1043 #endif
1044     } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
1045     break;
1046   case 2:
1047     if (!name || isCTetgen) {
1048 #if defined(PETSC_HAVE_CTETGEN)
1049       ierr = DMPlexGenerate_CTetgen(boundary, interpolate, mesh);CHKERRQ(ierr);
1050 #else
1051       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
1052 #endif
1053     } else if (isTetgen) {
1054 #if defined(PETSC_HAVE_TETGEN)
1055       ierr = DMPlexGenerate_Tetgen(boundary, interpolate, mesh);CHKERRQ(ierr);
1056 #else
1057       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --with-c-language=cxx --download-tetgen.");
1058 #endif
1059     } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
1060     break;
1061   default:
1062     SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation for a dimension %d boundary is not supported.", dim);
1063   }
1064   PetscFunctionReturn(0);
1065 }
1066 
1067 #undef __FUNCT__
1068 #define __FUNCT__ "DMRefine_Plex"
1069 PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined)
1070 {
1071   PetscReal      refinementLimit;
1072   PetscInt       dim, cStart, cEnd;
1073   char           genname[1024], *name = NULL;
1074   PetscBool      isUniform, isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;
1075   PetscErrorCode ierr;
1076 
1077   PetscFunctionBegin;
1078   ierr = DMPlexGetRefinementUniform(dm, &isUniform);CHKERRQ(ierr);
1079   if (isUniform) {
1080     CellRefiner cellRefiner;
1081 
1082     ierr = DMPlexGetCellRefiner_Internal(dm, &cellRefiner);CHKERRQ(ierr);
1083     ierr = DMPlexRefineUniform_Internal(dm, cellRefiner, dmRefined);CHKERRQ(ierr);
1084     ierr = DMPlexCopyBoundary(dm, *dmRefined);CHKERRQ(ierr);
1085     PetscFunctionReturn(0);
1086   }
1087   ierr = DMPlexGetRefinementLimit(dm, &refinementLimit);CHKERRQ(ierr);
1088   if (refinementLimit == 0.0) PetscFunctionReturn(0);
1089   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1090   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1091   ierr = PetscOptionsGetString(((PetscObject) dm)->prefix, "-dm_plex_generator", genname, 1024, &flg);CHKERRQ(ierr);
1092   if (flg) name = genname;
1093   if (name) {
1094     ierr = PetscStrcmp(name, "triangle", &isTriangle);CHKERRQ(ierr);
1095     ierr = PetscStrcmp(name, "tetgen",   &isTetgen);CHKERRQ(ierr);
1096     ierr = PetscStrcmp(name, "ctetgen",  &isCTetgen);CHKERRQ(ierr);
1097   }
1098   switch (dim) {
1099   case 2:
1100     if (!name || isTriangle) {
1101 #if defined(PETSC_HAVE_TRIANGLE)
1102       double  *maxVolumes;
1103       PetscInt c;
1104 
1105       ierr = PetscMalloc1((cEnd - cStart), &maxVolumes);CHKERRQ(ierr);
1106       for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1107       ierr = DMPlexRefine_Triangle(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
1108       ierr = PetscFree(maxVolumes);CHKERRQ(ierr);
1109 #else
1110       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement needs external package support.\nPlease reconfigure with --download-triangle.");
1111 #endif
1112     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
1113     break;
1114   case 3:
1115     if (!name || isCTetgen) {
1116 #if defined(PETSC_HAVE_CTETGEN)
1117       PetscReal *maxVolumes;
1118       PetscInt   c;
1119 
1120       ierr = PetscMalloc1((cEnd - cStart), &maxVolumes);CHKERRQ(ierr);
1121       for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1122       ierr = DMPlexRefine_CTetgen(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
1123 #else
1124       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
1125 #endif
1126     } else if (isTetgen) {
1127 #if defined(PETSC_HAVE_TETGEN)
1128       double  *maxVolumes;
1129       PetscInt c;
1130 
1131       ierr = PetscMalloc1((cEnd - cStart), &maxVolumes);CHKERRQ(ierr);
1132       for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1133       ierr = DMPlexRefine_Tetgen(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
1134 #else
1135       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --with-c-language=cxx --download-tetgen.");
1136 #endif
1137     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
1138     break;
1139   default:
1140     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %d is not supported.", dim);
1141   }
1142   ierr = DMPlexCopyBoundary(dm, *dmRefined);CHKERRQ(ierr);
1143   PetscFunctionReturn(0);
1144 }
1145 
1146 #undef __FUNCT__
1147 #define __FUNCT__ "DMRefineHierarchy_Plex"
1148 PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[])
1149 {
1150   DM             cdm = dm;
1151   PetscInt       r;
1152   PetscBool      isUniform;
1153   PetscErrorCode ierr;
1154 
1155   PetscFunctionBegin;
1156   ierr = DMPlexGetRefinementUniform(dm, &isUniform);CHKERRQ(ierr);
1157   if (!isUniform) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Non-uniform refinement is incompatible with the hierarchy");
1158   for (r = 0; r < nlevels; ++r) {
1159     CellRefiner cellRefiner;
1160 
1161     ierr = DMPlexGetCellRefiner_Internal(cdm, &cellRefiner);CHKERRQ(ierr);
1162     ierr = DMPlexRefineUniform_Internal(cdm, cellRefiner, &dmRefined[r]);CHKERRQ(ierr);
1163     ierr = DMPlexSetCoarseDM(dmRefined[r], cdm);CHKERRQ(ierr);
1164     cdm  = dmRefined[r];
1165   }
1166   PetscFunctionReturn(0);
1167 }
1168 
1169 #undef __FUNCT__
1170 #define __FUNCT__ "DMCoarsen_Plex"
1171 PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened)
1172 {
1173   DM_Plex       *mesh = (DM_Plex*) dm->data;
1174   PetscErrorCode ierr;
1175 
1176   PetscFunctionBegin;
1177   ierr = PetscObjectReference((PetscObject) mesh->coarseMesh);CHKERRQ(ierr);
1178   *dmCoarsened = mesh->coarseMesh;
1179   PetscFunctionReturn(0);
1180 }
1181