Actual source code: plexpartition.c

      1: #include <petsc/private/dmpleximpl.h>
      2: #include <petsc/private/partitionerimpl.h>
      3: #include <petsc/private/hashseti.h>

      5: const char *const DMPlexCSRAlgorithms[] = {"mat", "graph", "overlap", "DMPlexCSRAlgorithm", "DM_PLEX_CSR_", NULL};

      7: static inline PetscInt DMPlex_GlobalID(PetscInt point)
      8: {
      9:   return point >= 0 ? point : -(point + 1);
     10: }

     12: static PetscErrorCode DMPlexCreatePartitionerGraph_Overlap(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
     13: {
     14:   DM              ovdm;
     15:   PetscSF         sfPoint;
     16:   IS              cellNumbering;
     17:   const PetscInt *cellNum;
     18:   PetscInt       *adj = NULL, *vOffsets = NULL, *vAdj = NULL;
     19:   PetscBool       useCone, useClosure;
     20:   PetscInt        dim, depth, overlap, cStart, cEnd, c, v;
     21:   PetscMPIInt     rank, size;

     23:   PetscFunctionBegin;
     24:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
     25:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
     26:   PetscCall(DMGetDimension(dm, &dim));
     27:   PetscCall(DMPlexGetDepth(dm, &depth));
     28:   if (dim != depth) {
     29:     /* We do not handle the uninterpolated case here */
     30:     PetscCall(DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency));
     31:     /* DMPlexCreateNeighborCSR does not make a numbering */
     32:     if (globalNumbering) PetscCall(DMPlexCreateCellNumbering(dm, PETSC_TRUE, globalNumbering));
     33:     /* Different behavior for empty graphs */
     34:     if (!*numVertices) {
     35:       PetscCall(PetscMalloc1(1, offsets));
     36:       (*offsets)[0] = 0;
     37:     }
     38:     /* Broken in parallel */
     39:     if (rank) PetscCheck(!*numVertices, PETSC_COMM_SELF, PETSC_ERR_SUP, "Parallel partitioning of uninterpolated meshes not supported");
     40:     PetscFunctionReturn(PETSC_SUCCESS);
     41:   }
     42:   /* Always use FVM adjacency to create partitioner graph */
     43:   PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure));
     44:   PetscCall(DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE));
     45:   /* Need overlap >= 1 */
     46:   PetscCall(DMPlexGetOverlap(dm, &overlap));
     47:   if (size && overlap < 1) {
     48:     PetscCall(DMPlexDistributeOverlap(dm, 1, NULL, &ovdm));
     49:   } else {
     50:     PetscCall(PetscObjectReference((PetscObject)dm));
     51:     ovdm = dm;
     52:   }
     53:   PetscCall(DMGetPointSF(ovdm, &sfPoint));
     54:   PetscCall(DMPlexGetHeightStratum(ovdm, height, &cStart, &cEnd));
     55:   PetscCall(DMPlexCreateNumbering_Plex(ovdm, cStart, cEnd, 0, NULL, sfPoint, &cellNumbering));
     56:   if (globalNumbering) {
     57:     PetscCall(PetscObjectReference((PetscObject)cellNumbering));
     58:     *globalNumbering = cellNumbering;
     59:   }
     60:   PetscCall(ISGetIndices(cellNumbering, &cellNum));
     61:   /* Determine sizes */
     62:   for (*numVertices = 0, c = cStart; c < cEnd; ++c) {
     63:     /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
     64:     if (cellNum[c - cStart] < 0) continue;
     65:     (*numVertices)++;
     66:   }
     67:   PetscCall(PetscCalloc1(*numVertices + 1, &vOffsets));
     68:   for (c = cStart, v = 0; c < cEnd; ++c) {
     69:     PetscInt adjSize = PETSC_DETERMINE, a, vsize = 0;

     71:     if (cellNum[c - cStart] < 0) continue;
     72:     PetscCall(DMPlexGetAdjacency(ovdm, c, &adjSize, &adj));
     73:     for (a = 0; a < adjSize; ++a) {
     74:       const PetscInt point = adj[a];
     75:       if (point != c && cStart <= point && point < cEnd) ++vsize;
     76:     }
     77:     vOffsets[v + 1] = vOffsets[v] + vsize;
     78:     ++v;
     79:   }
     80:   /* Determine adjacency */
     81:   PetscCall(PetscMalloc1(vOffsets[*numVertices], &vAdj));
     82:   for (c = cStart, v = 0; c < cEnd; ++c) {
     83:     PetscInt adjSize = PETSC_DETERMINE, a, off = vOffsets[v];

     85:     /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
     86:     if (cellNum[c - cStart] < 0) continue;
     87:     PetscCall(DMPlexGetAdjacency(ovdm, c, &adjSize, &adj));
     88:     for (a = 0; a < adjSize; ++a) {
     89:       const PetscInt point = adj[a];
     90:       if (point != c && cStart <= point && point < cEnd) vAdj[off++] = DMPlex_GlobalID(cellNum[point - cStart]);
     91:     }
     92:     PetscCheck(off == vOffsets[v + 1], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Offsets %" PetscInt_FMT " should be %" PetscInt_FMT, off, vOffsets[v + 1]);
     93:     /* Sort adjacencies (not strictly necessary) */
     94:     PetscCall(PetscSortInt(off - vOffsets[v], &vAdj[vOffsets[v]]));
     95:     ++v;
     96:   }
     97:   PetscCall(PetscFree(adj));
     98:   PetscCall(ISRestoreIndices(cellNumbering, &cellNum));
     99:   PetscCall(ISDestroy(&cellNumbering));
    100:   PetscCall(DMSetBasicAdjacency(dm, useCone, useClosure));
    101:   PetscCall(DMDestroy(&ovdm));
    102:   if (offsets) {
    103:     *offsets = vOffsets;
    104:   } else PetscCall(PetscFree(vOffsets));
    105:   if (adjacency) {
    106:     *adjacency = vAdj;
    107:   } else PetscCall(PetscFree(vAdj));
    108:   PetscFunctionReturn(PETSC_SUCCESS);
    109: }

    111: static PetscErrorCode DMPlexCreatePartitionerGraph_Native(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
    112: {
    113:   PetscInt        dim, depth, p, pStart, pEnd, a, adjSize, idx, size;
    114:   PetscInt       *adj = NULL, *vOffsets = NULL, *graph = NULL;
    115:   IS              cellNumbering;
    116:   const PetscInt *cellNum;
    117:   PetscBool       useCone, useClosure;
    118:   PetscSection    section;
    119:   PetscSegBuffer  adjBuffer;
    120:   PetscSF         sfPoint;
    121:   PetscInt       *adjCells = NULL, *remoteCells = NULL;
    122:   const PetscInt *local;
    123:   PetscInt        nroots, nleaves, l;
    124:   PetscMPIInt     rank;

    126:   PetscFunctionBegin;
    127:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    128:   PetscCall(DMGetDimension(dm, &dim));
    129:   PetscCall(DMPlexGetDepth(dm, &depth));
    130:   if (dim != depth) {
    131:     /* We do not handle the uninterpolated case here */
    132:     PetscCall(DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency));
    133:     /* DMPlexCreateNeighborCSR does not make a numbering */
    134:     if (globalNumbering) PetscCall(DMPlexCreateCellNumbering(dm, PETSC_TRUE, globalNumbering));
    135:     /* Different behavior for empty graphs */
    136:     if (!*numVertices) {
    137:       PetscCall(PetscMalloc1(1, offsets));
    138:       (*offsets)[0] = 0;
    139:     }
    140:     /* Broken in parallel */
    141:     if (rank) PetscCheck(!*numVertices, PETSC_COMM_SELF, PETSC_ERR_SUP, "Parallel partitioning of uninterpolated meshes not supported");
    142:     PetscFunctionReturn(PETSC_SUCCESS);
    143:   }
    144:   PetscCall(DMGetPointSF(dm, &sfPoint));
    145:   PetscCall(DMPlexGetHeightStratum(dm, height, &pStart, &pEnd));
    146:   /* Build adjacency graph via a section/segbuffer */
    147:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section));
    148:   PetscCall(PetscSectionSetChart(section, pStart, pEnd));
    149:   PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 1000, &adjBuffer));
    150:   /* Always use FVM adjacency to create partitioner graph */
    151:   PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure));
    152:   PetscCall(DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE));
    153:   PetscCall(DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, sfPoint, &cellNumbering));
    154:   if (globalNumbering) {
    155:     PetscCall(PetscObjectReference((PetscObject)cellNumbering));
    156:     *globalNumbering = cellNumbering;
    157:   }
    158:   PetscCall(ISGetIndices(cellNumbering, &cellNum));
    159:   /* For all boundary faces (including faces adjacent to a ghost cell), record the local cell in adjCells
    160:      Broadcast adjCells to remoteCells (to get cells from roots) and Reduce adjCells to remoteCells (to get cells from leaves)
    161:    */
    162:   PetscCall(PetscSFGetGraph(sfPoint, &nroots, &nleaves, &local, NULL));
    163:   if (nroots >= 0) {
    164:     PetscInt fStart, fEnd, f;

    166:     PetscCall(PetscCalloc2(nroots, &adjCells, nroots, &remoteCells));
    167:     PetscCall(DMPlexGetHeightStratum(dm, height + 1, &fStart, &fEnd));
    168:     for (l = 0; l < nroots; ++l) adjCells[l] = -3;
    169:     for (f = fStart; f < fEnd; ++f) {
    170:       const PetscInt *support;
    171:       PetscInt        supportSize;

    173:       PetscCall(DMPlexGetSupport(dm, f, &support));
    174:       PetscCall(DMPlexGetSupportSize(dm, f, &supportSize));
    175:       if (supportSize == 1) adjCells[f] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
    176:       else if (supportSize == 2) {
    177:         PetscCall(PetscFindInt(support[0], nleaves, local, &p));
    178:         if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[1] - pStart]);
    179:         PetscCall(PetscFindInt(support[1], nleaves, local, &p));
    180:         if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
    181:       }
    182:       /* Handle non-conforming meshes */
    183:       if (supportSize > 2) {
    184:         PetscInt        numChildren, i;
    185:         const PetscInt *children;

    187:         PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, &children));
    188:         for (i = 0; i < numChildren; ++i) {
    189:           const PetscInt child = children[i];
    190:           if (fStart <= child && child < fEnd) {
    191:             PetscCall(DMPlexGetSupport(dm, child, &support));
    192:             PetscCall(DMPlexGetSupportSize(dm, child, &supportSize));
    193:             if (supportSize == 1) adjCells[child] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
    194:             else if (supportSize == 2) {
    195:               PetscCall(PetscFindInt(support[0], nleaves, local, &p));
    196:               if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[1] - pStart]);
    197:               PetscCall(PetscFindInt(support[1], nleaves, local, &p));
    198:               if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
    199:             }
    200:           }
    201:         }
    202:       }
    203:     }
    204:     for (l = 0; l < nroots; ++l) remoteCells[l] = -1;
    205:     PetscCall(PetscSFBcastBegin(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_REPLACE));
    206:     PetscCall(PetscSFBcastEnd(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_REPLACE));
    207:     PetscCall(PetscSFReduceBegin(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX));
    208:     PetscCall(PetscSFReduceEnd(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX));
    209:   }
    210:   /* Combine local and global adjacencies */
    211:   for (*numVertices = 0, p = pStart; p < pEnd; p++) {
    212:     /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
    213:     if (nroots > 0) {
    214:       if (cellNum[p - pStart] < 0) continue;
    215:     }
    216:     /* Add remote cells */
    217:     if (remoteCells) {
    218:       const PetscInt  gp = DMPlex_GlobalID(cellNum[p - pStart]);
    219:       PetscInt        coneSize, numChildren, c, i;
    220:       const PetscInt *cone, *children;

    222:       PetscCall(DMPlexGetCone(dm, p, &cone));
    223:       PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
    224:       for (c = 0; c < coneSize; ++c) {
    225:         const PetscInt point = cone[c];
    226:         if (remoteCells[point] >= 0 && remoteCells[point] != gp) {
    227:           PetscInt *PETSC_RESTRICT pBuf;
    228:           PetscCall(PetscSectionAddDof(section, p, 1));
    229:           PetscCall(PetscSegBufferGetInts(adjBuffer, 1, &pBuf));
    230:           *pBuf = remoteCells[point];
    231:         }
    232:         /* Handle non-conforming meshes */
    233:         PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, &children));
    234:         for (i = 0; i < numChildren; ++i) {
    235:           const PetscInt child = children[i];
    236:           if (remoteCells[child] >= 0 && remoteCells[child] != gp) {
    237:             PetscInt *PETSC_RESTRICT pBuf;
    238:             PetscCall(PetscSectionAddDof(section, p, 1));
    239:             PetscCall(PetscSegBufferGetInts(adjBuffer, 1, &pBuf));
    240:             *pBuf = remoteCells[child];
    241:           }
    242:         }
    243:       }
    244:     }
    245:     /* Add local cells */
    246:     adjSize = PETSC_DETERMINE;
    247:     PetscCall(DMPlexGetAdjacency(dm, p, &adjSize, &adj));
    248:     for (a = 0; a < adjSize; ++a) {
    249:       const PetscInt point = adj[a];
    250:       if (point != p && pStart <= point && point < pEnd) {
    251:         PetscInt *PETSC_RESTRICT pBuf;
    252:         PetscCall(PetscSectionAddDof(section, p, 1));
    253:         PetscCall(PetscSegBufferGetInts(adjBuffer, 1, &pBuf));
    254:         *pBuf = DMPlex_GlobalID(cellNum[point - pStart]);
    255:       }
    256:     }
    257:     (*numVertices)++;
    258:   }
    259:   PetscCall(PetscFree(adj));
    260:   PetscCall(PetscFree2(adjCells, remoteCells));
    261:   PetscCall(DMSetBasicAdjacency(dm, useCone, useClosure));

    263:   /* Derive CSR graph from section/segbuffer */
    264:   PetscCall(PetscSectionSetUp(section));
    265:   PetscCall(PetscSectionGetStorageSize(section, &size));
    266:   PetscCall(PetscMalloc1(*numVertices + 1, &vOffsets));
    267:   for (idx = 0, p = pStart; p < pEnd; p++) {
    268:     if (nroots > 0) {
    269:       if (cellNum[p - pStart] < 0) continue;
    270:     }
    271:     PetscCall(PetscSectionGetOffset(section, p, &vOffsets[idx++]));
    272:   }
    273:   vOffsets[*numVertices] = size;
    274:   PetscCall(PetscSegBufferExtractAlloc(adjBuffer, &graph));

    276:   if (nroots >= 0) {
    277:     /* Filter out duplicate edges using section/segbuffer */
    278:     PetscCall(PetscSectionReset(section));
    279:     PetscCall(PetscSectionSetChart(section, 0, *numVertices));
    280:     for (p = 0; p < *numVertices; p++) {
    281:       PetscInt start = vOffsets[p], end = vOffsets[p + 1];
    282:       PetscInt numEdges = end - start, *PETSC_RESTRICT edges;
    283:       PetscCall(PetscSortRemoveDupsInt(&numEdges, &graph[start]));
    284:       PetscCall(PetscSectionSetDof(section, p, numEdges));
    285:       PetscCall(PetscSegBufferGetInts(adjBuffer, numEdges, &edges));
    286:       PetscCall(PetscArraycpy(edges, &graph[start], numEdges));
    287:     }
    288:     PetscCall(PetscFree(vOffsets));
    289:     PetscCall(PetscFree(graph));
    290:     /* Derive CSR graph from section/segbuffer */
    291:     PetscCall(PetscSectionSetUp(section));
    292:     PetscCall(PetscSectionGetStorageSize(section, &size));
    293:     PetscCall(PetscMalloc1(*numVertices + 1, &vOffsets));
    294:     for (idx = 0, p = 0; p < *numVertices; p++) PetscCall(PetscSectionGetOffset(section, p, &vOffsets[idx++]));
    295:     vOffsets[*numVertices] = size;
    296:     PetscCall(PetscSegBufferExtractAlloc(adjBuffer, &graph));
    297:   } else {
    298:     /* Sort adjacencies (not strictly necessary) */
    299:     for (p = 0; p < *numVertices; p++) {
    300:       PetscInt start = vOffsets[p], end = vOffsets[p + 1];
    301:       PetscCall(PetscSortInt(end - start, &graph[start]));
    302:     }
    303:   }

    305:   if (offsets) *offsets = vOffsets;
    306:   if (adjacency) *adjacency = graph;

    308:   /* Cleanup */
    309:   PetscCall(PetscSegBufferDestroy(&adjBuffer));
    310:   PetscCall(PetscSectionDestroy(&section));
    311:   PetscCall(ISRestoreIndices(cellNumbering, &cellNum));
    312:   PetscCall(ISDestroy(&cellNumbering));
    313:   PetscFunctionReturn(PETSC_SUCCESS);
    314: }

    316: static PetscErrorCode DMPlexCreatePartitionerGraph_ViaMat(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
    317: {
    318:   Mat             conn, CSR;
    319:   IS              fis, cis, cis_own;
    320:   PetscSF         sfPoint;
    321:   const PetscInt *rows, *cols, *ii, *jj;
    322:   PetscInt       *idxs, *idxs2;
    323:   PetscInt        dim, depth, floc, cloc, i, M, N, c, lm, m, cStart, cEnd, fStart, fEnd;
    324:   PetscMPIInt     rank;
    325:   PetscBool       flg;

    327:   PetscFunctionBegin;
    328:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    329:   PetscCall(DMGetDimension(dm, &dim));
    330:   PetscCall(DMPlexGetDepth(dm, &depth));
    331:   if (dim != depth) {
    332:     /* We do not handle the uninterpolated case here */
    333:     PetscCall(DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency));
    334:     /* DMPlexCreateNeighborCSR does not make a numbering */
    335:     if (globalNumbering) PetscCall(DMPlexCreateCellNumbering(dm, PETSC_TRUE, globalNumbering));
    336:     /* Different behavior for empty graphs */
    337:     if (!*numVertices) {
    338:       PetscCall(PetscMalloc1(1, offsets));
    339:       (*offsets)[0] = 0;
    340:     }
    341:     /* Broken in parallel */
    342:     if (rank) PetscCheck(!*numVertices, PETSC_COMM_SELF, PETSC_ERR_SUP, "Parallel partitioning of uninterpolated meshes not supported");
    343:     PetscFunctionReturn(PETSC_SUCCESS);
    344:   }
    345:   /* Interpolated and parallel case */

    347:   /* numbering */
    348:   PetscCall(DMGetPointSF(dm, &sfPoint));
    349:   PetscCall(DMPlexGetHeightStratum(dm, height, &cStart, &cEnd));
    350:   PetscCall(DMPlexGetHeightStratum(dm, height + 1, &fStart, &fEnd));
    351:   PetscCall(DMPlexCreateNumbering_Plex(dm, cStart, cEnd, 0, &N, sfPoint, &cis));
    352:   PetscCall(DMPlexCreateNumbering_Plex(dm, fStart, fEnd, 0, &M, sfPoint, &fis));
    353:   if (globalNumbering) PetscCall(ISDuplicate(cis, globalNumbering));

    355:   /* get positive global ids and local sizes for facets and cells */
    356:   PetscCall(ISGetLocalSize(fis, &m));
    357:   PetscCall(ISGetIndices(fis, &rows));
    358:   PetscCall(PetscMalloc1(m, &idxs));
    359:   for (i = 0, floc = 0; i < m; i++) {
    360:     const PetscInt p = rows[i];

    362:     if (p < 0) {
    363:       idxs[i] = -(p + 1);
    364:     } else {
    365:       idxs[i] = p;
    366:       floc += 1;
    367:     }
    368:   }
    369:   PetscCall(ISRestoreIndices(fis, &rows));
    370:   PetscCall(ISDestroy(&fis));
    371:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, m, idxs, PETSC_OWN_POINTER, &fis));

    373:   PetscCall(ISGetLocalSize(cis, &m));
    374:   PetscCall(ISGetIndices(cis, &cols));
    375:   PetscCall(PetscMalloc1(m, &idxs));
    376:   PetscCall(PetscMalloc1(m, &idxs2));
    377:   for (i = 0, cloc = 0; i < m; i++) {
    378:     const PetscInt p = cols[i];

    380:     if (p < 0) {
    381:       idxs[i] = -(p + 1);
    382:     } else {
    383:       idxs[i]       = p;
    384:       idxs2[cloc++] = p;
    385:     }
    386:   }
    387:   PetscCall(ISRestoreIndices(cis, &cols));
    388:   PetscCall(ISDestroy(&cis));
    389:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), m, idxs, PETSC_OWN_POINTER, &cis));
    390:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), cloc, idxs2, PETSC_OWN_POINTER, &cis_own));

    392:   /* Create matrix to hold F-C connectivity (MatMatTranspose Mult not supported for MPIAIJ) */
    393:   PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), &conn));
    394:   PetscCall(MatSetSizes(conn, floc, cloc, M, N));
    395:   PetscCall(MatSetType(conn, MATMPIAIJ));
    396:   PetscCall(DMPlexGetMaxSizes(dm, NULL, &lm));
    397:   PetscCallMPI(MPIU_Allreduce(&lm, &m, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm)));
    398:   PetscCall(MatMPIAIJSetPreallocation(conn, m, NULL, m, NULL));

    400:   /* Assemble matrix */
    401:   PetscCall(ISGetIndices(fis, &rows));
    402:   PetscCall(ISGetIndices(cis, &cols));
    403:   for (c = cStart; c < cEnd; c++) {
    404:     const PetscInt *cone;
    405:     PetscInt        coneSize, row, col, f;

    407:     col = cols[c - cStart];
    408:     PetscCall(DMPlexGetCone(dm, c, &cone));
    409:     PetscCall(DMPlexGetConeSize(dm, c, &coneSize));
    410:     for (f = 0; f < coneSize; f++) {
    411:       const PetscScalar v = 1.0;
    412:       const PetscInt   *children;
    413:       PetscInt          numChildren, ch;

    415:       row = rows[cone[f] - fStart];
    416:       PetscCall(MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES));

    418:       /* non-conforming meshes */
    419:       PetscCall(DMPlexGetTreeChildren(dm, cone[f], &numChildren, &children));
    420:       for (ch = 0; ch < numChildren; ch++) {
    421:         const PetscInt child = children[ch];

    423:         if (child < fStart || child >= fEnd) continue;
    424:         row = rows[child - fStart];
    425:         PetscCall(MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES));
    426:       }
    427:     }
    428:   }
    429:   PetscCall(ISRestoreIndices(fis, &rows));
    430:   PetscCall(ISRestoreIndices(cis, &cols));
    431:   PetscCall(ISDestroy(&fis));
    432:   PetscCall(ISDestroy(&cis));
    433:   PetscCall(MatAssemblyBegin(conn, MAT_FINAL_ASSEMBLY));
    434:   PetscCall(MatAssemblyEnd(conn, MAT_FINAL_ASSEMBLY));

    436:   /* Get parallel CSR by doing conn^T * conn */
    437:   PetscCall(MatTransposeMatMult(conn, conn, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &CSR));
    438:   PetscCall(MatDestroy(&conn));

    440:   /* extract local part of the CSR */
    441:   PetscCall(MatMPIAIJGetLocalMat(CSR, MAT_INITIAL_MATRIX, &conn));
    442:   PetscCall(MatDestroy(&CSR));
    443:   PetscCall(MatGetRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg));
    444:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No IJ format");

    446:   /* get back requested output */
    447:   if (numVertices) *numVertices = m;
    448:   if (offsets) {
    449:     PetscCall(PetscCalloc1(m + 1, &idxs));
    450:     for (i = 1; i < m + 1; i++) idxs[i] = ii[i] - i; /* ParMetis does not like self-connectivity */
    451:     *offsets = idxs;
    452:   }
    453:   if (adjacency) {
    454:     PetscCall(PetscMalloc1(ii[m] - m, &idxs));
    455:     PetscCall(ISGetIndices(cis_own, &rows));
    456:     for (i = 0, c = 0; i < m; i++) {
    457:       PetscInt j, g = rows[i];

    459:       for (j = ii[i]; j < ii[i + 1]; j++) {
    460:         if (jj[j] == g) continue; /* again, self-connectivity */
    461:         idxs[c++] = jj[j];
    462:       }
    463:     }
    464:     PetscCheck(c == ii[m] - m, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected %" PetscInt_FMT " != %" PetscInt_FMT, c, ii[m] - m);
    465:     PetscCall(ISRestoreIndices(cis_own, &rows));
    466:     *adjacency = idxs;
    467:   }

    469:   /* cleanup */
    470:   PetscCall(ISDestroy(&cis_own));
    471:   PetscCall(MatRestoreRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg));
    472:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No IJ format");
    473:   PetscCall(MatDestroy(&conn));
    474:   PetscFunctionReturn(PETSC_SUCCESS);
    475: }

    477: /*@C
    478:   DMPlexCreatePartitionerGraph - Create a CSR graph of point connections for the partitioner

    480:   Collective

    482:   Input Parameters:
    483: + dm     - The mesh `DM`
    484: - height - Height of the strata from which to construct the graph

    486:   Output Parameters:
    487: + numVertices     - Number of vertices in the graph
    488: . offsets         - Point offsets in the graph
    489: . adjacency       - Point connectivity in the graph
    490: - globalNumbering - A map from the local cell numbering to the global numbering used in "adjacency".  Negative indicates that the cell is a duplicate from another process.

    492:   Options Database Key:
    493: . -dm_plex_csr_alg <mat,graph,overlap> - Choose the algorithm for computing the CSR graph

    495:   Level: developer

    497:   Note:
    498:   The user can control the definition of adjacency for the mesh using `DMSetAdjacency()`. They should choose the combination appropriate for the function
    499:   representation on the mesh. If requested, globalNumbering needs to be destroyed by the caller; offsets and adjacency need to be freed with PetscFree().

    501: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscPartitionerGetType()`, `PetscPartitionerCreate()`, `DMSetAdjacency()`
    502: @*/
    503: PetscErrorCode DMPlexCreatePartitionerGraph(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
    504: {
    505:   DMPlexCSRAlgorithm alg = DM_PLEX_CSR_GRAPH;

    507:   PetscFunctionBegin;
    508:   PetscCall(PetscOptionsGetEnum(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_csr_alg", DMPlexCSRAlgorithms, (PetscEnum *)&alg, NULL));
    509:   switch (alg) {
    510:   case DM_PLEX_CSR_MAT:
    511:     PetscCall(DMPlexCreatePartitionerGraph_ViaMat(dm, height, numVertices, offsets, adjacency, globalNumbering));
    512:     break;
    513:   case DM_PLEX_CSR_GRAPH:
    514:     PetscCall(DMPlexCreatePartitionerGraph_Native(dm, height, numVertices, offsets, adjacency, globalNumbering));
    515:     break;
    516:   case DM_PLEX_CSR_OVERLAP:
    517:     PetscCall(DMPlexCreatePartitionerGraph_Overlap(dm, height, numVertices, offsets, adjacency, globalNumbering));
    518:     break;
    519:   }
    520:   PetscFunctionReturn(PETSC_SUCCESS);
    521: }

    523: /*@C
    524:   DMPlexCreateNeighborCSR - Create a mesh graph (cell-cell adjacency) in parallel CSR format.

    526:   Collective

    528:   Input Parameters:
    529: + dm         - The `DMPLEX`
    530: - cellHeight - The height of mesh points to treat as cells (default should be 0)

    532:   Output Parameters:
    533: + numVertices - The number of local vertices in the graph, or cells in the mesh.
    534: . offsets     - The offset to the adjacency list for each cell
    535: - adjacency   - The adjacency list for all cells

    537:   Level: advanced

    539:   Note:
    540:   This is suitable for input to a mesh partitioner like ParMetis.

    542: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreate()`
    543: @*/
    544: PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
    545: {
    546:   const PetscInt maxFaceCases = 30;
    547:   PetscInt       numFaceCases = 0;
    548:   PetscInt       numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */
    549:   PetscInt      *off, *adj;
    550:   PetscInt      *neighborCells = NULL;
    551:   PetscInt       dim, cellDim, depth = 0, faceDepth, cStart, cEnd, c, numCells, cell;

    553:   PetscFunctionBegin;
    554:   /* For parallel partitioning, I think you have to communicate supports */
    555:   PetscCall(DMGetDimension(dm, &dim));
    556:   cellDim = dim - cellHeight;
    557:   PetscCall(DMPlexGetDepth(dm, &depth));
    558:   PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
    559:   if (cEnd - cStart == 0) {
    560:     if (numVertices) *numVertices = 0;
    561:     if (offsets) *offsets = NULL;
    562:     if (adjacency) *adjacency = NULL;
    563:     PetscFunctionReturn(PETSC_SUCCESS);
    564:   }
    565:   numCells  = cEnd - cStart;
    566:   faceDepth = depth - cellHeight;
    567:   if (dim == depth) {
    568:     PetscInt f, fStart, fEnd;

    570:     PetscCall(PetscCalloc1(numCells + 1, &off));
    571:     /* Count neighboring cells */
    572:     PetscCall(DMPlexGetHeightStratum(dm, cellHeight + 1, &fStart, &fEnd));
    573:     for (f = fStart; f < fEnd; ++f) {
    574:       const PetscInt *support;
    575:       PetscInt        supportSize;
    576:       PetscCall(DMPlexGetSupportSize(dm, f, &supportSize));
    577:       PetscCall(DMPlexGetSupport(dm, f, &support));
    578:       if (supportSize == 2) {
    579:         PetscInt numChildren;

    581:         PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
    582:         if (!numChildren) {
    583:           ++off[support[0] - cStart + 1];
    584:           ++off[support[1] - cStart + 1];
    585:         }
    586:       }
    587:     }
    588:     /* Prefix sum */
    589:     for (c = 1; c <= numCells; ++c) off[c] += off[c - 1];
    590:     if (adjacency) {
    591:       PetscInt *tmp;

    593:       PetscCall(PetscMalloc1(off[numCells], &adj));
    594:       PetscCall(PetscMalloc1(numCells + 1, &tmp));
    595:       PetscCall(PetscArraycpy(tmp, off, numCells + 1));
    596:       /* Get neighboring cells */
    597:       for (f = fStart; f < fEnd; ++f) {
    598:         const PetscInt *support;
    599:         PetscInt        supportSize;
    600:         PetscCall(DMPlexGetSupportSize(dm, f, &supportSize));
    601:         PetscCall(DMPlexGetSupport(dm, f, &support));
    602:         if (supportSize == 2) {
    603:           PetscInt numChildren;

    605:           PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
    606:           if (!numChildren) {
    607:             adj[tmp[support[0] - cStart]++] = support[1];
    608:             adj[tmp[support[1] - cStart]++] = support[0];
    609:           }
    610:         }
    611:       }
    612:       for (c = 0; c < cEnd - cStart; ++c) PetscAssert(tmp[c] == off[c + 1], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Offset %" PetscInt_FMT " != %" PetscInt_FMT " for cell %" PetscInt_FMT, tmp[c], off[c], c + cStart);
    613:       PetscCall(PetscFree(tmp));
    614:     }
    615:     if (numVertices) *numVertices = numCells;
    616:     if (offsets) *offsets = off;
    617:     if (adjacency) *adjacency = adj;
    618:     PetscFunctionReturn(PETSC_SUCCESS);
    619:   }
    620:   /* Setup face recognition */
    621:   if (faceDepth == 1) {
    622:     PetscInt cornersSeen[30] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; /* Could use PetscBT */

    624:     for (c = cStart; c < cEnd; ++c) {
    625:       PetscInt corners;

    627:       PetscCall(DMPlexGetConeSize(dm, c, &corners));
    628:       if (!cornersSeen[corners]) {
    629:         PetscInt nFV;

    631:         PetscCheck(numFaceCases < maxFaceCases, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Exceeded maximum number of face recognition cases");
    632:         cornersSeen[corners] = 1;

    634:         PetscCall(DMPlexGetNumFaceVertices(dm, cellDim, corners, &nFV));

    636:         numFaceVertices[numFaceCases++] = nFV;
    637:       }
    638:     }
    639:   }
    640:   PetscCall(PetscCalloc1(numCells + 1, &off));
    641:   /* Count neighboring cells */
    642:   for (cell = cStart; cell < cEnd; ++cell) {
    643:     PetscInt numNeighbors = PETSC_DETERMINE, n;

    645:     PetscCall(DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells));
    646:     /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
    647:     for (n = 0; n < numNeighbors; ++n) {
    648:       PetscInt        cellPair[2];
    649:       PetscBool       found    = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
    650:       PetscInt        meetSize = 0;
    651:       const PetscInt *meet     = NULL;

    653:       cellPair[0] = cell;
    654:       cellPair[1] = neighborCells[n];
    655:       if (cellPair[0] == cellPair[1]) continue;
    656:       if (!found) {
    657:         PetscCall(DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet));
    658:         if (meetSize) {
    659:           PetscInt f;

    661:           for (f = 0; f < numFaceCases; ++f) {
    662:             if (numFaceVertices[f] == meetSize) {
    663:               found = PETSC_TRUE;
    664:               break;
    665:             }
    666:           }
    667:         }
    668:         PetscCall(DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet));
    669:       }
    670:       if (found) ++off[cell - cStart + 1];
    671:     }
    672:   }
    673:   /* Prefix sum */
    674:   for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell - 1];

    676:   if (adjacency) {
    677:     PetscCall(PetscMalloc1(off[numCells], &adj));
    678:     /* Get neighboring cells */
    679:     for (cell = cStart; cell < cEnd; ++cell) {
    680:       PetscInt numNeighbors = PETSC_DETERMINE, n;
    681:       PetscInt cellOffset   = 0;

    683:       PetscCall(DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells));
    684:       /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
    685:       for (n = 0; n < numNeighbors; ++n) {
    686:         PetscInt        cellPair[2];
    687:         PetscBool       found    = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
    688:         PetscInt        meetSize = 0;
    689:         const PetscInt *meet     = NULL;

    691:         cellPair[0] = cell;
    692:         cellPair[1] = neighborCells[n];
    693:         if (cellPair[0] == cellPair[1]) continue;
    694:         if (!found) {
    695:           PetscCall(DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet));
    696:           if (meetSize) {
    697:             PetscInt f;

    699:             for (f = 0; f < numFaceCases; ++f) {
    700:               if (numFaceVertices[f] == meetSize) {
    701:                 found = PETSC_TRUE;
    702:                 break;
    703:               }
    704:             }
    705:           }
    706:           PetscCall(DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet));
    707:         }
    708:         if (found) {
    709:           adj[off[cell - cStart] + cellOffset] = neighborCells[n];
    710:           ++cellOffset;
    711:         }
    712:       }
    713:     }
    714:   }
    715:   PetscCall(PetscFree(neighborCells));
    716:   if (numVertices) *numVertices = numCells;
    717:   if (offsets) *offsets = off;
    718:   if (adjacency) *adjacency = adj;
    719:   PetscFunctionReturn(PETSC_SUCCESS);
    720: }

    722: /*@
    723:   PetscPartitionerDMPlexPartition - Create a non-overlapping partition of the cells in the mesh

    725:   Collective

    727:   Input Parameters:
    728: + part          - The `PetscPartitioner`
    729: . targetSection - The `PetscSection` describing the absolute weight of each partition (can be `NULL`)
    730: - dm            - The mesh `DM`

    732:   Output Parameters:
    733: + partSection - The `PetscSection` giving the division of points by partition
    734: - partition   - The list of points by partition

    736:   Level: developer

    738:   Note:
    739:   If the `DM` has a local section associated, each point to be partitioned will be weighted by the total number of dofs identified
    740:   by the section in the transitive closure of the point.

    742: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscPartitioner`, `PetscSection`, `DMPlexDistribute()`, `PetscPartitionerCreate()`, `PetscSectionCreate()`,
    743:          `PetscSectionSetChart()`, `PetscPartitionerPartition()`
    744: @*/
    745: PetscErrorCode PetscPartitionerDMPlexPartition(PetscPartitioner part, DM dm, PetscSection targetSection, PetscSection partSection, IS *partition)
    746: {
    747:   PetscMPIInt  size;
    748:   PetscBool    isplex;
    749:   PetscSection vertSection = NULL, edgeSection = NULL;

    751:   PetscFunctionBegin;
    756:   PetscAssertPointer(partition, 5);
    757:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex));
    758:   PetscCheck(isplex, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)dm)->type_name);
    759:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)part), &size));
    760:   if (size == 1) {
    761:     PetscInt *points;
    762:     PetscInt  cStart, cEnd, c;

    764:     PetscCall(DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd));
    765:     PetscCall(PetscSectionReset(partSection));
    766:     PetscCall(PetscSectionSetChart(partSection, 0, size));
    767:     PetscCall(PetscSectionSetDof(partSection, 0, cEnd - cStart));
    768:     PetscCall(PetscSectionSetUp(partSection));
    769:     PetscCall(PetscMalloc1(cEnd - cStart, &points));
    770:     for (c = cStart; c < cEnd; ++c) points[c] = c;
    771:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), cEnd - cStart, points, PETSC_OWN_POINTER, partition));
    772:     PetscFunctionReturn(PETSC_SUCCESS);
    773:   }
    774:   if (part->height == 0) {
    775:     PetscInt  numVertices = 0;
    776:     PetscInt *start       = NULL;
    777:     PetscInt *adjacency   = NULL;
    778:     IS        globalNumbering;

    780:     if (!part->noGraph || part->viewGraph) {
    781:       PetscCall(DMPlexCreatePartitionerGraph(dm, part->height, &numVertices, &start, &adjacency, &globalNumbering));
    782:     } else { /* only compute the number of owned local vertices */
    783:       const PetscInt *idxs;
    784:       PetscInt        p, pStart, pEnd;

    786:       PetscCall(DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd));
    787:       PetscCall(DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, dm->sf, &globalNumbering));
    788:       PetscCall(ISGetIndices(globalNumbering, &idxs));
    789:       for (p = 0; p < pEnd - pStart; p++) numVertices += idxs[p] < 0 ? 0 : 1;
    790:       PetscCall(ISRestoreIndices(globalNumbering, &idxs));
    791:     }
    792:     if (part->usevwgt) {
    793:       PetscSection    section = dm->localSection, clSection = NULL;
    794:       IS              clPoints = NULL;
    795:       const PetscInt *gid, *clIdx;
    796:       PetscInt        v, p, pStart, pEnd;

    798:       /* dm->localSection encodes degrees of freedom per point, not per cell. We need to get the closure index to properly specify cell weights (aka dofs) */
    799:       /* We do this only if the local section has been set */
    800:       if (section) {
    801:         PetscCall(PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, NULL));
    802:         if (!clSection) PetscCall(DMPlexCreateClosureIndex(dm, NULL));
    803:         PetscCall(PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, &clPoints));
    804:         PetscCall(ISGetIndices(clPoints, &clIdx));
    805:       }
    806:       PetscCall(DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd));
    807:       PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &vertSection));
    808:       PetscCall(PetscSectionSetChart(vertSection, 0, numVertices));
    809:       if (globalNumbering) {
    810:         PetscCall(ISGetIndices(globalNumbering, &gid));
    811:       } else gid = NULL;
    812:       for (p = pStart, v = 0; p < pEnd; ++p) {
    813:         PetscInt dof = 1;

    815:         /* skip cells in the overlap */
    816:         if (gid && gid[p - pStart] < 0) continue;

    818:         if (section) {
    819:           PetscInt cl, clSize, clOff;

    821:           dof = 0;
    822:           PetscCall(PetscSectionGetDof(clSection, p, &clSize));
    823:           PetscCall(PetscSectionGetOffset(clSection, p, &clOff));
    824:           for (cl = 0; cl < clSize; cl += 2) {
    825:             PetscInt clDof, clPoint = clIdx[clOff + cl]; /* odd indices are reserved for orientations */

    827:             PetscCall(PetscSectionGetDof(section, clPoint, &clDof));
    828:             dof += clDof;
    829:           }
    830:         }
    831:         PetscCheck(dof, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of dofs for point %" PetscInt_FMT " in the local section should be positive", p);
    832:         PetscCall(PetscSectionSetDof(vertSection, v, dof));
    833:         v++;
    834:       }
    835:       if (globalNumbering) PetscCall(ISRestoreIndices(globalNumbering, &gid));
    836:       if (clPoints) PetscCall(ISRestoreIndices(clPoints, &clIdx));
    837:       PetscCall(PetscSectionSetUp(vertSection));
    838:     }
    839:     if (part->useewgt) {
    840:       const PetscInt numEdges = start[numVertices];

    842:       PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &edgeSection));
    843:       PetscCall(PetscSectionSetChart(edgeSection, 0, numEdges));
    844:       for (PetscInt e = 0; e < start[numVertices]; ++e) PetscCall(PetscSectionSetDof(edgeSection, e, 1));
    845:       for (PetscInt v = 0; v < numVertices; ++v) {
    846:         DMPolytopeType ct;

    848:         // Assume v is the cell number
    849:         PetscCall(DMPlexGetCellType(dm, v, &ct));
    850:         if (ct != DM_POLYTOPE_POINT_PRISM_TENSOR && ct != DM_POLYTOPE_SEG_PRISM_TENSOR && ct != DM_POLYTOPE_TRI_PRISM_TENSOR && ct != DM_POLYTOPE_QUAD_PRISM_TENSOR) continue;

    852:         for (PetscInt e = start[v]; e < start[v + 1]; ++e) PetscCall(PetscSectionSetDof(edgeSection, e, 3));
    853:       }
    854:       PetscCall(PetscSectionSetUp(edgeSection));
    855:     }
    856:     PetscCall(PetscPartitionerPartition(part, size, numVertices, start, adjacency, vertSection, edgeSection, targetSection, partSection, partition));
    857:     PetscCall(PetscFree(start));
    858:     PetscCall(PetscFree(adjacency));
    859:     if (globalNumbering) { /* partition is wrt global unique numbering: change this to be wrt local numbering */
    860:       const PetscInt *globalNum;
    861:       const PetscInt *partIdx;
    862:       PetscInt       *map, cStart, cEnd;
    863:       PetscInt       *adjusted, i, localSize, offset;
    864:       IS              newPartition;

    866:       PetscCall(ISGetLocalSize(*partition, &localSize));
    867:       PetscCall(PetscMalloc1(localSize, &adjusted));
    868:       PetscCall(ISGetIndices(globalNumbering, &globalNum));
    869:       PetscCall(ISGetIndices(*partition, &partIdx));
    870:       PetscCall(PetscMalloc1(localSize, &map));
    871:       PetscCall(DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd));
    872:       for (i = cStart, offset = 0; i < cEnd; i++) {
    873:         if (globalNum[i - cStart] >= 0) map[offset++] = i;
    874:       }
    875:       for (i = 0; i < localSize; i++) adjusted[i] = map[partIdx[i]];
    876:       PetscCall(PetscFree(map));
    877:       PetscCall(ISRestoreIndices(*partition, &partIdx));
    878:       PetscCall(ISRestoreIndices(globalNumbering, &globalNum));
    879:       PetscCall(ISCreateGeneral(PETSC_COMM_SELF, localSize, adjusted, PETSC_OWN_POINTER, &newPartition));
    880:       PetscCall(ISDestroy(&globalNumbering));
    881:       PetscCall(ISDestroy(partition));
    882:       *partition = newPartition;
    883:     }
    884:   } else SETERRQ(PetscObjectComm((PetscObject)part), PETSC_ERR_ARG_OUTOFRANGE, "Invalid height %" PetscInt_FMT " for points to partition", part->height);
    885:   PetscCall(PetscSectionDestroy(&vertSection));
    886:   PetscCall(PetscSectionDestroy(&edgeSection));
    887:   PetscFunctionReturn(PETSC_SUCCESS);
    888: }

    890: /*@
    891:   DMPlexGetPartitioner - Get the mesh partitioner

    893:   Not Collective

    895:   Input Parameter:
    896: . dm - The `DM`

    898:   Output Parameter:
    899: . part - The `PetscPartitioner`

    901:   Level: developer

    903:   Note:
    904:   This gets a borrowed reference, so the user should not destroy this `PetscPartitioner`.

    906: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscPartitioner`, `PetscSection`, `DMPlexDistribute()`, `DMPlexSetPartitioner()`, `PetscPartitionerDMPlexPartition()`, `PetscPartitionerCreate()`
    907: @*/
    908: PetscErrorCode DMPlexGetPartitioner(DM dm, PetscPartitioner *part)
    909: {
    910:   DM_Plex *mesh = (DM_Plex *)dm->data;

    912:   PetscFunctionBegin;
    914:   PetscAssertPointer(part, 2);
    915:   *part = mesh->partitioner;
    916:   PetscFunctionReturn(PETSC_SUCCESS);
    917: }

    919: /*@
    920:   DMPlexSetPartitioner - Set the mesh partitioner

    922:   logically Collective

    924:   Input Parameters:
    925: + dm   - The `DM`
    926: - part - The partitioner

    928:   Level: developer

    930:   Note:
    931:   Any existing `PetscPartitioner` will be destroyed.

    933: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscPartitioner`,`DMPlexDistribute()`, `DMPlexGetPartitioner()`, `PetscPartitionerCreate()`
    934: @*/
    935: PetscErrorCode DMPlexSetPartitioner(DM dm, PetscPartitioner part)
    936: {
    937:   DM_Plex *mesh = (DM_Plex *)dm->data;

    939:   PetscFunctionBegin;
    942:   PetscCall(PetscObjectReference((PetscObject)part));
    943:   PetscCall(PetscPartitionerDestroy(&mesh->partitioner));
    944:   mesh->partitioner = part;
    945:   PetscFunctionReturn(PETSC_SUCCESS);
    946: }

    948: static PetscErrorCode DMPlexAddClosure_Private(DM dm, PetscHSetI ht, PetscInt point)
    949: {
    950:   const PetscInt *cone;
    951:   PetscInt        coneSize, c;
    952:   PetscBool       missing;

    954:   PetscFunctionBeginHot;
    955:   PetscCall(PetscHSetIQueryAdd(ht, point, &missing));
    956:   if (missing) {
    957:     PetscCall(DMPlexGetCone(dm, point, &cone));
    958:     PetscCall(DMPlexGetConeSize(dm, point, &coneSize));
    959:     for (c = 0; c < coneSize; c++) PetscCall(DMPlexAddClosure_Private(dm, ht, cone[c]));
    960:   }
    961:   PetscFunctionReturn(PETSC_SUCCESS);
    962: }

    964: PETSC_UNUSED static PetscErrorCode DMPlexAddClosure_Tree(DM dm, PetscHSetI ht, PetscInt point, PetscBool up, PetscBool down)
    965: {
    966:   PetscFunctionBegin;
    967:   if (up) {
    968:     PetscInt parent;

    970:     PetscCall(DMPlexGetTreeParent(dm, point, &parent, NULL));
    971:     if (parent != point) {
    972:       PetscInt closureSize, *closure = NULL, i;

    974:       PetscCall(DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
    975:       for (i = 0; i < closureSize; i++) {
    976:         PetscInt cpoint = closure[2 * i];

    978:         PetscCall(PetscHSetIAdd(ht, cpoint));
    979:         PetscCall(DMPlexAddClosure_Tree(dm, ht, cpoint, PETSC_TRUE, PETSC_FALSE));
    980:       }
    981:       PetscCall(DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
    982:     }
    983:   }
    984:   if (down) {
    985:     PetscInt        numChildren;
    986:     const PetscInt *children;

    988:     PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, &children));
    989:     if (numChildren) {
    990:       PetscInt i;

    992:       for (i = 0; i < numChildren; i++) {
    993:         PetscInt cpoint = children[i];

    995:         PetscCall(PetscHSetIAdd(ht, cpoint));
    996:         PetscCall(DMPlexAddClosure_Tree(dm, ht, cpoint, PETSC_FALSE, PETSC_TRUE));
    997:       }
    998:     }
    999:   }
    1000:   PetscFunctionReturn(PETSC_SUCCESS);
    1001: }

    1003: static PetscErrorCode DMPlexAddClosureTree_Up_Private(DM dm, PetscHSetI ht, PetscInt point)
    1004: {
    1005:   PetscInt parent;

    1007:   PetscFunctionBeginHot;
    1008:   PetscCall(DMPlexGetTreeParent(dm, point, &parent, NULL));
    1009:   if (point != parent) {
    1010:     const PetscInt *cone;
    1011:     PetscInt        coneSize, c;

    1013:     PetscCall(DMPlexAddClosureTree_Up_Private(dm, ht, parent));
    1014:     PetscCall(DMPlexAddClosure_Private(dm, ht, parent));
    1015:     PetscCall(DMPlexGetCone(dm, parent, &cone));
    1016:     PetscCall(DMPlexGetConeSize(dm, parent, &coneSize));
    1017:     for (c = 0; c < coneSize; c++) {
    1018:       const PetscInt cp = cone[c];

    1020:       PetscCall(DMPlexAddClosureTree_Up_Private(dm, ht, cp));
    1021:     }
    1022:   }
    1023:   PetscFunctionReturn(PETSC_SUCCESS);
    1024: }

    1026: static PetscErrorCode DMPlexAddClosureTree_Down_Private(DM dm, PetscHSetI ht, PetscInt point)
    1027: {
    1028:   PetscInt        i, numChildren;
    1029:   const PetscInt *children;

    1031:   PetscFunctionBeginHot;
    1032:   PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, &children));
    1033:   for (i = 0; i < numChildren; i++) PetscCall(PetscHSetIAdd(ht, children[i]));
    1034:   PetscFunctionReturn(PETSC_SUCCESS);
    1035: }

    1037: static PetscErrorCode DMPlexAddClosureTree_Private(DM dm, PetscHSetI ht, PetscInt point)
    1038: {
    1039:   const PetscInt *cone;
    1040:   PetscInt        coneSize, c;

    1042:   PetscFunctionBeginHot;
    1043:   PetscCall(PetscHSetIAdd(ht, point));
    1044:   PetscCall(DMPlexAddClosureTree_Up_Private(dm, ht, point));
    1045:   PetscCall(DMPlexAddClosureTree_Down_Private(dm, ht, point));
    1046:   PetscCall(DMPlexGetCone(dm, point, &cone));
    1047:   PetscCall(DMPlexGetConeSize(dm, point, &coneSize));
    1048:   for (c = 0; c < coneSize; c++) PetscCall(DMPlexAddClosureTree_Private(dm, ht, cone[c]));
    1049:   PetscFunctionReturn(PETSC_SUCCESS);
    1050: }

    1052: PetscErrorCode DMPlexClosurePoints_Private(DM dm, PetscInt numPoints, const PetscInt points[], IS *closureIS)
    1053: {
    1054:   DM_Plex        *mesh    = (DM_Plex *)dm->data;
    1055:   const PetscBool hasTree = (mesh->parentSection || mesh->childSection) ? PETSC_TRUE : PETSC_FALSE;
    1056:   PetscInt        nelems, *elems, off = 0, p;
    1057:   PetscHSetI      ht = NULL;

    1059:   PetscFunctionBegin;
    1060:   PetscCall(PetscHSetICreate(&ht));
    1061:   PetscCall(PetscHSetIResize(ht, numPoints * 16));
    1062:   if (!hasTree) {
    1063:     for (p = 0; p < numPoints; ++p) PetscCall(DMPlexAddClosure_Private(dm, ht, points[p]));
    1064:   } else {
    1065: #if 1
    1066:     for (p = 0; p < numPoints; ++p) PetscCall(DMPlexAddClosureTree_Private(dm, ht, points[p]));
    1067: #else
    1068:     PetscInt *closure = NULL, closureSize, c;
    1069:     for (p = 0; p < numPoints; ++p) {
    1070:       PetscCall(DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure));
    1071:       for (c = 0; c < closureSize * 2; c += 2) {
    1072:         PetscCall(PetscHSetIAdd(ht, closure[c]));
    1073:         if (hasTree) PetscCall(DMPlexAddClosure_Tree(dm, ht, closure[c], PETSC_TRUE, PETSC_TRUE));
    1074:       }
    1075:     }
    1076:     if (closure) PetscCall(DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure));
    1077: #endif
    1078:   }
    1079:   PetscCall(PetscHSetIGetSize(ht, &nelems));
    1080:   PetscCall(PetscMalloc1(nelems, &elems));
    1081:   PetscCall(PetscHSetIGetElems(ht, &off, elems));
    1082:   PetscCall(PetscHSetIDestroy(&ht));
    1083:   PetscCall(PetscSortInt(nelems, elems));
    1084:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nelems, elems, PETSC_OWN_POINTER, closureIS));
    1085:   PetscFunctionReturn(PETSC_SUCCESS);
    1086: }

    1088: /*@
    1089:   DMPlexPartitionLabelClosure - Add the closure of all points to the partition label

    1091:   Input Parameters:
    1092: + dm    - The `DM`
    1093: - label - `DMLabel` assigning ranks to remote roots

    1095:   Level: developer

    1097: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
    1098: @*/
    1099: PetscErrorCode DMPlexPartitionLabelClosure(DM dm, DMLabel label)
    1100: {
    1101:   IS              rankIS, pointIS, closureIS;
    1102:   const PetscInt *ranks, *points;
    1103:   PetscInt        numRanks, numPoints, r;

    1105:   PetscFunctionBegin;
    1106:   PetscCall(DMLabelGetValueIS(label, &rankIS));
    1107:   PetscCall(ISGetLocalSize(rankIS, &numRanks));
    1108:   PetscCall(ISGetIndices(rankIS, &ranks));
    1109:   for (r = 0; r < numRanks; ++r) {
    1110:     const PetscInt rank = ranks[r];
    1111:     PetscCall(DMLabelGetStratumIS(label, rank, &pointIS));
    1112:     PetscCall(ISGetLocalSize(pointIS, &numPoints));
    1113:     PetscCall(ISGetIndices(pointIS, &points));
    1114:     PetscCall(DMPlexClosurePoints_Private(dm, numPoints, points, &closureIS));
    1115:     PetscCall(ISRestoreIndices(pointIS, &points));
    1116:     PetscCall(ISDestroy(&pointIS));
    1117:     PetscCall(DMLabelSetStratumIS(label, rank, closureIS));
    1118:     PetscCall(ISDestroy(&closureIS));
    1119:   }
    1120:   PetscCall(ISRestoreIndices(rankIS, &ranks));
    1121:   PetscCall(ISDestroy(&rankIS));
    1122:   PetscFunctionReturn(PETSC_SUCCESS);
    1123: }

    1125: /*@
    1126:   DMPlexPartitionLabelAdjacency - Add one level of adjacent points to the partition label

    1128:   Input Parameters:
    1129: + dm    - The `DM`
    1130: - label - `DMLabel` assigning ranks to remote roots

    1132:   Level: developer

    1134: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
    1135: @*/
    1136: PetscErrorCode DMPlexPartitionLabelAdjacency(DM dm, DMLabel label)
    1137: {
    1138:   IS              rankIS, pointIS;
    1139:   const PetscInt *ranks, *points;
    1140:   PetscInt        numRanks, numPoints, r, p, a, adjSize;
    1141:   PetscInt       *adj = NULL;

    1143:   PetscFunctionBegin;
    1144:   PetscCall(DMLabelGetValueIS(label, &rankIS));
    1145:   PetscCall(ISGetLocalSize(rankIS, &numRanks));
    1146:   PetscCall(ISGetIndices(rankIS, &ranks));
    1147:   for (r = 0; r < numRanks; ++r) {
    1148:     const PetscInt rank = ranks[r];

    1150:     PetscCall(DMLabelGetStratumIS(label, rank, &pointIS));
    1151:     PetscCall(ISGetLocalSize(pointIS, &numPoints));
    1152:     PetscCall(ISGetIndices(pointIS, &points));
    1153:     for (p = 0; p < numPoints; ++p) {
    1154:       adjSize = PETSC_DETERMINE;
    1155:       PetscCall(DMPlexGetAdjacency(dm, points[p], &adjSize, &adj));
    1156:       for (a = 0; a < adjSize; ++a) PetscCall(DMLabelSetValue(label, adj[a], rank));
    1157:     }
    1158:     PetscCall(ISRestoreIndices(pointIS, &points));
    1159:     PetscCall(ISDestroy(&pointIS));
    1160:   }
    1161:   PetscCall(ISRestoreIndices(rankIS, &ranks));
    1162:   PetscCall(ISDestroy(&rankIS));
    1163:   PetscCall(PetscFree(adj));
    1164:   PetscFunctionReturn(PETSC_SUCCESS);
    1165: }

    1167: /*@
    1168:   DMPlexPartitionLabelPropagate - Propagate points in a partition label over the point `PetscSF`

    1170:   Input Parameters:
    1171: + dm    - The `DM`
    1172: - label - `DMLabel` assigning ranks to remote roots

    1174:   Level: developer

    1176:   Note:
    1177:   This is required when generating multi-level overlaps to capture
    1178:   overlap points from non-neighbouring partitions.

    1180: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
    1181: @*/
    1182: PetscErrorCode DMPlexPartitionLabelPropagate(DM dm, DMLabel label)
    1183: {
    1184:   MPI_Comm        comm;
    1185:   PetscMPIInt     rank;
    1186:   PetscSF         sfPoint;
    1187:   DMLabel         lblRoots, lblLeaves;
    1188:   IS              rankIS, pointIS;
    1189:   const PetscInt *ranks;
    1190:   PetscInt        numRanks, r;

    1192:   PetscFunctionBegin;
    1193:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    1194:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1195:   PetscCall(DMGetPointSF(dm, &sfPoint));
    1196:   /* Pull point contributions from remote leaves into local roots */
    1197:   PetscCall(DMLabelGather(label, sfPoint, &lblLeaves));
    1198:   PetscCall(DMLabelGetValueIS(lblLeaves, &rankIS));
    1199:   PetscCall(ISGetLocalSize(rankIS, &numRanks));
    1200:   PetscCall(ISGetIndices(rankIS, &ranks));
    1201:   for (r = 0; r < numRanks; ++r) {
    1202:     const PetscInt remoteRank = ranks[r];
    1203:     if (remoteRank == rank) continue;
    1204:     PetscCall(DMLabelGetStratumIS(lblLeaves, remoteRank, &pointIS));
    1205:     PetscCall(DMLabelInsertIS(label, pointIS, remoteRank));
    1206:     PetscCall(ISDestroy(&pointIS));
    1207:   }
    1208:   PetscCall(ISRestoreIndices(rankIS, &ranks));
    1209:   PetscCall(ISDestroy(&rankIS));
    1210:   PetscCall(DMLabelDestroy(&lblLeaves));
    1211:   /* Push point contributions from roots into remote leaves */
    1212:   PetscCall(DMLabelDistribute(label, sfPoint, &lblRoots));
    1213:   PetscCall(DMLabelGetValueIS(lblRoots, &rankIS));
    1214:   PetscCall(ISGetLocalSize(rankIS, &numRanks));
    1215:   PetscCall(ISGetIndices(rankIS, &ranks));
    1216:   for (r = 0; r < numRanks; ++r) {
    1217:     const PetscInt remoteRank = ranks[r];
    1218:     if (remoteRank == rank) continue;
    1219:     PetscCall(DMLabelGetStratumIS(lblRoots, remoteRank, &pointIS));
    1220:     PetscCall(DMLabelInsertIS(label, pointIS, remoteRank));
    1221:     PetscCall(ISDestroy(&pointIS));
    1222:   }
    1223:   PetscCall(ISRestoreIndices(rankIS, &ranks));
    1224:   PetscCall(ISDestroy(&rankIS));
    1225:   PetscCall(DMLabelDestroy(&lblRoots));
    1226:   PetscFunctionReturn(PETSC_SUCCESS);
    1227: }

    1229: /*@
    1230:   DMPlexPartitionLabelInvert - Create a partition label of remote roots from a local root label

    1232:   Input Parameters:
    1233: + dm        - The `DM`
    1234: . rootLabel - `DMLabel` assigning ranks to local roots
    1235: - processSF - A star forest mapping into the local index on each remote rank

    1237:   Output Parameter:
    1238: . leafLabel - `DMLabel` assigning ranks to remote roots

    1240:   Level: developer

    1242:   Note:
    1243:   The rootLabel defines a send pattern by mapping local points to remote target ranks. The
    1244:   resulting leafLabel is a receiver mapping of remote roots to their parent rank.

    1246: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
    1247: @*/
    1248: PetscErrorCode DMPlexPartitionLabelInvert(DM dm, DMLabel rootLabel, PetscSF processSF, DMLabel leafLabel)
    1249: {
    1250:   MPI_Comm           comm;
    1251:   PetscMPIInt        rank, size, r;
    1252:   PetscInt           p, n, numNeighbors, numPoints, dof, off, rootSize, l, nleaves, leafSize;
    1253:   PetscSF            sfPoint;
    1254:   PetscSection       rootSection;
    1255:   PetscSFNode       *rootPoints, *leafPoints;
    1256:   const PetscSFNode *remote;
    1257:   const PetscInt    *local, *neighbors;
    1258:   IS                 valueIS;
    1259:   PetscBool          mpiOverflow = PETSC_FALSE;

    1261:   PetscFunctionBegin;
    1262:   PetscCall(PetscLogEventBegin(DMPLEX_PartLabelInvert, dm, 0, 0, 0));
    1263:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    1264:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1265:   PetscCallMPI(MPI_Comm_size(comm, &size));
    1266:   PetscCall(DMGetPointSF(dm, &sfPoint));

    1268:   /* Convert to (point, rank) and use actual owners */
    1269:   PetscCall(PetscSectionCreate(comm, &rootSection));
    1270:   PetscCall(PetscSectionSetChart(rootSection, 0, size));
    1271:   PetscCall(DMLabelGetValueIS(rootLabel, &valueIS));
    1272:   PetscCall(ISGetLocalSize(valueIS, &numNeighbors));
    1273:   PetscCall(ISGetIndices(valueIS, &neighbors));
    1274:   for (n = 0; n < numNeighbors; ++n) {
    1275:     PetscCall(DMLabelGetStratumSize(rootLabel, neighbors[n], &numPoints));
    1276:     PetscCall(PetscSectionAddDof(rootSection, neighbors[n], numPoints));
    1277:   }
    1278:   PetscCall(PetscSectionSetUp(rootSection));
    1279:   PetscCall(PetscSectionGetStorageSize(rootSection, &rootSize));
    1280:   PetscCall(PetscMalloc1(rootSize, &rootPoints));
    1281:   PetscCall(PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote));
    1282:   for (n = 0; n < numNeighbors; ++n) {
    1283:     IS              pointIS;
    1284:     const PetscInt *points;

    1286:     PetscCall(PetscSectionGetOffset(rootSection, neighbors[n], &off));
    1287:     PetscCall(DMLabelGetStratumIS(rootLabel, neighbors[n], &pointIS));
    1288:     PetscCall(ISGetLocalSize(pointIS, &numPoints));
    1289:     PetscCall(ISGetIndices(pointIS, &points));
    1290:     for (p = 0; p < numPoints; ++p) {
    1291:       if (local) PetscCall(PetscFindInt(points[p], nleaves, local, &l));
    1292:       else l = -1;
    1293:       if (l >= 0) {
    1294:         rootPoints[off + p] = remote[l];
    1295:       } else {
    1296:         rootPoints[off + p].index = points[p];
    1297:         rootPoints[off + p].rank  = rank;
    1298:       }
    1299:     }
    1300:     PetscCall(ISRestoreIndices(pointIS, &points));
    1301:     PetscCall(ISDestroy(&pointIS));
    1302:   }

    1304:   /* Try to communicate overlap using All-to-All */
    1305:   if (!processSF) {
    1306:     PetscCount   counter     = 0;
    1307:     PetscBool    locOverflow = PETSC_FALSE;
    1308:     PetscMPIInt *scounts, *sdispls, *rcounts, *rdispls;

    1310:     PetscCall(PetscCalloc4(size, &scounts, size, &sdispls, size, &rcounts, size, &rdispls));
    1311:     for (n = 0; n < numNeighbors; ++n) {
    1312:       PetscCall(PetscSectionGetDof(rootSection, neighbors[n], &dof));
    1313:       PetscCall(PetscSectionGetOffset(rootSection, neighbors[n], &off));
    1314: #if defined(PETSC_USE_64BIT_INDICES)
    1315:       if (dof > PETSC_MPI_INT_MAX) {
    1316:         locOverflow = PETSC_TRUE;
    1317:         break;
    1318:       }
    1319:       if (off > PETSC_MPI_INT_MAX) {
    1320:         locOverflow = PETSC_TRUE;
    1321:         break;
    1322:       }
    1323: #endif
    1324:       scounts[neighbors[n]] = (PetscMPIInt)dof;
    1325:       sdispls[neighbors[n]] = (PetscMPIInt)off;
    1326:     }
    1327:     PetscCallMPI(MPI_Alltoall(scounts, 1, MPI_INT, rcounts, 1, MPI_INT, comm));
    1328:     for (r = 0; r < size; ++r) {
    1329:       rdispls[r] = (PetscMPIInt)counter;
    1330:       counter += rcounts[r];
    1331:     }
    1332:     if (counter > PETSC_MPI_INT_MAX) locOverflow = PETSC_TRUE;
    1333:     PetscCallMPI(MPIU_Allreduce(&locOverflow, &mpiOverflow, 1, MPIU_BOOL, MPI_LOR, comm));
    1334:     if (!mpiOverflow) {
    1335:       PetscCall(PetscInfo(dm, "Using MPI_Alltoallv() for mesh distribution\n"));
    1336:       leafSize = (PetscInt)counter;
    1337:       PetscCall(PetscMalloc1(leafSize, &leafPoints));
    1338:       PetscCallMPI(MPI_Alltoallv(rootPoints, scounts, sdispls, MPIU_SF_NODE, leafPoints, rcounts, rdispls, MPIU_SF_NODE, comm));
    1339:     }
    1340:     PetscCall(PetscFree4(scounts, sdispls, rcounts, rdispls));
    1341:   }

    1343:   /* Communicate overlap using process star forest */
    1344:   if (processSF || mpiOverflow) {
    1345:     PetscSF      procSF;
    1346:     PetscSection leafSection;

    1348:     if (processSF) {
    1349:       PetscCall(PetscInfo(dm, "Using processSF for mesh distribution\n"));
    1350:       PetscCall(PetscObjectReference((PetscObject)processSF));
    1351:       procSF = processSF;
    1352:     } else {
    1353:       PetscCall(PetscInfo(dm, "Using processSF for mesh distribution (MPI overflow)\n"));
    1354:       PetscCall(PetscSFCreate(comm, &procSF));
    1355:       PetscCall(PetscSFSetGraphWithPattern(procSF, NULL, PETSCSF_PATTERN_ALLTOALL));
    1356:     }

    1358:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &leafSection));
    1359:     PetscCall(DMPlexDistributeData(dm, procSF, rootSection, MPIU_SF_NODE, rootPoints, leafSection, (void **)&leafPoints));
    1360:     PetscCall(PetscSectionGetStorageSize(leafSection, &leafSize));
    1361:     PetscCall(PetscSectionDestroy(&leafSection));
    1362:     PetscCall(PetscSFDestroy(&procSF));
    1363:   }

    1365:   for (p = 0; p < leafSize; p++) PetscCall(DMLabelSetValue(leafLabel, leafPoints[p].index, leafPoints[p].rank));

    1367:   PetscCall(ISRestoreIndices(valueIS, &neighbors));
    1368:   PetscCall(ISDestroy(&valueIS));
    1369:   PetscCall(PetscSectionDestroy(&rootSection));
    1370:   PetscCall(PetscFree(rootPoints));
    1371:   PetscCall(PetscFree(leafPoints));
    1372:   PetscCall(PetscLogEventEnd(DMPLEX_PartLabelInvert, dm, 0, 0, 0));
    1373:   PetscFunctionReturn(PETSC_SUCCESS);
    1374: }

    1376: /*@
    1377:   DMPlexPartitionLabelCreateSF - Create a star forest from a label that assigns ranks to points

    1379:   Input Parameters:
    1380: + dm        - The `DM`
    1381: . label     - `DMLabel` assigning ranks to remote roots
    1382: - sortRanks - Whether or not to sort the `PetscSF` leaves by rank

    1384:   Output Parameter:
    1385: . sf - The star forest communication context encapsulating the defined mapping

    1387:   Level: developer

    1389:   Note:
    1390:   The incoming label is a receiver mapping of remote points to their parent rank.

    1392: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabel`, `PetscSF`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
    1393: @*/
    1394: PetscErrorCode DMPlexPartitionLabelCreateSF(DM dm, DMLabel label, PetscBool sortRanks, PetscSF *sf)
    1395: {
    1396:   PetscMPIInt     rank;
    1397:   PetscInt        n, numRemote, p, numPoints, pStart, pEnd, idx = 0, nNeighbors;
    1398:   PetscSFNode    *remotePoints;
    1399:   IS              remoteRootIS, neighborsIS;
    1400:   const PetscInt *remoteRoots, *neighbors;
    1401:   PetscMPIInt     myRank;

    1403:   PetscFunctionBegin;
    1404:   PetscCall(PetscLogEventBegin(DMPLEX_PartLabelCreateSF, dm, 0, 0, 0));
    1405:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    1406:   PetscCall(DMLabelGetValueIS(label, &neighborsIS));

    1408:   if (sortRanks) {
    1409:     IS is;

    1411:     PetscCall(ISDuplicate(neighborsIS, &is));
    1412:     PetscCall(ISSort(is));
    1413:     PetscCall(ISDestroy(&neighborsIS));
    1414:     neighborsIS = is;
    1415:   }
    1416:   myRank = sortRanks ? -1 : rank;

    1418:   PetscCall(ISGetLocalSize(neighborsIS, &nNeighbors));
    1419:   PetscCall(ISGetIndices(neighborsIS, &neighbors));
    1420:   for (numRemote = 0, n = 0; n < nNeighbors; ++n) {
    1421:     PetscCall(DMLabelGetStratumSize(label, neighbors[n], &numPoints));
    1422:     numRemote += numPoints;
    1423:   }
    1424:   PetscCall(PetscMalloc1(numRemote, &remotePoints));

    1426:   /* Put owned points first if not sorting the ranks */
    1427:   if (!sortRanks) {
    1428:     PetscCall(DMLabelGetStratumSize(label, rank, &numPoints));
    1429:     if (numPoints > 0) {
    1430:       PetscCall(DMLabelGetStratumIS(label, rank, &remoteRootIS));
    1431:       PetscCall(ISGetIndices(remoteRootIS, &remoteRoots));
    1432:       for (p = 0; p < numPoints; p++) {
    1433:         remotePoints[idx].index = remoteRoots[p];
    1434:         remotePoints[idx].rank  = rank;
    1435:         idx++;
    1436:       }
    1437:       PetscCall(ISRestoreIndices(remoteRootIS, &remoteRoots));
    1438:       PetscCall(ISDestroy(&remoteRootIS));
    1439:     }
    1440:   }

    1442:   /* Now add remote points */
    1443:   for (n = 0; n < nNeighbors; ++n) {
    1444:     PetscMPIInt nn;

    1446:     PetscCall(PetscMPIIntCast(neighbors[n], &nn));
    1447:     PetscCall(DMLabelGetStratumSize(label, nn, &numPoints));
    1448:     if (nn == myRank || numPoints <= 0) continue;
    1449:     PetscCall(DMLabelGetStratumIS(label, nn, &remoteRootIS));
    1450:     PetscCall(ISGetIndices(remoteRootIS, &remoteRoots));
    1451:     for (p = 0; p < numPoints; p++) {
    1452:       remotePoints[idx].index = remoteRoots[p];
    1453:       remotePoints[idx].rank  = nn;
    1454:       idx++;
    1455:     }
    1456:     PetscCall(ISRestoreIndices(remoteRootIS, &remoteRoots));
    1457:     PetscCall(ISDestroy(&remoteRootIS));
    1458:   }

    1460:   PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), sf));
    1461:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    1462:   PetscCall(PetscSFSetGraph(*sf, pEnd - pStart, numRemote, NULL, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER));
    1463:   PetscCall(PetscSFSetUp(*sf));
    1464:   PetscCall(ISDestroy(&neighborsIS));
    1465:   PetscCall(PetscLogEventEnd(DMPLEX_PartLabelCreateSF, dm, 0, 0, 0));
    1466:   PetscFunctionReturn(PETSC_SUCCESS);
    1467: }

    1469: #if defined(PETSC_HAVE_PARMETIS)
    1470:   #include <parmetis.h>
    1471: #endif

    1473: /* The two functions below are used by DMPlexRebalanceSharedPoints which errors
    1474:  * when PETSc is built without ParMETIS. To avoid -Wunused-function, we take
    1475:  * them out in that case. */
    1476: #if defined(PETSC_HAVE_PARMETIS)
    1477: /*
    1478:   DMPlexRewriteSF - Rewrites the ownership of the `PetscSF` of a `DM` (in place).

    1480:   Input parameters:b
    1481: + dm                - The `DMPLEX` object.
    1482: . n                 - The number of points.
    1483: . pointsToRewrite   - The points in the `PetscSF` whose ownership will change.
    1484: . targetOwners      - New owner for each element in pointsToRewrite.
    1485: - degrees           - Degrees of the points in the `PetscSF` as obtained by `PetscSFComputeDegreeBegin()`/`PetscSFComputeDegreeEnd()`.

    1487:   Level: developer

    1489: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabel`, `PetscSF`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
    1490: */
    1491: static PetscErrorCode DMPlexRewriteSF(DM dm, PetscInt n, PetscInt *pointsToRewrite, PetscInt *targetOwners, const PetscInt *degrees)
    1492: {
    1493:   PetscInt           pStart, pEnd, i, j, counter, leafCounter, sumDegrees, nroots, nleafs;
    1494:   PetscInt          *cumSumDegrees, *newOwners, *newNumbers, *rankOnLeafs, *locationsOfLeafs, *remoteLocalPointOfLeafs, *points, *leafsNew;
    1495:   PetscSFNode       *leafLocationsNew;
    1496:   const PetscSFNode *iremote;
    1497:   const PetscInt    *ilocal;
    1498:   PetscBool         *isLeaf;
    1499:   PetscSF            sf;
    1500:   MPI_Comm           comm;
    1501:   PetscMPIInt        rank, size;

    1503:   PetscFunctionBegin;
    1504:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    1505:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1506:   PetscCallMPI(MPI_Comm_size(comm, &size));
    1507:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));

    1509:   PetscCall(DMGetPointSF(dm, &sf));
    1510:   PetscCall(PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote));
    1511:   PetscCall(PetscMalloc1(pEnd - pStart, &isLeaf));
    1512:   for (i = 0; i < pEnd - pStart; i++) isLeaf[i] = PETSC_FALSE;
    1513:   for (i = 0; i < nleafs; i++) isLeaf[ilocal[i] - pStart] = PETSC_TRUE;

    1515:   PetscCall(PetscMalloc1(pEnd - pStart + 1, &cumSumDegrees));
    1516:   cumSumDegrees[0] = 0;
    1517:   for (i = 1; i <= pEnd - pStart; i++) cumSumDegrees[i] = cumSumDegrees[i - 1] + degrees[i - 1];
    1518:   sumDegrees = cumSumDegrees[pEnd - pStart];
    1519:   /* get the location of my leafs (we have sumDegrees many leafs pointing at our roots) */

    1521:   PetscCall(PetscMalloc1(sumDegrees, &locationsOfLeafs));
    1522:   PetscCall(PetscMalloc1(pEnd - pStart, &rankOnLeafs));
    1523:   for (i = 0; i < pEnd - pStart; i++) rankOnLeafs[i] = rank;
    1524:   PetscCall(PetscSFGatherBegin(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs));
    1525:   PetscCall(PetscSFGatherEnd(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs));
    1526:   PetscCall(PetscFree(rankOnLeafs));

    1528:   /* get the remote local points of my leaves */
    1529:   PetscCall(PetscMalloc1(sumDegrees, &remoteLocalPointOfLeafs));
    1530:   PetscCall(PetscMalloc1(pEnd - pStart, &points));
    1531:   for (i = 0; i < pEnd - pStart; i++) points[i] = pStart + i;
    1532:   PetscCall(PetscSFGatherBegin(sf, MPIU_INT, points, remoteLocalPointOfLeafs));
    1533:   PetscCall(PetscSFGatherEnd(sf, MPIU_INT, points, remoteLocalPointOfLeafs));
    1534:   PetscCall(PetscFree(points));
    1535:   /* Figure out the new owners of the vertices that are up for grabs and their numbers on the new owners */
    1536:   PetscCall(PetscMalloc1(pEnd - pStart, &newOwners));
    1537:   PetscCall(PetscMalloc1(pEnd - pStart, &newNumbers));
    1538:   for (i = 0; i < pEnd - pStart; i++) {
    1539:     newOwners[i]  = -1;
    1540:     newNumbers[i] = -1;
    1541:   }
    1542:   {
    1543:     PetscInt oldNumber, newNumber, oldOwner, newOwner;
    1544:     for (i = 0; i < n; i++) {
    1545:       oldNumber = pointsToRewrite[i];
    1546:       newNumber = -1;
    1547:       oldOwner  = rank;
    1548:       newOwner  = targetOwners[i];
    1549:       if (oldOwner == newOwner) {
    1550:         newNumber = oldNumber;
    1551:       } else {
    1552:         for (j = 0; j < degrees[oldNumber]; j++) {
    1553:           if (locationsOfLeafs[cumSumDegrees[oldNumber] + j] == newOwner) {
    1554:             newNumber = remoteLocalPointOfLeafs[cumSumDegrees[oldNumber] + j];
    1555:             break;
    1556:           }
    1557:         }
    1558:       }
    1559:       PetscCheck(newNumber != -1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Couldn't find the new owner of vertex.");

    1561:       newOwners[oldNumber]  = newOwner;
    1562:       newNumbers[oldNumber] = newNumber;
    1563:     }
    1564:   }
    1565:   PetscCall(PetscFree(cumSumDegrees));
    1566:   PetscCall(PetscFree(locationsOfLeafs));
    1567:   PetscCall(PetscFree(remoteLocalPointOfLeafs));

    1569:   PetscCall(PetscSFBcastBegin(sf, MPIU_INT, newOwners, newOwners, MPI_REPLACE));
    1570:   PetscCall(PetscSFBcastEnd(sf, MPIU_INT, newOwners, newOwners, MPI_REPLACE));
    1571:   PetscCall(PetscSFBcastBegin(sf, MPIU_INT, newNumbers, newNumbers, MPI_REPLACE));
    1572:   PetscCall(PetscSFBcastEnd(sf, MPIU_INT, newNumbers, newNumbers, MPI_REPLACE));

    1574:   /* Now count how many leafs we have on each processor. */
    1575:   leafCounter = 0;
    1576:   for (i = 0; i < pEnd - pStart; i++) {
    1577:     if (newOwners[i] >= 0) {
    1578:       if (newOwners[i] != rank) leafCounter++;
    1579:     } else {
    1580:       if (isLeaf[i]) leafCounter++;
    1581:     }
    1582:   }

    1584:   /* Now set up the new sf by creating the leaf arrays */
    1585:   PetscCall(PetscMalloc1(leafCounter, &leafsNew));
    1586:   PetscCall(PetscMalloc1(leafCounter, &leafLocationsNew));

    1588:   leafCounter = 0;
    1589:   counter     = 0;
    1590:   for (i = 0; i < pEnd - pStart; i++) {
    1591:     if (newOwners[i] >= 0) {
    1592:       if (newOwners[i] != rank) {
    1593:         leafsNew[leafCounter]               = i;
    1594:         leafLocationsNew[leafCounter].rank  = newOwners[i];
    1595:         leafLocationsNew[leafCounter].index = newNumbers[i];
    1596:         leafCounter++;
    1597:       }
    1598:     } else {
    1599:       if (isLeaf[i]) {
    1600:         leafsNew[leafCounter]               = i;
    1601:         leafLocationsNew[leafCounter].rank  = iremote[counter].rank;
    1602:         leafLocationsNew[leafCounter].index = iremote[counter].index;
    1603:         leafCounter++;
    1604:       }
    1605:     }
    1606:     if (isLeaf[i]) counter++;
    1607:   }

    1609:   PetscCall(PetscSFSetGraph(sf, nroots, leafCounter, leafsNew, PETSC_OWN_POINTER, leafLocationsNew, PETSC_OWN_POINTER));
    1610:   PetscCall(PetscFree(newOwners));
    1611:   PetscCall(PetscFree(newNumbers));
    1612:   PetscCall(PetscFree(isLeaf));
    1613:   PetscFunctionReturn(PETSC_SUCCESS);
    1614: }

    1616: static PetscErrorCode DMPlexViewDistribution(MPI_Comm comm, PetscInt n, PetscInt skip, PetscInt *vtxwgt, PetscInt *part, PetscViewer viewer)
    1617: {
    1618:   PetscInt   *distribution, min, max, sum;
    1619:   PetscMPIInt rank, size;

    1621:   PetscFunctionBegin;
    1622:   PetscCallMPI(MPI_Comm_size(comm, &size));
    1623:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1624:   PetscCall(PetscCalloc1(size, &distribution));
    1625:   for (PetscInt i = 0; i < n; i++) {
    1626:     if (part) distribution[part[i]] += vtxwgt[skip * i];
    1627:     else distribution[rank] += vtxwgt[skip * i];
    1628:   }
    1629:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, distribution, size, MPIU_INT, MPI_SUM, comm));
    1630:   min = distribution[0];
    1631:   max = distribution[0];
    1632:   sum = distribution[0];
    1633:   for (PetscInt i = 1; i < size; i++) {
    1634:     if (distribution[i] < min) min = distribution[i];
    1635:     if (distribution[i] > max) max = distribution[i];
    1636:     sum += distribution[i];
    1637:   }
    1638:   PetscCall(PetscViewerASCIIPrintf(viewer, "Min: %" PetscInt_FMT ", Avg: %" PetscInt_FMT ", Max: %" PetscInt_FMT ", Balance: %f\n", min, sum / size, max, (max * 1. * size) / sum));
    1639:   PetscCall(PetscFree(distribution));
    1640:   PetscFunctionReturn(PETSC_SUCCESS);
    1641: }

    1643: #endif

    1645: /*@
    1646:   DMPlexRebalanceSharedPoints - Redistribute points in the plex that are shared in order to achieve better balancing. This routine updates the `PointSF` of the `DM` inplace.

    1648:   Input Parameters:
    1649: + dm              - The `DMPLEX` object.
    1650: . entityDepth     - depth of the entity to balance (0 -> balance vertices).
    1651: . useInitialGuess - whether to use the current distribution as initial guess (only used by ParMETIS).
    1652: - parallel        - whether to use ParMETIS and do the partition in parallel or whether to gather the graph onto a single process and use METIS.

    1654:   Output Parameter:
    1655: . success - whether the graph partitioning was successful or not, optional. Unsuccessful simply means no change to the partitioning

    1657:   Options Database Keys:
    1658: + -dm_plex_rebalance_shared_points_parmetis             - Use ParMetis instead of Metis for the partitioner
    1659: . -dm_plex_rebalance_shared_points_use_initial_guess    - Use current partition to bootstrap ParMetis partition
    1660: . -dm_plex_rebalance_shared_points_use_mat_partitioning - Use the MatPartitioning object to perform the partition, the prefix for those operations is -dm_plex_rebalance_shared_points_
    1661: - -dm_plex_rebalance_shared_points_monitor              - Monitor the shared points rebalance process

    1663:   Level: intermediate

    1665: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
    1666: @*/
    1667: PetscErrorCode DMPlexRebalanceSharedPoints(DM dm, PetscInt entityDepth, PetscBool useInitialGuess, PetscBool parallel, PetscBool *success)
    1668: {
    1669: #if defined(PETSC_HAVE_PARMETIS)
    1670:   PetscSF            sf;
    1671:   PetscInt           ierr, i, j, idx, jdx;
    1672:   PetscInt           eBegin, eEnd, nroots, nleafs, pStart, pEnd;
    1673:   const PetscInt    *degrees, *ilocal;
    1674:   const PetscSFNode *iremote;
    1675:   PetscBool         *toBalance, *isLeaf, *isExclusivelyOwned, *isNonExclusivelyOwned;
    1676:   PetscInt           numExclusivelyOwned, numNonExclusivelyOwned;
    1677:   PetscMPIInt        rank, size;
    1678:   PetscInt          *globalNumbersOfLocalOwnedVertices, *leafGlobalNumbers;
    1679:   const PetscInt    *cumSumVertices;
    1680:   PetscInt           offset, counter;
    1681:   PetscInt          *vtxwgt;
    1682:   const PetscInt    *xadj, *adjncy;
    1683:   PetscInt          *part, *options;
    1684:   PetscInt           nparts, wgtflag, numflag, ncon, edgecut;
    1685:   real_t            *ubvec;
    1686:   PetscInt          *firstVertices, *renumbering;
    1687:   PetscInt           failed, failedGlobal;
    1688:   MPI_Comm           comm;
    1689:   Mat                A;
    1690:   PetscViewer        viewer;
    1691:   PetscViewerFormat  format;
    1692:   PetscLayout        layout;
    1693:   real_t            *tpwgts;
    1694:   PetscMPIInt       *counts, *mpiCumSumVertices;
    1695:   PetscInt          *pointsToRewrite;
    1696:   PetscInt           numRows;
    1697:   PetscBool          done, usematpartitioning = PETSC_FALSE;
    1698:   IS                 ispart = NULL;
    1699:   MatPartitioning    mp;
    1700:   const char        *prefix;

    1702:   PetscFunctionBegin;
    1703:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    1704:   PetscCallMPI(MPI_Comm_size(comm, &size));
    1705:   if (size == 1) {
    1706:     if (success) *success = PETSC_TRUE;
    1707:     PetscFunctionReturn(PETSC_SUCCESS);
    1708:   }
    1709:   if (success) *success = PETSC_FALSE;
    1710:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    1712:   parallel        = PETSC_FALSE;
    1713:   useInitialGuess = PETSC_FALSE;
    1714:   PetscObjectOptionsBegin((PetscObject)dm);
    1715:   PetscCall(PetscOptionsName("-dm_plex_rebalance_shared_points_parmetis", "Use ParMetis instead of Metis for the partitioner", "DMPlexRebalanceSharedPoints", &parallel));
    1716:   PetscCall(PetscOptionsBool("-dm_plex_rebalance_shared_points_use_initial_guess", "Use current partition to bootstrap ParMetis partition", "DMPlexRebalanceSharedPoints", useInitialGuess, &useInitialGuess, NULL));
    1717:   PetscCall(PetscOptionsBool("-dm_plex_rebalance_shared_points_use_mat_partitioning", "Use the MatPartitioning object to partition", "DMPlexRebalanceSharedPoints", usematpartitioning, &usematpartitioning, NULL));
    1718:   PetscCall(PetscOptionsViewer("-dm_plex_rebalance_shared_points_monitor", "Monitor the shared points rebalance process", "DMPlexRebalanceSharedPoints", &viewer, &format, NULL));
    1719:   PetscOptionsEnd();
    1720:   if (viewer) PetscCall(PetscViewerPushFormat(viewer, format));

    1722:   PetscCall(PetscLogEventBegin(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0));

    1724:   PetscCall(DMGetOptionsPrefix(dm, &prefix));
    1725:   PetscCall(PetscOptionsCreateViewer(comm, ((PetscObject)dm)->options, prefix, "-dm_rebalance_partition_view", &viewer, &format, NULL));
    1726:   if (viewer) PetscCall(PetscViewerPushFormat(viewer, format));

    1728:   /* Figure out all points in the plex that we are interested in balancing. */
    1729:   PetscCall(DMPlexGetDepthStratum(dm, entityDepth, &eBegin, &eEnd));
    1730:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    1731:   PetscCall(PetscMalloc1(pEnd - pStart, &toBalance));
    1732:   for (i = 0; i < pEnd - pStart; i++) toBalance[i] = (PetscBool)(i >= eBegin && i < eEnd);

    1734:   /* There are three types of points:
    1735:    * exclusivelyOwned: points that are owned by this process and only seen by this process
    1736:    * nonExclusivelyOwned: points that are owned by this process but seen by at least another process
    1737:    * leaf: a point that is seen by this process but owned by a different process
    1738:    */
    1739:   PetscCall(DMGetPointSF(dm, &sf));
    1740:   PetscCall(PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote));
    1741:   PetscCall(PetscMalloc1(pEnd - pStart, &isLeaf));
    1742:   PetscCall(PetscMalloc1(pEnd - pStart, &isNonExclusivelyOwned));
    1743:   PetscCall(PetscMalloc1(pEnd - pStart, &isExclusivelyOwned));
    1744:   for (i = 0; i < pEnd - pStart; i++) {
    1745:     isNonExclusivelyOwned[i] = PETSC_FALSE;
    1746:     isExclusivelyOwned[i]    = PETSC_FALSE;
    1747:     isLeaf[i]                = PETSC_FALSE;
    1748:   }

    1750:   /* mark all the leafs */
    1751:   for (i = 0; i < nleafs; i++) isLeaf[ilocal[i] - pStart] = PETSC_TRUE;

    1753:   /* for an owned point, we can figure out whether another processor sees it or
    1754:    * not by calculating its degree */
    1755:   PetscCall(PetscSFComputeDegreeBegin(sf, &degrees));
    1756:   PetscCall(PetscSFComputeDegreeEnd(sf, &degrees));
    1757:   numExclusivelyOwned    = 0;
    1758:   numNonExclusivelyOwned = 0;
    1759:   for (i = 0; i < pEnd - pStart; i++) {
    1760:     if (toBalance[i]) {
    1761:       if (degrees[i] > 0) {
    1762:         isNonExclusivelyOwned[i] = PETSC_TRUE;
    1763:         numNonExclusivelyOwned += 1;
    1764:       } else {
    1765:         if (!isLeaf[i]) {
    1766:           isExclusivelyOwned[i] = PETSC_TRUE;
    1767:           numExclusivelyOwned += 1;
    1768:         }
    1769:       }
    1770:     }
    1771:   }

    1773:   /* Build a graph with one vertex per core representing the
    1774:    * exclusively owned points and then one vertex per nonExclusively owned
    1775:    * point. */
    1776:   PetscCall(PetscLayoutCreate(comm, &layout));
    1777:   PetscCall(PetscLayoutSetLocalSize(layout, 1 + numNonExclusivelyOwned));
    1778:   PetscCall(PetscLayoutSetUp(layout));
    1779:   PetscCall(PetscLayoutGetRanges(layout, &cumSumVertices));
    1780:   PetscCall(PetscMalloc1(pEnd - pStart, &globalNumbersOfLocalOwnedVertices));
    1781:   for (i = 0; i < pEnd - pStart; i++) globalNumbersOfLocalOwnedVertices[i] = pStart - 1;
    1782:   offset  = cumSumVertices[rank];
    1783:   counter = 0;
    1784:   for (i = 0; i < pEnd - pStart; i++) {
    1785:     if (toBalance[i]) {
    1786:       if (degrees[i] > 0) {
    1787:         globalNumbersOfLocalOwnedVertices[i] = counter + 1 + offset;
    1788:         counter++;
    1789:       }
    1790:     }
    1791:   }

    1793:   /* send the global numbers of vertices I own to the leafs so that they know to connect to it */
    1794:   PetscCall(PetscMalloc1(pEnd - pStart, &leafGlobalNumbers));
    1795:   PetscCall(PetscSFBcastBegin(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers, MPI_REPLACE));
    1796:   PetscCall(PetscSFBcastEnd(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers, MPI_REPLACE));

    1798:   /* Build the graph for partitioning */
    1799:   numRows = 1 + numNonExclusivelyOwned;
    1800:   PetscCall(PetscLogEventBegin(DMPLEX_RebalBuildGraph, dm, 0, 0, 0));
    1801:   PetscCall(MatCreate(comm, &A));
    1802:   PetscCall(MatSetType(A, MATMPIADJ));
    1803:   PetscCall(MatSetSizes(A, numRows, numRows, cumSumVertices[size], cumSumVertices[size]));
    1804:   idx = cumSumVertices[rank];
    1805:   for (i = 0; i < pEnd - pStart; i++) {
    1806:     if (toBalance[i]) {
    1807:       if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
    1808:       else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
    1809:       else continue;
    1810:       PetscCall(MatSetValue(A, idx, jdx, 1, INSERT_VALUES));
    1811:       PetscCall(MatSetValue(A, jdx, idx, 1, INSERT_VALUES));
    1812:     }
    1813:   }
    1814:   PetscCall(PetscFree(globalNumbersOfLocalOwnedVertices));
    1815:   PetscCall(PetscFree(leafGlobalNumbers));
    1816:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    1817:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    1818:   PetscCall(PetscLogEventEnd(DMPLEX_RebalBuildGraph, dm, 0, 0, 0));

    1820:   nparts = size;
    1821:   ncon   = 1;
    1822:   PetscCall(PetscMalloc1(ncon * nparts, &tpwgts));
    1823:   for (i = 0; i < ncon * nparts; i++) tpwgts[i] = (real_t)(1. / (nparts));
    1824:   PetscCall(PetscMalloc1(ncon, &ubvec));
    1825:   for (i = 0; i < ncon; i++) ubvec[i] = (real_t)1.05;

    1827:   PetscCall(PetscMalloc1(ncon * (1 + numNonExclusivelyOwned), &vtxwgt));
    1828:   if (ncon == 2) {
    1829:     vtxwgt[0] = numExclusivelyOwned;
    1830:     vtxwgt[1] = 1;
    1831:     for (i = 0; i < numNonExclusivelyOwned; i++) {
    1832:       vtxwgt[ncon * (i + 1)]     = 1;
    1833:       vtxwgt[ncon * (i + 1) + 1] = 0;
    1834:     }
    1835:   } else {
    1836:     PetscInt base, ms;
    1837:     PetscCallMPI(MPIU_Allreduce(&numExclusivelyOwned, &base, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
    1838:     PetscCall(MatGetSize(A, &ms, NULL));
    1839:     ms -= size;
    1840:     base      = PetscMax(base, ms);
    1841:     vtxwgt[0] = base + numExclusivelyOwned;
    1842:     for (i = 0; i < numNonExclusivelyOwned; i++) vtxwgt[i + 1] = 1;
    1843:   }

    1845:   if (viewer) {
    1846:     PetscCall(PetscViewerASCIIPrintf(viewer, "Attempt rebalancing of shared points of depth %" PetscInt_FMT " on interface of mesh distribution.\n", entityDepth));
    1847:     PetscCall(PetscViewerASCIIPrintf(viewer, "Size of generated auxiliary graph: %" PetscInt_FMT "\n", cumSumVertices[size]));
    1848:   }
    1849:   /* TODO: Drop the parallel/sequential choice here and just use MatPartioner for much more flexibility */
    1850:   if (usematpartitioning) {
    1851:     const char *prefix;

    1853:     PetscCall(MatPartitioningCreate(PetscObjectComm((PetscObject)dm), &mp));
    1854:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mp, "dm_plex_rebalance_shared_points_"));
    1855:     PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
    1856:     PetscCall(PetscObjectPrependOptionsPrefix((PetscObject)mp, prefix));
    1857:     PetscCall(MatPartitioningSetAdjacency(mp, A));
    1858:     PetscCall(MatPartitioningSetNumberVertexWeights(mp, ncon));
    1859:     PetscCall(MatPartitioningSetVertexWeights(mp, vtxwgt));
    1860:     PetscCall(MatPartitioningSetFromOptions(mp));
    1861:     PetscCall(MatPartitioningApply(mp, &ispart));
    1862:     PetscCall(ISGetIndices(ispart, (const PetscInt **)&part));
    1863:   } else if (parallel) {
    1864:     if (viewer) PetscCall(PetscViewerASCIIPrintf(viewer, "Using ParMETIS to partition graph.\n"));
    1865:     PetscCall(PetscMalloc1(cumSumVertices[rank + 1] - cumSumVertices[rank], &part));
    1866:     PetscCall(MatGetRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, &xadj, &adjncy, &done));
    1867:     PetscCheck(done, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Could not get adjacency information");
    1868:     PetscCall(PetscMalloc1(4, &options));
    1869:     options[0] = 1;
    1870:     options[1] = 0; /* Verbosity */
    1871:     if (viewer) options[1] = 1;
    1872:     options[2] = 0;                    /* Seed */
    1873:     options[3] = PARMETIS_PSR_COUPLED; /* Seed */
    1874:     wgtflag    = 2;
    1875:     numflag    = 0;
    1876:     if (useInitialGuess) {
    1877:       PetscCall(PetscViewerASCIIPrintf(viewer, "THIS DOES NOT WORK! I don't know why. Using current distribution of points as initial guess.\n"));
    1878:       for (i = 0; i < numRows; i++) part[i] = rank;
    1879:       if (viewer) PetscCall(PetscViewerASCIIPrintf(viewer, "Using current distribution of points as initial guess.\n"));
    1880:       PetscStackPushExternal("ParMETIS_V3_RefineKway");
    1881:       PetscCall(PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0));
    1882:       ierr = ParMETIS_V3_RefineKway((PetscInt *)cumSumVertices, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt, NULL, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
    1883:       PetscCall(PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0));
    1884:       PetscStackPop;
    1885:       PetscCheck(ierr == METIS_OK, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ParMETIS_V3_RefineKway()");
    1886:     } else {
    1887:       PetscStackPushExternal("ParMETIS_V3_PartKway");
    1888:       PetscCall(PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0));
    1889:       ierr = ParMETIS_V3_PartKway((PetscInt *)cumSumVertices, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt, NULL, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
    1890:       PetscCall(PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0));
    1891:       PetscStackPop;
    1892:       PetscCheck(ierr == METIS_OK, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ParMETIS_V3_PartKway()");
    1893:     }
    1894:     PetscCall(MatRestoreRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, &xadj, &adjncy, &done));
    1895:     PetscCall(PetscFree(options));
    1896:   } else {
    1897:     if (viewer) PetscCall(PetscViewerASCIIPrintf(viewer, "Using METIS to partition graph.\n"));
    1898:     Mat       As;
    1899:     PetscInt *partGlobal;
    1900:     PetscInt *numExclusivelyOwnedAll;

    1902:     PetscCall(PetscMalloc1(cumSumVertices[rank + 1] - cumSumVertices[rank], &part));
    1903:     PetscCall(MatGetSize(A, &numRows, NULL));
    1904:     PetscCall(PetscLogEventBegin(DMPLEX_RebalGatherGraph, dm, 0, 0, 0));
    1905:     PetscCall(MatMPIAdjToSeqRankZero(A, &As));
    1906:     PetscCall(PetscLogEventEnd(DMPLEX_RebalGatherGraph, dm, 0, 0, 0));

    1908:     PetscCall(PetscMalloc1(size, &numExclusivelyOwnedAll));
    1909:     numExclusivelyOwnedAll[rank] = numExclusivelyOwned;
    1910:     PetscCallMPI(MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, numExclusivelyOwnedAll, 1, MPIU_INT, comm));

    1912:     PetscCall(PetscMalloc1(numRows, &partGlobal));
    1913:     PetscCall(PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0));
    1914:     if (rank == 0) {
    1915:       PetscInt *vtxwgt_g, numRows_g;

    1917:       PetscCall(MatGetRowIJ(As, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows_g, &xadj, &adjncy, &done));
    1918:       PetscCall(PetscMalloc1(2 * numRows_g, &vtxwgt_g));
    1919:       for (i = 0; i < size; i++) {
    1920:         vtxwgt_g[ncon * cumSumVertices[i]] = numExclusivelyOwnedAll[i];
    1921:         if (ncon > 1) vtxwgt_g[ncon * cumSumVertices[i] + 1] = 1;
    1922:         for (j = cumSumVertices[i] + 1; j < cumSumVertices[i + 1]; j++) {
    1923:           vtxwgt_g[ncon * j] = 1;
    1924:           if (ncon > 1) vtxwgt_g[2 * j + 1] = 0;
    1925:         }
    1926:       }

    1928:       PetscCall(PetscMalloc1(64, &options));
    1929:       ierr = METIS_SetDefaultOptions(options); /* initialize all defaults */
    1930:       PetscCheck(ierr == METIS_OK, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_SetDefaultOptions()");
    1931:       options[METIS_OPTION_CONTIG] = 1;
    1932:       PetscStackPushExternal("METIS_PartGraphKway");
    1933:       ierr = METIS_PartGraphKway(&numRows_g, &ncon, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt_g, NULL, NULL, &nparts, tpwgts, ubvec, options, &edgecut, partGlobal);
    1934:       PetscStackPop;
    1935:       PetscCheck(ierr == METIS_OK, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphKway()");
    1936:       PetscCall(PetscFree(options));
    1937:       PetscCall(PetscFree(vtxwgt_g));
    1938:       PetscCall(MatRestoreRowIJ(As, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows_g, &xadj, &adjncy, &done));
    1939:       PetscCall(MatDestroy(&As));
    1940:     }
    1941:     PetscCall(PetscBarrier((PetscObject)dm));
    1942:     PetscCall(PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0));
    1943:     PetscCall(PetscFree(numExclusivelyOwnedAll));

    1945:     /* scatter the partitioning information to ranks */
    1946:     PetscCall(PetscLogEventBegin(DMPLEX_RebalScatterPart, 0, 0, 0, 0));
    1947:     PetscCall(PetscMalloc1(size, &counts));
    1948:     PetscCall(PetscMalloc1(size + 1, &mpiCumSumVertices));
    1949:     for (i = 0; i < size; i++) PetscCall(PetscMPIIntCast(cumSumVertices[i + 1] - cumSumVertices[i], &counts[i]));
    1950:     for (i = 0; i <= size; i++) PetscCall(PetscMPIIntCast(cumSumVertices[i], &mpiCumSumVertices[i]));
    1951:     PetscCallMPI(MPI_Scatterv(partGlobal, counts, mpiCumSumVertices, MPIU_INT, part, counts[rank], MPIU_INT, 0, comm));
    1952:     PetscCall(PetscFree(counts));
    1953:     PetscCall(PetscFree(mpiCumSumVertices));
    1954:     PetscCall(PetscFree(partGlobal));
    1955:     PetscCall(PetscLogEventEnd(DMPLEX_RebalScatterPart, 0, 0, 0, 0));
    1956:   }
    1957:   PetscCall(PetscFree(ubvec));
    1958:   PetscCall(PetscFree(tpwgts));

    1960:   /* Rename the result so that the vertex resembling the exclusively owned points stays on the same rank */
    1961:   PetscCall(PetscMalloc2(size, &firstVertices, size, &renumbering));
    1962:   firstVertices[rank] = part[0];
    1963:   PetscCallMPI(MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, firstVertices, 1, MPIU_INT, comm));
    1964:   for (i = 0; i < size; i++) renumbering[firstVertices[i]] = i;
    1965:   for (i = 0; i < cumSumVertices[rank + 1] - cumSumVertices[rank]; i++) part[i] = renumbering[part[i]];
    1966:   PetscCall(PetscFree2(firstVertices, renumbering));

    1968:   /* Check if the renumbering worked (this can fail when ParMETIS gives fewer partitions than there are processes) */
    1969:   failed = (PetscInt)(part[0] != rank);
    1970:   PetscCallMPI(MPIU_Allreduce(&failed, &failedGlobal, 1, MPIU_INT, MPI_SUM, comm));
    1971:   if (failedGlobal > 0) {
    1972:     PetscCheck(failedGlobal <= 0, comm, PETSC_ERR_LIB, "Metis/Parmetis returned a bad partition");
    1973:     PetscCall(PetscFree(vtxwgt));
    1974:     PetscCall(PetscFree(toBalance));
    1975:     PetscCall(PetscFree(isLeaf));
    1976:     PetscCall(PetscFree(isNonExclusivelyOwned));
    1977:     PetscCall(PetscFree(isExclusivelyOwned));
    1978:     if (usematpartitioning) {
    1979:       PetscCall(ISRestoreIndices(ispart, (const PetscInt **)&part));
    1980:       PetscCall(ISDestroy(&ispart));
    1981:     } else PetscCall(PetscFree(part));
    1982:     if (viewer) {
    1983:       PetscCall(PetscViewerPopFormat(viewer));
    1984:       PetscCall(PetscViewerDestroy(&viewer));
    1985:     }
    1986:     PetscCall(PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0));
    1987:     PetscFunctionReturn(PETSC_SUCCESS);
    1988:   }

    1990:   /* Check how well we did distributing points*/
    1991:   if (viewer) {
    1992:     PetscCall(PetscViewerASCIIPrintf(viewer, "Number of owned entities of depth %" PetscInt_FMT ".\n", entityDepth));
    1993:     PetscCall(PetscViewerASCIIPrintf(viewer, "Initial      "));
    1994:     PetscCall(DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, NULL, viewer));
    1995:     PetscCall(PetscViewerASCIIPrintf(viewer, "Rebalanced   "));
    1996:     PetscCall(DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, part, viewer));
    1997:   }

    1999:   /* Check that every vertex is owned by a process that it is actually connected to. */
    2000:   PetscCall(MatGetRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done));
    2001:   for (i = 1; i <= numNonExclusivelyOwned; i++) {
    2002:     PetscInt loc = 0;
    2003:     PetscCall(PetscFindInt(cumSumVertices[part[i]], xadj[i + 1] - xadj[i], &adjncy[xadj[i]], &loc));
    2004:     /* If not, then just set the owner to the original owner (hopefully a rare event, it means that a vertex has been isolated) */
    2005:     if (loc < 0) part[i] = rank;
    2006:   }
    2007:   PetscCall(MatRestoreRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done));
    2008:   PetscCall(MatDestroy(&A));

    2010:   /* See how significant the influences of the previous fixing up step was.*/
    2011:   if (viewer) {
    2012:     PetscCall(PetscViewerASCIIPrintf(viewer, "After fix    "));
    2013:     PetscCall(DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, part, viewer));
    2014:   }
    2015:   if (!usematpartitioning) PetscCall(PetscFree(vtxwgt));
    2016:   else PetscCall(MatPartitioningDestroy(&mp));

    2018:   PetscCall(PetscLayoutDestroy(&layout));

    2020:   PetscCall(PetscLogEventBegin(DMPLEX_RebalRewriteSF, dm, 0, 0, 0));
    2021:   /* Rewrite the SF to reflect the new ownership. */
    2022:   PetscCall(PetscMalloc1(numNonExclusivelyOwned, &pointsToRewrite));
    2023:   counter = 0;
    2024:   for (i = 0; i < pEnd - pStart; i++) {
    2025:     if (toBalance[i]) {
    2026:       if (isNonExclusivelyOwned[i]) {
    2027:         pointsToRewrite[counter] = i + pStart;
    2028:         counter++;
    2029:       }
    2030:     }
    2031:   }
    2032:   PetscCall(DMPlexRewriteSF(dm, numNonExclusivelyOwned, pointsToRewrite, part + 1, degrees));
    2033:   PetscCall(PetscFree(pointsToRewrite));
    2034:   PetscCall(PetscLogEventEnd(DMPLEX_RebalRewriteSF, dm, 0, 0, 0));

    2036:   PetscCall(PetscFree(toBalance));
    2037:   PetscCall(PetscFree(isLeaf));
    2038:   PetscCall(PetscFree(isNonExclusivelyOwned));
    2039:   PetscCall(PetscFree(isExclusivelyOwned));
    2040:   if (usematpartitioning) {
    2041:     PetscCall(ISRestoreIndices(ispart, (const PetscInt **)&part));
    2042:     PetscCall(ISDestroy(&ispart));
    2043:   } else PetscCall(PetscFree(part));
    2044:   if (viewer) {
    2045:     PetscCall(PetscViewerPopFormat(viewer));
    2046:     PetscCall(PetscViewerDestroy(&viewer));
    2047:   }
    2048:   if (success) *success = PETSC_TRUE;
    2049:   PetscCall(PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0));
    2050:   PetscFunctionReturn(PETSC_SUCCESS);
    2051: #else
    2052:   SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
    2053: #endif
    2054: }