xref: /petsc/src/dm/impls/plex/plexgenerate.c (revision 966be33a19c9230d4aa438248a644248d45cc287) 
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 #undef __FUNCT__
1129 #define __FUNCT__ "DMRefine_Plex"
1130 PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined)
1131 {
1132   PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *);
1133   PetscReal        refinementLimit;
1134   PetscInt         dim, cStart, cEnd;
1135   char             genname[1024], *name = NULL;
1136   PetscBool        isUniform, isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg, localized;
1137   PetscErrorCode   ierr;
1138 
1139   PetscFunctionBegin;
1140   ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr);
1141   ierr = DMPlexGetRefinementUniform(dm, &isUniform);CHKERRQ(ierr);
1142   if (isUniform) {
1143     CellRefiner cellRefiner;
1144 
1145     ierr = DMPlexGetCellRefiner_Internal(dm, &cellRefiner);CHKERRQ(ierr);
1146     ierr = DMPlexRefineUniform_Internal(dm, cellRefiner, dmRefined);CHKERRQ(ierr);
1147     ierr = DMCopyBoundary(dm, *dmRefined);CHKERRQ(ierr);
1148     if (localized) {
1149       ierr = DMLocalizeCoordinates(*dmRefined);CHKERRQ(ierr);
1150     }
1151     PetscFunctionReturn(0);
1152   }
1153   ierr = DMPlexGetRefinementLimit(dm, &refinementLimit);CHKERRQ(ierr);
1154   ierr = DMPlexGetRefinementFunction(dm, &refinementFunc);CHKERRQ(ierr);
1155   if (refinementLimit == 0.0 && !refinementFunc) PetscFunctionReturn(0);
1156   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1157   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1158   ierr = PetscOptionsGetString(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_generator", genname, 1024, &flg);CHKERRQ(ierr);
1159   if (flg) name = genname;
1160   if (name) {
1161     ierr = PetscStrcmp(name, "triangle", &isTriangle);CHKERRQ(ierr);
1162     ierr = PetscStrcmp(name, "tetgen",   &isTetgen);CHKERRQ(ierr);
1163     ierr = PetscStrcmp(name, "ctetgen",  &isCTetgen);CHKERRQ(ierr);
1164   }
1165   switch (dim) {
1166   case 2:
1167     if (!name || isTriangle) {
1168 #if defined(PETSC_HAVE_TRIANGLE)
1169       double  *maxVolumes;
1170       PetscInt c;
1171 
1172       ierr = PetscMalloc1(cEnd - cStart, &maxVolumes);CHKERRQ(ierr);
1173       if (refinementFunc) {
1174         for (c = cStart; c < cEnd; ++c) {
1175           PetscReal vol, centroid[3];
1176           PetscReal maxVol;
1177 
1178           ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);CHKERRQ(ierr);
1179           ierr = (*refinementFunc)(centroid, &maxVol);CHKERRQ(ierr);
1180           maxVolumes[c - cStart] = (double) maxVol;
1181         }
1182       } else {
1183         for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1184       }
1185       ierr = DMPlexRefine_Triangle(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
1186       ierr = PetscFree(maxVolumes);CHKERRQ(ierr);
1187 #else
1188       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement needs external package support.\nPlease reconfigure with --download-triangle.");
1189 #endif
1190     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
1191     break;
1192   case 3:
1193     if (!name || isCTetgen) {
1194 #if defined(PETSC_HAVE_CTETGEN)
1195       PetscReal *maxVolumes;
1196       PetscInt   c;
1197 
1198       ierr = PetscMalloc1(cEnd - cStart, &maxVolumes);CHKERRQ(ierr);
1199       if (refinementFunc) {
1200         for (c = cStart; c < cEnd; ++c) {
1201           PetscReal vol, centroid[3];
1202 
1203           ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);CHKERRQ(ierr);
1204           ierr = (*refinementFunc)(centroid, &maxVolumes[c-cStart]);CHKERRQ(ierr);
1205         }
1206       } else {
1207         for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1208       }
1209       ierr = DMPlexRefine_CTetgen(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
1210 #else
1211       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
1212 #endif
1213     } else if (isTetgen) {
1214 #if defined(PETSC_HAVE_TETGEN)
1215       double  *maxVolumes;
1216       PetscInt c;
1217 
1218       ierr = PetscMalloc1(cEnd - cStart, &maxVolumes);CHKERRQ(ierr);
1219       if (refinementFunc) {
1220         for (c = cStart; c < cEnd; ++c) {
1221           PetscReal vol, centroid[3];
1222 
1223           ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);CHKERRQ(ierr);
1224           ierr = (*refinementFunc)(centroid, &maxVolumes[c-cStart]);CHKERRQ(ierr);
1225         }
1226       } else {
1227         for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1228       }
1229       ierr = DMPlexRefine_Tetgen(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
1230 #else
1231       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --with-c-language=cxx --download-tetgen.");
1232 #endif
1233     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
1234     break;
1235   default:
1236     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %d is not supported.", dim);
1237   }
1238   ierr = DMCopyBoundary(dm, *dmRefined);CHKERRQ(ierr);
1239   if (localized) {
1240     ierr = DMLocalizeCoordinates(*dmRefined);CHKERRQ(ierr);
1241   }
1242   PetscFunctionReturn(0);
1243 }
1244 
1245 #undef __FUNCT__
1246 #define __FUNCT__ "DMRefineHierarchy_Plex"
1247 PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[])
1248 {
1249   DM             cdm = dm;
1250   PetscInt       r;
1251   PetscBool      isUniform, localized;
1252   PetscErrorCode ierr;
1253 
1254   PetscFunctionBegin;
1255   ierr = DMPlexGetRefinementUniform(dm, &isUniform);CHKERRQ(ierr);
1256   ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr);
1257   if (!isUniform) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Non-uniform refinement is incompatible with the hierarchy");
1258   for (r = 0; r < nlevels; ++r) {
1259     CellRefiner cellRefiner;
1260 
1261     ierr = DMPlexGetCellRefiner_Internal(cdm, &cellRefiner);CHKERRQ(ierr);
1262     ierr = DMPlexRefineUniform_Internal(cdm, cellRefiner, &dmRefined[r]);CHKERRQ(ierr);
1263     ierr = DMCopyBoundary(cdm, dmRefined[r]);CHKERRQ(ierr);
1264     if (localized) {
1265       ierr = DMLocalizeCoordinates(dmRefined[r]);CHKERRQ(ierr);
1266     }
1267     ierr = DMSetCoarseDM(dmRefined[r], cdm);CHKERRQ(ierr);
1268     ierr = DMPlexSetRegularRefinement(dmRefined[r], PETSC_TRUE);CHKERRQ(ierr);
1269     cdm  = dmRefined[r];
1270   }
1271   PetscFunctionReturn(0);
1272 }
1273