xref: /petsc/src/dm/impls/plex/plexgenerate.c (revision fe998a80077c9ee0917a39496df43fc256e1b478) 
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 = DMPlexGetLabel(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 = DMPlexCreateLabel(*dm, labelName); ierr = DMPlexGetLabel(*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 = MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
355   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
356   ierr = DMPlexGetLabel(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 = DMPlexCreateLabel(*dmRefined, labelName); ierr = DMPlexGetLabel(*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 = DMPlexGetLabel(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 = DMPlexCreateLabel(*dm, labelName); ierr = DMPlexGetLabel(*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 = MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
627   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
628   ierr = DMPlexGetLabel(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 = DMPlexCreateLabel(*dmRefined, labelName); ierr = DMPlexGetLabel(*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(((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 = DMPlexGetLabel(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     const double   *meshCoords = out->pointlist;
842     int            *cells      = out->tetrahedronlist;
843 
844     ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr);
845     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr);
846     if (label) {ierr = DMPlexCreateLabel(*dm, labelName); ierr = DMPlexGetLabel(*dm, labelName, &glabel);}
847     /* Set labels */
848     for (v = 0; v < numVertices; ++v) {
849       if (out->pointmarkerlist[v]) {
850         if (glabel) {ierr = DMLabelSetValue(glabel, v+numCells, out->pointmarkerlist[v]);CHKERRQ(ierr);}
851       }
852     }
853     if (interpolate) {
854       PetscInt e;
855 
856       for (e = 0; e < out->numberofedges; e++) {
857         if (out->edgemarkerlist[e]) {
858           const PetscInt  vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
859           const PetscInt *edges;
860           PetscInt        numEdges;
861 
862           ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
863           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
864           if (glabel) {ierr = DMLabelSetValue(glabel, edges[0], out->edgemarkerlist[e]);CHKERRQ(ierr);}
865           ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
866         }
867       }
868       for (f = 0; f < out->numberoftrifaces; f++) {
869         if (out->trifacemarkerlist[f]) {
870           const PetscInt  vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
871           const PetscInt *faces;
872           PetscInt        numFaces;
873 
874           ierr = DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
875           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
876           if (glabel) {ierr = DMLabelSetValue(glabel, faces[0], out->trifacemarkerlist[f]);CHKERRQ(ierr);}
877           ierr = DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
878         }
879       }
880     }
881     ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
882   }
883 
884   ierr = PLCDestroy(&in);CHKERRQ(ierr);
885   ierr = PLCDestroy(&out);CHKERRQ(ierr);
886   PetscFunctionReturn(0);
887 }
888 
889 #undef __FUNCT__
890 #define __FUNCT__ "DMPlexRefine_CTetgen"
891 PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
892 {
893   MPI_Comm       comm;
894   const PetscInt dim       = 3;
895   const char    *labelName = "marker";
896   PLC           *in, *out;
897   DMLabel        label;
898   PetscInt       verbose = 0, vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
899   PetscMPIInt    rank;
900   PetscErrorCode ierr;
901 
902   PetscFunctionBegin;
903   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
904   ierr = PetscOptionsGetInt(((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);CHKERRQ(ierr);
905   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
906   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
907   ierr = MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
908   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
909   ierr = DMPlexGetLabel(dm, labelName, &label);CHKERRQ(ierr);
910   ierr = PLCCreate(&in);CHKERRQ(ierr);
911   ierr = PLCCreate(&out);CHKERRQ(ierr);
912 
913   in->numberofpoints = vEnd - vStart;
914   if (in->numberofpoints > 0) {
915     PetscSection coordSection;
916     Vec          coordinates;
917     PetscScalar *array;
918 
919     ierr = PetscMalloc1(in->numberofpoints*dim, &in->pointlist);CHKERRQ(ierr);
920     ierr = PetscMalloc1(in->numberofpoints,       &in->pointmarkerlist);CHKERRQ(ierr);
921     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
922     ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
923     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
924     for (v = vStart; v < vEnd; ++v) {
925       const PetscInt idx = v - vStart;
926       PetscInt       off, d, m = -1;
927 
928       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
929       for (d = 0; d < dim; ++d) {
930         in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
931       }
932       if (label) {ierr = DMLabelGetValue(label, v, &m);CHKERRQ(ierr);}
933 
934       in->pointmarkerlist[idx] = (int) m;
935     }
936     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
937   }
938   ierr  = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
939 
940   in->numberofcorners       = 4;
941   in->numberoftetrahedra    = cEnd - cStart;
942   in->tetrahedronvolumelist = maxVolumes;
943   if (in->numberoftetrahedra > 0) {
944     ierr = PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);CHKERRQ(ierr);
945     for (c = cStart; c < cEnd; ++c) {
946       const PetscInt idx      = c - cStart;
947       PetscInt      *closure = NULL;
948       PetscInt       closureSize;
949 
950       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
951       if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
952       for (v = 0; v < 4; ++v) {
953         in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
954       }
955       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
956     }
957   }
958   if (!rank) {
959     TetGenOpts t;
960 
961     ierr = TetGenOptsInitialize(&t);CHKERRQ(ierr);
962 
963     t.in        = dm; /* Should go away */
964     t.refine    = 1;
965     t.varvolume = 1;
966     t.quality   = 1;
967     t.edgesout  = 1;
968     t.zeroindex = 1;
969     t.quiet     = 1;
970     t.verbose   = verbose; /* Change this */
971 
972     ierr = TetGenCheckOpts(&t);CHKERRQ(ierr);
973     ierr = TetGenTetrahedralize(&t, in, out);CHKERRQ(ierr);
974   }
975   {
976     DMLabel        rlabel      = NULL;
977     const PetscInt numCorners  = 4;
978     const PetscInt numCells    = out->numberoftetrahedra;
979     const PetscInt numVertices = out->numberofpoints;
980     const double   *meshCoords = out->pointlist;
981     int            *cells      = out->tetrahedronlist;
982     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
983 
984     ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr);
985     ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr);
986     if (label) {ierr = DMPlexCreateLabel(*dmRefined, labelName); ierr = DMPlexGetLabel(*dmRefined, labelName, &rlabel);}
987     /* Set labels */
988     for (v = 0; v < numVertices; ++v) {
989       if (out->pointmarkerlist[v]) {
990         if (rlabel) {ierr = DMLabelSetValue(rlabel, v+numCells, out->pointmarkerlist[v]);CHKERRQ(ierr);}
991       }
992     }
993     if (interpolate) {
994       PetscInt e, f;
995 
996       for (e = 0; e < out->numberofedges; e++) {
997         if (out->edgemarkerlist[e]) {
998           const PetscInt  vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
999           const PetscInt *edges;
1000           PetscInt        numEdges;
1001 
1002           ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
1003           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
1004           if (rlabel) {ierr = DMLabelSetValue(rlabel, edges[0], out->edgemarkerlist[e]);CHKERRQ(ierr);}
1005           ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
1006         }
1007       }
1008       for (f = 0; f < out->numberoftrifaces; f++) {
1009         if (out->trifacemarkerlist[f]) {
1010           const PetscInt  vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
1011           const PetscInt *faces;
1012           PetscInt        numFaces;
1013 
1014           ierr = DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
1015           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
1016           if (rlabel) {ierr = DMLabelSetValue(rlabel, faces[0], out->trifacemarkerlist[f]);CHKERRQ(ierr);}
1017           ierr = DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
1018         }
1019       }
1020     }
1021     ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr);
1022   }
1023   ierr = PLCDestroy(&in);CHKERRQ(ierr);
1024   ierr = PLCDestroy(&out);CHKERRQ(ierr);
1025   PetscFunctionReturn(0);
1026 }
1027 #endif /* PETSC_HAVE_CTETGEN */
1028 
1029 #undef __FUNCT__
1030 #define __FUNCT__ "DMPlexGenerate"
1031 /*@C
1032   DMPlexGenerate - Generates a mesh.
1033 
1034   Not Collective
1035 
1036   Input Parameters:
1037 + boundary - The DMPlex boundary object
1038 . name - The mesh generation package name
1039 - interpolate - Flag to create intermediate mesh elements
1040 
1041   Output Parameter:
1042 . mesh - The DMPlex object
1043 
1044   Level: intermediate
1045 
1046 .keywords: mesh, elements
1047 .seealso: DMPlexCreate(), DMRefine()
1048 @*/
1049 PetscErrorCode DMPlexGenerate(DM boundary, const char name[], PetscBool interpolate, DM *mesh)
1050 {
1051   PetscInt       dim;
1052   char           genname[1024];
1053   PetscBool      isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;
1054   PetscErrorCode ierr;
1055 
1056   PetscFunctionBegin;
1057   PetscValidHeaderSpecific(boundary, DM_CLASSID, 1);
1058   PetscValidLogicalCollectiveBool(boundary, interpolate, 2);
1059   ierr = DMGetDimension(boundary, &dim);CHKERRQ(ierr);
1060   ierr = PetscOptionsGetString(((PetscObject) boundary)->prefix, "-dm_plex_generator", genname, 1024, &flg);CHKERRQ(ierr);
1061   if (flg) name = genname;
1062   if (name) {
1063     ierr = PetscStrcmp(name, "triangle", &isTriangle);CHKERRQ(ierr);
1064     ierr = PetscStrcmp(name, "tetgen",   &isTetgen);CHKERRQ(ierr);
1065     ierr = PetscStrcmp(name, "ctetgen",  &isCTetgen);CHKERRQ(ierr);
1066   }
1067   switch (dim) {
1068   case 1:
1069     if (!name || isTriangle) {
1070 #if defined(PETSC_HAVE_TRIANGLE)
1071       ierr = DMPlexGenerate_Triangle(boundary, interpolate, mesh);CHKERRQ(ierr);
1072 #else
1073       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation needs external package support.\nPlease reconfigure with --download-triangle.");
1074 #endif
1075     } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
1076     break;
1077   case 2:
1078     if (!name || isCTetgen) {
1079 #if defined(PETSC_HAVE_CTETGEN)
1080       ierr = DMPlexGenerate_CTetgen(boundary, interpolate, mesh);CHKERRQ(ierr);
1081 #else
1082       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
1083 #endif
1084     } else if (isTetgen) {
1085 #if defined(PETSC_HAVE_TETGEN)
1086       ierr = DMPlexGenerate_Tetgen(boundary, interpolate, mesh);CHKERRQ(ierr);
1087 #else
1088       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --with-c-language=cxx --download-tetgen.");
1089 #endif
1090     } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
1091     break;
1092   default:
1093     SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation for a dimension %d boundary is not supported.", dim);
1094   }
1095   PetscFunctionReturn(0);
1096 }
1097 
1098 #undef __FUNCT__
1099 #define __FUNCT__ "DMRefine_Plex"
1100 PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined)
1101 {
1102   PetscReal      refinementLimit;
1103   PetscInt       dim, cStart, cEnd;
1104   char           genname[1024], *name = NULL;
1105   PetscBool      isUniform, isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;
1106   PetscErrorCode ierr;
1107 
1108   PetscFunctionBegin;
1109   ierr = DMPlexGetRefinementUniform(dm, &isUniform);CHKERRQ(ierr);
1110   if (isUniform) {
1111     CellRefiner cellRefiner;
1112 
1113     ierr = DMPlexGetCellRefiner_Internal(dm, &cellRefiner);CHKERRQ(ierr);
1114     ierr = DMPlexRefineUniform_Internal(dm, cellRefiner, dmRefined);CHKERRQ(ierr);
1115     ierr = DMPlexCopyBoundary(dm, *dmRefined);CHKERRQ(ierr);
1116     PetscFunctionReturn(0);
1117   }
1118   ierr = DMPlexGetRefinementLimit(dm, &refinementLimit);CHKERRQ(ierr);
1119   if (refinementLimit == 0.0) PetscFunctionReturn(0);
1120   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1121   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1122   ierr = PetscOptionsGetString(((PetscObject) dm)->prefix, "-dm_plex_generator", genname, 1024, &flg);CHKERRQ(ierr);
1123   if (flg) name = genname;
1124   if (name) {
1125     ierr = PetscStrcmp(name, "triangle", &isTriangle);CHKERRQ(ierr);
1126     ierr = PetscStrcmp(name, "tetgen",   &isTetgen);CHKERRQ(ierr);
1127     ierr = PetscStrcmp(name, "ctetgen",  &isCTetgen);CHKERRQ(ierr);
1128   }
1129   switch (dim) {
1130   case 2:
1131     if (!name || isTriangle) {
1132 #if defined(PETSC_HAVE_TRIANGLE)
1133       double  *maxVolumes;
1134       PetscInt c;
1135 
1136       ierr = PetscMalloc1(cEnd - cStart, &maxVolumes);CHKERRQ(ierr);
1137       for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1138       ierr = DMPlexRefine_Triangle(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
1139       ierr = PetscFree(maxVolumes);CHKERRQ(ierr);
1140 #else
1141       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement needs external package support.\nPlease reconfigure with --download-triangle.");
1142 #endif
1143     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
1144     break;
1145   case 3:
1146     if (!name || isCTetgen) {
1147 #if defined(PETSC_HAVE_CTETGEN)
1148       PetscReal *maxVolumes;
1149       PetscInt   c;
1150 
1151       ierr = PetscMalloc1(cEnd - cStart, &maxVolumes);CHKERRQ(ierr);
1152       for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1153       ierr = DMPlexRefine_CTetgen(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
1154 #else
1155       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
1156 #endif
1157     } else if (isTetgen) {
1158 #if defined(PETSC_HAVE_TETGEN)
1159       double  *maxVolumes;
1160       PetscInt c;
1161 
1162       ierr = PetscMalloc1(cEnd - cStart, &maxVolumes);CHKERRQ(ierr);
1163       for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1164       ierr = DMPlexRefine_Tetgen(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
1165 #else
1166       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --with-c-language=cxx --download-tetgen.");
1167 #endif
1168     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
1169     break;
1170   default:
1171     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %d is not supported.", dim);
1172   }
1173   ierr = DMPlexCopyBoundary(dm, *dmRefined);CHKERRQ(ierr);
1174   PetscFunctionReturn(0);
1175 }
1176 
1177 #undef __FUNCT__
1178 #define __FUNCT__ "DMRefineHierarchy_Plex"
1179 PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[])
1180 {
1181   DM             cdm = dm;
1182   PetscInt       r;
1183   PetscBool      isUniform;
1184   PetscErrorCode ierr;
1185 
1186   PetscFunctionBegin;
1187   ierr = DMPlexGetRefinementUniform(dm, &isUniform);CHKERRQ(ierr);
1188   if (!isUniform) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Non-uniform refinement is incompatible with the hierarchy");
1189   for (r = 0; r < nlevels; ++r) {
1190     CellRefiner cellRefiner;
1191 
1192     ierr = DMPlexGetCellRefiner_Internal(cdm, &cellRefiner);CHKERRQ(ierr);
1193     ierr = DMPlexRefineUniform_Internal(cdm, cellRefiner, &dmRefined[r]);CHKERRQ(ierr);
1194     ierr = DMPlexSetCoarseDM(dmRefined[r], cdm);CHKERRQ(ierr);
1195     cdm  = dmRefined[r];
1196   }
1197   PetscFunctionReturn(0);
1198 }
1199 
1200 #undef __FUNCT__
1201 #define __FUNCT__ "DMCoarsen_Plex"
1202 PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened)
1203 {
1204   DM_Plex       *mesh = (DM_Plex*) dm->data;
1205   PetscErrorCode ierr;
1206 
1207   PetscFunctionBegin;
1208   ierr = PetscObjectReference((PetscObject) mesh->coarseMesh);CHKERRQ(ierr);
1209   *dmCoarsened = mesh->coarseMesh;
1210   PetscFunctionReturn(0);
1211 }
1212