Actual source code: ctetgenerate.c

      1: #include <petsc/private/dmpleximpl.h>

      3: #ifdef PETSC_HAVE_EGADS
      4:   #include <egads.h>
      5: #endif

      7: #include <ctetgen.h>

      9: /* This is to fix the tetrahedron orientation from TetGen */
     10: static PetscErrorCode DMPlexInvertCells_CTetgen(PetscInt numCells, PetscInt numCorners, PetscInt cells[])
     11: {
     12:   PetscInt bound = numCells * numCorners, coff;

     14:   PetscFunctionBegin;
     15: #define SWAP(a, b) \
     16:   do { \
     17:     PetscInt tmp = (a); \
     18:     (a)          = (b); \
     19:     (b)          = tmp; \
     20:   } while (0)
     21:   for (coff = 0; coff < bound; coff += numCorners) SWAP(cells[coff], cells[coff + 1]);
     22: #undef SWAP
     23:   PetscFunctionReturn(PETSC_SUCCESS);
     24: }

     26: PETSC_EXTERN PetscErrorCode DMPlexGenerate_CTetgen(DM boundary, PetscBool interpolate, DM *dm)
     27: {
     28:   MPI_Comm               comm;
     29:   const PetscInt         dim = 3;
     30:   PLC                   *in, *out;
     31:   DMUniversalLabel       universal;
     32:   PetscInt               vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, defVal, verbose = 0;
     33:   DMPlexInterpolatedFlag isInterpolated;
     34:   PetscMPIInt            rank;

     36:   PetscFunctionBegin;
     37:   PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)boundary)->prefix, "-ctetgen_verbose", &verbose, NULL));
     38:   PetscCall(PetscObjectGetComm((PetscObject)boundary, &comm));
     39:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
     40:   PetscCall(DMPlexIsInterpolatedCollective(boundary, &isInterpolated));
     41:   PetscCall(DMUniversalLabelCreate(boundary, &universal));
     42:   PetscCall(DMLabelGetDefaultValue(universal->label, &defVal));

     44:   PetscCall(PLCCreate(&in));
     45:   PetscCall(PLCCreate(&out));

     47:   PetscCall(DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd));
     48:   in->numberofpoints = vEnd - vStart;
     49:   if (in->numberofpoints > 0) {
     50:     PetscSection       coordSection;
     51:     Vec                coordinates;
     52:     const PetscScalar *array;

     54:     PetscCall(PetscMalloc1(in->numberofpoints * dim, &in->pointlist));
     55:     PetscCall(PetscCalloc1(in->numberofpoints, &in->pointmarkerlist));
     56:     PetscCall(DMGetCoordinatesLocal(boundary, &coordinates));
     57:     PetscCall(DMGetCoordinateSection(boundary, &coordSection));
     58:     PetscCall(VecGetArrayRead(coordinates, &array));
     59:     for (v = vStart; v < vEnd; ++v) {
     60:       const PetscInt idx = v - vStart;
     61:       PetscInt       off, d, m;

     63:       PetscCall(PetscSectionGetOffset(coordSection, v, &off));
     64:       for (d = 0; d < dim; ++d) in->pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
     65:       PetscCall(DMLabelGetValue(universal->label, v, &m));
     66:       if (m != defVal) in->pointmarkerlist[idx] = (int)m;
     67:     }
     68:     PetscCall(VecRestoreArrayRead(coordinates, &array));
     69:   }

     71:   PetscCall(DMPlexGetHeightStratum(boundary, 1, &eStart, &eEnd));
     72:   in->numberofedges = eEnd - eStart;
     73:   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberofedges > 0) {
     74:     PetscCall(PetscMalloc1(in->numberofedges * 2, &in->edgelist));
     75:     PetscCall(PetscMalloc1(in->numberofedges, &in->edgemarkerlist));
     76:     for (e = eStart; e < eEnd; ++e) {
     77:       const PetscInt  idx = e - eStart;
     78:       const PetscInt *cone;
     79:       PetscInt        coneSize, val;

     81:       PetscCall(DMPlexGetConeSize(boundary, e, &coneSize));
     82:       PetscCall(DMPlexGetCone(boundary, e, &cone));
     83:       in->edgelist[idx * 2]     = cone[0] - vStart;
     84:       in->edgelist[idx * 2 + 1] = cone[1] - vStart;

     86:       PetscCall(DMLabelGetValue(universal->label, e, &val));
     87:       if (val != defVal) in->edgemarkerlist[idx] = (int)val;
     88:     }
     89:   }

     91:   PetscCall(DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd));
     92:   in->numberoffacets = fEnd - fStart;
     93:   if (in->numberoffacets > 0) {
     94:     PetscCall(PetscMalloc1(in->numberoffacets, &in->facetlist));
     95:     PetscCall(PetscMalloc1(in->numberoffacets, &in->facetmarkerlist));
     96:     for (f = fStart; f < fEnd; ++f) {
     97:       const PetscInt idx    = f - fStart;
     98:       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v, m = -1;
     99:       polygon       *poly;

    101:       in->facetlist[idx].numberofpolygons = 1;
    102:       PetscCall(PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist));
    103:       in->facetlist[idx].numberofholes = 0;
    104:       in->facetlist[idx].holelist      = NULL;

    106:       PetscCall(DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points));
    107:       for (p = 0; p < numPoints * 2; p += 2) {
    108:         const PetscInt point = points[p];
    109:         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
    110:       }

    112:       poly                   = in->facetlist[idx].polygonlist;
    113:       poly->numberofvertices = numVertices;
    114:       PetscCall(PetscMalloc1(poly->numberofvertices, &poly->vertexlist));
    115:       for (v = 0; v < numVertices; ++v) {
    116:         const PetscInt vIdx = points[v] - vStart;
    117:         poly->vertexlist[v] = vIdx;
    118:       }
    119:       PetscCall(DMLabelGetValue(universal->label, f, &m));
    120:       if (m != defVal) in->facetmarkerlist[idx] = (int)m;
    121:       PetscCall(DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points));
    122:     }
    123:   }
    124:   if (rank == 0) {
    125:     TetGenOpts t;

    127:     PetscCall(TetGenOptsInitialize(&t));
    128:     t.in        = boundary; /* Should go away */
    129:     t.plc       = 1;
    130:     t.quality   = 1;
    131:     t.edgesout  = 1;
    132:     t.zeroindex = 1;
    133:     t.quiet     = 1;
    134:     t.verbose   = verbose;
    135: #if 0
    136:   #ifdef PETSC_HAVE_EGADS
    137:     /* Need to add in more TetGen code */
    138:     t.nobisect  = 1; /* Add Y to preserve Surface Mesh for EGADS */
    139:   #endif
    140: #endif

    142:     PetscCall(TetGenCheckOpts(&t));
    143:     PetscCall(TetGenTetrahedralize(&t, in, out));
    144:   }
    145:   {
    146:     const PetscInt numCorners  = 4;
    147:     const PetscInt numCells    = out->numberoftetrahedra;
    148:     const PetscInt numVertices = out->numberofpoints;
    149:     PetscReal     *meshCoords  = NULL;
    150:     PetscInt      *cells       = NULL;

    152:     if (sizeof(PetscReal) == sizeof(out->pointlist[0])) {
    153:       meshCoords = (PetscReal *)out->pointlist;
    154:     } else {
    155:       PetscInt i;

    157:       PetscCall(PetscMalloc1(dim * numVertices, &meshCoords));
    158:       for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal)out->pointlist[i];
    159:     }
    160:     if (sizeof(PetscInt) == sizeof(out->tetrahedronlist[0])) {
    161:       cells = (PetscInt *)out->tetrahedronlist;
    162:     } else {
    163:       PetscInt i;

    165:       PetscCall(PetscMalloc1(numCells * numCorners, &cells));
    166:       for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out->tetrahedronlist[i];
    167:     }

    169:     PetscCall(DMPlexInvertCells_CTetgen(numCells, numCorners, cells));
    170:     PetscCall(DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm));
    171:     if (sizeof(PetscReal) != sizeof(out->pointlist[0])) PetscCall(PetscFree(meshCoords));
    172:     if (sizeof(PetscInt) != sizeof(out->tetrahedronlist[0])) PetscCall(PetscFree(cells));

    174:     /* Set labels */
    175:     PetscCall(DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dm));
    176:     for (v = 0; v < numVertices; ++v) {
    177:       if (out->pointmarkerlist[v]) PetscCall(DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, v + numCells, out->pointmarkerlist[v]));
    178:     }
    179:     if (interpolate) {
    180:       PetscInt e;

    182:       for (e = 0; e < out->numberofedges; e++) {
    183:         if (out->edgemarkerlist[e]) {
    184:           const PetscInt  vertices[2] = {out->edgelist[e * 2 + 0] + numCells, out->edgelist[e * 2 + 1] + numCells};
    185:           const PetscInt *edges;
    186:           PetscInt        numEdges;

    188:           PetscCall(DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges));
    189:           PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %" PetscInt_FMT, numEdges);
    190:           PetscCall(DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, edges[0], out->edgemarkerlist[e]));
    191:           PetscCall(DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges));
    192:         }
    193:       }
    194:       for (f = 0; f < out->numberoftrifaces; f++) {
    195:         if (out->trifacemarkerlist[f]) {
    196:           const PetscInt  vertices[3] = {out->trifacelist[f * 3 + 0] + numCells, out->trifacelist[f * 3 + 1] + numCells, out->trifacelist[f * 3 + 2] + numCells};
    197:           const PetscInt *faces;
    198:           PetscInt        numFaces;

    200:           PetscCall(DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces));
    201:           PetscCheck(numFaces == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %" PetscInt_FMT, numFaces);
    202:           PetscCall(DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, faces[0], out->trifacemarkerlist[f]));
    203:           PetscCall(DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces));
    204:         }
    205:       }
    206:     }

    208: #ifdef PETSC_HAVE_EGADS
    209:     {
    210:       DMLabel        bodyLabel;
    211:       PetscContainer modelObj;
    212:       PetscInt       cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
    213:       ego           *bodies;
    214:       ego            model, geom;
    215:       int            Nb, oclass, mtype, *senses;

    217:       /* Get Attached EGADS Model from Original DMPlex */
    218:       PetscCall(PetscObjectQuery((PetscObject)boundary, "EGADS Model", (PetscObject *)&modelObj));
    219:       if (modelObj) {
    220:         PetscCall(PetscContainerGetPointer(modelObj, (void **)&model));
    221:         PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses));
    222:         /* Transfer EGADS Model to Volumetric Mesh */
    223:         PetscCall(PetscObjectCompose((PetscObject)*dm, "EGADS Model", (PetscObject)modelObj));

    225:         /* Set Cell Labels */
    226:         PetscCall(DMGetLabel(*dm, "EGADS Body ID", &bodyLabel));
    227:         PetscCall(DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd));
    228:         PetscCall(DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd));
    229:         PetscCall(DMPlexGetDepthStratum(*dm, 1, &eStart, &eEnd));

    231:         for (c = cStart; c < cEnd; ++c) {
    232:           PetscReal centroid[3] = {0., 0., 0.};
    233:           PetscInt  b;

    235:           /* Determine what body the cell's centroid is located in */
    236:           if (!interpolate) {
    237:             PetscSection coordSection;
    238:             Vec          coordinates;
    239:             PetscScalar *coords = NULL;
    240:             PetscInt     coordSize, s, d;

    242:             PetscCall(DMGetCoordinatesLocal(*dm, &coordinates));
    243:             PetscCall(DMGetCoordinateSection(*dm, &coordSection));
    244:             PetscCall(DMPlexVecGetClosure(*dm, coordSection, coordinates, c, &coordSize, &coords));
    245:             for (s = 0; s < coordSize; ++s)
    246:               for (d = 0; d < dim; ++d) centroid[d] += coords[s * dim + d];
    247:             PetscCall(DMPlexVecRestoreClosure(*dm, coordSection, coordinates, c, &coordSize, &coords));
    248:           } else PetscCall(DMPlexComputeCellGeometryFVM(*dm, c, NULL, centroid, NULL));
    249:           for (b = 0; b < Nb; ++b) {
    250:             if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
    251:           }
    252:           if (b < Nb) {
    253:             PetscInt  cval    = b, eVal, fVal;
    254:             PetscInt *closure = NULL, Ncl, cl;

    256:             PetscCall(DMLabelSetValue(bodyLabel, c, cval));
    257:             PetscCall(DMPlexGetTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure));
    258:             for (cl = 0; cl < Ncl; ++cl) {
    259:               const PetscInt p = closure[cl * 2];

    261:               if (p >= eStart && p < eEnd) {
    262:                 PetscCall(DMLabelGetValue(bodyLabel, p, &eVal));
    263:                 if (eVal < 0) PetscCall(DMLabelSetValue(bodyLabel, p, cval));
    264:               }
    265:               if (p >= fStart && p < fEnd) {
    266:                 PetscCall(DMLabelGetValue(bodyLabel, p, &fVal));
    267:                 if (fVal < 0) PetscCall(DMLabelSetValue(bodyLabel, p, cval));
    268:               }
    269:             }
    270:             PetscCall(DMPlexRestoreTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure));
    271:           }
    272:         }
    273:       }
    274:     }
    275: #endif
    276:     PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_FALSE));
    277:   }

    279:   PetscCall(DMUniversalLabelDestroy(&universal));
    280:   PetscCall(PLCDestroy(&in));
    281:   PetscCall(PLCDestroy(&out));
    282:   PetscFunctionReturn(PETSC_SUCCESS);
    283: }

    285: PETSC_EXTERN PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
    286: {
    287:   MPI_Comm               comm;
    288:   const PetscInt         dim = 3;
    289:   PLC                   *in, *out;
    290:   DMUniversalLabel       universal;
    291:   PetscInt               vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, cStart, cEnd, c, defVal, verbose = 0;
    292:   DMPlexInterpolatedFlag isInterpolated;
    293:   PetscMPIInt            rank;

    295:   PetscFunctionBegin;
    296:   PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)dm)->prefix, "-ctetgen_verbose", &verbose, NULL));
    297:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    298:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    299:   PetscCall(DMPlexIsInterpolatedCollective(dm, &isInterpolated));
    300:   PetscCall(DMUniversalLabelCreate(dm, &universal));
    301:   PetscCall(DMLabelGetDefaultValue(universal->label, &defVal));

    303:   PetscCall(PLCCreate(&in));
    304:   PetscCall(PLCCreate(&out));

    306:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
    307:   in->numberofpoints = vEnd - vStart;
    308:   if (in->numberofpoints > 0) {
    309:     PetscSection coordSection;
    310:     Vec          coordinates;
    311:     PetscScalar *array;

    313:     PetscCall(PetscMalloc1(in->numberofpoints * dim, &in->pointlist));
    314:     PetscCall(PetscCalloc1(in->numberofpoints, &in->pointmarkerlist));
    315:     PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
    316:     PetscCall(DMGetCoordinateSection(dm, &coordSection));
    317:     PetscCall(VecGetArray(coordinates, &array));
    318:     for (v = vStart; v < vEnd; ++v) {
    319:       const PetscInt idx = v - vStart;
    320:       PetscInt       off, d, m;

    322:       PetscCall(PetscSectionGetOffset(coordSection, v, &off));
    323:       for (d = 0; d < dim; ++d) in->pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
    324:       PetscCall(DMLabelGetValue(universal->label, v, &m));
    325:       if (m != defVal) in->pointmarkerlist[idx] = (int)m;
    326:     }
    327:     PetscCall(VecRestoreArray(coordinates, &array));
    328:   }

    330:   PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
    331:   in->numberofedges = eEnd - eStart;
    332:   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberofedges > 0) {
    333:     PetscCall(PetscMalloc1(in->numberofedges * 2, &in->edgelist));
    334:     PetscCall(PetscMalloc1(in->numberofedges, &in->edgemarkerlist));
    335:     for (e = eStart; e < eEnd; ++e) {
    336:       const PetscInt  idx = e - eStart;
    337:       const PetscInt *cone;
    338:       PetscInt        coneSize, val;

    340:       PetscCall(DMPlexGetConeSize(dm, e, &coneSize));
    341:       PetscCall(DMPlexGetCone(dm, e, &cone));
    342:       in->edgelist[idx * 2]     = cone[0] - vStart;
    343:       in->edgelist[idx * 2 + 1] = cone[1] - vStart;

    345:       PetscCall(DMLabelGetValue(universal->label, e, &val));
    346:       if (val != defVal) in->edgemarkerlist[idx] = (int)val;
    347:     }
    348:   }

    350:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
    351:   in->numberoftrifaces = 0;
    352:   for (f = fStart; f < fEnd; ++f) {
    353:     PetscInt supportSize;

    355:     PetscCall(DMPlexGetSupportSize(dm, f, &supportSize));
    356:     if (supportSize == 1) ++in->numberoftrifaces;
    357:   }
    358:   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberoftrifaces > 0) {
    359:     PetscInt tf = 0;

    361:     PetscCall(PetscMalloc1(in->numberoftrifaces * 3, &in->trifacelist));
    362:     PetscCall(PetscMalloc1(in->numberoftrifaces, &in->trifacemarkerlist));
    363:     for (f = fStart; f < fEnd; ++f) {
    364:       PetscInt *points = NULL;
    365:       PetscInt  supportSize, numPoints, p, Nv = 0, val;

    367:       PetscCall(DMPlexGetSupportSize(dm, f, &supportSize));
    368:       if (supportSize != 1) continue;
    369:       PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points));
    370:       for (p = 0; p < numPoints * 2; p += 2) {
    371:         const PetscInt point = points[p];
    372:         if ((point >= vStart) && (point < vEnd)) in->trifacelist[tf * 3 + Nv++] = point - vStart;
    373:       }
    374:       PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points));
    375:       PetscCheck(Nv == 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " has %" PetscInt_FMT " vertices, not 3", f, Nv);
    376:       PetscCall(DMLabelGetValue(universal->label, f, &val));
    377:       if (val != defVal) in->trifacemarkerlist[tf] = (int)val;
    378:       ++tf;
    379:     }
    380:   }

    382:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
    383:   in->numberofcorners       = 4;
    384:   in->numberoftetrahedra    = cEnd - cStart;
    385:   in->tetrahedronvolumelist = maxVolumes;
    386:   if (in->numberoftetrahedra > 0) {
    387:     PetscCall(PetscMalloc1(in->numberoftetrahedra * in->numberofcorners, &in->tetrahedronlist));
    388:     for (c = cStart; c < cEnd; ++c) {
    389:       const PetscInt idx     = c - cStart;
    390:       PetscInt      *closure = NULL;
    391:       PetscInt       closureSize;

    393:       PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
    394:       PetscCheck((closureSize == 5) || (closureSize == 15), comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %" PetscInt_FMT " vertices in closure", closureSize);
    395:       for (v = 0; v < 4; ++v) in->tetrahedronlist[idx * in->numberofcorners + v] = closure[(v + closureSize - 4) * 2] - vStart;
    396:       PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
    397:     }
    398:   }

    400:   if (rank == 0) {
    401:     TetGenOpts t;

    403:     PetscCall(TetGenOptsInitialize(&t));

    405:     t.in        = dm; /* Should go away */
    406:     t.refine    = 1;
    407:     t.varvolume = 1;
    408:     t.quality   = 1;
    409:     t.edgesout  = 1;
    410:     t.zeroindex = 1;
    411:     t.quiet     = 1;
    412:     t.verbose   = verbose; /* Change this */

    414:     PetscCall(TetGenCheckOpts(&t));
    415:     PetscCall(TetGenTetrahedralize(&t, in, out));
    416:   }

    418:   in->tetrahedronvolumelist = NULL;
    419:   {
    420:     const PetscInt numCorners  = 4;
    421:     const PetscInt numCells    = out->numberoftetrahedra;
    422:     const PetscInt numVertices = out->numberofpoints;
    423:     PetscReal     *meshCoords  = NULL;
    424:     PetscInt      *cells       = NULL;
    425:     PetscBool      interpolate = isInterpolated == DMPLEX_INTERPOLATED_FULL ? PETSC_TRUE : PETSC_FALSE;

    427:     if (sizeof(PetscReal) == sizeof(out->pointlist[0])) {
    428:       meshCoords = (PetscReal *)out->pointlist;
    429:     } else {
    430:       PetscInt i;

    432:       PetscCall(PetscMalloc1(dim * numVertices, &meshCoords));
    433:       for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal)out->pointlist[i];
    434:     }
    435:     if (sizeof(PetscInt) == sizeof(out->tetrahedronlist[0])) {
    436:       cells = (PetscInt *)out->tetrahedronlist;
    437:     } else {
    438:       PetscInt i;

    440:       PetscCall(PetscMalloc1(numCells * numCorners, &cells));
    441:       for (i = 0; i < numCells * numCorners; ++i) cells[i] = (PetscInt)out->tetrahedronlist[i];
    442:     }

    444:     PetscCall(DMPlexInvertCells_CTetgen(numCells, numCorners, cells));
    445:     PetscCall(DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined));
    446:     if (sizeof(PetscReal) != sizeof(out->pointlist[0])) PetscCall(PetscFree(meshCoords));
    447:     if (sizeof(PetscInt) != sizeof(out->tetrahedronlist[0])) PetscCall(PetscFree(cells));

    449:     /* Set labels */
    450:     PetscCall(DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dmRefined));
    451:     for (v = 0; v < numVertices; ++v) {
    452:       if (out->pointmarkerlist[v]) PetscCall(DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, v + numCells, out->pointmarkerlist[v]));
    453:     }
    454:     if (interpolate) {
    455:       PetscInt e, f;

    457:       for (e = 0; e < out->numberofedges; e++) {
    458:         if (out->edgemarkerlist[e]) {
    459:           const PetscInt  vertices[2] = {out->edgelist[e * 2 + 0] + numCells, out->edgelist[e * 2 + 1] + numCells};
    460:           const PetscInt *edges;
    461:           PetscInt        numEdges;

    463:           PetscCall(DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges));
    464:           PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %" PetscInt_FMT, numEdges);
    465:           PetscCall(DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out->edgemarkerlist[e]));
    466:           PetscCall(DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges));
    467:         }
    468:       }
    469:       for (f = 0; f < out->numberoftrifaces; f++) {
    470:         if (out->trifacemarkerlist[f]) {
    471:           const PetscInt  vertices[3] = {out->trifacelist[f * 3 + 0] + numCells, out->trifacelist[f * 3 + 1] + numCells, out->trifacelist[f * 3 + 2] + numCells};
    472:           const PetscInt *faces;
    473:           PetscInt        numFaces;

    475:           PetscCall(DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces));
    476:           PetscCheck(numFaces == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %" PetscInt_FMT, numFaces);
    477:           PetscCall(DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out->trifacemarkerlist[f]));
    478:           PetscCall(DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces));
    479:         }
    480:       }
    481:     }

    483: #ifdef PETSC_HAVE_EGADS
    484:     {
    485:       DMLabel        bodyLabel;
    486:       PetscContainer modelObj;
    487:       PetscInt       cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
    488:       ego           *bodies;
    489:       ego            model, geom;
    490:       int            Nb, oclass, mtype, *senses;

    492:       /* Get Attached EGADS Model from Original DMPlex */
    493:       PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj));
    494:       if (modelObj) {
    495:         PetscCall(PetscContainerGetPointer(modelObj, (void **)&model));
    496:         PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses));
    497:         /* Transfer EGADS Model to Volumetric Mesh */
    498:         PetscCall(PetscObjectCompose((PetscObject)*dmRefined, "EGADS Model", (PetscObject)modelObj));

    500:         /* Set Cell Labels */
    501:         PetscCall(DMGetLabel(*dmRefined, "EGADS Body ID", &bodyLabel));
    502:         PetscCall(DMPlexGetHeightStratum(*dmRefined, 0, &cStart, &cEnd));
    503:         PetscCall(DMPlexGetHeightStratum(*dmRefined, 1, &fStart, &fEnd));
    504:         PetscCall(DMPlexGetDepthStratum(*dmRefined, 1, &eStart, &eEnd));

    506:         for (c = cStart; c < cEnd; ++c) {
    507:           PetscReal centroid[3] = {0., 0., 0.};
    508:           PetscInt  b;

    510:           /* Determine what body the cell's centroid is located in */
    511:           if (!interpolate) {
    512:             PetscSection coordSection;
    513:             Vec          coordinates;
    514:             PetscScalar *coords = NULL;
    515:             PetscInt     coordSize, s, d;

    517:             PetscCall(DMGetCoordinatesLocal(*dmRefined, &coordinates));
    518:             PetscCall(DMGetCoordinateSection(*dmRefined, &coordSection));
    519:             PetscCall(DMPlexVecGetClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords));
    520:             for (s = 0; s < coordSize; ++s)
    521:               for (d = 0; d < dim; ++d) centroid[d] += coords[s * dim + d];
    522:             PetscCall(DMPlexVecRestoreClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords));
    523:           } else PetscCall(DMPlexComputeCellGeometryFVM(*dmRefined, c, NULL, centroid, NULL));
    524:           for (b = 0; b < Nb; ++b) {
    525:             if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
    526:           }
    527:           if (b < Nb) {
    528:             PetscInt  cval    = b, eVal, fVal;
    529:             PetscInt *closure = NULL, Ncl, cl;

    531:             PetscCall(DMLabelSetValue(bodyLabel, c, cval));
    532:             PetscCall(DMPlexGetTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure));
    533:             for (cl = 0; cl < Ncl; cl += 2) {
    534:               const PetscInt p = closure[cl];

    536:               if (p >= eStart && p < eEnd) {
    537:                 PetscCall(DMLabelGetValue(bodyLabel, p, &eVal));
    538:                 if (eVal < 0) PetscCall(DMLabelSetValue(bodyLabel, p, cval));
    539:               }
    540:               if (p >= fStart && p < fEnd) {
    541:                 PetscCall(DMLabelGetValue(bodyLabel, p, &fVal));
    542:                 if (fVal < 0) PetscCall(DMLabelSetValue(bodyLabel, p, cval));
    543:               }
    544:             }
    545:             PetscCall(DMPlexRestoreTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure));
    546:           }
    547:         }
    548:       }
    549:     }
    550: #endif
    551:     PetscCall(DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE));
    552:   }
    553:   PetscCall(DMUniversalLabelDestroy(&universal));
    554:   PetscCall(PLCDestroy(&in));
    555:   PetscCall(PLCDestroy(&out));
    556:   PetscFunctionReturn(PETSC_SUCCESS);
    557: }