Actual source code: plexpartition.c

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

      5: PETSC_STATIC_INLINE PetscInt DMPlex_GlobalID(PetscInt point) { return point >= 0 ? point : -(point+1); }

      7: static PetscErrorCode DMPlexCreatePartitionerGraph_Native(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
      8: {
      9:   PetscInt       dim, depth, p, pStart, pEnd, a, adjSize, idx, size;
     10:   PetscInt      *adj = NULL, *vOffsets = NULL, *graph = NULL;
     11:   IS             cellNumbering;
     12:   const PetscInt *cellNum;
     13:   PetscBool      useCone, useClosure;
     14:   PetscSection   section;
     15:   PetscSegBuffer adjBuffer;
     16:   PetscSF        sfPoint;
     17:   PetscInt       *adjCells = NULL, *remoteCells = NULL;
     18:   const PetscInt *local;
     19:   PetscInt       nroots, nleaves, l;
     20:   PetscMPIInt    rank;

     24:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
     25:   DMGetDimension(dm, &dim);
     26:   DMPlexGetDepth(dm, &depth);
     27:   if (dim != depth) {
     28:     /* We do not handle the uninterpolated case here */
     29:     DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
     30:     /* DMPlexCreateNeighborCSR does not make a numbering */
     31:     if (globalNumbering) {DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);}
     32:     /* Different behavior for empty graphs */
     33:     if (!*numVertices) {
     34:       PetscMalloc1(1, offsets);
     35:       (*offsets)[0] = 0;
     36:     }
     37:     /* Broken in parallel */
     38:     if (rank && *numVertices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parallel partitioning of uninterpolated meshes not supported");
     39:     return(0);
     40:   }
     41:   DMGetPointSF(dm, &sfPoint);
     42:   DMPlexGetHeightStratum(dm, height, &pStart, &pEnd);
     43:   /* Build adjacency graph via a section/segbuffer */
     44:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &section);
     45:   PetscSectionSetChart(section, pStart, pEnd);
     46:   PetscSegBufferCreate(sizeof(PetscInt),1000,&adjBuffer);
     47:   /* Always use FVM adjacency to create partitioner graph */
     48:   DMGetBasicAdjacency(dm, &useCone, &useClosure);
     49:   DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE);
     50:   DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, sfPoint, &cellNumbering);
     51:   if (globalNumbering) {
     52:     PetscObjectReference((PetscObject)cellNumbering);
     53:     *globalNumbering = cellNumbering;
     54:   }
     55:   ISGetIndices(cellNumbering, &cellNum);
     56:   /* For all boundary faces (including faces adjacent to a ghost cell), record the local cell in adjCells
     57:      Broadcast adjCells to remoteCells (to get cells from roots) and Reduce adjCells to remoteCells (to get cells from leaves)
     58:    */
     59:   PetscSFGetGraph(sfPoint, &nroots, &nleaves, &local, NULL);
     60:   if (nroots >= 0) {
     61:     PetscInt fStart, fEnd, f;

     63:     PetscCalloc2(nroots, &adjCells, nroots, &remoteCells);
     64:     DMPlexGetHeightStratum(dm, height+1, &fStart, &fEnd);
     65:     for (l = 0; l < nroots; ++l) adjCells[l] = -3;
     66:     for (f = fStart; f < fEnd; ++f) {
     67:       const PetscInt *support;
     68:       PetscInt        supportSize;

     70:       DMPlexGetSupport(dm, f, &support);
     71:       DMPlexGetSupportSize(dm, f, &supportSize);
     72:       if (supportSize == 1) adjCells[f] = DMPlex_GlobalID(cellNum[support[0]]);
     73:       else if (supportSize == 2) {
     74:         PetscFindInt(support[0], nleaves, local, &p);
     75:         if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[1]]);
     76:         PetscFindInt(support[1], nleaves, local, &p);
     77:         if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[0]]);
     78:       }
     79:       /* Handle non-conforming meshes */
     80:       if (supportSize > 2) {
     81:         PetscInt        numChildren, i;
     82:         const PetscInt *children;

     84:         DMPlexGetTreeChildren(dm, f, &numChildren, &children);
     85:         for (i = 0; i < numChildren; ++i) {
     86:           const PetscInt child = children[i];
     87:           if (fStart <= child && child < fEnd) {
     88:             DMPlexGetSupport(dm, child, &support);
     89:             DMPlexGetSupportSize(dm, child, &supportSize);
     90:             if (supportSize == 1) adjCells[child] = DMPlex_GlobalID(cellNum[support[0]]);
     91:             else if (supportSize == 2) {
     92:               PetscFindInt(support[0], nleaves, local, &p);
     93:               if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[1]]);
     94:               PetscFindInt(support[1], nleaves, local, &p);
     95:               if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[0]]);
     96:             }
     97:           }
     98:         }
     99:       }
    100:     }
    101:     for (l = 0; l < nroots; ++l) remoteCells[l] = -1;
    102:     PetscSFBcastBegin(dm->sf, MPIU_INT, adjCells, remoteCells,MPI_REPLACE);
    103:     PetscSFBcastEnd(dm->sf, MPIU_INT, adjCells, remoteCells,MPI_REPLACE);
    104:     PetscSFReduceBegin(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
    105:     PetscSFReduceEnd(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
    106:   }
    107:   /* Combine local and global adjacencies */
    108:   for (*numVertices = 0, p = pStart; p < pEnd; p++) {
    109:     /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
    110:     if (nroots > 0) {if (cellNum[p] < 0) continue;}
    111:     /* Add remote cells */
    112:     if (remoteCells) {
    113:       const PetscInt gp = DMPlex_GlobalID(cellNum[p]);
    114:       PetscInt       coneSize, numChildren, c, i;
    115:       const PetscInt *cone, *children;

    117:       DMPlexGetCone(dm, p, &cone);
    118:       DMPlexGetConeSize(dm, p, &coneSize);
    119:       for (c = 0; c < coneSize; ++c) {
    120:         const PetscInt point = cone[c];
    121:         if (remoteCells[point] >= 0 && remoteCells[point] != gp) {
    122:           PetscInt *PETSC_RESTRICT pBuf;
    123:           PetscSectionAddDof(section, p, 1);
    124:           PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
    125:           *pBuf = remoteCells[point];
    126:         }
    127:         /* Handle non-conforming meshes */
    128:         DMPlexGetTreeChildren(dm, point, &numChildren, &children);
    129:         for (i = 0; i < numChildren; ++i) {
    130:           const PetscInt child = children[i];
    131:           if (remoteCells[child] >= 0 && remoteCells[child] != gp) {
    132:             PetscInt *PETSC_RESTRICT pBuf;
    133:             PetscSectionAddDof(section, p, 1);
    134:             PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
    135:             *pBuf = remoteCells[child];
    136:           }
    137:         }
    138:       }
    139:     }
    140:     /* Add local cells */
    141:     adjSize = PETSC_DETERMINE;
    142:     DMPlexGetAdjacency(dm, p, &adjSize, &adj);
    143:     for (a = 0; a < adjSize; ++a) {
    144:       const PetscInt point = adj[a];
    145:       if (point != p && pStart <= point && point < pEnd) {
    146:         PetscInt *PETSC_RESTRICT pBuf;
    147:         PetscSectionAddDof(section, p, 1);
    148:         PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
    149:         *pBuf = DMPlex_GlobalID(cellNum[point]);
    150:       }
    151:     }
    152:     (*numVertices)++;
    153:   }
    154:   PetscFree(adj);
    155:   PetscFree2(adjCells, remoteCells);
    156:   DMSetBasicAdjacency(dm, useCone, useClosure);

    158:   /* Derive CSR graph from section/segbuffer */
    159:   PetscSectionSetUp(section);
    160:   PetscSectionGetStorageSize(section, &size);
    161:   PetscMalloc1(*numVertices+1, &vOffsets);
    162:   for (idx = 0, p = pStart; p < pEnd; p++) {
    163:     if (nroots > 0) {if (cellNum[p] < 0) continue;}
    164:     PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
    165:   }
    166:   vOffsets[*numVertices] = size;
    167:   PetscSegBufferExtractAlloc(adjBuffer, &graph);

    169:   if (nroots >= 0) {
    170:     /* Filter out duplicate edges using section/segbuffer */
    171:     PetscSectionReset(section);
    172:     PetscSectionSetChart(section, 0, *numVertices);
    173:     for (p = 0; p < *numVertices; p++) {
    174:       PetscInt start = vOffsets[p], end = vOffsets[p+1];
    175:       PetscInt numEdges = end-start, *PETSC_RESTRICT edges;
    176:       PetscSortRemoveDupsInt(&numEdges, &graph[start]);
    177:       PetscSectionSetDof(section, p, numEdges);
    178:       PetscSegBufferGetInts(adjBuffer, numEdges, &edges);
    179:       PetscArraycpy(edges, &graph[start], numEdges);
    180:     }
    181:     PetscFree(vOffsets);
    182:     PetscFree(graph);
    183:     /* Derive CSR graph from section/segbuffer */
    184:     PetscSectionSetUp(section);
    185:     PetscSectionGetStorageSize(section, &size);
    186:     PetscMalloc1(*numVertices+1, &vOffsets);
    187:     for (idx = 0, p = 0; p < *numVertices; p++) {
    188:       PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
    189:     }
    190:     vOffsets[*numVertices] = size;
    191:     PetscSegBufferExtractAlloc(adjBuffer, &graph);
    192:   } else {
    193:     /* Sort adjacencies (not strictly necessary) */
    194:     for (p = 0; p < *numVertices; p++) {
    195:       PetscInt start = vOffsets[p], end = vOffsets[p+1];
    196:       PetscSortInt(end-start, &graph[start]);
    197:     }
    198:   }

    200:   if (offsets) *offsets = vOffsets;
    201:   if (adjacency) *adjacency = graph;

    203:   /* Cleanup */
    204:   PetscSegBufferDestroy(&adjBuffer);
    205:   PetscSectionDestroy(&section);
    206:   ISRestoreIndices(cellNumbering, &cellNum);
    207:   ISDestroy(&cellNumbering);
    208:   return(0);
    209: }

    211: static PetscErrorCode DMPlexCreatePartitionerGraph_ViaMat(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
    212: {
    213:   Mat            conn, CSR;
    214:   IS             fis, cis, cis_own;
    215:   PetscSF        sfPoint;
    216:   const PetscInt *rows, *cols, *ii, *jj;
    217:   PetscInt       *idxs,*idxs2;
    218:   PetscInt       dim, depth, floc, cloc, i, M, N, c, lm, m, cStart, cEnd, fStart, fEnd;
    219:   PetscMPIInt    rank;
    220:   PetscBool      flg;

    224:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
    225:   DMGetDimension(dm, &dim);
    226:   DMPlexGetDepth(dm, &depth);
    227:   if (dim != depth) {
    228:     /* We do not handle the uninterpolated case here */
    229:     DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
    230:     /* DMPlexCreateNeighborCSR does not make a numbering */
    231:     if (globalNumbering) {DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);}
    232:     /* Different behavior for empty graphs */
    233:     if (!*numVertices) {
    234:       PetscMalloc1(1, offsets);
    235:       (*offsets)[0] = 0;
    236:     }
    237:     /* Broken in parallel */
    238:     if (rank && *numVertices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parallel partitioning of uninterpolated meshes not supported");
    239:     return(0);
    240:   }
    241:   /* Interpolated and parallel case */

    243:   /* numbering */
    244:   DMGetPointSF(dm, &sfPoint);
    245:   DMPlexGetHeightStratum(dm, height, &cStart, &cEnd);
    246:   DMPlexGetHeightStratum(dm, height+1, &fStart, &fEnd);
    247:   DMPlexCreateNumbering_Plex(dm, cStart, cEnd, 0, &N, sfPoint, &cis);
    248:   DMPlexCreateNumbering_Plex(dm, fStart, fEnd, 0, &M, sfPoint, &fis);
    249:   if (globalNumbering) {
    250:     ISDuplicate(cis, globalNumbering);
    251:   }

    253:   /* get positive global ids and local sizes for facets and cells */
    254:   ISGetLocalSize(fis, &m);
    255:   ISGetIndices(fis, &rows);
    256:   PetscMalloc1(m, &idxs);
    257:   for (i = 0, floc = 0; i < m; i++) {
    258:     const PetscInt p = rows[i];

    260:     if (p < 0) {
    261:       idxs[i] = -(p+1);
    262:     } else {
    263:       idxs[i] = p;
    264:       floc   += 1;
    265:     }
    266:   }
    267:   ISRestoreIndices(fis, &rows);
    268:   ISDestroy(&fis);
    269:   ISCreateGeneral(PETSC_COMM_SELF, m, idxs, PETSC_OWN_POINTER, &fis);

    271:   ISGetLocalSize(cis, &m);
    272:   ISGetIndices(cis, &cols);
    273:   PetscMalloc1(m, &idxs);
    274:   PetscMalloc1(m, &idxs2);
    275:   for (i = 0, cloc = 0; i < m; i++) {
    276:     const PetscInt p = cols[i];

    278:     if (p < 0) {
    279:       idxs[i] = -(p+1);
    280:     } else {
    281:       idxs[i]       = p;
    282:       idxs2[cloc++] = p;
    283:     }
    284:   }
    285:   ISRestoreIndices(cis, &cols);
    286:   ISDestroy(&cis);
    287:   ISCreateGeneral(PetscObjectComm((PetscObject)dm), m, idxs, PETSC_OWN_POINTER, &cis);
    288:   ISCreateGeneral(PetscObjectComm((PetscObject)dm), cloc, idxs2, PETSC_OWN_POINTER, &cis_own);

    290:   /* Create matrix to hold F-C connectivity (MatMatTranspose Mult not supported for MPIAIJ) */
    291:   MatCreate(PetscObjectComm((PetscObject)dm), &conn);
    292:   MatSetSizes(conn, floc, cloc, M, N);
    293:   MatSetType(conn, MATMPIAIJ);
    294:   DMPlexGetMaxSizes(dm, NULL, &lm);
    295:   MPI_Allreduce(&lm, &m, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) dm));
    296:   MatMPIAIJSetPreallocation(conn, m, NULL, m, NULL);

    298:   /* Assemble matrix */
    299:   ISGetIndices(fis, &rows);
    300:   ISGetIndices(cis, &cols);
    301:   for (c = cStart; c < cEnd; c++) {
    302:     const PetscInt *cone;
    303:     PetscInt        coneSize, row, col, f;

    305:     col  = cols[c-cStart];
    306:     DMPlexGetCone(dm, c, &cone);
    307:     DMPlexGetConeSize(dm, c, &coneSize);
    308:     for (f = 0; f < coneSize; f++) {
    309:       const PetscScalar v = 1.0;
    310:       const PetscInt *children;
    311:       PetscInt        numChildren, ch;

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

    316:       /* non-conforming meshes */
    317:       DMPlexGetTreeChildren(dm, cone[f], &numChildren, &children);
    318:       for (ch = 0; ch < numChildren; ch++) {
    319:         const PetscInt child = children[ch];

    321:         if (child < fStart || child >= fEnd) continue;
    322:         row  = rows[child-fStart];
    323:         MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES);
    324:       }
    325:     }
    326:   }
    327:   ISRestoreIndices(fis, &rows);
    328:   ISRestoreIndices(cis, &cols);
    329:   ISDestroy(&fis);
    330:   ISDestroy(&cis);
    331:   MatAssemblyBegin(conn, MAT_FINAL_ASSEMBLY);
    332:   MatAssemblyEnd(conn, MAT_FINAL_ASSEMBLY);

    334:   /* Get parallel CSR by doing conn^T * conn */
    335:   MatTransposeMatMult(conn, conn, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CSR);
    336:   MatDestroy(&conn);

    338:   /* extract local part of the CSR */
    339:   MatMPIAIJGetLocalMat(CSR, MAT_INITIAL_MATRIX, &conn);
    340:   MatDestroy(&CSR);
    341:   MatGetRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
    342:   if (!flg) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No IJ format");

    344:   /* get back requested output */
    345:   if (numVertices) *numVertices = m;
    346:   if (offsets) {
    347:     PetscCalloc1(m+1, &idxs);
    348:     for (i = 1; i < m+1; i++) idxs[i] = ii[i] - i; /* ParMetis does not like self-connectivity */
    349:     *offsets = idxs;
    350:   }
    351:   if (adjacency) {
    352:     PetscMalloc1(ii[m] - m, &idxs);
    353:     ISGetIndices(cis_own, &rows);
    354:     for (i = 0, c = 0; i < m; i++) {
    355:       PetscInt j, g = rows[i];

    357:       for (j = ii[i]; j < ii[i+1]; j++) {
    358:         if (jj[j] == g) continue; /* again, self-connectivity */
    359:         idxs[c++] = jj[j];
    360:       }
    361:     }
    362:     if (c != ii[m] - m) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected %D != %D",c,ii[m]-m);
    363:     ISRestoreIndices(cis_own, &rows);
    364:     *adjacency = idxs;
    365:   }

    367:   /* cleanup */
    368:   ISDestroy(&cis_own);
    369:   MatRestoreRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
    370:   if (!flg) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No IJ format");
    371:   MatDestroy(&conn);
    372:   return(0);
    373: }

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

    378:   Input Parameters:
    379: + dm      - The mesh DM dm
    380: - height  - Height of the strata from which to construct the graph

    382:   Output Parameter:
    383: + numVertices     - Number of vertices in the graph
    384: . offsets         - Point offsets in the graph
    385: . adjacency       - Point connectivity in the graph
    386: - 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.

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

    391:   Level: developer

    393: .seealso: PetscPartitionerGetType(), PetscPartitionerCreate(), DMSetAdjacency()
    394: @*/
    395: PetscErrorCode DMPlexCreatePartitionerGraph(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
    396: {
    398:   PetscBool      usemat = PETSC_FALSE;

    401:   PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_csr_via_mat", &usemat, NULL);
    402:   if (usemat) {
    403:     DMPlexCreatePartitionerGraph_ViaMat(dm, height, numVertices, offsets, adjacency, globalNumbering);
    404:   } else {
    405:     DMPlexCreatePartitionerGraph_Native(dm, height, numVertices, offsets, adjacency, globalNumbering);
    406:   }
    407:   return(0);
    408: }

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

    413:   Collective on DM

    415:   Input Arguments:
    416: + dm - The DMPlex
    417: - cellHeight - The height of mesh points to treat as cells (default should be 0)

    419:   Output Arguments:
    420: + numVertices - The number of local vertices in the graph, or cells in the mesh.
    421: . offsets     - The offset to the adjacency list for each cell
    422: - adjacency   - The adjacency list for all cells

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

    426:   Level: advanced

    428: .seealso: DMPlexCreate()
    429: @*/
    430: PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
    431: {
    432:   const PetscInt maxFaceCases = 30;
    433:   PetscInt       numFaceCases = 0;
    434:   PetscInt       numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */
    435:   PetscInt      *off, *adj;
    436:   PetscInt      *neighborCells = NULL;
    437:   PetscInt       dim, cellDim, depth = 0, faceDepth, cStart, cEnd, c, numCells, cell;

    441:   /* For parallel partitioning, I think you have to communicate supports */
    442:   DMGetDimension(dm, &dim);
    443:   cellDim = dim - cellHeight;
    444:   DMPlexGetDepth(dm, &depth);
    445:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
    446:   if (cEnd - cStart == 0) {
    447:     if (numVertices) *numVertices = 0;
    448:     if (offsets)   *offsets   = NULL;
    449:     if (adjacency) *adjacency = NULL;
    450:     return(0);
    451:   }
    452:   numCells  = cEnd - cStart;
    453:   faceDepth = depth - cellHeight;
    454:   if (dim == depth) {
    455:     PetscInt f, fStart, fEnd;

    457:     PetscCalloc1(numCells+1, &off);
    458:     /* Count neighboring cells */
    459:     DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
    460:     for (f = fStart; f < fEnd; ++f) {
    461:       const PetscInt *support;
    462:       PetscInt        supportSize;
    463:       DMPlexGetSupportSize(dm, f, &supportSize);
    464:       DMPlexGetSupport(dm, f, &support);
    465:       if (supportSize == 2) {
    466:         PetscInt numChildren;

    468:         DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
    469:         if (!numChildren) {
    470:           ++off[support[0]-cStart+1];
    471:           ++off[support[1]-cStart+1];
    472:         }
    473:       }
    474:     }
    475:     /* Prefix sum */
    476:     for (c = 1; c <= numCells; ++c) off[c] += off[c-1];
    477:     if (adjacency) {
    478:       PetscInt *tmp;

    480:       PetscMalloc1(off[numCells], &adj);
    481:       PetscMalloc1(numCells+1, &tmp);
    482:       PetscArraycpy(tmp, off, numCells+1);
    483:       /* Get neighboring cells */
    484:       for (f = fStart; f < fEnd; ++f) {
    485:         const PetscInt *support;
    486:         PetscInt        supportSize;
    487:         DMPlexGetSupportSize(dm, f, &supportSize);
    488:         DMPlexGetSupport(dm, f, &support);
    489:         if (supportSize == 2) {
    490:           PetscInt numChildren;

    492:           DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
    493:           if (!numChildren) {
    494:             adj[tmp[support[0]-cStart]++] = support[1];
    495:             adj[tmp[support[1]-cStart]++] = support[0];
    496:           }
    497:         }
    498:       }
    499:       if (PetscDefined(USE_DEBUG)) {
    500:         for (c = 0; c < cEnd-cStart; ++c) if (tmp[c] != off[c+1]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Offset %d != %d for cell %d", tmp[c], off[c], c+cStart);
    501:       }
    502:       PetscFree(tmp);
    503:     }
    504:     if (numVertices) *numVertices = numCells;
    505:     if (offsets)   *offsets   = off;
    506:     if (adjacency) *adjacency = adj;
    507:     return(0);
    508:   }
    509:   /* Setup face recognition */
    510:   if (faceDepth == 1) {
    511:     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 */

    513:     for (c = cStart; c < cEnd; ++c) {
    514:       PetscInt corners;

    516:       DMPlexGetConeSize(dm, c, &corners);
    517:       if (!cornersSeen[corners]) {
    518:         PetscInt nFV;

    520:         if (numFaceCases >= maxFaceCases) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Exceeded maximum number of face recognition cases");
    521:         cornersSeen[corners] = 1;

    523:         DMPlexGetNumFaceVertices(dm, cellDim, corners, &nFV);

    525:         numFaceVertices[numFaceCases++] = nFV;
    526:       }
    527:     }
    528:   }
    529:   PetscCalloc1(numCells+1, &off);
    530:   /* Count neighboring cells */
    531:   for (cell = cStart; cell < cEnd; ++cell) {
    532:     PetscInt numNeighbors = PETSC_DETERMINE, n;

    534:     DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
    535:     /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
    536:     for (n = 0; n < numNeighbors; ++n) {
    537:       PetscInt        cellPair[2];
    538:       PetscBool       found    = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
    539:       PetscInt        meetSize = 0;
    540:       const PetscInt *meet    = NULL;

    542:       cellPair[0] = cell; cellPair[1] = neighborCells[n];
    543:       if (cellPair[0] == cellPair[1]) continue;
    544:       if (!found) {
    545:         DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
    546:         if (meetSize) {
    547:           PetscInt f;

    549:           for (f = 0; f < numFaceCases; ++f) {
    550:             if (numFaceVertices[f] == meetSize) {
    551:               found = PETSC_TRUE;
    552:               break;
    553:             }
    554:           }
    555:         }
    556:         DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
    557:       }
    558:       if (found) ++off[cell-cStart+1];
    559:     }
    560:   }
    561:   /* Prefix sum */
    562:   for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell-1];

    564:   if (adjacency) {
    565:     PetscMalloc1(off[numCells], &adj);
    566:     /* Get neighboring cells */
    567:     for (cell = cStart; cell < cEnd; ++cell) {
    568:       PetscInt numNeighbors = PETSC_DETERMINE, n;
    569:       PetscInt cellOffset   = 0;

    571:       DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
    572:       /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
    573:       for (n = 0; n < numNeighbors; ++n) {
    574:         PetscInt        cellPair[2];
    575:         PetscBool       found    = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
    576:         PetscInt        meetSize = 0;
    577:         const PetscInt *meet    = NULL;

    579:         cellPair[0] = cell; cellPair[1] = neighborCells[n];
    580:         if (cellPair[0] == cellPair[1]) continue;
    581:         if (!found) {
    582:           DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
    583:           if (meetSize) {
    584:             PetscInt f;

    586:             for (f = 0; f < numFaceCases; ++f) {
    587:               if (numFaceVertices[f] == meetSize) {
    588:                 found = PETSC_TRUE;
    589:                 break;
    590:               }
    591:             }
    592:           }
    593:           DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
    594:         }
    595:         if (found) {
    596:           adj[off[cell-cStart]+cellOffset] = neighborCells[n];
    597:           ++cellOffset;
    598:         }
    599:       }
    600:     }
    601:   }
    602:   PetscFree(neighborCells);
    603:   if (numVertices) *numVertices = numCells;
    604:   if (offsets)   *offsets   = off;
    605:   if (adjacency) *adjacency = adj;
    606:   return(0);
    607: }

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

    612:   Collective on PetscPartitioner

    614:   Input Parameters:
    615: + part    - The PetscPartitioner
    616: . targetSection - The PetscSection describing the absolute weight of each partition (can be NULL)
    617: - dm      - The mesh DM

    619:   Output Parameters:
    620: + partSection     - The PetscSection giving the division of points by partition
    621: - partition       - The list of points by partition

    623:   Notes:
    624:     If the DM has a local section associated, each point to be partitioned will be weighted by the total number of dofs identified
    625:     by the section in the transitive closure of the point.

    627:   Level: developer

    629: .seealso DMPlexDistribute(), PetscPartitionerCreate(), PetscSectionCreate(), PetscSectionSetChart(), PetscPartitionerPartition()
    630: @*/
    631: PetscErrorCode PetscPartitionerDMPlexPartition(PetscPartitioner part, DM dm, PetscSection targetSection, PetscSection partSection, IS *partition)
    632: {
    633:   PetscMPIInt    size;
    634:   PetscBool      isplex;
    636:   PetscSection   vertSection = NULL;

    644:   PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isplex);
    645:   if (!isplex) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not for type %s",((PetscObject)dm)->type_name);
    646:   MPI_Comm_size(PetscObjectComm((PetscObject) part), &size);
    647:   if (size == 1) {
    648:     PetscInt *points;
    649:     PetscInt  cStart, cEnd, c;

    651:     DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
    652:     PetscSectionReset(partSection);
    653:     PetscSectionSetChart(partSection, 0, size);
    654:     PetscSectionSetDof(partSection, 0, cEnd-cStart);
    655:     PetscSectionSetUp(partSection);
    656:     PetscMalloc1(cEnd-cStart, &points);
    657:     for (c = cStart; c < cEnd; ++c) points[c] = c;
    658:     ISCreateGeneral(PetscObjectComm((PetscObject) part), cEnd-cStart, points, PETSC_OWN_POINTER, partition);
    659:     return(0);
    660:   }
    661:   if (part->height == 0) {
    662:     PetscInt numVertices = 0;
    663:     PetscInt *start     = NULL;
    664:     PetscInt *adjacency = NULL;
    665:     IS       globalNumbering;

    667:     if (!part->noGraph || part->viewGraph) {
    668:       DMPlexCreatePartitionerGraph(dm, part->height, &numVertices, &start, &adjacency, &globalNumbering);
    669:     } else { /* only compute the number of owned local vertices */
    670:       const PetscInt *idxs;
    671:       PetscInt       p, pStart, pEnd;

    673:       DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
    674:       DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, dm->sf, &globalNumbering);
    675:       ISGetIndices(globalNumbering, &idxs);
    676:       for (p = 0; p < pEnd - pStart; p++) numVertices += idxs[p] < 0 ? 0 : 1;
    677:       ISRestoreIndices(globalNumbering, &idxs);
    678:     }
    679:     if (part->usevwgt) {
    680:       PetscSection   section = dm->localSection, clSection = NULL;
    681:       IS             clPoints = NULL;
    682:       const PetscInt *gid,*clIdx;
    683:       PetscInt       v, p, pStart, pEnd;

    685:       /* 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) */
    686:       /* We do this only if the local section has been set */
    687:       if (section) {
    688:         PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, NULL);
    689:         if (!clSection) {
    690:           DMPlexCreateClosureIndex(dm,NULL);
    691:         }
    692:         PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, &clPoints);
    693:         ISGetIndices(clPoints,&clIdx);
    694:       }
    695:       DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
    696:       PetscSectionCreate(PETSC_COMM_SELF, &vertSection);
    697:       PetscSectionSetChart(vertSection, 0, numVertices);
    698:       if (globalNumbering) {
    699:         ISGetIndices(globalNumbering,&gid);
    700:       } else gid = NULL;
    701:       for (p = pStart, v = 0; p < pEnd; ++p) {
    702:         PetscInt dof = 1;

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

    707:         if (section) {
    708:           PetscInt cl, clSize, clOff;

    710:           dof  = 0;
    711:           PetscSectionGetDof(clSection, p, &clSize);
    712:           PetscSectionGetOffset(clSection, p, &clOff);
    713:           for (cl = 0; cl < clSize; cl+=2) {
    714:             PetscInt clDof, clPoint = clIdx[clOff + cl]; /* odd indices are reserved for orientations */

    716:             PetscSectionGetDof(section, clPoint, &clDof);
    717:             dof += clDof;
    718:           }
    719:         }
    720:         if (!dof) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of dofs for point %D in the local section should be positive",p);
    721:         PetscSectionSetDof(vertSection, v, dof);
    722:         v++;
    723:       }
    724:       if (globalNumbering) {
    725:         ISRestoreIndices(globalNumbering,&gid);
    726:       }
    727:       if (clPoints) {
    728:         ISRestoreIndices(clPoints,&clIdx);
    729:       }
    730:       PetscSectionSetUp(vertSection);
    731:     }
    732:     PetscPartitionerPartition(part, size, numVertices, start, adjacency, vertSection, targetSection, partSection, partition);
    733:     PetscFree(start);
    734:     PetscFree(adjacency);
    735:     if (globalNumbering) { /* partition is wrt global unique numbering: change this to be wrt local numbering */
    736:       const PetscInt *globalNum;
    737:       const PetscInt *partIdx;
    738:       PetscInt       *map, cStart, cEnd;
    739:       PetscInt       *adjusted, i, localSize, offset;
    740:       IS             newPartition;

    742:       ISGetLocalSize(*partition,&localSize);
    743:       PetscMalloc1(localSize,&adjusted);
    744:       ISGetIndices(globalNumbering,&globalNum);
    745:       ISGetIndices(*partition,&partIdx);
    746:       PetscMalloc1(localSize,&map);
    747:       DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
    748:       for (i = cStart, offset = 0; i < cEnd; i++) {
    749:         if (globalNum[i - cStart] >= 0) map[offset++] = i;
    750:       }
    751:       for (i = 0; i < localSize; i++) {
    752:         adjusted[i] = map[partIdx[i]];
    753:       }
    754:       PetscFree(map);
    755:       ISRestoreIndices(*partition,&partIdx);
    756:       ISRestoreIndices(globalNumbering,&globalNum);
    757:       ISCreateGeneral(PETSC_COMM_SELF,localSize,adjusted,PETSC_OWN_POINTER,&newPartition);
    758:       ISDestroy(&globalNumbering);
    759:       ISDestroy(partition);
    760:       *partition = newPartition;
    761:     }
    762:   } else SETERRQ1(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_OUTOFRANGE, "Invalid height %D for points to partition", part->height);
    763:   PetscSectionDestroy(&vertSection);
    764:   return(0);
    765: }

    767: /*@
    768:   DMPlexGetPartitioner - Get the mesh partitioner

    770:   Not collective

    772:   Input Parameter:
    773: . dm - The DM

    775:   Output Parameter:
    776: . part - The PetscPartitioner

    778:   Level: developer

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

    782: .seealso DMPlexDistribute(), DMPlexSetPartitioner(), PetscPartitionerDMPlexPartition(), PetscPartitionerCreate()
    783: @*/
    784: PetscErrorCode DMPlexGetPartitioner(DM dm, PetscPartitioner *part)
    785: {
    786:   DM_Plex       *mesh = (DM_Plex *) dm->data;

    791:   *part = mesh->partitioner;
    792:   return(0);
    793: }

    795: /*@
    796:   DMPlexSetPartitioner - Set the mesh partitioner

    798:   logically collective on DM

    800:   Input Parameters:
    801: + dm - The DM
    802: - part - The partitioner

    804:   Level: developer

    806:   Note: Any existing PetscPartitioner will be destroyed.

    808: .seealso DMPlexDistribute(), DMPlexGetPartitioner(), PetscPartitionerCreate()
    809: @*/
    810: PetscErrorCode DMPlexSetPartitioner(DM dm, PetscPartitioner part)
    811: {
    812:   DM_Plex       *mesh = (DM_Plex *) dm->data;

    818:   PetscObjectReference((PetscObject)part);
    819:   PetscPartitionerDestroy(&mesh->partitioner);
    820:   mesh->partitioner = part;
    821:   return(0);
    822: }

    824: static PetscErrorCode DMPlexAddClosure_Private(DM dm, PetscHSetI ht, PetscInt point)
    825: {
    826:   const PetscInt *cone;
    827:   PetscInt       coneSize, c;
    828:   PetscBool      missing;

    832:   PetscHSetIQueryAdd(ht, point, &missing);
    833:   if (missing) {
    834:     DMPlexGetCone(dm, point, &cone);
    835:     DMPlexGetConeSize(dm, point, &coneSize);
    836:     for (c = 0; c < coneSize; c++) {
    837:       DMPlexAddClosure_Private(dm, ht, cone[c]);
    838:     }
    839:   }
    840:   return(0);
    841: }

    843: PETSC_UNUSED static PetscErrorCode DMPlexAddClosure_Tree(DM dm, PetscHSetI ht, PetscInt point, PetscBool up, PetscBool down)
    844: {

    848:   if (up) {
    849:     PetscInt parent;

    851:     DMPlexGetTreeParent(dm,point,&parent,NULL);
    852:     if (parent != point) {
    853:       PetscInt closureSize, *closure = NULL, i;

    855:       DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
    856:       for (i = 0; i < closureSize; i++) {
    857:         PetscInt cpoint = closure[2*i];

    859:         PetscHSetIAdd(ht, cpoint);
    860:         DMPlexAddClosure_Tree(dm,ht,cpoint,PETSC_TRUE,PETSC_FALSE);
    861:       }
    862:       DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
    863:     }
    864:   }
    865:   if (down) {
    866:     PetscInt numChildren;
    867:     const PetscInt *children;

    869:     DMPlexGetTreeChildren(dm,point,&numChildren,&children);
    870:     if (numChildren) {
    871:       PetscInt i;

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

    876:         PetscHSetIAdd(ht, cpoint);
    877:         DMPlexAddClosure_Tree(dm,ht,cpoint,PETSC_FALSE,PETSC_TRUE);
    878:       }
    879:     }
    880:   }
    881:   return(0);
    882: }

    884: static PetscErrorCode DMPlexAddClosureTree_Up_Private(DM dm, PetscHSetI ht, PetscInt point)
    885: {
    886:   PetscInt       parent;

    890:   DMPlexGetTreeParent(dm, point, &parent,NULL);
    891:   if (point != parent) {
    892:     const PetscInt *cone;
    893:     PetscInt       coneSize, c;

    895:     DMPlexAddClosureTree_Up_Private(dm, ht, parent);
    896:     DMPlexAddClosure_Private(dm, ht, parent);
    897:     DMPlexGetCone(dm, parent, &cone);
    898:     DMPlexGetConeSize(dm, parent, &coneSize);
    899:     for (c = 0; c < coneSize; c++) {
    900:       const PetscInt cp = cone[c];

    902:       DMPlexAddClosureTree_Up_Private(dm, ht, cp);
    903:     }
    904:   }
    905:   return(0);
    906: }

    908: static PetscErrorCode DMPlexAddClosureTree_Down_Private(DM dm, PetscHSetI ht, PetscInt point)
    909: {
    910:   PetscInt       i, numChildren;
    911:   const PetscInt *children;

    915:   DMPlexGetTreeChildren(dm, point, &numChildren, &children);
    916:   for (i = 0; i < numChildren; i++) {
    917:     PetscHSetIAdd(ht, children[i]);
    918:   }
    919:   return(0);
    920: }

    922: static PetscErrorCode DMPlexAddClosureTree_Private(DM dm, PetscHSetI ht, PetscInt point)
    923: {
    924:   const PetscInt *cone;
    925:   PetscInt       coneSize, c;

    929:   PetscHSetIAdd(ht, point);
    930:   DMPlexAddClosureTree_Up_Private(dm, ht, point);
    931:   DMPlexAddClosureTree_Down_Private(dm, ht, point);
    932:   DMPlexGetCone(dm, point, &cone);
    933:   DMPlexGetConeSize(dm, point, &coneSize);
    934:   for (c = 0; c < coneSize; c++) {
    935:     DMPlexAddClosureTree_Private(dm, ht, cone[c]);
    936:   }
    937:   return(0);
    938: }

    940: PetscErrorCode DMPlexClosurePoints_Private(DM dm, PetscInt numPoints, const PetscInt points[], IS *closureIS)
    941: {
    942:   DM_Plex         *mesh = (DM_Plex *)dm->data;
    943:   const PetscBool hasTree = (mesh->parentSection || mesh->childSection) ? PETSC_TRUE : PETSC_FALSE;
    944:   PetscInt        nelems, *elems, off = 0, p;
    945:   PetscHSetI      ht;
    946:   PetscErrorCode  ierr;

    949:   PetscHSetICreate(&ht);
    950:   PetscHSetIResize(ht, numPoints*16);
    951:   if (!hasTree) {
    952:     for (p = 0; p < numPoints; ++p) {
    953:       DMPlexAddClosure_Private(dm, ht, points[p]);
    954:     }
    955:   } else {
    956: #if 1
    957:     for (p = 0; p < numPoints; ++p) {
    958:       DMPlexAddClosureTree_Private(dm, ht, points[p]);
    959:     }
    960: #else
    961:     PetscInt  *closure = NULL, closureSize, c;
    962:     for (p = 0; p < numPoints; ++p) {
    963:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
    964:       for (c = 0; c < closureSize*2; c += 2) {
    965:         PetscHSetIAdd(ht, closure[c]);
    966:         if (hasTree) {DMPlexAddClosure_Tree(dm, ht, closure[c], PETSC_TRUE, PETSC_TRUE);}
    967:       }
    968:     }
    969:     if (closure) {DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure);}
    970: #endif
    971:   }
    972:   PetscHSetIGetSize(ht, &nelems);
    973:   PetscMalloc1(nelems, &elems);
    974:   PetscHSetIGetElems(ht, &off, elems);
    975:   PetscHSetIDestroy(&ht);
    976:   PetscSortInt(nelems, elems);
    977:   ISCreateGeneral(PETSC_COMM_SELF, nelems, elems, PETSC_OWN_POINTER, closureIS);
    978:   return(0);
    979: }

    981: /*@
    982:   DMPlexPartitionLabelClosure - Add the closure of all points to the partition label

    984:   Input Parameters:
    985: + dm     - The DM
    986: - label  - DMLabel assinging ranks to remote roots

    988:   Level: developer

    990: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
    991: @*/
    992: PetscErrorCode DMPlexPartitionLabelClosure(DM dm, DMLabel label)
    993: {
    994:   IS              rankIS,   pointIS, closureIS;
    995:   const PetscInt *ranks,   *points;
    996:   PetscInt        numRanks, numPoints, r;
    997:   PetscErrorCode  ierr;

    1000:   DMLabelGetValueIS(label, &rankIS);
    1001:   ISGetLocalSize(rankIS, &numRanks);
    1002:   ISGetIndices(rankIS, &ranks);
    1003:   for (r = 0; r < numRanks; ++r) {
    1004:     const PetscInt rank = ranks[r];
    1005:     DMLabelGetStratumIS(label, rank, &pointIS);
    1006:     ISGetLocalSize(pointIS, &numPoints);
    1007:     ISGetIndices(pointIS, &points);
    1008:     DMPlexClosurePoints_Private(dm, numPoints, points, &closureIS);
    1009:     ISRestoreIndices(pointIS, &points);
    1010:     ISDestroy(&pointIS);
    1011:     DMLabelSetStratumIS(label, rank, closureIS);
    1012:     ISDestroy(&closureIS);
    1013:   }
    1014:   ISRestoreIndices(rankIS, &ranks);
    1015:   ISDestroy(&rankIS);
    1016:   return(0);
    1017: }

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

    1022:   Input Parameters:
    1023: + dm     - The DM
    1024: - label  - DMLabel assinging ranks to remote roots

    1026:   Level: developer

    1028: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
    1029: @*/
    1030: PetscErrorCode DMPlexPartitionLabelAdjacency(DM dm, DMLabel label)
    1031: {
    1032:   IS              rankIS,   pointIS;
    1033:   const PetscInt *ranks,   *points;
    1034:   PetscInt        numRanks, numPoints, r, p, a, adjSize;
    1035:   PetscInt       *adj = NULL;
    1036:   PetscErrorCode  ierr;

    1039:   DMLabelGetValueIS(label, &rankIS);
    1040:   ISGetLocalSize(rankIS, &numRanks);
    1041:   ISGetIndices(rankIS, &ranks);
    1042:   for (r = 0; r < numRanks; ++r) {
    1043:     const PetscInt rank = ranks[r];

    1045:     DMLabelGetStratumIS(label, rank, &pointIS);
    1046:     ISGetLocalSize(pointIS, &numPoints);
    1047:     ISGetIndices(pointIS, &points);
    1048:     for (p = 0; p < numPoints; ++p) {
    1049:       adjSize = PETSC_DETERMINE;
    1050:       DMPlexGetAdjacency(dm, points[p], &adjSize, &adj);
    1051:       for (a = 0; a < adjSize; ++a) {DMLabelSetValue(label, adj[a], rank);}
    1052:     }
    1053:     ISRestoreIndices(pointIS, &points);
    1054:     ISDestroy(&pointIS);
    1055:   }
    1056:   ISRestoreIndices(rankIS, &ranks);
    1057:   ISDestroy(&rankIS);
    1058:   PetscFree(adj);
    1059:   return(0);
    1060: }

    1062: /*@
    1063:   DMPlexPartitionLabelPropagate - Propagate points in a partition label over the point SF

    1065:   Input Parameters:
    1066: + dm     - The DM
    1067: - label  - DMLabel assinging ranks to remote roots

    1069:   Level: developer

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

    1074: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
    1075: @*/
    1076: PetscErrorCode DMPlexPartitionLabelPropagate(DM dm, DMLabel label)
    1077: {
    1078:   MPI_Comm        comm;
    1079:   PetscMPIInt     rank;
    1080:   PetscSF         sfPoint;
    1081:   DMLabel         lblRoots, lblLeaves;
    1082:   IS              rankIS, pointIS;
    1083:   const PetscInt *ranks;
    1084:   PetscInt        numRanks, r;
    1085:   PetscErrorCode  ierr;

    1088:   PetscObjectGetComm((PetscObject) dm, &comm);
    1089:   MPI_Comm_rank(comm, &rank);
    1090:   DMGetPointSF(dm, &sfPoint);
    1091:   /* Pull point contributions from remote leaves into local roots */
    1092:   DMLabelGather(label, sfPoint, &lblLeaves);
    1093:   DMLabelGetValueIS(lblLeaves, &rankIS);
    1094:   ISGetLocalSize(rankIS, &numRanks);
    1095:   ISGetIndices(rankIS, &ranks);
    1096:   for (r = 0; r < numRanks; ++r) {
    1097:     const PetscInt remoteRank = ranks[r];
    1098:     if (remoteRank == rank) continue;
    1099:     DMLabelGetStratumIS(lblLeaves, remoteRank, &pointIS);
    1100:     DMLabelInsertIS(label, pointIS, remoteRank);
    1101:     ISDestroy(&pointIS);
    1102:   }
    1103:   ISRestoreIndices(rankIS, &ranks);
    1104:   ISDestroy(&rankIS);
    1105:   DMLabelDestroy(&lblLeaves);
    1106:   /* Push point contributions from roots into remote leaves */
    1107:   DMLabelDistribute(label, sfPoint, &lblRoots);
    1108:   DMLabelGetValueIS(lblRoots, &rankIS);
    1109:   ISGetLocalSize(rankIS, &numRanks);
    1110:   ISGetIndices(rankIS, &ranks);
    1111:   for (r = 0; r < numRanks; ++r) {
    1112:     const PetscInt remoteRank = ranks[r];
    1113:     if (remoteRank == rank) continue;
    1114:     DMLabelGetStratumIS(lblRoots, remoteRank, &pointIS);
    1115:     DMLabelInsertIS(label, pointIS, remoteRank);
    1116:     ISDestroy(&pointIS);
    1117:   }
    1118:   ISRestoreIndices(rankIS, &ranks);
    1119:   ISDestroy(&rankIS);
    1120:   DMLabelDestroy(&lblRoots);
    1121:   return(0);
    1122: }

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

    1127:   Input Parameters:
    1128: + dm        - The DM
    1129: . rootLabel - DMLabel assinging ranks to local roots
    1130: - processSF - A star forest mapping into the local index on each remote rank

    1132:   Output Parameter:
    1133: . leafLabel - DMLabel assinging ranks to remote roots

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

    1138:   Level: developer

    1140: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
    1141: @*/
    1142: PetscErrorCode DMPlexPartitionLabelInvert(DM dm, DMLabel rootLabel, PetscSF processSF, DMLabel leafLabel)
    1143: {
    1144:   MPI_Comm           comm;
    1145:   PetscMPIInt        rank, size, r;
    1146:   PetscInt           p, n, numNeighbors, numPoints, dof, off, rootSize, l, nleaves, leafSize;
    1147:   PetscSF            sfPoint;
    1148:   PetscSection       rootSection;
    1149:   PetscSFNode       *rootPoints, *leafPoints;
    1150:   const PetscSFNode *remote;
    1151:   const PetscInt    *local, *neighbors;
    1152:   IS                 valueIS;
    1153:   PetscBool          mpiOverflow = PETSC_FALSE;
    1154:   PetscErrorCode     ierr;

    1157:   PetscLogEventBegin(DMPLEX_PartLabelInvert,dm,0,0,0);
    1158:   PetscObjectGetComm((PetscObject) dm, &comm);
    1159:   MPI_Comm_rank(comm, &rank);
    1160:   MPI_Comm_size(comm, &size);
    1161:   DMGetPointSF(dm, &sfPoint);

    1163:   /* Convert to (point, rank) and use actual owners */
    1164:   PetscSectionCreate(comm, &rootSection);
    1165:   PetscSectionSetChart(rootSection, 0, size);
    1166:   DMLabelGetValueIS(rootLabel, &valueIS);
    1167:   ISGetLocalSize(valueIS, &numNeighbors);
    1168:   ISGetIndices(valueIS, &neighbors);
    1169:   for (n = 0; n < numNeighbors; ++n) {
    1170:     DMLabelGetStratumSize(rootLabel, neighbors[n], &numPoints);
    1171:     PetscSectionAddDof(rootSection, neighbors[n], numPoints);
    1172:   }
    1173:   PetscSectionSetUp(rootSection);
    1174:   PetscSectionGetStorageSize(rootSection, &rootSize);
    1175:   PetscMalloc1(rootSize, &rootPoints);
    1176:   PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote);
    1177:   for (n = 0; n < numNeighbors; ++n) {
    1178:     IS              pointIS;
    1179:     const PetscInt *points;

    1181:     PetscSectionGetOffset(rootSection, neighbors[n], &off);
    1182:     DMLabelGetStratumIS(rootLabel, neighbors[n], &pointIS);
    1183:     ISGetLocalSize(pointIS, &numPoints);
    1184:     ISGetIndices(pointIS, &points);
    1185:     for (p = 0; p < numPoints; ++p) {
    1186:       if (local) {PetscFindInt(points[p], nleaves, local, &l);}
    1187:       else       {l = -1;}
    1188:       if (l >= 0) {rootPoints[off+p] = remote[l];}
    1189:       else        {rootPoints[off+p].index = points[p]; rootPoints[off+p].rank = rank;}
    1190:     }
    1191:     ISRestoreIndices(pointIS, &points);
    1192:     ISDestroy(&pointIS);
    1193:   }

    1195:   /* Try to communicate overlap using All-to-All */
    1196:   if (!processSF) {
    1197:     PetscInt64  counter = 0;
    1198:     PetscBool   locOverflow = PETSC_FALSE;
    1199:     PetscMPIInt *scounts, *sdispls, *rcounts, *rdispls;

    1201:     PetscCalloc4(size, &scounts, size, &sdispls, size, &rcounts, size, &rdispls);
    1202:     for (n = 0; n < numNeighbors; ++n) {
    1203:       PetscSectionGetDof(rootSection, neighbors[n], &dof);
    1204:       PetscSectionGetOffset(rootSection, neighbors[n], &off);
    1205: #if defined(PETSC_USE_64BIT_INDICES)
    1206:       if (dof > PETSC_MPI_INT_MAX) {locOverflow = PETSC_TRUE; break;}
    1207:       if (off > PETSC_MPI_INT_MAX) {locOverflow = PETSC_TRUE; break;}
    1208: #endif
    1209:       scounts[neighbors[n]] = (PetscMPIInt) dof;
    1210:       sdispls[neighbors[n]] = (PetscMPIInt) off;
    1211:     }
    1212:     MPI_Alltoall(scounts, 1, MPI_INT, rcounts, 1, MPI_INT, comm);
    1213:     for (r = 0; r < size; ++r) { rdispls[r] = (int)counter; counter += rcounts[r]; }
    1214:     if (counter > PETSC_MPI_INT_MAX) locOverflow = PETSC_TRUE;
    1215:     MPI_Allreduce(&locOverflow, &mpiOverflow, 1, MPIU_BOOL, MPI_LOR, comm);
    1216:     if (!mpiOverflow) {
    1217:       PetscInfo(dm,"Using Alltoallv for mesh distribution\n");
    1218:       leafSize = (PetscInt) counter;
    1219:       PetscMalloc1(leafSize, &leafPoints);
    1220:       MPI_Alltoallv(rootPoints, scounts, sdispls, MPIU_2INT, leafPoints, rcounts, rdispls, MPIU_2INT, comm);
    1221:     }
    1222:     PetscFree4(scounts, sdispls, rcounts, rdispls);
    1223:   }

    1225:   /* Communicate overlap using process star forest */
    1226:   if (processSF || mpiOverflow) {
    1227:     PetscSF      procSF;
    1228:     PetscSection leafSection;

    1230:     if (processSF) {
    1231:       PetscInfo(dm,"Using processSF for mesh distribution\n");
    1232:       PetscObjectReference((PetscObject)processSF);
    1233:       procSF = processSF;
    1234:     } else {
    1235:       PetscInfo(dm,"Using processSF for mesh distribution (MPI overflow)\n");
    1236:       PetscSFCreate(comm,&procSF);
    1237:       PetscSFSetGraphWithPattern(procSF,NULL,PETSCSF_PATTERN_ALLTOALL);
    1238:     }

    1240:     PetscSectionCreate(PetscObjectComm((PetscObject)dm), &leafSection);
    1241:     DMPlexDistributeData(dm, procSF, rootSection, MPIU_2INT, rootPoints, leafSection, (void**) &leafPoints);
    1242:     PetscSectionGetStorageSize(leafSection, &leafSize);
    1243:     PetscSectionDestroy(&leafSection);
    1244:     PetscSFDestroy(&procSF);
    1245:   }

    1247:   for (p = 0; p < leafSize; p++) {
    1248:     DMLabelSetValue(leafLabel, leafPoints[p].index, leafPoints[p].rank);
    1249:   }

    1251:   ISRestoreIndices(valueIS, &neighbors);
    1252:   ISDestroy(&valueIS);
    1253:   PetscSectionDestroy(&rootSection);
    1254:   PetscFree(rootPoints);
    1255:   PetscFree(leafPoints);
    1256:   PetscLogEventEnd(DMPLEX_PartLabelInvert,dm,0,0,0);
    1257:   return(0);
    1258: }

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

    1263:   Input Parameters:
    1264: + dm    - The DM
    1265: - label - DMLabel assinging ranks to remote roots

    1267:   Output Parameter:
    1268: . sf    - The star forest communication context encapsulating the defined mapping

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

    1272:   Level: developer

    1274: .seealso: DMPlexDistribute(), DMPlexCreateOverlap()
    1275: @*/
    1276: PetscErrorCode DMPlexPartitionLabelCreateSF(DM dm, DMLabel label, PetscSF *sf)
    1277: {
    1278:   PetscMPIInt     rank;
    1279:   PetscInt        n, numRemote, p, numPoints, pStart, pEnd, idx = 0, nNeighbors;
    1280:   PetscSFNode    *remotePoints;
    1281:   IS              remoteRootIS, neighborsIS;
    1282:   const PetscInt *remoteRoots, *neighbors;

    1286:   PetscLogEventBegin(DMPLEX_PartLabelCreateSF,dm,0,0,0);
    1287:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);

    1289:   DMLabelGetValueIS(label, &neighborsIS);
    1290: #if 0
    1291:   {
    1292:     IS is;
    1293:     ISDuplicate(neighborsIS, &is);
    1294:     ISSort(is);
    1295:     ISDestroy(&neighborsIS);
    1296:     neighborsIS = is;
    1297:   }
    1298: #endif
    1299:   ISGetLocalSize(neighborsIS, &nNeighbors);
    1300:   ISGetIndices(neighborsIS, &neighbors);
    1301:   for (numRemote = 0, n = 0; n < nNeighbors; ++n) {
    1302:     DMLabelGetStratumSize(label, neighbors[n], &numPoints);
    1303:     numRemote += numPoints;
    1304:   }
    1305:   PetscMalloc1(numRemote, &remotePoints);
    1306:   /* Put owned points first */
    1307:   DMLabelGetStratumSize(label, rank, &numPoints);
    1308:   if (numPoints > 0) {
    1309:     DMLabelGetStratumIS(label, rank, &remoteRootIS);
    1310:     ISGetIndices(remoteRootIS, &remoteRoots);
    1311:     for (p = 0; p < numPoints; p++) {
    1312:       remotePoints[idx].index = remoteRoots[p];
    1313:       remotePoints[idx].rank = rank;
    1314:       idx++;
    1315:     }
    1316:     ISRestoreIndices(remoteRootIS, &remoteRoots);
    1317:     ISDestroy(&remoteRootIS);
    1318:   }
    1319:   /* Now add remote points */
    1320:   for (n = 0; n < nNeighbors; ++n) {
    1321:     const PetscInt nn = neighbors[n];

    1323:     DMLabelGetStratumSize(label, nn, &numPoints);
    1324:     if (nn == rank || numPoints <= 0) continue;
    1325:     DMLabelGetStratumIS(label, nn, &remoteRootIS);
    1326:     ISGetIndices(remoteRootIS, &remoteRoots);
    1327:     for (p = 0; p < numPoints; p++) {
    1328:       remotePoints[idx].index = remoteRoots[p];
    1329:       remotePoints[idx].rank = nn;
    1330:       idx++;
    1331:     }
    1332:     ISRestoreIndices(remoteRootIS, &remoteRoots);
    1333:     ISDestroy(&remoteRootIS);
    1334:   }
    1335:   PetscSFCreate(PetscObjectComm((PetscObject) dm), sf);
    1336:   DMPlexGetChart(dm, &pStart, &pEnd);
    1337:   PetscSFSetGraph(*sf, pEnd-pStart, numRemote, NULL, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
    1338:   PetscSFSetUp(*sf);
    1339:   ISDestroy(&neighborsIS);
    1340:   PetscLogEventEnd(DMPLEX_PartLabelCreateSF,dm,0,0,0);
    1341:   return(0);
    1342: }

    1344: #if defined(PETSC_HAVE_PARMETIS)
    1345: #include <parmetis.h>
    1346: #endif

    1348: /* The two functions below are used by DMPlexRebalanceSharedPoints which errors
    1349:  * when PETSc is built without ParMETIS. To avoid -Wunused-function, we take
    1350:  * them out in that case. */
    1351: #if defined(PETSC_HAVE_PARMETIS)
    1352: /*@C

    1354:   DMPlexRewriteSF - Rewrites the ownership of the SF of a DM (in place).

    1356:   Input parameters:
    1357: + dm                - The DMPlex object.
    1358: . n                 - The number of points.
    1359: . pointsToRewrite   - The points in the SF whose ownership will change.
    1360: . targetOwners      - New owner for each element in pointsToRewrite.
    1361: - degrees           - Degrees of the points in the SF as obtained by PetscSFComputeDegreeBegin/PetscSFComputeDegreeEnd.

    1363:   Level: developer

    1365: @*/
    1366: static PetscErrorCode DMPlexRewriteSF(DM dm, PetscInt n, PetscInt *pointsToRewrite, PetscInt *targetOwners, const PetscInt *degrees)
    1367: {
    1368:   PetscInt      ierr, pStart, pEnd, i, j, counter, leafCounter, sumDegrees, nroots, nleafs;
    1369:   PetscInt     *cumSumDegrees, *newOwners, *newNumbers, *rankOnLeafs, *locationsOfLeafs, *remoteLocalPointOfLeafs, *points, *leafsNew;
    1370:   PetscSFNode  *leafLocationsNew;
    1371:   const         PetscSFNode *iremote;
    1372:   const         PetscInt *ilocal;
    1373:   PetscBool    *isLeaf;
    1374:   PetscSF       sf;
    1375:   MPI_Comm      comm;
    1376:   PetscMPIInt   rank, size;

    1379:   PetscObjectGetComm((PetscObject) dm, &comm);
    1380:   MPI_Comm_rank(comm, &rank);
    1381:   MPI_Comm_size(comm, &size);
    1382:   DMPlexGetChart(dm, &pStart, &pEnd);

    1384:   DMGetPointSF(dm, &sf);
    1385:   PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
    1386:   PetscMalloc1(pEnd-pStart, &isLeaf);
    1387:   for (i=0; i<pEnd-pStart; i++) {
    1388:     isLeaf[i] = PETSC_FALSE;
    1389:   }
    1390:   for (i=0; i<nleafs; i++) {
    1391:     isLeaf[ilocal[i]-pStart] = PETSC_TRUE;
    1392:   }

    1394:   PetscMalloc1(pEnd-pStart+1, &cumSumDegrees);
    1395:   cumSumDegrees[0] = 0;
    1396:   for (i=1; i<=pEnd-pStart; i++) {
    1397:     cumSumDegrees[i] = cumSumDegrees[i-1] + degrees[i-1];
    1398:   }
    1399:   sumDegrees = cumSumDegrees[pEnd-pStart];
    1400:   /* get the location of my leafs (we have sumDegrees many leafs pointing at our roots) */

    1402:   PetscMalloc1(sumDegrees, &locationsOfLeafs);
    1403:   PetscMalloc1(pEnd-pStart, &rankOnLeafs);
    1404:   for (i=0; i<pEnd-pStart; i++) {
    1405:     rankOnLeafs[i] = rank;
    1406:   }
    1407:   PetscSFGatherBegin(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
    1408:   PetscSFGatherEnd(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
    1409:   PetscFree(rankOnLeafs);

    1411:   /* get the remote local points of my leaves */
    1412:   PetscMalloc1(sumDegrees, &remoteLocalPointOfLeafs);
    1413:   PetscMalloc1(pEnd-pStart, &points);
    1414:   for (i=0; i<pEnd-pStart; i++) {
    1415:     points[i] = pStart+i;
    1416:   }
    1417:   PetscSFGatherBegin(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
    1418:   PetscSFGatherEnd(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
    1419:   PetscFree(points);
    1420:   /* Figure out the new owners of the vertices that are up for grabs and their numbers on the new owners */
    1421:   PetscMalloc1(pEnd-pStart, &newOwners);
    1422:   PetscMalloc1(pEnd-pStart, &newNumbers);
    1423:   for (i=0; i<pEnd-pStart; i++) {
    1424:     newOwners[i] = -1;
    1425:     newNumbers[i] = -1;
    1426:   }
    1427:   {
    1428:     PetscInt oldNumber, newNumber, oldOwner, newOwner;
    1429:     for (i=0; i<n; i++) {
    1430:       oldNumber = pointsToRewrite[i];
    1431:       newNumber = -1;
    1432:       oldOwner = rank;
    1433:       newOwner = targetOwners[i];
    1434:       if (oldOwner == newOwner) {
    1435:         newNumber = oldNumber;
    1436:       } else {
    1437:         for (j=0; j<degrees[oldNumber]; j++) {
    1438:           if (locationsOfLeafs[cumSumDegrees[oldNumber]+j] == newOwner) {
    1439:             newNumber = remoteLocalPointOfLeafs[cumSumDegrees[oldNumber]+j];
    1440:             break;
    1441:           }
    1442:         }
    1443:       }
    1444:       if (newNumber == -1) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Couldn't find the new owner of vertex.");

    1446:       newOwners[oldNumber] = newOwner;
    1447:       newNumbers[oldNumber] = newNumber;
    1448:     }
    1449:   }
    1450:   PetscFree(cumSumDegrees);
    1451:   PetscFree(locationsOfLeafs);
    1452:   PetscFree(remoteLocalPointOfLeafs);

    1454:   PetscSFBcastBegin(sf, MPIU_INT, newOwners, newOwners,MPI_REPLACE);
    1455:   PetscSFBcastEnd(sf, MPIU_INT, newOwners, newOwners,MPI_REPLACE);
    1456:   PetscSFBcastBegin(sf, MPIU_INT, newNumbers, newNumbers,MPI_REPLACE);
    1457:   PetscSFBcastEnd(sf, MPIU_INT, newNumbers, newNumbers,MPI_REPLACE);

    1459:   /* Now count how many leafs we have on each processor. */
    1460:   leafCounter=0;
    1461:   for (i=0; i<pEnd-pStart; i++) {
    1462:     if (newOwners[i] >= 0) {
    1463:       if (newOwners[i] != rank) {
    1464:         leafCounter++;
    1465:       }
    1466:     } else {
    1467:       if (isLeaf[i]) {
    1468:         leafCounter++;
    1469:       }
    1470:     }
    1471:   }

    1473:   /* Now set up the new sf by creating the leaf arrays */
    1474:   PetscMalloc1(leafCounter, &leafsNew);
    1475:   PetscMalloc1(leafCounter, &leafLocationsNew);

    1477:   leafCounter = 0;
    1478:   counter = 0;
    1479:   for (i=0; i<pEnd-pStart; i++) {
    1480:     if (newOwners[i] >= 0) {
    1481:       if (newOwners[i] != rank) {
    1482:         leafsNew[leafCounter] = i;
    1483:         leafLocationsNew[leafCounter].rank = newOwners[i];
    1484:         leafLocationsNew[leafCounter].index = newNumbers[i];
    1485:         leafCounter++;
    1486:       }
    1487:     } else {
    1488:       if (isLeaf[i]) {
    1489:         leafsNew[leafCounter] = i;
    1490:         leafLocationsNew[leafCounter].rank = iremote[counter].rank;
    1491:         leafLocationsNew[leafCounter].index = iremote[counter].index;
    1492:         leafCounter++;
    1493:       }
    1494:     }
    1495:     if (isLeaf[i]) {
    1496:       counter++;
    1497:     }
    1498:   }

    1500:   PetscSFSetGraph(sf, nroots, leafCounter, leafsNew, PETSC_OWN_POINTER, leafLocationsNew, PETSC_OWN_POINTER);
    1501:   PetscFree(newOwners);
    1502:   PetscFree(newNumbers);
    1503:   PetscFree(isLeaf);
    1504:   return(0);
    1505: }

    1507: static PetscErrorCode DMPlexViewDistribution(MPI_Comm comm, PetscInt n, PetscInt skip, PetscInt *vtxwgt, PetscInt *part, PetscViewer viewer)
    1508: {
    1509:   PetscInt *distribution, min, max, sum, i, ierr;
    1510:   PetscMPIInt rank, size;
    1512:   MPI_Comm_size(comm, &size);
    1513:   MPI_Comm_rank(comm, &rank);
    1514:   PetscCalloc1(size, &distribution);
    1515:   for (i=0; i<n; i++) {
    1516:     if (part) distribution[part[i]] += vtxwgt[skip*i];
    1517:     else distribution[rank] += vtxwgt[skip*i];
    1518:   }
    1519:   MPI_Allreduce(MPI_IN_PLACE, distribution, size, MPIU_INT, MPI_SUM, comm);
    1520:   min = distribution[0];
    1521:   max = distribution[0];
    1522:   sum = distribution[0];
    1523:   for (i=1; i<size; i++) {
    1524:     if (distribution[i]<min) min=distribution[i];
    1525:     if (distribution[i]>max) max=distribution[i];
    1526:     sum += distribution[i];
    1527:   }
    1528:   PetscViewerASCIIPrintf(viewer, "Min: %D, Avg: %D, Max: %D, Balance: %f\n", min, sum/size, max, (max*1.*size)/sum);
    1529:   PetscFree(distribution);
    1530:   return(0);
    1531: }

    1533: #endif

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

    1538:   Input parameters:
    1539: + dm               - The DMPlex object.
    1540: . entityDepth      - depth of the entity to balance (0 -> balance vertices).
    1541: . useInitialGuess  - whether to use the current distribution as initial guess (only used by ParMETIS).
    1542: - parallel         - whether to use ParMETIS and do the partition in parallel or whether to gather the graph onto a single process and use METIS.

    1544:   Output parameters:
    1545: . success          - whether the graph partitioning was successful or not. If not, try useInitialGuess=True and parallel=True.

    1547:   Level: intermediate

    1549: @*/

    1551: PetscErrorCode DMPlexRebalanceSharedPoints(DM dm, PetscInt entityDepth, PetscBool useInitialGuess, PetscBool parallel, PetscBool *success)
    1552: {
    1553: #if defined(PETSC_HAVE_PARMETIS)
    1554:   PetscSF     sf;
    1555:   PetscInt    ierr, i, j, idx, jdx;
    1556:   PetscInt    eBegin, eEnd, nroots, nleafs, pStart, pEnd;
    1557:   const       PetscInt *degrees, *ilocal;
    1558:   const       PetscSFNode *iremote;
    1559:   PetscBool   *toBalance, *isLeaf, *isExclusivelyOwned, *isNonExclusivelyOwned;
    1560:   PetscInt    numExclusivelyOwned, numNonExclusivelyOwned;
    1561:   PetscMPIInt rank, size;
    1562:   PetscInt    *globalNumbersOfLocalOwnedVertices, *leafGlobalNumbers;
    1563:   const       PetscInt *cumSumVertices;
    1564:   PetscInt    offset, counter;
    1565:   PetscInt    lenadjncy;
    1566:   PetscInt    *xadj, *adjncy, *vtxwgt;
    1567:   PetscInt    lenxadj;
    1568:   PetscInt    *adjwgt = NULL;
    1569:   PetscInt    *part, *options;
    1570:   PetscInt    nparts, wgtflag, numflag, ncon, edgecut;
    1571:   real_t      *ubvec;
    1572:   PetscInt    *firstVertices, *renumbering;
    1573:   PetscInt    failed, failedGlobal;
    1574:   MPI_Comm    comm;
    1575:   Mat         A, Apre;
    1576:   const char *prefix = NULL;
    1577:   PetscViewer       viewer;
    1578:   PetscViewerFormat format;
    1579:   PetscLayout layout;

    1582:   if (success) *success = PETSC_FALSE;
    1583:   PetscObjectGetComm((PetscObject) dm, &comm);
    1584:   MPI_Comm_rank(comm, &rank);
    1585:   MPI_Comm_size(comm, &size);
    1586:   if (size==1) return(0);

    1588:   PetscLogEventBegin(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);

    1590:   PetscOptionsGetViewer(comm,((PetscObject)dm)->options, prefix,"-dm_rebalance_partition_view",&viewer,&format,NULL);
    1591:   if (viewer) {
    1592:     PetscViewerPushFormat(viewer,format);
    1593:   }

    1595:   /* Figure out all points in the plex that we are interested in balancing. */
    1596:   DMPlexGetDepthStratum(dm, entityDepth, &eBegin, &eEnd);
    1597:   DMPlexGetChart(dm, &pStart, &pEnd);
    1598:   PetscMalloc1(pEnd-pStart, &toBalance);

    1600:   for (i=0; i<pEnd-pStart; i++) {
    1601:     toBalance[i] = (PetscBool)(i-pStart>=eBegin && i-pStart<eEnd);
    1602:   }

    1604:   /* There are three types of points:
    1605:    * exclusivelyOwned: points that are owned by this process and only seen by this process
    1606:    * nonExclusivelyOwned: points that are owned by this process but seen by at least another process
    1607:    * leaf: a point that is seen by this process but owned by a different process
    1608:    */
    1609:   DMGetPointSF(dm, &sf);
    1610:   PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
    1611:   PetscMalloc1(pEnd-pStart, &isLeaf);
    1612:   PetscMalloc1(pEnd-pStart, &isNonExclusivelyOwned);
    1613:   PetscMalloc1(pEnd-pStart, &isExclusivelyOwned);
    1614:   for (i=0; i<pEnd-pStart; i++) {
    1615:     isNonExclusivelyOwned[i] = PETSC_FALSE;
    1616:     isExclusivelyOwned[i] = PETSC_FALSE;
    1617:     isLeaf[i] = PETSC_FALSE;
    1618:   }

    1620:   /* start by marking all the leafs */
    1621:   for (i=0; i<nleafs; i++) {
    1622:     isLeaf[ilocal[i]-pStart] = PETSC_TRUE;
    1623:   }

    1625:   /* for an owned point, we can figure out whether another processor sees it or
    1626:    * not by calculating its degree */
    1627:   PetscSFComputeDegreeBegin(sf, &degrees);
    1628:   PetscSFComputeDegreeEnd(sf, &degrees);

    1630:   numExclusivelyOwned = 0;
    1631:   numNonExclusivelyOwned = 0;
    1632:   for (i=0; i<pEnd-pStart; i++) {
    1633:     if (toBalance[i]) {
    1634:       if (degrees[i] > 0) {
    1635:         isNonExclusivelyOwned[i] = PETSC_TRUE;
    1636:         numNonExclusivelyOwned += 1;
    1637:       } else {
    1638:         if (!isLeaf[i]) {
    1639:           isExclusivelyOwned[i] = PETSC_TRUE;
    1640:           numExclusivelyOwned += 1;
    1641:         }
    1642:       }
    1643:     }
    1644:   }

    1646:   /* We are going to build a graph with one vertex per core representing the
    1647:    * exclusively owned points and then one vertex per nonExclusively owned
    1648:    * point. */

    1650:   PetscLayoutCreate(comm, &layout);
    1651:   PetscLayoutSetLocalSize(layout, 1 + numNonExclusivelyOwned);
    1652:   PetscLayoutSetUp(layout);
    1653:   PetscLayoutGetRanges(layout, &cumSumVertices);

    1655:   PetscMalloc1(pEnd-pStart, &globalNumbersOfLocalOwnedVertices);
    1656:   for (i=0; i<pEnd-pStart; i++) {globalNumbersOfLocalOwnedVertices[i] = pStart - 1;}
    1657:   offset = cumSumVertices[rank];
    1658:   counter = 0;
    1659:   for (i=0; i<pEnd-pStart; i++) {
    1660:     if (toBalance[i]) {
    1661:       if (degrees[i] > 0) {
    1662:         globalNumbersOfLocalOwnedVertices[i] = counter + 1 + offset;
    1663:         counter++;
    1664:       }
    1665:     }
    1666:   }

    1668:   /* send the global numbers of vertices I own to the leafs so that they know to connect to it */
    1669:   PetscMalloc1(pEnd-pStart, &leafGlobalNumbers);
    1670:   PetscSFBcastBegin(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers,MPI_REPLACE);
    1671:   PetscSFBcastEnd(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers,MPI_REPLACE);

    1673:   /* Now start building the data structure for ParMETIS */

    1675:   MatCreate(comm, &Apre);
    1676:   MatSetType(Apre, MATPREALLOCATOR);
    1677:   MatSetSizes(Apre, 1+numNonExclusivelyOwned, 1+numNonExclusivelyOwned, cumSumVertices[size], cumSumVertices[size]);
    1678:   MatSetUp(Apre);

    1680:   for (i=0; i<pEnd-pStart; i++) {
    1681:     if (toBalance[i]) {
    1682:       idx = cumSumVertices[rank];
    1683:       if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
    1684:       else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
    1685:       else continue;
    1686:       MatSetValue(Apre, idx, jdx, 1, INSERT_VALUES);
    1687:       MatSetValue(Apre, jdx, idx, 1, INSERT_VALUES);
    1688:     }
    1689:   }

    1691:   MatAssemblyBegin(Apre, MAT_FINAL_ASSEMBLY);
    1692:   MatAssemblyEnd(Apre, MAT_FINAL_ASSEMBLY);

    1694:   MatCreate(comm, &A);
    1695:   MatSetType(A, MATMPIAIJ);
    1696:   MatSetSizes(A, 1+numNonExclusivelyOwned, 1+numNonExclusivelyOwned, cumSumVertices[size], cumSumVertices[size]);
    1697:   MatPreallocatorPreallocate(Apre, PETSC_FALSE, A);
    1698:   MatDestroy(&Apre);

    1700:   for (i=0; i<pEnd-pStart; i++) {
    1701:     if (toBalance[i]) {
    1702:       idx = cumSumVertices[rank];
    1703:       if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
    1704:       else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
    1705:       else continue;
    1706:       MatSetValue(A, idx, jdx, 1, INSERT_VALUES);
    1707:       MatSetValue(A, jdx, idx, 1, INSERT_VALUES);
    1708:     }
    1709:   }

    1711:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    1712:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
    1713:   PetscFree(leafGlobalNumbers);
    1714:   PetscFree(globalNumbersOfLocalOwnedVertices);

    1716:   nparts = size;
    1717:   wgtflag = 2;
    1718:   numflag = 0;
    1719:   ncon = 2;
    1720:   real_t *tpwgts;
    1721:   PetscMalloc1(ncon * nparts, &tpwgts);
    1722:   for (i=0; i<ncon*nparts; i++) {
    1723:     tpwgts[i] = 1./(nparts);
    1724:   }

    1726:   PetscMalloc1(ncon, &ubvec);
    1727:   ubvec[0] = 1.01;
    1728:   ubvec[1] = 1.01;
    1729:   lenadjncy = 0;
    1730:   for (i=0; i<1+numNonExclusivelyOwned; i++) {
    1731:     PetscInt temp=0;
    1732:     MatGetRow(A, cumSumVertices[rank] + i, &temp, NULL, NULL);
    1733:     lenadjncy += temp;
    1734:     MatRestoreRow(A, cumSumVertices[rank] + i, &temp, NULL, NULL);
    1735:   }
    1736:   PetscMalloc1(lenadjncy, &adjncy);
    1737:   lenxadj = 2 + numNonExclusivelyOwned;
    1738:   PetscMalloc1(lenxadj, &xadj);
    1739:   xadj[0] = 0;
    1740:   counter = 0;
    1741:   for (i=0; i<1+numNonExclusivelyOwned; i++) {
    1742:     PetscInt        temp=0;
    1743:     const PetscInt *cols;
    1744:     MatGetRow(A, cumSumVertices[rank] + i, &temp, &cols, NULL);
    1745:     PetscArraycpy(&adjncy[counter], cols, temp);
    1746:     counter += temp;
    1747:     xadj[i+1] = counter;
    1748:     MatRestoreRow(A, cumSumVertices[rank] + i, &temp, &cols, NULL);
    1749:   }

    1751:   PetscMalloc1(cumSumVertices[rank+1]-cumSumVertices[rank], &part);
    1752:   PetscMalloc1(ncon*(1 + numNonExclusivelyOwned), &vtxwgt);
    1753:   vtxwgt[0] = numExclusivelyOwned;
    1754:   if (ncon>1) vtxwgt[1] = 1;
    1755:   for (i=0; i<numNonExclusivelyOwned; i++) {
    1756:     vtxwgt[ncon*(i+1)] = 1;
    1757:     if (ncon>1) vtxwgt[ncon*(i+1)+1] = 0;
    1758:   }

    1760:   if (viewer) {
    1761:     PetscViewerASCIIPrintf(viewer, "Attempt rebalancing of shared points of depth %D on interface of mesh distribution.\n", entityDepth);
    1762:     PetscViewerASCIIPrintf(viewer, "Size of generated auxiliary graph: %D\n", cumSumVertices[size]);
    1763:   }
    1764:   if (parallel) {
    1765:     PetscMalloc1(4, &options);
    1766:     options[0] = 1;
    1767:     options[1] = 0; /* Verbosity */
    1768:     options[2] = 0; /* Seed */
    1769:     options[3] = PARMETIS_PSR_COUPLED; /* Seed */
    1770:     if (viewer) { PetscViewerASCIIPrintf(viewer, "Using ParMETIS to partition graph.\n"); }
    1771:     if (useInitialGuess) {
    1772:       if (viewer) { PetscViewerASCIIPrintf(viewer, "Using current distribution of points as initial guess.\n"); }
    1773:       PetscStackPush("ParMETIS_V3_RefineKway");
    1774:       ParMETIS_V3_RefineKway((PetscInt*)cumSumVertices, xadj, adjncy, vtxwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
    1775:       if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ParMETIS_V3_RefineKway()");
    1776:       PetscStackPop;
    1777:     } else {
    1778:       PetscStackPush("ParMETIS_V3_PartKway");
    1779:       ParMETIS_V3_PartKway((PetscInt*)cumSumVertices, xadj, adjncy, vtxwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
    1780:       PetscStackPop;
    1781:       if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ParMETIS_V3_PartKway()");
    1782:     }
    1783:     PetscFree(options);
    1784:   } else {
    1785:     if (viewer) { PetscViewerASCIIPrintf(viewer, "Using METIS to partition graph.\n"); }
    1786:     Mat As;
    1787:     PetscInt numRows;
    1788:     PetscInt *partGlobal;
    1789:     MatCreateRedundantMatrix(A, size, MPI_COMM_NULL, MAT_INITIAL_MATRIX, &As);

    1791:     PetscInt *numExclusivelyOwnedAll;
    1792:     PetscMalloc1(size, &numExclusivelyOwnedAll);
    1793:     numExclusivelyOwnedAll[rank] = numExclusivelyOwned;
    1794:     MPI_Allgather(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,numExclusivelyOwnedAll,1,MPIU_INT,comm);

    1796:     MatGetSize(As, &numRows, NULL);
    1797:     PetscMalloc1(numRows, &partGlobal);
    1798:     if (!rank) {
    1799:       PetscInt *adjncy_g, *xadj_g, *vtxwgt_g;
    1800:       lenadjncy = 0;

    1802:       for (i=0; i<numRows; i++) {
    1803:         PetscInt temp=0;
    1804:         MatGetRow(As, i, &temp, NULL, NULL);
    1805:         lenadjncy += temp;
    1806:         MatRestoreRow(As, i, &temp, NULL, NULL);
    1807:       }
    1808:       PetscMalloc1(lenadjncy, &adjncy_g);
    1809:       lenxadj = 1 + numRows;
    1810:       PetscMalloc1(lenxadj, &xadj_g);
    1811:       xadj_g[0] = 0;
    1812:       counter = 0;
    1813:       for (i=0; i<numRows; i++) {
    1814:         PetscInt        temp=0;
    1815:         const PetscInt *cols;
    1816:         MatGetRow(As, i, &temp, &cols, NULL);
    1817:         PetscArraycpy(&adjncy_g[counter], cols, temp);
    1818:         counter += temp;
    1819:         xadj_g[i+1] = counter;
    1820:         MatRestoreRow(As, i, &temp, &cols, NULL);
    1821:       }
    1822:       PetscMalloc1(2*numRows, &vtxwgt_g);
    1823:       for (i=0; i<size; i++){
    1824:         vtxwgt_g[ncon*cumSumVertices[i]] = numExclusivelyOwnedAll[i];
    1825:         if (ncon>1) vtxwgt_g[ncon*cumSumVertices[i]+1] = 1;
    1826:         for (j=cumSumVertices[i]+1; j<cumSumVertices[i+1]; j++) {
    1827:           vtxwgt_g[ncon*j] = 1;
    1828:           if (ncon>1) vtxwgt_g[2*j+1] = 0;
    1829:         }
    1830:       }
    1831:       PetscMalloc1(64, &options);
    1832:       METIS_SetDefaultOptions(options); /* initialize all defaults */
    1833:       if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_SetDefaultOptions()");
    1834:       options[METIS_OPTION_CONTIG] = 1;
    1835:       PetscStackPush("METIS_PartGraphKway");
    1836:       METIS_PartGraphKway(&numRows, &ncon, xadj_g, adjncy_g, vtxwgt_g, NULL, NULL, &nparts, tpwgts, ubvec, options, &edgecut, partGlobal);
    1837:       PetscStackPop;
    1838:       if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphKway()");
    1839:       PetscFree(options);
    1840:       PetscFree(xadj_g);
    1841:       PetscFree(adjncy_g);
    1842:       PetscFree(vtxwgt_g);
    1843:     }
    1844:     PetscFree(numExclusivelyOwnedAll);

    1846:     /* Now scatter the parts array. */
    1847:     {
    1848:       PetscMPIInt *counts, *mpiCumSumVertices;
    1849:       PetscMalloc1(size, &counts);
    1850:       PetscMalloc1(size+1, &mpiCumSumVertices);
    1851:       for (i=0; i<size; i++) {
    1852:         PetscMPIIntCast(cumSumVertices[i+1] - cumSumVertices[i], &(counts[i]));
    1853:       }
    1854:       for (i=0; i<=size; i++) {
    1855:         PetscMPIIntCast(cumSumVertices[i], &(mpiCumSumVertices[i]));
    1856:       }
    1857:       MPI_Scatterv(partGlobal, counts, mpiCumSumVertices, MPIU_INT, part, counts[rank], MPIU_INT, 0, comm);
    1858:       PetscFree(counts);
    1859:       PetscFree(mpiCumSumVertices);
    1860:     }

    1862:     PetscFree(partGlobal);
    1863:     MatDestroy(&As);
    1864:   }

    1866:   MatDestroy(&A);
    1867:   PetscFree(ubvec);
    1868:   PetscFree(tpwgts);

    1870:   /* Now rename the result so that the vertex resembling the exclusively owned points stays on the same rank */

    1872:   PetscMalloc1(size, &firstVertices);
    1873:   PetscMalloc1(size, &renumbering);
    1874:   firstVertices[rank] = part[0];
    1875:   MPI_Allgather(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,firstVertices,1,MPIU_INT,comm);
    1876:   for (i=0; i<size; i++) {
    1877:     renumbering[firstVertices[i]] = i;
    1878:   }
    1879:   for (i=0; i<cumSumVertices[rank+1]-cumSumVertices[rank]; i++) {
    1880:     part[i] = renumbering[part[i]];
    1881:   }
    1882:   /* Check if the renumbering worked (this can fail when ParMETIS gives fewer partitions than there are processes) */
    1883:   failed = (PetscInt)(part[0] != rank);
    1884:   MPI_Allreduce(&failed, &failedGlobal, 1, MPIU_INT, MPI_SUM, comm);

    1886:   PetscFree(firstVertices);
    1887:   PetscFree(renumbering);

    1889:   if (failedGlobal > 0) {
    1890:     PetscLayoutDestroy(&layout);
    1891:     PetscFree(xadj);
    1892:     PetscFree(adjncy);
    1893:     PetscFree(vtxwgt);
    1894:     PetscFree(toBalance);
    1895:     PetscFree(isLeaf);
    1896:     PetscFree(isNonExclusivelyOwned);
    1897:     PetscFree(isExclusivelyOwned);
    1898:     PetscFree(part);
    1899:     if (viewer) {
    1900:       PetscViewerPopFormat(viewer);
    1901:       PetscViewerDestroy(&viewer);
    1902:     }
    1903:     PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
    1904:     return(0);
    1905:   }

    1907:   /*Let's check how well we did distributing points*/
    1908:   if (viewer) {
    1909:     PetscViewerASCIIPrintf(viewer, "Comparing number of owned entities of depth %D on each process before rebalancing, after rebalancing, and after consistency checks.\n", entityDepth);
    1910:     PetscViewerASCIIPrintf(viewer, "Initial.     ");
    1911:     DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, NULL, viewer);
    1912:     PetscViewerASCIIPrintf(viewer, "Rebalanced.  ");
    1913:     DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, part, viewer);
    1914:   }

    1916:   /* Now check that every vertex is owned by a process that it is actually connected to. */
    1917:   for (i=1; i<=numNonExclusivelyOwned; i++) {
    1918:     PetscInt loc = 0;
    1919:     PetscFindInt(cumSumVertices[part[i]], xadj[i+1]-xadj[i], &adjncy[xadj[i]], &loc);
    1920:     /* If not, then just set the owner to the original owner (hopefully a rare event, it means that a vertex has been isolated) */
    1921:     if (loc<0) {
    1922:       part[i] = rank;
    1923:     }
    1924:   }

    1926:   /* Let's see how significant the influences of the previous fixing up step was.*/
    1927:   if (viewer) {
    1928:     PetscViewerASCIIPrintf(viewer, "After.       ");
    1929:     DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, part, viewer);
    1930:   }

    1932:   PetscLayoutDestroy(&layout);
    1933:   PetscFree(xadj);
    1934:   PetscFree(adjncy);
    1935:   PetscFree(vtxwgt);

    1937:   /* Almost done, now rewrite the SF to reflect the new ownership. */
    1938:   {
    1939:     PetscInt *pointsToRewrite;
    1940:     PetscMalloc1(numNonExclusivelyOwned, &pointsToRewrite);
    1941:     counter = 0;
    1942:     for (i=0; i<pEnd-pStart; i++) {
    1943:       if (toBalance[i]) {
    1944:         if (isNonExclusivelyOwned[i]) {
    1945:           pointsToRewrite[counter] = i + pStart;
    1946:           counter++;
    1947:         }
    1948:       }
    1949:     }
    1950:     DMPlexRewriteSF(dm, numNonExclusivelyOwned, pointsToRewrite, part+1, degrees);
    1951:     PetscFree(pointsToRewrite);
    1952:   }

    1954:   PetscFree(toBalance);
    1955:   PetscFree(isLeaf);
    1956:   PetscFree(isNonExclusivelyOwned);
    1957:   PetscFree(isExclusivelyOwned);
    1958:   PetscFree(part);
    1959:   if (viewer) {
    1960:     PetscViewerPopFormat(viewer);
    1961:     PetscViewerDestroy(&viewer);
    1962:   }
    1963:   if (success) *success = PETSC_TRUE;
    1964:   PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
    1965:   return(0);
    1966: #else
    1967:   SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
    1968: #endif
    1969: }