Actual source code: plexgenerate.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/

  5: PetscErrorCode DMPlexInvertCell_Internal(PetscInt dim, PetscInt numCorners, PetscInt cone[])
  6: {
  7:   int tmpc;

 10:   if (dim != 3) return(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:   return(0);
 25: }

 29: /*@C
 30:   DMPlexInvertCell - This flips tetrahedron and hexahedron orientation since Plex stores them internally with outward normals. Other cells are left untouched.

 32:   Input Parameters:
 33: + numCorners - The number of vertices in a cell
 34: - cone - The incoming cone

 36:   Output Parameter:
 37: . cone - The inverted cone (in-place)

 39:   Level: developer

 41: .seealso: DMPlexGenerate()
 42: @*/
 43: PetscErrorCode DMPlexInvertCell(PetscInt dim, PetscInt numCorners, int cone[])
 44: {
 45:   int tmpc;

 48:   if (dim != 3) return(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:   return(0);
 63: }

 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;

 74:   for (coff = 0; coff < bound; coff += numCorners) {
 75:     DMPlexInvertCell(dim, numCorners, &cells[coff]);
 76:   }
 77:   return(0);
 78: }

 82: /*@
 83:   DMPlexTriangleSetOptions - Set the options used for the Triangle mesh generator

 85:   Not Collective

 87:   Inputs Parameters:
 88: + dm - The DMPlex object
 89: - opts - The command line options

 91:   Level: developer

 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;

104:   PetscFree(mesh->triangleOpts);
105:   PetscStrallocpy(opts, &mesh->triangleOpts);
106:   return(0);
107: }

111: /*@
112:   DMPlexTetgenSetOptions - Set the options used for the Tetgen mesh generator

114:   Not Collective

116:   Inputs Parameters:
117: + dm - The DMPlex object
118: - opts - The command line options

120:   Level: developer

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;

133:   PetscFree(mesh->tetgenOpts);
134:   PetscStrallocpy(opts, &mesh->tetgenOpts);
135:   return(0);
136: }

138: #if defined(PETSC_HAVE_TRIANGLE)
139: #include <triangle.h>

143: PetscErrorCode InitInput_Triangle(struct triangulateio *inputCtx)
144: {
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:   return(0);
161: }

165: PetscErrorCode InitOutput_Triangle(struct triangulateio *outputCtx)
166: {
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:   return(0);
182: }

186: PetscErrorCode FiniOutput_Triangle(struct triangulateio *outputCtx)
187: {
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:   return(0);
198: }

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;

218:   PetscObjectGetComm((PetscObject)boundary,&comm);
219:   MPI_Comm_rank(comm, &rank);
220:   InitInput_Triangle(&in);
221:   InitOutput_Triangle(&out);
222:   DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
223:   DMGetLabel(boundary, labelName, &label);

225:   in.numberofpoints = vEnd - vStart;
226:   if (in.numberofpoints > 0) {
227:     PetscSection coordSection;
228:     Vec          coordinates;
229:     PetscScalar *array;

231:     PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
232:     PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
233:     DMGetCoordinatesLocal(boundary, &coordinates);
234:     DMGetCoordinateSection(boundary, &coordSection);
235:     VecGetArray(coordinates, &array);
236:     for (v = vStart; v < vEnd; ++v) {
237:       const PetscInt idx = v - vStart;
238:       PetscInt       off, d;

240:       PetscSectionGetOffset(coordSection, v, &off);
241:       for (d = 0; d < dim; ++d) {
242:         in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
243:       }
244:       if (label) {DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);}
245:     }
246:     VecRestoreArray(coordinates, &array);
247:   }
248:   DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);
249:   in.numberofsegments = eEnd - eStart;
250:   if (in.numberofsegments > 0) {
251:     PetscMalloc1(in.numberofsegments*2, &in.segmentlist);
252:     PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist);
253:     for (e = eStart; e < eEnd; ++e) {
254:       const PetscInt  idx = e - eStart;
255:       const PetscInt *cone;

257:       DMPlexGetCone(boundary, e, &cone);

259:       in.segmentlist[idx*2+0] = cone[0] - vStart;
260:       in.segmentlist[idx*2+1] = cone[1] - vStart;

262:       if (label) {DMLabelGetValue(label, e, &in.segmentmarkerlist[idx]);}
263:     }
264:   }
265: #if 0 /* Do not currently support holes */
266:   PetscReal *holeCoords;
267:   PetscInt   h, d;

269:   DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
270:   if (in.numberofholes > 0) {
271:     PetscMalloc1(in.numberofholes*dim, &in.holelist);
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];

282:     /* Take away 'Q' for verbose output */
283:     PetscStrcpy(args, "pqezQ");
284:     if (createConvexHull)   {PetscStrcat(args, "c");}
285:     if (constrained)        {PetscStrcpy(args, "zepDQ");}
286:     if (mesh->triangleOpts) {triangulate(mesh->triangleOpts, &in, &out, NULL);}
287:     else                    {triangulate(args, &in, &out, NULL);}
288:   }
289:   PetscFree(in.pointlist);
290:   PetscFree(in.pointmarkerlist);
291:   PetscFree(in.segmentlist);
292:   PetscFree(in.segmentmarkerlist);
293:   PetscFree(in.holelist);

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;

303:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
304:     if (label) {DMCreateLabel(*dm, labelName); DMGetLabel(*dm, labelName, &glabel);}
305:     /* Set labels */
306:     for (v = 0; v < numVertices; ++v) {
307:       if (out.pointmarkerlist[v]) {
308:         if (glabel) {DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);}
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;

318:           DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
319:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
320:           if (glabel) {DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);}
321:           DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
322:         }
323:       }
324:     }
325:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
326:   }
327: #if 0 /* Do not currently support holes */
328:   DMPlexCopyHoles(*dm, boundary);
329: #endif
330:   FiniOutput_Triangle(&out);
331:   return(0);
332: }

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;

349:   PetscObjectGetComm((PetscObject)dm,&comm);
350:   MPI_Comm_rank(comm, &rank);
351:   InitInput_Triangle(&in);
352:   InitOutput_Triangle(&out);
353:   DMPlexGetDepth(dm, &depth);
354:   MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
355:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
356:   DMGetLabel(dm, labelName, &label);

358:   in.numberofpoints = vEnd - vStart;
359:   if (in.numberofpoints > 0) {
360:     PetscSection coordSection;
361:     Vec          coordinates;
362:     PetscScalar *array;

364:     PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
365:     PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
366:     DMGetCoordinatesLocal(dm, &coordinates);
367:     DMGetCoordinateSection(dm, &coordSection);
368:     VecGetArray(coordinates, &array);
369:     for (v = vStart; v < vEnd; ++v) {
370:       const PetscInt idx = v - vStart;
371:       PetscInt       off, d;

373:       PetscSectionGetOffset(coordSection, v, &off);
374:       for (d = 0; d < dim; ++d) {
375:         in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
376:       }
377:       if (label) {DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);}
378:     }
379:     VecRestoreArray(coordinates, &array);
380:   }
381:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);

383:   in.numberofcorners   = 3;
384:   in.numberoftriangles = cEnd - cStart;

386:   in.trianglearealist  = (double*) maxVolumes;
387:   if (in.numberoftriangles > 0) {
388:     PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);
389:     for (c = cStart; c < cEnd; ++c) {
390:       const PetscInt idx      = c - cStart;
391:       PetscInt      *closure = NULL;
392:       PetscInt       closureSize;

394:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
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:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
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;

407:   DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
408:   if (in.numberofholes > 0) {
409:     PetscMalloc1(in.numberofholes*dim, &in.holelist);
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];

420:     /* Take away 'Q' for verbose output */
421:     PetscStrcpy(args, "pqezQra");
422:     triangulate(args, &in, &out, NULL);
423:   }
424:   PetscFree(in.pointlist);
425:   PetscFree(in.pointmarkerlist);
426:   PetscFree(in.segmentlist);
427:   PetscFree(in.segmentmarkerlist);
428:   PetscFree(in.trianglelist);

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;

439:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
440:     if (label) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*dmRefined, labelName, &rlabel);}
441:     /* Set labels */
442:     for (v = 0; v < numVertices; ++v) {
443:       if (out.pointmarkerlist[v]) {
444:         if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);}
445:       }
446:     }
447:     if (interpolate) {
448:       PetscInt e;

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;

456:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
457:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
458:           if (rlabel) {DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);}
459:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
460:         }
461:       }
462:     }
463:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
464:   }
465: #if 0 /* Do not currently support holes */
466:   DMPlexCopyHoles(*dm, boundary);
467: #endif
468:   FiniOutput_Triangle(&out);
469:   return(0);
470: }
471: #endif /* PETSC_HAVE_TRIANGLE */

473: #if defined(PETSC_HAVE_TETGEN)
474: #include <tetgen.h>
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;

491:   PetscObjectGetComm((PetscObject)boundary,&comm);
492:   MPI_Comm_rank(comm, &rank);
493:   DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
494:   DMGetLabel(boundary, labelName, &label);

496:   in.numberofpoints = vEnd - vStart;
497:   if (in.numberofpoints > 0) {
498:     PetscSection coordSection;
499:     Vec          coordinates;
500:     PetscScalar *array;

502:     in.pointlist       = new double[in.numberofpoints*dim];
503:     in.pointmarkerlist = new int[in.numberofpoints];

505:     DMGetCoordinatesLocal(boundary, &coordinates);
506:     DMGetCoordinateSection(boundary, &coordSection);
507:     VecGetArray(coordinates, &array);
508:     for (v = vStart; v < vEnd; ++v) {
509:       const PetscInt idx = v - vStart;
510:       PetscInt       off, d;

512:       PetscSectionGetOffset(coordSection, v, &off);
513:       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
514:       if (label) {DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);}
515:     }
516:     VecRestoreArray(coordinates, &array);
517:   }
518:   DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);

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;

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;

533:       DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
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:       }

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) {DMLabelGetValue(label, f, &in.facetmarkerlist[idx]);}
547:       DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
548:     }
549:   }
550:   if (!rank) {
551:     char args[32];

553:     /* Take away 'Q' for verbose output */
554:     PetscStrcpy(args, "pqezQ");
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;

566:     DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
567:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
568:     if (label) {DMCreateLabel(*dm, labelName); DMGetLabel(*dm, labelName, &glabel);}
569:     /* Set labels */
570:     for (v = 0; v < numVertices; ++v) {
571:       if (out.pointmarkerlist[v]) {
572:         if (glabel) {DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);}
573:       }
574:     }
575:     if (interpolate) {
576:       PetscInt e;

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;

584:           DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
585:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
586:           if (glabel) {DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);}
587:           DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
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;

596:           DMPlexGetJoin(*dm, 3, vertices, &numFaces, &faces);
597:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
598:           if (glabel) {DMLabelSetValue(glabel, faces[0], out.trifacemarkerlist[f]);}
599:           DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
600:         }
601:       }
602:     }
603:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
604:   }
605:   return(0);
606: }

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;

623:   PetscObjectGetComm((PetscObject)dm,&comm);
624:   MPI_Comm_rank(comm, &rank);
625:   DMPlexGetDepth(dm, &depth);
626:   MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
627:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
628:   DMGetLabel(dm, labelName, &label);

630:   in.numberofpoints = vEnd - vStart;
631:   if (in.numberofpoints > 0) {
632:     PetscSection coordSection;
633:     Vec          coordinates;
634:     PetscScalar *array;

636:     in.pointlist       = new double[in.numberofpoints*dim];
637:     in.pointmarkerlist = new int[in.numberofpoints];

639:     DMGetCoordinatesLocal(dm, &coordinates);
640:     DMGetCoordinateSection(dm, &coordSection);
641:     VecGetArray(coordinates, &array);
642:     for (v = vStart; v < vEnd; ++v) {
643:       const PetscInt idx = v - vStart;
644:       PetscInt       off, d;

646:       PetscSectionGetOffset(coordSection, v, &off);
647:       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
648:       if (label) {DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);}
649:     }
650:     VecRestoreArray(coordinates, &array);
651:   }
652:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);

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;

664:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
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:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
670:     }
671:   }
672:   /* TODO: Put in boundary faces with markers */
673:   if (!rank) {
674:     char args[32];

676:     /* Take away 'Q' for verbose output */
677:     /*PetscStrcpy(args, "qezQra"); */
678:     PetscStrcpy(args, "qezraVVVV");
679:     ::tetrahedralize(args, &in, &out);
680:   }
681:   in.tetrahedronvolumelist = NULL;

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;

691:     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;

693:     DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
694:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
695:     if (label) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*dmRefined, labelName, &rlabel);}
696:     /* Set labels */
697:     for (v = 0; v < numVertices; ++v) {
698:       if (out.pointmarkerlist[v]) {
699:         if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);}
700:       }
701:     }
702:     if (interpolate) {
703:       PetscInt e, f;

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;

711:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
712:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
713:           if (rlabel) {DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);}
714:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
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;

723:           DMPlexGetJoin(*dmRefined, 3, vertices, &numFaces, &faces);
724:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
725:           if (rlabel) {DMLabelSetValue(rlabel, faces[0], out.trifacemarkerlist[f]);}
726:           DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
727:         }
728:       }
729:     }
730:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
731:   }
732:   return(0);
733: }
734: #endif /* PETSC_HAVE_TETGEN */

736: #if defined(PETSC_HAVE_CTETGEN)
737: #include <ctetgen.h>

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;

753:   PetscObjectGetComm((PetscObject)boundary,&comm);
754:   PetscOptionsGetInt(NULL,((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);
755:   MPI_Comm_rank(comm, &rank);
756:   DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
757:   DMGetLabel(boundary, labelName, &label);
758:   PLCCreate(&in);
759:   PLCCreate(&out);

761:   in->numberofpoints = vEnd - vStart;
762:   if (in->numberofpoints > 0) {
763:     PetscSection coordSection;
764:     Vec          coordinates;
765:     PetscScalar *array;

767:     PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
768:     PetscMalloc1(in->numberofpoints,       &in->pointmarkerlist);
769:     DMGetCoordinatesLocal(boundary, &coordinates);
770:     DMGetCoordinateSection(boundary, &coordSection);
771:     VecGetArray(coordinates, &array);
772:     for (v = vStart; v < vEnd; ++v) {
773:       const PetscInt idx = v - vStart;
774:       PetscInt       off, d, m = -1;

776:       PetscSectionGetOffset(coordSection, v, &off);
777:       for (d = 0; d < dim; ++d) {
778:         in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
779:       }
780:       if (label) {DMLabelGetValue(label, v, &m);}

782:       in->pointmarkerlist[idx] = (int) m;
783:     }
784:     VecRestoreArray(coordinates, &array);
785:   }
786:   DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);

788:   in->numberoffacets = fEnd - fStart;
789:   if (in->numberoffacets > 0) {
790:     PetscMalloc1(in->numberoffacets, &in->facetlist);
791:     PetscMalloc1(in->numberoffacets,   &in->facetmarkerlist);
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;

797:       in->facetlist[idx].numberofpolygons = 1;

799:       PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist);

801:       in->facetlist[idx].numberofholes    = 0;
802:       in->facetlist[idx].holelist         = NULL;

804:       DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
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:       }

810:       poly                   = in->facetlist[idx].polygonlist;
811:       poly->numberofvertices = numVertices;
812:       PetscMalloc1(poly->numberofvertices, &poly->vertexlist);
813:       for (v = 0; v < numVertices; ++v) {
814:         const PetscInt vIdx = points[v] - vStart;
815:         poly->vertexlist[v] = vIdx;
816:       }
817:       if (label) {DMLabelGetValue(label, f, &m);}
818:       in->facetmarkerlist[idx] = (int) m;
819:       DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
820:     }
821:   }
822:   if (!rank) {
823:     TetGenOpts t;

825:     TetGenOptsInitialize(&t);
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:     TetGenCheckOpts(&t);
834:     TetGenTetrahedralize(&t, in, out);
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;

844:     if (sizeof (PetscReal) == sizeof (double)) {
845:       meshCoords = (double *) out->pointlist;
846:     }
847:     else {
848:       PetscInt i;

850:       PetscMalloc1(3 * numVertices,&meshCoords);
851:       for (i = 0; i < 3 * numVertices; i++) {
852:         meshCoords[i] = (PetscReal) out->pointlist[i];
853:       }
854:     }

856:     DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
857:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
858:     if (sizeof (PetscReal) != sizeof (double)) {
859:       PetscFree(meshCoords);
860:     }
861:     if (label) {DMCreateLabel(*dm, labelName); DMGetLabel(*dm, labelName, &glabel);}
862:     /* Set labels */
863:     for (v = 0; v < numVertices; ++v) {
864:       if (out->pointmarkerlist[v]) {
865:         if (glabel) {DMLabelSetValue(glabel, v+numCells, out->pointmarkerlist[v]);}
866:       }
867:     }
868:     if (interpolate) {
869:       PetscInt e;

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;

877:           DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
878:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
879:           if (glabel) {DMLabelSetValue(glabel, edges[0], out->edgemarkerlist[e]);}
880:           DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
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;

889:           DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
890:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
891:           if (glabel) {DMLabelSetValue(glabel, faces[0], out->trifacemarkerlist[f]);}
892:           DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
893:         }
894:       }
895:     }
896:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
897:   }

899:   PLCDestroy(&in);
900:   PLCDestroy(&out);
901:   return(0);
902: }

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;

918:   PetscObjectGetComm((PetscObject)dm,&comm);
919:   PetscOptionsGetInt(NULL,((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);
920:   MPI_Comm_rank(comm, &rank);
921:   DMPlexGetDepth(dm, &depth);
922:   MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
923:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
924:   DMGetLabel(dm, labelName, &label);
925:   PLCCreate(&in);
926:   PLCCreate(&out);

928:   in->numberofpoints = vEnd - vStart;
929:   if (in->numberofpoints > 0) {
930:     PetscSection coordSection;
931:     Vec          coordinates;
932:     PetscScalar *array;

934:     PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
935:     PetscMalloc1(in->numberofpoints,       &in->pointmarkerlist);
936:     DMGetCoordinatesLocal(dm, &coordinates);
937:     DMGetCoordinateSection(dm, &coordSection);
938:     VecGetArray(coordinates, &array);
939:     for (v = vStart; v < vEnd; ++v) {
940:       const PetscInt idx = v - vStart;
941:       PetscInt       off, d, m = -1;

943:       PetscSectionGetOffset(coordSection, v, &off);
944:       for (d = 0; d < dim; ++d) {
945:         in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
946:       }
947:       if (label) {DMLabelGetValue(label, v, &m);}

949:       in->pointmarkerlist[idx] = (int) m;
950:     }
951:     VecRestoreArray(coordinates, &array);
952:   }
953:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);

955:   in->numberofcorners       = 4;
956:   in->numberoftetrahedra    = cEnd - cStart;
957:   in->tetrahedronvolumelist = maxVolumes;
958:   if (in->numberoftetrahedra > 0) {
959:     PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);
960:     for (c = cStart; c < cEnd; ++c) {
961:       const PetscInt idx      = c - cStart;
962:       PetscInt      *closure = NULL;
963:       PetscInt       closureSize;

965:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
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:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
971:     }
972:   }
973:   if (!rank) {
974:     TetGenOpts t;

976:     TetGenOptsInitialize(&t);

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 */

987:     TetGenCheckOpts(&t);
988:     TetGenTetrahedralize(&t, in, out);
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;

999:     if (sizeof (PetscReal) == sizeof (double)) {
1000:       meshCoords = (double *) out->pointlist;
1001:     }
1002:     else {
1003:       PetscInt i;

1005:       PetscMalloc1(3 * numVertices,&meshCoords);
1006:       for (i = 0; i < 3 * numVertices; i++) {
1007:         meshCoords[i] = (PetscReal) out->pointlist[i];
1008:       }
1009:     }

1011:     DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
1012:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
1013:     if (sizeof (PetscReal) != sizeof (double)) {
1014:       PetscFree(meshCoords);
1015:     }
1016:     if (label) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*dmRefined, labelName, &rlabel);}
1017:     /* Set labels */
1018:     for (v = 0; v < numVertices; ++v) {
1019:       if (out->pointmarkerlist[v]) {
1020:         if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out->pointmarkerlist[v]);}
1021:       }
1022:     }
1023:     if (interpolate) {
1024:       PetscInt e, f;

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;

1032:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
1033:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
1034:           if (rlabel) {DMLabelSetValue(rlabel, edges[0], out->edgemarkerlist[e]);}
1035:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
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;

1044:           DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
1045:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
1046:           if (rlabel) {DMLabelSetValue(rlabel, faces[0], out->trifacemarkerlist[f]);}
1047:           DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
1048:         }
1049:       }
1050:     }
1051:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
1052:   }
1053:   PLCDestroy(&in);
1054:   PLCDestroy(&out);
1055:   return(0);
1056: }
1057: #endif /* PETSC_HAVE_CTETGEN */

1061: /*@C
1062:   DMPlexGenerate - Generates a mesh.

1064:   Not Collective

1066:   Input Parameters:
1067: + boundary - The DMPlex boundary object
1068: . name - The mesh generation package name
1069: - interpolate - Flag to create intermediate mesh elements

1071:   Output Parameter:
1072: . mesh - The DMPlex object

1074:   Level: intermediate

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;

1089:   DMGetDimension(boundary, &dim);
1090:   PetscOptionsGetString(((PetscObject) boundary)->options,((PetscObject) boundary)->prefix, "-dm_plex_generator", genname, 1024, &flg);
1091:   if (flg) name = genname;
1092:   if (name) {
1093:     PetscStrcmp(name, "triangle", &isTriangle);
1094:     PetscStrcmp(name, "tetgen",   &isTetgen);
1095:     PetscStrcmp(name, "ctetgen",  &isCTetgen);
1096:   }
1097:   switch (dim) {
1098:   case 1:
1099:     if (!name || isTriangle) {
1100: #if defined(PETSC_HAVE_TRIANGLE)
1101:       DMPlexGenerate_Triangle(boundary, interpolate, mesh);
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:       DMPlexGenerate_CTetgen(boundary, interpolate, mesh);
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:       DMPlexGenerate_Tetgen(boundary, interpolate, mesh);
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:   return(0);
1126: }

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;

1140:   DMGetCoordinatesLocalized(dm, &localized);
1141:   DMPlexGetRefinementUniform(dm, &isUniform);
1142:   if (isUniform) {
1143:     CellRefiner cellRefiner;

1145:     DMPlexGetCellRefiner_Internal(dm, &cellRefiner);
1146:     DMPlexRefineUniform_Internal(dm, cellRefiner, dmRefined);
1147:     DMCopyBoundary(dm, *dmRefined);
1148:     if (localized) {
1149:       DMLocalizeCoordinates(*dmRefined);
1150:     }
1151:     return(0);
1152:   }
1153:   DMPlexGetRefinementLimit(dm, &refinementLimit);
1154:   DMPlexGetRefinementFunction(dm, &refinementFunc);
1155:   if (refinementLimit == 0.0 && !refinementFunc) return(0);
1156:   DMGetDimension(dm, &dim);
1157:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1158:   PetscOptionsGetString(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_generator", genname, 1024, &flg);
1159:   if (flg) name = genname;
1160:   if (name) {
1161:     PetscStrcmp(name, "triangle", &isTriangle);
1162:     PetscStrcmp(name, "tetgen",   &isTetgen);
1163:     PetscStrcmp(name, "ctetgen",  &isCTetgen);
1164:   }
1165:   switch (dim) {
1166:   case 2:
1167:     if (!name || isTriangle) {
1168: #if defined(PETSC_HAVE_TRIANGLE)
1169:       double  *maxVolumes;
1170:       PetscInt c;

1172:       PetscMalloc1(cEnd - cStart, &maxVolumes);
1173:       if (refinementFunc) {
1174:         for (c = cStart; c < cEnd; ++c) {
1175:           PetscReal vol, centroid[3];
1176:           PetscReal maxVol;

1178:           DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);
1179:           (*refinementFunc)(centroid, &maxVol);
1180:           maxVolumes[c - cStart] = (double) maxVol;
1181:         }
1182:       } else {
1183:         for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1184:       }
1185:       DMPlexRefine_Triangle(dm, maxVolumes, dmRefined);
1186:       PetscFree(maxVolumes);
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;

1198:       PetscMalloc1(cEnd - cStart, &maxVolumes);
1199:       if (refinementFunc) {
1200:         for (c = cStart; c < cEnd; ++c) {
1201:           PetscReal vol, centroid[3];

1203:           DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);
1204:           (*refinementFunc)(centroid, &maxVolumes[c-cStart]);
1205:         }
1206:       } else {
1207:         for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1208:       }
1209:       DMPlexRefine_CTetgen(dm, maxVolumes, dmRefined);
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;

1218:       PetscMalloc1(cEnd - cStart, &maxVolumes);
1219:       if (refinementFunc) {
1220:         for (c = cStart; c < cEnd; ++c) {
1221:           PetscReal vol, centroid[3];

1223:           DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);
1224:           (*refinementFunc)(centroid, &maxVolumes[c-cStart]);
1225:         }
1226:       } else {
1227:         for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1228:       }
1229:       DMPlexRefine_Tetgen(dm, maxVolumes, dmRefined);
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:   DMCopyBoundary(dm, *dmRefined);
1239:   if (localized) {
1240:     DMLocalizeCoordinates(*dmRefined);
1241:   }
1242:   return(0);
1243: }

1247: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[])
1248: {
1249:   DM             cdm = dm;
1250:   PetscInt       r;
1251:   PetscBool      isUniform, localized;

1255:   DMPlexGetRefinementUniform(dm, &isUniform);
1256:   DMGetCoordinatesLocalized(dm, &localized);
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;

1261:     DMPlexGetCellRefiner_Internal(cdm, &cellRefiner);
1262:     DMPlexRefineUniform_Internal(cdm, cellRefiner, &dmRefined[r]);
1263:     DMCopyBoundary(cdm, dmRefined[r]);
1264:     if (localized) {
1265:       DMLocalizeCoordinates(dmRefined[r]);
1266:     }
1267:     DMSetCoarseDM(dmRefined[r], cdm);
1268:     DMPlexSetRegularRefinement(dmRefined[r], PETSC_TRUE);
1269:     cdm  = dmRefined[r];
1270:   }
1271:   return(0);
1272: }