Actual source code: plexinterpolate.c

      1: #include <petsc/private/dmpleximpl.h>
      2: #include <petsc/private/hashmapi.h>
      3: #include <petsc/private/hashmapij.h>

      5: const char *const DMPlexInterpolatedFlags[] = {"none", "partial", "mixed", "full", "DMPlexInterpolatedFlag", "DMPLEX_INTERPOLATED_", NULL};

      7: /* HMapIJKL */

      9: #include <petsc/private/hashmapijkl.h>

     11: static PetscSFNode _PetscInvalidSFNode = {-1, -1};

     13: typedef struct _PetscHMapIJKLRemoteKey {
     14:   PetscSFNode i, j, k, l;
     15: } PetscHMapIJKLRemoteKey;

     17: #define PetscHMapIJKLRemoteKeyHash(key) \
     18:   PetscHashCombine(PetscHashCombine(PetscHashInt((key).i.rank + (key).i.index), PetscHashInt((key).j.rank + (key).j.index)), PetscHashCombine(PetscHashInt((key).k.rank + (key).k.index), PetscHashInt((key).l.rank + (key).l.index)))

     20: #define PetscHMapIJKLRemoteKeyEqual(k1, k2) \
     21:   (((k1).i.rank == (k2).i.rank) ? ((k1).i.index == (k2).i.index) ? ((k1).j.rank == (k2).j.rank) ? ((k1).j.index == (k2).j.index) ? ((k1).k.rank == (k2).k.rank) ? ((k1).k.index == (k2).k.index) ? ((k1).l.rank == (k2).l.rank) ? ((k1).l.index == (k2).l.index) : 0 : 0 : 0 : 0 : 0 : 0 : 0)

     23: PetscDisableStaticAnalyzerForExpressionUnderstandingThatThisIsDangerousAndBugprone(PETSC_HASH_MAP(HMapIJKLRemote, PetscHMapIJKLRemoteKey, PetscSFNode, PetscHMapIJKLRemoteKeyHash, PetscHMapIJKLRemoteKeyEqual, _PetscInvalidSFNode))

     25:   static PetscErrorCode PetscSortSFNode(PetscInt n, PetscSFNode A[])
     26: {
     27:   PetscInt i;

     29:   PetscFunctionBegin;
     30:   for (i = 1; i < n; ++i) {
     31:     PetscSFNode x = A[i];
     32:     PetscInt    j;

     34:     for (j = i - 1; j >= 0; --j) {
     35:       if ((A[j].rank > x.rank) || (A[j].rank == x.rank && A[j].index > x.index)) break;
     36:       A[j + 1] = A[j];
     37:     }
     38:     A[j + 1] = x;
     39:   }
     40:   PetscFunctionReturn(PETSC_SUCCESS);
     41: }

     43: /*
     44:   DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone
     45: */
     46: PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[])
     47: {
     48:   DMPolytopeType *typesTmp = NULL;
     49:   PetscInt       *sizesTmp = NULL, *facesTmp = NULL;
     50:   PetscInt       *tmp;
     51:   PetscInt        maxConeSize, maxSupportSize, maxSize;
     52:   PetscInt        getSize = 0;

     54:   PetscFunctionBegin;
     56:   if (cone) PetscAssertPointer(cone, 3);
     57:   PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize));
     58:   maxSize = PetscMax(maxConeSize, maxSupportSize);
     59:   if (faceTypes) getSize += maxSize;
     60:   if (faceSizes) getSize += maxSize;
     61:   if (faces) getSize += PetscSqr(maxSize);
     62:   PetscCall(DMGetWorkArray(dm, getSize, MPIU_INT, &tmp));
     63:   if (faceTypes) {
     64:     typesTmp = (DMPolytopeType *)tmp;
     65:     tmp += maxSize;
     66:   }
     67:   if (faceSizes) {
     68:     sizesTmp = tmp;
     69:     tmp += maxSize;
     70:   }
     71:   if (faces) facesTmp = tmp;
     72:   switch (ct) {
     73:   case DM_POLYTOPE_POINT:
     74:     if (numFaces) *numFaces = 0;
     75:     if (faceTypes) *faceTypes = typesTmp;
     76:     if (faceSizes) *faceSizes = sizesTmp;
     77:     if (faces) *faces = facesTmp;
     78:     break;
     79:   case DM_POLYTOPE_SEGMENT:
     80:     if (numFaces) *numFaces = 2;
     81:     if (faceTypes) {
     82:       typesTmp[0] = DM_POLYTOPE_POINT;
     83:       typesTmp[1] = DM_POLYTOPE_POINT;
     84:       *faceTypes  = typesTmp;
     85:     }
     86:     if (faceSizes) {
     87:       sizesTmp[0] = 1;
     88:       sizesTmp[1] = 1;
     89:       *faceSizes  = sizesTmp;
     90:     }
     91:     if (faces) {
     92:       facesTmp[0] = cone[0];
     93:       facesTmp[1] = cone[1];
     94:       *faces      = facesTmp;
     95:     }
     96:     break;
     97:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
     98:     if (numFaces) *numFaces = 2;
     99:     if (faceTypes) {
    100:       typesTmp[0] = DM_POLYTOPE_POINT;
    101:       typesTmp[1] = DM_POLYTOPE_POINT;
    102:       *faceTypes  = typesTmp;
    103:     }
    104:     if (faceSizes) {
    105:       sizesTmp[0] = 1;
    106:       sizesTmp[1] = 1;
    107:       *faceSizes  = sizesTmp;
    108:     }
    109:     if (faces) {
    110:       facesTmp[0] = cone[0];
    111:       facesTmp[1] = cone[1];
    112:       *faces      = facesTmp;
    113:     }
    114:     break;
    115:   case DM_POLYTOPE_TRIANGLE:
    116:     if (numFaces) *numFaces = 3;
    117:     if (faceTypes) {
    118:       typesTmp[0] = DM_POLYTOPE_SEGMENT;
    119:       typesTmp[1] = DM_POLYTOPE_SEGMENT;
    120:       typesTmp[2] = DM_POLYTOPE_SEGMENT;
    121:       *faceTypes  = typesTmp;
    122:     }
    123:     if (faceSizes) {
    124:       sizesTmp[0] = 2;
    125:       sizesTmp[1] = 2;
    126:       sizesTmp[2] = 2;
    127:       *faceSizes  = sizesTmp;
    128:     }
    129:     if (faces) {
    130:       facesTmp[0] = cone[0];
    131:       facesTmp[1] = cone[1];
    132:       facesTmp[2] = cone[1];
    133:       facesTmp[3] = cone[2];
    134:       facesTmp[4] = cone[2];
    135:       facesTmp[5] = cone[0];
    136:       *faces      = facesTmp;
    137:     }
    138:     break;
    139:   case DM_POLYTOPE_QUADRILATERAL:
    140:     /* Vertices follow right hand rule */
    141:     if (numFaces) *numFaces = 4;
    142:     if (faceTypes) {
    143:       typesTmp[0] = DM_POLYTOPE_SEGMENT;
    144:       typesTmp[1] = DM_POLYTOPE_SEGMENT;
    145:       typesTmp[2] = DM_POLYTOPE_SEGMENT;
    146:       typesTmp[3] = DM_POLYTOPE_SEGMENT;
    147:       *faceTypes  = typesTmp;
    148:     }
    149:     if (faceSizes) {
    150:       sizesTmp[0] = 2;
    151:       sizesTmp[1] = 2;
    152:       sizesTmp[2] = 2;
    153:       sizesTmp[3] = 2;
    154:       *faceSizes  = sizesTmp;
    155:     }
    156:     if (faces) {
    157:       facesTmp[0] = cone[0];
    158:       facesTmp[1] = cone[1];
    159:       facesTmp[2] = cone[1];
    160:       facesTmp[3] = cone[2];
    161:       facesTmp[4] = cone[2];
    162:       facesTmp[5] = cone[3];
    163:       facesTmp[6] = cone[3];
    164:       facesTmp[7] = cone[0];
    165:       *faces      = facesTmp;
    166:     }
    167:     break;
    168:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
    169:     if (numFaces) *numFaces = 4;
    170:     if (faceTypes) {
    171:       typesTmp[0] = DM_POLYTOPE_SEGMENT;
    172:       typesTmp[1] = DM_POLYTOPE_SEGMENT;
    173:       typesTmp[2] = DM_POLYTOPE_POINT_PRISM_TENSOR;
    174:       typesTmp[3] = DM_POLYTOPE_POINT_PRISM_TENSOR;
    175:       *faceTypes  = typesTmp;
    176:     }
    177:     if (faceSizes) {
    178:       sizesTmp[0] = 2;
    179:       sizesTmp[1] = 2;
    180:       sizesTmp[2] = 2;
    181:       sizesTmp[3] = 2;
    182:       *faceSizes  = sizesTmp;
    183:     }
    184:     if (faces) {
    185:       facesTmp[0] = cone[0];
    186:       facesTmp[1] = cone[1];
    187:       facesTmp[2] = cone[2];
    188:       facesTmp[3] = cone[3];
    189:       facesTmp[4] = cone[0];
    190:       facesTmp[5] = cone[2];
    191:       facesTmp[6] = cone[1];
    192:       facesTmp[7] = cone[3];
    193:       *faces      = facesTmp;
    194:     }
    195:     break;
    196:   case DM_POLYTOPE_TETRAHEDRON:
    197:     /* Vertices of first face follow right hand rule and normal points away from last vertex */
    198:     if (numFaces) *numFaces = 4;
    199:     if (faceTypes) {
    200:       typesTmp[0] = DM_POLYTOPE_TRIANGLE;
    201:       typesTmp[1] = DM_POLYTOPE_TRIANGLE;
    202:       typesTmp[2] = DM_POLYTOPE_TRIANGLE;
    203:       typesTmp[3] = DM_POLYTOPE_TRIANGLE;
    204:       *faceTypes  = typesTmp;
    205:     }
    206:     if (faceSizes) {
    207:       sizesTmp[0] = 3;
    208:       sizesTmp[1] = 3;
    209:       sizesTmp[2] = 3;
    210:       sizesTmp[3] = 3;
    211:       *faceSizes  = sizesTmp;
    212:     }
    213:     if (faces) {
    214:       facesTmp[0]  = cone[0];
    215:       facesTmp[1]  = cone[1];
    216:       facesTmp[2]  = cone[2];
    217:       facesTmp[3]  = cone[0];
    218:       facesTmp[4]  = cone[3];
    219:       facesTmp[5]  = cone[1];
    220:       facesTmp[6]  = cone[0];
    221:       facesTmp[7]  = cone[2];
    222:       facesTmp[8]  = cone[3];
    223:       facesTmp[9]  = cone[2];
    224:       facesTmp[10] = cone[1];
    225:       facesTmp[11] = cone[3];
    226:       *faces       = facesTmp;
    227:     }
    228:     break;
    229:   case DM_POLYTOPE_HEXAHEDRON:
    230:     /*  7--------6
    231:          /|       /|
    232:         / |      / |
    233:        4--------5  |
    234:        |  |     |  |
    235:        |  |     |  |
    236:        |  1--------2
    237:        | /      | /
    238:        |/       |/
    239:        0--------3
    240:        */
    241:     if (numFaces) *numFaces = 6;
    242:     if (faceTypes) {
    243:       typesTmp[0] = DM_POLYTOPE_QUADRILATERAL;
    244:       typesTmp[1] = DM_POLYTOPE_QUADRILATERAL;
    245:       typesTmp[2] = DM_POLYTOPE_QUADRILATERAL;
    246:       typesTmp[3] = DM_POLYTOPE_QUADRILATERAL;
    247:       typesTmp[4] = DM_POLYTOPE_QUADRILATERAL;
    248:       typesTmp[5] = DM_POLYTOPE_QUADRILATERAL;
    249:       *faceTypes  = typesTmp;
    250:     }
    251:     if (faceSizes) {
    252:       sizesTmp[0] = 4;
    253:       sizesTmp[1] = 4;
    254:       sizesTmp[2] = 4;
    255:       sizesTmp[3] = 4;
    256:       sizesTmp[4] = 4;
    257:       sizesTmp[5] = 4;
    258:       *faceSizes  = sizesTmp;
    259:     }
    260:     if (faces) {
    261:       facesTmp[0]  = cone[0];
    262:       facesTmp[1]  = cone[1];
    263:       facesTmp[2]  = cone[2];
    264:       facesTmp[3]  = cone[3]; /* Bottom */
    265:       facesTmp[4]  = cone[4];
    266:       facesTmp[5]  = cone[5];
    267:       facesTmp[6]  = cone[6];
    268:       facesTmp[7]  = cone[7]; /* Top */
    269:       facesTmp[8]  = cone[0];
    270:       facesTmp[9]  = cone[3];
    271:       facesTmp[10] = cone[5];
    272:       facesTmp[11] = cone[4]; /* Front */
    273:       facesTmp[12] = cone[2];
    274:       facesTmp[13] = cone[1];
    275:       facesTmp[14] = cone[7];
    276:       facesTmp[15] = cone[6]; /* Back */
    277:       facesTmp[16] = cone[3];
    278:       facesTmp[17] = cone[2];
    279:       facesTmp[18] = cone[6];
    280:       facesTmp[19] = cone[5]; /* Right */
    281:       facesTmp[20] = cone[0];
    282:       facesTmp[21] = cone[4];
    283:       facesTmp[22] = cone[7];
    284:       facesTmp[23] = cone[1]; /* Left */
    285:       *faces       = facesTmp;
    286:     }
    287:     break;
    288:   case DM_POLYTOPE_TRI_PRISM:
    289:     if (numFaces) *numFaces = 5;
    290:     if (faceTypes) {
    291:       typesTmp[0] = DM_POLYTOPE_TRIANGLE;
    292:       typesTmp[1] = DM_POLYTOPE_TRIANGLE;
    293:       typesTmp[2] = DM_POLYTOPE_QUADRILATERAL;
    294:       typesTmp[3] = DM_POLYTOPE_QUADRILATERAL;
    295:       typesTmp[4] = DM_POLYTOPE_QUADRILATERAL;
    296:       *faceTypes  = typesTmp;
    297:     }
    298:     if (faceSizes) {
    299:       sizesTmp[0] = 3;
    300:       sizesTmp[1] = 3;
    301:       sizesTmp[2] = 4;
    302:       sizesTmp[3] = 4;
    303:       sizesTmp[4] = 4;
    304:       *faceSizes  = sizesTmp;
    305:     }
    306:     if (faces) {
    307:       facesTmp[0]  = cone[0];
    308:       facesTmp[1]  = cone[1];
    309:       facesTmp[2]  = cone[2]; /* Bottom */
    310:       facesTmp[3]  = cone[3];
    311:       facesTmp[4]  = cone[4];
    312:       facesTmp[5]  = cone[5]; /* Top */
    313:       facesTmp[6]  = cone[0];
    314:       facesTmp[7]  = cone[2];
    315:       facesTmp[8]  = cone[4];
    316:       facesTmp[9]  = cone[3]; /* Back left */
    317:       facesTmp[10] = cone[2];
    318:       facesTmp[11] = cone[1];
    319:       facesTmp[12] = cone[5];
    320:       facesTmp[13] = cone[4]; /* Front */
    321:       facesTmp[14] = cone[1];
    322:       facesTmp[15] = cone[0];
    323:       facesTmp[16] = cone[3];
    324:       facesTmp[17] = cone[5]; /* Back right */
    325:       *faces       = facesTmp;
    326:     }
    327:     break;
    328:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
    329:     if (numFaces) *numFaces = 5;
    330:     if (faceTypes) {
    331:       typesTmp[0] = DM_POLYTOPE_TRIANGLE;
    332:       typesTmp[1] = DM_POLYTOPE_TRIANGLE;
    333:       typesTmp[2] = DM_POLYTOPE_SEG_PRISM_TENSOR;
    334:       typesTmp[3] = DM_POLYTOPE_SEG_PRISM_TENSOR;
    335:       typesTmp[4] = DM_POLYTOPE_SEG_PRISM_TENSOR;
    336:       *faceTypes  = typesTmp;
    337:     }
    338:     if (faceSizes) {
    339:       sizesTmp[0] = 3;
    340:       sizesTmp[1] = 3;
    341:       sizesTmp[2] = 4;
    342:       sizesTmp[3] = 4;
    343:       sizesTmp[4] = 4;
    344:       *faceSizes  = sizesTmp;
    345:     }
    346:     if (faces) {
    347:       facesTmp[0]  = cone[0];
    348:       facesTmp[1]  = cone[1];
    349:       facesTmp[2]  = cone[2]; /* Bottom */
    350:       facesTmp[3]  = cone[3];
    351:       facesTmp[4]  = cone[4];
    352:       facesTmp[5]  = cone[5]; /* Top */
    353:       facesTmp[6]  = cone[0];
    354:       facesTmp[7]  = cone[1];
    355:       facesTmp[8]  = cone[3];
    356:       facesTmp[9]  = cone[4]; /* Back left */
    357:       facesTmp[10] = cone[1];
    358:       facesTmp[11] = cone[2];
    359:       facesTmp[12] = cone[4];
    360:       facesTmp[13] = cone[5]; /* Back right */
    361:       facesTmp[14] = cone[2];
    362:       facesTmp[15] = cone[0];
    363:       facesTmp[16] = cone[5];
    364:       facesTmp[17] = cone[3]; /* Front */
    365:       *faces       = facesTmp;
    366:     }
    367:     break;
    368:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
    369:     /*  7--------6
    370:          /|       /|
    371:         / |      / |
    372:        4--------5  |
    373:        |  |     |  |
    374:        |  |     |  |
    375:        |  3--------2
    376:        | /      | /
    377:        |/       |/
    378:        0--------1
    379:        */
    380:     if (numFaces) *numFaces = 6;
    381:     if (faceTypes) {
    382:       typesTmp[0] = DM_POLYTOPE_QUADRILATERAL;
    383:       typesTmp[1] = DM_POLYTOPE_QUADRILATERAL;
    384:       typesTmp[2] = DM_POLYTOPE_SEG_PRISM_TENSOR;
    385:       typesTmp[3] = DM_POLYTOPE_SEG_PRISM_TENSOR;
    386:       typesTmp[4] = DM_POLYTOPE_SEG_PRISM_TENSOR;
    387:       typesTmp[5] = DM_POLYTOPE_SEG_PRISM_TENSOR;
    388:       *faceTypes  = typesTmp;
    389:     }
    390:     if (faceSizes) {
    391:       sizesTmp[0] = 4;
    392:       sizesTmp[1] = 4;
    393:       sizesTmp[2] = 4;
    394:       sizesTmp[3] = 4;
    395:       sizesTmp[4] = 4;
    396:       sizesTmp[5] = 4;
    397:       *faceSizes  = sizesTmp;
    398:     }
    399:     if (faces) {
    400:       facesTmp[0]  = cone[0];
    401:       facesTmp[1]  = cone[1];
    402:       facesTmp[2]  = cone[2];
    403:       facesTmp[3]  = cone[3]; /* Bottom */
    404:       facesTmp[4]  = cone[4];
    405:       facesTmp[5]  = cone[5];
    406:       facesTmp[6]  = cone[6];
    407:       facesTmp[7]  = cone[7]; /* Top */
    408:       facesTmp[8]  = cone[0];
    409:       facesTmp[9]  = cone[1];
    410:       facesTmp[10] = cone[4];
    411:       facesTmp[11] = cone[5]; /* Front */
    412:       facesTmp[12] = cone[1];
    413:       facesTmp[13] = cone[2];
    414:       facesTmp[14] = cone[5];
    415:       facesTmp[15] = cone[6]; /* Right */
    416:       facesTmp[16] = cone[2];
    417:       facesTmp[17] = cone[3];
    418:       facesTmp[18] = cone[6];
    419:       facesTmp[19] = cone[7]; /* Back */
    420:       facesTmp[20] = cone[3];
    421:       facesTmp[21] = cone[0];
    422:       facesTmp[22] = cone[7];
    423:       facesTmp[23] = cone[4]; /* Left */
    424:       *faces       = facesTmp;
    425:     }
    426:     break;
    427:   case DM_POLYTOPE_PYRAMID:
    428:     /*
    429:        4----
    430:        |\-\ \-----
    431:        | \ -\     \
    432:        |  1--\-----2
    433:        | /    \   /
    434:        |/      \ /
    435:        0--------3
    436:        */
    437:     if (numFaces) *numFaces = 5;
    438:     if (faceTypes) {
    439:       typesTmp[0] = DM_POLYTOPE_QUADRILATERAL;
    440:       typesTmp[1] = DM_POLYTOPE_TRIANGLE;
    441:       typesTmp[2] = DM_POLYTOPE_TRIANGLE;
    442:       typesTmp[3] = DM_POLYTOPE_TRIANGLE;
    443:       typesTmp[4] = DM_POLYTOPE_TRIANGLE;
    444:       *faceTypes  = typesTmp;
    445:     }
    446:     if (faceSizes) {
    447:       sizesTmp[0] = 4;
    448:       sizesTmp[1] = 3;
    449:       sizesTmp[2] = 3;
    450:       sizesTmp[3] = 3;
    451:       sizesTmp[4] = 3;
    452:       *faceSizes  = sizesTmp;
    453:     }
    454:     if (faces) {
    455:       facesTmp[0]  = cone[0];
    456:       facesTmp[1]  = cone[1];
    457:       facesTmp[2]  = cone[2];
    458:       facesTmp[3]  = cone[3]; /* Bottom */
    459:       facesTmp[4]  = cone[0];
    460:       facesTmp[5]  = cone[3];
    461:       facesTmp[6]  = cone[4]; /* Front */
    462:       facesTmp[7]  = cone[3];
    463:       facesTmp[8]  = cone[2];
    464:       facesTmp[9]  = cone[4]; /* Right */
    465:       facesTmp[10] = cone[2];
    466:       facesTmp[11] = cone[1];
    467:       facesTmp[12] = cone[4]; /* Back */
    468:       facesTmp[13] = cone[1];
    469:       facesTmp[14] = cone[0];
    470:       facesTmp[15] = cone[4]; /* Left */
    471:       *faces       = facesTmp;
    472:     }
    473:     break;
    474:   default:
    475:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No face description for cell type %s", DMPolytopeTypes[ct]);
    476:   }
    477:   PetscFunctionReturn(PETSC_SUCCESS);
    478: }

    480: PetscErrorCode DMPlexRestoreRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[])
    481: {
    482:   PetscFunctionBegin;
    483:   if (faceTypes) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faceTypes));
    484:   else if (faceSizes) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faceSizes));
    485:   else if (faces) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faces));
    486:   if (faceTypes) *faceTypes = NULL;
    487:   if (faceSizes) *faceSizes = NULL;
    488:   if (faces) *faces = NULL;
    489:   PetscFunctionReturn(PETSC_SUCCESS);
    490: }

    492: /* This interpolates faces for cells at some stratum */
    493: static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
    494: {
    495:   DMLabel       ctLabel;
    496:   PetscHMapIJKL faceTable;
    497:   PetscInt      faceTypeNum[DM_NUM_POLYTOPES];
    498:   PetscInt      depth, pStart, Np, cStart, cEnd, fStart, fEnd, vStart, vEnd;
    499:   PetscInt      cntFaces, *facesId, minCone;

    501:   PetscFunctionBegin;
    502:   PetscCall(DMPlexGetDepth(dm, &depth));
    503:   PetscCall(PetscHMapIJKLCreate(&faceTable));
    504:   PetscCall(PetscArrayzero(faceTypeNum, DM_NUM_POLYTOPES));
    505:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
    506:   PetscCall(DMPlexGetDepthStratum(dm, cellDepth, &cStart, &cEnd));
    507:   // If the range incorporates the vertices, it means we have a non-manifold topology, so choose just cells
    508:   if (cStart <= vStart && cEnd >= vEnd) cEnd = vStart;
    509:   // Number new faces and save face vertices in hash table
    510:   //   If depth > cellDepth, meaning we are interpolating faces, put the new (d-1)-faces after them
    511:   //   otherwise, we are interpolating cells, so put the faces after the vertices
    512:   PetscCall(DMPlexGetDepthStratum(dm, depth > cellDepth ? cellDepth : 0, NULL, &fStart));
    513:   fEnd = fStart;

    515:   minCone  = PETSC_INT_MAX;
    516:   cntFaces = 0;
    517:   for (PetscInt c = cStart; c < cEnd; ++c) {
    518:     const PetscInt *cone;
    519:     DMPolytopeType  ct;
    520:     PetscInt        numFaces = 0, coneSize;

    522:     PetscCall(DMPlexGetCellType(dm, c, &ct));
    523:     PetscCall(DMPlexGetCone(dm, c, &cone));
    524:     PetscCall(DMPlexGetConeSize(dm, c, &coneSize));
    525:     for (PetscInt j = 0; j < coneSize; j++) minCone = PetscMin(cone[j], minCone);
    526:     // Ignore faces since they are interpolated
    527:     if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, NULL, NULL, NULL));
    528:     cntFaces += numFaces;
    529:   }
    530:   // Encode so that we can use 0 as an excluded value, instead of PETSC_INT_MAX
    531:   minCone = -(minCone - 1);

    533:   PetscCall(PetscMalloc1(cntFaces, &facesId));

    535:   cntFaces = 0;
    536:   for (PetscInt c = cStart; c < cEnd; ++c) {
    537:     const PetscInt       *cone, *faceSizes, *faces;
    538:     const DMPolytopeType *faceTypes;
    539:     DMPolytopeType        ct;
    540:     PetscInt              numFaces, foff = 0;

    542:     PetscCall(DMPlexGetCellType(dm, c, &ct));
    543:     PetscCall(DMPlexGetCone(dm, c, &cone));
    544:     // Ignore faces since they are interpolated
    545:     if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) {
    546:       PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
    547:     } else {
    548:       numFaces = 0;
    549:     }
    550:     for (PetscInt cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
    551:       const PetscInt       faceSize = faceSizes[cf];
    552:       const DMPolytopeType faceType = faceTypes[cf];
    553:       const PetscInt      *face     = &faces[foff];
    554:       PetscHashIJKLKey     key;
    555:       PetscHashIter        iter;
    556:       PetscBool            missing;

    558:       PetscCheck(faceSize <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %" PetscInt_FMT " > 4", faceSize);
    559:       key.i = face[0] + minCone;
    560:       key.j = faceSize > 1 ? face[1] + minCone : 0;
    561:       key.k = faceSize > 2 ? face[2] + minCone : 0;
    562:       key.l = faceSize > 3 ? face[3] + minCone : 0;
    563:       PetscCall(PetscSortInt(faceSize, (PetscInt *)&key));
    564:       PetscCall(PetscHMapIJKLPut(faceTable, key, &iter, &missing));
    565:       if (missing) {
    566:         facesId[cntFaces] = fEnd;
    567:         PetscCall(PetscHMapIJKLIterSet(faceTable, iter, fEnd++));
    568:         ++faceTypeNum[faceType];
    569:       } else PetscCall(PetscHMapIJKLIterGet(faceTable, iter, &facesId[cntFaces]));
    570:       cntFaces++;
    571:     }
    572:     if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
    573:   }
    574:   /* We need to number faces contiguously among types */
    575:   {
    576:     PetscInt faceTypeStart[DM_NUM_POLYTOPES], ct, numFT = 0;

    578:     for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) {
    579:       if (faceTypeNum[ct]) ++numFT;
    580:       faceTypeStart[ct] = 0;
    581:     }
    582:     if (numFT > 1) {
    583:       PetscCall(PetscHMapIJKLClear(faceTable));
    584:       faceTypeStart[0] = fStart;
    585:       for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) faceTypeStart[ct] = faceTypeStart[ct - 1] + faceTypeNum[ct - 1];
    586:       cntFaces = 0;
    587:       for (PetscInt c = cStart; c < cEnd; ++c) {
    588:         const PetscInt       *cone, *faceSizes, *faces;
    589:         const DMPolytopeType *faceTypes;
    590:         DMPolytopeType        ct;
    591:         PetscInt              numFaces, foff = 0;

    593:         PetscCall(DMPlexGetCellType(dm, c, &ct));
    594:         PetscCall(DMPlexGetCone(dm, c, &cone));
    595:         if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) {
    596:           PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
    597:         } else {
    598:           numFaces = 0;
    599:         }
    600:         for (PetscInt cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
    601:           const PetscInt       faceSize = faceSizes[cf];
    602:           const DMPolytopeType faceType = faceTypes[cf];
    603:           const PetscInt      *face     = &faces[foff];
    604:           PetscHashIJKLKey     key;
    605:           PetscHashIter        iter;
    606:           PetscBool            missing;

    608:           key.i = face[0] + minCone;
    609:           key.j = faceSize > 1 ? face[1] + minCone : 0;
    610:           key.k = faceSize > 2 ? face[2] + minCone : 0;
    611:           key.l = faceSize > 3 ? face[3] + minCone : 0;
    612:           PetscCall(PetscSortInt(faceSize, (PetscInt *)&key));
    613:           PetscCall(PetscHMapIJKLPut(faceTable, key, &iter, &missing));
    614:           if (missing) {
    615:             facesId[cntFaces] = faceTypeStart[faceType];
    616:             PetscCall(PetscHMapIJKLIterSet(faceTable, iter, faceTypeStart[faceType]++));
    617:           } else PetscCall(PetscHMapIJKLIterGet(faceTable, iter, &facesId[cntFaces]));
    618:           cntFaces++;
    619:         }
    620:         if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
    621:       }
    622:       for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) {
    623:         PetscCheck(faceTypeStart[ct] == faceTypeStart[ct - 1] + faceTypeNum[ct], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent numbering for cell type %s, %" PetscInt_FMT " != %" PetscInt_FMT " + %" PetscInt_FMT, DMPolytopeTypes[ct], faceTypeStart[ct], faceTypeStart[ct - 1], faceTypeNum[ct]);
    624:       }
    625:     }
    626:   }
    627:   PetscCall(PetscHMapIJKLDestroy(&faceTable));

    629:   // Add new points, perhaps inserting into the numbering
    630:   PetscCall(DMPlexGetChart(dm, &pStart, &Np));
    631:   PetscCall(DMPlexSetChart(idm, pStart, Np + (fEnd - fStart)));
    632:   // Set cone sizes
    633:   //   Must create the celltype label here so that we do not automatically try to compute the types
    634:   PetscCall(DMCreateLabel(idm, "celltype"));
    635:   PetscCall(DMPlexGetCellTypeLabel(idm, &ctLabel));
    636:   for (PetscInt d = 0; d <= depth; ++d) {
    637:     DMPolytopeType ct;
    638:     PetscInt       coneSize, pStart, pEnd, poff = 0;

    640:     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
    641:     // Check for non-manifold condition
    642:     if (d == cellDepth) {
    643:       if (pEnd == cEnd) continue;
    644:       else pStart = vEnd;
    645:     }
    646:     // Account for insertion
    647:     if (pStart >= fStart) poff = fEnd - fStart;
    648:     for (PetscInt p = pStart; p < pEnd; ++p) {
    649:       PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
    650:       PetscCall(DMPlexSetConeSize(idm, p + poff, coneSize));
    651:       PetscCall(DMPlexGetCellType(dm, p, &ct));
    652:       PetscCall(DMPlexSetCellType(idm, p + poff, ct));
    653:     }
    654:   }
    655:   cntFaces = 0;
    656:   for (PetscInt c = cStart; c < cEnd; ++c) {
    657:     const PetscInt       *cone, *faceSizes;
    658:     const DMPolytopeType *faceTypes;
    659:     DMPolytopeType        ct;
    660:     PetscInt              numFaces, poff = 0;

    662:     PetscCall(DMPlexGetCellType(dm, c, &ct));
    663:     PetscCall(DMPlexGetCone(dm, c, &cone));
    664:     if (c >= fStart) poff = fEnd - fStart;
    665:     if (ct == DM_POLYTOPE_SEGMENT || ct == DM_POLYTOPE_POINT_PRISM_TENSOR) {
    666:       PetscCall(DMPlexSetCellType(idm, c + poff, ct));
    667:       PetscCall(DMPlexSetConeSize(idm, c + poff, 2));
    668:       continue;
    669:     }
    670:     PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, NULL));
    671:     PetscCall(DMPlexSetCellType(idm, c + poff, ct));
    672:     PetscCall(DMPlexSetConeSize(idm, c + poff, numFaces));
    673:     for (PetscInt cf = 0; cf < numFaces; ++cf) {
    674:       const PetscInt f        = facesId[cntFaces];
    675:       DMPolytopeType faceType = faceTypes[cf];
    676:       const PetscInt faceSize = faceSizes[cf];
    677:       PetscCall(DMPlexSetConeSize(idm, f, faceSize));
    678:       PetscCall(DMPlexSetCellType(idm, f, faceType));
    679:       cntFaces++;
    680:     }
    681:     PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, NULL));
    682:   }
    683:   PetscCall(DMSetUp(idm));
    684:   // Initialize cones so we do not need the bash table to tell us that a cone has been set
    685:   {
    686:     PetscSection cs;
    687:     PetscInt    *cones, csize;

    689:     PetscCall(DMPlexGetConeSection(idm, &cs));
    690:     PetscCall(DMPlexGetCones(idm, &cones));
    691:     PetscCall(PetscSectionGetStorageSize(cs, &csize));
    692:     for (PetscInt c = 0; c < csize; ++c) cones[c] = -1;
    693:   }
    694:   // Set cones
    695:   {
    696:     PetscInt *icone;
    697:     PetscInt  maxConeSize;

    699:     PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL));
    700:     PetscCall(PetscMalloc1(maxConeSize, &icone));
    701:     for (PetscInt d = 0; d <= depth; ++d) {
    702:       const PetscInt *cone;
    703:       PetscInt        pStart, pEnd, poff = 0, coneSize;

    705:       PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
    706:       // Check for non-manifold condition
    707:       if (d == cellDepth) {
    708:         if (pEnd == cEnd) continue;
    709:         else pStart = vEnd;
    710:       }
    711:       // Account for insertion
    712:       if (pStart >= fStart) poff = fEnd - fStart;
    713:       for (PetscInt p = pStart; p < pEnd; ++p) {
    714:         PetscCall(DMPlexGetCone(dm, p, &cone));
    715:         PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
    716:         for (PetscInt cp = 0; cp < coneSize; ++cp) icone[cp] = cone[cp] + (cone[cp] >= fStart ? fEnd - fStart : 0);
    717:         PetscCall(DMPlexSetCone(idm, p + poff, icone));
    718:         PetscCall(DMPlexGetConeOrientation(dm, p, &cone));
    719:         PetscCall(DMPlexSetConeOrientation(idm, p + poff, cone));
    720:       }
    721:     }
    722:     cntFaces = 0;
    723:     for (PetscInt c = cStart; c < cEnd; ++c) {
    724:       const PetscInt       *cone, *faceSizes, *faces;
    725:       const DMPolytopeType *faceTypes;
    726:       DMPolytopeType        ct;
    727:       PetscInt              coneSize, numFaces, foff = 0, poff = 0;

    729:       PetscCall(DMPlexGetCellType(dm, c, &ct));
    730:       PetscCall(DMPlexGetCone(dm, c, &cone));
    731:       PetscCall(DMPlexGetConeSize(dm, c, &coneSize));
    732:       if (c >= fStart) poff = fEnd - fStart;
    733:       if (ct == DM_POLYTOPE_SEGMENT || ct == DM_POLYTOPE_POINT_PRISM_TENSOR) {
    734:         for (PetscInt cp = 0; cp < coneSize; ++cp) icone[cp] = cone[cp] + (cone[cp] >= fStart ? fEnd - fStart : 0);
    735:         PetscCall(DMPlexSetCone(idm, c + poff, icone));
    736:         PetscCall(DMPlexGetConeOrientation(dm, c, &cone));
    737:         PetscCall(DMPlexSetConeOrientation(idm, c + poff, cone));
    738:         continue;
    739:       }
    740:       PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
    741:       for (PetscInt cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
    742:         DMPolytopeType  faceType = faceTypes[cf];
    743:         const PetscInt  faceSize = faceSizes[cf];
    744:         const PetscInt  f        = facesId[cntFaces];
    745:         const PetscInt *face     = &faces[foff];
    746:         const PetscInt *fcone;

    748:         PetscCall(DMPlexInsertCone(idm, c, cf, f));
    749:         PetscCall(DMPlexGetCone(idm, f, &fcone));
    750:         if (fcone[0] < 0) PetscCall(DMPlexSetCone(idm, f, face));
    751:         {
    752:           const PetscInt *fcone2;
    753:           PetscInt        ornt;

    755:           PetscCall(DMPlexGetConeSize(idm, f, &coneSize));
    756:           PetscCall(DMPlexGetCone(idm, f, &fcone2));
    757:           PetscCheck(coneSize == faceSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of face vertices %" PetscInt_FMT " for face %" PetscInt_FMT " should be %" PetscInt_FMT, coneSize, f, faceSize);
    758:           /* Notice that we have to use vertices here because the lower dimensional faces have not been created yet */
    759:           PetscCall(DMPolytopeGetVertexOrientation(faceType, fcone2, face, &ornt));
    760:           PetscCall(DMPlexInsertConeOrientation(idm, c + poff, cf, ornt));
    761:         }
    762:         cntFaces++;
    763:       }
    764:       PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces));
    765:     }
    766:     PetscCall(PetscFree(icone));
    767:   }
    768:   PetscCall(PetscFree(facesId));
    769:   PetscCall(DMPlexSymmetrize(idm));
    770:   PetscCall(DMPlexStratify(idm));
    771:   PetscFunctionReturn(PETSC_SUCCESS);
    772: }

    774: static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[])
    775: {
    776:   PetscInt           nleaves;
    777:   PetscMPIInt        nranks;
    778:   const PetscMPIInt *ranks   = NULL;
    779:   const PetscInt    *roffset = NULL, *rmine = NULL, *rremote = NULL;
    780:   PetscInt           n, o;

    782:   PetscFunctionBegin;
    783:   PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote));
    784:   nleaves = roffset[nranks];
    785:   PetscCall(PetscMalloc2(nleaves, rmine1, nleaves, rremote1));
    786:   for (PetscMPIInt r = 0; r < nranks; r++) {
    787:     /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote
    788:        - to unify order with the other side */
    789:     o = roffset[r];
    790:     n = roffset[r + 1] - o;
    791:     PetscCall(PetscArraycpy(&(*rmine1)[o], &rmine[o], n));
    792:     PetscCall(PetscArraycpy(&(*rremote1)[o], &rremote[o], n));
    793:     PetscCall(PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]));
    794:   }
    795:   PetscFunctionReturn(PETSC_SUCCESS);
    796: }

    798: PetscErrorCode DMPlexOrientInterface_Internal(DM dm)
    799: {
    800:   PetscSF            sf;
    801:   const PetscInt    *locals;
    802:   const PetscSFNode *remotes;
    803:   const PetscMPIInt *ranks;
    804:   const PetscInt    *roffset;
    805:   PetscInt          *rmine1, *rremote1; /* rmine and rremote copies simultaneously sorted by rank and rremote */
    806:   PetscInt           nroots, p, nleaves, maxConeSize = 0;
    807:   PetscMPIInt        nranks, r;
    808:   PetscInt(*roots)[4], (*leaves)[4], mainCone[4];
    809:   PetscMPIInt(*rootsRanks)[4], (*leavesRanks)[4];
    810:   MPI_Comm    comm;
    811:   PetscMPIInt rank, size;
    812:   PetscInt    debug = 0;

    814:   PetscFunctionBegin;
    815:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    816:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    817:   PetscCallMPI(MPI_Comm_size(comm, &size));
    818:   PetscCall(DMGetPointSF(dm, &sf));
    819:   PetscCall(DMViewFromOptions(dm, NULL, "-before_orient_interface_dm_view"));
    820:   if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sf, PETSC_FALSE));
    821:   PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes));
    822:   if (nroots < 0) PetscFunctionReturn(PETSC_SUCCESS);
    823:   PetscCall(PetscSFSetUp(sf));
    824:   PetscCall(SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1));
    825:   for (p = 0; p < nleaves; ++p) {
    826:     PetscInt coneSize;
    827:     PetscCall(DMPlexGetConeSize(dm, locals[p], &coneSize));
    828:     maxConeSize = PetscMax(maxConeSize, coneSize);
    829:   }
    830:   PetscCheck(maxConeSize <= 4, comm, PETSC_ERR_SUP, "This method does not support cones of size %" PetscInt_FMT, maxConeSize);
    831:   PetscCall(PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks));
    832:   for (p = 0; p < nroots; ++p) {
    833:     const PetscInt *cone;
    834:     PetscInt        coneSize, c, ind0;

    836:     PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
    837:     PetscCall(DMPlexGetCone(dm, p, &cone));
    838:     /* Ignore vertices */
    839:     if (coneSize < 2) {
    840:       for (c = 0; c < 4; c++) {
    841:         roots[p][c]      = -1;
    842:         rootsRanks[p][c] = -1;
    843:       }
    844:       continue;
    845:     }
    846:     /* Translate all points to root numbering */
    847:     for (c = 0; c < PetscMin(coneSize, 4); c++) {
    848:       PetscCall(PetscFindInt(cone[c], nleaves, locals, &ind0));
    849:       if (ind0 < 0) {
    850:         roots[p][c]      = cone[c];
    851:         rootsRanks[p][c] = rank;
    852:       } else {
    853:         roots[p][c]      = remotes[ind0].index;
    854:         rootsRanks[p][c] = (PetscMPIInt)remotes[ind0].rank;
    855:       }
    856:     }
    857:     for (c = coneSize; c < 4; c++) {
    858:       roots[p][c]      = -1;
    859:       rootsRanks[p][c] = -1;
    860:     }
    861:   }
    862:   for (p = 0; p < nroots; ++p) {
    863:     PetscInt c;
    864:     for (c = 0; c < 4; c++) {
    865:       leaves[p][c]      = -2;
    866:       leavesRanks[p][c] = -2;
    867:     }
    868:   }
    869:   PetscCall(PetscSFBcastBegin(sf, MPIU_4INT, roots, leaves, MPI_REPLACE));
    870:   PetscCall(PetscSFBcastBegin(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE));
    871:   PetscCall(PetscSFBcastEnd(sf, MPIU_4INT, roots, leaves, MPI_REPLACE));
    872:   PetscCall(PetscSFBcastEnd(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE));
    873:   if (debug) {
    874:     PetscCall(PetscSynchronizedFlush(comm, NULL));
    875:     if (rank == 0) PetscCall(PetscSynchronizedPrintf(comm, "Referenced roots\n"));
    876:   }
    877:   PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL));
    878:   for (p = 0; p < nroots; ++p) {
    879:     DMPolytopeType  ct;
    880:     const PetscInt *cone;
    881:     PetscInt        coneSize, c, ind0, o, ir;

    883:     if (leaves[p][0] < 0) continue; /* Ignore vertices */
    884:     PetscCall(DMPlexGetCellType(dm, p, &ct));
    885:     PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
    886:     PetscCall(DMPlexGetCone(dm, p, &cone));
    887:     if (debug) {
    888:       PetscCall(PetscSynchronizedPrintf(comm, "[%d]  %4" PetscInt_FMT ": cone=[%4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT "] roots=[(%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ")] leaves=[(%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ")]", rank, p, cone[0], cone[1], cone[2], cone[3], rootsRanks[p][0], roots[p][0], rootsRanks[p][1], roots[p][1], rootsRanks[p][2], roots[p][2], rootsRanks[p][3], roots[p][3], leavesRanks[p][0], leaves[p][0], leavesRanks[p][1], leaves[p][1], leavesRanks[p][2], leaves[p][2], leavesRanks[p][3], leaves[p][3]));
    889:     }
    890:     if (leavesRanks[p][0] != rootsRanks[p][0] || leaves[p][0] != roots[p][0] || leavesRanks[p][1] != rootsRanks[p][1] || leaves[p][1] != roots[p][1] || leavesRanks[p][2] != rootsRanks[p][2] || leaves[p][2] != roots[p][2] || leavesRanks[p][3] != rootsRanks[p][3] || leaves[p][3] != roots[p][3]) {
    891:       /* Translate these leaves to my cone points; mainCone means desired order p's cone points */
    892:       for (c = 0; c < PetscMin(coneSize, 4); ++c) {
    893:         PetscInt rS, rN;

    895:         if (leavesRanks[p][c] == rank) {
    896:           /* A local leaf is just taken as it is */
    897:           mainCone[c] = leaves[p][c];
    898:           continue;
    899:         }
    900:         /* Find index of rank leavesRanks[p][c] among remote ranks */
    901:         /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
    902:         PetscCall(PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &ir));
    903:         PetscCall(PetscMPIIntCast(ir, &r));
    904:         PetscCheck(r >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " cone[%" PetscInt_FMT "]=%" PetscInt_FMT " root (%d,%" PetscInt_FMT ") leaf (%d,%" PetscInt_FMT "): leaf rank not found among remote ranks", p, c, cone[c], rootsRanks[p][c], roots[p][c], leavesRanks[p][c], leaves[p][c]);
    905:         PetscCheck(ranks[r] >= 0 && ranks[r] < size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "p=%" PetscInt_FMT " c=%" PetscInt_FMT " commsize=%d: ranks[%d] = %d makes no sense", p, c, size, r, ranks[r]);
    906:         /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */
    907:         rS = roffset[r];
    908:         rN = roffset[r + 1] - rS;
    909:         PetscCall(PetscFindInt(leaves[p][c], rN, &rremote1[rS], &ind0));
    910:         PetscCheck(ind0 >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " cone[%" PetscInt_FMT "]=%" PetscInt_FMT " root (%d,%" PetscInt_FMT ") leave (%d,%" PetscInt_FMT "): corresponding remote point not found - it seems there is missing connection in point SF!", p, c, cone[c], rootsRanks[p][c], roots[p][c], leavesRanks[p][c], leaves[p][c]);
    911:         /* Get the corresponding local point */
    912:         mainCone[c] = rmine1[rS + ind0];
    913:       }
    914:       if (debug) PetscCall(PetscSynchronizedPrintf(comm, " mainCone=[%4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT "]\n", mainCone[0], mainCone[1], mainCone[2], mainCone[3]));
    915:       /* Set the desired order of p's cone points and fix orientations accordingly */
    916:       PetscCall(DMPolytopeGetOrientation(ct, cone, mainCone, &o));
    917:       PetscCall(DMPlexOrientPoint(dm, p, o));
    918:     } else if (debug) PetscCall(PetscSynchronizedPrintf(comm, " ==\n"));
    919:   }
    920:   if (debug) {
    921:     PetscCall(PetscSynchronizedFlush(comm, NULL));
    922:     PetscCallMPI(MPI_Barrier(comm));
    923:   }
    924:   PetscCall(DMViewFromOptions(dm, NULL, "-after_orient_interface_dm_view"));
    925:   PetscCall(PetscFree4(roots, leaves, rootsRanks, leavesRanks));
    926:   PetscCall(PetscFree2(rmine1, rremote1));
    927:   PetscFunctionReturn(PETSC_SUCCESS);
    928: }

    930: static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[])
    931: {
    932:   PetscInt    idx;
    933:   PetscMPIInt rank;
    934:   PetscBool   flg;

    936:   PetscFunctionBegin;
    937:   PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg));
    938:   if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
    939:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    940:   PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name));
    941:   for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s %" PetscInt_FMT " %s %" PetscInt_FMT "\n", rank, idxname, idx, valname, a[idx]));
    942:   PetscCall(PetscSynchronizedFlush(comm, NULL));
    943:   PetscFunctionReturn(PETSC_SUCCESS);
    944: }

    946: static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[])
    947: {
    948:   PetscInt    idx;
    949:   PetscMPIInt rank;
    950:   PetscBool   flg;

    952:   PetscFunctionBegin;
    953:   PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg));
    954:   if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
    955:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    956:   PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name));
    957:   if (idxname) {
    958:     for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s %" PetscInt_FMT " rank %d index %" PetscInt_FMT "\n", rank, idxname, idx, (PetscMPIInt)a[idx].rank, a[idx].index));
    959:   } else {
    960:     for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]rank %d index %" PetscInt_FMT "\n", rank, (PetscMPIInt)a[idx].rank, a[idx].index));
    961:   }
    962:   PetscCall(PetscSynchronizedFlush(comm, NULL));
    963:   PetscFunctionReturn(PETSC_SUCCESS);
    964: }

    966: static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint, PetscBool *mapFailed)
    967: {
    968:   PetscSF         sf;
    969:   const PetscInt *locals;
    970:   PetscMPIInt     rank;

    972:   PetscFunctionBegin;
    973:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    974:   PetscCall(DMGetPointSF(dm, &sf));
    975:   PetscCall(PetscSFGetGraph(sf, NULL, NULL, &locals, NULL));
    976:   if (mapFailed) *mapFailed = PETSC_FALSE;
    977:   if (remotePoint.rank == rank) {
    978:     *localPoint = remotePoint.index;
    979:   } else {
    980:     PetscHashIJKey key;
    981:     PetscInt       l;

    983:     key.i = remotePoint.index;
    984:     key.j = remotePoint.rank;
    985:     PetscCall(PetscHMapIJGet(remotehash, key, &l));
    986:     if (l >= 0) {
    987:       *localPoint = locals[l];
    988:     } else if (mapFailed) *mapFailed = PETSC_TRUE;
    989:   }
    990:   PetscFunctionReturn(PETSC_SUCCESS);
    991: }

    993: static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint, PetscBool *mapFailed)
    994: {
    995:   PetscSF            sf;
    996:   const PetscInt    *locals, *rootdegree;
    997:   const PetscSFNode *remotes;
    998:   PetscInt           Nl, l;
    999:   PetscMPIInt        rank;

    1001:   PetscFunctionBegin;
    1002:   if (mapFailed) *mapFailed = PETSC_FALSE;
    1003:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    1004:   PetscCall(DMGetPointSF(dm, &sf));
    1005:   PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes));
    1006:   if (Nl < 0) goto owned;
    1007:   PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
    1008:   PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
    1009:   if (rootdegree[localPoint]) goto owned;
    1010:   PetscCall(PetscFindInt(localPoint, Nl, locals, &l));
    1011:   if (l < 0) {
    1012:     if (mapFailed) *mapFailed = PETSC_TRUE;
    1013:   } else *remotePoint = remotes[l];
    1014:   PetscFunctionReturn(PETSC_SUCCESS);
    1015: owned:
    1016:   remotePoint->rank  = rank;
    1017:   remotePoint->index = localPoint;
    1018:   PetscFunctionReturn(PETSC_SUCCESS);
    1019: }

    1021: static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared)
    1022: {
    1023:   PetscSF         sf;
    1024:   const PetscInt *locals, *rootdegree;
    1025:   PetscInt        Nl, idx;

    1027:   PetscFunctionBegin;
    1028:   *isShared = PETSC_FALSE;
    1029:   PetscCall(DMGetPointSF(dm, &sf));
    1030:   PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL));
    1031:   if (Nl < 0) PetscFunctionReturn(PETSC_SUCCESS);
    1032:   PetscCall(PetscFindInt(p, Nl, locals, &idx));
    1033:   if (idx >= 0) {
    1034:     *isShared = PETSC_TRUE;
    1035:     PetscFunctionReturn(PETSC_SUCCESS);
    1036:   }
    1037:   PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
    1038:   PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
    1039:   if (rootdegree[p] > 0) *isShared = PETSC_TRUE;
    1040:   PetscFunctionReturn(PETSC_SUCCESS);
    1041: }

    1043: static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared)
    1044: {
    1045:   const PetscInt *cone;
    1046:   PetscInt        coneSize, c;
    1047:   PetscBool       cShared = PETSC_TRUE;

    1049:   PetscFunctionBegin;
    1050:   PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
    1051:   PetscCall(DMPlexGetCone(dm, p, &cone));
    1052:   for (c = 0; c < coneSize; ++c) {
    1053:     PetscBool pointShared;

    1055:     PetscCall(DMPlexPointIsShared(dm, cone[c], &pointShared));
    1056:     cShared = (PetscBool)(cShared && pointShared);
    1057:   }
    1058:   *isShared = coneSize ? cShared : PETSC_FALSE;
    1059:   PetscFunctionReturn(PETSC_SUCCESS);
    1060: }

    1062: static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin)
    1063: {
    1064:   const PetscInt *cone;
    1065:   PetscInt        coneSize, c;
    1066:   PetscSFNode     cmin = {PETSC_INT_MAX, PETSC_MPI_INT_MAX}, missing = {-1, -1};

    1068:   PetscFunctionBegin;
    1069:   PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
    1070:   PetscCall(DMPlexGetCone(dm, p, &cone));
    1071:   for (c = 0; c < coneSize; ++c) {
    1072:     PetscSFNode rcp;
    1073:     PetscBool   mapFailed;

    1075:     PetscCall(DMPlexMapToGlobalPoint(dm, cone[c], &rcp, &mapFailed));
    1076:     if (mapFailed) {
    1077:       cmin = missing;
    1078:     } else {
    1079:       cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin;
    1080:     }
    1081:   }
    1082:   *cpmin = coneSize ? cmin : missing;
    1083:   PetscFunctionReturn(PETSC_SUCCESS);
    1084: }

    1086: /*
    1087:   Each shared face has an entry in the candidates array:
    1088:     (-1, coneSize-1), {(global cone point)}
    1089:   where the set is missing the point p which we use as the key for the face
    1090: */
    1091: static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug)
    1092: {
    1093:   MPI_Comm        comm;
    1094:   const PetscInt *support;
    1095:   PetscInt        supportSize, s, off = 0, idx = 0, overlap, cellHeight, height;
    1096:   PetscMPIInt     rank;

    1098:   PetscFunctionBegin;
    1099:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    1100:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1101:   PetscCall(DMPlexGetOverlap(dm, &overlap));
    1102:   PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
    1103:   PetscCall(DMPlexGetPointHeight(dm, p, &height));
    1104:   if (!overlap && height <= cellHeight + 1) {
    1105:     /* cells can't be shared for non-overlapping meshes */
    1106:     if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Skipping face %" PetscInt_FMT " to avoid adding cell to hashmap since this is nonoverlapping mesh\n", rank, p));
    1107:     PetscFunctionReturn(PETSC_SUCCESS);
    1108:   }
    1109:   PetscCall(DMPlexGetSupportSize(dm, p, &supportSize));
    1110:   PetscCall(DMPlexGetSupport(dm, p, &support));
    1111:   if (candidates) PetscCall(PetscSectionGetOffset(candidateSection, p, &off));
    1112:   for (s = 0; s < supportSize; ++s) {
    1113:     const PetscInt  face = support[s];
    1114:     const PetscInt *cone;
    1115:     PetscSFNode     cpmin = {-1, -1}, rp = {-1, -1};
    1116:     PetscInt        coneSize, c, f;
    1117:     PetscBool       isShared = PETSC_FALSE;
    1118:     PetscHashIJKey  key;

    1120:     /* Only add point once */
    1121:     if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Support face %" PetscInt_FMT "\n", rank, face));
    1122:     key.i = p;
    1123:     key.j = face;
    1124:     PetscCall(PetscHMapIJGet(faceHash, key, &f));
    1125:     if (f >= 0) continue;
    1126:     PetscCall(DMPlexConeIsShared(dm, face, &isShared));
    1127:     PetscCall(DMPlexGetConeMinimum(dm, face, &cpmin));
    1128:     PetscCall(DMPlexMapToGlobalPoint(dm, p, &rp, NULL));
    1129:     if (debug) {
    1130:       PetscCall(PetscSynchronizedPrintf(comm, "[%d]      Face point %" PetscInt_FMT " is shared: %d\n", rank, face, (int)isShared));
    1131:       PetscCall(PetscSynchronizedPrintf(comm, "[%d]      Global point (%d, %" PetscInt_FMT ") Min Cone Point (%d, %" PetscInt_FMT ")\n", rank, (PetscMPIInt)rp.rank, rp.index, (PetscMPIInt)cpmin.rank, cpmin.index));
    1132:     }
    1133:     if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) {
    1134:       PetscCall(PetscHMapIJSet(faceHash, key, p));
    1135:       if (candidates) {
    1136:         if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Adding shared face %" PetscInt_FMT " at idx %" PetscInt_FMT "\n[%d]     ", rank, face, idx, rank));
    1137:         PetscCall(DMPlexGetConeSize(dm, face, &coneSize));
    1138:         PetscCall(DMPlexGetCone(dm, face, &cone));
    1139:         candidates[off + idx].rank    = -1;
    1140:         candidates[off + idx++].index = coneSize - 1;
    1141:         candidates[off + idx].rank    = rank;
    1142:         candidates[off + idx++].index = face;
    1143:         for (c = 0; c < coneSize; ++c) {
    1144:           const PetscInt cp = cone[c];

    1146:           if (cp == p) continue;
    1147:           PetscCall(DMPlexMapToGlobalPoint(dm, cp, &candidates[off + idx], NULL));
    1148:           if (debug) PetscCall(PetscSynchronizedPrintf(comm, " (%d,%" PetscInt_FMT ")", (PetscMPIInt)candidates[off + idx].rank, candidates[off + idx].index));
    1149:           ++idx;
    1150:         }
    1151:         if (debug) PetscCall(PetscSynchronizedPrintf(comm, "\n"));
    1152:       } else {
    1153:         /* Add cone size to section */
    1154:         if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Scheduling shared face %" PetscInt_FMT "\n", rank, face));
    1155:         PetscCall(DMPlexGetConeSize(dm, face, &coneSize));
    1156:         PetscCall(PetscHMapIJSet(faceHash, key, p));
    1157:         PetscCall(PetscSectionAddDof(candidateSection, p, coneSize + 1));
    1158:       }
    1159:     }
    1160:   }
    1161:   PetscFunctionReturn(PETSC_SUCCESS);
    1162: }

    1164: /*@
    1165:   DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the `PointSF` in parallel, following local interpolation

    1167:   Collective

    1169:   Input Parameters:
    1170: + dm      - The interpolated `DMPLEX`
    1171: - pointSF - The initial `PetscSF` without interpolated points

    1173:   Level: developer

    1175:   Note:
    1176:   Debugging for this process can be turned on with the options: `-dm_interp_pre_view` `-petscsf_interp_pre_view` `-petscsection_interp_candidate_view` `-petscsection_interp_candidate_remote_view` `-petscsection_interp_claim_view` `-petscsf_interp_pre_view` `-dmplex_interp_debug`

    1178: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexUninterpolate()`
    1179: @*/
    1180: PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
    1181: {
    1182:   MPI_Comm           comm;
    1183:   PetscHMapIJ        remoteHash;
    1184:   PetscHMapI         claimshash;
    1185:   PetscSection       candidateSection, candidateRemoteSection, claimSection;
    1186:   PetscSFNode       *candidates, *candidatesRemote, *claims;
    1187:   const PetscInt    *localPoints, *rootdegree;
    1188:   const PetscSFNode *remotePoints;
    1189:   PetscInt           ov, Nr, r, Nl, l;
    1190:   PetscInt           candidatesSize, candidatesRemoteSize, claimsSize;
    1191:   PetscBool          flg, debug = PETSC_FALSE;
    1192:   PetscMPIInt        rank;

    1194:   PetscFunctionBegin;
    1197:   PetscCall(DMPlexIsDistributed(dm, &flg));
    1198:   if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
    1199:   /* Set initial SF so that lower level queries work */
    1200:   PetscCall(DMSetPointSF(dm, pointSF));
    1201:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    1202:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1203:   PetscCall(DMPlexGetOverlap(dm, &ov));
    1204:   PetscCheck(!ov, comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet");
    1205:   PetscCall(PetscOptionsHasName(NULL, ((PetscObject)dm)->prefix, "-dmplex_interp_debug", &debug));
    1206:   PetscCall(PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_interp_pre_view"));
    1207:   PetscCall(PetscObjectViewFromOptions((PetscObject)pointSF, NULL, "-petscsf_interp_pre_view"));
    1208:   PetscCall(PetscLogEventBegin(DMPLEX_InterpolateSF, dm, 0, 0, 0));
    1209:   /* Step 0: Precalculations */
    1210:   PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints));
    1211:   PetscCheck(Nr >= 0, comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set");
    1212:   PetscCall(PetscHMapIJCreate(&remoteHash));
    1213:   for (l = 0; l < Nl; ++l) {
    1214:     PetscHashIJKey key;

    1216:     key.i = remotePoints[l].index;
    1217:     key.j = remotePoints[l].rank;
    1218:     PetscCall(PetscHMapIJSet(remoteHash, key, l));
    1219:   }
    1220:   /*   Compute root degree to identify shared points */
    1221:   PetscCall(PetscSFComputeDegreeBegin(pointSF, &rootdegree));
    1222:   PetscCall(PetscSFComputeDegreeEnd(pointSF, &rootdegree));
    1223:   PetscCall(IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree));
    1224:   /*
    1225:   1) Loop over each leaf point $p$ at depth $d$ in the SF
    1226:   \item Get set $F(p)$ of faces $f$ in the support of $p$ for which
    1227:   \begin{itemize}
    1228:     \item all cone points of $f$ are shared
    1229:     \item $p$ is the cone point with smallest canonical number
    1230:   \end{itemize}
    1231:   \item Send $F(p)$ and the cone of each face to the active root point $r(p)$
    1232:   \item At the root, if at least two faces with a given cone are present, including a local face, mark the face as shared \label{alg:rootStep} and choose the root face
    1233:   \item Send the root face from the root back to all leaf process
    1234:   \item Leaf processes add the shared face to the SF
    1235:   */
    1236:   /* Step 1: Construct section+SFNode array
    1237:        The section has entries for all shared faces for which we have a leaf point in the cone
    1238:        The array holds candidate shared faces, each face is referred to by the leaf point */
    1239:   PetscCall(PetscSectionCreate(comm, &candidateSection));
    1240:   PetscCall(PetscSectionSetChart(candidateSection, 0, Nr));
    1241:   {
    1242:     PetscHMapIJ faceHash;

    1244:     PetscCall(PetscHMapIJCreate(&faceHash));
    1245:     for (l = 0; l < Nl; ++l) {
    1246:       const PetscInt p = localPoints[l];

    1248:       if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]  First pass leaf point %" PetscInt_FMT "\n", rank, p));
    1249:       PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug));
    1250:     }
    1251:     PetscCall(PetscHMapIJClear(faceHash));
    1252:     PetscCall(PetscSectionSetUp(candidateSection));
    1253:     PetscCall(PetscSectionGetStorageSize(candidateSection, &candidatesSize));
    1254:     PetscCall(PetscMalloc1(candidatesSize, &candidates));
    1255:     for (l = 0; l < Nl; ++l) {
    1256:       const PetscInt p = localPoints[l];

    1258:       if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]  Second pass leaf point %" PetscInt_FMT "\n", rank, p));
    1259:       PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug));
    1260:     }
    1261:     PetscCall(PetscHMapIJDestroy(&faceHash));
    1262:     if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL));
    1263:   }
    1264:   PetscCall(PetscObjectSetName((PetscObject)candidateSection, "Candidate Section"));
    1265:   PetscCall(PetscObjectViewFromOptions((PetscObject)candidateSection, NULL, "-petscsection_interp_candidate_view"));
    1266:   PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates));
    1267:   /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
    1268:   /*   Note that this section is indexed by offsets into leaves, not by point number */
    1269:   {
    1270:     PetscSF   sfMulti, sfInverse, sfCandidates;
    1271:     PetscInt *remoteOffsets;

    1273:     PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti));
    1274:     PetscCall(PetscSFCreateInverseSF(sfMulti, &sfInverse));
    1275:     PetscCall(PetscSectionCreate(comm, &candidateRemoteSection));
    1276:     PetscCall(PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection));
    1277:     PetscCall(PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates));
    1278:     PetscCall(PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize));
    1279:     PetscCall(PetscMalloc1(candidatesRemoteSize, &candidatesRemote));
    1280:     PetscCall(PetscSFBcastBegin(sfCandidates, MPIU_SF_NODE, candidates, candidatesRemote, MPI_REPLACE));
    1281:     PetscCall(PetscSFBcastEnd(sfCandidates, MPIU_SF_NODE, candidates, candidatesRemote, MPI_REPLACE));
    1282:     PetscCall(PetscSFDestroy(&sfInverse));
    1283:     PetscCall(PetscSFDestroy(&sfCandidates));
    1284:     PetscCall(PetscFree(remoteOffsets));

    1286:     PetscCall(PetscObjectSetName((PetscObject)candidateRemoteSection, "Remote Candidate Section"));
    1287:     PetscCall(PetscObjectViewFromOptions((PetscObject)candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view"));
    1288:     PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote));
    1289:   }
    1290:   /* Step 3: At the root, if at least two faces with a given cone are present, including a local face, mark the face as shared and choose the root face */
    1291:   {
    1292:     PetscHMapIJKLRemote faceTable;
    1293:     PetscInt            idx, idx2;

    1295:     PetscCall(PetscHMapIJKLRemoteCreate(&faceTable));
    1296:     /* There is a section point for every leaf attached to a given root point */
    1297:     for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) {
    1298:       PetscInt deg;

    1300:       for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) {
    1301:         PetscInt offset, dof, d;

    1303:         PetscCall(PetscSectionGetDof(candidateRemoteSection, idx, &dof));
    1304:         PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx, &offset));
    1305:         /* dof may include many faces from the remote process */
    1306:         for (d = 0; d < dof; ++d) {
    1307:           const PetscInt         hidx  = offset + d;
    1308:           const PetscInt         Np    = candidatesRemote[hidx].index + 1;
    1309:           const PetscSFNode      rface = candidatesRemote[hidx + 1];
    1310:           const PetscSFNode     *fcone = &candidatesRemote[hidx + 2];
    1311:           PetscSFNode            fcp0;
    1312:           const PetscSFNode      pmax = {-1, -1};
    1313:           const PetscInt        *join = NULL;
    1314:           PetscHMapIJKLRemoteKey key;
    1315:           PetscHashIter          iter;
    1316:           PetscBool              missing, mapToLocalPointFailed = PETSC_FALSE;
    1317:           PetscInt               points[1024], p, joinSize;

    1319:           if (debug)
    1320:             PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Checking face (%d, %" PetscInt_FMT ") at (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") with cone size %" PetscInt_FMT "\n", rank, (PetscMPIInt)rface.rank,
    1321:                                               rface.index, r, idx, d, Np));

    1323:           PetscCheck(Np <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle face (%d, %" PetscInt_FMT ") at (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") with %" PetscInt_FMT " cone points", (PetscMPIInt)rface.rank, rface.index, r, idx, d, Np);
    1324:           fcp0.rank  = rank;
    1325:           fcp0.index = r;
    1326:           d += Np;
    1327:           /* Put remote face in hash table */
    1328:           key.i = fcp0;
    1329:           key.j = fcone[0];
    1330:           key.k = Np > 2 ? fcone[1] : pmax;
    1331:           key.l = Np > 3 ? fcone[2] : pmax;
    1332:           PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key));
    1333:           PetscCall(PetscHMapIJKLRemotePut(faceTable, key, &iter, &missing));
    1334:           if (missing) {
    1335:             if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Setting remote face (%" PetscInt_FMT ", %d)\n", rank, rface.index, (PetscMPIInt)rface.rank));
    1336:             PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, rface));
    1337:           } else {
    1338:             PetscSFNode oface;

    1340:             PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &oface));
    1341:             if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) {
    1342:               if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Replacing with remote face (%" PetscInt_FMT ", %d)\n", rank, rface.index, (PetscMPIInt)rface.rank));
    1343:               PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, rface));
    1344:             }
    1345:           }
    1346:           /* Check for local face */
    1347:           points[0] = r;
    1348:           for (p = 1; p < Np; ++p) {
    1349:             PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, fcone[p - 1], &points[p], &mapToLocalPointFailed));
    1350:             if (mapToLocalPointFailed) break; /* We got a point not in our overlap */
    1351:             if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Checking local candidate %" PetscInt_FMT "\n", rank, points[p]));
    1352:           }
    1353:           if (mapToLocalPointFailed) continue;
    1354:           PetscCall(DMPlexGetJoin(dm, Np, points, &joinSize, &join));
    1355:           if (joinSize == 1) {
    1356:             PetscSFNode lface;
    1357:             PetscSFNode oface;

    1359:             /* Always replace with local face */
    1360:             lface.rank  = rank;
    1361:             lface.index = join[0];
    1362:             PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &oface));
    1363:             if (debug)
    1364:               PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Replacing (%" PetscInt_FMT ", %d) with local face (%" PetscInt_FMT ", %d)\n", rank, oface.index, (PetscMPIInt)oface.rank, lface.index, (PetscMPIInt)lface.rank));
    1365:             PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, lface));
    1366:           }
    1367:           PetscCall(DMPlexRestoreJoin(dm, Np, points, &joinSize, &join));
    1368:         }
    1369:       }
    1370:       /* Put back faces for this root */
    1371:       for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) {
    1372:         PetscInt offset, dof, d;

    1374:         PetscCall(PetscSectionGetDof(candidateRemoteSection, idx2, &dof));
    1375:         PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx2, &offset));
    1376:         /* dof may include many faces from the remote process */
    1377:         for (d = 0; d < dof; ++d) {
    1378:           const PetscInt         hidx  = offset + d;
    1379:           const PetscInt         Np    = candidatesRemote[hidx].index + 1;
    1380:           const PetscSFNode     *fcone = &candidatesRemote[hidx + 2];
    1381:           PetscSFNode            fcp0;
    1382:           const PetscSFNode      pmax = {-1, -1};
    1383:           PetscHMapIJKLRemoteKey key;
    1384:           PetscHashIter          iter;
    1385:           PetscBool              missing;

    1387:           if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]  Entering face at (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, r, idx));
    1388:           PetscCheck(Np <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %" PetscInt_FMT " cone points", Np);
    1389:           fcp0.rank  = rank;
    1390:           fcp0.index = r;
    1391:           d += Np;
    1392:           /* Find remote face in hash table */
    1393:           key.i = fcp0;
    1394:           key.j = fcone[0];
    1395:           key.k = Np > 2 ? fcone[1] : pmax;
    1396:           key.l = Np > 3 ? fcone[2] : pmax;
    1397:           PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key));
    1398:           if (debug)
    1399:             PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]    key (%d, %" PetscInt_FMT ") (%d, %" PetscInt_FMT ") (%d, %" PetscInt_FMT ") (%d, %" PetscInt_FMT ")\n", rank, (PetscMPIInt)key.i.rank, key.i.index,
    1400:                                               (PetscMPIInt)key.j.rank, key.j.index, (PetscMPIInt)key.k.rank, key.k.index, (PetscMPIInt)key.l.rank, key.l.index));
    1401:           PetscCall(PetscHMapIJKLRemotePut(faceTable, key, &iter, &missing));
    1402:           PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %" PetscInt_FMT " Idx %" PetscInt_FMT " ought to have an associated face", r, idx2);
    1403:           PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx]));
    1404:         }
    1405:       }
    1406:     }
    1407:     if (debug) PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)dm), NULL));
    1408:     PetscCall(PetscHMapIJKLRemoteDestroy(&faceTable));
    1409:   }
    1410:   /* Step 4: Push back owned faces */
    1411:   {
    1412:     PetscSF      sfMulti, sfClaims, sfPointNew;
    1413:     PetscSFNode *remotePointsNew;
    1414:     PetscInt    *remoteOffsets, *localPointsNew;
    1415:     PetscInt     pStart, pEnd, r, NlNew, p;

    1417:     /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
    1418:     PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti));
    1419:     PetscCall(PetscSectionCreate(comm, &claimSection));
    1420:     PetscCall(PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection));
    1421:     PetscCall(PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims));
    1422:     PetscCall(PetscSectionGetStorageSize(claimSection, &claimsSize));
    1423:     PetscCall(PetscMalloc1(claimsSize, &claims));
    1424:     for (p = 0; p < claimsSize; ++p) claims[p].rank = -1;
    1425:     PetscCall(PetscSFBcastBegin(sfClaims, MPIU_SF_NODE, candidatesRemote, claims, MPI_REPLACE));
    1426:     PetscCall(PetscSFBcastEnd(sfClaims, MPIU_SF_NODE, candidatesRemote, claims, MPI_REPLACE));
    1427:     PetscCall(PetscSFDestroy(&sfClaims));
    1428:     PetscCall(PetscFree(remoteOffsets));
    1429:     PetscCall(PetscObjectSetName((PetscObject)claimSection, "Claim Section"));
    1430:     PetscCall(PetscObjectViewFromOptions((PetscObject)claimSection, NULL, "-petscsection_interp_claim_view"));
    1431:     PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims));
    1432:     /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */
    1433:     /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */
    1434:     PetscCall(PetscHMapICreate(&claimshash));
    1435:     for (r = 0; r < Nr; ++r) {
    1436:       PetscInt dof, off, d;

    1438:       if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]  Checking root for claims %" PetscInt_FMT "\n", rank, r));
    1439:       PetscCall(PetscSectionGetDof(candidateSection, r, &dof));
    1440:       PetscCall(PetscSectionGetOffset(candidateSection, r, &off));
    1441:       for (d = 0; d < dof;) {
    1442:         if (claims[off + d].rank >= 0) {
    1443:           const PetscInt  faceInd = off + d;
    1444:           const PetscInt  Np      = candidates[off + d].index;
    1445:           const PetscInt *join    = NULL;
    1446:           PetscInt        joinSize, points[1024], c;

    1448:           if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Found claim for remote point (%d, %" PetscInt_FMT ")\n", rank, (PetscMPIInt)claims[faceInd].rank, claims[faceInd].index));
    1449:           points[0] = r;
    1450:           if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]      point %" PetscInt_FMT "\n", rank, points[0]));
    1451:           for (c = 0, d += 2; c < Np; ++c, ++d) {
    1452:             PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, candidates[off + d], &points[c + 1], NULL));
    1453:             if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]      point %" PetscInt_FMT "\n", rank, points[c + 1]));
    1454:           }
    1455:           PetscCall(DMPlexGetJoin(dm, Np + 1, points, &joinSize, &join));
    1456:           if (joinSize == 1) {
    1457:             if (claims[faceInd].rank == rank) {
    1458:               if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Ignoring local face %" PetscInt_FMT " for non-remote partner\n", rank, join[0]));
    1459:             } else {
    1460:               if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Found local face %" PetscInt_FMT "\n", rank, join[0]));
    1461:               PetscCall(PetscHMapISet(claimshash, join[0], faceInd));
    1462:             }
    1463:           } else {
    1464:             if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    Failed to find face\n", rank));
    1465:           }
    1466:           PetscCall(DMPlexRestoreJoin(dm, Np + 1, points, &joinSize, &join));
    1467:         } else {
    1468:           if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]    No claim for point %" PetscInt_FMT "\n", rank, r));
    1469:           d += claims[off + d].index + 1;
    1470:         }
    1471:       }
    1472:     }
    1473:     if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL));
    1474:     /* Step 6) Create new pointSF from hashed claims */
    1475:     PetscCall(PetscHMapIGetSize(claimshash, &NlNew));
    1476:     PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    1477:     PetscCall(PetscMalloc1(Nl + NlNew, &localPointsNew));
    1478:     PetscCall(PetscMalloc1(Nl + NlNew, &remotePointsNew));
    1479:     for (l = 0; l < Nl; ++l) {
    1480:       localPointsNew[l]        = localPoints[l];
    1481:       remotePointsNew[l].index = remotePoints[l].index;
    1482:       remotePointsNew[l].rank  = remotePoints[l].rank;
    1483:     }
    1484:     p = Nl;
    1485:     PetscCall(PetscHMapIGetKeys(claimshash, &p, localPointsNew));
    1486:     /* We sort new points, and assume they are numbered after all existing points */
    1487:     PetscCall(PetscSortInt(NlNew, PetscSafePointerPlusOffset(localPointsNew, Nl)));
    1488:     for (p = Nl; p < Nl + NlNew; ++p) {
    1489:       PetscInt off;
    1490:       PetscCall(PetscHMapIGet(claimshash, localPointsNew[p], &off));
    1491:       PetscCheck(claims[off].rank >= 0 && claims[off].index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid claim for local point %" PetscInt_FMT ", (%d, %" PetscInt_FMT ")", localPointsNew[p], (PetscMPIInt)claims[off].rank, claims[off].index);
    1492:       remotePointsNew[p] = claims[off];
    1493:     }
    1494:     PetscCall(PetscSFCreate(comm, &sfPointNew));
    1495:     PetscCall(PetscSFSetGraph(sfPointNew, pEnd - pStart, Nl + NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
    1496:     PetscCall(PetscSFSetUp(sfPointNew));
    1497:     PetscCall(DMSetPointSF(dm, sfPointNew));
    1498:     PetscCall(PetscObjectViewFromOptions((PetscObject)sfPointNew, NULL, "-petscsf_interp_view"));
    1499:     if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPointNew, PETSC_FALSE));
    1500:     PetscCall(PetscSFDestroy(&sfPointNew));
    1501:     PetscCall(PetscHMapIDestroy(&claimshash));
    1502:   }
    1503:   PetscCall(PetscHMapIJDestroy(&remoteHash));
    1504:   PetscCall(PetscSectionDestroy(&candidateSection));
    1505:   PetscCall(PetscSectionDestroy(&candidateRemoteSection));
    1506:   PetscCall(PetscSectionDestroy(&claimSection));
    1507:   PetscCall(PetscFree(candidates));
    1508:   PetscCall(PetscFree(candidatesRemote));
    1509:   PetscCall(PetscFree(claims));
    1510:   PetscCall(PetscLogEventEnd(DMPLEX_InterpolateSF, dm, 0, 0, 0));
    1511:   PetscFunctionReturn(PETSC_SUCCESS);
    1512: }

    1514: /*@
    1515:   DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc.

    1517:   Collective

    1519:   Input Parameter:
    1520: . dm - The `DMPLEX` object with only cells and vertices

    1522:   Output Parameter:
    1523: . dmInt - The complete `DMPLEX` object

    1525:   Level: intermediate

    1527:   Note:
    1528:   Labels and coordinates are copied.

    1530:   Developer Notes:
    1531:   It sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`.

    1533: .seealso: `DMPLEX`, `DMPlexUninterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()`
    1534: @*/
    1535: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
    1536: {
    1537:   DMPlexInterpolatedFlag interpolated;
    1538:   DM                     idm, odm = dm;
    1539:   PetscSF                sfPoint;
    1540:   PetscInt               depth, dim, d;
    1541:   const char            *name;
    1542:   PetscBool              flg = PETSC_TRUE;

    1544:   PetscFunctionBegin;
    1546:   PetscAssertPointer(dmInt, 2);
    1547:   PetscCall(PetscLogEventBegin(DMPLEX_Interpolate, dm, 0, 0, 0));
    1548:   PetscCall(DMPlexGetDepth(dm, &depth));
    1549:   PetscCall(DMGetDimension(dm, &dim));
    1550:   PetscCall(DMPlexIsInterpolated(dm, &interpolated));
    1551:   PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
    1552:   if (interpolated == DMPLEX_INTERPOLATED_FULL) {
    1553:     PetscCall(PetscObjectReference((PetscObject)dm));
    1554:     idm = dm;
    1555:   } else {
    1556:     PetscBool nonmanifold = PETSC_FALSE;

    1558:     PetscCall(PetscOptionsGetBool(NULL, dm->hdr.prefix, "-dm_plex_stratify_celltype", &nonmanifold, NULL));
    1559:     if (nonmanifold) {
    1560:       do {
    1561:         const char *prefix;
    1562:         PetscInt    pStart, pEnd, pdepth;
    1563:         PetscBool   done = PETSC_TRUE;

    1565:         // Find a point which is not correctly interpolated
    1566:         PetscCall(DMPlexGetChart(odm, &pStart, &pEnd));
    1567:         for (PetscInt p = pStart; p < pEnd; ++p) {
    1568:           DMPolytopeType  ct;
    1569:           const PetscInt *cone;
    1570:           PetscInt        coneSize, cdepth;

    1572:           PetscCall(DMPlexGetPointDepth(odm, p, &pdepth));
    1573:           PetscCall(DMPlexGetCellType(odm, p, &ct));
    1574:           // Check against celltype
    1575:           if (pdepth != DMPolytopeTypeGetDim(ct)) {
    1576:             done = PETSC_FALSE;
    1577:             break;
    1578:           }
    1579:           // Check against boundary
    1580:           PetscCall(DMPlexGetCone(odm, p, &cone));
    1581:           PetscCall(DMPlexGetConeSize(odm, p, &coneSize));
    1582:           for (PetscInt c = 0; c < coneSize; ++c) {
    1583:             PetscCall(DMPlexGetPointDepth(odm, cone[c], &cdepth));
    1584:             if (cdepth != pdepth - 1) {
    1585:               done = PETSC_FALSE;
    1586:               p    = pEnd;
    1587:               break;
    1588:             }
    1589:           }
    1590:         }
    1591:         if (done) break;
    1592:         /* Create interpolated mesh */
    1593:         PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &idm));
    1594:         PetscCall(DMSetType(idm, DMPLEX));
    1595:         PetscCall(DMSetDimension(idm, dim));
    1596:         PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
    1597:         PetscCall(PetscObjectSetOptionsPrefix((PetscObject)idm, prefix));
    1598:         if (depth > 0) {
    1599:           PetscCall(DMPlexInterpolateFaces_Internal(odm, pdepth, idm));
    1600:           PetscCall(DMGetPointSF(odm, &sfPoint));
    1601:           if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(odm, sfPoint, PETSC_FALSE));
    1602:           {
    1603:             /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */
    1604:             PetscInt nroots;
    1605:             PetscCall(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL));
    1606:             if (nroots >= 0) PetscCall(DMPlexInterpolatePointSF(idm, sfPoint));
    1607:           }
    1608:         }
    1609:         if (odm != dm) PetscCall(DMDestroy(&odm));
    1610:         odm = idm;
    1611:       } while (1);
    1612:     } else {
    1613:       for (d = 1; d < dim; ++d) {
    1614:         const char *prefix;

    1616:         /* Create interpolated mesh */
    1617:         PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &idm));
    1618:         PetscCall(DMSetType(idm, DMPLEX));
    1619:         PetscCall(DMSetDimension(idm, dim));
    1620:         PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
    1621:         PetscCall(PetscObjectSetOptionsPrefix((PetscObject)idm, prefix));
    1622:         if (depth > 0) {
    1623:           PetscCall(DMPlexInterpolateFaces_Internal(odm, 1, idm));
    1624:           PetscCall(DMGetPointSF(odm, &sfPoint));
    1625:           if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(odm, sfPoint, PETSC_FALSE));
    1626:           {
    1627:             /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */
    1628:             PetscInt nroots;
    1629:             PetscCall(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL));
    1630:             if (nroots >= 0) PetscCall(DMPlexInterpolatePointSF(idm, sfPoint));
    1631:           }
    1632:         }
    1633:         if (odm != dm) PetscCall(DMDestroy(&odm));
    1634:         odm = idm;
    1635:       }
    1636:     }
    1637:     PetscCall(PetscObjectGetName((PetscObject)dm, &name));
    1638:     PetscCall(PetscObjectSetName((PetscObject)idm, name));
    1639:     PetscCall(DMPlexCopyCoordinates(dm, idm));
    1640:     PetscCall(DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL));
    1641:     PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL));
    1642:     if (flg) PetscCall(DMPlexOrientInterface_Internal(idm));
    1643:   }
    1644:   /* This function makes the mesh fully interpolated on all ranks */
    1645:   {
    1646:     DM_Plex *plex      = (DM_Plex *)idm->data;
    1647:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL;
    1648:   }
    1649:   PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, idm));
    1650:   *dmInt = idm;
    1651:   PetscCall(PetscLogEventEnd(DMPLEX_Interpolate, dm, 0, 0, 0));
    1652:   PetscFunctionReturn(PETSC_SUCCESS);
    1653: }

    1655: /*@
    1656:   DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices

    1658:   Collective

    1660:   Input Parameter:
    1661: . dmA - The `DMPLEX` object with initial coordinates

    1663:   Output Parameter:
    1664: . dmB - The `DMPLEX` object with copied coordinates

    1666:   Level: intermediate

    1668:   Notes:
    1669:   This is typically used when adding pieces other than vertices to a mesh

    1671:   This function does not copy localized coordinates.

    1673: .seealso: `DMPLEX`, `DMCopyLabels()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMGetCoordinateSection()`
    1674: @*/
    1675: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
    1676: {
    1677:   Vec          coordinatesA, coordinatesB;
    1678:   VecType      vtype;
    1679:   PetscSection coordSectionA, coordSectionB;
    1680:   PetscScalar *coordsA, *coordsB;
    1681:   PetscInt     spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
    1682:   PetscInt     cStartA, cEndA, cStartB, cEndB, cS, cE, cdim;

    1684:   PetscFunctionBegin;
    1687:   if (dmA == dmB) PetscFunctionReturn(PETSC_SUCCESS);
    1688:   PetscCall(DMGetCoordinateDim(dmA, &cdim));
    1689:   PetscCall(DMSetCoordinateDim(dmB, cdim));
    1690:   PetscCall(DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA));
    1691:   PetscCall(DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB));
    1692:   PetscCheck((vEndA - vStartA) == (vEndB - vStartB), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of vertices in first DM %" PetscInt_FMT " != %" PetscInt_FMT " in the second DM", vEndA - vStartA, vEndB - vStartB);
    1693:   /* Copy over discretization if it exists */
    1694:   {
    1695:     DM                 cdmA, cdmB;
    1696:     PetscDS            dsA, dsB;
    1697:     PetscObject        objA, objB;
    1698:     PetscClassId       idA, idB;
    1699:     const PetscScalar *constants;
    1700:     PetscInt           cdim, Nc;

    1702:     PetscCall(DMGetCoordinateDM(dmA, &cdmA));
    1703:     PetscCall(DMGetCoordinateDM(dmB, &cdmB));
    1704:     PetscCall(DMGetField(cdmA, 0, NULL, &objA));
    1705:     PetscCall(DMGetField(cdmB, 0, NULL, &objB));
    1706:     PetscCall(PetscObjectGetClassId(objA, &idA));
    1707:     PetscCall(PetscObjectGetClassId(objB, &idB));
    1708:     if ((idA == PETSCFE_CLASSID) && (idA != idB)) {
    1709:       PetscCall(DMSetField(cdmB, 0, NULL, objA));
    1710:       PetscCall(DMCreateDS(cdmB));
    1711:       PetscCall(DMGetDS(cdmA, &dsA));
    1712:       PetscCall(DMGetDS(cdmB, &dsB));
    1713:       PetscCall(PetscDSGetCoordinateDimension(dsA, &cdim));
    1714:       PetscCall(PetscDSSetCoordinateDimension(dsB, cdim));
    1715:       PetscCall(PetscDSGetConstants(dsA, &Nc, &constants));
    1716:       PetscCall(PetscDSSetConstants(dsB, Nc, (PetscScalar *)constants));
    1717:     }
    1718:   }
    1719:   PetscCall(DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA));
    1720:   PetscCall(DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB));
    1721:   PetscCall(DMGetCoordinateSection(dmA, &coordSectionA));
    1722:   PetscCall(DMGetCoordinateSection(dmB, &coordSectionB));
    1723:   if (coordSectionA == coordSectionB) PetscFunctionReturn(PETSC_SUCCESS);
    1724:   PetscCall(PetscSectionGetNumFields(coordSectionA, &Nf));
    1725:   if (!Nf) PetscFunctionReturn(PETSC_SUCCESS);
    1726:   PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %" PetscInt_FMT, Nf);
    1727:   if (!coordSectionB) {
    1728:     PetscInt dim;

    1730:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coordSectionA), &coordSectionB));
    1731:     PetscCall(DMGetCoordinateDim(dmA, &dim));
    1732:     PetscCall(DMSetCoordinateSection(dmB, dim, coordSectionB));
    1733:     PetscCall(PetscObjectDereference((PetscObject)coordSectionB));
    1734:   }
    1735:   PetscCall(PetscSectionSetNumFields(coordSectionB, 1));
    1736:   PetscCall(PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim));
    1737:   PetscCall(PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim));
    1738:   PetscCall(PetscSectionGetChart(coordSectionA, &cS, &cE));
    1739:   cS = vStartB;
    1740:   cE = vEndB;
    1741:   PetscCall(PetscSectionSetChart(coordSectionB, cS, cE));
    1742:   for (v = vStartB; v < vEndB; ++v) {
    1743:     PetscCall(PetscSectionSetDof(coordSectionB, v, spaceDim));
    1744:     PetscCall(PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim));
    1745:   }
    1746:   PetscCall(PetscSectionSetUp(coordSectionB));
    1747:   PetscCall(PetscSectionGetStorageSize(coordSectionB, &coordSizeB));
    1748:   PetscCall(DMGetCoordinatesLocal(dmA, &coordinatesA));
    1749:   PetscCall(VecCreate(PETSC_COMM_SELF, &coordinatesB));
    1750:   PetscCall(PetscObjectSetName((PetscObject)coordinatesB, "coordinates"));
    1751:   PetscCall(VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE));
    1752:   PetscCall(VecGetBlockSize(coordinatesA, &d));
    1753:   PetscCall(VecSetBlockSize(coordinatesB, d));
    1754:   PetscCall(VecGetType(coordinatesA, &vtype));
    1755:   PetscCall(VecSetType(coordinatesB, vtype));
    1756:   PetscCall(VecGetArray(coordinatesA, &coordsA));
    1757:   PetscCall(VecGetArray(coordinatesB, &coordsB));
    1758:   for (v = 0; v < vEndB - vStartB; ++v) {
    1759:     PetscInt offA, offB;

    1761:     PetscCall(PetscSectionGetOffset(coordSectionA, v + vStartA, &offA));
    1762:     PetscCall(PetscSectionGetOffset(coordSectionB, v + vStartB, &offB));
    1763:     for (d = 0; d < spaceDim; ++d) coordsB[offB + d] = coordsA[offA + d];
    1764:   }
    1765:   PetscCall(VecRestoreArray(coordinatesA, &coordsA));
    1766:   PetscCall(VecRestoreArray(coordinatesB, &coordsB));
    1767:   PetscCall(DMSetCoordinatesLocal(dmB, coordinatesB));
    1768:   PetscCall(VecDestroy(&coordinatesB));
    1769:   PetscFunctionReturn(PETSC_SUCCESS);
    1770: }

    1772: /*@
    1773:   DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh

    1775:   Collective

    1777:   Input Parameter:
    1778: . dm - The complete `DMPLEX` object

    1780:   Output Parameter:
    1781: . dmUnint - The `DMPLEX` object with only cells and vertices

    1783:   Level: intermediate

    1785:   Note:
    1786:   It does not copy over the coordinates.

    1788:   Developer Notes:
    1789:   Sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`.

    1791: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()`
    1792: @*/
    1793: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
    1794: {
    1795:   DMPlexInterpolatedFlag interpolated;
    1796:   DM                     udm;
    1797:   PetscInt               dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;

    1799:   PetscFunctionBegin;
    1801:   PetscAssertPointer(dmUnint, 2);
    1802:   PetscCall(PetscLogEventBegin(DMPLEX_Uninterpolate, dm, 0, 0, 0));
    1803:   PetscCall(DMGetDimension(dm, &dim));
    1804:   PetscCall(DMPlexIsInterpolated(dm, &interpolated));
    1805:   PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
    1806:   if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) {
    1807:     /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */
    1808:     PetscCall(PetscObjectReference((PetscObject)dm));
    1809:     *dmUnint = dm;
    1810:     PetscFunctionReturn(PETSC_SUCCESS);
    1811:   }
    1812:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
    1813:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
    1814:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &udm));
    1815:   PetscCall(DMSetType(udm, DMPLEX));
    1816:   PetscCall(DMSetDimension(udm, dim));
    1817:   PetscCall(DMPlexSetChart(udm, cStart, vEnd));
    1818:   for (c = cStart; c < cEnd; ++c) {
    1819:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

    1821:     PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
    1822:     for (cl = 0; cl < closureSize * 2; cl += 2) {
    1823:       const PetscInt p = closure[cl];

    1825:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
    1826:     }
    1827:     PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
    1828:     PetscCall(DMPlexSetConeSize(udm, c, coneSize));
    1829:     maxConeSize = PetscMax(maxConeSize, coneSize);
    1830:   }
    1831:   PetscCall(DMSetUp(udm));
    1832:   PetscCall(PetscMalloc1(maxConeSize, &cone));
    1833:   for (c = cStart; c < cEnd; ++c) {
    1834:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

    1836:     PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
    1837:     for (cl = 0; cl < closureSize * 2; cl += 2) {
    1838:       const PetscInt p = closure[cl];

    1840:       if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
    1841:     }
    1842:     PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
    1843:     PetscCall(DMPlexSetCone(udm, c, cone));
    1844:   }
    1845:   PetscCall(PetscFree(cone));
    1846:   PetscCall(DMPlexSymmetrize(udm));
    1847:   PetscCall(DMPlexStratify(udm));
    1848:   /* Reduce SF */
    1849:   {
    1850:     PetscSF            sfPoint, sfPointUn;
    1851:     const PetscSFNode *remotePoints;
    1852:     const PetscInt    *localPoints;
    1853:     PetscSFNode       *remotePointsUn;
    1854:     PetscInt          *localPointsUn;
    1855:     PetscInt           numRoots, numLeaves, l;
    1856:     PetscInt           numLeavesUn = 0, n = 0;

    1858:     /* Get original SF information */
    1859:     PetscCall(DMGetPointSF(dm, &sfPoint));
    1860:     if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPoint, PETSC_FALSE));
    1861:     PetscCall(DMGetPointSF(udm, &sfPointUn));
    1862:     PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints));
    1863:     if (numRoots >= 0) {
    1864:       /* Allocate space for cells and vertices */
    1865:       for (l = 0; l < numLeaves; ++l) {
    1866:         const PetscInt p = localPoints[l];

    1868:         if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) numLeavesUn++;
    1869:       }
    1870:       /* Fill in leaves */
    1871:       PetscCall(PetscMalloc1(numLeavesUn, &remotePointsUn));
    1872:       PetscCall(PetscMalloc1(numLeavesUn, &localPointsUn));
    1873:       for (l = 0; l < numLeaves; l++) {
    1874:         const PetscInt p = localPoints[l];

    1876:         if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) {
    1877:           localPointsUn[n]        = p;
    1878:           remotePointsUn[n].rank  = remotePoints[l].rank;
    1879:           remotePointsUn[n].index = remotePoints[l].index;
    1880:           ++n;
    1881:         }
    1882:       }
    1883:       PetscCheck(n == numLeavesUn, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %" PetscInt_FMT " != %" PetscInt_FMT, n, numLeavesUn);
    1884:       PetscCall(PetscSFSetGraph(sfPointUn, cEnd - cStart + vEnd - vStart, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER));
    1885:     }
    1886:   }
    1887:   /* This function makes the mesh fully uninterpolated on all ranks */
    1888:   {
    1889:     DM_Plex *plex      = (DM_Plex *)udm->data;
    1890:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE;
    1891:   }
    1892:   PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, udm));
    1893:   if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(udm, NULL, PETSC_FALSE));
    1894:   *dmUnint = udm;
    1895:   PetscCall(PetscLogEventEnd(DMPLEX_Uninterpolate, dm, 0, 0, 0));
    1896:   PetscFunctionReturn(PETSC_SUCCESS);
    1897: }

    1899: static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated)
    1900: {
    1901:   PetscInt coneSize, depth, dim, h, p, pStart, pEnd;
    1902:   MPI_Comm comm;

    1904:   PetscFunctionBegin;
    1905:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    1906:   PetscCall(DMPlexGetDepth(dm, &depth));
    1907:   PetscCall(DMGetDimension(dm, &dim));

    1909:   if (depth == dim) {
    1910:     *interpolated = DMPLEX_INTERPOLATED_FULL;
    1911:     if (!dim) goto finish;

    1913:     /* Check points at height = dim are vertices (have no cones) */
    1914:     PetscCall(DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd));
    1915:     for (p = pStart; p < pEnd; p++) {
    1916:       PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
    1917:       if (coneSize) {
    1918:         *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
    1919:         goto finish;
    1920:       }
    1921:     }

    1923:     /* Check points at height < dim have cones */
    1924:     for (h = 0; h < dim; h++) {
    1925:       PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd));
    1926:       for (p = pStart; p < pEnd; p++) {
    1927:         PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
    1928:         if (!coneSize) {
    1929:           *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
    1930:           goto finish;
    1931:         }
    1932:       }
    1933:     }
    1934:   } else if (depth == 1) {
    1935:     *interpolated = DMPLEX_INTERPOLATED_NONE;
    1936:   } else {
    1937:     *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
    1938:   }
    1939: finish:
    1940:   PetscFunctionReturn(PETSC_SUCCESS);
    1941: }

    1943: /*@
    1944:   DMPlexIsInterpolated - Find out to what extent the `DMPLEX` is topologically interpolated.

    1946:   Not Collective

    1948:   Input Parameter:
    1949: . dm - The `DMPLEX` object

    1951:   Output Parameter:
    1952: . interpolated - Flag whether the `DM` is interpolated

    1954:   Level: intermediate

    1956:   Notes:
    1957:   Unlike `DMPlexIsInterpolatedCollective()`, this is NOT collective
    1958:   so the results can be different on different ranks in special cases.
    1959:   However, `DMPlexInterpolate()` guarantees the result is the same on all.

    1961:   Unlike `DMPlexIsInterpolatedCollective()`, this cannot return `DMPLEX_INTERPOLATED_MIXED`.

    1963:   Developer Notes:
    1964:   Initially, plex->interpolated = `DMPLEX_INTERPOLATED_INVALID`.

    1966:   If plex->interpolated == `DMPLEX_INTERPOLATED_INVALID`, `DMPlexIsInterpolated_Internal()` is called.
    1967:   It checks the actual topology and sets plex->interpolated on each rank separately to one of
    1968:   `DMPLEX_INTERPOLATED_NONE`, `DMPLEX_INTERPOLATED_PARTIAL` or `DMPLEX_INTERPOLATED_FULL`.

    1970:   If plex->interpolated != `DMPLEX_INTERPOLATED_INVALID`, this function just returns plex->interpolated.

    1972:   `DMPlexInterpolate()` sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`,
    1973:   and DMPlexUninterpolate() sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`.

    1975: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolatedCollective()`
    1976: @*/
    1977: PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated)
    1978: {
    1979:   DM_Plex *plex = (DM_Plex *)dm->data;

    1981:   PetscFunctionBegin;
    1983:   PetscAssertPointer(interpolated, 2);
    1984:   if (plex->interpolated < 0) {
    1985:     PetscCall(DMPlexIsInterpolated_Internal(dm, &plex->interpolated));
    1986:   } else if (PetscDefined(USE_DEBUG)) {
    1987:     DMPlexInterpolatedFlag flg;

    1989:     PetscCall(DMPlexIsInterpolated_Internal(dm, &flg));
    1990:     PetscCheck(plex->tr || flg == plex->interpolated, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stashed DMPlexInterpolatedFlag %s is inconsistent with current %s", DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[flg]);
    1991:   }
    1992:   *interpolated = plex->interpolated;
    1993:   PetscFunctionReturn(PETSC_SUCCESS);
    1994: }

    1996: /*@
    1997:   DMPlexIsInterpolatedCollective - Find out to what extent the `DMPLEX` is topologically interpolated (in collective manner).

    1999:   Collective

    2001:   Input Parameter:
    2002: . dm - The `DMPLEX` object

    2004:   Output Parameter:
    2005: . interpolated - Flag whether the `DM` is interpolated

    2007:   Level: intermediate

    2009:   Notes:
    2010:   Unlike `DMPlexIsInterpolated()`, this is collective so the results are guaranteed to be the same on all ranks.

    2012:   This function will return `DMPLEX_INTERPOLATED_MIXED` if the results of `DMPlexIsInterpolated()` are different on different ranks.

    2014:   Developer Notes:
    2015:   Initially, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_INVALID`.

    2017:   If plex->interpolatedCollective == `DMPLEX_INTERPOLATED_INVALID`, this function calls `DMPlexIsInterpolated()` which sets plex->interpolated.
    2018:   `MPI_Allreduce()` is then called and collectively consistent flag plex->interpolatedCollective is set and returned;
    2019:   if plex->interpolated varies on different ranks, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_MIXED`,
    2020:   otherwise sets plex->interpolatedCollective = plex->interpolated.

    2022:   If plex->interpolatedCollective != `DMPLEX_INTERPOLATED_INVALID`, this function just returns plex->interpolatedCollective.

    2024: .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolated()`
    2025: @*/
    2026: PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated)
    2027: {
    2028:   DM_Plex  *plex  = (DM_Plex *)dm->data;
    2029:   PetscBool debug = PETSC_FALSE;

    2031:   PetscFunctionBegin;
    2033:   PetscAssertPointer(interpolated, 2);
    2034:   PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL));
    2035:   if (plex->interpolatedCollective < 0) {
    2036:     DMPlexInterpolatedFlag min, max;
    2037:     MPI_Comm               comm;

    2039:     PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    2040:     PetscCall(DMPlexIsInterpolated(dm, &plex->interpolatedCollective));
    2041:     PetscCallMPI(MPIU_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm));
    2042:     PetscCallMPI(MPIU_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm));
    2043:     if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED;
    2044:     if (debug) {
    2045:       PetscMPIInt rank;

    2047:       PetscCallMPI(MPI_Comm_rank(comm, &rank));
    2048:       PetscCall(PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]));
    2049:       PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
    2050:     }
    2051:   }
    2052:   *interpolated = plex->interpolatedCollective;
    2053:   PetscFunctionReturn(PETSC_SUCCESS);
    2054: }