Actual source code: plexsubmesh.c

      1: #include <petsc/private/dmpleximpl.h>
      2: #include <petsc/private/dmlabelimpl.h>
      3: #include <petscsf.h>

      5: static PetscErrorCode DMPlexCellIsHybrid_Internal(DM dm, PetscInt p, PetscBool *isHybrid)
      6: {
      7:   DMPolytopeType ct;

     11:   DMPlexGetCellType(dm, p, &ct);
     12:   switch (ct) {
     13:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
     14:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
     15:     case DM_POLYTOPE_TRI_PRISM_TENSOR:
     16:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
     17:       *isHybrid = PETSC_TRUE;
     18:     default: *isHybrid = PETSC_FALSE;
     19:   }
     20:   return(0);
     21: }

     23: static PetscErrorCode DMPlexGetTensorPrismBounds_Internal(DM dm, PetscInt dim, PetscInt *cStart, PetscInt *cEnd)
     24: {
     25:   DMLabel        ctLabel;

     29:   if (cStart) *cStart = -1;
     30:   if (cEnd)   *cEnd   = -1;
     31:   DMPlexGetCellTypeLabel(dm, &ctLabel);
     32:   switch (dim) {
     33:     case 1: DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_POINT_PRISM_TENSOR, cStart, cEnd);break;
     34:     case 2: DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_SEG_PRISM_TENSOR, cStart, cEnd);break;
     35:     case 3:
     36:       DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_TRI_PRISM_TENSOR, cStart, cEnd);
     37:       if (*cStart < 0) {DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_QUAD_PRISM_TENSOR, cStart, cEnd);}
     38:       break;
     39:     default: return(0);
     40:   }
     41:   return(0);
     42: }

     44: static PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM dm, PetscInt val, PetscInt cellHeight, DMLabel label)
     45: {
     46:   PetscInt       fStart, fEnd, f;

     50:   DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
     51:   for (f = fStart; f < fEnd; ++f) {
     52:     PetscInt supportSize;

     54:     DMPlexGetSupportSize(dm, f, &supportSize);
     55:     if (supportSize == 1) {
     56:       if (val < 0) {
     57:         PetscInt *closure = NULL;
     58:         PetscInt  clSize, cl, cval;

     60:         DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure);
     61:         for (cl = 0; cl < clSize*2; cl += 2) {
     62:           DMLabelGetValue(label, closure[cl], &cval);
     63:           if (cval < 0) continue;
     64:           DMLabelSetValue(label, f, cval);
     65:           break;
     66:         }
     67:         if (cl == clSize*2) {DMLabelSetValue(label, f, 1);}
     68:         DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure);
     69:       } else {
     70:         DMLabelSetValue(label, f, val);
     71:       }
     72:     }
     73:   }
     74:   return(0);
     75: }

     77: /*@
     78:   DMPlexMarkBoundaryFaces - Mark all faces on the boundary

     80:   Not Collective

     82:   Input Parameter:
     83: + dm - The original DM
     84: - val - The marker value, or PETSC_DETERMINE to use some value in the closure (or 1 if none are found)

     86:   Output Parameter:
     87: . label - The DMLabel marking boundary faces with the given value

     89:   Level: developer

     91: .seealso: DMLabelCreate(), DMCreateLabel()
     92: @*/
     93: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, PetscInt val, DMLabel label)
     94: {
     95:   DMPlexInterpolatedFlag  flg;
     96:   PetscErrorCode          ierr;

    100:   DMPlexIsInterpolated(dm, &flg);
    101:   if (flg != DMPLEX_INTERPOLATED_FULL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not fully interpolated on this rank");
    102:   DMPlexMarkBoundaryFaces_Internal(dm, val, 0, label);
    103:   return(0);
    104: }

    106: static PetscErrorCode DMPlexLabelComplete_Internal(DM dm, DMLabel label, PetscBool completeCells)
    107: {
    108:   IS              valueIS;
    109:   PetscSF         sfPoint;
    110:   const PetscInt *values;
    111:   PetscInt        numValues, v, cStart, cEnd, nroots;
    112:   PetscErrorCode  ierr;

    115:   DMLabelGetNumValues(label, &numValues);
    116:   DMLabelGetValueIS(label, &valueIS);
    117:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
    118:   ISGetIndices(valueIS, &values);
    119:   for (v = 0; v < numValues; ++v) {
    120:     IS              pointIS;
    121:     const PetscInt *points;
    122:     PetscInt        numPoints, p;

    124:     DMLabelGetStratumSize(label, values[v], &numPoints);
    125:     DMLabelGetStratumIS(label, values[v], &pointIS);
    126:     ISGetIndices(pointIS, &points);
    127:     for (p = 0; p < numPoints; ++p) {
    128:       PetscInt  q = points[p];
    129:       PetscInt *closure = NULL;
    130:       PetscInt  closureSize, c;

    132:       if (cStart <= q && q < cEnd && !completeCells) { /* skip cells */
    133:         continue;
    134:       }
    135:       DMPlexGetTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
    136:       for (c = 0; c < closureSize*2; c += 2) {
    137:         DMLabelSetValue(label, closure[c], values[v]);
    138:       }
    139:       DMPlexRestoreTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
    140:     }
    141:     ISRestoreIndices(pointIS, &points);
    142:     ISDestroy(&pointIS);
    143:   }
    144:   ISRestoreIndices(valueIS, &values);
    145:   ISDestroy(&valueIS);
    146:   DMGetPointSF(dm, &sfPoint);
    147:   PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
    148:   if (nroots >= 0) {
    149:     DMLabel         lblRoots, lblLeaves;
    150:     IS              valueIS, pointIS;
    151:     const PetscInt *values;
    152:     PetscInt        numValues, v;
    153:     PetscErrorCode  ierr;

    155:     /* Pull point contributions from remote leaves into local roots */
    156:     DMLabelGather(label, sfPoint, &lblLeaves);
    157:     DMLabelGetValueIS(lblLeaves, &valueIS);
    158:     ISGetLocalSize(valueIS, &numValues);
    159:     ISGetIndices(valueIS, &values);
    160:     for (v = 0; v < numValues; ++v) {
    161:       const PetscInt value = values[v];

    163:       DMLabelGetStratumIS(lblLeaves, value, &pointIS);
    164:       DMLabelInsertIS(label, pointIS, value);
    165:       ISDestroy(&pointIS);
    166:     }
    167:     ISRestoreIndices(valueIS, &values);
    168:     ISDestroy(&valueIS);
    169:     DMLabelDestroy(&lblLeaves);
    170:     /* Push point contributions from roots into remote leaves */
    171:     DMLabelDistribute(label, sfPoint, &lblRoots);
    172:     DMLabelGetValueIS(lblRoots, &valueIS);
    173:     ISGetLocalSize(valueIS, &numValues);
    174:     ISGetIndices(valueIS, &values);
    175:     for (v = 0; v < numValues; ++v) {
    176:       const PetscInt value = values[v];

    178:       DMLabelGetStratumIS(lblRoots, value, &pointIS);
    179:       DMLabelInsertIS(label, pointIS, value);
    180:       ISDestroy(&pointIS);
    181:     }
    182:     ISRestoreIndices(valueIS, &values);
    183:     ISDestroy(&valueIS);
    184:     DMLabelDestroy(&lblRoots);
    185:   }
    186:   return(0);
    187: }

    189: /*@
    190:   DMPlexLabelComplete - Starting with a label marking points on a surface, we add the transitive closure to the surface

    192:   Input Parameters:
    193: + dm - The DM
    194: - label - A DMLabel marking the surface points

    196:   Output Parameter:
    197: . label - A DMLabel marking all surface points in the transitive closure

    199:   Level: developer

    201: .seealso: DMPlexLabelCohesiveComplete()
    202: @*/
    203: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
    204: {

    208:   DMPlexLabelComplete_Internal(dm, label, PETSC_TRUE);
    209:   return(0);
    210: }

    212: /*@
    213:   DMPlexLabelAddCells - Starting with a label marking points on a surface, we add a cell for each point

    215:   Input Parameters:
    216: + dm - The DM
    217: - label - A DMLabel marking the surface points

    219:   Output Parameter:
    220: . label - A DMLabel incorporating cells

    222:   Level: developer

    224:   Note: The cells allow FEM boundary conditions to be applied using the cell geometry

    226: .seealso: DMPlexLabelAddFaceCells(), DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
    227: @*/
    228: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
    229: {
    230:   IS              valueIS;
    231:   const PetscInt *values;
    232:   PetscInt        numValues, v, cStart, cEnd;
    233:   PetscErrorCode  ierr;

    236:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
    237:   DMLabelGetNumValues(label, &numValues);
    238:   DMLabelGetValueIS(label, &valueIS);
    239:   ISGetIndices(valueIS, &values);
    240:   for (v = 0; v < numValues; ++v) {
    241:     IS              pointIS;
    242:     const PetscInt *points;
    243:     PetscInt        numPoints, p;

    245:     DMLabelGetStratumSize(label, values[v], &numPoints);
    246:     DMLabelGetStratumIS(label, values[v], &pointIS);
    247:     ISGetIndices(pointIS, &points);
    248:     for (p = 0; p < numPoints; ++p) {
    249:       PetscInt *closure = NULL;
    250:       PetscInt  closureSize, cl;

    252:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
    253:       for (cl = closureSize-1; cl > 0; --cl) {
    254:         const PetscInt cell = closure[cl*2];
    255:         if ((cell >= cStart) && (cell < cEnd)) {DMLabelSetValue(label, cell, values[v]); break;}
    256:       }
    257:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
    258:     }
    259:     ISRestoreIndices(pointIS, &points);
    260:     ISDestroy(&pointIS);
    261:   }
    262:   ISRestoreIndices(valueIS, &values);
    263:   ISDestroy(&valueIS);
    264:   return(0);
    265: }

    267: /*@
    268:   DMPlexLabelAddFaceCells - Starting with a label marking faces on a surface, we add a cell for each face

    270:   Input Parameters:
    271: + dm - The DM
    272: - label - A DMLabel marking the surface points

    274:   Output Parameter:
    275: . label - A DMLabel incorporating cells

    277:   Level: developer

    279:   Note: The cells allow FEM boundary conditions to be applied using the cell geometry

    281: .seealso: DMPlexLabelAddCells(), DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
    282: @*/
    283: PetscErrorCode DMPlexLabelAddFaceCells(DM dm, DMLabel label)
    284: {
    285:   IS              valueIS;
    286:   const PetscInt *values;
    287:   PetscInt        numValues, v, cStart, cEnd, fStart, fEnd;
    288:   PetscErrorCode  ierr;

    291:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
    292:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
    293:   DMLabelGetNumValues(label, &numValues);
    294:   DMLabelGetValueIS(label, &valueIS);
    295:   ISGetIndices(valueIS, &values);
    296:   for (v = 0; v < numValues; ++v) {
    297:     IS              pointIS;
    298:     const PetscInt *points;
    299:     PetscInt        numPoints, p;

    301:     DMLabelGetStratumSize(label, values[v], &numPoints);
    302:     DMLabelGetStratumIS(label, values[v], &pointIS);
    303:     ISGetIndices(pointIS, &points);
    304:     for (p = 0; p < numPoints; ++p) {
    305:       const PetscInt face = points[p];
    306:       PetscInt      *closure = NULL;
    307:       PetscInt       closureSize, cl;

    309:       if ((face < fStart) || (face >= fEnd)) continue;
    310:       DMPlexGetTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure);
    311:       for (cl = closureSize-1; cl > 0; --cl) {
    312:         const PetscInt cell = closure[cl*2];
    313:         if ((cell >= cStart) && (cell < cEnd)) {DMLabelSetValue(label, cell, values[v]); break;}
    314:       }
    315:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure);
    316:     }
    317:     ISRestoreIndices(pointIS, &points);
    318:     ISDestroy(&pointIS);
    319:   }
    320:   ISRestoreIndices(valueIS, &values);
    321:   ISDestroy(&valueIS);
    322:   return(0);
    323: }

    325: /*@
    326:   DMPlexLabelClearCells - Remove cells from a label

    328:   Input Parameters:
    329: + dm - The DM
    330: - label - A DMLabel marking surface points and their adjacent cells

    332:   Output Parameter:
    333: . label - A DMLabel without cells

    335:   Level: developer

    337:   Note: This undoes DMPlexLabelAddCells() or DMPlexLabelAddFaceCells()

    339: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete(), DMPlexLabelAddCells()
    340: @*/
    341: PetscErrorCode DMPlexLabelClearCells(DM dm, DMLabel label)
    342: {
    343:   IS              valueIS;
    344:   const PetscInt *values;
    345:   PetscInt        numValues, v, cStart, cEnd;
    346:   PetscErrorCode  ierr;

    349:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
    350:   DMLabelGetNumValues(label, &numValues);
    351:   DMLabelGetValueIS(label, &valueIS);
    352:   ISGetIndices(valueIS, &values);
    353:   for (v = 0; v < numValues; ++v) {
    354:     IS              pointIS;
    355:     const PetscInt *points;
    356:     PetscInt        numPoints, p;

    358:     DMLabelGetStratumSize(label, values[v], &numPoints);
    359:     DMLabelGetStratumIS(label, values[v], &pointIS);
    360:     ISGetIndices(pointIS, &points);
    361:     for (p = 0; p < numPoints; ++p) {
    362:       PetscInt point = points[p];

    364:       if (point >= cStart && point < cEnd) {
    365:         DMLabelClearValue(label,point,values[v]);
    366:       }
    367:     }
    368:     ISRestoreIndices(pointIS, &points);
    369:     ISDestroy(&pointIS);
    370:   }
    371:   ISRestoreIndices(valueIS, &values);
    372:   ISDestroy(&valueIS);
    373:   return(0);
    374: }

    376: /* take (oldEnd, added) pairs, ordered by height and convert them to (oldstart, newstart) pairs, ordered by ascending
    377:  * index (skipping first, which is (0,0)) */
    378: PETSC_STATIC_INLINE PetscErrorCode DMPlexShiftPointSetUp_Internal(PetscInt depth, PetscInt depthShift[])
    379: {
    380:   PetscInt d, off = 0;

    383:   /* sort by (oldend): yes this is an O(n^2) sort, we expect depth <= 3 */
    384:   for (d = 0; d < depth; d++) {
    385:     PetscInt firstd = d;
    386:     PetscInt firstStart = depthShift[2*d];
    387:     PetscInt e;

    389:     for (e = d+1; e <= depth; e++) {
    390:       if (depthShift[2*e] < firstStart) {
    391:         firstd = e;
    392:         firstStart = depthShift[2*d];
    393:       }
    394:     }
    395:     if (firstd != d) {
    396:       PetscInt swap[2];

    398:       e = firstd;
    399:       swap[0] = depthShift[2*d];
    400:       swap[1] = depthShift[2*d+1];
    401:       depthShift[2*d]   = depthShift[2*e];
    402:       depthShift[2*d+1] = depthShift[2*e+1];
    403:       depthShift[2*e]   = swap[0];
    404:       depthShift[2*e+1] = swap[1];
    405:     }
    406:   }
    407:   /* convert (oldstart, added) to (oldstart, newstart) */
    408:   for (d = 0; d <= depth; d++) {
    409:     off += depthShift[2*d+1];
    410:     depthShift[2*d+1] = depthShift[2*d] + off;
    411:   }
    412:   return(0);
    413: }

    415: /* depthShift is a list of (old, new) pairs */
    416: PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
    417: {
    418:   PetscInt d;
    419:   PetscInt newOff = 0;

    421:   for (d = 0; d <= depth; d++) {
    422:     if (p < depthShift[2*d]) return p + newOff;
    423:     else newOff = depthShift[2*d+1] - depthShift[2*d];
    424:   }
    425:   return p + newOff;
    426: }

    428: /* depthShift is a list of (old, new) pairs */
    429: PETSC_STATIC_INLINE PetscInt DMPlexShiftPointInverse_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
    430: {
    431:   PetscInt d;
    432:   PetscInt newOff = 0;

    434:   for (d = 0; d <= depth; d++) {
    435:     if (p < depthShift[2*d+1]) return p + newOff;
    436:     else newOff = depthShift[2*d] - depthShift[2*d+1];
    437:   }
    438:   return p + newOff;
    439: }

    441: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
    442: {
    443:   PetscInt       depth = 0, d, pStart, pEnd, p;
    444:   DMLabel        depthLabel;

    448:   DMPlexGetDepth(dm, &depth);
    449:   if (depth < 0) return(0);
    450:   /* Step 1: Expand chart */
    451:   DMPlexGetChart(dm, &pStart, &pEnd);
    452:   pEnd = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
    453:   DMPlexSetChart(dmNew, pStart, pEnd);
    454:   DMCreateLabel(dmNew,"depth");
    455:   DMPlexGetDepthLabel(dmNew,&depthLabel);
    456:   DMCreateLabel(dmNew, "celltype");
    457:   /* Step 2: Set cone and support sizes */
    458:   for (d = 0; d <= depth; ++d) {
    459:     PetscInt pStartNew, pEndNew;
    460:     IS pIS;

    462:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
    463:     pStartNew = DMPlexShiftPoint_Internal(pStart, depth, depthShift);
    464:     pEndNew = DMPlexShiftPoint_Internal(pEnd, depth, depthShift);
    465:     ISCreateStride(PETSC_COMM_SELF, pEndNew - pStartNew, pStartNew, 1, &pIS);
    466:     DMLabelSetStratumIS(depthLabel, d, pIS);
    467:     ISDestroy(&pIS);
    468:     for (p = pStart; p < pEnd; ++p) {
    469:       PetscInt       newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
    470:       PetscInt       size;
    471:       DMPolytopeType ct;

    473:       DMPlexGetConeSize(dm, p, &size);
    474:       DMPlexSetConeSize(dmNew, newp, size);
    475:       DMPlexGetSupportSize(dm, p, &size);
    476:       DMPlexSetSupportSize(dmNew, newp, size);
    477:       DMPlexGetCellType(dm, p, &ct);
    478:       DMPlexSetCellType(dmNew, newp, ct);
    479:     }
    480:   }
    481:   return(0);
    482: }

    484: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
    485: {
    486:   PetscInt      *newpoints;
    487:   PetscInt       depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;

    491:   DMPlexGetDepth(dm, &depth);
    492:   if (depth < 0) return(0);
    493:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
    494:   DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
    495:   PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
    496:   /* Step 5: Set cones and supports */
    497:   DMPlexGetChart(dm, &pStart, &pEnd);
    498:   for (p = pStart; p < pEnd; ++p) {
    499:     const PetscInt *points = NULL, *orientations = NULL;
    500:     PetscInt        size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);

    502:     DMPlexGetConeSize(dm, p, &size);
    503:     DMPlexGetCone(dm, p, &points);
    504:     DMPlexGetConeOrientation(dm, p, &orientations);
    505:     for (i = 0; i < size; ++i) {
    506:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
    507:     }
    508:     DMPlexSetCone(dmNew, newp, newpoints);
    509:     DMPlexSetConeOrientation(dmNew, newp, orientations);
    510:     DMPlexGetSupportSize(dm, p, &size);
    511:     DMPlexGetSupportSize(dmNew, newp, &sizeNew);
    512:     DMPlexGetSupport(dm, p, &points);
    513:     for (i = 0; i < size; ++i) {
    514:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
    515:     }
    516:     for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
    517:     DMPlexSetSupport(dmNew, newp, newpoints);
    518:   }
    519:   PetscFree(newpoints);
    520:   return(0);
    521: }

    523: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
    524: {
    525:   PetscSection   coordSection, newCoordSection;
    526:   Vec            coordinates, newCoordinates;
    527:   PetscScalar   *coords, *newCoords;
    528:   PetscInt       coordSize, sStart, sEnd;
    529:   PetscInt       dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v;
    530:   PetscBool      hasCells;

    534:   DMGetCoordinateDim(dm, &dim);
    535:   DMSetCoordinateDim(dmNew, dim);
    536:   DMPlexGetDepth(dm, &depth);
    537:   /* Step 8: Convert coordinates */
    538:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
    539:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
    540:   DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
    541:   DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);
    542:   DMGetCoordinateSection(dm, &coordSection);
    543:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
    544:   PetscSectionSetNumFields(newCoordSection, 1);
    545:   PetscSectionSetFieldComponents(newCoordSection, 0, dim);
    546:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
    547:   hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE;
    548:   PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);
    549:   if (hasCells) {
    550:     for (c = cStart; c < cEnd; ++c) {
    551:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof;

    553:       PetscSectionGetDof(coordSection, c, &dof);
    554:       PetscSectionSetDof(newCoordSection, cNew, dof);
    555:       PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);
    556:     }
    557:   }
    558:   for (v = vStartNew; v < vEndNew; ++v) {
    559:     PetscSectionSetDof(newCoordSection, v, dim);
    560:     PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
    561:   }
    562:   PetscSectionSetUp(newCoordSection);
    563:   DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
    564:   PetscSectionGetStorageSize(newCoordSection, &coordSize);
    565:   VecCreate(PETSC_COMM_SELF, &newCoordinates);
    566:   PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
    567:   VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
    568:   VecSetBlockSize(newCoordinates, dim);
    569:   VecSetType(newCoordinates,VECSTANDARD);
    570:   DMSetCoordinatesLocal(dmNew, newCoordinates);
    571:   DMGetCoordinatesLocal(dm, &coordinates);
    572:   VecGetArray(coordinates, &coords);
    573:   VecGetArray(newCoordinates, &newCoords);
    574:   if (hasCells) {
    575:     for (c = cStart; c < cEnd; ++c) {
    576:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d;

    578:       PetscSectionGetDof(coordSection, c, &dof);
    579:       PetscSectionGetOffset(coordSection, c, &off);
    580:       PetscSectionGetOffset(newCoordSection, cNew, &noff);
    581:       for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
    582:     }
    583:   }
    584:   for (v = vStart; v < vEnd; ++v) {
    585:     PetscInt dof, off, noff, d;

    587:     PetscSectionGetDof(coordSection, v, &dof);
    588:     PetscSectionGetOffset(coordSection, v, &off);
    589:     PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);
    590:     for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
    591:   }
    592:   VecRestoreArray(coordinates, &coords);
    593:   VecRestoreArray(newCoordinates, &newCoords);
    594:   VecDestroy(&newCoordinates);
    595:   PetscSectionDestroy(&newCoordSection);
    596:   return(0);
    597: }

    599: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
    600: {
    601:   PetscInt           depth = 0;
    602:   PetscSF            sfPoint, sfPointNew;
    603:   const PetscSFNode *remotePoints;
    604:   PetscSFNode       *gremotePoints;
    605:   const PetscInt    *localPoints;
    606:   PetscInt          *glocalPoints, *newLocation, *newRemoteLocation;
    607:   PetscInt           numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
    608:   PetscErrorCode     ierr;

    611:   DMPlexGetDepth(dm, &depth);
    612:   /* Step 9: Convert pointSF */
    613:   DMGetPointSF(dm, &sfPoint);
    614:   DMGetPointSF(dmNew, &sfPointNew);
    615:   DMPlexGetChart(dm, &pStart, &pEnd);
    616:   PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
    617:   totShift = DMPlexShiftPoint_Internal(pEnd,depth,depthShift) - pEnd;
    618:   if (numRoots >= 0) {
    619:     PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
    620:     for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
    621:     PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation,MPI_REPLACE);
    622:     PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation,MPI_REPLACE);
    623:     PetscMalloc1(numLeaves,    &glocalPoints);
    624:     PetscMalloc1(numLeaves, &gremotePoints);
    625:     for (l = 0; l < numLeaves; ++l) {
    626:       glocalPoints[l]        = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
    627:       gremotePoints[l].rank  = remotePoints[l].rank;
    628:       gremotePoints[l].index = newRemoteLocation[localPoints[l]];
    629:     }
    630:     PetscFree2(newLocation,newRemoteLocation);
    631:     PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
    632:   }
    633:   return(0);
    634: }

    636: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
    637: {
    638:   PetscSF            sfPoint;
    639:   DMLabel            vtkLabel, ghostLabel;
    640:   const PetscSFNode *leafRemote;
    641:   const PetscInt    *leafLocal;
    642:   PetscInt           depth = 0, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
    643:   PetscMPIInt        rank;
    644:   PetscErrorCode     ierr;

    647:   DMPlexGetDepth(dm, &depth);
    648:   /* Step 10: Convert labels */
    649:   DMGetNumLabels(dm, &numLabels);
    650:   for (l = 0; l < numLabels; ++l) {
    651:     DMLabel         label, newlabel;
    652:     const char     *lname;
    653:     PetscBool       isDepth, isDim;
    654:     IS              valueIS;
    655:     const PetscInt *values;
    656:     PetscInt        numValues, val;

    658:     DMGetLabelName(dm, l, &lname);
    659:     PetscStrcmp(lname, "depth", &isDepth);
    660:     if (isDepth) continue;
    661:     PetscStrcmp(lname, "dim", &isDim);
    662:     if (isDim) continue;
    663:     DMCreateLabel(dmNew, lname);
    664:     DMGetLabel(dm, lname, &label);
    665:     DMGetLabel(dmNew, lname, &newlabel);
    666:     DMLabelGetDefaultValue(label,&val);
    667:     DMLabelSetDefaultValue(newlabel,val);
    668:     DMLabelGetValueIS(label, &valueIS);
    669:     ISGetLocalSize(valueIS, &numValues);
    670:     ISGetIndices(valueIS, &values);
    671:     for (val = 0; val < numValues; ++val) {
    672:       IS              pointIS;
    673:       const PetscInt *points;
    674:       PetscInt        numPoints, p;

    676:       DMLabelGetStratumIS(label, values[val], &pointIS);
    677:       ISGetLocalSize(pointIS, &numPoints);
    678:       ISGetIndices(pointIS, &points);
    679:       for (p = 0; p < numPoints; ++p) {
    680:         const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);

    682:         DMLabelSetValue(newlabel, newpoint, values[val]);
    683:       }
    684:       ISRestoreIndices(pointIS, &points);
    685:       ISDestroy(&pointIS);
    686:     }
    687:     ISRestoreIndices(valueIS, &values);
    688:     ISDestroy(&valueIS);
    689:   }
    690:   /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
    691:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
    692:   DMGetPointSF(dm, &sfPoint);
    693:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
    694:   PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
    695:   DMCreateLabel(dmNew, "vtk");
    696:   DMCreateLabel(dmNew, "ghost");
    697:   DMGetLabel(dmNew, "vtk", &vtkLabel);
    698:   DMGetLabel(dmNew, "ghost", &ghostLabel);
    699:   for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
    700:     for (; c < leafLocal[l] && c < cEnd; ++c) {
    701:       DMLabelSetValue(vtkLabel, c, 1);
    702:     }
    703:     if (leafLocal[l] >= cEnd) break;
    704:     if (leafRemote[l].rank == rank) {
    705:       DMLabelSetValue(vtkLabel, c, 1);
    706:     } else {
    707:       DMLabelSetValue(ghostLabel, c, 2);
    708:     }
    709:   }
    710:   for (; c < cEnd; ++c) {
    711:     DMLabelSetValue(vtkLabel, c, 1);
    712:   }
    713:   DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
    714:   for (f = fStart; f < fEnd; ++f) {
    715:     PetscInt numCells;

    717:     DMPlexGetSupportSize(dmNew, f, &numCells);
    718:     if (numCells < 2) {
    719:       DMLabelSetValue(ghostLabel, f, 1);
    720:     } else {
    721:       const PetscInt *cells = NULL;
    722:       PetscInt        vA, vB;

    724:       DMPlexGetSupport(dmNew, f, &cells);
    725:       DMLabelGetValue(vtkLabel, cells[0], &vA);
    726:       DMLabelGetValue(vtkLabel, cells[1], &vB);
    727:       if (vA != 1 && vB != 1) {DMLabelSetValue(ghostLabel, f, 1);}
    728:     }
    729:   }
    730:   return(0);
    731: }

    733: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
    734: {
    735:   DM             refTree;
    736:   PetscSection   pSec;
    737:   PetscInt       *parents, *childIDs;

    741:   DMPlexGetReferenceTree(dm,&refTree);
    742:   DMPlexSetReferenceTree(dmNew,refTree);
    743:   DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);
    744:   if (pSec) {
    745:     PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
    746:     PetscInt *childIDsShifted;
    747:     PetscSection pSecShifted;

    749:     PetscSectionGetChart(pSec,&pStart,&pEnd);
    750:     DMPlexGetDepth(dm,&depth);
    751:     pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift);
    752:     pEndShifted   = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
    753:     PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);
    754:     PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);
    755:     PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);
    756:     for (p = pStartShifted; p < pEndShifted; p++) {
    757:       /* start off assuming no children */
    758:       PetscSectionSetDof(pSecShifted,p,0);
    759:     }
    760:     for (p = pStart; p < pEnd; p++) {
    761:       PetscInt dof;
    762:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

    764:       PetscSectionGetDof(pSec,p,&dof);
    765:       PetscSectionSetDof(pSecShifted,pNew,dof);
    766:     }
    767:     PetscSectionSetUp(pSecShifted);
    768:     for (p = pStart; p < pEnd; p++) {
    769:       PetscInt dof;
    770:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

    772:       PetscSectionGetDof(pSec,p,&dof);
    773:       if (dof) {
    774:         PetscInt off, offNew;

    776:         PetscSectionGetOffset(pSec,p,&off);
    777:         PetscSectionGetOffset(pSecShifted,pNew,&offNew);
    778:         parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift);
    779:         childIDsShifted[offNew] = childIDs[off];
    780:       }
    781:     }
    782:     DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);
    783:     PetscFree2(parentsShifted,childIDsShifted);
    784:     PetscSectionDestroy(&pSecShifted);
    785:   }
    786:   return(0);
    787: }

    789: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
    790: {
    791:   PetscSF               sf;
    792:   IS                    valueIS;
    793:   const PetscInt       *values, *leaves;
    794:   PetscInt             *depthShift;
    795:   PetscInt              d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
    796:   PetscBool             isper;
    797:   const PetscReal      *maxCell, *L;
    798:   const DMBoundaryType *bd;
    799:   PetscErrorCode        ierr;

    802:   DMGetPointSF(dm, &sf);
    803:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
    804:   nleaves = PetscMax(0, nleaves);
    805:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
    806:   /* Count ghost cells */
    807:   DMLabelGetValueIS(label, &valueIS);
    808:   ISGetLocalSize(valueIS, &numFS);
    809:   ISGetIndices(valueIS, &values);
    810:   Ng   = 0;
    811:   for (fs = 0; fs < numFS; ++fs) {
    812:     IS              faceIS;
    813:     const PetscInt *faces;
    814:     PetscInt        numFaces, f, numBdFaces = 0;

    816:     DMLabelGetStratumIS(label, values[fs], &faceIS);
    817:     ISGetLocalSize(faceIS, &numFaces);
    818:     ISGetIndices(faceIS, &faces);
    819:     for (f = 0; f < numFaces; ++f) {
    820:       PetscInt numChildren;

    822:       PetscFindInt(faces[f], nleaves, leaves, &loc);
    823:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
    824:       /* non-local and ancestors points don't get to register ghosts */
    825:       if (loc >= 0 || numChildren) continue;
    826:       if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
    827:     }
    828:     Ng += numBdFaces;
    829:     ISRestoreIndices(faceIS, &faces);
    830:     ISDestroy(&faceIS);
    831:   }
    832:   DMPlexGetDepth(dm, &depth);
    833:   PetscMalloc1(2*(depth+1), &depthShift);
    834:   for (d = 0; d <= depth; d++) {
    835:     PetscInt dEnd;

    837:     DMPlexGetDepthStratum(dm,d,NULL,&dEnd);
    838:     depthShift[2*d]   = dEnd;
    839:     depthShift[2*d+1] = 0;
    840:   }
    841:   if (depth >= 0) depthShift[2*depth+1] = Ng;
    842:   DMPlexShiftPointSetUp_Internal(depth,depthShift);
    843:   DMPlexShiftSizes_Internal(dm, depthShift, gdm);
    844:   /* Step 3: Set cone/support sizes for new points */
    845:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
    846:   for (c = cEnd; c < cEnd + Ng; ++c) {
    847:     DMPlexSetConeSize(gdm, c, 1);
    848:   }
    849:   for (fs = 0; fs < numFS; ++fs) {
    850:     IS              faceIS;
    851:     const PetscInt *faces;
    852:     PetscInt        numFaces, f;

    854:     DMLabelGetStratumIS(label, values[fs], &faceIS);
    855:     ISGetLocalSize(faceIS, &numFaces);
    856:     ISGetIndices(faceIS, &faces);
    857:     for (f = 0; f < numFaces; ++f) {
    858:       PetscInt size, numChildren;

    860:       PetscFindInt(faces[f], nleaves, leaves, &loc);
    861:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
    862:       if (loc >= 0 || numChildren) continue;
    863:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
    864:       DMPlexGetSupportSize(dm, faces[f], &size);
    865:       if (size != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
    866:       DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
    867:     }
    868:     ISRestoreIndices(faceIS, &faces);
    869:     ISDestroy(&faceIS);
    870:   }
    871:   /* Step 4: Setup ghosted DM */
    872:   DMSetUp(gdm);
    873:   DMPlexShiftPoints_Internal(dm, depthShift, gdm);
    874:   /* Step 6: Set cones and supports for new points */
    875:   ghostCell = cEnd;
    876:   for (fs = 0; fs < numFS; ++fs) {
    877:     IS              faceIS;
    878:     const PetscInt *faces;
    879:     PetscInt        numFaces, f;

    881:     DMLabelGetStratumIS(label, values[fs], &faceIS);
    882:     ISGetLocalSize(faceIS, &numFaces);
    883:     ISGetIndices(faceIS, &faces);
    884:     for (f = 0; f < numFaces; ++f) {
    885:       PetscInt newFace = faces[f] + Ng, numChildren;

    887:       PetscFindInt(faces[f], nleaves, leaves, &loc);
    888:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
    889:       if (loc >= 0 || numChildren) continue;
    890:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
    891:       DMPlexSetCone(gdm, ghostCell, &newFace);
    892:       DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
    893:       ++ghostCell;
    894:     }
    895:     ISRestoreIndices(faceIS, &faces);
    896:     ISDestroy(&faceIS);
    897:   }
    898:   ISRestoreIndices(valueIS, &values);
    899:   ISDestroy(&valueIS);
    900:   DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
    901:   DMPlexShiftSF_Internal(dm, depthShift, gdm);
    902:   DMPlexShiftLabels_Internal(dm, depthShift, gdm);
    903:   DMPlexShiftTree_Internal(dm, depthShift, gdm);
    904:   PetscFree(depthShift);
    905:   for (c = cEnd; c < cEnd + Ng; ++c) {
    906:     DMPlexSetCellType(gdm, c, DM_POLYTOPE_FV_GHOST);
    907:   }
    908:   /* Step 7: Periodicity */
    909:   DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
    910:   DMSetPeriodicity(gdm, isper, maxCell,  L,  bd);
    911:   if (numGhostCells) *numGhostCells = Ng;
    912:   return(0);
    913: }

    915: /*@C
    916:   DMPlexConstructGhostCells - Construct ghost cells which connect to every boundary face

    918:   Collective on dm

    920:   Input Parameters:
    921: + dm - The original DM
    922: - labelName - The label specifying the boundary faces, or "Face Sets" if this is NULL

    924:   Output Parameters:
    925: + numGhostCells - The number of ghost cells added to the DM
    926: - dmGhosted - The new DM

    928:   Note: If no label exists of that name, one will be created marking all boundary faces

    930:   Level: developer

    932: .seealso: DMCreate()
    933: @*/
    934: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
    935: {
    936:   DM             gdm;
    937:   DMLabel        label;
    938:   const char    *name = labelName ? labelName : "Face Sets";
    939:   PetscInt       dim, Ng = 0;
    940:   PetscBool      useCone, useClosure;

    947:   DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
    948:   DMSetType(gdm, DMPLEX);
    949:   DMGetDimension(dm, &dim);
    950:   DMSetDimension(gdm, dim);
    951:   DMGetBasicAdjacency(dm, &useCone, &useClosure);
    952:   DMSetBasicAdjacency(gdm, useCone,  useClosure);
    953:   DMGetLabel(dm, name, &label);
    954:   if (!label) {
    955:     /* Get label for boundary faces */
    956:     DMCreateLabel(dm, name);
    957:     DMGetLabel(dm, name, &label);
    958:     DMPlexMarkBoundaryFaces(dm, 1, label);
    959:   }
    960:   DMPlexConstructGhostCells_Internal(dm, label, &Ng, gdm);
    961:   DMCopyBoundary(dm, gdm);
    962:   DMCopyDisc(dm, gdm);
    963:   gdm->setfromoptionscalled = dm->setfromoptionscalled;
    964:   if (numGhostCells) *numGhostCells = Ng;
    965:   *dmGhosted = gdm;
    966:   return(0);
    967: }

    969: /*
    970:   We are adding three kinds of points here:
    971:     Replicated:     Copies of points which exist in the mesh, such as vertices identified across a fault
    972:     Non-replicated: Points which exist on the fault, but are not replicated
    973:     Hybrid:         Entirely new points, such as cohesive cells

    975:   When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
    976:   each depth so that the new split/hybrid points can be inserted as a block.
    977: */
    978: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DMLabel splitLabel, DM sdm)
    979: {
    980:   MPI_Comm         comm;
    981:   IS               valueIS;
    982:   PetscInt         numSP = 0;       /* The number of depths for which we have replicated points */
    983:   const PetscInt  *values;          /* List of depths for which we have replicated points */
    984:   IS              *splitIS;
    985:   IS              *unsplitIS;
    986:   PetscInt        *numSplitPoints;     /* The number of replicated points at each depth */
    987:   PetscInt        *numUnsplitPoints;   /* The number of non-replicated points at each depth which still give rise to hybrid points */
    988:   PetscInt        *numHybridPoints;    /* The number of new hybrid points at each depth */
    989:   PetscInt        *numHybridPointsOld; /* The number of existing hybrid points at each depth */
    990:   const PetscInt **splitPoints;        /* Replicated points for each depth */
    991:   const PetscInt **unsplitPoints;      /* Non-replicated points for each depth */
    992:   PetscSection     coordSection;
    993:   Vec              coordinates;
    994:   PetscScalar     *coords;
    995:   PetscInt        *depthMax;           /* The first hybrid point at each depth in the original mesh */
    996:   PetscInt        *depthEnd;           /* The point limit at each depth in the original mesh */
    997:   PetscInt        *depthShift;         /* Number of replicated+hybrid points at each depth */
    998:   PetscInt        *pMaxNew;            /* The first replicated point at each depth in the new mesh, hybrids come after this */
    999:   PetscInt        *coneNew, *coneONew, *supportNew;
    1000:   PetscInt         shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
    1001:   PetscErrorCode   ierr;

    1004:   PetscObjectGetComm((PetscObject)dm,&comm);
    1005:   DMGetDimension(dm, &dim);
    1006:   DMPlexGetDepth(dm, &depth);
    1007:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
    1008:   /* We do not want this label automatically computed, instead we compute it here */
    1009:   DMCreateLabel(sdm, "celltype");
    1010:   /* Count split points and add cohesive cells */
    1011:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
    1012:   PetscMalloc5(depth+1,&depthMax,depth+1,&depthEnd,2*(depth+1),&depthShift,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
    1013:   PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
    1014:   for (d = 0; d <= depth; ++d) {
    1015:     DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
    1016:     DMPlexGetTensorPrismBounds_Internal(dm, d, &depthMax[d], NULL);
    1017:     depthEnd[d]           = pMaxNew[d];
    1018:     depthMax[d]           = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
    1019:     numSplitPoints[d]     = 0;
    1020:     numUnsplitPoints[d]   = 0;
    1021:     numHybridPoints[d]    = 0;
    1022:     numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
    1023:     splitPoints[d]        = NULL;
    1024:     unsplitPoints[d]      = NULL;
    1025:     splitIS[d]            = NULL;
    1026:     unsplitIS[d]          = NULL;
    1027:     /* we are shifting the existing hybrid points with the stratum behind them, so
    1028:      * the split comes at the end of the normal points, i.e., at depthMax[d] */
    1029:     depthShift[2*d]       = depthMax[d];
    1030:     depthShift[2*d+1]     = 0;
    1031:   }
    1032:   if (label) {
    1033:     DMLabelGetValueIS(label, &valueIS);
    1034:     ISGetLocalSize(valueIS, &numSP);
    1035:     ISGetIndices(valueIS, &values);
    1036:   }
    1037:   for (sp = 0; sp < numSP; ++sp) {
    1038:     const PetscInt dep = values[sp];

    1040:     if ((dep < 0) || (dep > depth)) continue;
    1041:     DMLabelGetStratumIS(label, dep, &splitIS[dep]);
    1042:     if (splitIS[dep]) {
    1043:       ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
    1044:       ISGetIndices(splitIS[dep], &splitPoints[dep]);
    1045:     }
    1046:     DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
    1047:     if (unsplitIS[dep]) {
    1048:       ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
    1049:       ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
    1050:     }
    1051:   }
    1052:   /* Calculate number of hybrid points */
    1053:   for (d = 1; d <= depth; ++d) numHybridPoints[d]     = numSplitPoints[d-1] + numUnsplitPoints[d-1]; /* There is a hybrid cell/face/edge for every split face/edge/vertex   */
    1054:   for (d = 0; d <= depth; ++d) depthShift[2*d+1]      = numSplitPoints[d] + numHybridPoints[d];
    1055:   DMPlexShiftPointSetUp_Internal(depth,depthShift);
    1056:   /* the end of the points in this stratum that come before the new points:
    1057:    * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
    1058:    * added points */
    1059:   for (d = 0; d <= depth; ++d) pMaxNew[d]             = DMPlexShiftPoint_Internal(pMaxNew[d],depth,depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
    1060:   DMPlexShiftSizes_Internal(dm, depthShift, sdm);
    1061:   /* Step 3: Set cone/support sizes for new points */
    1062:   for (dep = 0; dep <= depth; ++dep) {
    1063:     for (p = 0; p < numSplitPoints[dep]; ++p) {
    1064:       const PetscInt  oldp   = splitPoints[dep][p];
    1065:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
    1066:       const PetscInt  splitp = p    + pMaxNew[dep];
    1067:       const PetscInt *support;
    1068:       PetscInt        coneSize, supportSize, qf, qn, qp, e;

    1070:       DMPlexGetConeSize(dm, oldp, &coneSize);
    1071:       DMPlexSetConeSize(sdm, splitp, coneSize);
    1072:       DMPlexGetSupportSize(dm, oldp, &supportSize);
    1073:       DMPlexSetSupportSize(sdm, splitp, supportSize);
    1074:       if (dep == depth-1) {
    1075:         const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

    1077:         /* Add cohesive cells, they are prisms */
    1078:         DMPlexSetConeSize(sdm, hybcell, 2 + coneSize);
    1079:         switch (coneSize) {
    1080:           case 2: DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_SEG_PRISM_TENSOR);break;
    1081:           case 3: DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_TRI_PRISM_TENSOR);break;
    1082:           case 4: DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_QUAD_PRISM_TENSOR);break;
    1083:         }
    1084:       } else if (dep == 0) {
    1085:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

    1087:         DMPlexGetSupport(dm, oldp, &support);
    1088:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
    1089:           PetscInt val;

    1091:           DMLabelGetValue(label, support[e], &val);
    1092:           if (val == 1) ++qf;
    1093:           if ((val == 1) || (val ==  (shift + 1))) ++qn;
    1094:           if ((val == 1) || (val == -(shift + 1))) ++qp;
    1095:         }
    1096:         /* Split old vertex: Edges into original vertex and new cohesive edge */
    1097:         DMPlexSetSupportSize(sdm, newp, qn+1);
    1098:         /* Split new vertex: Edges into split vertex and new cohesive edge */
    1099:         DMPlexSetSupportSize(sdm, splitp, qp+1);
    1100:         /* Add hybrid edge */
    1101:         DMPlexSetConeSize(sdm, hybedge, 2);
    1102:         DMPlexSetSupportSize(sdm, hybedge, qf);
    1103:         DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR);
    1104:       } else if (dep == dim-2) {
    1105:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

    1107:         DMPlexGetSupport(dm, oldp, &support);
    1108:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
    1109:           PetscInt val;

    1111:           DMLabelGetValue(label, support[e], &val);
    1112:           if (val == dim-1) ++qf;
    1113:           if ((val == dim-1) || (val ==  (shift + dim-1))) ++qn;
    1114:           if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
    1115:         }
    1116:         /* Split old edge: Faces into original edge and cohesive face (positive side?) */
    1117:         DMPlexSetSupportSize(sdm, newp, qn+1);
    1118:         /* Split new edge: Faces into split edge and cohesive face (negative side?) */
    1119:         DMPlexSetSupportSize(sdm, splitp, qp+1);
    1120:         /* Add hybrid face */
    1121:         DMPlexSetConeSize(sdm, hybface, 4);
    1122:         DMPlexSetSupportSize(sdm, hybface, qf);
    1123:         DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR);
    1124:       }
    1125:     }
    1126:   }
    1127:   for (dep = 0; dep <= depth; ++dep) {
    1128:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
    1129:       const PetscInt  oldp   = unsplitPoints[dep][p];
    1130:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
    1131:       const PetscInt *support;
    1132:       PetscInt        coneSize, supportSize, qf, e, s;

    1134:       DMPlexGetConeSize(dm, oldp, &coneSize);
    1135:       DMPlexGetSupportSize(dm, oldp, &supportSize);
    1136:       DMPlexGetSupport(dm, oldp, &support);
    1137:       if (dep == 0) {
    1138:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

    1140:         /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
    1141:         for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
    1142:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
    1143:           if (e >= 0) ++qf;
    1144:         }
    1145:         DMPlexSetSupportSize(sdm, newp, qf+2);
    1146:         /* Add hybrid edge */
    1147:         DMPlexSetConeSize(sdm, hybedge, 2);
    1148:         for (e = 0, qf = 0; e < supportSize; ++e) {
    1149:           PetscInt val;

    1151:           DMLabelGetValue(label, support[e], &val);
    1152:           /* Split and unsplit edges produce hybrid faces */
    1153:           if (val == 1) ++qf;
    1154:           if (val == (shift2 + 1)) ++qf;
    1155:         }
    1156:         DMPlexSetSupportSize(sdm, hybedge, qf);
    1157:         DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR);
    1158:       } else if (dep == dim-2) {
    1159:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
    1160:         PetscInt       val;

    1162:         for (e = 0, qf = 0; e < supportSize; ++e) {
    1163:           DMLabelGetValue(label, support[e], &val);
    1164:           if (val == dim-1) qf += 2;
    1165:           else              ++qf;
    1166:         }
    1167:         /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
    1168:         DMPlexSetSupportSize(sdm, newp, qf+2);
    1169:         /* Add hybrid face */
    1170:         for (e = 0, qf = 0; e < supportSize; ++e) {
    1171:           DMLabelGetValue(label, support[e], &val);
    1172:           if (val == dim-1) ++qf;
    1173:         }
    1174:         DMPlexSetConeSize(sdm, hybface, 4);
    1175:         DMPlexSetSupportSize(sdm, hybface, qf);
    1176:         DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR);
    1177:       }
    1178:     }
    1179:   }
    1180:   /* Step 4: Setup split DM */
    1181:   DMSetUp(sdm);
    1182:   DMPlexShiftPoints_Internal(dm, depthShift, sdm);
    1183:   DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
    1184:   PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
    1185:   /* Step 6: Set cones and supports for new points */
    1186:   for (dep = 0; dep <= depth; ++dep) {
    1187:     for (p = 0; p < numSplitPoints[dep]; ++p) {
    1188:       const PetscInt  oldp   = splitPoints[dep][p];
    1189:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
    1190:       const PetscInt  splitp = p    + pMaxNew[dep];
    1191:       const PetscInt *cone, *support, *ornt;
    1192:       PetscInt        coneSize, supportSize, q, qf, qn, qp, v, e, s;

    1194:       DMPlexGetConeSize(dm, oldp, &coneSize);
    1195:       DMPlexGetCone(dm, oldp, &cone);
    1196:       DMPlexGetConeOrientation(dm, oldp, &ornt);
    1197:       DMPlexGetSupportSize(dm, oldp, &supportSize);
    1198:       DMPlexGetSupport(dm, oldp, &support);
    1199:       if (dep == depth-1) {
    1200:         PetscBool       hasUnsplit = PETSC_FALSE;
    1201:         const PetscInt  hybcell    = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
    1202:         const PetscInt *supportF;

    1204:         /* Split face:       copy in old face to new face to start */
    1205:         DMPlexGetSupport(sdm, newp,  &supportF);
    1206:         DMPlexSetSupport(sdm, splitp, supportF);
    1207:         /* Split old face:   old vertices/edges in cone so no change */
    1208:         /* Split new face:   new vertices/edges in cone */
    1209:         for (q = 0; q < coneSize; ++q) {
    1210:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
    1211:           if (v < 0) {
    1212:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
    1213:             if (v < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", cone[q], dep-1);
    1214:             coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
    1215:             hasUnsplit   = PETSC_TRUE;
    1216:           } else {
    1217:             coneNew[2+q] = v + pMaxNew[dep-1];
    1218:             if (dep > 1) {
    1219:               const PetscInt *econe;
    1220:               PetscInt        econeSize, r, vs, vu;

    1222:               DMPlexGetConeSize(dm, cone[q], &econeSize);
    1223:               DMPlexGetCone(dm, cone[q], &econe);
    1224:               for (r = 0; r < econeSize; ++r) {
    1225:                 PetscFindInt(econe[r], numSplitPoints[dep-2],   splitPoints[dep-2],   &vs);
    1226:                 PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
    1227:                 if (vs >= 0) continue;
    1228:                 if (vu < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", econe[r], dep-2);
    1229:                 hasUnsplit   = PETSC_TRUE;
    1230:               }
    1231:             }
    1232:           }
    1233:         }
    1234:         DMPlexSetCone(sdm, splitp, &coneNew[2]);
    1235:         DMPlexSetConeOrientation(sdm, splitp, ornt);
    1236:         /* Face support */
    1237:         for (s = 0; s < supportSize; ++s) {
    1238:           PetscInt val;

    1240:           DMLabelGetValue(label, support[s], &val);
    1241:           if (val < 0) {
    1242:             /* Split old face:   Replace negative side cell with cohesive cell */
    1243:              DMPlexInsertSupport(sdm, newp, s, hybcell);
    1244:           } else {
    1245:             /* Split new face:   Replace positive side cell with cohesive cell */
    1246:             DMPlexInsertSupport(sdm, splitp, s, hybcell);
    1247:             /* Get orientation for cohesive face */
    1248:             {
    1249:               const PetscInt *ncone, *nconeO;
    1250:               PetscInt        nconeSize, nc;

    1252:               DMPlexGetConeSize(dm, support[s], &nconeSize);
    1253:               DMPlexGetCone(dm, support[s], &ncone);
    1254:               DMPlexGetConeOrientation(dm, support[s], &nconeO);
    1255:               for (nc = 0; nc < nconeSize; ++nc) {
    1256:                 if (ncone[nc] == oldp) {
    1257:                   coneONew[0] = nconeO[nc];
    1258:                   break;
    1259:                 }
    1260:               }
    1261:               if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
    1262:             }
    1263:           }
    1264:         }
    1265:         /* Cohesive cell:    Old and new split face, then new cohesive faces */
    1266:         coneNew[0]  = newp;   /* Extracted negative side orientation above */
    1267:         coneNew[1]  = splitp;
    1268:         coneONew[1] = coneONew[0];
    1269:         for (q = 0; q < coneSize; ++q) {
    1270:           /* Hybrid faces must follow order from oriented end face */
    1271:           const PetscInt o  = coneONew[0];
    1272:           const PetscInt qo = o < 0 ? (-(o+1)+coneSize-q)%coneSize : (q+o)%coneSize;

    1274:           PetscFindInt(cone[qo], numSplitPoints[dep-1], splitPoints[dep-1], &v);
    1275:           if (v < 0) {
    1276:             PetscFindInt(cone[qo], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
    1277:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
    1278:           } else {
    1279:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep];
    1280:           }
    1281:           coneONew[2+q] = ((o < 0) + (ornt[qo] < 0))%2 ? -1 : 0;
    1282:         }
    1283:         DMPlexSetCone(sdm, hybcell, coneNew);
    1284:         DMPlexSetConeOrientation(sdm, hybcell, coneONew);
    1285:         /* Label the hybrid cells on the boundary of the split */
    1286:         if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
    1287:       } else if (dep == 0) {
    1288:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

    1290:         /* Split old vertex: Edges in old split faces and new cohesive edge */
    1291:         for (e = 0, qn = 0; e < supportSize; ++e) {
    1292:           PetscInt val;

    1294:           DMLabelGetValue(label, support[e], &val);
    1295:           if ((val == 1) || (val == (shift + 1))) {
    1296:             supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
    1297:           }
    1298:         }
    1299:         supportNew[qn] = hybedge;
    1300:         DMPlexSetSupport(sdm, newp, supportNew);
    1301:         /* Split new vertex: Edges in new split faces and new cohesive edge */
    1302:         for (e = 0, qp = 0; e < supportSize; ++e) {
    1303:           PetscInt val, edge;

    1305:           DMLabelGetValue(label, support[e], &val);
    1306:           if (val == 1) {
    1307:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
    1308:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
    1309:             supportNew[qp++] = edge + pMaxNew[dep+1];
    1310:           } else if (val == -(shift + 1)) {
    1311:             supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
    1312:           }
    1313:         }
    1314:         supportNew[qp] = hybedge;
    1315:         DMPlexSetSupport(sdm, splitp, supportNew);
    1316:         /* Hybrid edge:    Old and new split vertex */
    1317:         coneNew[0] = newp;
    1318:         coneNew[1] = splitp;
    1319:         DMPlexSetCone(sdm, hybedge, coneNew);
    1320:         for (e = 0, qf = 0; e < supportSize; ++e) {
    1321:           PetscInt val, edge;

    1323:           DMLabelGetValue(label, support[e], &val);
    1324:           if (val == 1) {
    1325:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
    1326:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
    1327:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
    1328:           }
    1329:         }
    1330:         DMPlexSetSupport(sdm, hybedge, supportNew);
    1331:       } else if (dep == dim-2) {
    1332:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

    1334:         /* Split old edge:   old vertices in cone so no change */
    1335:         /* Split new edge:   new vertices in cone */
    1336:         for (q = 0; q < coneSize; ++q) {
    1337:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
    1338:           if (v < 0) {
    1339:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
    1340:             if (v < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", cone[q], dep-1);
    1341:             coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
    1342:           } else {
    1343:             coneNew[q] = v + pMaxNew[dep-1];
    1344:           }
    1345:         }
    1346:         DMPlexSetCone(sdm, splitp, coneNew);
    1347:         /* Split old edge: Faces in positive side cells and old split faces */
    1348:         for (e = 0, q = 0; e < supportSize; ++e) {
    1349:           PetscInt val;

    1351:           DMLabelGetValue(label, support[e], &val);
    1352:           if (val == dim-1) {
    1353:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
    1354:           } else if (val == (shift + dim-1)) {
    1355:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
    1356:           }
    1357:         }
    1358:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
    1359:         DMPlexSetSupport(sdm, newp, supportNew);
    1360:         /* Split new edge: Faces in negative side cells and new split faces */
    1361:         for (e = 0, q = 0; e < supportSize; ++e) {
    1362:           PetscInt val, face;

    1364:           DMLabelGetValue(label, support[e], &val);
    1365:           if (val == dim-1) {
    1366:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
    1367:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
    1368:             supportNew[q++] = face + pMaxNew[dep+1];
    1369:           } else if (val == -(shift + dim-1)) {
    1370:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
    1371:           }
    1372:         }
    1373:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
    1374:         DMPlexSetSupport(sdm, splitp, supportNew);
    1375:         /* Hybrid face */
    1376:         coneNew[0] = newp;
    1377:         coneNew[1] = splitp;
    1378:         for (v = 0; v < coneSize; ++v) {
    1379:           PetscInt vertex;
    1380:           PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
    1381:           if (vertex < 0) {
    1382:             PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
    1383:             if (vertex < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", cone[v], dep-1);
    1384:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
    1385:           } else {
    1386:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
    1387:           }
    1388:         }
    1389:         DMPlexSetCone(sdm, hybface, coneNew);
    1390:         for (e = 0, qf = 0; e < supportSize; ++e) {
    1391:           PetscInt val, face;

    1393:           DMLabelGetValue(label, support[e], &val);
    1394:           if (val == dim-1) {
    1395:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
    1396:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
    1397:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
    1398:           }
    1399:         }
    1400:         DMPlexSetSupport(sdm, hybface, supportNew);
    1401:       }
    1402:     }
    1403:   }
    1404:   for (dep = 0; dep <= depth; ++dep) {
    1405:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
    1406:       const PetscInt  oldp   = unsplitPoints[dep][p];
    1407:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
    1408:       const PetscInt *cone, *support, *ornt;
    1409:       PetscInt        coneSize, supportSize, supportSizeNew, q, qf, e, f, s;

    1411:       DMPlexGetConeSize(dm, oldp, &coneSize);
    1412:       DMPlexGetCone(dm, oldp, &cone);
    1413:       DMPlexGetConeOrientation(dm, oldp, &ornt);
    1414:       DMPlexGetSupportSize(dm, oldp, &supportSize);
    1415:       DMPlexGetSupport(dm, oldp, &support);
    1416:       if (dep == 0) {
    1417:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

    1419:         /* Unsplit vertex */
    1420:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
    1421:         for (s = 0, q = 0; s < supportSize; ++s) {
    1422:           supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
    1423:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
    1424:           if (e >= 0) {
    1425:             supportNew[q++] = e + pMaxNew[dep+1];
    1426:           }
    1427:         }
    1428:         supportNew[q++] = hybedge;
    1429:         supportNew[q++] = hybedge;
    1430:         if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
    1431:         DMPlexSetSupport(sdm, newp, supportNew);
    1432:         /* Hybrid edge */
    1433:         coneNew[0] = newp;
    1434:         coneNew[1] = newp;
    1435:         DMPlexSetCone(sdm, hybedge, coneNew);
    1436:         for (e = 0, qf = 0; e < supportSize; ++e) {
    1437:           PetscInt val, edge;

    1439:           DMLabelGetValue(label, support[e], &val);
    1440:           if (val == 1) {
    1441:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
    1442:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
    1443:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
    1444:           } else if  (val ==  (shift2 + 1)) {
    1445:             PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
    1446:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
    1447:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
    1448:           }
    1449:         }
    1450:         DMPlexSetSupport(sdm, hybedge, supportNew);
    1451:       } else if (dep == dim-2) {
    1452:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

    1454:         /* Unsplit edge: Faces into original edge, split face, and hybrid face twice */
    1455:         for (f = 0, qf = 0; f < supportSize; ++f) {
    1456:           PetscInt val, face;

    1458:           DMLabelGetValue(label, support[f], &val);
    1459:           if (val == dim-1) {
    1460:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
    1461:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
    1462:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
    1463:             supportNew[qf++] = face + pMaxNew[dep+1];
    1464:           } else {
    1465:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
    1466:           }
    1467:         }
    1468:         supportNew[qf++] = hybface;
    1469:         supportNew[qf++] = hybface;
    1470:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
    1471:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
    1472:         DMPlexSetSupport(sdm, newp, supportNew);
    1473:         /* Add hybrid face */
    1474:         coneNew[0] = newp;
    1475:         coneNew[1] = newp;
    1476:         PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
    1477:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
    1478:         coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
    1479:         PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
    1480:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
    1481:         coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
    1482:         DMPlexSetCone(sdm, hybface, coneNew);
    1483:         for (f = 0, qf = 0; f < supportSize; ++f) {
    1484:           PetscInt val, face;

    1486:           DMLabelGetValue(label, support[f], &val);
    1487:           if (val == dim-1) {
    1488:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
    1489:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
    1490:           }
    1491:         }
    1492:         DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
    1493:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
    1494:         DMPlexSetSupport(sdm, hybface, supportNew);
    1495:       }
    1496:     }
    1497:   }
    1498:   /* Step 6b: Replace split points in negative side cones */
    1499:   for (sp = 0; sp < numSP; ++sp) {
    1500:     PetscInt        dep = values[sp];
    1501:     IS              pIS;
    1502:     PetscInt        numPoints;
    1503:     const PetscInt *points;

    1505:     if (dep >= 0) continue;
    1506:     DMLabelGetStratumIS(label, dep, &pIS);
    1507:     if (!pIS) continue;
    1508:     dep  = -dep - shift;
    1509:     ISGetLocalSize(pIS, &numPoints);
    1510:     ISGetIndices(pIS, &points);
    1511:     for (p = 0; p < numPoints; ++p) {
    1512:       const PetscInt  oldp = points[p];
    1513:       const PetscInt  newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
    1514:       const PetscInt *cone;
    1515:       PetscInt        coneSize, c;
    1516:       /* PetscBool       replaced = PETSC_FALSE; */

    1518:       /* Negative edge: replace split vertex */
    1519:       /* Negative cell: replace split face */
    1520:       DMPlexGetConeSize(sdm, newp, &coneSize);
    1521:       DMPlexGetCone(sdm, newp, &cone);
    1522:       for (c = 0; c < coneSize; ++c) {
    1523:         const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c],depth,depthShift);
    1524:         PetscInt       csplitp, cp, val;

    1526:         DMLabelGetValue(label, coldp, &val);
    1527:         if (val == dep-1) {
    1528:           PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
    1529:           if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
    1530:           csplitp  = pMaxNew[dep-1] + cp;
    1531:           DMPlexInsertCone(sdm, newp, c, csplitp);
    1532:           /* replaced = PETSC_TRUE; */
    1533:         }
    1534:       }
    1535:       /* Cells with only a vertex or edge on the submesh have no replacement */
    1536:       /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
    1537:     }
    1538:     ISRestoreIndices(pIS, &points);
    1539:     ISDestroy(&pIS);
    1540:   }
    1541:   /* Step 7: Coordinates */
    1542:   DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
    1543:   DMGetCoordinateSection(sdm, &coordSection);
    1544:   DMGetCoordinatesLocal(sdm, &coordinates);
    1545:   VecGetArray(coordinates, &coords);
    1546:   for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
    1547:     const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
    1548:     const PetscInt splitp = pMaxNew[0] + v;
    1549:     PetscInt       dof, off, soff, d;

    1551:     PetscSectionGetDof(coordSection, newp, &dof);
    1552:     PetscSectionGetOffset(coordSection, newp, &off);
    1553:     PetscSectionGetOffset(coordSection, splitp, &soff);
    1554:     for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
    1555:   }
    1556:   VecRestoreArray(coordinates, &coords);
    1557:   /* Step 8: SF, if I can figure this out we can split the mesh in parallel */
    1558:   DMPlexShiftSF_Internal(dm, depthShift, sdm);
    1559:   /* Step 9: Labels */
    1560:   DMPlexShiftLabels_Internal(dm, depthShift, sdm);
    1561:   DMGetNumLabels(sdm, &numLabels);
    1562:   for (dep = 0; dep <= depth; ++dep) {
    1563:     for (p = 0; p < numSplitPoints[dep]; ++p) {
    1564:       const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
    1565:       const PetscInt splitp = pMaxNew[dep] + p;
    1566:       PetscInt       l;

    1568:       if (splitLabel) {
    1569:         const PetscInt val = 100 + dep;

    1571:         DMLabelSetValue(splitLabel, newp,    val);
    1572:         DMLabelSetValue(splitLabel, splitp, -val);
    1573:       }
    1574:       for (l = 0; l < numLabels; ++l) {
    1575:         DMLabel     mlabel;
    1576:         const char *lname;
    1577:         PetscInt    val;
    1578:         PetscBool   isDepth;

    1580:         DMGetLabelName(sdm, l, &lname);
    1581:         PetscStrcmp(lname, "depth", &isDepth);
    1582:         if (isDepth) continue;
    1583:         DMGetLabel(sdm, lname, &mlabel);
    1584:         DMLabelGetValue(mlabel, newp, &val);
    1585:         if (val >= 0) {
    1586:           DMLabelSetValue(mlabel, splitp, val);
    1587:         }
    1588:       }
    1589:     }
    1590:   }
    1591:   for (sp = 0; sp < numSP; ++sp) {
    1592:     const PetscInt dep = values[sp];

    1594:     if ((dep < 0) || (dep > depth)) continue;
    1595:     if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
    1596:     ISDestroy(&splitIS[dep]);
    1597:     if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
    1598:     ISDestroy(&unsplitIS[dep]);
    1599:   }
    1600:   if (label) {
    1601:     ISRestoreIndices(valueIS, &values);
    1602:     ISDestroy(&valueIS);
    1603:   }
    1604:   for (d = 0; d <= depth; ++d) {
    1605:     DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
    1606:     pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
    1607:   }
    1608:   PetscFree3(coneNew, coneONew, supportNew);
    1609:   PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
    1610:   PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
    1611:   return(0);
    1612: }

    1614: /*@C
    1615:   DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface

    1617:   Collective on dm

    1619:   Input Parameters:
    1620: + dm - The original DM
    1621: - label - The label specifying the boundary faces (this could be auto-generated)

    1623:   Output Parameters:
    1624: + splitLabel - The label containing the split points, or NULL if no output is desired
    1625: - dmSplit - The new DM

    1627:   Level: developer

    1629: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
    1630: @*/
    1631: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DMLabel splitLabel, DM *dmSplit)
    1632: {
    1633:   DM             sdm;
    1634:   PetscInt       dim;

    1640:   DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
    1641:   DMSetType(sdm, DMPLEX);
    1642:   DMGetDimension(dm, &dim);
    1643:   DMSetDimension(sdm, dim);
    1644:   switch (dim) {
    1645:   case 2:
    1646:   case 3:
    1647:     DMPlexConstructCohesiveCells_Internal(dm, label, splitLabel, sdm);
    1648:     break;
    1649:   default:
    1650:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
    1651:   }
    1652:   *dmSplit = sdm;
    1653:   return(0);
    1654: }

    1656: /* Returns the side of the surface for a given cell with a face on the surface */
    1657: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
    1658: {
    1659:   const PetscInt *cone, *ornt;
    1660:   PetscInt        dim, coneSize, c;
    1661:   PetscErrorCode  ierr;

    1664:   *pos = PETSC_TRUE;
    1665:   DMGetDimension(dm, &dim);
    1666:   DMPlexGetConeSize(dm, cell, &coneSize);
    1667:   DMPlexGetCone(dm, cell, &cone);
    1668:   DMPlexGetConeOrientation(dm, cell, &ornt);
    1669:   for (c = 0; c < coneSize; ++c) {
    1670:     if (cone[c] == face) {
    1671:       PetscInt o = ornt[c];

    1673:       if (subdm) {
    1674:         const PetscInt *subcone, *subornt;
    1675:         PetscInt        subpoint, subface, subconeSize, sc;

    1677:         PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
    1678:         PetscFindInt(face, numSubpoints, subpoints, &subface);
    1679:         DMPlexGetConeSize(subdm, subpoint, &subconeSize);
    1680:         DMPlexGetCone(subdm, subpoint, &subcone);
    1681:         DMPlexGetConeOrientation(subdm, subpoint, &subornt);
    1682:         for (sc = 0; sc < subconeSize; ++sc) {
    1683:           if (subcone[sc] == subface) {
    1684:             o = subornt[0];
    1685:             break;
    1686:           }
    1687:         }
    1688:         if (sc >= subconeSize) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find subpoint %d (%d) in cone for subpoint %d (%d)", subface, face, subpoint, cell);
    1689:       }
    1690:       if (o >= 0) *pos = PETSC_TRUE;
    1691:       else        *pos = PETSC_FALSE;
    1692:       break;
    1693:     }
    1694:   }
    1695:   if (c == coneSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cell %d in split face %d support does not have it in the cone", cell, face);
    1696:   return(0);
    1697: }

    1699: /*@
    1700:   DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces
    1701:   to complete the surface

    1703:   Input Parameters:
    1704: + dm     - The DM
    1705: . label  - A DMLabel marking the surface
    1706: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
    1707: . flip   - Flag to flip the submesh normal and replace points on the other side
    1708: - subdm  - The subDM associated with the label, or NULL

    1710:   Output Parameter:
    1711: . label - A DMLabel marking all surface points

    1713:   Note: The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation.

    1715:   Level: developer

    1717: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
    1718: @*/
    1719: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
    1720: {
    1721:   DMLabel         depthLabel;
    1722:   IS              dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL;
    1723:   const PetscInt *points, *subpoints;
    1724:   const PetscInt  rev   = flip ? -1 : 1;
    1725:   PetscInt        shift = 100, shift2 = 200, dim, depth, dep, cStart, cEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
    1726:   PetscErrorCode  ierr;

    1729:   DMPlexGetDepth(dm, &depth);
    1730:   DMGetDimension(dm, &dim);
    1731:   DMPlexGetDepthLabel(dm, &depthLabel);
    1732:   if (subdm) {
    1733:     DMPlexGetSubpointIS(subdm, &subpointIS);
    1734:     if (subpointIS) {
    1735:       ISGetLocalSize(subpointIS, &numSubpoints);
    1736:       ISGetIndices(subpointIS, &subpoints);
    1737:     }
    1738:   }
    1739:   /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
    1740:   DMLabelGetStratumIS(label, dim-1, &dimIS);
    1741:   if (!dimIS) return(0);
    1742:   ISGetLocalSize(dimIS, &numPoints);
    1743:   ISGetIndices(dimIS, &points);
    1744:   for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
    1745:     const PetscInt *support;
    1746:     PetscInt        supportSize, s;

    1748:     DMPlexGetSupportSize(dm, points[p], &supportSize);
    1749: #if 0
    1750:     if (supportSize != 2) {
    1751:       const PetscInt *lp;
    1752:       PetscInt        Nlp, pind;

    1754:       /* Check that for a cell with a single support face, that face is in the SF */
    1755:       /*   THis check only works for the remote side. We would need root side information */
    1756:       PetscSFGetGraph(dm->sf, NULL, &Nlp, &lp, NULL);
    1757:       PetscFindInt(points[p], Nlp, lp, &pind);
    1758:       if (pind < 0) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports, and the face is not shared with another process", points[p], supportSize);
    1759:     }
    1760: #endif
    1761:     DMPlexGetSupport(dm, points[p], &support);
    1762:     for (s = 0; s < supportSize; ++s) {
    1763:       const PetscInt *cone;
    1764:       PetscInt        coneSize, c;
    1765:       PetscBool       pos;

    1767:       GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
    1768:       if (pos) {DMLabelSetValue(label, support[s],  rev*(shift+dim));}
    1769:       else     {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
    1770:       if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
    1771:       /* Put faces touching the fault in the label */
    1772:       DMPlexGetConeSize(dm, support[s], &coneSize);
    1773:       DMPlexGetCone(dm, support[s], &cone);
    1774:       for (c = 0; c < coneSize; ++c) {
    1775:         const PetscInt point = cone[c];

    1777:         DMLabelGetValue(label, point, &val);
    1778:         if (val == -1) {
    1779:           PetscInt *closure = NULL;
    1780:           PetscInt  closureSize, cl;

    1782:           DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
    1783:           for (cl = 0; cl < closureSize*2; cl += 2) {
    1784:             const PetscInt clp  = closure[cl];
    1785:             PetscInt       bval = -1;

    1787:             DMLabelGetValue(label, clp, &val);
    1788:             if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
    1789:             if ((val >= 0) && (val < dim-1) && (bval < 0)) {
    1790:               DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
    1791:               break;
    1792:             }
    1793:           }
    1794:           DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
    1795:         }
    1796:       }
    1797:     }
    1798:   }
    1799:   ISRestoreIndices(dimIS, &points);
    1800:   ISDestroy(&dimIS);
    1801:   if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
    1802:   /* Mark boundary points as unsplit */
    1803:   if (blabel) {
    1804:     DMLabelGetStratumIS(blabel, 1, &dimIS);
    1805:     ISGetLocalSize(dimIS, &numPoints);
    1806:     ISGetIndices(dimIS, &points);
    1807:     for (p = 0; p < numPoints; ++p) {
    1808:       const PetscInt point = points[p];
    1809:       PetscInt       val, bval;

    1811:       DMLabelGetValue(blabel, point, &bval);
    1812:       if (bval >= 0) {
    1813:         DMLabelGetValue(label, point, &val);
    1814:         if ((val < 0) || (val > dim)) {
    1815:           /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
    1816:           DMLabelClearValue(blabel, point, bval);
    1817:         }
    1818:       }
    1819:     }
    1820:     for (p = 0; p < numPoints; ++p) {
    1821:       const PetscInt point = points[p];
    1822:       PetscInt       val, bval;

    1824:       DMLabelGetValue(blabel, point, &bval);
    1825:       if (bval >= 0) {
    1826:         const PetscInt *cone,    *support;
    1827:         PetscInt        coneSize, supportSize, s, valA, valB, valE;

    1829:         /* Mark as unsplit */
    1830:         DMLabelGetValue(label, point, &val);
    1831:         if ((val < 0) || (val > dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d has label value %d, should be part of the fault", point, val);
    1832:         DMLabelClearValue(label, point, val);
    1833:         DMLabelSetValue(label, point, shift2+val);
    1834:         /* Check for cross-edge
    1835:              A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
    1836:         if (val != 0) continue;
    1837:         DMPlexGetSupport(dm, point, &support);
    1838:         DMPlexGetSupportSize(dm, point, &supportSize);
    1839:         for (s = 0; s < supportSize; ++s) {
    1840:           DMPlexGetCone(dm, support[s], &cone);
    1841:           DMPlexGetConeSize(dm, support[s], &coneSize);
    1842:           if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
    1843:           DMLabelGetValue(blabel, cone[0], &valA);
    1844:           DMLabelGetValue(blabel, cone[1], &valB);
    1845:           DMLabelGetValue(blabel, support[s], &valE);
    1846:           if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
    1847:         }
    1848:       }
    1849:     }
    1850:     ISRestoreIndices(dimIS, &points);
    1851:     ISDestroy(&dimIS);
    1852:   }
    1853:   /* Search for other cells/faces/edges connected to the fault by a vertex */
    1854:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
    1855:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
    1856:   DMLabelGetStratumIS(label, 0, &dimIS);
    1857:   if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
    1858:   if (dimIS && crossEdgeIS) {
    1859:     IS vertIS = dimIS;

    1861:     ISExpand(vertIS, crossEdgeIS, &dimIS);
    1862:     ISDestroy(&crossEdgeIS);
    1863:     ISDestroy(&vertIS);
    1864:   }
    1865:   if (!dimIS) {
    1866:     return(0);
    1867:   }
    1868:   ISGetLocalSize(dimIS, &numPoints);
    1869:   ISGetIndices(dimIS, &points);
    1870:   for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
    1871:     PetscInt *star = NULL;
    1872:     PetscInt  starSize, s;
    1873:     PetscInt  again = 1;  /* 0: Finished 1: Keep iterating after a change 2: No change */

    1875:     /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
    1876:     DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
    1877:     while (again) {
    1878:       if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
    1879:       again = 0;
    1880:       for (s = 0; s < starSize*2; s += 2) {
    1881:         const PetscInt  point = star[s];
    1882:         const PetscInt *cone;
    1883:         PetscInt        coneSize, c;

    1885:         if ((point < cStart) || (point >= cEnd)) continue;
    1886:         DMLabelGetValue(label, point, &val);
    1887:         if (val != -1) continue;
    1888:         again = again == 1 ? 1 : 2;
    1889:         DMPlexGetConeSize(dm, point, &coneSize);
    1890:         DMPlexGetCone(dm, point, &cone);
    1891:         for (c = 0; c < coneSize; ++c) {
    1892:           DMLabelGetValue(label, cone[c], &val);
    1893:           if (val != -1) {
    1894:             const PetscInt *ccone;
    1895:             PetscInt        cconeSize, cc, side;

    1897:             if (PetscAbs(val) < shift) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d on cell %d has an invalid label %d", cone[c], point, val);
    1898:             if (val > 0) side =  1;
    1899:             else         side = -1;
    1900:             DMLabelSetValue(label, point, side*(shift+dim));
    1901:             /* Mark cell faces which touch the fault */
    1902:             DMPlexGetConeSize(dm, point, &cconeSize);
    1903:             DMPlexGetCone(dm, point, &ccone);
    1904:             for (cc = 0; cc < cconeSize; ++cc) {
    1905:               PetscInt *closure = NULL;
    1906:               PetscInt  closureSize, cl;

    1908:               DMLabelGetValue(label, ccone[cc], &val);
    1909:               if (val != -1) continue;
    1910:               DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
    1911:               for (cl = 0; cl < closureSize*2; cl += 2) {
    1912:                 const PetscInt clp = closure[cl];

    1914:                 DMLabelGetValue(label, clp, &val);
    1915:                 if (val == -1) continue;
    1916:                 DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
    1917:                 break;
    1918:               }
    1919:               DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
    1920:             }
    1921:             again = 1;
    1922:             break;
    1923:           }
    1924:         }
    1925:       }
    1926:     }
    1927:     /* Classify the rest by cell membership */
    1928:     for (s = 0; s < starSize*2; s += 2) {
    1929:       const PetscInt point = star[s];

    1931:       DMLabelGetValue(label, point, &val);
    1932:       if (val == -1) {
    1933:         PetscInt      *sstar = NULL;
    1934:         PetscInt       sstarSize, ss;
    1935:         PetscBool      marked = PETSC_FALSE, isHybrid;

    1937:         DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
    1938:         for (ss = 0; ss < sstarSize*2; ss += 2) {
    1939:           const PetscInt spoint = sstar[ss];

    1941:           if ((spoint < cStart) || (spoint >= cEnd)) continue;
    1942:           DMLabelGetValue(label, spoint, &val);
    1943:           if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
    1944:           DMLabelGetValue(depthLabel, point, &dep);
    1945:           if (val > 0) {
    1946:             DMLabelSetValue(label, point,   shift+dep);
    1947:           } else {
    1948:             DMLabelSetValue(label, point, -(shift+dep));
    1949:           }
    1950:           marked = PETSC_TRUE;
    1951:           break;
    1952:         }
    1953:         DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
    1954:         DMPlexCellIsHybrid_Internal(dm, point, &isHybrid);
    1955:         if (!isHybrid && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
    1956:       }
    1957:     }
    1958:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
    1959:   }
    1960:   ISRestoreIndices(dimIS, &points);
    1961:   ISDestroy(&dimIS);
    1962:   /* If any faces touching the fault divide cells on either side, split them */
    1963:   DMLabelGetStratumIS(label,   shift+dim-1,  &facePosIS);
    1964:   DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
    1965:   ISExpand(facePosIS, faceNegIS, &dimIS);
    1966:   ISDestroy(&facePosIS);
    1967:   ISDestroy(&faceNegIS);
    1968:   ISGetLocalSize(dimIS, &numPoints);
    1969:   ISGetIndices(dimIS, &points);
    1970:   for (p = 0; p < numPoints; ++p) {
    1971:     const PetscInt  point = points[p];
    1972:     const PetscInt *support;
    1973:     PetscInt        supportSize, valA, valB;

    1975:     DMPlexGetSupportSize(dm, point, &supportSize);
    1976:     if (supportSize != 2) continue;
    1977:     DMPlexGetSupport(dm, point, &support);
    1978:     DMLabelGetValue(label, support[0], &valA);
    1979:     DMLabelGetValue(label, support[1], &valB);
    1980:     if ((valA == -1) || (valB == -1)) continue;
    1981:     if (valA*valB > 0) continue;
    1982:     /* Split the face */
    1983:     DMLabelGetValue(label, point, &valA);
    1984:     DMLabelClearValue(label, point, valA);
    1985:     DMLabelSetValue(label, point, dim-1);
    1986:     /* Label its closure:
    1987:       unmarked: label as unsplit
    1988:       incident: relabel as split
    1989:       split:    do nothing
    1990:     */
    1991:     {
    1992:       PetscInt *closure = NULL;
    1993:       PetscInt  closureSize, cl;

    1995:       DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
    1996:       for (cl = 0; cl < closureSize*2; cl += 2) {
    1997:         DMLabelGetValue(label, closure[cl], &valA);
    1998:         if (valA == -1) { /* Mark as unsplit */
    1999:           DMLabelGetValue(depthLabel, closure[cl], &dep);
    2000:           DMLabelSetValue(label, closure[cl], shift2+dep);
    2001:         } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
    2002:           DMLabelGetValue(depthLabel, closure[cl], &dep);
    2003:           DMLabelClearValue(label, closure[cl], valA);
    2004:           DMLabelSetValue(label, closure[cl], dep);
    2005:         }
    2006:       }
    2007:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
    2008:     }
    2009:   }
    2010:   ISRestoreIndices(dimIS, &points);
    2011:   ISDestroy(&dimIS);
    2012:   return(0);
    2013: }

    2015: /* Check that no cell have all vertices on the fault */
    2016: PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm)
    2017: {
    2018:   IS              subpointIS;
    2019:   const PetscInt *dmpoints;
    2020:   PetscInt        defaultValue, cStart, cEnd, c, vStart, vEnd;
    2021:   PetscErrorCode  ierr;

    2024:   if (!label) return(0);
    2025:   DMLabelGetDefaultValue(label, &defaultValue);
    2026:   DMPlexGetSubpointIS(subdm, &subpointIS);
    2027:   if (!subpointIS) return(0);
    2028:   DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd);
    2029:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
    2030:   ISGetIndices(subpointIS, &dmpoints);
    2031:   for (c = cStart; c < cEnd; ++c) {
    2032:     PetscBool invalidCell = PETSC_TRUE;
    2033:     PetscInt *closure     = NULL;
    2034:     PetscInt  closureSize, cl;

    2036:     DMPlexGetTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
    2037:     for (cl = 0; cl < closureSize*2; cl += 2) {
    2038:       PetscInt value = 0;

    2040:       if ((closure[cl] < vStart) || (closure[cl] >= vEnd)) continue;
    2041:       DMLabelGetValue(label, closure[cl], &value);
    2042:       if (value == defaultValue) {invalidCell = PETSC_FALSE; break;}
    2043:     }
    2044:     DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
    2045:     if (invalidCell) {
    2046:       ISRestoreIndices(subpointIS, &dmpoints);
    2047:       ISDestroy(&subpointIS);
    2048:       DMDestroy(&subdm);
    2049:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ambiguous submesh. Cell %D has all of its vertices on the submesh.", dmpoints[c]);
    2050:     }
    2051:   }
    2052:   ISRestoreIndices(subpointIS, &dmpoints);
    2053:   return(0);
    2054: }

    2057: /*@
    2058:   DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface

    2060:   Collective on dm

    2062:   Input Parameters:
    2063: + dm - The original DM
    2064: . label - The label specifying the interface vertices
    2065: - bdlabel - The optional label specifying the interface boundary vertices

    2067:   Output Parameters:
    2068: + hybridLabel - The label fully marking the interface, or NULL if no output is desired
    2069: . splitLabel - The label containing the split points, or NULL if no output is desired
    2070: . dmInterface - The new interface DM, or NULL
    2071: - dmHybrid - The new DM with cohesive cells

    2073:   Note: The hybridLabel indicates what parts of the original mesh impinged on the on division surface. For points
    2074:   directly on the division surface, they are labeled with their dimension, so an edge 7 on the division surface would be
    2075:   7 (1) in hybridLabel. For points that impinge from the positive side, they are labeled with 100+dim, so an edge 6 with
    2076:   one vertex 3 on the surface would be 6 (101) and 3 (0) in hybridLabel. If an edge 9 from the negative side of the
    2077:   surface also hits vertex 3, it would be 9 (-101) in hybridLabel.

    2079:   The splitLabel indicates what points in the new hybrid mesh were the result of splitting points in the original
    2080:   mesh. The label value is +=100+dim for each point. For example, if two edges 10 and 14 in the hybrid resulting from
    2081:   splitting an edge in the original mesh, you would have 10 (101) and 14 (-101) in the splitLabel.

    2083:   The dmInterface is a DM built from the original division surface. It has a label which can be retrieved using
    2084:   DMPlexGetSubpointMap() which maps each point back to the point in the surface of the original mesh.

    2086:   Level: developer

    2088: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMPlexGetSubpointMap(), DMCreate()
    2089: @*/
    2090: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid)
    2091: {
    2092:   DM             idm;
    2093:   DMLabel        subpointMap, hlabel, slabel = NULL;
    2094:   PetscInt       dim;

    2104:   DMGetDimension(dm, &dim);
    2105:   DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm);
    2106:   DMPlexCheckValidSubmesh_Private(dm, label, idm);
    2107:   DMPlexOrient(idm);
    2108:   DMPlexGetSubpointMap(idm, &subpointMap);
    2109:   DMLabelDuplicate(subpointMap, &hlabel);
    2110:   DMLabelClearStratum(hlabel, dim);
    2111:   if (splitLabel) {
    2112:     const char *name;
    2113:     char        sname[PETSC_MAX_PATH_LEN];

    2115:     PetscObjectGetName((PetscObject) hlabel, &name);
    2116:     PetscStrncpy(sname, name, PETSC_MAX_PATH_LEN);
    2117:     PetscStrcat(sname, " split");
    2118:     DMLabelCreate(PETSC_COMM_SELF, sname, &slabel);
    2119:   }
    2120:   DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, PETSC_FALSE, idm);
    2121:   if (dmInterface) {*dmInterface = idm;}
    2122:   else             {DMDestroy(&idm);}
    2123:   DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid);
    2124:   if (hybridLabel) *hybridLabel = hlabel;
    2125:   else             {DMLabelDestroy(&hlabel);}
    2126:   if (splitLabel)  *splitLabel  = slabel;
    2127:   return(0);
    2128: }

    2130: /* Here we need the explicit assumption that:

    2132:      For any marked cell, the marked vertices constitute a single face
    2133: */
    2134: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
    2135: {
    2136:   IS               subvertexIS = NULL;
    2137:   const PetscInt  *subvertices;
    2138:   PetscInt        *pStart, *pEnd, pSize;
    2139:   PetscInt         depth, dim, d, numSubVerticesInitial = 0, v;
    2140:   PetscErrorCode   ierr;

    2143:   *numFaces = 0;
    2144:   *nFV      = 0;
    2145:   DMPlexGetDepth(dm, &depth);
    2146:   DMGetDimension(dm, &dim);
    2147:   pSize = PetscMax(depth, dim) + 1;
    2148:   PetscMalloc2(pSize, &pStart, pSize, &pEnd);
    2149:   for (d = 0; d <= depth; ++d) {
    2150:     DMPlexGetSimplexOrBoxCells(dm, depth-d, &pStart[d], &pEnd[d]);
    2151:   }
    2152:   /* Loop over initial vertices and mark all faces in the collective star() */
    2153:   if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
    2154:   if (subvertexIS) {
    2155:     ISGetSize(subvertexIS, &numSubVerticesInitial);
    2156:     ISGetIndices(subvertexIS, &subvertices);
    2157:   }
    2158:   for (v = 0; v < numSubVerticesInitial; ++v) {
    2159:     const PetscInt vertex = subvertices[v];
    2160:     PetscInt      *star   = NULL;
    2161:     PetscInt       starSize, s, numCells = 0, c;

    2163:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
    2164:     for (s = 0; s < starSize*2; s += 2) {
    2165:       const PetscInt point = star[s];
    2166:       if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
    2167:     }
    2168:     for (c = 0; c < numCells; ++c) {
    2169:       const PetscInt cell    = star[c];
    2170:       PetscInt      *closure = NULL;
    2171:       PetscInt       closureSize, cl;
    2172:       PetscInt       cellLoc, numCorners = 0, faceSize = 0;

    2174:       DMLabelGetValue(subpointMap, cell, &cellLoc);
    2175:       if (cellLoc == 2) continue;
    2176:       if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
    2177:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
    2178:       for (cl = 0; cl < closureSize*2; cl += 2) {
    2179:         const PetscInt point = closure[cl];
    2180:         PetscInt       vertexLoc;

    2182:         if ((point >= pStart[0]) && (point < pEnd[0])) {
    2183:           ++numCorners;
    2184:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
    2185:           if (vertexLoc == value) closure[faceSize++] = point;
    2186:         }
    2187:       }
    2188:       if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
    2189:       if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
    2190:       if (faceSize == *nFV) {
    2191:         const PetscInt *cells = NULL;
    2192:         PetscInt        numCells, nc;

    2194:         ++(*numFaces);
    2195:         for (cl = 0; cl < faceSize; ++cl) {
    2196:           DMLabelSetValue(subpointMap, closure[cl], 0);
    2197:         }
    2198:         DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
    2199:         for (nc = 0; nc < numCells; ++nc) {
    2200:           DMLabelSetValue(subpointMap, cells[nc], 2);
    2201:         }
    2202:         DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
    2203:       }
    2204:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
    2205:     }
    2206:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
    2207:   }
    2208:   if (subvertexIS) {
    2209:     ISRestoreIndices(subvertexIS, &subvertices);
    2210:   }
    2211:   ISDestroy(&subvertexIS);
    2212:   PetscFree2(pStart, pEnd);
    2213:   return(0);
    2214: }

    2216: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DMLabel subpointMap, DM subdm)
    2217: {
    2218:   IS               subvertexIS = NULL;
    2219:   const PetscInt  *subvertices;
    2220:   PetscInt        *pStart, *pEnd;
    2221:   PetscInt         dim, d, numSubVerticesInitial = 0, v;
    2222:   PetscErrorCode   ierr;

    2225:   DMGetDimension(dm, &dim);
    2226:   PetscMalloc2(dim+1, &pStart, dim+1, &pEnd);
    2227:   for (d = 0; d <= dim; ++d) {
    2228:     DMPlexGetSimplexOrBoxCells(dm, dim-d, &pStart[d], &pEnd[d]);
    2229:   }
    2230:   /* Loop over initial vertices and mark all faces in the collective star() */
    2231:   if (vertexLabel) {
    2232:     DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
    2233:     if (subvertexIS) {
    2234:       ISGetSize(subvertexIS, &numSubVerticesInitial);
    2235:       ISGetIndices(subvertexIS, &subvertices);
    2236:     }
    2237:   }
    2238:   for (v = 0; v < numSubVerticesInitial; ++v) {
    2239:     const PetscInt vertex = subvertices[v];
    2240:     PetscInt      *star   = NULL;
    2241:     PetscInt       starSize, s, numFaces = 0, f;

    2243:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
    2244:     for (s = 0; s < starSize*2; s += 2) {
    2245:       const PetscInt point = star[s];
    2246:       PetscInt       faceLoc;

    2248:       if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) {
    2249:         if (markedFaces) {
    2250:           DMLabelGetValue(vertexLabel, point, &faceLoc);
    2251:           if (faceLoc < 0) continue;
    2252:         }
    2253:         star[numFaces++] = point;
    2254:       }
    2255:     }
    2256:     for (f = 0; f < numFaces; ++f) {
    2257:       const PetscInt face    = star[f];
    2258:       PetscInt      *closure = NULL;
    2259:       PetscInt       closureSize, c;
    2260:       PetscInt       faceLoc;

    2262:       DMLabelGetValue(subpointMap, face, &faceLoc);
    2263:       if (faceLoc == dim-1) continue;
    2264:       if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
    2265:       DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
    2266:       for (c = 0; c < closureSize*2; c += 2) {
    2267:         const PetscInt point = closure[c];
    2268:         PetscInt       vertexLoc;

    2270:         if ((point >= pStart[0]) && (point < pEnd[0])) {
    2271:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
    2272:           if (vertexLoc != value) break;
    2273:         }
    2274:       }
    2275:       if (c == closureSize*2) {
    2276:         const PetscInt *support;
    2277:         PetscInt        supportSize, s;

    2279:         for (c = 0; c < closureSize*2; c += 2) {
    2280:           const PetscInt point = closure[c];

    2282:           for (d = 0; d < dim; ++d) {
    2283:             if ((point >= pStart[d]) && (point < pEnd[d])) {
    2284:               DMLabelSetValue(subpointMap, point, d);
    2285:               break;
    2286:             }
    2287:           }
    2288:         }
    2289:         DMPlexGetSupportSize(dm, face, &supportSize);
    2290:         DMPlexGetSupport(dm, face, &support);
    2291:         for (s = 0; s < supportSize; ++s) {
    2292:           DMLabelSetValue(subpointMap, support[s], dim);
    2293:         }
    2294:       }
    2295:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
    2296:     }
    2297:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
    2298:   }
    2299:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
    2300:   ISDestroy(&subvertexIS);
    2301:   PetscFree2(pStart, pEnd);
    2302:   return(0);
    2303: }

    2305: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
    2306: {
    2307:   DMLabel         label = NULL;
    2308:   const PetscInt *cone;
    2309:   PetscInt        dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
    2310:   PetscErrorCode  ierr;

    2313:   *numFaces = 0;
    2314:   *nFV = 0;
    2315:   if (labelname) {DMGetLabel(dm, labelname, &label);}
    2316:   *subCells = NULL;
    2317:   DMGetDimension(dm, &dim);
    2318:   DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
    2319:   if (cMax < 0) return(0);
    2320:   if (label) {
    2321:     for (c = cMax; c < cEnd; ++c) {
    2322:       PetscInt val;

    2324:       DMLabelGetValue(label, c, &val);
    2325:       if (val == value) {
    2326:         ++(*numFaces);
    2327:         DMPlexGetConeSize(dm, c, &coneSize);
    2328:       }
    2329:     }
    2330:   } else {
    2331:     *numFaces = cEnd - cMax;
    2332:     DMPlexGetConeSize(dm, cMax, &coneSize);
    2333:   }
    2334:   PetscMalloc1(*numFaces *2, subCells);
    2335:   if (!(*numFaces)) return(0);
    2336:   *nFV = hasLagrange ? coneSize/3 : coneSize/2;
    2337:   for (c = cMax; c < cEnd; ++c) {
    2338:     const PetscInt *cells;
    2339:     PetscInt        numCells;

    2341:     if (label) {
    2342:       PetscInt val;

    2344:       DMLabelGetValue(label, c, &val);
    2345:       if (val != value) continue;
    2346:     }
    2347:     DMPlexGetCone(dm, c, &cone);
    2348:     for (p = 0; p < *nFV; ++p) {
    2349:       DMLabelSetValue(subpointMap, cone[p], 0);
    2350:     }
    2351:     /* Negative face */
    2352:     DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
    2353:     /* Not true in parallel
    2354:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
    2355:     for (p = 0; p < numCells; ++p) {
    2356:       DMLabelSetValue(subpointMap, cells[p], 2);
    2357:       (*subCells)[subc++] = cells[p];
    2358:     }
    2359:     DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
    2360:     /* Positive face is not included */
    2361:   }
    2362:   return(0);
    2363: }

    2365: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
    2366: {
    2367:   PetscInt      *pStart, *pEnd;
    2368:   PetscInt       dim, cMax, cEnd, c, d;

    2372:   DMGetDimension(dm, &dim);
    2373:   DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
    2374:   if (cMax < 0) return(0);
    2375:   PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
    2376:   for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
    2377:   for (c = cMax; c < cEnd; ++c) {
    2378:     const PetscInt *cone;
    2379:     PetscInt       *closure = NULL;
    2380:     PetscInt        fconeSize, coneSize, closureSize, cl, val;

    2382:     if (label) {
    2383:       DMLabelGetValue(label, c, &val);
    2384:       if (val != value) continue;
    2385:     }
    2386:     DMPlexGetConeSize(dm, c, &coneSize);
    2387:     DMPlexGetCone(dm, c, &cone);
    2388:     DMPlexGetConeSize(dm, cone[0], &fconeSize);
    2389:     if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
    2390:     /* Negative face */
    2391:     DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
    2392:     for (cl = 0; cl < closureSize*2; cl += 2) {
    2393:       const PetscInt point = closure[cl];

    2395:       for (d = 0; d <= dim; ++d) {
    2396:         if ((point >= pStart[d]) && (point < pEnd[d])) {
    2397:           DMLabelSetValue(subpointMap, point, d);
    2398:           break;
    2399:         }
    2400:       }
    2401:     }
    2402:     DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
    2403:     /* Cells -- positive face is not included */
    2404:     for (cl = 0; cl < 1; ++cl) {
    2405:       const PetscInt *support;
    2406:       PetscInt        supportSize, s;

    2408:       DMPlexGetSupportSize(dm, cone[cl], &supportSize);
    2409:       /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
    2410:       DMPlexGetSupport(dm, cone[cl], &support);
    2411:       for (s = 0; s < supportSize; ++s) {
    2412:         DMLabelSetValue(subpointMap, support[s], dim);
    2413:       }
    2414:     }
    2415:   }
    2416:   PetscFree2(pStart, pEnd);
    2417:   return(0);
    2418: }

    2420: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
    2421: {
    2422:   MPI_Comm       comm;
    2423:   PetscBool      posOrient = PETSC_FALSE;
    2424:   const PetscInt debug     = 0;
    2425:   PetscInt       cellDim, faceSize, f;

    2429:   PetscObjectGetComm((PetscObject)dm,&comm);
    2430:   DMGetDimension(dm, &cellDim);
    2431:   if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}

    2433:   if (cellDim == 1 && numCorners == 2) {
    2434:     /* Triangle */
    2435:     faceSize  = numCorners-1;
    2436:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
    2437:   } else if (cellDim == 2 && numCorners == 3) {
    2438:     /* Triangle */
    2439:     faceSize  = numCorners-1;
    2440:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
    2441:   } else if (cellDim == 3 && numCorners == 4) {
    2442:     /* Tetrahedron */
    2443:     faceSize  = numCorners-1;
    2444:     posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
    2445:   } else if (cellDim == 1 && numCorners == 3) {
    2446:     /* Quadratic line */
    2447:     faceSize  = 1;
    2448:     posOrient = PETSC_TRUE;
    2449:   } else if (cellDim == 2 && numCorners == 4) {
    2450:     /* Quads */
    2451:     faceSize = 2;
    2452:     if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
    2453:       posOrient = PETSC_TRUE;
    2454:     } else if ((indices[0] == 3) && (indices[1] == 0)) {
    2455:       posOrient = PETSC_TRUE;
    2456:     } else {
    2457:       if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
    2458:         posOrient = PETSC_FALSE;
    2459:       } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
    2460:     }
    2461:   } else if (cellDim == 2 && numCorners == 6) {
    2462:     /* Quadratic triangle (I hate this) */
    2463:     /* Edges are determined by the first 2 vertices (corners of edges) */
    2464:     const PetscInt faceSizeTri = 3;
    2465:     PetscInt       sortedIndices[3], i, iFace;
    2466:     PetscBool      found                    = PETSC_FALSE;
    2467:     PetscInt       faceVerticesTriSorted[9] = {
    2468:       0, 3,  4, /* bottom */
    2469:       1, 4,  5, /* right */
    2470:       2, 3,  5, /* left */
    2471:     };
    2472:     PetscInt       faceVerticesTri[9] = {
    2473:       0, 3,  4, /* bottom */
    2474:       1, 4,  5, /* right */
    2475:       2, 5,  3, /* left */
    2476:     };

    2478:     for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
    2479:     PetscSortInt(faceSizeTri, sortedIndices);
    2480:     for (iFace = 0; iFace < 3; ++iFace) {
    2481:       const PetscInt ii = iFace*faceSizeTri;
    2482:       PetscInt       fVertex, cVertex;

    2484:       if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
    2485:           (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
    2486:         for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
    2487:           for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
    2488:             if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
    2489:               faceVertices[fVertex] = origVertices[cVertex];
    2490:               break;
    2491:             }
    2492:           }
    2493:         }
    2494:         found = PETSC_TRUE;
    2495:         break;
    2496:       }
    2497:     }
    2498:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
    2499:     if (posOriented) *posOriented = PETSC_TRUE;
    2500:     return(0);
    2501:   } else if (cellDim == 2 && numCorners == 9) {
    2502:     /* Quadratic quad (I hate this) */
    2503:     /* Edges are determined by the first 2 vertices (corners of edges) */
    2504:     const PetscInt faceSizeQuad = 3;
    2505:     PetscInt       sortedIndices[3], i, iFace;
    2506:     PetscBool      found                      = PETSC_FALSE;
    2507:     PetscInt       faceVerticesQuadSorted[12] = {
    2508:       0, 1,  4, /* bottom */
    2509:       1, 2,  5, /* right */
    2510:       2, 3,  6, /* top */
    2511:       0, 3,  7, /* left */
    2512:     };
    2513:     PetscInt       faceVerticesQuad[12] = {
    2514:       0, 1,  4, /* bottom */
    2515:       1, 2,  5, /* right */
    2516:       2, 3,  6, /* top */
    2517:       3, 0,  7, /* left */
    2518:     };

    2520:     for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
    2521:     PetscSortInt(faceSizeQuad, sortedIndices);
    2522:     for (iFace = 0; iFace < 4; ++iFace) {
    2523:       const PetscInt ii = iFace*faceSizeQuad;
    2524:       PetscInt       fVertex, cVertex;

    2526:       if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
    2527:           (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
    2528:         for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
    2529:           for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
    2530:             if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
    2531:               faceVertices[fVertex] = origVertices[cVertex];
    2532:               break;
    2533:             }
    2534:           }
    2535:         }
    2536:         found = PETSC_TRUE;
    2537:         break;
    2538:       }
    2539:     }
    2540:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
    2541:     if (posOriented) *posOriented = PETSC_TRUE;
    2542:     return(0);
    2543:   } else if (cellDim == 3 && numCorners == 8) {
    2544:     /* Hexes
    2545:        A hex is two oriented quads with the normal of the first
    2546:        pointing up at the second.

    2548:           7---6
    2549:          /|  /|
    2550:         4---5 |
    2551:         | 1-|-2
    2552:         |/  |/
    2553:         0---3

    2555:         Faces are determined by the first 4 vertices (corners of faces) */
    2556:     const PetscInt faceSizeHex = 4;
    2557:     PetscInt       sortedIndices[4], i, iFace;
    2558:     PetscBool      found                     = PETSC_FALSE;
    2559:     PetscInt       faceVerticesHexSorted[24] = {
    2560:       0, 1, 2, 3,  /* bottom */
    2561:       4, 5, 6, 7,  /* top */
    2562:       0, 3, 4, 5,  /* front */
    2563:       2, 3, 5, 6,  /* right */
    2564:       1, 2, 6, 7,  /* back */
    2565:       0, 1, 4, 7,  /* left */
    2566:     };
    2567:     PetscInt       faceVerticesHex[24] = {
    2568:       1, 2, 3, 0,  /* bottom */
    2569:       4, 5, 6, 7,  /* top */
    2570:       0, 3, 5, 4,  /* front */
    2571:       3, 2, 6, 5,  /* right */
    2572:       2, 1, 7, 6,  /* back */
    2573:       1, 0, 4, 7,  /* left */
    2574:     };

    2576:     for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
    2577:     PetscSortInt(faceSizeHex, sortedIndices);
    2578:     for (iFace = 0; iFace < 6; ++iFace) {
    2579:       const PetscInt ii = iFace*faceSizeHex;
    2580:       PetscInt       fVertex, cVertex;

    2582:       if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
    2583:           (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
    2584:           (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
    2585:           (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
    2586:         for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
    2587:           for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
    2588:             if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
    2589:               faceVertices[fVertex] = origVertices[cVertex];
    2590:               break;
    2591:             }
    2592:           }
    2593:         }
    2594:         found = PETSC_TRUE;
    2595:         break;
    2596:       }
    2597:     }
    2598:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
    2599:     if (posOriented) *posOriented = PETSC_TRUE;
    2600:     return(0);
    2601:   } else if (cellDim == 3 && numCorners == 10) {
    2602:     /* Quadratic tet */
    2603:     /* Faces are determined by the first 3 vertices (corners of faces) */
    2604:     const PetscInt faceSizeTet = 6;
    2605:     PetscInt       sortedIndices[6], i, iFace;
    2606:     PetscBool      found                     = PETSC_FALSE;
    2607:     PetscInt       faceVerticesTetSorted[24] = {
    2608:       0, 1, 2,  6, 7, 8, /* bottom */
    2609:       0, 3, 4,  6, 7, 9,  /* front */
    2610:       1, 4, 5,  7, 8, 9,  /* right */
    2611:       2, 3, 5,  6, 8, 9,  /* left */
    2612:     };
    2613:     PetscInt       faceVerticesTet[24] = {
    2614:       0, 1, 2,  6, 7, 8, /* bottom */
    2615:       0, 4, 3,  6, 7, 9,  /* front */
    2616:       1, 5, 4,  7, 8, 9,  /* right */
    2617:       2, 3, 5,  8, 6, 9,  /* left */
    2618:     };

    2620:     for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
    2621:     PetscSortInt(faceSizeTet, sortedIndices);
    2622:     for (iFace=0; iFace < 4; ++iFace) {
    2623:       const PetscInt ii = iFace*faceSizeTet;
    2624:       PetscInt       fVertex, cVertex;

    2626:       if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
    2627:           (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
    2628:           (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
    2629:           (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
    2630:         for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
    2631:           for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
    2632:             if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
    2633:               faceVertices[fVertex] = origVertices[cVertex];
    2634:               break;
    2635:             }
    2636:           }
    2637:         }
    2638:         found = PETSC_TRUE;
    2639:         break;
    2640:       }
    2641:     }
    2642:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
    2643:     if (posOriented) *posOriented = PETSC_TRUE;
    2644:     return(0);
    2645:   } else if (cellDim == 3 && numCorners == 27) {
    2646:     /* Quadratic hexes (I hate this)
    2647:        A hex is two oriented quads with the normal of the first
    2648:        pointing up at the second.

    2650:          7---6
    2651:         /|  /|
    2652:        4---5 |
    2653:        | 3-|-2
    2654:        |/  |/
    2655:        0---1

    2657:        Faces are determined by the first 4 vertices (corners of faces) */
    2658:     const PetscInt faceSizeQuadHex = 9;
    2659:     PetscInt       sortedIndices[9], i, iFace;
    2660:     PetscBool      found                         = PETSC_FALSE;
    2661:     PetscInt       faceVerticesQuadHexSorted[54] = {
    2662:       0, 1, 2, 3,  8, 9, 10, 11,  24, /* bottom */
    2663:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
    2664:       0, 1, 4, 5,  8, 12, 16, 17,  22, /* front */
    2665:       1, 2, 5, 6,  9, 13, 17, 18,  21, /* right */
    2666:       2, 3, 6, 7,  10, 14, 18, 19,  23, /* back */
    2667:       0, 3, 4, 7,  11, 15, 16, 19,  20, /* left */
    2668:     };
    2669:     PetscInt       faceVerticesQuadHex[54] = {
    2670:       3, 2, 1, 0,  10, 9, 8, 11,  24, /* bottom */
    2671:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
    2672:       0, 1, 5, 4,  8, 17, 12, 16,  22, /* front */
    2673:       1, 2, 6, 5,  9, 18, 13, 17,  21, /* right */
    2674:       2, 3, 7, 6,  10, 19, 14, 18,  23, /* back */
    2675:       3, 0, 4, 7,  11, 16, 15, 19,  20 /* left */
    2676:     };

    2678:     for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
    2679:     PetscSortInt(faceSizeQuadHex, sortedIndices);
    2680:     for (iFace = 0; iFace < 6; ++iFace) {
    2681:       const PetscInt ii = iFace*faceSizeQuadHex;
    2682:       PetscInt       fVertex, cVertex;

    2684:       if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
    2685:           (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
    2686:           (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
    2687:           (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
    2688:         for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
    2689:           for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
    2690:             if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
    2691:               faceVertices[fVertex] = origVertices[cVertex];
    2692:               break;
    2693:             }
    2694:           }
    2695:         }
    2696:         found = PETSC_TRUE;
    2697:         break;
    2698:       }
    2699:     }
    2700:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
    2701:     if (posOriented) *posOriented = PETSC_TRUE;
    2702:     return(0);
    2703:   } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
    2704:   if (!posOrient) {
    2705:     if (debug) {PetscPrintf(comm, "  Reversing initial face orientation\n");}
    2706:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
    2707:   } else {
    2708:     if (debug) {PetscPrintf(comm, "  Keeping initial face orientation\n");}
    2709:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
    2710:   }
    2711:   if (posOriented) *posOriented = posOrient;
    2712:   return(0);
    2713: }

    2715: /*@
    2716:   DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices,
    2717:   in faceVertices. The orientation is such that the face normal points out of the cell

    2719:   Not collective

    2721:   Input Parameters:
    2722: + dm           - The original mesh
    2723: . cell         - The cell mesh point
    2724: . faceSize     - The number of vertices on the face
    2725: . face         - The face vertices
    2726: . numCorners   - The number of vertices on the cell
    2727: . indices      - Local numbering of face vertices in cell cone
    2728: - origVertices - Original face vertices

    2730:   Output Parameter:
    2731: + faceVertices - The face vertices properly oriented
    2732: - posOriented  - PETSC_TRUE if the face was oriented with outward normal

    2734:   Level: developer

    2736: .seealso: DMPlexGetCone()
    2737: @*/
    2738: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
    2739: {
    2740:   const PetscInt *cone = NULL;
    2741:   PetscInt        coneSize, v, f, v2;
    2742:   PetscInt        oppositeVertex = -1;
    2743:   PetscErrorCode  ierr;

    2746:   DMPlexGetConeSize(dm, cell, &coneSize);
    2747:   DMPlexGetCone(dm, cell, &cone);
    2748:   for (v = 0, v2 = 0; v < coneSize; ++v) {
    2749:     PetscBool found = PETSC_FALSE;

    2751:     for (f = 0; f < faceSize; ++f) {
    2752:       if (face[f] == cone[v]) {
    2753:         found = PETSC_TRUE; break;
    2754:       }
    2755:     }
    2756:     if (found) {
    2757:       indices[v2]      = v;
    2758:       origVertices[v2] = cone[v];
    2759:       ++v2;
    2760:     } else {
    2761:       oppositeVertex = v;
    2762:     }
    2763:   }
    2764:   DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
    2765:   return(0);
    2766: }

    2768: /*
    2769:   DMPlexInsertFace_Internal - Puts a face into the mesh

    2771:   Not collective

    2773:   Input Parameters:
    2774:   + dm              - The DMPlex
    2775:   . numFaceVertex   - The number of vertices in the face
    2776:   . faceVertices    - The vertices in the face for dm
    2777:   . subfaceVertices - The vertices in the face for subdm
    2778:   . numCorners      - The number of vertices in the cell
    2779:   . cell            - A cell in dm containing the face
    2780:   . subcell         - A cell in subdm containing the face
    2781:   . firstFace       - First face in the mesh
    2782:   - newFacePoint    - Next face in the mesh

    2784:   Output Parameters:
    2785:   . newFacePoint - Contains next face point number on input, updated on output

    2787:   Level: developer
    2788: */
    2789: static PetscErrorCode DMPlexInsertFace_Internal(DM dm, DM subdm, PetscInt numFaceVertices, const PetscInt faceVertices[], const PetscInt subfaceVertices[], PetscInt numCorners, PetscInt cell, PetscInt subcell, PetscInt firstFace, PetscInt *newFacePoint)
    2790: {
    2791:   MPI_Comm        comm;
    2792:   DM_Plex        *submesh = (DM_Plex*) subdm->data;
    2793:   const PetscInt *faces;
    2794:   PetscInt        numFaces, coneSize;
    2795:   PetscErrorCode  ierr;

    2798:   PetscObjectGetComm((PetscObject)dm,&comm);
    2799:   DMPlexGetConeSize(subdm, subcell, &coneSize);
    2800:   if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
    2801: #if 0
    2802:   /* Cannot use this because support() has not been constructed yet */
    2803:   DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
    2804: #else
    2805:   {
    2806:     PetscInt f;

    2808:     numFaces = 0;
    2809:     DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
    2810:     for (f = firstFace; f < *newFacePoint; ++f) {
    2811:       PetscInt dof, off, d;

    2813:       PetscSectionGetDof(submesh->coneSection, f, &dof);
    2814:       PetscSectionGetOffset(submesh->coneSection, f, &off);
    2815:       /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
    2816:       for (d = 0; d < dof; ++d) {
    2817:         const PetscInt p = submesh->cones[off+d];
    2818:         PetscInt       v;

    2820:         for (v = 0; v < numFaceVertices; ++v) {
    2821:           if (subfaceVertices[v] == p) break;
    2822:         }
    2823:         if (v == numFaceVertices) break;
    2824:       }
    2825:       if (d == dof) {
    2826:         numFaces               = 1;
    2827:         ((PetscInt*) faces)[0] = f;
    2828:       }
    2829:     }
    2830:   }
    2831: #endif
    2832:   if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
    2833:   else if (numFaces == 1) {
    2834:     /* Add the other cell neighbor for this face */
    2835:     DMPlexSetCone(subdm, subcell, faces);
    2836:   } else {
    2837:     PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
    2838:     PetscBool posOriented;

    2840:     DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
    2841:     origVertices        = &orientedVertices[numFaceVertices];
    2842:     indices             = &orientedVertices[numFaceVertices*2];
    2843:     orientedSubVertices = &orientedVertices[numFaceVertices*3];
    2844:     DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
    2845:     /* TODO: I know that routine should return a permutation, not the indices */
    2846:     for (v = 0; v < numFaceVertices; ++v) {
    2847:       const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
    2848:       for (ov = 0; ov < numFaceVertices; ++ov) {
    2849:         if (orientedVertices[ov] == vertex) {
    2850:           orientedSubVertices[ov] = subvertex;
    2851:           break;
    2852:         }
    2853:       }
    2854:       if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
    2855:     }
    2856:     DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
    2857:     DMPlexSetCone(subdm, subcell, newFacePoint);
    2858:     DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
    2859:     ++(*newFacePoint);
    2860:   }
    2861: #if 0
    2862:   DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
    2863: #else
    2864:   DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
    2865: #endif
    2866:   return(0);
    2867: }

    2869: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
    2870: {
    2871:   MPI_Comm        comm;
    2872:   DMLabel         subpointMap;
    2873:   IS              subvertexIS,  subcellIS;
    2874:   const PetscInt *subVertices, *subCells;
    2875:   PetscInt        numSubVertices, firstSubVertex, numSubCells;
    2876:   PetscInt       *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
    2877:   PetscInt        vStart, vEnd, c, f;
    2878:   PetscErrorCode  ierr;

    2881:   PetscObjectGetComm((PetscObject)dm,&comm);
    2882:   /* Create subpointMap which marks the submesh */
    2883:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
    2884:   DMPlexSetSubpointMap(subdm, subpointMap);
    2885:   DMLabelDestroy(&subpointMap);
    2886:   if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
    2887:   /* Setup chart */
    2888:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
    2889:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
    2890:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
    2891:   DMPlexSetVTKCellHeight(subdm, 1);
    2892:   /* Set cone sizes */
    2893:   firstSubVertex = numSubCells;
    2894:   firstSubFace   = numSubCells+numSubVertices;
    2895:   newFacePoint   = firstSubFace;
    2896:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
    2897:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
    2898:   DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
    2899:   if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
    2900:   for (c = 0; c < numSubCells; ++c) {
    2901:     DMPlexSetConeSize(subdm, c, 1);
    2902:   }
    2903:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
    2904:     DMPlexSetConeSize(subdm, f, nFV);
    2905:   }
    2906:   DMSetUp(subdm);
    2907:   /* Create face cones */
    2908:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
    2909:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
    2910:   DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
    2911:   for (c = 0; c < numSubCells; ++c) {
    2912:     const PetscInt cell    = subCells[c];
    2913:     const PetscInt subcell = c;
    2914:     PetscInt      *closure = NULL;
    2915:     PetscInt       closureSize, cl, numCorners = 0, faceSize = 0;

    2917:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
    2918:     for (cl = 0; cl < closureSize*2; cl += 2) {
    2919:       const PetscInt point = closure[cl];
    2920:       PetscInt       subVertex;

    2922:       if ((point >= vStart) && (point < vEnd)) {
    2923:         ++numCorners;
    2924:         PetscFindInt(point, numSubVertices, subVertices, &subVertex);
    2925:         if (subVertex >= 0) {
    2926:           closure[faceSize] = point;
    2927:           subface[faceSize] = firstSubVertex+subVertex;
    2928:           ++faceSize;
    2929:         }
    2930:       }
    2931:     }
    2932:     if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
    2933:     if (faceSize == nFV) {
    2934:       DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
    2935:     }
    2936:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
    2937:   }
    2938:   DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
    2939:   DMPlexSymmetrize(subdm);
    2940:   DMPlexStratify(subdm);
    2941:   /* Build coordinates */
    2942:   {
    2943:     PetscSection coordSection, subCoordSection;
    2944:     Vec          coordinates, subCoordinates;
    2945:     PetscScalar *coords, *subCoords;
    2946:     PetscInt     numComp, coordSize, v;
    2947:     const char  *name;

    2949:     DMGetCoordinateSection(dm, &coordSection);
    2950:     DMGetCoordinatesLocal(dm, &coordinates);
    2951:     DMGetCoordinateSection(subdm, &subCoordSection);
    2952:     PetscSectionSetNumFields(subCoordSection, 1);
    2953:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
    2954:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
    2955:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
    2956:     for (v = 0; v < numSubVertices; ++v) {
    2957:       const PetscInt vertex    = subVertices[v];
    2958:       const PetscInt subvertex = firstSubVertex+v;
    2959:       PetscInt       dof;

    2961:       PetscSectionGetDof(coordSection, vertex, &dof);
    2962:       PetscSectionSetDof(subCoordSection, subvertex, dof);
    2963:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
    2964:     }
    2965:     PetscSectionSetUp(subCoordSection);
    2966:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
    2967:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
    2968:     PetscObjectGetName((PetscObject)coordinates,&name);
    2969:     PetscObjectSetName((PetscObject)subCoordinates,name);
    2970:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
    2971:     VecSetType(subCoordinates,VECSTANDARD);
    2972:     if (coordSize) {
    2973:       VecGetArray(coordinates,    &coords);
    2974:       VecGetArray(subCoordinates, &subCoords);
    2975:       for (v = 0; v < numSubVertices; ++v) {
    2976:         const PetscInt vertex    = subVertices[v];
    2977:         const PetscInt subvertex = firstSubVertex+v;
    2978:         PetscInt       dof, off, sdof, soff, d;

    2980:         PetscSectionGetDof(coordSection, vertex, &dof);
    2981:         PetscSectionGetOffset(coordSection, vertex, &off);
    2982:         PetscSectionGetDof(subCoordSection, subvertex, &sdof);
    2983:         PetscSectionGetOffset(subCoordSection, subvertex, &soff);
    2984:         if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
    2985:         for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
    2986:       }
    2987:       VecRestoreArray(coordinates,    &coords);
    2988:       VecRestoreArray(subCoordinates, &subCoords);
    2989:     }
    2990:     DMSetCoordinatesLocal(subdm, subCoordinates);
    2991:     VecDestroy(&subCoordinates);
    2992:   }
    2993:   /* Cleanup */
    2994:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
    2995:   ISDestroy(&subvertexIS);
    2996:   if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
    2997:   ISDestroy(&subcellIS);
    2998:   return(0);
    2999: }

    3001: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
    3002: {
    3003:   PetscInt       subPoint;

    3006:   PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
    3007:   return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
    3008: }

    3010: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
    3011: {
    3012:   MPI_Comm         comm;
    3013:   DMLabel          subpointMap;
    3014:   IS              *subpointIS;
    3015:   const PetscInt **subpoints;
    3016:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
    3017:   PetscInt         totSubPoints = 0, maxConeSize, dim, p, d, v;
    3018:   PetscMPIInt      rank;
    3019:   PetscErrorCode   ierr;

    3022:   PetscObjectGetComm((PetscObject)dm,&comm);
    3023:   MPI_Comm_rank(comm, &rank);
    3024:   /* Create subpointMap which marks the submesh */
    3025:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
    3026:   DMPlexSetSubpointMap(subdm, subpointMap);
    3027:   if (cellHeight) {
    3028:     if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
    3029:     else            {DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm);}
    3030:   } else {
    3031:     DMLabel         depth;
    3032:     IS              pointIS;
    3033:     const PetscInt *points;
    3034:     PetscInt        numPoints=0;

    3036:     DMPlexGetDepthLabel(dm, &depth);
    3037:     DMLabelGetStratumIS(label, value, &pointIS);
    3038:     if (pointIS) {
    3039:       ISGetIndices(pointIS, &points);
    3040:       ISGetLocalSize(pointIS, &numPoints);
    3041:     }
    3042:     for (p = 0; p < numPoints; ++p) {
    3043:       PetscInt *closure = NULL;
    3044:       PetscInt  closureSize, c, pdim;

    3046:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
    3047:       for (c = 0; c < closureSize*2; c += 2) {
    3048:         DMLabelGetValue(depth, closure[c], &pdim);
    3049:         DMLabelSetValue(subpointMap, closure[c], pdim);
    3050:       }
    3051:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
    3052:     }
    3053:     if (pointIS) {ISRestoreIndices(pointIS, &points);}
    3054:     ISDestroy(&pointIS);
    3055:   }
    3056:   /* Setup chart */
    3057:   DMGetDimension(dm, &dim);
    3058:   PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
    3059:   for (d = 0; d <= dim; ++d) {
    3060:     DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
    3061:     totSubPoints += numSubPoints[d];
    3062:   }
    3063:   DMPlexSetChart(subdm, 0, totSubPoints);
    3064:   DMPlexSetVTKCellHeight(subdm, cellHeight);
    3065:   /* Set cone sizes */
    3066:   firstSubPoint[dim] = 0;
    3067:   firstSubPoint[0]   = firstSubPoint[dim] + numSubPoints[dim];
    3068:   if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0]     + numSubPoints[0];}
    3069:   if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
    3070:   for (d = 0; d <= dim; ++d) {
    3071:     DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
    3072:     if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
    3073:   }
    3074:   /* We do not want this label automatically computed, instead we compute it here */
    3075:   DMCreateLabel(subdm, "celltype");
    3076:   for (d = 0; d <= dim; ++d) {
    3077:     for (p = 0; p < numSubPoints[d]; ++p) {
    3078:       const PetscInt  point    = subpoints[d][p];
    3079:       const PetscInt  subpoint = firstSubPoint[d] + p;
    3080:       const PetscInt *cone;
    3081:       PetscInt        coneSize, coneSizeNew, c, val;
    3082:       DMPolytopeType  ct;

    3084:       DMPlexGetConeSize(dm, point, &coneSize);
    3085:       DMPlexSetConeSize(subdm, subpoint, coneSize);
    3086:       DMPlexGetCellType(dm, point, &ct);
    3087:       DMPlexSetCellType(subdm, subpoint, ct);
    3088:       if (cellHeight && (d == dim)) {
    3089:         DMPlexGetCone(dm, point, &cone);
    3090:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
    3091:           DMLabelGetValue(subpointMap, cone[c], &val);
    3092:           if (val >= 0) coneSizeNew++;
    3093:         }
    3094:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
    3095:         DMPlexSetCellType(subdm, subpoint, DM_POLYTOPE_FV_GHOST);
    3096:       }
    3097:     }
    3098:   }
    3099:   DMLabelDestroy(&subpointMap);
    3100:   DMSetUp(subdm);
    3101:   /* Set cones */
    3102:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
    3103:   PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
    3104:   for (d = 0; d <= dim; ++d) {
    3105:     for (p = 0; p < numSubPoints[d]; ++p) {
    3106:       const PetscInt  point    = subpoints[d][p];
    3107:       const PetscInt  subpoint = firstSubPoint[d] + p;
    3108:       const PetscInt *cone, *ornt;
    3109:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;

    3111:       if (d == dim-1) {
    3112:         const PetscInt *support, *cone, *ornt;
    3113:         PetscInt        supportSize, coneSize, s, subc;

    3115:         DMPlexGetSupport(dm, point, &support);
    3116:         DMPlexGetSupportSize(dm, point, &supportSize);
    3117:         for (s = 0; s < supportSize; ++s) {
    3118:           PetscBool isHybrid;

    3120:           DMPlexCellIsHybrid_Internal(dm, support[s], &isHybrid);
    3121:           if (!isHybrid) continue;
    3122:           PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
    3123:           if (subc >= 0) {
    3124:             const PetscInt ccell = subpoints[d+1][subc];

    3126:             DMPlexGetCone(dm, ccell, &cone);
    3127:             DMPlexGetConeSize(dm, ccell, &coneSize);
    3128:             DMPlexGetConeOrientation(dm, ccell, &ornt);
    3129:             for (c = 0; c < coneSize; ++c) {
    3130:               if (cone[c] == point) {
    3131:                 fornt = ornt[c];
    3132:                 break;
    3133:               }
    3134:             }
    3135:             break;
    3136:           }
    3137:         }
    3138:       }
    3139:       DMPlexGetConeSize(dm, point, &coneSize);
    3140:       DMPlexGetConeSize(subdm, subpoint, &subconeSize);
    3141:       DMPlexGetCone(dm, point, &cone);
    3142:       DMPlexGetConeOrientation(dm, point, &ornt);
    3143:       for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
    3144:         PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
    3145:         if (subc >= 0) {
    3146:           coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
    3147:           orntNew[coneSizeNew] = ornt[c];
    3148:           ++coneSizeNew;
    3149:         }
    3150:       }
    3151:       if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
    3152:       if (fornt < 0) {
    3153:         /* This should be replaced by a call to DMPlexReverseCell() */
    3154: #if 0
    3155:         DMPlexReverseCell(subdm, subpoint);
    3156: #else
    3157:         for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) {
    3158:           PetscInt faceSize, tmp;

    3160:           tmp        = coneNew[c];
    3161:           coneNew[c] = coneNew[coneSizeNew-1-c];
    3162:           coneNew[coneSizeNew-1-c] = tmp;
    3163:           DMPlexGetConeSize(dm, cone[c], &faceSize);
    3164:           tmp        = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c];
    3165:           orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c];
    3166:           orntNew[coneSizeNew-1-c] = tmp;
    3167:         }
    3168:       }
    3169:       DMPlexSetCone(subdm, subpoint, coneNew);
    3170:       DMPlexSetConeOrientation(subdm, subpoint, orntNew);
    3171: #endif
    3172:     }
    3173:   }
    3174:   PetscFree2(coneNew,orntNew);
    3175:   DMPlexSymmetrize(subdm);
    3176:   DMPlexStratify(subdm);
    3177:   /* Build coordinates */
    3178:   {
    3179:     PetscSection coordSection, subCoordSection;
    3180:     Vec          coordinates, subCoordinates;
    3181:     PetscScalar *coords, *subCoords;
    3182:     PetscInt     cdim, numComp, coordSize;
    3183:     const char  *name;

    3185:     DMGetCoordinateDim(dm, &cdim);
    3186:     DMGetCoordinateSection(dm, &coordSection);
    3187:     DMGetCoordinatesLocal(dm, &coordinates);
    3188:     DMGetCoordinateSection(subdm, &subCoordSection);
    3189:     PetscSectionSetNumFields(subCoordSection, 1);
    3190:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
    3191:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
    3192:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
    3193:     for (v = 0; v < numSubPoints[0]; ++v) {
    3194:       const PetscInt vertex    = subpoints[0][v];
    3195:       const PetscInt subvertex = firstSubPoint[0]+v;
    3196:       PetscInt       dof;

    3198:       PetscSectionGetDof(coordSection, vertex, &dof);
    3199:       PetscSectionSetDof(subCoordSection, subvertex, dof);
    3200:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
    3201:     }
    3202:     PetscSectionSetUp(subCoordSection);
    3203:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
    3204:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
    3205:     PetscObjectGetName((PetscObject)coordinates,&name);
    3206:     PetscObjectSetName((PetscObject)subCoordinates,name);
    3207:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
    3208:     VecSetBlockSize(subCoordinates, cdim);
    3209:     VecSetType(subCoordinates,VECSTANDARD);
    3210:     VecGetArray(coordinates,    &coords);
    3211:     VecGetArray(subCoordinates, &subCoords);
    3212:     for (v = 0; v < numSubPoints[0]; ++v) {
    3213:       const PetscInt vertex    = subpoints[0][v];
    3214:       const PetscInt subvertex = firstSubPoint[0]+v;
    3215:       PetscInt dof, off, sdof, soff, d;

    3217:       PetscSectionGetDof(coordSection, vertex, &dof);
    3218:       PetscSectionGetOffset(coordSection, vertex, &off);
    3219:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
    3220:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
    3221:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
    3222:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
    3223:     }
    3224:     VecRestoreArray(coordinates,    &coords);
    3225:     VecRestoreArray(subCoordinates, &subCoords);
    3226:     DMSetCoordinatesLocal(subdm, subCoordinates);
    3227:     VecDestroy(&subCoordinates);
    3228:   }
    3229:   /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
    3230:   {
    3231:     PetscSF            sfPoint, sfPointSub;
    3232:     IS                 subpIS;
    3233:     const PetscSFNode *remotePoints;
    3234:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
    3235:     const PetscInt    *localPoints, *subpoints;
    3236:     PetscInt          *slocalPoints;
    3237:     PetscInt           numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
    3238:     PetscMPIInt        rank;

    3240:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
    3241:     DMGetPointSF(dm, &sfPoint);
    3242:     DMGetPointSF(subdm, &sfPointSub);
    3243:     DMPlexGetChart(dm, &pStart, &pEnd);
    3244:     DMPlexGetChart(subdm, NULL, &numSubroots);
    3245:     DMPlexGetSubpointIS(subdm, &subpIS);
    3246:     if (subpIS) {
    3247:       ISGetIndices(subpIS, &subpoints);
    3248:       ISGetLocalSize(subpIS, &numSubpoints);
    3249:     }
    3250:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
    3251:     if (numRoots >= 0) {
    3252:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
    3253:       for (p = 0; p < pEnd-pStart; ++p) {
    3254:         newLocalPoints[p].rank  = -2;
    3255:         newLocalPoints[p].index = -2;
    3256:       }
    3257:       /* Set subleaves */
    3258:       for (l = 0; l < numLeaves; ++l) {
    3259:         const PetscInt point    = localPoints[l];
    3260:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

    3262:         if (subpoint < 0) continue;
    3263:         newLocalPoints[point-pStart].rank  = rank;
    3264:         newLocalPoints[point-pStart].index = subpoint;
    3265:         ++numSubleaves;
    3266:       }
    3267:       /* Must put in owned subpoints */
    3268:       for (p = pStart; p < pEnd; ++p) {
    3269:         const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);

    3271:         if (subpoint < 0) {
    3272:           newOwners[p-pStart].rank  = -3;
    3273:           newOwners[p-pStart].index = -3;
    3274:         } else {
    3275:           newOwners[p-pStart].rank  = rank;
    3276:           newOwners[p-pStart].index = subpoint;
    3277:         }
    3278:       }
    3279:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
    3280:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
    3281:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);
    3282:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);
    3283:       PetscMalloc1(numSubleaves, &slocalPoints);
    3284:       PetscMalloc1(numSubleaves, &sremotePoints);
    3285:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
    3286:         const PetscInt point    = localPoints[l];
    3287:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

    3289:         if (subpoint < 0) continue;
    3290:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
    3291:         slocalPoints[sl]        = subpoint;
    3292:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
    3293:         sremotePoints[sl].index = newLocalPoints[point].index;
    3294:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
    3295:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
    3296:         ++sl;
    3297:       }
    3298:       if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
    3299:       PetscFree2(newLocalPoints,newOwners);
    3300:       PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
    3301:     }
    3302:     if (subpIS) {
    3303:       ISRestoreIndices(subpIS, &subpoints);
    3304:     }
    3305:   }
    3306:   /* Cleanup */
    3307:   for (d = 0; d <= dim; ++d) {
    3308:     if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
    3309:     ISDestroy(&subpointIS[d]);
    3310:   }
    3311:   PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
    3312:   return(0);
    3313: }

    3315: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm)
    3316: {

    3320:   DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm);
    3321:   return(0);
    3322: }

    3324: /*@
    3325:   DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label

    3327:   Input Parameters:
    3328: + dm           - The original mesh
    3329: . vertexLabel  - The DMLabel marking points contained in the surface
    3330: . value        - The label value to use
    3331: - markedFaces  - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked

    3333:   Output Parameter:
    3334: . subdm - The surface mesh

    3336:   Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().

    3338:   Level: developer

    3340: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
    3341: @*/
    3342: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm)
    3343: {
    3344:   DMPlexInterpolatedFlag interpolated;
    3345:   PetscInt       dim, cdim;

    3351:   DMGetDimension(dm, &dim);
    3352:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
    3353:   DMSetType(*subdm, DMPLEX);
    3354:   DMSetDimension(*subdm, dim-1);
    3355:   DMGetCoordinateDim(dm, &cdim);
    3356:   DMSetCoordinateDim(*subdm, cdim);
    3357:   DMPlexIsInterpolated(dm, &interpolated);
    3358:   if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
    3359:   if (interpolated) {
    3360:     DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm);
    3361:   } else {
    3362:     DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
    3363:   }
    3364:   return(0);
    3365: }

    3367: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
    3368: {
    3369:   MPI_Comm        comm;
    3370:   DMLabel         subpointMap;
    3371:   IS              subvertexIS;
    3372:   const PetscInt *subVertices;
    3373:   PetscInt        numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
    3374:   PetscInt       *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
    3375:   PetscInt        c, f;
    3376:   PetscErrorCode  ierr;

    3379:   PetscObjectGetComm((PetscObject)dm, &comm);
    3380:   /* Create subpointMap which marks the submesh */
    3381:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
    3382:   DMPlexSetSubpointMap(subdm, subpointMap);
    3383:   DMLabelDestroy(&subpointMap);
    3384:   DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
    3385:   /* Setup chart */
    3386:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
    3387:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
    3388:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
    3389:   DMPlexSetVTKCellHeight(subdm, 1);
    3390:   /* Set cone sizes */
    3391:   firstSubVertex = numSubCells;
    3392:   firstSubFace   = numSubCells+numSubVertices;
    3393:   newFacePoint   = firstSubFace;
    3394:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
    3395:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
    3396:   for (c = 0; c < numSubCells; ++c) {
    3397:     DMPlexSetConeSize(subdm, c, 1);
    3398:   }
    3399:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
    3400:     DMPlexSetConeSize(subdm, f, nFV);
    3401:   }
    3402:   DMSetUp(subdm);
    3403:   /* Create face cones */
    3404:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
    3405:   DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
    3406:   for (c = 0; c < numSubCells; ++c) {
    3407:     const PetscInt  cell    = subCells[c];
    3408:     const PetscInt  subcell = c;
    3409:     const PetscInt *cone, *cells;
    3410:     PetscBool       isHybrid;
    3411:     PetscInt        numCells, subVertex, p, v;

    3413:     DMPlexCellIsHybrid_Internal(dm, cell, &isHybrid);
    3414:     if (!isHybrid) continue;
    3415:     DMPlexGetCone(dm, cell, &cone);
    3416:     for (v = 0; v < nFV; ++v) {
    3417:       PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
    3418:       subface[v] = firstSubVertex+subVertex;
    3419:     }
    3420:     DMPlexSetCone(subdm, newFacePoint, subface);
    3421:     DMPlexSetCone(subdm, subcell, &newFacePoint);
    3422:     DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
    3423:     /* Not true in parallel
    3424:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
    3425:     for (p = 0; p < numCells; ++p) {
    3426:       PetscInt  negsubcell;
    3427:       PetscBool isHybrid;

    3429:       DMPlexCellIsHybrid_Internal(dm, cells[p], &isHybrid);
    3430:       if (isHybrid) continue;
    3431:       /* I know this is a crap search */
    3432:       for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
    3433:         if (subCells[negsubcell] == cells[p]) break;
    3434:       }
    3435:       if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
    3436:       DMPlexSetCone(subdm, negsubcell, &newFacePoint);
    3437:     }
    3438:     DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
    3439:     ++newFacePoint;
    3440:   }
    3441:   DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
    3442:   DMPlexSymmetrize(subdm);
    3443:   DMPlexStratify(subdm);
    3444:   /* Build coordinates */
    3445:   {
    3446:     PetscSection coordSection, subCoordSection;
    3447:     Vec          coordinates, subCoordinates;
    3448:     PetscScalar *coords, *subCoords;
    3449:     PetscInt     cdim, numComp, coordSize, v;
    3450:     const char  *name;

    3452:     DMGetCoordinateDim(dm, &cdim);
    3453:     DMGetCoordinateSection(dm, &coordSection);
    3454:     DMGetCoordinatesLocal(dm, &coordinates);
    3455:     DMGetCoordinateSection(subdm, &subCoordSection);
    3456:     PetscSectionSetNumFields(subCoordSection, 1);
    3457:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
    3458:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
    3459:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
    3460:     for (v = 0; v < numSubVertices; ++v) {
    3461:       const PetscInt vertex    = subVertices[v];
    3462:       const PetscInt subvertex = firstSubVertex+v;
    3463:       PetscInt       dof;

    3465:       PetscSectionGetDof(coordSection, vertex, &dof);
    3466:       PetscSectionSetDof(subCoordSection, subvertex, dof);
    3467:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
    3468:     }
    3469:     PetscSectionSetUp(subCoordSection);
    3470:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
    3471:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
    3472:     PetscObjectGetName((PetscObject)coordinates,&name);
    3473:     PetscObjectSetName((PetscObject)subCoordinates,name);
    3474:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
    3475:     VecSetBlockSize(subCoordinates, cdim);
    3476:     VecSetType(subCoordinates,VECSTANDARD);
    3477:     VecGetArray(coordinates,    &coords);
    3478:     VecGetArray(subCoordinates, &subCoords);
    3479:     for (v = 0; v < numSubVertices; ++v) {
    3480:       const PetscInt vertex    = subVertices[v];
    3481:       const PetscInt subvertex = firstSubVertex+v;
    3482:       PetscInt       dof, off, sdof, soff, d;

    3484:       PetscSectionGetDof(coordSection, vertex, &dof);
    3485:       PetscSectionGetOffset(coordSection, vertex, &off);
    3486:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
    3487:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
    3488:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
    3489:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
    3490:     }
    3491:     VecRestoreArray(coordinates,    &coords);
    3492:     VecRestoreArray(subCoordinates, &subCoords);
    3493:     DMSetCoordinatesLocal(subdm, subCoordinates);
    3494:     VecDestroy(&subCoordinates);
    3495:   }
    3496:   /* Build SF */
    3497:   CHKMEMQ;
    3498:   {
    3499:     PetscSF            sfPoint, sfPointSub;
    3500:     const PetscSFNode *remotePoints;
    3501:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
    3502:     const PetscInt    *localPoints;
    3503:     PetscInt          *slocalPoints;
    3504:     PetscInt           numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
    3505:     PetscMPIInt        rank;

    3507:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
    3508:     DMGetPointSF(dm, &sfPoint);
    3509:     DMGetPointSF(subdm, &sfPointSub);
    3510:     DMPlexGetChart(dm, &pStart, &pEnd);
    3511:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
    3512:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
    3513:     if (numRoots >= 0) {
    3514:       /* Only vertices should be shared */
    3515:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
    3516:       for (p = 0; p < pEnd-pStart; ++p) {
    3517:         newLocalPoints[p].rank  = -2;
    3518:         newLocalPoints[p].index = -2;
    3519:       }
    3520:       /* Set subleaves */
    3521:       for (l = 0; l < numLeaves; ++l) {
    3522:         const PetscInt point    = localPoints[l];
    3523:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

    3525:         if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
    3526:         if (subPoint < 0) continue;
    3527:         newLocalPoints[point-pStart].rank  = rank;
    3528:         newLocalPoints[point-pStart].index = subPoint;
    3529:         ++numSubLeaves;
    3530:       }
    3531:       /* Must put in owned subpoints */
    3532:       for (p = pStart; p < pEnd; ++p) {
    3533:         const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);

    3535:         if (subPoint < 0) {
    3536:           newOwners[p-pStart].rank  = -3;
    3537:           newOwners[p-pStart].index = -3;
    3538:         } else {
    3539:           newOwners[p-pStart].rank  = rank;
    3540:           newOwners[p-pStart].index = subPoint;
    3541:         }
    3542:       }
    3543:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
    3544:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
    3545:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);
    3546:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);
    3547:       PetscMalloc1(numSubLeaves,    &slocalPoints);
    3548:       PetscMalloc1(numSubLeaves, &sremotePoints);
    3549:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
    3550:         const PetscInt point    = localPoints[l];
    3551:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

    3553:         if (subPoint < 0) continue;
    3554:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
    3555:         slocalPoints[sl]        = subPoint;
    3556:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
    3557:         sremotePoints[sl].index = newLocalPoints[point].index;
    3558:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
    3559:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
    3560:         ++sl;
    3561:       }
    3562:       PetscFree2(newLocalPoints,newOwners);
    3563:       if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
    3564:       PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
    3565:     }
    3566:   }
    3567:   CHKMEMQ;
    3568:   /* Cleanup */
    3569:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
    3570:   ISDestroy(&subvertexIS);
    3571:   PetscFree(subCells);
    3572:   return(0);
    3573: }

    3575: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
    3576: {
    3577:   DMLabel        label = NULL;

    3581:   if (labelname) {DMGetLabel(dm, labelname, &label);}
    3582:   DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm);
    3583:   return(0);
    3584: }

    3586: /*@C
    3587:   DMPlexCreateCohesiveSubmesh - Extract from a mesh with cohesive cells the hypersurface defined by one face of the cells. Optionally, a Label an be given to restrict the cells.

    3589:   Input Parameters:
    3590: + dm          - The original mesh
    3591: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
    3592: . label       - A label name, or NULL
    3593: - value  - A label value

    3595:   Output Parameter:
    3596: . subdm - The surface mesh

    3598:   Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().

    3600:   Level: developer

    3602: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
    3603: @*/
    3604: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
    3605: {
    3606:   PetscInt       dim, cdim, depth;

    3612:   DMGetDimension(dm, &dim);
    3613:   DMPlexGetDepth(dm, &depth);
    3614:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
    3615:   DMSetType(*subdm, DMPLEX);
    3616:   DMSetDimension(*subdm, dim-1);
    3617:   DMGetCoordinateDim(dm, &cdim);
    3618:   DMSetCoordinateDim(*subdm, cdim);
    3619:   if (depth == dim) {
    3620:     DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
    3621:   } else {
    3622:     DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
    3623:   }
    3624:   return(0);
    3625: }

    3627: /*@
    3628:   DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh

    3630:   Input Parameters:
    3631: + dm        - The original mesh
    3632: . cellLabel - The DMLabel marking cells contained in the new mesh
    3633: - value     - The label value to use

    3635:   Output Parameter:
    3636: . subdm - The new mesh

    3638:   Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().

    3640:   Level: developer

    3642: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
    3643: @*/
    3644: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
    3645: {
    3646:   PetscInt       dim;

    3652:   DMGetDimension(dm, &dim);
    3653:   DMCreate(PetscObjectComm((PetscObject) dm), subdm);
    3654:   DMSetType(*subdm, DMPLEX);
    3655:   DMSetDimension(*subdm, dim);
    3656:   /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
    3657:   DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm);
    3658:   return(0);
    3659: }

    3661: /*@
    3662:   DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values

    3664:   Input Parameter:
    3665: . dm - The submesh DM

    3667:   Output Parameter:
    3668: . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh

    3670:   Level: developer

    3672: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointIS()
    3673: @*/
    3674: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
    3675: {
    3679:   *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
    3680:   return(0);
    3681: }

    3683: /*@
    3684:   DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values

    3686:   Input Parameters:
    3687: + dm - The submesh DM
    3688: - subpointMap - The DMLabel of all the points from the original mesh in this submesh

    3690:   Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh()

    3692:   Level: developer

    3694: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointIS()
    3695: @*/
    3696: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
    3697: {
    3698:   DM_Plex       *mesh = (DM_Plex *) dm->data;
    3699:   DMLabel        tmp;

    3704:   tmp  = mesh->subpointMap;
    3705:   mesh->subpointMap = subpointMap;
    3706:   PetscObjectReference((PetscObject) mesh->subpointMap);
    3707:   DMLabelDestroy(&tmp);
    3708:   return(0);
    3709: }

    3711: static PetscErrorCode DMPlexCreateSubpointIS_Internal(DM dm, IS *subpointIS)
    3712: {
    3713:   DMLabel        spmap;
    3714:   PetscInt       depth, d;

    3718:   DMPlexGetSubpointMap(dm, &spmap);
    3719:   DMPlexGetDepth(dm, &depth);
    3720:   if (spmap && depth >= 0) {
    3721:     DM_Plex  *mesh = (DM_Plex *) dm->data;
    3722:     PetscInt *points, *depths;
    3723:     PetscInt  pStart, pEnd, p, off;

    3725:     DMPlexGetChart(dm, &pStart, &pEnd);
    3726:     if (pStart) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
    3727:     PetscMalloc1(pEnd, &points);
    3728:     DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);
    3729:     depths[0] = depth;
    3730:     depths[1] = 0;
    3731:     for (d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
    3732:     for (d = 0, off = 0; d <= depth; ++d) {
    3733:       const PetscInt dep = depths[d];
    3734:       PetscInt       depStart, depEnd, n;

    3736:       DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
    3737:       DMLabelGetStratumSize(spmap, dep, &n);
    3738:       if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
    3739:         if (n != depEnd-depStart) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d at depth %d should be %d", n, dep, depEnd-depStart);
    3740:       } else {
    3741:         if (!n) {
    3742:           if (d == 0) {
    3743:             /* Missing cells */
    3744:             for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
    3745:           } else {
    3746:             /* Missing faces */
    3747:             for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
    3748:           }
    3749:         }
    3750:       }
    3751:       if (n) {
    3752:         IS              is;
    3753:         const PetscInt *opoints;

    3755:         DMLabelGetStratumIS(spmap, dep, &is);
    3756:         ISGetIndices(is, &opoints);
    3757:         for (p = 0; p < n; ++p, ++off) points[off] = opoints[p];
    3758:         ISRestoreIndices(is, &opoints);
    3759:         ISDestroy(&is);
    3760:       }
    3761:     }
    3762:     DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);
    3763:     if (off != pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
    3764:     ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
    3765:     PetscObjectStateGet((PetscObject) spmap, &mesh->subpointState);
    3766:   }
    3767:   return(0);
    3768: }

    3770: /*@
    3771:   DMPlexGetSubpointIS - Returns an IS covering the entire subdm chart with the original points as data

    3773:   Input Parameter:
    3774: . dm - The submesh DM

    3776:   Output Parameter:
    3777: . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh

    3779:   Note: This IS is guaranteed to be sorted by the construction of the submesh

    3781:   Level: developer

    3783: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
    3784: @*/
    3785: PetscErrorCode DMPlexGetSubpointIS(DM dm, IS *subpointIS)
    3786: {
    3787:   DM_Plex         *mesh = (DM_Plex *) dm->data;
    3788:   DMLabel          spmap;
    3789:   PetscObjectState state;
    3790:   PetscErrorCode   ierr;

    3795:   DMPlexGetSubpointMap(dm, &spmap);
    3796:   PetscObjectStateGet((PetscObject) spmap, &state);
    3797:   if (state != mesh->subpointState || !mesh->subpointIS) {DMPlexCreateSubpointIS_Internal(dm, &mesh->subpointIS);}
    3798:   *subpointIS = mesh->subpointIS;
    3799:   return(0);
    3800: }

    3802: /*@
    3803:   DMGetEnclosureRelation - Get the relationship between dmA and dmB

    3805:   Input Parameters:
    3806: + dmA - The first DM
    3807: - dmB - The second DM

    3809:   Output Parameter:
    3810: . rel - The relation of dmA to dmB

    3812:   Level: intermediate

    3814: .seealso: DMGetEnclosurePoint()
    3815: @*/
    3816: PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel)
    3817: {
    3818:   DM             plexA, plexB, sdm;
    3819:   DMLabel        spmap;
    3820:   PetscInt       pStartA, pEndA, pStartB, pEndB, NpA, NpB;

    3825:   *rel = DM_ENC_NONE;
    3826:   if (!dmA || !dmB) return(0);
    3829:   if (dmA == dmB) {*rel = DM_ENC_EQUALITY; return(0);}
    3830:   DMConvert(dmA, DMPLEX, &plexA);
    3831:   DMConvert(dmB, DMPLEX, &plexB);
    3832:   DMPlexGetChart(plexA, &pStartA, &pEndA);
    3833:   DMPlexGetChart(plexB, &pStartB, &pEndB);
    3834:   /* Assumption 1: subDMs have smaller charts than the DMs that they originate from
    3835:     - The degenerate case of a subdomain which includes all of the domain on some process can be treated as equality */
    3836:   if ((pStartA == pStartB) && (pEndA == pEndB)) {
    3837:     *rel = DM_ENC_EQUALITY;
    3838:     goto end;
    3839:   }
    3840:   NpA = pEndA - pStartA;
    3841:   NpB = pEndB - pStartB;
    3842:   if (NpA == NpB) goto end;
    3843:   sdm = NpA > NpB ? plexB : plexA; /* The other is the original, enclosing dm */
    3844:   DMPlexGetSubpointMap(sdm, &spmap);
    3845:   if (!spmap) goto end;
    3846:   /* TODO Check the space mapped to by subpointMap is same size as dm */
    3847:   if (NpA > NpB) {
    3848:     *rel = DM_ENC_SUPERMESH;
    3849:   } else {
    3850:     *rel = DM_ENC_SUBMESH;
    3851:   }
    3852:   end:
    3853:   DMDestroy(&plexA);
    3854:   DMDestroy(&plexB);
    3855:   return(0);
    3856: }

    3858: /*@
    3859:   DMGetEnclosurePoint - Get the point pA in dmA which corresponds to the point pB in dmB

    3861:   Input Parameters:
    3862: + dmA   - The first DM
    3863: . dmB   - The second DM
    3864: . etype - The type of enclosure relation that dmA has to dmB
    3865: - pB    - A point of dmB

    3867:   Output Parameter:
    3868: . pA    - The corresponding point of dmA

    3870:   Level: intermediate

    3872: .seealso: DMGetEnclosureRelation()
    3873: @*/
    3874: PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA)
    3875: {
    3876:   DM              sdm;
    3877:   IS              subpointIS;
    3878:   const PetscInt *subpoints;
    3879:   PetscInt        numSubpoints;
    3880:   PetscErrorCode  ierr;

    3883:   /* TODO Cache the IS, making it look like an index */
    3884:   switch (etype) {
    3885:     case DM_ENC_SUPERMESH:
    3886:     sdm  = dmB;
    3887:     DMPlexGetSubpointIS(sdm, &subpointIS);
    3888:     ISGetIndices(subpointIS, &subpoints);
    3889:     *pA  = subpoints[pB];
    3890:     ISRestoreIndices(subpointIS, &subpoints);
    3891:     break;
    3892:     case DM_ENC_SUBMESH:
    3893:     sdm  = dmA;
    3894:     DMPlexGetSubpointIS(sdm, &subpointIS);
    3895:     ISGetLocalSize(subpointIS, &numSubpoints);
    3896:     ISGetIndices(subpointIS, &subpoints);
    3897:     PetscFindInt(pB, numSubpoints, subpoints, pA);
    3898:     if (*pA < 0) {
    3899:       DMViewFromOptions(dmA, NULL, "-dm_enc_A_view");
    3900:       DMViewFromOptions(dmB, NULL, "-dm_enc_B_view");
    3901:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", pB);
    3902:     }
    3903:     ISRestoreIndices(subpointIS, &subpoints);
    3904:     break;
    3905:     case DM_ENC_EQUALITY:
    3906:     case DM_ENC_NONE:
    3907:     *pA = pB;break;
    3908:     case DM_ENC_UNKNOWN:
    3909:     {
    3910:       DMEnclosureType enc;

    3912:       DMGetEnclosureRelation(dmA, dmB, &enc);
    3913:       DMGetEnclosurePoint(dmA, dmB, enc, pB, pA);
    3914:     }
    3915:     break;
    3916:     default: SETERRQ1(PetscObjectComm((PetscObject) dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int) etype);
    3917:   }
    3918:   return(0);
    3919: }