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

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

     11: typedef struct _PetscHashIJKLKey { PetscInt i, j, k, l; } PetscHashIJKLKey;

     13: #define PetscHashIJKLKeyHash(key) \
     14:   PetscHashCombine(PetscHashCombine(PetscHashInt((key).i),PetscHashInt((key).j)), \
     15:                    PetscHashCombine(PetscHashInt((key).k),PetscHashInt((key).l)))

     17: #define PetscHashIJKLKeyEqual(k1,k2) \
     18:   (((k1).i==(k2).i) ? ((k1).j==(k2).j) ? ((k1).k==(k2).k) ? ((k1).l==(k2).l) : 0 : 0 : 0)

     20: PETSC_HASH_MAP(HashIJKL, PetscHashIJKLKey, PetscInt, PetscHashIJKLKeyHash, PetscHashIJKLKeyEqual, -1)

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

     24: typedef struct _PetscHashIJKLRemoteKey { PetscSFNode i, j, k, l; } PetscHashIJKLRemoteKey;

     26: #define PetscHashIJKLRemoteKeyHash(key) \
     27:   PetscHashCombine(PetscHashCombine(PetscHashInt((key).i.rank + (key).i.index),PetscHashInt((key).j.rank + (key).j.index)), \
     28:                    PetscHashCombine(PetscHashInt((key).k.rank + (key).k.index),PetscHashInt((key).l.rank + (key).l.index)))

     30: #define PetscHashIJKLRemoteKeyEqual(k1,k2) \
     31:   (((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)

     33: PETSC_HASH_MAP(HashIJKLRemote, PetscHashIJKLRemoteKey, PetscSFNode, PetscHashIJKLRemoteKeyHash, PetscHashIJKLRemoteKeyEqual, _PetscInvalidSFNode)

     35: static PetscErrorCode PetscSortSFNode(PetscInt n, PetscSFNode A[])
     36: {
     37:   PetscInt i;

     40:   for (i = 1; i < n; ++i) {
     41:     PetscSFNode x = A[i];
     42:     PetscInt    j;

     44:     for (j = i-1; j >= 0; --j) {
     45:       if ((A[j].rank > x.rank) || (A[j].rank == x.rank && A[j].index > x.index)) break;
     46:       A[j+1] = A[j];
     47:     }
     48:     A[j+1] = x;
     49:   }
     50:   return(0);
     51: }

     53: /*
     54:   DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone
     55: */
     56: PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[])
     57: {
     58:   DMPolytopeType *typesTmp;
     59:   PetscInt       *sizesTmp, *facesTmp;
     60:   PetscInt        maxConeSize, maxSupportSize;
     61:   PetscErrorCode  ierr;

     66:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
     67:   if (faceTypes) {DMGetWorkArray(dm, PetscMax(maxConeSize, maxSupportSize),           MPIU_INT, &typesTmp);}
     68:   if (faceSizes) {DMGetWorkArray(dm, PetscMax(maxConeSize, maxSupportSize),           MPIU_INT, &sizesTmp);}
     69:   if (faces)     {DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), MPIU_INT, &facesTmp);}
     70:   switch (ct) {
     71:     case DM_POLYTOPE_POINT:
     72:       if (numFaces) *numFaces = 0;
     73:       break;
     74:     case DM_POLYTOPE_SEGMENT:
     75:       if (numFaces) *numFaces = 2;
     76:       if (faceTypes) {
     77:         typesTmp[0] = DM_POLYTOPE_POINT; typesTmp[1] = DM_POLYTOPE_POINT;
     78:         *faceTypes = typesTmp;
     79:       }
     80:       if (faceSizes) {
     81:         sizesTmp[0] = 1; sizesTmp[1] = 1;
     82:         *faceSizes = sizesTmp;
     83:       }
     84:       if (faces) {
     85:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
     86:         *faces = facesTmp;
     87:       }
     88:       break;
     89:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
     90:       if (numFaces) *numFaces = 2;
     91:       if (faceTypes) {
     92:         typesTmp[0] = DM_POLYTOPE_POINT; typesTmp[1] = DM_POLYTOPE_POINT;
     93:         *faceTypes = typesTmp;
     94:       }
     95:       if (faceSizes) {
     96:         sizesTmp[0] = 1; sizesTmp[1] = 1;
     97:         *faceSizes = sizesTmp;
     98:       }
     99:       if (faces) {
    100:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
    101:         *faces = facesTmp;
    102:       }
    103:       break;
    104:     case DM_POLYTOPE_TRIANGLE:
    105:       if (numFaces) *numFaces = 3;
    106:       if (faceTypes) {
    107:         typesTmp[0] = DM_POLYTOPE_SEGMENT; typesTmp[1] = DM_POLYTOPE_SEGMENT; typesTmp[2] = DM_POLYTOPE_SEGMENT;
    108:         *faceTypes = typesTmp;
    109:       }
    110:       if (faceSizes) {
    111:         sizesTmp[0] = 2; sizesTmp[1] = 2; sizesTmp[2] = 2;
    112:         *faceSizes = sizesTmp;
    113:       }
    114:       if (faces) {
    115:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
    116:         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
    117:         facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
    118:         *faces = facesTmp;
    119:       }
    120:       break;
    121:     case DM_POLYTOPE_QUADRILATERAL:
    122:       /* Vertices follow right hand rule */
    123:       if (numFaces) *numFaces = 4;
    124:       if (faceTypes) {
    125:         typesTmp[0] = DM_POLYTOPE_SEGMENT; typesTmp[1] = DM_POLYTOPE_SEGMENT; typesTmp[2] = DM_POLYTOPE_SEGMENT; typesTmp[3] = DM_POLYTOPE_SEGMENT;
    126:         *faceTypes = typesTmp;
    127:       }
    128:       if (faceSizes) {
    129:         sizesTmp[0] = 2; sizesTmp[1] = 2; sizesTmp[2] = 2; sizesTmp[3] = 2;
    130:         *faceSizes = sizesTmp;
    131:       }
    132:       if (faces) {
    133:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
    134:         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
    135:         facesTmp[4] = cone[2]; facesTmp[5] = cone[3];
    136:         facesTmp[6] = cone[3]; facesTmp[7] = cone[0];
    137:         *faces = facesTmp;
    138:       }
    139:       break;
    140:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
    141:       if (numFaces) *numFaces = 4;
    142:       if (faceTypes) {
    143:         typesTmp[0] = DM_POLYTOPE_SEGMENT; typesTmp[1] = DM_POLYTOPE_SEGMENT; typesTmp[2] = DM_POLYTOPE_POINT_PRISM_TENSOR; typesTmp[3] = DM_POLYTOPE_POINT_PRISM_TENSOR;
    144:         *faceTypes = typesTmp;
    145:       }
    146:       if (faceSizes) {
    147:         sizesTmp[0] = 2; sizesTmp[1] = 2; sizesTmp[2] = 2; sizesTmp[3] = 2;
    148:         *faceSizes = sizesTmp;
    149:       }
    150:       if (faces) {
    151:         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
    152:         facesTmp[2] = cone[2]; facesTmp[3] = cone[3];
    153:         facesTmp[4] = cone[0]; facesTmp[5] = cone[2];
    154:         facesTmp[6] = cone[1]; facesTmp[7] = cone[3];
    155:         *faces = facesTmp;
    156:       }
    157:       break;
    158:     case DM_POLYTOPE_TETRAHEDRON:
    159:       /* Vertices of first face follow right hand rule and normal points away from last vertex */
    160:       if (numFaces) *numFaces = 4;
    161:       if (faceTypes) {
    162:         typesTmp[0] = DM_POLYTOPE_TRIANGLE; typesTmp[1] = DM_POLYTOPE_TRIANGLE; typesTmp[2] = DM_POLYTOPE_TRIANGLE; typesTmp[3] = DM_POLYTOPE_TRIANGLE;
    163:         *faceTypes = typesTmp;
    164:       }
    165:       if (faceSizes) {
    166:         sizesTmp[0] = 3; sizesTmp[1] = 3; sizesTmp[2] = 3; sizesTmp[3] = 3;
    167:         *faceSizes = sizesTmp;
    168:       }
    169:       if (faces) {
    170:         facesTmp[0] = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2];
    171:         facesTmp[3] = cone[0]; facesTmp[4]  = cone[3]; facesTmp[5]  = cone[1];
    172:         facesTmp[6] = cone[0]; facesTmp[7]  = cone[2]; facesTmp[8]  = cone[3];
    173:         facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3];
    174:         *faces = facesTmp;
    175:       }
    176:       break;
    177:     case DM_POLYTOPE_HEXAHEDRON:
    178:       /*  7--------6
    179:          /|       /|
    180:         / |      / |
    181:        4--------5  |
    182:        |  |     |  |
    183:        |  |     |  |
    184:        |  1--------2
    185:        | /      | /
    186:        |/       |/
    187:        0--------3
    188:        */
    189:       if (numFaces) *numFaces = 6;
    190:       if (faceTypes) {
    191:         typesTmp[0] = DM_POLYTOPE_QUADRILATERAL; typesTmp[1] = DM_POLYTOPE_QUADRILATERAL; typesTmp[2] = DM_POLYTOPE_QUADRILATERAL;
    192:         typesTmp[3] = DM_POLYTOPE_QUADRILATERAL; typesTmp[4] = DM_POLYTOPE_QUADRILATERAL; typesTmp[5] = DM_POLYTOPE_QUADRILATERAL;
    193:         *faceTypes = typesTmp;
    194:       }
    195:       if (faceSizes) {
    196:         sizesTmp[0] = 4; sizesTmp[1] = 4; sizesTmp[2] = 4; sizesTmp[3] = 4; sizesTmp[4] = 4; sizesTmp[5] = 4;
    197:         *faceSizes = sizesTmp;
    198:       }
    199:       if (faces) {
    200:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2]; facesTmp[3]  = cone[3]; /* Bottom */
    201:         facesTmp[4]  = cone[4]; facesTmp[5]  = cone[5]; facesTmp[6]  = cone[6]; facesTmp[7]  = cone[7]; /* Top */
    202:         facesTmp[8]  = cone[0]; facesTmp[9]  = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4]; /* Front */
    203:         facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6]; /* Back */
    204:         facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5]; /* Right */
    205:         facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1]; /* Left */
    206:         *faces = facesTmp;
    207:       }
    208:       break;
    209:     case DM_POLYTOPE_TRI_PRISM:
    210:       if (numFaces) *numFaces = 5;
    211:       if (faceTypes) {
    212:         typesTmp[0] = DM_POLYTOPE_TRIANGLE; typesTmp[1] = DM_POLYTOPE_TRIANGLE;
    213:         typesTmp[2] = DM_POLYTOPE_QUADRILATERAL; typesTmp[3] = DM_POLYTOPE_QUADRILATERAL; typesTmp[4] = DM_POLYTOPE_QUADRILATERAL;
    214:         *faceTypes = typesTmp;
    215:       }
    216:       if (faceSizes) {
    217:         sizesTmp[0] = 3; sizesTmp[1] = 3;
    218:         sizesTmp[2] = 4; sizesTmp[3] = 4; sizesTmp[4] = 4;
    219:         *faceSizes = sizesTmp;
    220:       }
    221:       if (faces) {
    222:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2];                         /* Bottom */
    223:         facesTmp[3]  = cone[3]; facesTmp[4]  = cone[4]; facesTmp[5]  = cone[5];                         /* Top */
    224:         facesTmp[6]  = cone[0]; facesTmp[7]  = cone[2]; facesTmp[8]  = cone[4]; facesTmp[9]  = cone[3]; /* Back left */
    225:         facesTmp[10] = cone[2]; facesTmp[11] = cone[1]; facesTmp[12] = cone[5]; facesTmp[13] = cone[4]; /* Front */
    226:         facesTmp[14] = cone[1]; facesTmp[15] = cone[0]; facesTmp[16] = cone[3]; facesTmp[17] = cone[5]; /* Back right */
    227:         *faces = facesTmp;
    228:       }
    229:       break;
    230:     case DM_POLYTOPE_TRI_PRISM_TENSOR:
    231:       if (numFaces)     *numFaces = 5;
    232:       if (faceTypes) {
    233:         typesTmp[0] = DM_POLYTOPE_TRIANGLE; typesTmp[1] = DM_POLYTOPE_TRIANGLE;
    234:         typesTmp[2] = DM_POLYTOPE_SEG_PRISM_TENSOR; typesTmp[3] = DM_POLYTOPE_SEG_PRISM_TENSOR; typesTmp[4] = DM_POLYTOPE_SEG_PRISM_TENSOR;
    235:         *faceTypes = typesTmp;
    236:       }
    237:       if (faceSizes) {
    238:         sizesTmp[0] = 3; sizesTmp[1] = 3;
    239:         sizesTmp[2] = 4; sizesTmp[3] = 4; sizesTmp[4] = 4;
    240:         *faceSizes = sizesTmp;
    241:       }
    242:       if (faces) {
    243:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2];                         /* Bottom */
    244:         facesTmp[3]  = cone[3]; facesTmp[4]  = cone[4]; facesTmp[5]  = cone[5];                         /* Top */
    245:         facesTmp[6]  = cone[0]; facesTmp[7]  = cone[1]; facesTmp[8]  = cone[3]; facesTmp[9]  = cone[4]; /* Back left */
    246:         facesTmp[10] = cone[1]; facesTmp[11] = cone[2]; facesTmp[12] = cone[4]; facesTmp[13] = cone[5]; /* Back right */
    247:         facesTmp[14] = cone[2]; facesTmp[15] = cone[0]; facesTmp[16] = cone[5]; facesTmp[17] = cone[3]; /* Front */
    248:         *faces = facesTmp;
    249:       }
    250:       break;
    251:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
    252:       /*  7--------6
    253:          /|       /|
    254:         / |      / |
    255:        4--------5  |
    256:        |  |     |  |
    257:        |  |     |  |
    258:        |  3--------2
    259:        | /      | /
    260:        |/       |/
    261:        0--------1
    262:        */
    263:       if (numFaces) *numFaces = 6;
    264:       if (faceTypes) {
    265:         typesTmp[0] = DM_POLYTOPE_QUADRILATERAL;    typesTmp[1] = DM_POLYTOPE_QUADRILATERAL;    typesTmp[2] = DM_POLYTOPE_SEG_PRISM_TENSOR;
    266:         typesTmp[3] = DM_POLYTOPE_SEG_PRISM_TENSOR; typesTmp[4] = DM_POLYTOPE_SEG_PRISM_TENSOR; typesTmp[5] = DM_POLYTOPE_SEG_PRISM_TENSOR;
    267:         *faceTypes = typesTmp;
    268:       }
    269:       if (faceSizes) {
    270:         sizesTmp[0] = 4; sizesTmp[1] = 4; sizesTmp[2] = 4; sizesTmp[3] = 4; sizesTmp[4] = 4; sizesTmp[5] = 4;
    271:         *faceSizes = sizesTmp;
    272:       }
    273:       if (faces) {
    274:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2]; facesTmp[3]  = cone[3]; /* Bottom */
    275:         facesTmp[4]  = cone[4]; facesTmp[5]  = cone[5]; facesTmp[6]  = cone[6]; facesTmp[7]  = cone[7]; /* Top */
    276:         facesTmp[8]  = cone[0]; facesTmp[9]  = cone[1]; facesTmp[10] = cone[4]; facesTmp[11] = cone[5]; /* Front */
    277:         facesTmp[12] = cone[1]; facesTmp[13] = cone[2]; facesTmp[14] = cone[5]; facesTmp[15] = cone[6]; /* Right */
    278:         facesTmp[16] = cone[2]; facesTmp[17] = cone[3]; facesTmp[18] = cone[6]; facesTmp[19] = cone[7]; /* Back */
    279:         facesTmp[20] = cone[3]; facesTmp[21] = cone[0]; facesTmp[22] = cone[7]; facesTmp[23] = cone[4]; /* Left */
    280:         *faces = facesTmp;
    281:       }
    282:       break;
    283:     case DM_POLYTOPE_PYRAMID:
    284:       /*
    285:        4----
    286:        |\-\ \-----
    287:        | \ -\     \
    288:        |  1--\-----2
    289:        | /    \   /
    290:        |/      \ /
    291:        0--------3
    292:        */
    293:       if (numFaces) *numFaces = 5;
    294:       if (faceTypes) {
    295:         typesTmp[0] = DM_POLYTOPE_QUADRILATERAL;
    296:         typesTmp[1] = DM_POLYTOPE_TRIANGLE; typesTmp[2] = DM_POLYTOPE_TRIANGLE; typesTmp[3] = DM_POLYTOPE_TRIANGLE; typesTmp[4] = DM_POLYTOPE_TRIANGLE;
    297:         *faceTypes = typesTmp;
    298:       }
    299:       if (faceSizes) {
    300:         sizesTmp[0] = 4;
    301:         sizesTmp[1] = 3; sizesTmp[2] = 3; sizesTmp[3] = 3; sizesTmp[4] = 3;
    302:         *faceSizes = sizesTmp;
    303:       }
    304:       if (faces) {
    305:         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2]; facesTmp[3]  = cone[3]; /* Bottom */
    306:         facesTmp[4]  = cone[0]; facesTmp[5]  = cone[3]; facesTmp[6]  = cone[4];                         /* Front */
    307:         facesTmp[7]  = cone[3]; facesTmp[8]  = cone[2]; facesTmp[9]  = cone[4];                         /* Right */
    308:         facesTmp[10] = cone[2]; facesTmp[11] = cone[1]; facesTmp[12] = cone[4];                         /* Back */
    309:         facesTmp[13] = cone[1]; facesTmp[14] = cone[0]; facesTmp[15] = cone[4];                         /* Left */
    310:         *faces = facesTmp;
    311:       }
    312:       break;
    313:     default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No face description for cell type %s", DMPolytopeTypes[ct]);
    314:   }
    315:   return(0);
    316: }

    318: PetscErrorCode DMPlexRestoreRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[])
    319: {
    320:   PetscErrorCode  ierr;

    323:   if (faceTypes) {DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faceTypes);}
    324:   if (faceSizes) {DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faceSizes);}
    325:   if (faces)     {DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faces);}
    326:   return(0);
    327: }

    329: /* This interpolates faces for cells at some stratum */
    330: static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
    331: {
    332:   DMLabel        ctLabel;
    333:   PetscHashIJKL  faceTable;
    334:   PetscInt       faceTypeNum[DM_NUM_POLYTOPES];
    335:   PetscInt       depth, d, Np, cStart, cEnd, c, fStart, fEnd;

    339:   DMPlexGetDepth(dm, &depth);
    340:   PetscHashIJKLCreate(&faceTable);
    341:   PetscArrayzero(faceTypeNum, DM_NUM_POLYTOPES);
    342:   DMPlexGetDepthStratum(dm, cellDepth, &cStart, &cEnd);
    343:   /* Number new faces and save face vertices in hash table */
    344:   DMPlexGetDepthStratum(dm, depth > cellDepth ? cellDepth : 0, NULL, &fStart);
    345:   fEnd = fStart;
    346:   for (c = cStart; c < cEnd; ++c) {
    347:     const PetscInt       *cone, *faceSizes, *faces;
    348:     const DMPolytopeType *faceTypes;
    349:     DMPolytopeType        ct;
    350:     PetscInt              numFaces, cf, foff = 0;

    352:     DMPlexGetCellType(dm, c, &ct);
    353:     DMPlexGetCone(dm, c, &cone);
    354:     DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
    355:     for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
    356:       const PetscInt       faceSize = faceSizes[cf];
    357:       const DMPolytopeType faceType = faceTypes[cf];
    358:       const PetscInt      *face     = &faces[foff];
    359:       PetscHashIJKLKey     key;
    360:       PetscHashIter        iter;
    361:       PetscBool            missing;

    363:       if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %D > 4", faceSize);
    364:       key.i = face[0];
    365:       key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
    366:       key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
    367:       key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
    368:       PetscSortInt(faceSize, (PetscInt *) &key);
    369:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
    370:       if (missing) {
    371:         PetscHashIJKLIterSet(faceTable, iter, fEnd++);
    372:         ++faceTypeNum[faceType];
    373:       }
    374:     }
    375:     DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
    376:   }
    377:   /* We need to number faces contiguously among types */
    378:   {
    379:     PetscInt faceTypeStart[DM_NUM_POLYTOPES], ct, numFT = 0;

    381:     for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) {if (faceTypeNum[ct]) ++numFT; faceTypeStart[ct] = 0;}
    382:     if (numFT > 1) {
    383:       PetscHashIJKLClear(faceTable);
    384:       faceTypeStart[0] = fStart;
    385:       for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) faceTypeStart[ct] = faceTypeStart[ct-1] + faceTypeNum[ct-1];
    386:       for (c = cStart; c < cEnd; ++c) {
    387:         const PetscInt       *cone, *faceSizes, *faces;
    388:         const DMPolytopeType *faceTypes;
    389:         DMPolytopeType        ct;
    390:         PetscInt              numFaces, cf, foff = 0;

    392:         DMPlexGetCellType(dm, c, &ct);
    393:         DMPlexGetCone(dm, c, &cone);
    394:         DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
    395:         for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
    396:           const PetscInt       faceSize = faceSizes[cf];
    397:           const DMPolytopeType faceType = faceTypes[cf];
    398:           const PetscInt      *face     = &faces[foff];
    399:           PetscHashIJKLKey     key;
    400:           PetscHashIter        iter;
    401:           PetscBool            missing;

    403:           if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %D > 4", faceSize);
    404:           key.i = face[0];
    405:           key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
    406:           key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
    407:           key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
    408:           PetscSortInt(faceSize, (PetscInt *) &key);
    409:           PetscHashIJKLPut(faceTable, key, &iter, &missing);
    410:           if (missing) {PetscHashIJKLIterSet(faceTable, iter, faceTypeStart[faceType]++);}
    411:         }
    412:         DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
    413:       }
    414:       for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) {
    415:         if (faceTypeStart[ct] != faceTypeStart[ct-1] + faceTypeNum[ct]) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent numbering for cell type %s, %D != %D + %D", DMPolytopeTypes[ct], faceTypeStart[ct], faceTypeStart[ct-1], faceTypeNum[ct]);
    416:       }
    417:     }
    418:   }
    419:   /* Add new points, always at the end of the numbering */
    420:   DMPlexGetChart(dm, NULL, &Np);
    421:   DMPlexSetChart(idm, 0, Np + (fEnd - fStart));
    422:   /* Set cone sizes */
    423:   /*   Must create the celltype label here so that we do not automatically try to compute the types */
    424:   DMCreateLabel(idm, "celltype");
    425:   DMPlexGetCellTypeLabel(idm, &ctLabel);
    426:   for (d = 0; d <= depth; ++d) {
    427:     DMPolytopeType ct;
    428:     PetscInt       coneSize, pStart, pEnd, p;

    430:     if (d == cellDepth) continue;
    431:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
    432:     for (p = pStart; p < pEnd; ++p) {
    433:       DMPlexGetConeSize(dm, p, &coneSize);
    434:       DMPlexSetConeSize(idm, p, coneSize);
    435:       DMPlexGetCellType(dm, p, &ct);
    436:       DMPlexSetCellType(idm, p, ct);
    437:     }
    438:   }
    439:   for (c = cStart; c < cEnd; ++c) {
    440:     const PetscInt       *cone, *faceSizes, *faces;
    441:     const DMPolytopeType *faceTypes;
    442:     DMPolytopeType        ct;
    443:     PetscInt              numFaces, cf, foff = 0;

    445:     DMPlexGetCellType(dm, c, &ct);
    446:     DMPlexGetCone(dm, c, &cone);
    447:     DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
    448:     DMPlexSetCellType(idm, c, ct);
    449:     DMPlexSetConeSize(idm, c, numFaces);
    450:     for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
    451:       const PetscInt       faceSize = faceSizes[cf];
    452:       const DMPolytopeType faceType = faceTypes[cf];
    453:       const PetscInt      *face     = &faces[foff];
    454:       PetscHashIJKLKey     key;
    455:       PetscHashIter        iter;
    456:       PetscBool            missing;
    457:       PetscInt             f;

    459:       if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %D > 4", faceSize);
    460:       key.i = face[0];
    461:       key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
    462:       key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
    463:       key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
    464:       PetscSortInt(faceSize, (PetscInt *) &key);
    465:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
    466:       if (missing) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing face (cell %D, lf %D)", c, cf);
    467:       PetscHashIJKLIterGet(faceTable, iter, &f);
    468:       DMPlexSetConeSize(idm, f, faceSize);
    469:       DMPlexSetCellType(idm, f, faceType);
    470:     }
    471:     DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
    472:   }
    473:   DMSetUp(idm);
    474:   /* Initialize cones so we do not need the bash table to tell us that a cone has been set */
    475:   {
    476:     PetscSection cs;
    477:     PetscInt    *cones, csize;

    479:     DMPlexGetConeSection(idm, &cs);
    480:     DMPlexGetCones(idm, &cones);
    481:     PetscSectionGetStorageSize(cs, &csize);
    482:     for (c = 0; c < csize; ++c) cones[c] = -1;
    483:   }
    484:   /* Set cones */
    485:   for (d = 0; d <= depth; ++d) {
    486:     const PetscInt *cone;
    487:     PetscInt        pStart, pEnd, p;

    489:     if (d == cellDepth) continue;
    490:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
    491:     for (p = pStart; p < pEnd; ++p) {
    492:       DMPlexGetCone(dm, p, &cone);
    493:       DMPlexSetCone(idm, p, cone);
    494:       DMPlexGetConeOrientation(dm, p, &cone);
    495:       DMPlexSetConeOrientation(idm, p, cone);
    496:     }
    497:   }
    498:   for (c = cStart; c < cEnd; ++c) {
    499:     const PetscInt       *cone, *faceSizes, *faces;
    500:     const DMPolytopeType *faceTypes;
    501:     DMPolytopeType        ct;
    502:     PetscInt              numFaces, cf, foff = 0;

    504:     DMPlexGetCellType(dm, c, &ct);
    505:     DMPlexGetCone(dm, c, &cone);
    506:     DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
    507:     for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
    508:       DMPolytopeType   faceType = faceTypes[cf];
    509:       const PetscInt   faceSize = faceSizes[cf];
    510:       const PetscInt  *face     = &faces[foff];
    511:       const PetscInt  *fcone;
    512:       PetscHashIJKLKey key;
    513:       PetscHashIter    iter;
    514:       PetscBool        missing;
    515:       PetscInt         f;

    517:       if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %D > 4", faceSize);
    518:       key.i = face[0];
    519:       key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
    520:       key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
    521:       key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
    522:       PetscSortInt(faceSize, (PetscInt *) &key);
    523:       PetscHashIJKLPut(faceTable, key, &iter, &missing);
    524:       PetscHashIJKLIterGet(faceTable, iter, &f);
    525:       DMPlexInsertCone(idm, c, cf, f);
    526:       DMPlexGetCone(idm, f, &fcone);
    527:       if (fcone[0] < 0) {DMPlexSetCone(idm, f, face);}
    528:       /* TODO This should be unnecessary since we have autoamtic orientation */
    529:       {
    530:         /* when matching hybrid faces in 3D, only few cases are possible.
    531:            Face traversal however can no longer follow the usual convention, this seems a serious issue to me */
    532:         PetscInt        tquad_map[4][4] = { {PETSC_MIN_INT,            0,PETSC_MIN_INT,PETSC_MIN_INT},
    533:                                             {           -1,PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT},
    534:                                             {PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT,            1},
    535:                                             {PETSC_MIN_INT,PETSC_MIN_INT,           -2,PETSC_MIN_INT} };
    536:         PetscInt        i, i2, j;
    537:         const PetscInt *cone;
    538:         PetscInt        coneSize, ornt;

    540:         /* Orient face: Do not allow reverse orientation at the first vertex */
    541:         DMPlexGetConeSize(idm, f, &coneSize);
    542:         DMPlexGetCone(idm, f, &cone);
    543:         if (coneSize != faceSize) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of face vertices %D for face %D should be %D", coneSize, f, faceSize);
    544:         /* - First find the initial vertex */
    545:         for (i = 0; i < faceSize; ++i) if (face[0] == cone[i]) break;
    546:         /* If we want to compare tensor faces to regular faces, we have to flip them and take the else branch here */
    547:         if (faceType == DM_POLYTOPE_SEG_PRISM_TENSOR) {
    548:           /* find the second vertex */
    549:           for (i2 = 0; i2 < faceSize; ++i2) if (face[1] == cone[i2]) break;
    550:           switch (faceSize) {
    551:           case 2:
    552:             ornt = i ? -2 : 0;
    553:             break;
    554:           case 4:
    555:             ornt = tquad_map[i][i2];
    556:             break;
    557:           default: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unhandled face size %D for face %D in cell %D", faceSize, f, c);
    558:           }
    559:         } else {
    560:           /* Try forward comparison */
    561:           for (j = 0; j < faceSize; ++j) if (face[j] != cone[(i+j)%faceSize]) break;
    562:           if (j == faceSize) {
    563:             if ((faceSize == 2) && (i == 1)) ornt = -2;
    564:             else                             ornt = i;
    565:           } else {
    566:             /* Try backward comparison */
    567:             for (j = 0; j < faceSize; ++j) if (face[j] != cone[(i+faceSize-j)%faceSize]) break;
    568:             if (j == faceSize) {
    569:               if (i == 0) ornt = -faceSize;
    570:               else        ornt = -i;
    571:             } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine orientation of face %D in cell %D", f, c);
    572:           }
    573:         }
    574:         DMPlexInsertConeOrientation(idm, c, cf, ornt);
    575:       }
    576:     }
    577:     DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);
    578:   }
    579:   PetscHashIJKLDestroy(&faceTable);
    580:   DMPlexSymmetrize(idm);
    581:   DMPlexStratify(idm);
    582:   return(0);
    583: }

    585: static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[])
    586: {
    587:   PetscInt            nleaves;
    588:   PetscInt            nranks;
    589:   const PetscMPIInt  *ranks=NULL;
    590:   const PetscInt     *roffset=NULL, *rmine=NULL, *rremote=NULL;
    591:   PetscInt            n, o, r;
    592:   PetscErrorCode      ierr;

    595:   PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);
    596:   nleaves = roffset[nranks];
    597:   PetscMalloc2(nleaves, rmine1, nleaves, rremote1);
    598:   for (r=0; r<nranks; r++) {
    599:     /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote
    600:        - to unify order with the other side */
    601:     o = roffset[r];
    602:     n = roffset[r+1] - o;
    603:     PetscArraycpy(&(*rmine1)[o], &rmine[o], n);
    604:     PetscArraycpy(&(*rremote1)[o], &rremote[o], n);
    605:     PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);
    606:   }
    607:   return(0);
    608: }

    610: PetscErrorCode DMPlexOrientInterface_Internal(DM dm)
    611: {
    612:   /* Here we only compare first 2 points of the cone. Full cone size would lead to stronger self-checking. */
    613:   PetscInt          mainCone[2];
    614:   PetscInt          (*roots)[2], (*leaves)[2];
    615:   PetscMPIInt       (*rootsRanks)[2], (*leavesRanks)[2];

    617:   PetscSF           sf=NULL;
    618:   const PetscInt    *locals=NULL;
    619:   const PetscSFNode *remotes=NULL;
    620:   PetscInt           nroots, nleaves, p, c;
    621:   PetscInt           nranks, n, o, r;
    622:   const PetscMPIInt *ranks=NULL;
    623:   const PetscInt    *roffset=NULL;
    624:   PetscInt          *rmine1=NULL, *rremote1=NULL; /* rmine and rremote copies simultaneously sorted by rank and rremote */
    625:   const PetscInt    *cone=NULL;
    626:   PetscInt           coneSize, ind0;
    627:   MPI_Comm           comm;
    628:   PetscMPIInt        rank,size;
    629:   PetscInt           debug = 0;
    630:   PetscErrorCode     ierr;

    633:   DMGetPointSF(dm, &sf);
    634:   PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes);
    635:   if (nroots < 0) return(0);
    636:   PetscSFSetUp(sf);
    637:   PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL);
    638:   DMViewFromOptions(dm, NULL, "-before_fix_dm_view");
    639:   if (PetscDefined(USE_DEBUG)) {DMPlexCheckPointSF(dm);}
    640:   SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1);
    641:   PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks);
    642:   PetscObjectGetComm((PetscObject) dm, &comm);
    643:   MPI_Comm_rank(comm, &rank);
    644:   MPI_Comm_size(comm, &size);
    645:   for (p = 0; p < nroots; ++p) {
    646:     DMPlexGetConeSize(dm, p, &coneSize);
    647:     DMPlexGetCone(dm, p, &cone);
    648:     if (coneSize < 2) {
    649:       for (c = 0; c < 2; c++) {
    650:         roots[p][c] = -1;
    651:         rootsRanks[p][c] = -1;
    652:       }
    653:       continue;
    654:     }
    655:     /* Translate all points to root numbering */
    656:     for (c = 0; c < 2; c++) {
    657:       PetscFindInt(cone[c], nleaves, locals, &ind0);
    658:       if (ind0 < 0) {
    659:         roots[p][c] = cone[c];
    660:         rootsRanks[p][c] = rank;
    661:       } else {
    662:         roots[p][c] = remotes[ind0].index;
    663:         rootsRanks[p][c] = remotes[ind0].rank;
    664:       }
    665:     }
    666:   }
    667:   for (p = 0; p < nroots; ++p) {
    668:     for (c = 0; c < 2; c++) {
    669:       leaves[p][c] = -2;
    670:       leavesRanks[p][c] = -2;
    671:     }
    672:   }
    673:   PetscSFBcastBegin(sf, MPIU_2INT, roots, leaves,MPI_REPLACE);
    674:   PetscSFBcastBegin(sf, MPI_2INT, rootsRanks, leavesRanks,MPI_REPLACE);
    675:   PetscSFBcastEnd(sf, MPIU_2INT, roots, leaves,MPI_REPLACE);
    676:   PetscSFBcastEnd(sf, MPI_2INT, rootsRanks, leavesRanks,MPI_REPLACE);
    677:   if (debug) {PetscSynchronizedFlush(comm, NULL);}
    678:   if (debug && rank == 0) {PetscSynchronizedPrintf(comm, "Referenced roots\n");}
    679:   for (p = 0; p < nroots; ++p) {
    680:     if (leaves[p][0] < 0) continue;
    681:     DMPlexGetConeSize(dm, p, &coneSize);
    682:     DMPlexGetCone(dm, p, &cone);
    683:     if (debug) {PetscSynchronizedPrintf(comm, "[%d]  %4D: cone=[%4D %4D] roots=[(%d,%4D) (%d,%4D)] leaves=[(%d,%4D) (%d,%4D)]", rank, p, cone[0], cone[1], rootsRanks[p][0], roots[p][0], rootsRanks[p][1], roots[p][1], leavesRanks[p][0], leaves[p][0], leavesRanks[p][1], leaves[p][1]);}
    684:     if ((leaves[p][0] != roots[p][0]) || (leaves[p][1] != roots[p][1]) || (leavesRanks[p][0] != rootsRanks[p][0]) || (leavesRanks[p][1] != rootsRanks[p][1])) {
    685:       /* Translate these two leaves to my cone points; mainCone means desired order p's cone points */
    686:       for (c = 0; c < 2; c++) {
    687:         if (leavesRanks[p][c] == rank) {
    688:           /* A local leave is just taken as it is */
    689:           mainCone[c] = leaves[p][c];
    690:           continue;
    691:         }
    692:         /* Find index of rank leavesRanks[p][c] among remote ranks */
    693:         /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
    694:         PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &r);
    695:         if (PetscUnlikely(r < 0)) SETERRQ7(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D cone[%D]=%D root (%d,%D) leave (%d,%D): leave rank not found among remote ranks",p,c,cone[c],rootsRanks[p][c],roots[p][c],leavesRanks[p][c],leaves[p][c]);
    696:         if (PetscUnlikely(ranks[r] < 0 || ranks[r] >= size)) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_PLIB, "p=%D c=%D commsize=%d: ranks[%D] = %d makes no sense",p,c,size,r,ranks[r]);
    697:         /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */
    698:         o = roffset[r];
    699:         n = roffset[r+1] - o;
    700:         PetscFindInt(leaves[p][c], n, &rremote1[o], &ind0);
    701:         if (PetscUnlikely(ind0 < 0)) SETERRQ7(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D cone[%D]=%D root (%d,%D) leave (%d,%D): 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]);
    702:         /* Get the corresponding local point */
    703:         mainCone[c] = rmine1[o+ind0];
    704:       }
    705:       if (debug) {PetscSynchronizedPrintf(comm, " mainCone=[%4D %4D]\n", mainCone[0], mainCone[1]);}
    706:       /* Set the desired order of p's cone points and fix orientations accordingly */
    707:       /* Vaclav's note: Here we only compare first 2 points of the cone. Full cone size would lead to stronger self-checking. */
    708:       DMPlexOrientCell(dm, p, 2, mainCone);
    709:     } else if (debug) {PetscSynchronizedPrintf(comm, " ==\n");}
    710:   }
    711:   if (debug) {
    712:     PetscSynchronizedFlush(comm, NULL);
    713:     MPI_Barrier(comm);
    714:   }
    715:   DMViewFromOptions(dm, NULL, "-after_fix_dm_view");
    716:   PetscFree4(roots, leaves, rootsRanks, leavesRanks);
    717:   PetscFree2(rmine1, rremote1);
    718:   return(0);
    719: }

    721: static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[])
    722: {
    723:   PetscInt       idx;
    724:   PetscMPIInt    rank;
    725:   PetscBool      flg;

    729:   PetscOptionsHasName(NULL, NULL, opt, &flg);
    730:   if (!flg) return(0);
    731:   MPI_Comm_rank(comm, &rank);
    732:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
    733:   for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]%s %D %s %D\n", rank, idxname, idx, valname, a[idx]);}
    734:   PetscSynchronizedFlush(comm, NULL);
    735:   return(0);
    736: }

    738: static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[])
    739: {
    740:   PetscInt       idx;
    741:   PetscMPIInt    rank;
    742:   PetscBool      flg;

    746:   PetscOptionsHasName(NULL, NULL, opt, &flg);
    747:   if (!flg) return(0);
    748:   MPI_Comm_rank(comm, &rank);
    749:   PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);
    750:   if (idxname) {
    751:     for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]%s %D rank %D index %D\n", rank, idxname, idx, a[idx].rank, a[idx].index);}
    752:   } else {
    753:     for (idx = 0; idx < n; ++idx) {PetscSynchronizedPrintf(comm, "[%d]rank %D index %D\n", rank, a[idx].rank, a[idx].index);}
    754:   }
    755:   PetscSynchronizedFlush(comm, NULL);
    756:   return(0);
    757: }

    759: static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint)
    760: {
    761:   PetscSF         sf;
    762:   const PetscInt *locals;
    763:   PetscMPIInt     rank;
    764:   PetscErrorCode  ierr;

    767:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
    768:   DMGetPointSF(dm, &sf);
    769:   PetscSFGetGraph(sf, NULL, NULL, &locals, NULL);
    770:   if (remotePoint.rank == rank) {
    771:     *localPoint = remotePoint.index;
    772:   } else {
    773:     PetscHashIJKey key;
    774:     PetscInt       l;

    776:     key.i = remotePoint.index;
    777:     key.j = remotePoint.rank;
    778:     PetscHMapIJGet(remotehash, key, &l);
    779:     if (l >= 0) {
    780:       *localPoint = locals[l];
    781:     } else PetscFunctionReturn(1);
    782:   }
    783:   return(0);
    784: }

    786: static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint)
    787: {
    788:   PetscSF            sf;
    789:   const PetscInt    *locals, *rootdegree;
    790:   const PetscSFNode *remotes;
    791:   PetscInt           Nl, l;
    792:   PetscMPIInt        rank;
    793:   PetscErrorCode     ierr;

    796:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
    797:   DMGetPointSF(dm, &sf);
    798:   PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes);
    799:   if (Nl < 0) goto owned;
    800:   PetscSFComputeDegreeBegin(sf, &rootdegree);
    801:   PetscSFComputeDegreeEnd(sf, &rootdegree);
    802:   if (rootdegree[localPoint]) goto owned;
    803:   PetscFindInt(localPoint, Nl, locals, &l);
    804:   if (l < 0) PetscFunctionReturn(1);
    805:   *remotePoint = remotes[l];
    806:   return(0);
    807:   owned:
    808:   remotePoint->rank  = rank;
    809:   remotePoint->index = localPoint;
    810:   return(0);
    811: }

    814: static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared)
    815: {
    816:   PetscSF         sf;
    817:   const PetscInt *locals, *rootdegree;
    818:   PetscInt        Nl, idx;
    819:   PetscErrorCode  ierr;

    822:   *isShared = PETSC_FALSE;
    823:   DMGetPointSF(dm, &sf);
    824:   PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL);
    825:   if (Nl < 0) return(0);
    826:   PetscFindInt(p, Nl, locals, &idx);
    827:   if (idx >= 0) {*isShared = PETSC_TRUE; return(0);}
    828:   PetscSFComputeDegreeBegin(sf, &rootdegree);
    829:   PetscSFComputeDegreeEnd(sf, &rootdegree);
    830:   if (rootdegree[p] > 0) *isShared = PETSC_TRUE;
    831:   return(0);
    832: }

    834: static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared)
    835: {
    836:   const PetscInt *cone;
    837:   PetscInt        coneSize, c;
    838:   PetscBool       cShared = PETSC_TRUE;
    839:   PetscErrorCode  ierr;

    842:   DMPlexGetConeSize(dm, p, &coneSize);
    843:   DMPlexGetCone(dm, p, &cone);
    844:   for (c = 0; c < coneSize; ++c) {
    845:     PetscBool pointShared;

    847:     DMPlexPointIsShared(dm, cone[c], &pointShared);
    848:     cShared = (PetscBool) (cShared && pointShared);
    849:   }
    850:   *isShared = coneSize ? cShared : PETSC_FALSE;
    851:   return(0);
    852: }

    854: static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin)
    855: {
    856:   const PetscInt *cone;
    857:   PetscInt        coneSize, c;
    858:   PetscSFNode     cmin = {PETSC_MAX_INT, PETSC_MAX_INT}, missing = {-1, -1};
    859:   PetscErrorCode  ierr;

    862:   DMPlexGetConeSize(dm, p, &coneSize);
    863:   DMPlexGetCone(dm, p, &cone);
    864:   for (c = 0; c < coneSize; ++c) {
    865:     PetscSFNode rcp;

    867:     DMPlexMapToGlobalPoint(dm, cone[c], &rcp);
    868:     if (ierr) {
    869:       cmin = missing;
    870:     } else {
    871:       cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin;
    872:     }
    873:   }
    874:   *cpmin = coneSize ? cmin : missing;
    875:   return(0);
    876: }

    878: /*
    879:   Each shared face has an entry in the candidates array:
    880:     (-1, coneSize-1), {(global cone point)}
    881:   where the set is missing the point p which we use as the key for the face
    882: */
    883: static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug)
    884: {
    885:   MPI_Comm        comm;
    886:   const PetscInt *support;
    887:   PetscInt        supportSize, s, off = 0, idx = 0, overlap, cellHeight, height;
    888:   PetscMPIInt     rank;
    889:   PetscErrorCode  ierr;

    892:   PetscObjectGetComm((PetscObject) dm, &comm);
    893:   MPI_Comm_rank(comm, &rank);
    894:   DMPlexGetOverlap(dm, &overlap);
    895:   DMPlexGetVTKCellHeight(dm, &cellHeight);
    896:   DMPlexGetPointHeight(dm, p, &height);
    897:   if (!overlap && height <= cellHeight+1) {
    898:     /* cells can't be shared for non-overlapping meshes */
    899:     if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Skipping face %D to avoid adding cell to hashmap since this is nonoverlapping mesh\n", rank, p);}
    900:     return(0);
    901:   }
    902:   DMPlexGetSupportSize(dm, p, &supportSize);
    903:   DMPlexGetSupport(dm, p, &support);
    904:   if (candidates) {PetscSectionGetOffset(candidateSection, p, &off);}
    905:   for (s = 0; s < supportSize; ++s) {
    906:     const PetscInt  face = support[s];
    907:     const PetscInt *cone;
    908:     PetscSFNode     cpmin={-1,-1}, rp={-1,-1};
    909:     PetscInt        coneSize, c, f;
    910:     PetscBool       isShared = PETSC_FALSE;
    911:     PetscHashIJKey  key;

    913:     /* Only add point once */
    914:     if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Support face %D\n", rank, face);}
    915:     key.i = p;
    916:     key.j = face;
    917:     PetscHMapIJGet(faceHash, key, &f);
    918:     if (f >= 0) continue;
    919:     DMPlexConeIsShared(dm, face, &isShared);
    920:     DMPlexGetConeMinimum(dm, face, &cpmin);
    921:     DMPlexMapToGlobalPoint(dm, p, &rp);
    922:     if (debug) {
    923:       PetscSynchronizedPrintf(comm, "[%d]      Face point %D is shared: %d\n", rank, face, (int) isShared);
    924:       PetscSynchronizedPrintf(comm, "[%d]      Global point (%D, %D) Min Cone Point (%D, %D)\n", rank, rp.rank, rp.index, cpmin.rank, cpmin.index);
    925:     }
    926:     if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) {
    927:       PetscHMapIJSet(faceHash, key, p);
    928:       if (candidates) {
    929:         if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Adding shared face %D at idx %D\n[%d]     ", rank, face, idx, rank);}
    930:         DMPlexGetConeSize(dm, face, &coneSize);
    931:         DMPlexGetCone(dm, face, &cone);
    932:         candidates[off+idx].rank    = -1;
    933:         candidates[off+idx++].index = coneSize-1;
    934:         candidates[off+idx].rank    = rank;
    935:         candidates[off+idx++].index = face;
    936:         for (c = 0; c < coneSize; ++c) {
    937:           const PetscInt cp = cone[c];

    939:           if (cp == p) continue;
    940:           DMPlexMapToGlobalPoint(dm, cp, &candidates[off+idx]);
    941:           if (debug) {PetscSynchronizedPrintf(comm, " (%D,%D)", candidates[off+idx].rank, candidates[off+idx].index);}
    942:           ++idx;
    943:         }
    944:         if (debug) {PetscSynchronizedPrintf(comm, "\n");}
    945:       } else {
    946:         /* Add cone size to section */
    947:         if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Scheduling shared face %D\n", rank, face);}
    948:         DMPlexGetConeSize(dm, face, &coneSize);
    949:         PetscHMapIJSet(faceHash, key, p);
    950:         PetscSectionAddDof(candidateSection, p, coneSize+1);
    951:       }
    952:     }
    953:   }
    954:   return(0);
    955: }

    957: /*@
    958:   DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the PointSF in parallel, following local interpolation

    960:   Collective on dm

    962:   Input Parameters:
    963: + dm      - The interpolated DM
    964: - pointSF - The initial SF without interpolated points

    966:   Output Parameter:
    967: . pointSF - The SF including interpolated points

    969:   Level: developer

    971:    Note: All 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

    973: .seealso: DMPlexInterpolate(), DMPlexUninterpolate()
    974: @*/
    975: PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
    976: {
    977:   MPI_Comm           comm;
    978:   PetscHMapIJ        remoteHash;
    979:   PetscHMapI         claimshash;
    980:   PetscSection       candidateSection, candidateRemoteSection, claimSection;
    981:   PetscSFNode       *candidates, *candidatesRemote, *claims;
    982:   const PetscInt    *localPoints, *rootdegree;
    983:   const PetscSFNode *remotePoints;
    984:   PetscInt           ov, Nr, r, Nl, l;
    985:   PetscInt           candidatesSize, candidatesRemoteSize, claimsSize;
    986:   PetscBool          flg, debug = PETSC_FALSE;
    987:   PetscMPIInt        rank;
    988:   PetscErrorCode     ierr;

    993:   DMPlexIsDistributed(dm, &flg);
    994:   if (!flg) return(0);
    995:   /* Set initial SF so that lower level queries work */
    996:   DMSetPointSF(dm, pointSF);
    997:   PetscObjectGetComm((PetscObject) dm, &comm);
    998:   MPI_Comm_rank(comm, &rank);
    999:   DMPlexGetOverlap(dm, &ov);
    1000:   if (ov) SETERRQ(comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet");
    1001:   PetscOptionsHasName(NULL, ((PetscObject) dm)->prefix, "-dmplex_interp_debug", &debug);
    1002:   PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_interp_pre_view");
    1003:   PetscObjectViewFromOptions((PetscObject) pointSF, NULL, "-petscsf_interp_pre_view");
    1004:   PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);
    1005:   /* Step 0: Precalculations */
    1006:   PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints);
    1007:   if (Nr < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set");
    1008:   PetscHMapIJCreate(&remoteHash);
    1009:   for (l = 0; l < Nl; ++l) {
    1010:     PetscHashIJKey key;
    1011:     key.i = remotePoints[l].index;
    1012:     key.j = remotePoints[l].rank;
    1013:     PetscHMapIJSet(remoteHash, key, l);
    1014:   }
    1015:   /*   Compute root degree to identify shared points */
    1016:   PetscSFComputeDegreeBegin(pointSF, &rootdegree);
    1017:   PetscSFComputeDegreeEnd(pointSF, &rootdegree);
    1018:   IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree);
    1019:   /*
    1020:   1) Loop over each leaf point $p$ at depth $d$ in the SF
    1021:   \item Get set $F(p)$ of faces $f$ in the support of $p$ for which
    1022:   \begin{itemize}
    1023:     \item all cone points of $f$ are shared
    1024:     \item $p$ is the cone point with smallest canonical number
    1025:   \end{itemize}
    1026:   \item Send $F(p)$ and the cone of each face to the active root point $r(p)$
    1027:   \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
    1028:   \item Send the root face from the root back to all leaf process
    1029:   \item Leaf processes add the shared face to the SF
    1030:   */
    1031:   /* Step 1: Construct section+SFNode array
    1032:        The section has entries for all shared faces for which we have a leaf point in the cone
    1033:        The array holds candidate shared faces, each face is refered to by the leaf point */
    1034:   PetscSectionCreate(comm, &candidateSection);
    1035:   PetscSectionSetChart(candidateSection, 0, Nr);
    1036:   {
    1037:     PetscHMapIJ faceHash;

    1039:     PetscHMapIJCreate(&faceHash);
    1040:     for (l = 0; l < Nl; ++l) {
    1041:       const PetscInt p = localPoints[l];

    1043:       if (debug) {PetscSynchronizedPrintf(comm, "[%d]  First pass leaf point %D\n", rank, p);}
    1044:       DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug);
    1045:     }
    1046:     PetscHMapIJClear(faceHash);
    1047:     PetscSectionSetUp(candidateSection);
    1048:     PetscSectionGetStorageSize(candidateSection, &candidatesSize);
    1049:     PetscMalloc1(candidatesSize, &candidates);
    1050:     for (l = 0; l < Nl; ++l) {
    1051:       const PetscInt p = localPoints[l];

    1053:       if (debug) {PetscSynchronizedPrintf(comm, "[%d]  Second pass leaf point %D\n", rank, p);}
    1054:       DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug);
    1055:     }
    1056:     PetscHMapIJDestroy(&faceHash);
    1057:     if (debug) {PetscSynchronizedFlush(comm, NULL);}
    1058:   }
    1059:   PetscObjectSetName((PetscObject) candidateSection, "Candidate Section");
    1060:   PetscObjectViewFromOptions((PetscObject) candidateSection, NULL, "-petscsection_interp_candidate_view");
    1061:   SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates);
    1062:   /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
    1063:   /*   Note that this section is indexed by offsets into leaves, not by point number */
    1064:   {
    1065:     PetscSF   sfMulti, sfInverse, sfCandidates;
    1066:     PetscInt *remoteOffsets;

    1068:     PetscSFGetMultiSF(pointSF, &sfMulti);
    1069:     PetscSFCreateInverseSF(sfMulti, &sfInverse);
    1070:     PetscSectionCreate(comm, &candidateRemoteSection);
    1071:     PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection);
    1072:     PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates);
    1073:     PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize);
    1074:     PetscMalloc1(candidatesRemoteSize, &candidatesRemote);
    1075:     PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote,MPI_REPLACE);
    1076:     PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote,MPI_REPLACE);
    1077:     PetscSFDestroy(&sfInverse);
    1078:     PetscSFDestroy(&sfCandidates);
    1079:     PetscFree(remoteOffsets);

    1081:     PetscObjectSetName((PetscObject) candidateRemoteSection, "Remote Candidate Section");
    1082:     PetscObjectViewFromOptions((PetscObject) candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view");
    1083:     SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote);
    1084:   }
    1085:   /* 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 */
    1086:   {
    1087:     PetscHashIJKLRemote faceTable;
    1088:     PetscInt            idx, idx2;

    1090:     PetscHashIJKLRemoteCreate(&faceTable);
    1091:     /* There is a section point for every leaf attached to a given root point */
    1092:     for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) {
    1093:       PetscInt deg;

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

    1098:         PetscSectionGetDof(candidateRemoteSection, idx, &dof);
    1099:         PetscSectionGetOffset(candidateRemoteSection, idx, &offset);
    1100:         /* dof may include many faces from the remote process */
    1101:         for (d = 0; d < dof; ++d) {
    1102:           const PetscInt         hidx  = offset+d;
    1103:           const PetscInt         Np    = candidatesRemote[hidx].index+1;
    1104:           const PetscSFNode      rface = candidatesRemote[hidx+1];
    1105:           const PetscSFNode     *fcone = &candidatesRemote[hidx+2];
    1106:           PetscSFNode            fcp0;
    1107:           const PetscSFNode      pmax  = {PETSC_MAX_INT, PETSC_MAX_INT};
    1108:           const PetscInt        *join  = NULL;
    1109:           PetscHashIJKLRemoteKey key;
    1110:           PetscHashIter          iter;
    1111:           PetscBool              missing;
    1112:           PetscInt               points[1024], p, joinSize;

    1114:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking face (%D, %D) at (%D, %D, %D) with cone size %D\n", rank, rface.rank, rface.index, r, idx, d, Np);}
    1115:           if (Np > 4) SETERRQ6(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle face (%D, %D) at (%D, %D, %D) with %D cone points", rface.rank, rface.index, r, idx, d, Np);
    1116:           fcp0.rank  = rank;
    1117:           fcp0.index = r;
    1118:           d += Np;
    1119:           /* Put remote face in hash table */
    1120:           key.i = fcp0;
    1121:           key.j = fcone[0];
    1122:           key.k = Np > 2 ? fcone[1] : pmax;
    1123:           key.l = Np > 3 ? fcone[2] : pmax;
    1124:           PetscSortSFNode(Np, (PetscSFNode *) &key);
    1125:           PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);
    1126:           if (missing) {
    1127:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Setting remote face (%D, %D)\n", rank, rface.index, rface.rank);}
    1128:             PetscHashIJKLRemoteIterSet(faceTable, iter, rface);
    1129:           } else {
    1130:             PetscSFNode oface;

    1132:             PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);
    1133:             if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) {
    1134:               if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Replacing with remote face (%D, %D)\n", rank, rface.index, rface.rank);}
    1135:               PetscHashIJKLRemoteIterSet(faceTable, iter, rface);
    1136:             }
    1137:           }
    1138:           /* Check for local face */
    1139:           points[0] = r;
    1140:           for (p = 1; p < Np; ++p) {
    1141:             DMPlexMapToLocalPoint(dm, remoteHash, fcone[p-1], &points[p]);
    1142:             if (ierr) break; /* We got a point not in our overlap */
    1143:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Checking local candidate %D\n", rank, points[p]);}
    1144:           }
    1145:           if (ierr) continue;
    1146:           DMPlexGetJoin(dm, Np, points, &joinSize, &join);
    1147:           if (joinSize == 1) {
    1148:             PetscSFNode lface;
    1149:             PetscSFNode oface;

    1151:             /* Always replace with local face */
    1152:             lface.rank  = rank;
    1153:             lface.index = join[0];
    1154:             PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);
    1155:             if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Replacing (%D, %D) with local face (%D, %D)\n", rank, oface.index, oface.rank, lface.index, lface.rank);}
    1156:             PetscHashIJKLRemoteIterSet(faceTable, iter, lface);
    1157:           }
    1158:           DMPlexRestoreJoin(dm, Np, points, &joinSize, &join);
    1159:         }
    1160:       }
    1161:       /* Put back faces for this root */
    1162:       for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) {
    1163:         PetscInt offset, dof, d;

    1165:         PetscSectionGetDof(candidateRemoteSection, idx2, &dof);
    1166:         PetscSectionGetOffset(candidateRemoteSection, idx2, &offset);
    1167:         /* dof may include many faces from the remote process */
    1168:         for (d = 0; d < dof; ++d) {
    1169:           const PetscInt         hidx  = offset+d;
    1170:           const PetscInt         Np    = candidatesRemote[hidx].index+1;
    1171:           const PetscSFNode     *fcone = &candidatesRemote[hidx+2];
    1172:           PetscSFNode            fcp0;
    1173:           const PetscSFNode      pmax  = {PETSC_MAX_INT, PETSC_MAX_INT};
    1174:           PetscHashIJKLRemoteKey key;
    1175:           PetscHashIter          iter;
    1176:           PetscBool              missing;

    1178:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]  Entering face at (%D, %D)\n", rank, r, idx);}
    1179:           if (Np > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %D cone points", Np);
    1180:           fcp0.rank  = rank;
    1181:           fcp0.index = r;
    1182:           d += Np;
    1183:           /* Find remote face in hash table */
    1184:           key.i = fcp0;
    1185:           key.j = fcone[0];
    1186:           key.k = Np > 2 ? fcone[1] : pmax;
    1187:           key.l = Np > 3 ? fcone[2] : pmax;
    1188:           PetscSortSFNode(Np, (PetscSFNode *) &key);
    1189:           if (debug) {PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d]    key (%D, %D) (%D, %D) (%D, %D) (%D, %D)\n", rank, key.i.rank, key.i.index, key.j.rank, key.j.index, key.k.rank, key.k.index, key.l.rank, key.l.index);}
    1190:           PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);
    1191:           if (missing) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %D Idx %D ought to have an associated face", r, idx2);
    1192:           else        {PetscHashIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx]);}
    1193:         }
    1194:       }
    1195:     }
    1196:     if (debug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);}
    1197:     PetscHashIJKLRemoteDestroy(&faceTable);
    1198:   }
    1199:   /* Step 4: Push back owned faces */
    1200:   {
    1201:     PetscSF      sfMulti, sfClaims, sfPointNew;
    1202:     PetscSFNode *remotePointsNew;
    1203:     PetscInt    *remoteOffsets, *localPointsNew;
    1204:     PetscInt     pStart, pEnd, r, NlNew, p;

    1206:     /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
    1207:     PetscSFGetMultiSF(pointSF, &sfMulti);
    1208:     PetscSectionCreate(comm, &claimSection);
    1209:     PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection);
    1210:     PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims);
    1211:     PetscSectionGetStorageSize(claimSection, &claimsSize);
    1212:     PetscMalloc1(claimsSize, &claims);
    1213:     for (p = 0; p < claimsSize; ++p) claims[p].rank = -1;
    1214:     PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims,MPI_REPLACE);
    1215:     PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims,MPI_REPLACE);
    1216:     PetscSFDestroy(&sfClaims);
    1217:     PetscFree(remoteOffsets);
    1218:     PetscObjectSetName((PetscObject) claimSection, "Claim Section");
    1219:     PetscObjectViewFromOptions((PetscObject) claimSection, NULL, "-petscsection_interp_claim_view");
    1220:     SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims);
    1221:     /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */
    1222:     /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */
    1223:     PetscHMapICreate(&claimshash);
    1224:     for (r = 0; r < Nr; ++r) {
    1225:       PetscInt dof, off, d;

    1227:       if (debug) {PetscSynchronizedPrintf(comm, "[%d]  Checking root for claims %D\n", rank, r);}
    1228:       PetscSectionGetDof(candidateSection, r, &dof);
    1229:       PetscSectionGetOffset(candidateSection, r, &off);
    1230:       for (d = 0; d < dof;) {
    1231:         if (claims[off+d].rank >= 0) {
    1232:           const PetscInt  faceInd = off+d;
    1233:           const PetscInt  Np      = candidates[off+d].index;
    1234:           const PetscInt *join    = NULL;
    1235:           PetscInt        joinSize, points[1024], c;

    1237:           if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Found claim for remote point (%D, %D)\n", rank, claims[faceInd].rank, claims[faceInd].index);}
    1238:           points[0] = r;
    1239:           if (debug) {PetscSynchronizedPrintf(comm, "[%d]      point %D\n", rank, points[0]);}
    1240:           for (c = 0, d += 2; c < Np; ++c, ++d) {
    1241:             DMPlexMapToLocalPoint(dm, remoteHash, candidates[off+d], &points[c+1]);
    1242:             if (debug) {PetscSynchronizedPrintf(comm, "[%d]      point %D\n", rank, points[c+1]);}
    1243:           }
    1244:           DMPlexGetJoin(dm, Np+1, points, &joinSize, &join);
    1245:           if (joinSize == 1) {
    1246:             if (claims[faceInd].rank == rank) {
    1247:               if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Ignoring local face %D for non-remote partner\n", rank, join[0]);}
    1248:             } else {
    1249:               if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Found local face %D\n", rank, join[0]);}
    1250:               PetscHMapISet(claimshash, join[0], faceInd);
    1251:             }
    1252:           } else {
    1253:             if (debug) {PetscSynchronizedPrintf(comm, "[%d]    Failed to find face\n", rank);}
    1254:           }
    1255:           DMPlexRestoreJoin(dm, Np+1, points, &joinSize, &join);
    1256:         } else {
    1257:           if (debug) {PetscSynchronizedPrintf(comm, "[%d]    No claim for point %D\n", rank, r);}
    1258:           d += claims[off+d].index+1;
    1259:         }
    1260:       }
    1261:     }
    1262:     if (debug) {PetscSynchronizedFlush(comm, NULL);}
    1263:     /* Step 6) Create new pointSF from hashed claims */
    1264:     PetscHMapIGetSize(claimshash, &NlNew);
    1265:     DMPlexGetChart(dm, &pStart, &pEnd);
    1266:     PetscMalloc1(Nl + NlNew, &localPointsNew);
    1267:     PetscMalloc1(Nl + NlNew, &remotePointsNew);
    1268:     for (l = 0; l < Nl; ++l) {
    1269:       localPointsNew[l] = localPoints[l];
    1270:       remotePointsNew[l].index = remotePoints[l].index;
    1271:       remotePointsNew[l].rank  = remotePoints[l].rank;
    1272:     }
    1273:     p = Nl;
    1274:     PetscHMapIGetKeys(claimshash, &p, localPointsNew);
    1275:     /* We sort new points, and assume they are numbered after all existing points */
    1276:     PetscSortInt(NlNew, &localPointsNew[Nl]);
    1277:     for (p = Nl; p < Nl + NlNew; ++p) {
    1278:       PetscInt off;
    1279:       PetscHMapIGet(claimshash, localPointsNew[p], &off);
    1280:       if (claims[off].rank < 0 || claims[off].index < 0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid claim for local point %D, (%D, %D)", localPointsNew[p], claims[off].rank, claims[off].index);
    1281:       remotePointsNew[p] = claims[off];
    1282:     }
    1283:     PetscSFCreate(comm, &sfPointNew);
    1284:     PetscSFSetGraph(sfPointNew, pEnd-pStart, Nl+NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
    1285:     PetscSFSetUp(sfPointNew);
    1286:     DMSetPointSF(dm, sfPointNew);
    1287:     PetscObjectViewFromOptions((PetscObject) sfPointNew, NULL, "-petscsf_interp_view");
    1288:     PetscSFDestroy(&sfPointNew);
    1289:     PetscHMapIDestroy(&claimshash);
    1290:   }
    1291:   PetscHMapIJDestroy(&remoteHash);
    1292:   PetscSectionDestroy(&candidateSection);
    1293:   PetscSectionDestroy(&candidateRemoteSection);
    1294:   PetscSectionDestroy(&claimSection);
    1295:   PetscFree(candidates);
    1296:   PetscFree(candidatesRemote);
    1297:   PetscFree(claims);
    1298:   PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);
    1299:   return(0);
    1300: }

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

    1305:   Collective on dm

    1307:   Input Parameters:
    1308: + dm - The DMPlex object with only cells and vertices
    1309: - dmInt - The interpolated DM

    1311:   Output Parameter:
    1312: . dmInt - The complete DMPlex object

    1314:   Level: intermediate

    1316:   Notes:
    1317:     It does not copy over the coordinates.

    1319:   Developer Notes:
    1320:     It sets plex->interpolated = DMPLEX_INTERPOLATED_FULL.

    1322: .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates()
    1323: @*/
    1324: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
    1325: {
    1326:   DMPlexInterpolatedFlag interpolated;
    1327:   DM             idm, odm = dm;
    1328:   PetscSF        sfPoint;
    1329:   PetscInt       depth, dim, d;
    1330:   const char    *name;
    1331:   PetscBool      flg=PETSC_TRUE;

    1337:   PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);
    1338:   DMPlexGetDepth(dm, &depth);
    1339:   DMGetDimension(dm, &dim);
    1340:   DMPlexIsInterpolated(dm, &interpolated);
    1341:   if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
    1342:   if (interpolated == DMPLEX_INTERPOLATED_FULL) {
    1343:     PetscObjectReference((PetscObject) dm);
    1344:     idm  = dm;
    1345:   } else {
    1346:     for (d = 1; d < dim; ++d) {
    1347:       /* Create interpolated mesh */
    1348:       DMCreate(PetscObjectComm((PetscObject)dm), &idm);
    1349:       DMSetType(idm, DMPLEX);
    1350:       DMSetDimension(idm, dim);
    1351:       if (depth > 0) {
    1352:         DMPlexInterpolateFaces_Internal(odm, 1, idm);
    1353:         DMGetPointSF(odm, &sfPoint);
    1354:         {
    1355:           /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */
    1356:           PetscInt nroots;
    1357:           PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
    1358:           if (nroots >= 0) {DMPlexInterpolatePointSF(idm, sfPoint);}
    1359:         }
    1360:       }
    1361:       if (odm != dm) {DMDestroy(&odm);}
    1362:       odm = idm;
    1363:     }
    1364:     PetscObjectGetName((PetscObject) dm,  &name);
    1365:     PetscObjectSetName((PetscObject) idm,  name);
    1366:     DMPlexCopyCoordinates(dm, idm);
    1367:     DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE);
    1368:     PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);
    1369:     if (flg) {DMPlexOrientInterface_Internal(idm);}
    1370:   }
    1371:   {
    1372:     PetscBool            isper;
    1373:     const PetscReal      *maxCell, *L;
    1374:     const DMBoundaryType *bd;

    1376:     DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
    1377:     DMSetPeriodicity(idm,isper,maxCell,L,bd);
    1378:   }
    1379:   /* This function makes the mesh fully interpolated on all ranks */
    1380:   {
    1381:     DM_Plex *plex = (DM_Plex *) idm->data;
    1382:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL;
    1383:   }
    1384:   *dmInt = idm;
    1385:   PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);
    1386:   return(0);
    1387: }

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

    1392:   Collective on dmA

    1394:   Input Parameter:
    1395: . dmA - The DMPlex object with initial coordinates

    1397:   Output Parameter:
    1398: . dmB - The DMPlex object with copied coordinates

    1400:   Level: intermediate

    1402:   Note: This is typically used when adding pieces other than vertices to a mesh

    1404: .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
    1405: @*/
    1406: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
    1407: {
    1408:   Vec            coordinatesA, coordinatesB;
    1409:   VecType        vtype;
    1410:   PetscSection   coordSectionA, coordSectionB;
    1411:   PetscScalar   *coordsA, *coordsB;
    1412:   PetscInt       spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
    1413:   PetscInt       cStartA, cEndA, cStartB, cEndB, cS, cE, cdim;
    1414:   PetscBool      lc = PETSC_FALSE;

    1420:   if (dmA == dmB) return(0);
    1421:   DMGetCoordinateDim(dmA, &cdim);
    1422:   DMSetCoordinateDim(dmB, cdim);
    1423:   DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
    1424:   DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
    1425:   if ((vEndA-vStartA) != (vEndB-vStartB)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of vertices in first DM %d != %d in the second DM", vEndA-vStartA, vEndB-vStartB);
    1426:   DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);
    1427:   DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);
    1428:   DMGetCoordinateSection(dmA, &coordSectionA);
    1429:   DMGetCoordinateSection(dmB, &coordSectionB);
    1430:   if (coordSectionA == coordSectionB) return(0);
    1431:   PetscSectionGetNumFields(coordSectionA, &Nf);
    1432:   if (!Nf) return(0);
    1433:   if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf);
    1434:   if (!coordSectionB) {
    1435:     PetscInt dim;

    1437:     PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);
    1438:     DMGetCoordinateDim(dmA, &dim);
    1439:     DMSetCoordinateSection(dmB, dim, coordSectionB);
    1440:     PetscObjectDereference((PetscObject) coordSectionB);
    1441:   }
    1442:   PetscSectionSetNumFields(coordSectionB, 1);
    1443:   PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
    1444:   PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
    1445:   PetscSectionGetChart(coordSectionA, &cS, &cE);
    1446:   if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
    1447:     if ((cEndA-cStartA) != (cEndB-cStartB)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cells in first DM %D != %D in the second DM", cEndA-cStartA, cEndB-cStartB);
    1448:     cS = cS - cStartA + cStartB;
    1449:     cE = vEndB;
    1450:     lc = PETSC_TRUE;
    1451:   } else {
    1452:     cS = vStartB;
    1453:     cE = vEndB;
    1454:   }
    1455:   PetscSectionSetChart(coordSectionB, cS, cE);
    1456:   for (v = vStartB; v < vEndB; ++v) {
    1457:     PetscSectionSetDof(coordSectionB, v, spaceDim);
    1458:     PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
    1459:   }
    1460:   if (lc) { /* localized coordinates */
    1461:     PetscInt c;

    1463:     for (c = cS-cStartB; c < cEndB-cStartB; c++) {
    1464:       PetscInt dof;

    1466:       PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
    1467:       PetscSectionSetDof(coordSectionB, c + cStartB, dof);
    1468:       PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);
    1469:     }
    1470:   }
    1471:   PetscSectionSetUp(coordSectionB);
    1472:   PetscSectionGetStorageSize(coordSectionB, &coordSizeB);
    1473:   DMGetCoordinatesLocal(dmA, &coordinatesA);
    1474:   VecCreate(PETSC_COMM_SELF, &coordinatesB);
    1475:   PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
    1476:   VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);
    1477:   VecGetBlockSize(coordinatesA, &d);
    1478:   VecSetBlockSize(coordinatesB, d);
    1479:   VecGetType(coordinatesA, &vtype);
    1480:   VecSetType(coordinatesB, vtype);
    1481:   VecGetArray(coordinatesA, &coordsA);
    1482:   VecGetArray(coordinatesB, &coordsB);
    1483:   for (v = 0; v < vEndB-vStartB; ++v) {
    1484:     PetscInt offA, offB;

    1486:     PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);
    1487:     PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);
    1488:     for (d = 0; d < spaceDim; ++d) {
    1489:       coordsB[offB+d] = coordsA[offA+d];
    1490:     }
    1491:   }
    1492:   if (lc) { /* localized coordinates */
    1493:     PetscInt c;

    1495:     for (c = cS-cStartB; c < cEndB-cStartB; c++) {
    1496:       PetscInt dof, offA, offB;

    1498:       PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);
    1499:       PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);
    1500:       PetscSectionGetDof(coordSectionA, c + cStartA, &dof);
    1501:       PetscArraycpy(coordsB + offB,coordsA + offA,dof);
    1502:     }
    1503:   }
    1504:   VecRestoreArray(coordinatesA, &coordsA);
    1505:   VecRestoreArray(coordinatesB, &coordsB);
    1506:   DMSetCoordinatesLocal(dmB, coordinatesB);
    1507:   VecDestroy(&coordinatesB);
    1508:   return(0);
    1509: }

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

    1514:   Collective on dm

    1516:   Input Parameter:
    1517: . dm - The complete DMPlex object

    1519:   Output Parameter:
    1520: . dmUnint - The DMPlex object with only cells and vertices

    1522:   Level: intermediate

    1524:   Notes:
    1525:     It does not copy over the coordinates.

    1527:   Developer Notes:
    1528:     It sets plex->interpolated = DMPLEX_INTERPOLATED_NONE.

    1530: .seealso: DMPlexInterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates()
    1531: @*/
    1532: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
    1533: {
    1534:   DMPlexInterpolatedFlag interpolated;
    1535:   DM             udm;
    1536:   PetscInt       dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;

    1542:   DMGetDimension(dm, &dim);
    1543:   DMPlexIsInterpolated(dm, &interpolated);
    1544:   if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
    1545:   if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) {
    1546:     /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */
    1547:     PetscObjectReference((PetscObject) dm);
    1548:     *dmUnint = dm;
    1549:     return(0);
    1550:   }
    1551:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
    1552:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
    1553:   DMCreate(PetscObjectComm((PetscObject) dm), &udm);
    1554:   DMSetType(udm, DMPLEX);
    1555:   DMSetDimension(udm, dim);
    1556:   DMPlexSetChart(udm, cStart, vEnd);
    1557:   for (c = cStart; c < cEnd; ++c) {
    1558:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

    1560:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
    1561:     for (cl = 0; cl < closureSize*2; cl += 2) {
    1562:       const PetscInt p = closure[cl];

    1564:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
    1565:     }
    1566:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
    1567:     DMPlexSetConeSize(udm, c, coneSize);
    1568:     maxConeSize = PetscMax(maxConeSize, coneSize);
    1569:   }
    1570:   DMSetUp(udm);
    1571:   PetscMalloc1(maxConeSize, &cone);
    1572:   for (c = cStart; c < cEnd; ++c) {
    1573:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

    1575:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
    1576:     for (cl = 0; cl < closureSize*2; cl += 2) {
    1577:       const PetscInt p = closure[cl];

    1579:       if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
    1580:     }
    1581:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
    1582:     DMPlexSetCone(udm, c, cone);
    1583:   }
    1584:   PetscFree(cone);
    1585:   DMPlexSymmetrize(udm);
    1586:   DMPlexStratify(udm);
    1587:   /* Reduce SF */
    1588:   {
    1589:     PetscSF            sfPoint, sfPointUn;
    1590:     const PetscSFNode *remotePoints;
    1591:     const PetscInt    *localPoints;
    1592:     PetscSFNode       *remotePointsUn;
    1593:     PetscInt          *localPointsUn;
    1594:     PetscInt           vEnd, numRoots, numLeaves, l;
    1595:     PetscInt           numLeavesUn = 0, n = 0;
    1596:     PetscErrorCode     ierr;

    1598:     /* Get original SF information */
    1599:     DMGetPointSF(dm, &sfPoint);
    1600:     DMGetPointSF(udm, &sfPointUn);
    1601:     DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);
    1602:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
    1603:     /* Allocate space for cells and vertices */
    1604:     for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++;
    1605:     /* Fill in leaves */
    1606:     if (vEnd >= 0) {
    1607:       PetscMalloc1(numLeavesUn, &remotePointsUn);
    1608:       PetscMalloc1(numLeavesUn, &localPointsUn);
    1609:       for (l = 0; l < numLeaves; l++) {
    1610:         if (localPoints[l] < vEnd) {
    1611:           localPointsUn[n]        = localPoints[l];
    1612:           remotePointsUn[n].rank  = remotePoints[l].rank;
    1613:           remotePointsUn[n].index = remotePoints[l].index;
    1614:           ++n;
    1615:         }
    1616:       }
    1617:       if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn);
    1618:       PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);
    1619:     }
    1620:   }
    1621:   {
    1622:     PetscBool            isper;
    1623:     const PetscReal      *maxCell, *L;
    1624:     const DMBoundaryType *bd;

    1626:     DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);
    1627:     DMSetPeriodicity(udm,isper,maxCell,L,bd);
    1628:   }
    1629:   /* This function makes the mesh fully uninterpolated on all ranks */
    1630:   {
    1631:     DM_Plex *plex = (DM_Plex *) udm->data;
    1632:     plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE;
    1633:   }
    1634:   *dmUnint = udm;
    1635:   return(0);
    1636: }

    1638: static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated)
    1639: {
    1640:   PetscInt       coneSize, depth, dim, h, p, pStart, pEnd;
    1641:   MPI_Comm       comm;

    1645:   PetscObjectGetComm((PetscObject)dm, &comm);
    1646:   DMPlexGetDepth(dm, &depth);
    1647:   DMGetDimension(dm, &dim);

    1649:   if (depth == dim) {
    1650:     *interpolated = DMPLEX_INTERPOLATED_FULL;
    1651:     if (!dim) goto finish;

    1653:     /* Check points at height = dim are vertices (have no cones) */
    1654:     DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd);
    1655:     for (p=pStart; p<pEnd; p++) {
    1656:       DMPlexGetConeSize(dm, p, &coneSize);
    1657:       if (coneSize) {
    1658:         *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
    1659:         goto finish;
    1660:       }
    1661:     }

    1663:     /* Check points at height < dim have cones */
    1664:     for (h=0; h<dim; h++) {
    1665:       DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
    1666:       for (p=pStart; p<pEnd; p++) {
    1667:         DMPlexGetConeSize(dm, p, &coneSize);
    1668:         if (!coneSize) {
    1669:           *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
    1670:           goto finish;
    1671:         }
    1672:       }
    1673:     }
    1674:   } else if (depth == 1) {
    1675:     *interpolated = DMPLEX_INTERPOLATED_NONE;
    1676:   } else {
    1677:     *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
    1678:   }
    1679: finish:
    1680:   return(0);
    1681: }

    1683: /*@
    1684:   DMPlexIsInterpolated - Find out to what extent the DMPlex is topologically interpolated.

    1686:   Not Collective

    1688:   Input Parameter:
    1689: . dm      - The DM object

    1691:   Output Parameter:
    1692: . interpolated - Flag whether the DM is interpolated

    1694:   Level: intermediate

    1696:   Notes:
    1697:   Unlike DMPlexIsInterpolatedCollective(), this is NOT collective
    1698:   so the results can be different on different ranks in special cases.
    1699:   However, DMPlexInterpolate() guarantees the result is the same on all.

    1701:   Unlike DMPlexIsInterpolatedCollective(), this cannot return DMPLEX_INTERPOLATED_MIXED.

    1703:   Developer Notes:
    1704:   Initially, plex->interpolated = DMPLEX_INTERPOLATED_INVALID.

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

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

    1712:   DMPlexInterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_FULL,
    1713:   and DMPlexUninterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_NONE.

    1715: .seealso: DMPlexInterpolate(), DMPlexIsInterpolatedCollective()
    1716: @*/
    1717: PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated)
    1718: {
    1719:   DM_Plex        *plex = (DM_Plex *) dm->data;
    1720:   PetscErrorCode  ierr;

    1725:   if (plex->interpolated < 0) {
    1726:     DMPlexIsInterpolated_Internal(dm, &plex->interpolated);
    1727:   } else if (PetscDefined (USE_DEBUG)) {
    1728:     DMPlexInterpolatedFlag flg;

    1730:     DMPlexIsInterpolated_Internal(dm, &flg);
    1731:     if (flg != plex->interpolated) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stashed DMPlexInterpolatedFlag %s is inconsistent with current %s", DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[flg]);
    1732:   }
    1733:   *interpolated = plex->interpolated;
    1734:   return(0);
    1735: }

    1737: /*@
    1738:   DMPlexIsInterpolatedCollective - Find out to what extent the DMPlex is topologically interpolated (in collective manner).

    1740:   Collective

    1742:   Input Parameter:
    1743: . dm      - The DM object

    1745:   Output Parameter:
    1746: . interpolated - Flag whether the DM is interpolated

    1748:   Level: intermediate

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

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

    1755:   Developer Notes:
    1756:   Initially, plex->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID.

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

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

    1765: .seealso: DMPlexInterpolate(), DMPlexIsInterpolated()
    1766: @*/
    1767: PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated)
    1768: {
    1769:   DM_Plex        *plex = (DM_Plex *) dm->data;
    1770:   PetscBool       debug=PETSC_FALSE;
    1771:   PetscErrorCode  ierr;

    1776:   PetscOptionsGetBool(((PetscObject) dm)->options, ((PetscObject) dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL);
    1777:   if (plex->interpolatedCollective < 0) {
    1778:     DMPlexInterpolatedFlag  min, max;
    1779:     MPI_Comm                comm;

    1781:     PetscObjectGetComm((PetscObject)dm, &comm);
    1782:     DMPlexIsInterpolated(dm, &plex->interpolatedCollective);
    1783:     MPI_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm);
    1784:     MPI_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm);
    1785:     if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED;
    1786:     if (debug) {
    1787:       PetscMPIInt rank;

    1789:       MPI_Comm_rank(comm, &rank);
    1790:       PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]);
    1791:       PetscSynchronizedFlush(comm, PETSC_STDOUT);
    1792:     }
    1793:   }
    1794:   *interpolated = plex->interpolatedCollective;
    1795:   return(0);
    1796: }