Actual source code: plexpartition.c

petsc-master 2020-09-26
Report Typos and Errors
  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);
103:     PetscSFBcastEnd(dm->sf, MPIU_INT, adjCells, remoteCells);
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);
1455:   PetscSFBcastEnd(sf, MPIU_INT, newOwners, newOwners);
1456:   PetscSFBcastBegin(sf, MPIU_INT, newNumbers, newNumbers);
1457:   PetscSFBcastEnd(sf, MPIU_INT, newNumbers, newNumbers);

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);
1671:   PetscSFBcastEnd(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers);

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: #else
1966:   SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
1967: #endif
1968:   return(0);
1969: }