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:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
 24:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
 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 */
 39:     return 0;
 40:   }
 41:   /* Always use FVM adjacency to create partitioner graph */
 42:   DMGetBasicAdjacency(dm, &useCone, &useClosure);
 43:   DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE);
 44:   /* Need overlap >= 1 */
 45:   DMPlexGetOverlap(dm, &overlap);
 46:   if (size && overlap < 1) {
 47:     DMPlexDistributeOverlap(dm, 1, NULL, &ovdm);
 48:   } else {
 49:     PetscObjectReference((PetscObject)dm);
 50:     ovdm = dm;
 51:   }
 52:   DMGetPointSF(ovdm, &sfPoint);
 53:   DMPlexGetHeightStratum(ovdm, height, &cStart, &cEnd);
 54:   DMPlexCreateNumbering_Plex(ovdm, cStart, cEnd, 0, NULL, sfPoint, &cellNumbering);
 55:   if (globalNumbering) {
 56:     PetscObjectReference((PetscObject)cellNumbering);
 57:     *globalNumbering = cellNumbering;
 58:   }
 59:   ISGetIndices(cellNumbering, &cellNum);
 60:   /* Determine sizes */
 61:   for (*numVertices = 0, c = cStart; c < cEnd; ++c) {
 62:     /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
 63:     if (cellNum[c - cStart] < 0) continue;
 64:     (*numVertices)++;
 65:   }
 66:   PetscCalloc1(*numVertices + 1, &vOffsets);
 67:   for (c = cStart, v = 0; c < cEnd; ++c) {
 68:     PetscInt adjSize = PETSC_DETERMINE, a, vsize = 0;

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

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

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

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

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

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

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

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

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

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

303:   if (offsets) *offsets = vOffsets;
304:   if (adjacency) *adjacency = graph;

306:   /* Cleanup */
307:   PetscSegBufferDestroy(&adjBuffer);
308:   PetscSectionDestroy(&section);
309:   ISRestoreIndices(cellNumbering, &cellNum);
310:   ISDestroy(&cellNumbering);
311:   return 0;
312: }

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

325:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
326:   DMGetDimension(dm, &dim);
327:   DMPlexGetDepth(dm, &depth);
328:   if (dim != depth) {
329:     /* We do not handle the uninterpolated case here */
330:     DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
331:     /* DMPlexCreateNeighborCSR does not make a numbering */
332:     if (globalNumbering) DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);
333:     /* Different behavior for empty graphs */
334:     if (!*numVertices) {
335:       PetscMalloc1(1, offsets);
336:       (*offsets)[0] = 0;
337:     }
338:     /* Broken in parallel */
340:     return 0;
341:   }
342:   /* Interpolated and parallel case */

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

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

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

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

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

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

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

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

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

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

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

433:   /* Get parallel CSR by doing conn^T * conn */
434:   MatTransposeMatMult(conn, conn, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CSR);
435:   MatDestroy(&conn);

437:   /* extract local part of the CSR */
438:   MatMPIAIJGetLocalMat(CSR, MAT_INITIAL_MATRIX, &conn);
439:   MatDestroy(&CSR);
440:   MatGetRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);

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

456:       for (j = ii[i]; j < ii[i + 1]; j++) {
457:         if (jj[j] == g) continue; /* again, self-connectivity */
458:         idxs[c++] = jj[j];
459:       }
460:     }
462:     ISRestoreIndices(cis_own, &rows);
463:     *adjacency = idxs;
464:   }

466:   /* cleanup */
467:   ISDestroy(&cis_own);
468:   MatRestoreRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
470:   MatDestroy(&conn);
471:   return 0;
472: }

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

477:   Collective on dm

479:   Input Parameters:
480: + dm      - The mesh `DM`
481: - height  - Height of the strata from which to construct the graph

483:   Output Parameters:
484: + numVertices     - Number of vertices in the graph
485: . offsets         - Point offsets in the graph
486: . adjacency       - Point connectivity in the graph
487: - 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.

489:   Options Database Keys:
490: . -dm_plex_csr_alg <mat,graph,overlap> - Choose the algorithm for computing the CSR graph

492:   Level: developer

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

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

504:   PetscOptionsGetEnum(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_csr_alg", DMPlexCSRAlgorithms, (PetscEnum *)&alg, NULL);
505:   switch (alg) {
506:   case DM_PLEX_CSR_MAT:
507:     DMPlexCreatePartitionerGraph_ViaMat(dm, height, numVertices, offsets, adjacency, globalNumbering);
508:     break;
509:   case DM_PLEX_CSR_GRAPH:
510:     DMPlexCreatePartitionerGraph_Native(dm, height, numVertices, offsets, adjacency, globalNumbering);
511:     break;
512:   case DM_PLEX_CSR_OVERLAP:
513:     DMPlexCreatePartitionerGraph_Overlap(dm, height, numVertices, offsets, adjacency, globalNumbering);
514:     break;
515:   }
516:   return 0;
517: }

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

522:   Collective on dm

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

528:   Output Parameters:
529: + numVertices - The number of local vertices in the graph, or cells in the mesh.
530: . offsets     - The offset to the adjacency list for each cell
531: - adjacency   - The adjacency list for all cells

533:   Level: advanced

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

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

549:   /* For parallel partitioning, I think you have to communicate supports */
550:   DMGetDimension(dm, &dim);
551:   cellDim = dim - cellHeight;
552:   DMPlexGetDepth(dm, &depth);
553:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
554:   if (cEnd - cStart == 0) {
555:     if (numVertices) *numVertices = 0;
556:     if (offsets) *offsets = NULL;
557:     if (adjacency) *adjacency = NULL;
558:     return 0;
559:   }
560:   numCells  = cEnd - cStart;
561:   faceDepth = depth - cellHeight;
562:   if (dim == depth) {
563:     PetscInt f, fStart, fEnd;

565:     PetscCalloc1(numCells + 1, &off);
566:     /* Count neighboring cells */
567:     DMPlexGetHeightStratum(dm, cellHeight + 1, &fStart, &fEnd);
568:     for (f = fStart; f < fEnd; ++f) {
569:       const PetscInt *support;
570:       PetscInt        supportSize;
571:       DMPlexGetSupportSize(dm, f, &supportSize);
572:       DMPlexGetSupport(dm, f, &support);
573:       if (supportSize == 2) {
574:         PetscInt numChildren;

576:         DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
577:         if (!numChildren) {
578:           ++off[support[0] - cStart + 1];
579:           ++off[support[1] - cStart + 1];
580:         }
581:       }
582:     }
583:     /* Prefix sum */
584:     for (c = 1; c <= numCells; ++c) off[c] += off[c - 1];
585:     if (adjacency) {
586:       PetscInt *tmp;

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

600:           DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
601:           if (!numChildren) {
602:             adj[tmp[support[0] - cStart]++] = support[1];
603:             adj[tmp[support[1] - cStart]++] = support[0];
604:           }
605:         }
606:       }
607:       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);
608:       PetscFree(tmp);
609:     }
610:     if (numVertices) *numVertices = numCells;
611:     if (offsets) *offsets = off;
612:     if (adjacency) *adjacency = adj;
613:     return 0;
614:   }
615:   /* Setup face recognition */
616:   if (faceDepth == 1) {
617:     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 */

619:     for (c = cStart; c < cEnd; ++c) {
620:       PetscInt corners;

622:       DMPlexGetConeSize(dm, c, &corners);
623:       if (!cornersSeen[corners]) {
624:         PetscInt nFV;

627:         cornersSeen[corners] = 1;

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

631:         numFaceVertices[numFaceCases++] = nFV;
632:       }
633:     }
634:   }
635:   PetscCalloc1(numCells + 1, &off);
636:   /* Count neighboring cells */
637:   for (cell = cStart; cell < cEnd; ++cell) {
638:     PetscInt numNeighbors = PETSC_DETERMINE, n;

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

648:       cellPair[0] = cell;
649:       cellPair[1] = neighborCells[n];
650:       if (cellPair[0] == cellPair[1]) continue;
651:       if (!found) {
652:         DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
653:         if (meetSize) {
654:           PetscInt f;

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

671:   if (adjacency) {
672:     PetscMalloc1(off[numCells], &adj);
673:     /* Get neighboring cells */
674:     for (cell = cStart; cell < cEnd; ++cell) {
675:       PetscInt numNeighbors = PETSC_DETERMINE, n;
676:       PetscInt cellOffset   = 0;

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

686:         cellPair[0] = cell;
687:         cellPair[1] = neighborCells[n];
688:         if (cellPair[0] == cellPair[1]) continue;
689:         if (!found) {
690:           DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
691:           if (meetSize) {
692:             PetscInt f;

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

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

720:   Collective on part

722:   Input Parameters:
723: + part    - The `PetscPartitioner`
724: . targetSection - The `PetscSection` describing the absolute weight of each partition (can be NULL)
725: - dm      - The mesh `DM`

727:   Output Parameters:
728: + partSection     - The `PetscSection` giving the division of points by partition
729: - partition       - The list of points by partition

731:   Level: developer

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

737: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `PetscPartitioner`, `PetscSection`, `DMPlexDistribute()`, `PetscPartitionerCreate()`, `PetscSectionCreate()`,
738:          `PetscSectionSetChart()`, `PetscPartitionerPartition()`
739: @*/
740: PetscErrorCode PetscPartitionerDMPlexPartition(PetscPartitioner part, DM dm, PetscSection targetSection, PetscSection partSection, IS *partition)
741: {
742:   PetscMPIInt  size;
743:   PetscBool    isplex;
744:   PetscSection vertSection = NULL;

751:   PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex);
753:   MPI_Comm_size(PetscObjectComm((PetscObject)part), &size);
754:   if (size == 1) {
755:     PetscInt *points;
756:     PetscInt  cStart, cEnd, c;

758:     DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
759:     PetscSectionReset(partSection);
760:     PetscSectionSetChart(partSection, 0, size);
761:     PetscSectionSetDof(partSection, 0, cEnd - cStart);
762:     PetscSectionSetUp(partSection);
763:     PetscMalloc1(cEnd - cStart, &points);
764:     for (c = cStart; c < cEnd; ++c) points[c] = c;
765:     ISCreateGeneral(PetscObjectComm((PetscObject)part), cEnd - cStart, points, PETSC_OWN_POINTER, partition);
766:     return 0;
767:   }
768:   if (part->height == 0) {
769:     PetscInt  numVertices = 0;
770:     PetscInt *start       = NULL;
771:     PetscInt *adjacency   = NULL;
772:     IS        globalNumbering;

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

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

792:       /* 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) */
793:       /* We do this only if the local section has been set */
794:       if (section) {
795:         PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, NULL);
796:         if (!clSection) DMPlexCreateClosureIndex(dm, NULL);
797:         PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, &clPoints);
798:         ISGetIndices(clPoints, &clIdx);
799:       }
800:       DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
801:       PetscSectionCreate(PETSC_COMM_SELF, &vertSection);
802:       PetscSectionSetChart(vertSection, 0, numVertices);
803:       if (globalNumbering) {
804:         ISGetIndices(globalNumbering, &gid);
805:       } else gid = NULL;
806:       for (p = pStart, v = 0; p < pEnd; ++p) {
807:         PetscInt dof = 1;

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

812:         if (section) {
813:           PetscInt cl, clSize, clOff;

815:           dof = 0;
816:           PetscSectionGetDof(clSection, p, &clSize);
817:           PetscSectionGetOffset(clSection, p, &clOff);
818:           for (cl = 0; cl < clSize; cl += 2) {
819:             PetscInt clDof, clPoint = clIdx[clOff + cl]; /* odd indices are reserved for orientations */

821:             PetscSectionGetDof(section, clPoint, &clDof);
822:             dof += clDof;
823:           }
824:         }
826:         PetscSectionSetDof(vertSection, v, dof);
827:         v++;
828:       }
829:       if (globalNumbering) ISRestoreIndices(globalNumbering, &gid);
830:       if (clPoints) ISRestoreIndices(clPoints, &clIdx);
831:       PetscSectionSetUp(vertSection);
832:     }
833:     PetscPartitionerPartition(part, size, numVertices, start, adjacency, vertSection, targetSection, partSection, partition);
834:     PetscFree(start);
835:     PetscFree(adjacency);
836:     if (globalNumbering) { /* partition is wrt global unique numbering: change this to be wrt local numbering */
837:       const PetscInt *globalNum;
838:       const PetscInt *partIdx;
839:       PetscInt       *map, cStart, cEnd;
840:       PetscInt       *adjusted, i, localSize, offset;
841:       IS              newPartition;

843:       ISGetLocalSize(*partition, &localSize);
844:       PetscMalloc1(localSize, &adjusted);
845:       ISGetIndices(globalNumbering, &globalNum);
846:       ISGetIndices(*partition, &partIdx);
847:       PetscMalloc1(localSize, &map);
848:       DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
849:       for (i = cStart, offset = 0; i < cEnd; i++) {
850:         if (globalNum[i - cStart] >= 0) map[offset++] = i;
851:       }
852:       for (i = 0; i < localSize; i++) adjusted[i] = map[partIdx[i]];
853:       PetscFree(map);
854:       ISRestoreIndices(*partition, &partIdx);
855:       ISRestoreIndices(globalNumbering, &globalNum);
856:       ISCreateGeneral(PETSC_COMM_SELF, localSize, adjusted, PETSC_OWN_POINTER, &newPartition);
857:       ISDestroy(&globalNumbering);
858:       ISDestroy(partition);
859:       *partition = newPartition;
860:     }
861:   } else SETERRQ(PetscObjectComm((PetscObject)part), PETSC_ERR_ARG_OUTOFRANGE, "Invalid height %" PetscInt_FMT " for points to partition", part->height);
862:   PetscSectionDestroy(&vertSection);
863:   return 0;
864: }

866: /*@
867:   DMPlexGetPartitioner - Get the mesh partitioner

869:   Not collective

871:   Input Parameter:
872: . dm - The `DM`

874:   Output Parameter:
875: . part - The `PetscPartitioner`

877:   Level: developer

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

882: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `PetscPartitioner`, `PetscSection`, `DMPlexDistribute()`, `DMPlexSetPartitioner()`, `PetscPartitionerDMPlexPartition()`, `PetscPartitionerCreate()`
883: @*/
884: PetscErrorCode DMPlexGetPartitioner(DM dm, PetscPartitioner *part)
885: {
886:   DM_Plex *mesh = (DM_Plex *)dm->data;

890:   *part = mesh->partitioner;
891:   return 0;
892: }

894: /*@
895:   DMPlexSetPartitioner - Set the mesh partitioner

897:   logically collective on dm

899:   Input Parameters:
900: + dm - The `DM`
901: - part - The partitioner

903:   Level: developer

905:   Note:
906:   Any existing `PetscPartitioner` will be destroyed.

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

916:   PetscObjectReference((PetscObject)part);
917:   PetscPartitionerDestroy(&mesh->partitioner);
918:   mesh->partitioner = part;
919:   return 0;
920: }

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

929:   PetscHSetIQueryAdd(ht, point, &missing);
930:   if (missing) {
931:     DMPlexGetCone(dm, point, &cone);
932:     DMPlexGetConeSize(dm, point, &coneSize);
933:     for (c = 0; c < coneSize; c++) DMPlexAddClosure_Private(dm, ht, cone[c]);
934:   }
935:   return 0;
936: }

938: PETSC_UNUSED static PetscErrorCode DMPlexAddClosure_Tree(DM dm, PetscHSetI ht, PetscInt point, PetscBool up, PetscBool down)
939: {
940:   if (up) {
941:     PetscInt parent;

943:     DMPlexGetTreeParent(dm, point, &parent, NULL);
944:     if (parent != point) {
945:       PetscInt closureSize, *closure = NULL, i;

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

951:         PetscHSetIAdd(ht, cpoint);
952:         DMPlexAddClosure_Tree(dm, ht, cpoint, PETSC_TRUE, PETSC_FALSE);
953:       }
954:       DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure);
955:     }
956:   }
957:   if (down) {
958:     PetscInt        numChildren;
959:     const PetscInt *children;

961:     DMPlexGetTreeChildren(dm, point, &numChildren, &children);
962:     if (numChildren) {
963:       PetscInt i;

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

968:         PetscHSetIAdd(ht, cpoint);
969:         DMPlexAddClosure_Tree(dm, ht, cpoint, PETSC_FALSE, PETSC_TRUE);
970:       }
971:     }
972:   }
973:   return 0;
974: }

976: static PetscErrorCode DMPlexAddClosureTree_Up_Private(DM dm, PetscHSetI ht, PetscInt point)
977: {
978:   PetscInt parent;

981:   DMPlexGetTreeParent(dm, point, &parent, NULL);
982:   if (point != parent) {
983:     const PetscInt *cone;
984:     PetscInt        coneSize, c;

986:     DMPlexAddClosureTree_Up_Private(dm, ht, parent);
987:     DMPlexAddClosure_Private(dm, ht, parent);
988:     DMPlexGetCone(dm, parent, &cone);
989:     DMPlexGetConeSize(dm, parent, &coneSize);
990:     for (c = 0; c < coneSize; c++) {
991:       const PetscInt cp = cone[c];

993:       DMPlexAddClosureTree_Up_Private(dm, ht, cp);
994:     }
995:   }
996:   return 0;
997: }

999: static PetscErrorCode DMPlexAddClosureTree_Down_Private(DM dm, PetscHSetI ht, PetscInt point)
1000: {
1001:   PetscInt        i, numChildren;
1002:   const PetscInt *children;

1005:   DMPlexGetTreeChildren(dm, point, &numChildren, &children);
1006:   for (i = 0; i < numChildren; i++) PetscHSetIAdd(ht, children[i]);
1007:   return 0;
1008: }

1010: static PetscErrorCode DMPlexAddClosureTree_Private(DM dm, PetscHSetI ht, PetscInt point)
1011: {
1012:   const PetscInt *cone;
1013:   PetscInt        coneSize, c;

1016:   PetscHSetIAdd(ht, point);
1017:   DMPlexAddClosureTree_Up_Private(dm, ht, point);
1018:   DMPlexAddClosureTree_Down_Private(dm, ht, point);
1019:   DMPlexGetCone(dm, point, &cone);
1020:   DMPlexGetConeSize(dm, point, &coneSize);
1021:   for (c = 0; c < coneSize; c++) DMPlexAddClosureTree_Private(dm, ht, cone[c]);
1022:   return 0;
1023: }

1025: PetscErrorCode DMPlexClosurePoints_Private(DM dm, PetscInt numPoints, const PetscInt points[], IS *closureIS)
1026: {
1027:   DM_Plex        *mesh    = (DM_Plex *)dm->data;
1028:   const PetscBool hasTree = (mesh->parentSection || mesh->childSection) ? PETSC_TRUE : PETSC_FALSE;
1029:   PetscInt        nelems, *elems, off = 0, p;
1030:   PetscHSetI      ht = NULL;

1032:   PetscHSetICreate(&ht);
1033:   PetscHSetIResize(ht, numPoints * 16);
1034:   if (!hasTree) {
1035:     for (p = 0; p < numPoints; ++p) DMPlexAddClosure_Private(dm, ht, points[p]);
1036:   } else {
1037: #if 1
1038:     for (p = 0; p < numPoints; ++p) DMPlexAddClosureTree_Private(dm, ht, points[p]);
1039: #else
1040:     PetscInt *closure = NULL, closureSize, c;
1041:     for (p = 0; p < numPoints; ++p) {
1042:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
1043:       for (c = 0; c < closureSize * 2; c += 2) {
1044:         PetscHSetIAdd(ht, closure[c]);
1045:         if (hasTree) DMPlexAddClosure_Tree(dm, ht, closure[c], PETSC_TRUE, PETSC_TRUE);
1046:       }
1047:     }
1048:     if (closure) DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure);
1049: #endif
1050:   }
1051:   PetscHSetIGetSize(ht, &nelems);
1052:   PetscMalloc1(nelems, &elems);
1053:   PetscHSetIGetElems(ht, &off, elems);
1054:   PetscHSetIDestroy(&ht);
1055:   PetscSortInt(nelems, elems);
1056:   ISCreateGeneral(PETSC_COMM_SELF, nelems, elems, PETSC_OWN_POINTER, closureIS);
1057:   return 0;
1058: }

1060: /*@
1061:   DMPlexPartitionLabelClosure - Add the closure of all points to the partition label

1063:   Input Parameters:
1064: + dm     - The `DM`
1065: - label  - `DMLabel` assigning ranks to remote roots

1067:   Level: developer

1069: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1070: @*/
1071: PetscErrorCode DMPlexPartitionLabelClosure(DM dm, DMLabel label)
1072: {
1073:   IS              rankIS, pointIS, closureIS;
1074:   const PetscInt *ranks, *points;
1075:   PetscInt        numRanks, numPoints, r;

1077:   DMLabelGetValueIS(label, &rankIS);
1078:   ISGetLocalSize(rankIS, &numRanks);
1079:   ISGetIndices(rankIS, &ranks);
1080:   for (r = 0; r < numRanks; ++r) {
1081:     const PetscInt rank = ranks[r];
1082:     DMLabelGetStratumIS(label, rank, &pointIS);
1083:     ISGetLocalSize(pointIS, &numPoints);
1084:     ISGetIndices(pointIS, &points);
1085:     DMPlexClosurePoints_Private(dm, numPoints, points, &closureIS);
1086:     ISRestoreIndices(pointIS, &points);
1087:     ISDestroy(&pointIS);
1088:     DMLabelSetStratumIS(label, rank, closureIS);
1089:     ISDestroy(&closureIS);
1090:   }
1091:   ISRestoreIndices(rankIS, &ranks);
1092:   ISDestroy(&rankIS);
1093:   return 0;
1094: }

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

1099:   Input Parameters:
1100: + dm     - The `DM`
1101: - label  - `DMLabel` assigning ranks to remote roots

1103:   Level: developer

1105: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1106: @*/
1107: PetscErrorCode DMPlexPartitionLabelAdjacency(DM dm, DMLabel label)
1108: {
1109:   IS              rankIS, pointIS;
1110:   const PetscInt *ranks, *points;
1111:   PetscInt        numRanks, numPoints, r, p, a, adjSize;
1112:   PetscInt       *adj = NULL;

1114:   DMLabelGetValueIS(label, &rankIS);
1115:   ISGetLocalSize(rankIS, &numRanks);
1116:   ISGetIndices(rankIS, &ranks);
1117:   for (r = 0; r < numRanks; ++r) {
1118:     const PetscInt rank = ranks[r];

1120:     DMLabelGetStratumIS(label, rank, &pointIS);
1121:     ISGetLocalSize(pointIS, &numPoints);
1122:     ISGetIndices(pointIS, &points);
1123:     for (p = 0; p < numPoints; ++p) {
1124:       adjSize = PETSC_DETERMINE;
1125:       DMPlexGetAdjacency(dm, points[p], &adjSize, &adj);
1126:       for (a = 0; a < adjSize; ++a) DMLabelSetValue(label, adj[a], rank);
1127:     }
1128:     ISRestoreIndices(pointIS, &points);
1129:     ISDestroy(&pointIS);
1130:   }
1131:   ISRestoreIndices(rankIS, &ranks);
1132:   ISDestroy(&rankIS);
1133:   PetscFree(adj);
1134:   return 0;
1135: }

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

1140:   Input Parameters:
1141: + dm     - The `DM`
1142: - label  - `DMLabel` assigning ranks to remote roots

1144:   Level: developer

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

1150: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1151: @*/
1152: PetscErrorCode DMPlexPartitionLabelPropagate(DM dm, DMLabel label)
1153: {
1154:   MPI_Comm        comm;
1155:   PetscMPIInt     rank;
1156:   PetscSF         sfPoint;
1157:   DMLabel         lblRoots, lblLeaves;
1158:   IS              rankIS, pointIS;
1159:   const PetscInt *ranks;
1160:   PetscInt        numRanks, r;

1162:   PetscObjectGetComm((PetscObject)dm, &comm);
1163:   MPI_Comm_rank(comm, &rank);
1164:   DMGetPointSF(dm, &sfPoint);
1165:   /* Pull point contributions from remote leaves into local roots */
1166:   DMLabelGather(label, sfPoint, &lblLeaves);
1167:   DMLabelGetValueIS(lblLeaves, &rankIS);
1168:   ISGetLocalSize(rankIS, &numRanks);
1169:   ISGetIndices(rankIS, &ranks);
1170:   for (r = 0; r < numRanks; ++r) {
1171:     const PetscInt remoteRank = ranks[r];
1172:     if (remoteRank == rank) continue;
1173:     DMLabelGetStratumIS(lblLeaves, remoteRank, &pointIS);
1174:     DMLabelInsertIS(label, pointIS, remoteRank);
1175:     ISDestroy(&pointIS);
1176:   }
1177:   ISRestoreIndices(rankIS, &ranks);
1178:   ISDestroy(&rankIS);
1179:   DMLabelDestroy(&lblLeaves);
1180:   /* Push point contributions from roots into remote leaves */
1181:   DMLabelDistribute(label, sfPoint, &lblRoots);
1182:   DMLabelGetValueIS(lblRoots, &rankIS);
1183:   ISGetLocalSize(rankIS, &numRanks);
1184:   ISGetIndices(rankIS, &ranks);
1185:   for (r = 0; r < numRanks; ++r) {
1186:     const PetscInt remoteRank = ranks[r];
1187:     if (remoteRank == rank) continue;
1188:     DMLabelGetStratumIS(lblRoots, remoteRank, &pointIS);
1189:     DMLabelInsertIS(label, pointIS, remoteRank);
1190:     ISDestroy(&pointIS);
1191:   }
1192:   ISRestoreIndices(rankIS, &ranks);
1193:   ISDestroy(&rankIS);
1194:   DMLabelDestroy(&lblRoots);
1195:   return 0;
1196: }

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

1201:   Input Parameters:
1202: + dm        - The `DM`
1203: . rootLabel - `DMLabel` assigning ranks to local roots
1204: - processSF - A star forest mapping into the local index on each remote rank

1206:   Output Parameter:
1207: . leafLabel - `DMLabel `assigning ranks to remote roots

1209:   Level: developer

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

1215: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1216: @*/
1217: PetscErrorCode DMPlexPartitionLabelInvert(DM dm, DMLabel rootLabel, PetscSF processSF, DMLabel leafLabel)
1218: {
1219:   MPI_Comm           comm;
1220:   PetscMPIInt        rank, size, r;
1221:   PetscInt           p, n, numNeighbors, numPoints, dof, off, rootSize, l, nleaves, leafSize;
1222:   PetscSF            sfPoint;
1223:   PetscSection       rootSection;
1224:   PetscSFNode       *rootPoints, *leafPoints;
1225:   const PetscSFNode *remote;
1226:   const PetscInt    *local, *neighbors;
1227:   IS                 valueIS;
1228:   PetscBool          mpiOverflow = PETSC_FALSE;

1230:   PetscLogEventBegin(DMPLEX_PartLabelInvert, dm, 0, 0, 0);
1231:   PetscObjectGetComm((PetscObject)dm, &comm);
1232:   MPI_Comm_rank(comm, &rank);
1233:   MPI_Comm_size(comm, &size);
1234:   DMGetPointSF(dm, &sfPoint);

1236:   /* Convert to (point, rank) and use actual owners */
1237:   PetscSectionCreate(comm, &rootSection);
1238:   PetscSectionSetChart(rootSection, 0, size);
1239:   DMLabelGetValueIS(rootLabel, &valueIS);
1240:   ISGetLocalSize(valueIS, &numNeighbors);
1241:   ISGetIndices(valueIS, &neighbors);
1242:   for (n = 0; n < numNeighbors; ++n) {
1243:     DMLabelGetStratumSize(rootLabel, neighbors[n], &numPoints);
1244:     PetscSectionAddDof(rootSection, neighbors[n], numPoints);
1245:   }
1246:   PetscSectionSetUp(rootSection);
1247:   PetscSectionGetStorageSize(rootSection, &rootSize);
1248:   PetscMalloc1(rootSize, &rootPoints);
1249:   PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote);
1250:   for (n = 0; n < numNeighbors; ++n) {
1251:     IS              pointIS;
1252:     const PetscInt *points;

1254:     PetscSectionGetOffset(rootSection, neighbors[n], &off);
1255:     DMLabelGetStratumIS(rootLabel, neighbors[n], &pointIS);
1256:     ISGetLocalSize(pointIS, &numPoints);
1257:     ISGetIndices(pointIS, &points);
1258:     for (p = 0; p < numPoints; ++p) {
1259:       if (local) PetscFindInt(points[p], nleaves, local, &l);
1260:       else l = -1;
1261:       if (l >= 0) {
1262:         rootPoints[off + p] = remote[l];
1263:       } else {
1264:         rootPoints[off + p].index = points[p];
1265:         rootPoints[off + p].rank  = rank;
1266:       }
1267:     }
1268:     ISRestoreIndices(pointIS, &points);
1269:     ISDestroy(&pointIS);
1270:   }

1272:   /* Try to communicate overlap using All-to-All */
1273:   if (!processSF) {
1274:     PetscInt64   counter     = 0;
1275:     PetscBool    locOverflow = PETSC_FALSE;
1276:     PetscMPIInt *scounts, *sdispls, *rcounts, *rdispls;

1278:     PetscCalloc4(size, &scounts, size, &sdispls, size, &rcounts, size, &rdispls);
1279:     for (n = 0; n < numNeighbors; ++n) {
1280:       PetscSectionGetDof(rootSection, neighbors[n], &dof);
1281:       PetscSectionGetOffset(rootSection, neighbors[n], &off);
1282: #if defined(PETSC_USE_64BIT_INDICES)
1283:       if (dof > PETSC_MPI_INT_MAX) {
1284:         locOverflow = PETSC_TRUE;
1285:         break;
1286:       }
1287:       if (off > PETSC_MPI_INT_MAX) {
1288:         locOverflow = PETSC_TRUE;
1289:         break;
1290:       }
1291: #endif
1292:       scounts[neighbors[n]] = (PetscMPIInt)dof;
1293:       sdispls[neighbors[n]] = (PetscMPIInt)off;
1294:     }
1295:     MPI_Alltoall(scounts, 1, MPI_INT, rcounts, 1, MPI_INT, comm);
1296:     for (r = 0; r < size; ++r) {
1297:       rdispls[r] = (int)counter;
1298:       counter += rcounts[r];
1299:     }
1300:     if (counter > PETSC_MPI_INT_MAX) locOverflow = PETSC_TRUE;
1301:     MPI_Allreduce(&locOverflow, &mpiOverflow, 1, MPIU_BOOL, MPI_LOR, comm);
1302:     if (!mpiOverflow) {
1303:       PetscInfo(dm, "Using Alltoallv for mesh distribution\n");
1304:       leafSize = (PetscInt)counter;
1305:       PetscMalloc1(leafSize, &leafPoints);
1306:       MPI_Alltoallv(rootPoints, scounts, sdispls, MPIU_2INT, leafPoints, rcounts, rdispls, MPIU_2INT, comm);
1307:     }
1308:     PetscFree4(scounts, sdispls, rcounts, rdispls);
1309:   }

1311:   /* Communicate overlap using process star forest */
1312:   if (processSF || mpiOverflow) {
1313:     PetscSF      procSF;
1314:     PetscSection leafSection;

1316:     if (processSF) {
1317:       PetscInfo(dm, "Using processSF for mesh distribution\n");
1318:       PetscObjectReference((PetscObject)processSF);
1319:       procSF = processSF;
1320:     } else {
1321:       PetscInfo(dm, "Using processSF for mesh distribution (MPI overflow)\n");
1322:       PetscSFCreate(comm, &procSF);
1323:       PetscSFSetGraphWithPattern(procSF, NULL, PETSCSF_PATTERN_ALLTOALL);
1324:     }

1326:     PetscSectionCreate(PetscObjectComm((PetscObject)dm), &leafSection);
1327:     DMPlexDistributeData(dm, procSF, rootSection, MPIU_2INT, rootPoints, leafSection, (void **)&leafPoints);
1328:     PetscSectionGetStorageSize(leafSection, &leafSize);
1329:     PetscSectionDestroy(&leafSection);
1330:     PetscSFDestroy(&procSF);
1331:   }

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

1335:   ISRestoreIndices(valueIS, &neighbors);
1336:   ISDestroy(&valueIS);
1337:   PetscSectionDestroy(&rootSection);
1338:   PetscFree(rootPoints);
1339:   PetscFree(leafPoints);
1340:   PetscLogEventEnd(DMPLEX_PartLabelInvert, dm, 0, 0, 0);
1341:   return 0;
1342: }

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

1347:   Input Parameters:
1348: + dm    - The `DM`
1349: - label - `DMLabel` assigning ranks to remote roots

1351:   Output Parameter:
1352: . sf    - The star forest communication context encapsulating the defined mapping

1354:   Level: developer

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

1359: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMLabel`, `PetscSF`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1360: @*/
1361: PetscErrorCode DMPlexPartitionLabelCreateSF(DM dm, DMLabel label, PetscSF *sf)
1362: {
1363:   PetscMPIInt     rank;
1364:   PetscInt        n, numRemote, p, numPoints, pStart, pEnd, idx = 0, nNeighbors;
1365:   PetscSFNode    *remotePoints;
1366:   IS              remoteRootIS, neighborsIS;
1367:   const PetscInt *remoteRoots, *neighbors;

1369:   PetscLogEventBegin(DMPLEX_PartLabelCreateSF, dm, 0, 0, 0);
1370:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);

1372:   DMLabelGetValueIS(label, &neighborsIS);
1373: #if 0
1374:   {
1375:     IS is;
1376:     ISDuplicate(neighborsIS, &is);
1377:     ISSort(is);
1378:     ISDestroy(&neighborsIS);
1379:     neighborsIS = is;
1380:   }
1381: #endif
1382:   ISGetLocalSize(neighborsIS, &nNeighbors);
1383:   ISGetIndices(neighborsIS, &neighbors);
1384:   for (numRemote = 0, n = 0; n < nNeighbors; ++n) {
1385:     DMLabelGetStratumSize(label, neighbors[n], &numPoints);
1386:     numRemote += numPoints;
1387:   }
1388:   PetscMalloc1(numRemote, &remotePoints);
1389:   /* Put owned points first */
1390:   DMLabelGetStratumSize(label, rank, &numPoints);
1391:   if (numPoints > 0) {
1392:     DMLabelGetStratumIS(label, rank, &remoteRootIS);
1393:     ISGetIndices(remoteRootIS, &remoteRoots);
1394:     for (p = 0; p < numPoints; p++) {
1395:       remotePoints[idx].index = remoteRoots[p];
1396:       remotePoints[idx].rank  = rank;
1397:       idx++;
1398:     }
1399:     ISRestoreIndices(remoteRootIS, &remoteRoots);
1400:     ISDestroy(&remoteRootIS);
1401:   }
1402:   /* Now add remote points */
1403:   for (n = 0; n < nNeighbors; ++n) {
1404:     const PetscInt nn = neighbors[n];

1406:     DMLabelGetStratumSize(label, nn, &numPoints);
1407:     if (nn == rank || numPoints <= 0) continue;
1408:     DMLabelGetStratumIS(label, nn, &remoteRootIS);
1409:     ISGetIndices(remoteRootIS, &remoteRoots);
1410:     for (p = 0; p < numPoints; p++) {
1411:       remotePoints[idx].index = remoteRoots[p];
1412:       remotePoints[idx].rank  = nn;
1413:       idx++;
1414:     }
1415:     ISRestoreIndices(remoteRootIS, &remoteRoots);
1416:     ISDestroy(&remoteRootIS);
1417:   }
1418:   PetscSFCreate(PetscObjectComm((PetscObject)dm), sf);
1419:   DMPlexGetChart(dm, &pStart, &pEnd);
1420:   PetscSFSetGraph(*sf, pEnd - pStart, numRemote, NULL, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
1421:   PetscSFSetUp(*sf);
1422:   ISDestroy(&neighborsIS);
1423:   PetscLogEventEnd(DMPLEX_PartLabelCreateSF, dm, 0, 0, 0);
1424:   return 0;
1425: }

1427: #if defined(PETSC_HAVE_PARMETIS)
1428:   #include <parmetis.h>
1429: #endif

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

1437:   DMPlexRewriteSF - Rewrites the ownership of the `PetsSF` of a `DM` (in place).

1439:   Input parameters:
1440: + dm                - The `DMPLEX` object.
1441: . n                 - The number of points.
1442: . pointsToRewrite   - The points in the `PetscSF` whose ownership will change.
1443: . targetOwners      - New owner for each element in pointsToRewrite.
1444: - degrees           - Degrees of the points in the `PetscSF` as obtained by `PetscSFComputeDegreeBegin()`/`PetscSFComputeDegreeEnd()`.

1446:   Level: developer

1448: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMLabel`, `PetscSF`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1449: @*/
1450: static PetscErrorCode DMPlexRewriteSF(DM dm, PetscInt n, PetscInt *pointsToRewrite, PetscInt *targetOwners, const PetscInt *degrees)
1451: {
1452:   PetscInt           pStart, pEnd, i, j, counter, leafCounter, sumDegrees, nroots, nleafs;
1453:   PetscInt          *cumSumDegrees, *newOwners, *newNumbers, *rankOnLeafs, *locationsOfLeafs, *remoteLocalPointOfLeafs, *points, *leafsNew;
1454:   PetscSFNode       *leafLocationsNew;
1455:   const PetscSFNode *iremote;
1456:   const PetscInt    *ilocal;
1457:   PetscBool         *isLeaf;
1458:   PetscSF            sf;
1459:   MPI_Comm           comm;
1460:   PetscMPIInt        rank, size;

1462:   PetscObjectGetComm((PetscObject)dm, &comm);
1463:   MPI_Comm_rank(comm, &rank);
1464:   MPI_Comm_size(comm, &size);
1465:   DMPlexGetChart(dm, &pStart, &pEnd);

1467:   DMGetPointSF(dm, &sf);
1468:   PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
1469:   PetscMalloc1(pEnd - pStart, &isLeaf);
1470:   for (i = 0; i < pEnd - pStart; i++) isLeaf[i] = PETSC_FALSE;
1471:   for (i = 0; i < nleafs; i++) isLeaf[ilocal[i] - pStart] = PETSC_TRUE;

1473:   PetscMalloc1(pEnd - pStart + 1, &cumSumDegrees);
1474:   cumSumDegrees[0] = 0;
1475:   for (i = 1; i <= pEnd - pStart; i++) cumSumDegrees[i] = cumSumDegrees[i - 1] + degrees[i - 1];
1476:   sumDegrees = cumSumDegrees[pEnd - pStart];
1477:   /* get the location of my leafs (we have sumDegrees many leafs pointing at our roots) */

1479:   PetscMalloc1(sumDegrees, &locationsOfLeafs);
1480:   PetscMalloc1(pEnd - pStart, &rankOnLeafs);
1481:   for (i = 0; i < pEnd - pStart; i++) rankOnLeafs[i] = rank;
1482:   PetscSFGatherBegin(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
1483:   PetscSFGatherEnd(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
1484:   PetscFree(rankOnLeafs);

1486:   /* get the remote local points of my leaves */
1487:   PetscMalloc1(sumDegrees, &remoteLocalPointOfLeafs);
1488:   PetscMalloc1(pEnd - pStart, &points);
1489:   for (i = 0; i < pEnd - pStart; i++) points[i] = pStart + i;
1490:   PetscSFGatherBegin(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
1491:   PetscSFGatherEnd(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
1492:   PetscFree(points);
1493:   /* Figure out the new owners of the vertices that are up for grabs and their numbers on the new owners */
1494:   PetscMalloc1(pEnd - pStart, &newOwners);
1495:   PetscMalloc1(pEnd - pStart, &newNumbers);
1496:   for (i = 0; i < pEnd - pStart; i++) {
1497:     newOwners[i]  = -1;
1498:     newNumbers[i] = -1;
1499:   }
1500:   {
1501:     PetscInt oldNumber, newNumber, oldOwner, newOwner;
1502:     for (i = 0; i < n; i++) {
1503:       oldNumber = pointsToRewrite[i];
1504:       newNumber = -1;
1505:       oldOwner  = rank;
1506:       newOwner  = targetOwners[i];
1507:       if (oldOwner == newOwner) {
1508:         newNumber = oldNumber;
1509:       } else {
1510:         for (j = 0; j < degrees[oldNumber]; j++) {
1511:           if (locationsOfLeafs[cumSumDegrees[oldNumber] + j] == newOwner) {
1512:             newNumber = remoteLocalPointOfLeafs[cumSumDegrees[oldNumber] + j];
1513:             break;
1514:           }
1515:         }
1516:       }

1519:       newOwners[oldNumber]  = newOwner;
1520:       newNumbers[oldNumber] = newNumber;
1521:     }
1522:   }
1523:   PetscFree(cumSumDegrees);
1524:   PetscFree(locationsOfLeafs);
1525:   PetscFree(remoteLocalPointOfLeafs);

1527:   PetscSFBcastBegin(sf, MPIU_INT, newOwners, newOwners, MPI_REPLACE);
1528:   PetscSFBcastEnd(sf, MPIU_INT, newOwners, newOwners, MPI_REPLACE);
1529:   PetscSFBcastBegin(sf, MPIU_INT, newNumbers, newNumbers, MPI_REPLACE);
1530:   PetscSFBcastEnd(sf, MPIU_INT, newNumbers, newNumbers, MPI_REPLACE);

1532:   /* Now count how many leafs we have on each processor. */
1533:   leafCounter = 0;
1534:   for (i = 0; i < pEnd - pStart; i++) {
1535:     if (newOwners[i] >= 0) {
1536:       if (newOwners[i] != rank) leafCounter++;
1537:     } else {
1538:       if (isLeaf[i]) leafCounter++;
1539:     }
1540:   }

1542:   /* Now set up the new sf by creating the leaf arrays */
1543:   PetscMalloc1(leafCounter, &leafsNew);
1544:   PetscMalloc1(leafCounter, &leafLocationsNew);

1546:   leafCounter = 0;
1547:   counter     = 0;
1548:   for (i = 0; i < pEnd - pStart; i++) {
1549:     if (newOwners[i] >= 0) {
1550:       if (newOwners[i] != rank) {
1551:         leafsNew[leafCounter]               = i;
1552:         leafLocationsNew[leafCounter].rank  = newOwners[i];
1553:         leafLocationsNew[leafCounter].index = newNumbers[i];
1554:         leafCounter++;
1555:       }
1556:     } else {
1557:       if (isLeaf[i]) {
1558:         leafsNew[leafCounter]               = i;
1559:         leafLocationsNew[leafCounter].rank  = iremote[counter].rank;
1560:         leafLocationsNew[leafCounter].index = iremote[counter].index;
1561:         leafCounter++;
1562:       }
1563:     }
1564:     if (isLeaf[i]) counter++;
1565:   }

1567:   PetscSFSetGraph(sf, nroots, leafCounter, leafsNew, PETSC_OWN_POINTER, leafLocationsNew, PETSC_OWN_POINTER);
1568:   PetscFree(newOwners);
1569:   PetscFree(newNumbers);
1570:   PetscFree(isLeaf);
1571:   return 0;
1572: }

1574: static PetscErrorCode DMPlexViewDistribution(MPI_Comm comm, PetscInt n, PetscInt skip, PetscInt *vtxwgt, PetscInt *part, PetscViewer viewer)
1575: {
1576:   PetscInt   *distribution, min, max, sum;
1577:   PetscMPIInt rank, size;

1579:   MPI_Comm_size(comm, &size);
1580:   MPI_Comm_rank(comm, &rank);
1581:   PetscCalloc1(size, &distribution);
1582:   for (PetscInt i = 0; i < n; i++) {
1583:     if (part) distribution[part[i]] += vtxwgt[skip * i];
1584:     else distribution[rank] += vtxwgt[skip * i];
1585:   }
1586:   MPI_Allreduce(MPI_IN_PLACE, distribution, size, MPIU_INT, MPI_SUM, comm);
1587:   min = distribution[0];
1588:   max = distribution[0];
1589:   sum = distribution[0];
1590:   for (PetscInt i = 1; i < size; i++) {
1591:     if (distribution[i] < min) min = distribution[i];
1592:     if (distribution[i] > max) max = distribution[i];
1593:     sum += distribution[i];
1594:   }
1595:   PetscViewerASCIIPrintf(viewer, "Min: %" PetscInt_FMT ", Avg: %" PetscInt_FMT ", Max: %" PetscInt_FMT ", Balance: %f\n", min, sum / size, max, (max * 1. * size) / sum);
1596:   PetscFree(distribution);
1597:   return 0;
1598: }

1600: #endif

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

1605:   Input parameters:
1606: + dm               - The `DMPLEX` object.
1607: . entityDepth      - depth of the entity to balance (0 -> balance vertices).
1608: . useInitialGuess  - whether to use the current distribution as initial guess (only used by ParMETIS).
1609: - parallel         - whether to use ParMETIS and do the partition in parallel or whether to gather the graph onto a single process and use METIS.

1611:   Output parameters:
1612: . success          - whether the graph partitioning was successful or not, optional. Unsuccessful simply means no change to the partitioning

1614:   Options Database Keys:
1615: +  -dm_plex_rebalance_shared_points_parmetis - Use ParMetis instead of Metis for the partitioner
1616: .  -dm_plex_rebalance_shared_points_use_initial_guess - Use current partition to bootstrap ParMetis partition
1617: .  -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_
1618: -  -dm_plex_rebalance_shared_points_monitor - Monitor the shared points rebalance process

1620:   Level: intermediate

1622: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1623: @*/
1624: PetscErrorCode DMPlexRebalanceSharedPoints(DM dm, PetscInt entityDepth, PetscBool useInitialGuess, PetscBool parallel, PetscBool *success)
1625: {
1626: #if defined(PETSC_HAVE_PARMETIS)
1627:   PetscSF            sf;
1628:   PetscInt           ierr, i, j, idx, jdx;
1629:   PetscInt           eBegin, eEnd, nroots, nleafs, pStart, pEnd;
1630:   const PetscInt    *degrees, *ilocal;
1631:   const PetscSFNode *iremote;
1632:   PetscBool         *toBalance, *isLeaf, *isExclusivelyOwned, *isNonExclusivelyOwned;
1633:   PetscInt           numExclusivelyOwned, numNonExclusivelyOwned;
1634:   PetscMPIInt        rank, size;
1635:   PetscInt          *globalNumbersOfLocalOwnedVertices, *leafGlobalNumbers;
1636:   const PetscInt    *cumSumVertices;
1637:   PetscInt           offset, counter;
1638:   PetscInt          *vtxwgt;
1639:   const PetscInt    *xadj, *adjncy;
1640:   PetscInt          *part, *options;
1641:   PetscInt           nparts, wgtflag, numflag, ncon, edgecut;
1642:   real_t            *ubvec;
1643:   PetscInt          *firstVertices, *renumbering;
1644:   PetscInt           failed, failedGlobal;
1645:   MPI_Comm           comm;
1646:   Mat                A;
1647:   PetscViewer        viewer;
1648:   PetscViewerFormat  format;
1649:   PetscLayout        layout;
1650:   real_t            *tpwgts;
1651:   PetscMPIInt       *counts, *mpiCumSumVertices;
1652:   PetscInt          *pointsToRewrite;
1653:   PetscInt           numRows;
1654:   PetscBool          done, usematpartitioning = PETSC_FALSE;
1655:   IS                 ispart = NULL;
1656:   MatPartitioning    mp;
1657:   const char        *prefix;

1659:   PetscObjectGetComm((PetscObject)dm, &comm);
1660:   MPI_Comm_size(comm, &size);
1661:   if (size == 1) {
1662:     if (success) *success = PETSC_TRUE;
1663:     return 0;
1664:   }
1665:   if (success) *success = PETSC_FALSE;
1666:   MPI_Comm_rank(comm, &rank);

1668:   parallel        = PETSC_FALSE;
1669:   useInitialGuess = PETSC_FALSE;
1670:   PetscObjectOptionsBegin((PetscObject)dm);
1671:   PetscOptionsName("-dm_plex_rebalance_shared_points_parmetis", "Use ParMetis instead of Metis for the partitioner", "DMPlexRebalanceSharedPoints", &parallel);
1672:   PetscOptionsBool("-dm_plex_rebalance_shared_points_use_initial_guess", "Use current partition to bootstrap ParMetis partition", "DMPlexRebalanceSharedPoints", useInitialGuess, &useInitialGuess, NULL);
1673:   PetscOptionsBool("-dm_plex_rebalance_shared_points_use_mat_partitioning", "Use the MatPartitioning object to partition", "DMPlexRebalanceSharedPoints", usematpartitioning, &usematpartitioning, NULL);
1674:   PetscOptionsViewer("-dm_plex_rebalance_shared_points_monitor", "Monitor the shared points rebalance process", "DMPlexRebalanceSharedPoints", &viewer, &format, NULL);
1675:   PetscOptionsEnd();
1676:   if (viewer) PetscViewerPushFormat(viewer, format);

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

1680:   DMGetOptionsPrefix(dm, &prefix);
1681:   PetscOptionsGetViewer(comm, ((PetscObject)dm)->options, prefix, "-dm_rebalance_partition_view", &viewer, &format, NULL);
1682:   if (viewer) PetscViewerPushFormat(viewer, format);

1684:   /* Figure out all points in the plex that we are interested in balancing. */
1685:   DMPlexGetDepthStratum(dm, entityDepth, &eBegin, &eEnd);
1686:   DMPlexGetChart(dm, &pStart, &pEnd);
1687:   PetscMalloc1(pEnd - pStart, &toBalance);
1688:   for (i = 0; i < pEnd - pStart; i++) toBalance[i] = (PetscBool)(i >= eBegin && i < eEnd);

1690:   /* There are three types of points:
1691:    * exclusivelyOwned: points that are owned by this process and only seen by this process
1692:    * nonExclusivelyOwned: points that are owned by this process but seen by at least another process
1693:    * leaf: a point that is seen by this process but owned by a different process
1694:    */
1695:   DMGetPointSF(dm, &sf);
1696:   PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
1697:   PetscMalloc1(pEnd - pStart, &isLeaf);
1698:   PetscMalloc1(pEnd - pStart, &isNonExclusivelyOwned);
1699:   PetscMalloc1(pEnd - pStart, &isExclusivelyOwned);
1700:   for (i = 0; i < pEnd - pStart; i++) {
1701:     isNonExclusivelyOwned[i] = PETSC_FALSE;
1702:     isExclusivelyOwned[i]    = PETSC_FALSE;
1703:     isLeaf[i]                = PETSC_FALSE;
1704:   }

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

1709:   /* for an owned point, we can figure out whether another processor sees it or
1710:    * not by calculating its degree */
1711:   PetscSFComputeDegreeBegin(sf, &degrees);
1712:   PetscSFComputeDegreeEnd(sf, &degrees);
1713:   numExclusivelyOwned    = 0;
1714:   numNonExclusivelyOwned = 0;
1715:   for (i = 0; i < pEnd - pStart; i++) {
1716:     if (toBalance[i]) {
1717:       if (degrees[i] > 0) {
1718:         isNonExclusivelyOwned[i] = PETSC_TRUE;
1719:         numNonExclusivelyOwned += 1;
1720:       } else {
1721:         if (!isLeaf[i]) {
1722:           isExclusivelyOwned[i] = PETSC_TRUE;
1723:           numExclusivelyOwned += 1;
1724:         }
1725:       }
1726:     }
1727:   }

1729:   /* Build a graph with one vertex per core representing the
1730:    * exclusively owned points and then one vertex per nonExclusively owned
1731:    * point. */
1732:   PetscLayoutCreate(comm, &layout);
1733:   PetscLayoutSetLocalSize(layout, 1 + numNonExclusivelyOwned);
1734:   PetscLayoutSetUp(layout);
1735:   PetscLayoutGetRanges(layout, &cumSumVertices);
1736:   PetscMalloc1(pEnd - pStart, &globalNumbersOfLocalOwnedVertices);
1737:   for (i = 0; i < pEnd - pStart; i++) globalNumbersOfLocalOwnedVertices[i] = pStart - 1;
1738:   offset  = cumSumVertices[rank];
1739:   counter = 0;
1740:   for (i = 0; i < pEnd - pStart; i++) {
1741:     if (toBalance[i]) {
1742:       if (degrees[i] > 0) {
1743:         globalNumbersOfLocalOwnedVertices[i] = counter + 1 + offset;
1744:         counter++;
1745:       }
1746:     }
1747:   }

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

1754:   /* Build the graph for partitioning */
1755:   numRows = 1 + numNonExclusivelyOwned;
1756:   PetscLogEventBegin(DMPLEX_RebalBuildGraph, dm, 0, 0, 0);
1757:   MatCreate(comm, &A);
1758:   MatSetType(A, MATMPIADJ);
1759:   MatSetSizes(A, numRows, numRows, cumSumVertices[size], cumSumVertices[size]);
1760:   idx = cumSumVertices[rank];
1761:   for (i = 0; i < pEnd - pStart; i++) {
1762:     if (toBalance[i]) {
1763:       if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
1764:       else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
1765:       else continue;
1766:       MatSetValue(A, idx, jdx, 1, INSERT_VALUES);
1767:       MatSetValue(A, jdx, idx, 1, INSERT_VALUES);
1768:     }
1769:   }
1770:   PetscFree(globalNumbersOfLocalOwnedVertices);
1771:   PetscFree(leafGlobalNumbers);
1772:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1773:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1774:   PetscLogEventEnd(DMPLEX_RebalBuildGraph, dm, 0, 0, 0);

1776:   nparts = size;
1777:   ncon   = 1;
1778:   PetscMalloc1(ncon * nparts, &tpwgts);
1779:   for (i = 0; i < ncon * nparts; i++) tpwgts[i] = 1. / (nparts);
1780:   PetscMalloc1(ncon, &ubvec);
1781:   ubvec[0] = 1.05;
1782:   ubvec[1] = 1.05;

1784:   PetscMalloc1(ncon * (1 + numNonExclusivelyOwned), &vtxwgt);
1785:   if (ncon == 2) {
1786:     vtxwgt[0] = numExclusivelyOwned;
1787:     vtxwgt[1] = 1;
1788:     for (i = 0; i < numNonExclusivelyOwned; i++) {
1789:       vtxwgt[ncon * (i + 1)]     = 1;
1790:       vtxwgt[ncon * (i + 1) + 1] = 0;
1791:     }
1792:   } else {
1793:     PetscInt base, ms;
1794:     MPI_Allreduce(&numExclusivelyOwned, &base, 1, MPIU_INT, MPIU_MAX, PetscObjectComm((PetscObject)dm));
1795:     MatGetSize(A, &ms, NULL);
1796:     ms -= size;
1797:     base      = PetscMax(base, ms);
1798:     vtxwgt[0] = base + numExclusivelyOwned;
1799:     for (i = 0; i < numNonExclusivelyOwned; i++) vtxwgt[i + 1] = 1;
1800:   }

1802:   if (viewer) {
1803:     PetscViewerASCIIPrintf(viewer, "Attempt rebalancing of shared points of depth %" PetscInt_FMT " on interface of mesh distribution.\n", entityDepth);
1804:     PetscViewerASCIIPrintf(viewer, "Size of generated auxiliary graph: %" PetscInt_FMT "\n", cumSumVertices[size]);
1805:   }
1806:   /* TODO: Drop the parallel/sequential choice here and just use MatPartioner for much more flexibility */
1807:   if (usematpartitioning) {
1808:     const char *prefix;

1810:     MatPartitioningCreate(PetscObjectComm((PetscObject)dm), &mp);
1811:     PetscObjectSetOptionsPrefix((PetscObject)mp, "dm_plex_rebalance_shared_points_");
1812:     PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
1813:     PetscObjectPrependOptionsPrefix((PetscObject)mp, prefix);
1814:     MatPartitioningSetAdjacency(mp, A);
1815:     MatPartitioningSetNumberVertexWeights(mp, ncon);
1816:     MatPartitioningSetVertexWeights(mp, vtxwgt);
1817:     MatPartitioningSetFromOptions(mp);
1818:     MatPartitioningApply(mp, &ispart);
1819:     ISGetIndices(ispart, (const PetscInt **)&part);
1820:   } else if (parallel) {
1821:     if (viewer) PetscViewerASCIIPrintf(viewer, "Using ParMETIS to partition graph.\n");
1822:     PetscMalloc1(cumSumVertices[rank + 1] - cumSumVertices[rank], &part);
1823:     MatGetRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, &xadj, &adjncy, &done);
1825:     PetscMalloc1(4, &options);
1826:     options[0] = 1;
1827:     options[1] = 0; /* Verbosity */
1828:     if (viewer) options[1] = 1;
1829:     options[2] = 0;                    /* Seed */
1830:     options[3] = PARMETIS_PSR_COUPLED; /* Seed */
1831:     wgtflag    = 2;
1832:     numflag    = 0;
1833:     if (useInitialGuess) {
1834:       PetscViewerASCIIPrintf(viewer, "THIS DOES NOT WORK! I don't know why. Using current distribution of points as initial guess.\n");
1835:       for (i = 0; i < numRows; i++) part[i] = rank;
1836:       if (viewer) PetscViewerASCIIPrintf(viewer, "Using current distribution of points as initial guess.\n");
1837:       PetscStackPushExternal("ParMETIS_V3_RefineKway");
1838:       PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0);
1839:       ParMETIS_V3_RefineKway((PetscInt *)cumSumVertices, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt, NULL, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
1840:       PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0);
1841:       PetscStackPop;
1843:     } else {
1844:       PetscStackPushExternal("ParMETIS_V3_PartKway");
1845:       PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0);
1846:       ParMETIS_V3_PartKway((PetscInt *)cumSumVertices, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt, NULL, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
1847:       PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0);
1848:       PetscStackPop;
1850:     }
1851:     MatRestoreRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, &xadj, &adjncy, &done);
1852:     PetscFree(options);
1853:   } else {
1854:     if (viewer) PetscViewerASCIIPrintf(viewer, "Using METIS to partition graph.\n");
1855:     Mat       As;
1856:     PetscInt *partGlobal;
1857:     PetscInt *numExclusivelyOwnedAll;

1859:     PetscMalloc1(cumSumVertices[rank + 1] - cumSumVertices[rank], &part);
1860:     MatGetSize(A, &numRows, NULL);
1861:     PetscLogEventBegin(DMPLEX_RebalGatherGraph, dm, 0, 0, 0);
1862:     MatMPIAdjToSeqRankZero(A, &As);
1863:     PetscLogEventEnd(DMPLEX_RebalGatherGraph, dm, 0, 0, 0);

1865:     PetscMalloc1(size, &numExclusivelyOwnedAll);
1866:     numExclusivelyOwnedAll[rank] = numExclusivelyOwned;
1867:     MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, numExclusivelyOwnedAll, 1, MPIU_INT, comm);

1869:     PetscMalloc1(numRows, &partGlobal);
1870:     PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0);
1871:     if (rank == 0) {
1872:       PetscInt *vtxwgt_g, numRows_g;

1874:       MatGetRowIJ(As, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows_g, &xadj, &adjncy, &done);
1875:       PetscMalloc1(2 * numRows_g, &vtxwgt_g);
1876:       for (i = 0; i < size; i++) {
1877:         vtxwgt_g[ncon * cumSumVertices[i]] = numExclusivelyOwnedAll[i];
1878:         if (ncon > 1) vtxwgt_g[ncon * cumSumVertices[i] + 1] = 1;
1879:         for (j = cumSumVertices[i] + 1; j < cumSumVertices[i + 1]; j++) {
1880:           vtxwgt_g[ncon * j] = 1;
1881:           if (ncon > 1) vtxwgt_g[2 * j + 1] = 0;
1882:         }
1883:       }

1885:       PetscMalloc1(64, &options);
1886:       METIS_SetDefaultOptions(options); /* initialize all defaults */
1888:       options[METIS_OPTION_CONTIG] = 1;
1889:       PetscStackPushExternal("METIS_PartGraphKway");
1890:       METIS_PartGraphKway(&numRows_g, &ncon, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt_g, NULL, NULL, &nparts, tpwgts, ubvec, options, &edgecut, partGlobal);
1891:       PetscStackPop;
1893:       PetscFree(options);
1894:       PetscFree(vtxwgt_g);
1895:       MatRestoreRowIJ(As, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows_g, &xadj, &adjncy, &done);
1896:       MatDestroy(&As);
1897:     }
1898:     PetscBarrier((PetscObject)dm);
1899:     PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0);
1900:     PetscFree(numExclusivelyOwnedAll);

1902:     /* scatter the partitioning information to ranks */
1903:     PetscLogEventBegin(DMPLEX_RebalScatterPart, 0, 0, 0, 0);
1904:     PetscMalloc1(size, &counts);
1905:     PetscMalloc1(size + 1, &mpiCumSumVertices);
1906:     for (i = 0; i < size; i++) PetscMPIIntCast(cumSumVertices[i + 1] - cumSumVertices[i], &(counts[i]));
1907:     for (i = 0; i <= size; i++) PetscMPIIntCast(cumSumVertices[i], &(mpiCumSumVertices[i]));
1908:     MPI_Scatterv(partGlobal, counts, mpiCumSumVertices, MPIU_INT, part, counts[rank], MPIU_INT, 0, comm);
1909:     PetscFree(counts);
1910:     PetscFree(mpiCumSumVertices);
1911:     PetscFree(partGlobal);
1912:     PetscLogEventEnd(DMPLEX_RebalScatterPart, 0, 0, 0, 0);
1913:   }
1914:   PetscFree(ubvec);
1915:   PetscFree(tpwgts);

1917:   /* Rename the result so that the vertex resembling the exclusively owned points stays on the same rank */
1918:   PetscMalloc2(size, &firstVertices, size, &renumbering);
1919:   firstVertices[rank] = part[0];
1920:   MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, firstVertices, 1, MPIU_INT, comm);
1921:   for (i = 0; i < size; i++) renumbering[firstVertices[i]] = i;
1922:   for (i = 0; i < cumSumVertices[rank + 1] - cumSumVertices[rank]; i++) part[i] = renumbering[part[i]];
1923:   PetscFree2(firstVertices, renumbering);

1925:   /* Check if the renumbering worked (this can fail when ParMETIS gives fewer partitions than there are processes) */
1926:   failed = (PetscInt)(part[0] != rank);
1927:   MPI_Allreduce(&failed, &failedGlobal, 1, MPIU_INT, MPI_SUM, comm);
1928:   if (failedGlobal > 0) {
1930:     PetscFree(vtxwgt);
1931:     PetscFree(toBalance);
1932:     PetscFree(isLeaf);
1933:     PetscFree(isNonExclusivelyOwned);
1934:     PetscFree(isExclusivelyOwned);
1935:     if (usematpartitioning) {
1936:       ISRestoreIndices(ispart, (const PetscInt **)&part);
1937:       ISDestroy(&ispart);
1938:     } else PetscFree(part);
1939:     if (viewer) {
1940:       PetscViewerPopFormat(viewer);
1941:       PetscViewerDestroy(&viewer);
1942:     }
1943:     PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
1944:     return 0;
1945:   }

1947:   /* Check how well we did distributing points*/
1948:   if (viewer) {
1949:     PetscViewerASCIIPrintf(viewer, "Number of owned entities of depth %" PetscInt_FMT ".\n", entityDepth);
1950:     PetscViewerASCIIPrintf(viewer, "Initial      ");
1951:     DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, NULL, viewer);
1952:     PetscViewerASCIIPrintf(viewer, "Rebalanced   ");
1953:     DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, part, viewer);
1954:   }

1956:   /* Check that every vertex is owned by a process that it is actually connected to. */
1957:   MatGetRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done);
1958:   for (i = 1; i <= numNonExclusivelyOwned; i++) {
1959:     PetscInt loc = 0;
1960:     PetscFindInt(cumSumVertices[part[i]], xadj[i + 1] - xadj[i], &adjncy[xadj[i]], &loc);
1961:     /* If not, then just set the owner to the original owner (hopefully a rare event, it means that a vertex has been isolated) */
1962:     if (loc < 0) part[i] = rank;
1963:   }
1964:   MatRestoreRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done);
1965:   MatDestroy(&A);

1967:   /* See how significant the influences of the previous fixing up step was.*/
1968:   if (viewer) {
1969:     PetscViewerASCIIPrintf(viewer, "After fix    ");
1970:     DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, part, viewer);
1971:   }
1972:   if (!usematpartitioning) PetscFree(vtxwgt);
1973:   else MatPartitioningDestroy(&mp);

1975:   PetscLayoutDestroy(&layout);

1977:   PetscLogEventBegin(DMPLEX_RebalRewriteSF, dm, 0, 0, 0);
1978:   /* Rewrite the SF to reflect the new ownership. */
1979:   PetscMalloc1(numNonExclusivelyOwned, &pointsToRewrite);
1980:   counter = 0;
1981:   for (i = 0; i < pEnd - pStart; i++) {
1982:     if (toBalance[i]) {
1983:       if (isNonExclusivelyOwned[i]) {
1984:         pointsToRewrite[counter] = i + pStart;
1985:         counter++;
1986:       }
1987:     }
1988:   }
1989:   DMPlexRewriteSF(dm, numNonExclusivelyOwned, pointsToRewrite, part + 1, degrees);
1990:   PetscFree(pointsToRewrite);
1991:   PetscLogEventEnd(DMPLEX_RebalRewriteSF, dm, 0, 0, 0);

1993:   PetscFree(toBalance);
1994:   PetscFree(isLeaf);
1995:   PetscFree(isNonExclusivelyOwned);
1996:   PetscFree(isExclusivelyOwned);
1997:   if (usematpartitioning) {
1998:     ISRestoreIndices(ispart, (const PetscInt **)&part);
1999:     ISDestroy(&ispart);
2000:   } else PetscFree(part);
2001:   if (viewer) {
2002:     PetscViewerPopFormat(viewer);
2003:     PetscViewerDestroy(&viewer);
2004:   }
2005:   if (success) *success = PETSC_TRUE;
2006:   PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
2007:   return 0;
2008: #else
2009:   SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
2010: #endif
2011: }