Actual source code: plexgenerate.c

petsc-3.6.1 2015-08-06
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:   DMPlexGetLabel(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) {DMPlexCreateLabel(*dm, labelName); DMPlexGetLabel(*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:   MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
355:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
356:   DMPlexGetLabel(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) {DMPlexCreateLabel(*dmRefined, labelName); DMPlexGetLabel(*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:   DMPlexGetLabel(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) {DMPlexCreateLabel(*dm, labelName); DMPlexGetLabel(*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:   MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
627:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
628:   DMPlexGetLabel(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) {DMPlexCreateLabel(*dmRefined, labelName); DMPlexGetLabel(*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(((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);
755:   MPI_Comm_rank(comm, &rank);
756:   DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
757:   DMPlexGetLabel(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:     const double   *meshCoords = out->pointlist;
842:     int            *cells      = out->tetrahedronlist;

844:     DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
845:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
846:     if (label) {DMPlexCreateLabel(*dm, labelName); DMPlexGetLabel(*dm, labelName, &glabel);}
847:     /* Set labels */
848:     for (v = 0; v < numVertices; ++v) {
849:       if (out->pointmarkerlist[v]) {
850:         if (glabel) {DMLabelSetValue(glabel, v+numCells, out->pointmarkerlist[v]);}
851:       }
852:     }
853:     if (interpolate) {
854:       PetscInt e;

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;

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

874:           DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
875:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
876:           if (glabel) {DMLabelSetValue(glabel, faces[0], out->trifacemarkerlist[f]);}
877:           DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
878:         }
879:       }
880:     }
881:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
882:   }

884:   PLCDestroy(&in);
885:   PLCDestroy(&out);
886:   return(0);
887: }

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;

903:   PetscObjectGetComm((PetscObject)dm,&comm);
904:   PetscOptionsGetInt(((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);
905:   MPI_Comm_rank(comm, &rank);
906:   DMPlexGetDepth(dm, &depth);
907:   MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
908:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
909:   DMPlexGetLabel(dm, labelName, &label);
910:   PLCCreate(&in);
911:   PLCCreate(&out);

913:   in->numberofpoints = vEnd - vStart;
914:   if (in->numberofpoints > 0) {
915:     PetscSection coordSection;
916:     Vec          coordinates;
917:     PetscScalar *array;

919:     PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
920:     PetscMalloc1(in->numberofpoints,       &in->pointmarkerlist);
921:     DMGetCoordinatesLocal(dm, &coordinates);
922:     DMGetCoordinateSection(dm, &coordSection);
923:     VecGetArray(coordinates, &array);
924:     for (v = vStart; v < vEnd; ++v) {
925:       const PetscInt idx = v - vStart;
926:       PetscInt       off, d, m = -1;

928:       PetscSectionGetOffset(coordSection, v, &off);
929:       for (d = 0; d < dim; ++d) {
930:         in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
931:       }
932:       if (label) {DMLabelGetValue(label, v, &m);}

934:       in->pointmarkerlist[idx] = (int) m;
935:     }
936:     VecRestoreArray(coordinates, &array);
937:   }
938:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);

940:   in->numberofcorners       = 4;
941:   in->numberoftetrahedra    = cEnd - cStart;
942:   in->tetrahedronvolumelist = maxVolumes;
943:   if (in->numberoftetrahedra > 0) {
944:     PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);
945:     for (c = cStart; c < cEnd; ++c) {
946:       const PetscInt idx      = c - cStart;
947:       PetscInt      *closure = NULL;
948:       PetscInt       closureSize;

950:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
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:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
956:     }
957:   }
958:   if (!rank) {
959:     TetGenOpts t;

961:     TetGenOptsInitialize(&t);

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

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

984:     DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
985:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
986:     if (label) {DMPlexCreateLabel(*dmRefined, labelName); DMPlexGetLabel(*dmRefined, labelName, &rlabel);}
987:     /* Set labels */
988:     for (v = 0; v < numVertices; ++v) {
989:       if (out->pointmarkerlist[v]) {
990:         if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out->pointmarkerlist[v]);}
991:       }
992:     }
993:     if (interpolate) {
994:       PetscInt e, f;

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;

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

1014:           DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
1015:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
1016:           if (rlabel) {DMLabelSetValue(rlabel, faces[0], out->trifacemarkerlist[f]);}
1017:           DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
1018:         }
1019:       }
1020:     }
1021:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
1022:   }
1023:   PLCDestroy(&in);
1024:   PLCDestroy(&out);
1025:   return(0);
1026: }
1027: #endif /* PETSC_HAVE_CTETGEN */

1031: /*@C
1032:   DMPlexGenerate - Generates a mesh.

1034:   Not Collective

1036:   Input Parameters:
1037: + boundary - The DMPlex boundary object
1038: . name - The mesh generation package name
1039: - interpolate - Flag to create intermediate mesh elements

1041:   Output Parameter:
1042: . mesh - The DMPlex object

1044:   Level: intermediate

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;

1059:   DMGetDimension(boundary, &dim);
1060:   PetscOptionsGetString(((PetscObject) boundary)->prefix, "-dm_plex_generator", genname, 1024, &flg);
1061:   if (flg) name = genname;
1062:   if (name) {
1063:     PetscStrcmp(name, "triangle", &isTriangle);
1064:     PetscStrcmp(name, "tetgen",   &isTetgen);
1065:     PetscStrcmp(name, "ctetgen",  &isCTetgen);
1066:   }
1067:   switch (dim) {
1068:   case 1:
1069:     if (!name || isTriangle) {
1070: #if defined(PETSC_HAVE_TRIANGLE)
1071:       DMPlexGenerate_Triangle(boundary, interpolate, mesh);
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:       DMPlexGenerate_CTetgen(boundary, interpolate, mesh);
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:       DMPlexGenerate_Tetgen(boundary, interpolate, mesh);
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:   return(0);
1096: }

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;

1109:   DMPlexGetRefinementUniform(dm, &isUniform);
1110:   if (isUniform) {
1111:     CellRefiner cellRefiner;

1113:     DMPlexGetCellRefiner_Internal(dm, &cellRefiner);
1114:     DMPlexRefineUniform_Internal(dm, cellRefiner, dmRefined);
1115:     DMPlexCopyBoundary(dm, *dmRefined);
1116:     return(0);
1117:   }
1118:   DMPlexGetRefinementLimit(dm, &refinementLimit);
1119:   if (refinementLimit == 0.0) return(0);
1120:   DMGetDimension(dm, &dim);
1121:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1122:   PetscOptionsGetString(((PetscObject) dm)->prefix, "-dm_plex_generator", genname, 1024, &flg);
1123:   if (flg) name = genname;
1124:   if (name) {
1125:     PetscStrcmp(name, "triangle", &isTriangle);
1126:     PetscStrcmp(name, "tetgen",   &isTetgen);
1127:     PetscStrcmp(name, "ctetgen",  &isCTetgen);
1128:   }
1129:   switch (dim) {
1130:   case 2:
1131:     if (!name || isTriangle) {
1132: #if defined(PETSC_HAVE_TRIANGLE)
1133:       double  *maxVolumes;
1134:       PetscInt c;

1136:       PetscMalloc1(cEnd - cStart, &maxVolumes);
1137:       for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1138:       DMPlexRefine_Triangle(dm, maxVolumes, dmRefined);
1139:       PetscFree(maxVolumes);
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;

1151:       PetscMalloc1(cEnd - cStart, &maxVolumes);
1152:       for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1153:       DMPlexRefine_CTetgen(dm, maxVolumes, dmRefined);
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;

1162:       PetscMalloc1(cEnd - cStart, &maxVolumes);
1163:       for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1164:       DMPlexRefine_Tetgen(dm, maxVolumes, dmRefined);
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:   DMPlexCopyBoundary(dm, *dmRefined);
1174:   return(0);
1175: }

1179: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[])
1180: {
1181:   DM             cdm = dm;
1182:   PetscInt       r;
1183:   PetscBool      isUniform;

1187:   DMPlexGetRefinementUniform(dm, &isUniform);
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;

1192:     DMPlexGetCellRefiner_Internal(cdm, &cellRefiner);
1193:     DMPlexRefineUniform_Internal(cdm, cellRefiner, &dmRefined[r]);
1194:     DMPlexCopyBoundary(cdm, dmRefined[r]);
1195:     DMPlexSetCoarseDM(dmRefined[r], cdm);
1196:     cdm  = dmRefined[r];
1197:   }
1198:   return(0);
1199: }

1203: PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened)
1204: {
1205:   DM_Plex       *mesh = (DM_Plex*) dm->data;

1209:   PetscObjectReference((PetscObject) mesh->coarseMesh);
1210:   *dmCoarsened = mesh->coarseMesh;
1211:   return(0);
1212: }