Actual source code: plexdistribute.c

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

  5: /*@C
  6:   DMPlexSetAdjacencyUser - Define adjacency in the mesh using a user-provided callback

  8:   Input Parameters:
  9: + dm   - The DM object
 10: . user - The user callback, may be `NULL` (to clear the callback)
 11: - ctx  - context for callback evaluation, may be `NULL`

 13:   Level: advanced

 15:   Notes:
 16:   The caller of `DMPlexGetAdjacency()` may need to arrange that a large enough array is available for the adjacency.

 18:   Any setting here overrides other configuration of `DMPLEX` adjacency determination.

 20: .seealso: `DMPLEX`, `DMSetAdjacency()`, `DMPlexDistribute()`, `DMPlexPreallocateOperator()`, `DMPlexGetAdjacency()`, `DMPlexGetAdjacencyUser()`
 21: @*/
 22: PetscErrorCode DMPlexSetAdjacencyUser(DM dm, PetscErrorCode (*user)(DM, PetscInt, PetscInt *, PetscInt[], void *), void *ctx)
 23: {
 24:   DM_Plex *mesh = (DM_Plex *)dm->data;

 26:   PetscFunctionBegin;
 28:   mesh->useradjacency    = user;
 29:   mesh->useradjacencyctx = ctx;
 30:   PetscFunctionReturn(PETSC_SUCCESS);
 31: }

 33: /*@C
 34:   DMPlexGetAdjacencyUser - get the user-defined adjacency callback

 36:   Input Parameter:
 37: . dm - The `DM` object

 39:   Output Parameters:
 40: + user - The callback
 41: - ctx  - context for callback evaluation

 43:   Level: advanced

 45: .seealso: `DMPLEX`, `DMSetAdjacency()`, `DMPlexDistribute()`, `DMPlexPreallocateOperator()`, `DMPlexGetAdjacency()`, `DMPlexSetAdjacencyUser()`
 46: @*/
 47: PetscErrorCode DMPlexGetAdjacencyUser(DM dm, PetscErrorCode (**user)(DM, PetscInt, PetscInt *, PetscInt[], void *), void **ctx)
 48: {
 49:   DM_Plex *mesh = (DM_Plex *)dm->data;

 51:   PetscFunctionBegin;
 53:   if (user) *user = mesh->useradjacency;
 54:   if (ctx) *ctx = mesh->useradjacencyctx;
 55:   PetscFunctionReturn(PETSC_SUCCESS);
 56: }

 58: /*@
 59:   DMPlexSetAdjacencyUseAnchors - Define adjacency in the mesh using the point-to-point constraints.

 61:   Input Parameters:
 62: + dm         - The `DM` object
 63: - useAnchors - Flag to use the constraints.  If PETSC_TRUE, then constrained points are omitted from DMPlexGetAdjacency(), and their anchor points appear in their place.

 65:   Level: intermediate

 67: .seealso: `DMPLEX`, `DMGetAdjacency()`, `DMSetAdjacency()`, `DMPlexDistribute()`, `DMPlexPreallocateOperator()`, `DMPlexSetAnchors()`
 68: @*/
 69: PetscErrorCode DMPlexSetAdjacencyUseAnchors(DM dm, PetscBool useAnchors)
 70: {
 71:   DM_Plex *mesh = (DM_Plex *)dm->data;

 73:   PetscFunctionBegin;
 75:   mesh->useAnchors = useAnchors;
 76:   PetscFunctionReturn(PETSC_SUCCESS);
 77: }

 79: /*@
 80:   DMPlexGetAdjacencyUseAnchors - Query whether adjacency in the mesh uses the point-to-point constraints.

 82:   Input Parameter:
 83: . dm - The `DM` object

 85:   Output Parameter:
 86: . useAnchors - Flag to use the closure.  If PETSC_TRUE, then constrained points are omitted from DMPlexGetAdjacency(), and their anchor points appear in their place.

 88:   Level: intermediate

 90: .seealso: `DMPLEX`, `DMPlexSetAdjacencyUseAnchors()`, `DMSetAdjacency()`, `DMGetAdjacency()`, `DMPlexDistribute()`, `DMPlexPreallocateOperator()`, `DMPlexSetAnchors()`
 91: @*/
 92: PetscErrorCode DMPlexGetAdjacencyUseAnchors(DM dm, PetscBool *useAnchors)
 93: {
 94:   DM_Plex *mesh = (DM_Plex *)dm->data;

 96:   PetscFunctionBegin;
 98:   PetscAssertPointer(useAnchors, 2);
 99:   *useAnchors = mesh->useAnchors;
100:   PetscFunctionReturn(PETSC_SUCCESS);
101: }

103: static PetscErrorCode DMPlexGetAdjacency_Cone_Internal(DM dm, PetscInt p, PetscInt *adjSize, PetscInt adj[])
104: {
105:   const PetscInt *cone   = NULL;
106:   PetscInt        numAdj = 0, maxAdjSize = *adjSize, coneSize, c;

108:   PetscFunctionBeginHot;
109:   PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
110:   PetscCall(DMPlexGetCone(dm, p, &cone));
111:   for (c = 0; c <= coneSize; ++c) {
112:     const PetscInt  point   = !c ? p : cone[c - 1];
113:     const PetscInt *support = NULL;
114:     PetscInt        supportSize, s, q;

116:     PetscCall(DMPlexGetSupportSize(dm, point, &supportSize));
117:     PetscCall(DMPlexGetSupport(dm, point, &support));
118:     for (s = 0; s < supportSize; ++s) {
119:       for (q = 0; q < numAdj || ((void)(adj[numAdj++] = support[s]), 0); ++q) {
120:         if (support[s] == adj[q]) break;
121:       }
122:       PetscCheck(numAdj <= maxAdjSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%" PetscInt_FMT ")", maxAdjSize);
123:     }
124:   }
125:   *adjSize = numAdj;
126:   PetscFunctionReturn(PETSC_SUCCESS);
127: }

129: static PetscErrorCode DMPlexGetAdjacency_Support_Internal(DM dm, PetscInt p, PetscInt *adjSize, PetscInt adj[])
130: {
131:   const PetscInt *support = NULL;
132:   PetscInt        numAdj = 0, maxAdjSize = *adjSize, supportSize, s;

134:   PetscFunctionBeginHot;
135:   PetscCall(DMPlexGetSupportSize(dm, p, &supportSize));
136:   PetscCall(DMPlexGetSupport(dm, p, &support));
137:   for (s = 0; s <= supportSize; ++s) {
138:     const PetscInt  point = !s ? p : support[s - 1];
139:     const PetscInt *cone  = NULL;
140:     PetscInt        coneSize, c, q;

142:     PetscCall(DMPlexGetConeSize(dm, point, &coneSize));
143:     PetscCall(DMPlexGetCone(dm, point, &cone));
144:     for (c = 0; c < coneSize; ++c) {
145:       for (q = 0; q < numAdj || ((void)(adj[numAdj++] = cone[c]), 0); ++q) {
146:         if (cone[c] == adj[q]) break;
147:       }
148:       PetscCheck(numAdj <= maxAdjSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%" PetscInt_FMT ")", maxAdjSize);
149:     }
150:   }
151:   *adjSize = numAdj;
152:   PetscFunctionReturn(PETSC_SUCCESS);
153: }

155: static PetscErrorCode DMPlexGetAdjacency_Transitive_Internal(DM dm, PetscInt p, PetscBool useClosure, PetscInt *adjSize, PetscInt adj[])
156: {
157:   PetscInt *star   = NULL;
158:   PetscInt  numAdj = 0, maxAdjSize = *adjSize, starSize, s;

160:   PetscFunctionBeginHot;
161:   PetscCall(DMPlexGetTransitiveClosure(dm, p, useClosure, &starSize, &star));
162:   for (s = 0; s < starSize * 2; s += 2) {
163:     const PetscInt *closure = NULL;
164:     PetscInt        closureSize, c, q;

166:     PetscCall(DMPlexGetTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt **)&closure));
167:     for (c = 0; c < closureSize * 2; c += 2) {
168:       for (q = 0; q < numAdj || ((void)(adj[numAdj++] = closure[c]), 0); ++q) {
169:         if (closure[c] == adj[q]) break;
170:       }
171:       PetscCheck(numAdj <= maxAdjSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%" PetscInt_FMT ")", maxAdjSize);
172:     }
173:     PetscCall(DMPlexRestoreTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt **)&closure));
174:   }
175:   PetscCall(DMPlexRestoreTransitiveClosure(dm, p, useClosure, &starSize, &star));
176:   *adjSize = numAdj;
177:   PetscFunctionReturn(PETSC_SUCCESS);
178: }

180: // Returns the maximum number of adjacent points in the DMPlex
181: PetscErrorCode DMPlexGetMaxAdjacencySize_Internal(DM dm, PetscBool useAnchors, PetscInt *max_adjacency_size)
182: {
183:   PetscInt depth, maxC, maxS, maxP, pStart, pEnd, asiz, maxAnchors = 1;

185:   PetscFunctionBeginUser;
186:   if (useAnchors) {
187:     PetscSection aSec = NULL;
188:     IS           aIS  = NULL;
189:     PetscInt     aStart, aEnd;
190:     PetscCall(DMPlexGetAnchors(dm, &aSec, &aIS));
191:     if (aSec) {
192:       PetscCall(PetscSectionGetMaxDof(aSec, &maxAnchors));
193:       maxAnchors = PetscMax(1, maxAnchors);
194:       PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));
195:     }
196:   }

198:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
199:   PetscCall(DMPlexGetDepth(dm, &depth));
200:   depth = PetscMax(depth, -depth);
201:   PetscCall(DMPlexGetMaxSizes(dm, &maxC, &maxS));
202:   maxP = maxS * maxC;
203:   /* Adjacency can be as large as supp(cl(cell)) or cl(supp(vertex)),
204:           supp(cell) + supp(maxC faces) + supp(maxC^2 edges) + supp(maxC^3 vertices)
205:         = 0 + maxS*maxC + maxS^2*maxC^2 + maxS^3*maxC^3
206:         = \sum^d_{i=0} (maxS*maxC)^i - 1
207:         = (maxS*maxC)^{d+1} - 1 / (maxS*maxC - 1) - 1
208:     We could improve this by getting the max by strata:
209:           supp[d](cell) + supp[d-1](maxC[d] faces) + supp[1](maxC[d]*maxC[d-1] edges) + supp[0](maxC[d]*maxC[d-1]*maxC[d-2] vertices)
210:         = 0 + maxS[d-1]*maxC[d] + maxS[1]*maxC[d]*maxC[d-1] + maxS[0]*maxC[d]*maxC[d-1]*maxC[d-2]
211:     and the same with S and C reversed
212:   */
213:   if ((depth == 3 && maxP > 200) || (depth == 2 && maxP > 580)) asiz = pEnd - pStart;
214:   else asiz = (maxP > 1) ? ((PetscPowInt(maxP, depth + 1) - 1) / (maxP - 1)) : depth + 1;
215:   asiz *= maxAnchors;
216:   *max_adjacency_size = PetscMin(asiz, pEnd - pStart);
217:   PetscFunctionReturn(PETSC_SUCCESS);
218: }

220: // Returns Adjacent mesh points to the selected point given specific criteria
221: //
222: // + adjSize - Number of adjacent points
223: // - adj - Array of the adjacent points
224: PetscErrorCode DMPlexGetAdjacency_Internal(DM dm, PetscInt p, PetscBool useCone, PetscBool useTransitiveClosure, PetscBool useAnchors, PetscInt *adjSize, PetscInt *adj[])
225: {
226:   static PetscInt asiz   = 0;
227:   PetscInt        aStart = -1, aEnd = -1;
228:   PetscInt        maxAdjSize;
229:   PetscSection    aSec = NULL;
230:   IS              aIS  = NULL;
231:   const PetscInt *anchors;
232:   DM_Plex        *mesh = (DM_Plex *)dm->data;

234:   PetscFunctionBeginHot;
235:   if (useAnchors) {
236:     PetscCall(DMPlexGetAnchors(dm, &aSec, &aIS));
237:     if (aSec) {
238:       PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));
239:       PetscCall(ISGetIndices(aIS, &anchors));
240:     }
241:   }
242:   if (!*adj) {
243:     PetscCall(DMPlexGetMaxAdjacencySize_Internal(dm, useAnchors, &asiz));
244:     PetscCall(PetscMalloc1(asiz, adj));
245:   }
246:   if (*adjSize < 0) *adjSize = asiz;
247:   maxAdjSize = *adjSize;
248:   if (mesh->useradjacency) {
249:     PetscCall((*mesh->useradjacency)(dm, p, adjSize, *adj, mesh->useradjacencyctx));
250:   } else if (useTransitiveClosure) {
251:     PetscCall(DMPlexGetAdjacency_Transitive_Internal(dm, p, useCone, adjSize, *adj));
252:   } else if (useCone) {
253:     PetscCall(DMPlexGetAdjacency_Cone_Internal(dm, p, adjSize, *adj));
254:   } else {
255:     PetscCall(DMPlexGetAdjacency_Support_Internal(dm, p, adjSize, *adj));
256:   }
257:   if (useAnchors && aSec) {
258:     PetscInt  origSize = *adjSize;
259:     PetscInt  numAdj   = origSize;
260:     PetscInt  i        = 0, j;
261:     PetscInt *orig     = *adj;

263:     while (i < origSize) {
264:       PetscInt p    = orig[i];
265:       PetscInt aDof = 0;

267:       if (p >= aStart && p < aEnd) PetscCall(PetscSectionGetDof(aSec, p, &aDof));
268:       if (aDof) {
269:         PetscInt aOff;
270:         PetscInt s, q;

272:         for (j = i + 1; j < numAdj; j++) orig[j - 1] = orig[j];
273:         origSize--;
274:         numAdj--;
275:         PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
276:         for (s = 0; s < aDof; ++s) {
277:           for (q = 0; q < numAdj || ((void)(orig[numAdj++] = anchors[aOff + s]), 0); ++q) {
278:             if (anchors[aOff + s] == orig[q]) break;
279:           }
280:           PetscCheck(numAdj <= maxAdjSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%" PetscInt_FMT ")", maxAdjSize);
281:         }
282:       } else {
283:         i++;
284:       }
285:     }
286:     *adjSize = numAdj;
287:     PetscCall(ISRestoreIndices(aIS, &anchors));
288:   }
289:   PetscFunctionReturn(PETSC_SUCCESS);
290: }

292: /*@
293:   DMPlexGetAdjacency - Return all points adjacent to the given point

295:   Input Parameters:
296: + dm - The `DM` object
297: - p  - The point

299:   Input/Output Parameters:
300: + adjSize - The maximum size of `adj` if it is non-`NULL`, or `PETSC_DETERMINE`;
301:             on output the number of adjacent points
302: - adj     - Either `NULL` so that the array is allocated, or an existing array with size `adjSize`;
303:             on output contains the adjacent points

305:   Level: advanced

307:   Notes:
308:   The user must `PetscFree()` the `adj` array if it was not passed in.

310: .seealso: `DMPLEX`, `DMSetAdjacency()`, `DMPlexDistribute()`, `DMCreateMatrix()`, `DMPlexPreallocateOperator()`
311: @*/
312: PetscErrorCode DMPlexGetAdjacency(DM dm, PetscInt p, PetscInt *adjSize, PetscInt *adj[])
313: {
314:   PetscBool useCone, useClosure, useAnchors;

316:   PetscFunctionBeginHot;
318:   PetscAssertPointer(adjSize, 3);
319:   PetscAssertPointer(adj, 4);
320:   PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure));
321:   PetscCall(DMPlexGetAdjacencyUseAnchors(dm, &useAnchors));
322:   PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, adjSize, adj));
323:   PetscFunctionReturn(PETSC_SUCCESS);
324: }

326: /*@
327:   DMPlexCreateTwoSidedProcessSF - Create an `PetscSF` which just has process connectivity

329:   Collective

331:   Input Parameters:
332: + dm              - The `DM`
333: . sfPoint         - The `PetscSF` which encodes point connectivity
334: . rootRankSection - to be documented
335: . rootRanks       - to be documented
336: . leafRankSection - to be documented
337: - leafRanks       - to be documented

339:   Output Parameters:
340: + processRanks - A list of process neighbors, or `NULL`
341: - sfProcess    - An `PetscSF` encoding the two-sided process connectivity, or `NULL`

343:   Level: developer

345: .seealso: `DMPLEX`, `PetscSFCreate()`, `DMPlexCreateProcessSF()`
346: @*/
347: PetscErrorCode DMPlexCreateTwoSidedProcessSF(DM dm, PetscSF sfPoint, PetscSection rootRankSection, IS rootRanks, PetscSection leafRankSection, IS leafRanks, PeOp IS *processRanks, PeOp PetscSF *sfProcess)
348: {
349:   const PetscSFNode *remotePoints;
350:   PetscInt          *localPointsNew;
351:   PetscSFNode       *remotePointsNew;
352:   const PetscInt    *nranks;
353:   PetscInt          *ranksNew;
354:   PetscBT            neighbors;
355:   PetscInt           pStart, pEnd, p, numLeaves, l, numNeighbors, n;
356:   PetscMPIInt        size, proc, rank;

358:   PetscFunctionBegin;
361:   if (processRanks) PetscAssertPointer(processRanks, 7);
362:   if (sfProcess) PetscAssertPointer(sfProcess, 8);
363:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
364:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
365:   PetscCall(PetscSFGetGraph(sfPoint, NULL, &numLeaves, NULL, &remotePoints));
366:   PetscCall(PetscBTCreate(size, &neighbors));
367:   PetscCall(PetscBTMemzero(size, neighbors));
368:   /* Compute root-to-leaf process connectivity */
369:   PetscCall(PetscSectionGetChart(rootRankSection, &pStart, &pEnd));
370:   PetscCall(ISGetIndices(rootRanks, &nranks));
371:   for (p = pStart; p < pEnd; ++p) {
372:     PetscInt ndof, noff, n;

374:     PetscCall(PetscSectionGetDof(rootRankSection, p, &ndof));
375:     PetscCall(PetscSectionGetOffset(rootRankSection, p, &noff));
376:     for (n = 0; n < ndof; ++n) PetscCall(PetscBTSet(neighbors, nranks[noff + n]));
377:   }
378:   PetscCall(ISRestoreIndices(rootRanks, &nranks));
379:   /* Compute leaf-to-neighbor process connectivity */
380:   PetscCall(PetscSectionGetChart(leafRankSection, &pStart, &pEnd));
381:   PetscCall(ISGetIndices(leafRanks, &nranks));
382:   for (p = pStart; p < pEnd; ++p) {
383:     PetscInt ndof, noff, n;

385:     PetscCall(PetscSectionGetDof(leafRankSection, p, &ndof));
386:     PetscCall(PetscSectionGetOffset(leafRankSection, p, &noff));
387:     for (n = 0; n < ndof; ++n) PetscCall(PetscBTSet(neighbors, nranks[noff + n]));
388:   }
389:   PetscCall(ISRestoreIndices(leafRanks, &nranks));
390:   /* Compute leaf-to-root process connectivity */
391:   for (l = 0; l < numLeaves; ++l) PetscCall(PetscBTSet(neighbors, remotePoints[l].rank));
392:   /* Calculate edges */
393:   PetscCall(PetscBTClear(neighbors, rank));
394:   for (proc = 0, numNeighbors = 0; proc < size; ++proc) {
395:     if (PetscBTLookup(neighbors, proc)) ++numNeighbors;
396:   }
397:   PetscCall(PetscMalloc1(numNeighbors, &ranksNew));
398:   PetscCall(PetscMalloc1(numNeighbors, &localPointsNew));
399:   PetscCall(PetscMalloc1(numNeighbors, &remotePointsNew));
400:   for (proc = 0, n = 0; proc < size; ++proc) {
401:     if (PetscBTLookup(neighbors, proc)) {
402:       ranksNew[n]              = proc;
403:       localPointsNew[n]        = proc;
404:       remotePointsNew[n].index = rank;
405:       remotePointsNew[n].rank  = proc;
406:       ++n;
407:     }
408:   }
409:   PetscCall(PetscBTDestroy(&neighbors));
410:   if (processRanks) PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), numNeighbors, ranksNew, PETSC_OWN_POINTER, processRanks));
411:   else PetscCall(PetscFree(ranksNew));
412:   if (sfProcess) {
413:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess));
414:     PetscCall(PetscObjectSetName((PetscObject)*sfProcess, "Two-Sided Process SF"));
415:     PetscCall(PetscSFSetFromOptions(*sfProcess));
416:     PetscCall(PetscSFSetGraph(*sfProcess, size, numNeighbors, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
417:   }
418:   PetscFunctionReturn(PETSC_SUCCESS);
419: }

421: /*@
422:   DMPlexDistributeOwnership - Compute owner information for shared points. This basically gets two-sided for an SF.

424:   Collective

426:   Input Parameter:
427: . dm - The `DM`

429:   Output Parameters:
430: + rootSection - The number of leaves for a given root point
431: . rootrank    - The rank of each edge into the root point
432: . leafSection - The number of processes sharing a given leaf point
433: - leafrank    - The rank of each process sharing a leaf point

435:   Level: developer

437: .seealso: `DMPLEX`, `DMPlexCreateOverlapLabel()`, `DMPlexDistribute()`, `DMPlexDistributeOverlap()`
438: @*/
439: PetscErrorCode DMPlexDistributeOwnership(DM dm, PetscSection rootSection, IS *rootrank, PetscSection leafSection, IS *leafrank)
440: {
441:   MPI_Comm        comm;
442:   PetscSF         sfPoint;
443:   const PetscInt *rootdegree;
444:   PetscInt       *myrank, *remoterank;
445:   PetscInt        pStart, pEnd, p, nedges;
446:   PetscMPIInt     rank;

448:   PetscFunctionBegin;
449:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
450:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
451:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
452:   PetscCall(DMGetPointSF(dm, &sfPoint));
453:   /* Compute number of leaves for each root */
454:   PetscCall(PetscObjectSetName((PetscObject)rootSection, "Root Section"));
455:   PetscCall(PetscSectionSetChart(rootSection, pStart, pEnd));
456:   PetscCall(PetscSFComputeDegreeBegin(sfPoint, &rootdegree));
457:   PetscCall(PetscSFComputeDegreeEnd(sfPoint, &rootdegree));
458:   for (p = pStart; p < pEnd; ++p) PetscCall(PetscSectionSetDof(rootSection, p, rootdegree[p - pStart]));
459:   PetscCall(PetscSectionSetUp(rootSection));
460:   /* Gather rank of each leaf to root */
461:   PetscCall(PetscSectionGetStorageSize(rootSection, &nedges));
462:   PetscCall(PetscMalloc1(pEnd - pStart, &myrank));
463:   PetscCall(PetscMalloc1(nedges, &remoterank));
464:   for (p = 0; p < pEnd - pStart; ++p) myrank[p] = rank;
465:   PetscCall(PetscSFGatherBegin(sfPoint, MPIU_INT, myrank, remoterank));
466:   PetscCall(PetscSFGatherEnd(sfPoint, MPIU_INT, myrank, remoterank));
467:   PetscCall(PetscFree(myrank));
468:   PetscCall(ISCreateGeneral(comm, nedges, remoterank, PETSC_OWN_POINTER, rootrank));
469:   /* Distribute remote ranks to leaves */
470:   PetscCall(PetscObjectSetName((PetscObject)leafSection, "Leaf Section"));
471:   PetscCall(DMPlexDistributeFieldIS(dm, sfPoint, rootSection, *rootrank, leafSection, leafrank));
472:   PetscFunctionReturn(PETSC_SUCCESS);
473: }

475: /*@
476:   DMPlexCreateOverlapLabel - Compute a label indicating what overlap points should be sent to new processes

478:   Collective

480:   Input Parameters:
481: + dm          - The `DM`
482: . levels      - Number of overlap levels
483: . rootSection - The number of leaves for a given root point
484: . rootrank    - The rank of each edge into the root point
485: . leafSection - The number of processes sharing a given leaf point
486: - leafrank    - The rank of each process sharing a leaf point

488:   Output Parameter:
489: . ovLabel - `DMLabel` containing remote overlap contributions as point/rank pairings

491:   Level: developer

493: .seealso: `DMPLEX`, `DMPlexCreateOverlapLabelFromLabels()`, `DMPlexGetAdjacency()`, `DMPlexDistributeOwnership()`, `DMPlexDistribute()`
494: @*/
495: PetscErrorCode DMPlexCreateOverlapLabel(DM dm, PetscInt levels, PetscSection rootSection, IS rootrank, PetscSection leafSection, IS leafrank, DMLabel *ovLabel)
496: {
497:   MPI_Comm           comm;
498:   DMLabel            ovAdjByRank; /* A DMLabel containing all points adjacent to shared points, separated by rank (value in label) */
499:   PetscSF            sfPoint;
500:   const PetscSFNode *remote;
501:   const PetscInt    *local;
502:   const PetscInt    *nrank, *rrank;
503:   PetscInt          *adj = NULL;
504:   PetscInt           pStart, pEnd, p, sStart, sEnd, nleaves, l;
505:   PetscMPIInt        rank, size;
506:   PetscBool          flg;

508:   PetscFunctionBegin;
509:   *ovLabel = NULL;
510:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
511:   PetscCallMPI(MPI_Comm_size(comm, &size));
512:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
513:   if (size == 1) PetscFunctionReturn(PETSC_SUCCESS);
514:   PetscCall(DMGetPointSF(dm, &sfPoint));
515:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
516:   if (!levels) {
517:     /* Add owned points */
518:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Overlap label", ovLabel));
519:     for (p = pStart; p < pEnd; ++p) PetscCall(DMLabelSetValue(*ovLabel, p, rank));
520:     PetscFunctionReturn(PETSC_SUCCESS);
521:   }
522:   PetscCall(PetscSectionGetChart(leafSection, &sStart, &sEnd));
523:   PetscCall(PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote));
524:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Overlap adjacency", &ovAdjByRank));
525:   /* Handle leaves: shared with the root point */
526:   for (l = 0; l < nleaves; ++l) {
527:     PetscInt adjSize = PETSC_DETERMINE, a;

529:     PetscCall(DMPlexGetAdjacency(dm, local ? local[l] : l, &adjSize, &adj));
530:     for (a = 0; a < adjSize; ++a) PetscCall(DMLabelSetValue(ovAdjByRank, adj[a], remote[l].rank));
531:   }
532:   PetscCall(ISGetIndices(rootrank, &rrank));
533:   PetscCall(ISGetIndices(leafrank, &nrank));
534:   /* Handle roots */
535:   for (p = pStart; p < pEnd; ++p) {
536:     PetscInt adjSize = PETSC_DETERMINE, neighbors = 0, noff, n, a;

538:     if ((p >= sStart) && (p < sEnd)) {
539:       /* Some leaves share a root with other leaves on different processes */
540:       PetscCall(PetscSectionGetDof(leafSection, p, &neighbors));
541:       if (neighbors) {
542:         PetscCall(PetscSectionGetOffset(leafSection, p, &noff));
543:         PetscCall(DMPlexGetAdjacency(dm, p, &adjSize, &adj));
544:         for (n = 0; n < neighbors; ++n) {
545:           const PetscInt remoteRank = nrank[noff + n];

547:           if (remoteRank == rank) continue;
548:           for (a = 0; a < adjSize; ++a) PetscCall(DMLabelSetValue(ovAdjByRank, adj[a], remoteRank));
549:         }
550:       }
551:     }
552:     /* Roots are shared with leaves */
553:     PetscCall(PetscSectionGetDof(rootSection, p, &neighbors));
554:     if (!neighbors) continue;
555:     PetscCall(PetscSectionGetOffset(rootSection, p, &noff));
556:     PetscCall(DMPlexGetAdjacency(dm, p, &adjSize, &adj));
557:     for (n = 0; n < neighbors; ++n) {
558:       const PetscInt remoteRank = rrank[noff + n];

560:       if (remoteRank == rank) continue;
561:       for (a = 0; a < adjSize; ++a) PetscCall(DMLabelSetValue(ovAdjByRank, adj[a], remoteRank));
562:     }
563:   }
564:   PetscCall(PetscFree(adj));
565:   PetscCall(ISRestoreIndices(rootrank, &rrank));
566:   PetscCall(ISRestoreIndices(leafrank, &nrank));
567:   /* Add additional overlap levels */
568:   for (l = 1; l < levels; l++) {
569:     /* Propagate point donations over SF to capture remote connections */
570:     PetscCall(DMPlexPartitionLabelPropagate(dm, ovAdjByRank));
571:     /* Add next level of point donations to the label */
572:     PetscCall(DMPlexPartitionLabelAdjacency(dm, ovAdjByRank));
573:   }
574:   /* We require the closure in the overlap */
575:   PetscCall(DMPlexPartitionLabelClosure(dm, ovAdjByRank));
576:   PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-overlap_view", &flg));
577:   if (flg) {
578:     PetscViewer viewer;
579:     PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm), &viewer));
580:     PetscCall(DMLabelView(ovAdjByRank, viewer));
581:   }
582:   /* Invert sender to receiver label */
583:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Overlap label", ovLabel));
584:   PetscCall(DMPlexPartitionLabelInvert(dm, ovAdjByRank, NULL, *ovLabel));
585:   /* Add owned points, except for shared local points */
586:   for (p = pStart; p < pEnd; ++p) PetscCall(DMLabelSetValue(*ovLabel, p, rank));
587:   for (l = 0; l < nleaves; ++l) {
588:     PetscCall(DMLabelClearValue(*ovLabel, local[l], rank));
589:     PetscCall(DMLabelSetValue(*ovLabel, remote[l].index, remote[l].rank));
590:   }
591:   /* Clean up */
592:   PetscCall(DMLabelDestroy(&ovAdjByRank));
593:   PetscFunctionReturn(PETSC_SUCCESS);
594: }

596: static PetscErrorCode HandlePoint_Private(DM dm, PetscInt p, PetscSection section, const PetscInt ranks[], PetscInt numExLabels, const DMLabel exLabel[], const PetscInt exValue[], DMLabel ovAdjByRank)
597: {
598:   PetscInt neighbors, el;

600:   PetscFunctionBegin;
601:   PetscCall(PetscSectionGetDof(section, p, &neighbors));
602:   if (neighbors) {
603:     PetscInt   *adj     = NULL;
604:     PetscInt    adjSize = PETSC_DETERMINE, noff, n, a;
605:     PetscMPIInt rank;

607:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
608:     PetscCall(PetscSectionGetOffset(section, p, &noff));
609:     PetscCall(DMPlexGetAdjacency(dm, p, &adjSize, &adj));
610:     for (n = 0; n < neighbors; ++n) {
611:       const PetscInt remoteRank = ranks[noff + n];

613:       if (remoteRank == rank) continue;
614:       for (a = 0; a < adjSize; ++a) {
615:         PetscBool insert = PETSC_TRUE;

617:         for (el = 0; el < numExLabels; ++el) {
618:           PetscInt exVal;
619:           PetscCall(DMLabelGetValue(exLabel[el], adj[a], &exVal));
620:           if (exVal == exValue[el]) {
621:             insert = PETSC_FALSE;
622:             break;
623:           }
624:         }
625:         if (insert) PetscCall(DMLabelSetValue(ovAdjByRank, adj[a], remoteRank));
626:       }
627:     }
628:     PetscCall(PetscFree(adj));
629:   }
630:   PetscFunctionReturn(PETSC_SUCCESS);
631: }

633: /*@C
634:   DMPlexCreateOverlapLabelFromLabels - Compute a label indicating what overlap points should be sent to new processes

636:   Collective

638:   Input Parameters:
639: + dm          - The `DM`
640: . numLabels   - The number of labels to draw candidate points from
641: . label       - An array of labels containing candidate points
642: . value       - An array of label values marking the candidate points
643: . numExLabels - The number of labels to use for exclusion
644: . exLabel     - An array of labels indicating points to be excluded, or `NULL`
645: . exValue     - An array of label values to be excluded, or `NULL`
646: . rootSection - The number of leaves for a given root point
647: . rootrank    - The rank of each edge into the root point
648: . leafSection - The number of processes sharing a given leaf point
649: - leafrank    - The rank of each process sharing a leaf point

651:   Output Parameter:
652: . ovLabel - `DMLabel` containing remote overlap contributions as point/rank pairings

654:   Level: developer

656:   Note:
657:   The candidate points are only added to the overlap if they are adjacent to a shared point

659: .seealso: `DMPLEX`, `DMPlexCreateOverlapLabel()`, `DMPlexGetAdjacency()`, `DMPlexDistributeOwnership()`, `DMPlexDistribute()`
660: @*/
661: PetscErrorCode DMPlexCreateOverlapLabelFromLabels(DM dm, PetscInt numLabels, const DMLabel label[], const PetscInt value[], PetscInt numExLabels, const DMLabel exLabel[], const PetscInt exValue[], PetscSection rootSection, IS rootrank, PetscSection leafSection, IS leafrank, DMLabel *ovLabel)
662: {
663:   MPI_Comm           comm;
664:   DMLabel            ovAdjByRank; /* A DMLabel containing all points adjacent to shared points, separated by rank (value in label) */
665:   PetscSF            sfPoint;
666:   const PetscSFNode *remote;
667:   const PetscInt    *local;
668:   const PetscInt    *nrank, *rrank;
669:   PetscInt          *adj = NULL;
670:   PetscInt           pStart, pEnd, p, sStart, sEnd, nleaves, l, el;
671:   PetscMPIInt        rank, size;
672:   PetscBool          flg;

674:   PetscFunctionBegin;
675:   *ovLabel = NULL;
676:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
677:   PetscCallMPI(MPI_Comm_size(comm, &size));
678:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
679:   if (size == 1) PetscFunctionReturn(PETSC_SUCCESS);
680:   PetscCall(DMGetPointSF(dm, &sfPoint));
681:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
682:   PetscCall(PetscSectionGetChart(leafSection, &sStart, &sEnd));
683:   PetscCall(PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote));
684:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Overlap adjacency", &ovAdjByRank));
685:   PetscCall(ISGetIndices(rootrank, &rrank));
686:   PetscCall(ISGetIndices(leafrank, &nrank));
687:   for (l = 0; l < numLabels; ++l) {
688:     IS              valIS;
689:     const PetscInt *points;
690:     PetscInt        n;

692:     PetscCall(DMLabelGetStratumIS(label[l], value[l], &valIS));
693:     if (!valIS) continue;
694:     PetscCall(ISGetIndices(valIS, &points));
695:     PetscCall(ISGetLocalSize(valIS, &n));
696:     for (PetscInt i = 0; i < n; ++i) {
697:       const PetscInt p = points[i];

699:       if ((p >= sStart) && (p < sEnd)) {
700:         PetscInt loc, adjSize = PETSC_DETERMINE;

702:         /* Handle leaves: shared with the root point */
703:         if (local) PetscCall(PetscFindInt(p, nleaves, local, &loc));
704:         else loc = (p >= 0 && p < nleaves) ? p : -1;
705:         if (loc >= 0) {
706:           const PetscInt remoteRank = remote[loc].rank;

708:           PetscCall(DMPlexGetAdjacency(dm, p, &adjSize, &adj));
709:           for (PetscInt a = 0; a < adjSize; ++a) {
710:             PetscBool insert = PETSC_TRUE;

712:             for (el = 0; el < numExLabels; ++el) {
713:               PetscInt exVal;
714:               PetscCall(DMLabelGetValue(exLabel[el], adj[a], &exVal));
715:               if (exVal == exValue[el]) {
716:                 insert = PETSC_FALSE;
717:                 break;
718:               }
719:             }
720:             if (insert) PetscCall(DMLabelSetValue(ovAdjByRank, adj[a], remoteRank));
721:           }
722:         }
723:         /* Some leaves share a root with other leaves on different processes */
724:         PetscCall(HandlePoint_Private(dm, p, leafSection, nrank, numExLabels, exLabel, exValue, ovAdjByRank));
725:       }
726:       /* Roots are shared with leaves */
727:       PetscCall(HandlePoint_Private(dm, p, rootSection, rrank, numExLabels, exLabel, exValue, ovAdjByRank));
728:     }
729:     PetscCall(ISRestoreIndices(valIS, &points));
730:     PetscCall(ISDestroy(&valIS));
731:   }
732:   PetscCall(PetscFree(adj));
733:   PetscCall(ISRestoreIndices(rootrank, &rrank));
734:   PetscCall(ISRestoreIndices(leafrank, &nrank));
735:   /* We require the closure in the overlap */
736:   PetscCall(DMPlexPartitionLabelClosure(dm, ovAdjByRank));
737:   PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-overlap_view", &flg));
738:   if (flg) {
739:     PetscViewer viewer;
740:     PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm), &viewer));
741:     PetscCall(DMLabelView(ovAdjByRank, viewer));
742:   }
743:   /* Invert sender to receiver label */
744:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Overlap label", ovLabel));
745:   PetscCall(DMPlexPartitionLabelInvert(dm, ovAdjByRank, NULL, *ovLabel));
746:   /* Add owned points, except for shared local points */
747:   for (p = pStart; p < pEnd; ++p) PetscCall(DMLabelSetValue(*ovLabel, p, rank));
748:   for (l = 0; l < nleaves; ++l) {
749:     PetscCall(DMLabelClearValue(*ovLabel, local[l], rank));
750:     PetscCall(DMLabelSetValue(*ovLabel, remote[l].index, remote[l].rank));
751:   }
752:   /* Clean up */
753:   PetscCall(DMLabelDestroy(&ovAdjByRank));
754:   PetscFunctionReturn(PETSC_SUCCESS);
755: }

757: /*@
758:   DMPlexCreateOverlapMigrationSF - Create a `PetscSF` describing the new mesh distribution to make the overlap described by the input `PetscSF`

760:   Collective

762:   Input Parameters:
763: + dm        - The `DM`
764: - overlapSF - The `PetscSF` mapping ghost points in overlap to owner points on other processes

766:   Output Parameter:
767: . migrationSF - A `PetscSF` that maps original points in old locations to points in new locations

769:   Level: developer

771: .seealso: `DMPLEX`, `DMPlexCreateOverlapLabel()`, `DMPlexDistributeOverlap()`, `DMPlexDistribute()`
772: @*/
773: PetscErrorCode DMPlexCreateOverlapMigrationSF(DM dm, PetscSF overlapSF, PetscSF *migrationSF)
774: {
775:   MPI_Comm           comm;
776:   PetscMPIInt        rank, size;
777:   PetscInt           d, dim, p, pStart, pEnd, nroots, nleaves, newLeaves, point, numSharedPoints;
778:   PetscInt          *pointDepths, *remoteDepths, *ilocal;
779:   PetscInt          *depthRecv, *depthShift, *depthIdx;
780:   PetscSFNode       *iremote;
781:   PetscSF            pointSF;
782:   const PetscInt    *sharedLocal;
783:   const PetscSFNode *overlapRemote, *sharedRemote;

785:   PetscFunctionBegin;
787:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
788:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
789:   PetscCallMPI(MPI_Comm_size(comm, &size));
790:   PetscCall(DMGetDimension(dm, &dim));

792:   /* Before building the migration SF we need to know the new stratum offsets */
793:   PetscCall(PetscSFGetGraph(overlapSF, &nroots, &nleaves, NULL, &overlapRemote));
794:   PetscCall(PetscMalloc2(nroots, &pointDepths, nleaves, &remoteDepths));
795:   for (d = 0; d < dim + 1; d++) {
796:     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
797:     for (p = pStart; p < pEnd; p++) pointDepths[p] = d;
798:   }
799:   for (p = 0; p < nleaves; p++) remoteDepths[p] = -1;
800:   PetscCall(PetscSFBcastBegin(overlapSF, MPIU_INT, pointDepths, remoteDepths, MPI_REPLACE));
801:   PetscCall(PetscSFBcastEnd(overlapSF, MPIU_INT, pointDepths, remoteDepths, MPI_REPLACE));

803:   /* Count received points in each stratum and compute the internal strata shift */
804:   PetscCall(PetscMalloc3(dim + 1, &depthRecv, dim + 1, &depthShift, dim + 1, &depthIdx));
805:   for (d = 0; d < dim + 1; d++) depthRecv[d] = 0;
806:   for (p = 0; p < nleaves; p++) depthRecv[remoteDepths[p]]++;
807:   depthShift[dim] = 0;
808:   for (d = 0; d < dim; d++) depthShift[d] = depthRecv[dim];
809:   for (d = 1; d < dim; d++) depthShift[d] += depthRecv[0];
810:   for (d = dim - 2; d > 0; d--) depthShift[d] += depthRecv[d + 1];
811:   for (d = 0; d < dim + 1; d++) {
812:     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
813:     depthIdx[d] = pStart + depthShift[d];
814:   }

816:   /* Form the overlap SF build an SF that describes the full overlap migration SF */
817:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
818:   newLeaves = pEnd - pStart + nleaves;
819:   PetscCall(PetscMalloc1(newLeaves, &ilocal));
820:   PetscCall(PetscMalloc1(newLeaves, &iremote));
821:   /* First map local points to themselves */
822:   for (d = 0; d < dim + 1; d++) {
823:     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
824:     for (p = pStart; p < pEnd; p++) {
825:       point                = p + depthShift[d];
826:       ilocal[point]        = point;
827:       iremote[point].index = p;
828:       iremote[point].rank  = rank;
829:       depthIdx[d]++;
830:     }
831:   }

833:   /* Add in the remote roots for currently shared points */
834:   PetscCall(DMGetPointSF(dm, &pointSF));
835:   PetscCall(PetscSFGetGraph(pointSF, NULL, &numSharedPoints, &sharedLocal, &sharedRemote));
836:   for (d = 0; d < dim + 1; d++) {
837:     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
838:     for (p = 0; p < numSharedPoints; p++) {
839:       if (pStart <= sharedLocal[p] && sharedLocal[p] < pEnd) {
840:         point                = sharedLocal[p] + depthShift[d];
841:         iremote[point].index = sharedRemote[p].index;
842:         iremote[point].rank  = sharedRemote[p].rank;
843:       }
844:     }
845:   }

847:   /* Now add the incoming overlap points */
848:   for (p = 0; p < nleaves; p++) {
849:     point                = depthIdx[remoteDepths[p]];
850:     ilocal[point]        = point;
851:     iremote[point].index = overlapRemote[p].index;
852:     iremote[point].rank  = overlapRemote[p].rank;
853:     depthIdx[remoteDepths[p]]++;
854:   }
855:   PetscCall(PetscFree2(pointDepths, remoteDepths));

857:   PetscCall(PetscSFCreate(comm, migrationSF));
858:   PetscCall(PetscObjectSetName((PetscObject)*migrationSF, "Overlap Migration SF"));
859:   PetscCall(PetscSFSetFromOptions(*migrationSF));
860:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
861:   PetscCall(PetscSFSetGraph(*migrationSF, pEnd - pStart, newLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));

863:   PetscCall(PetscFree3(depthRecv, depthShift, depthIdx));
864:   PetscFunctionReturn(PETSC_SUCCESS);
865: }

867: /*@
868:   DMPlexStratifyMigrationSF - Rearrange the leaves of a migration sf for stratification.

870:   Input Parameters:
871: + dm - The DM
872: - sf - A star forest with non-ordered leaves, usually defining a DM point migration

874:   Output Parameter:
875: . migrationSF - A star forest with added leaf indirection that ensures the resulting DM is stratified

877:   Level: developer

879:   Note:
880:   This lexicographically sorts by (depth, cellType)

882: .seealso: `DMPLEX`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexDistributeOverlap()`
883: @*/
884: PetscErrorCode DMPlexStratifyMigrationSF(DM dm, PetscSF sf, PetscSF *migrationSF)
885: {
886:   MPI_Comm           comm;
887:   PetscMPIInt        rank, size;
888:   PetscInt           d, ldepth, depth, dim, p, pStart, pEnd, nroots, nleaves;
889:   PetscSFNode       *pointDepths, *remoteDepths;
890:   PetscInt          *ilocal;
891:   PetscInt          *depthRecv, *depthShift, *depthIdx;
892:   PetscInt          *ctRecv, *ctShift, *ctIdx;
893:   const PetscSFNode *iremote;

895:   PetscFunctionBegin;
897:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
898:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
899:   PetscCallMPI(MPI_Comm_size(comm, &size));
900:   PetscCall(DMPlexGetDepth(dm, &ldepth));
901:   PetscCall(DMGetDimension(dm, &dim));
902:   PetscCallMPI(MPIU_Allreduce(&ldepth, &depth, 1, MPIU_INT, MPI_MAX, comm));
903:   PetscCheck(!(ldepth >= 0) || !(depth != ldepth), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent Plex depth %" PetscInt_FMT " != %" PetscInt_FMT, ldepth, depth);
904:   PetscCall(PetscLogEventBegin(DMPLEX_PartStratSF, dm, 0, 0, 0));

906:   /* Before building the migration SF we need to know the new stratum offsets */
907:   PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, &iremote));
908:   PetscCall(PetscMalloc2(nroots, &pointDepths, nleaves, &remoteDepths));
909:   for (d = 0; d < depth + 1; ++d) {
910:     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
911:     for (p = pStart; p < pEnd; ++p) {
912:       DMPolytopeType ct;

914:       PetscCall(DMPlexGetCellType(dm, p, &ct));
915:       pointDepths[p].index = d;
916:       pointDepths[p].rank  = ct;
917:     }
918:   }
919:   for (p = 0; p < nleaves; ++p) {
920:     remoteDepths[p].index = -1;
921:     remoteDepths[p].rank  = -1;
922:   }
923:   PetscCall(PetscSFBcastBegin(sf, MPIU_SF_NODE, pointDepths, remoteDepths, MPI_REPLACE));
924:   PetscCall(PetscSFBcastEnd(sf, MPIU_SF_NODE, pointDepths, remoteDepths, MPI_REPLACE));
925:   /* Count received points in each stratum and compute the internal strata shift */
926:   PetscCall(PetscCalloc6(depth + 1, &depthRecv, depth + 1, &depthShift, depth + 1, &depthIdx, DM_NUM_POLYTOPES, &ctRecv, DM_NUM_POLYTOPES, &ctShift, DM_NUM_POLYTOPES, &ctIdx));
927:   for (p = 0; p < nleaves; ++p) {
928:     if (remoteDepths[p].rank < 0) {
929:       ++depthRecv[remoteDepths[p].index];
930:     } else {
931:       ++ctRecv[remoteDepths[p].rank];
932:     }
933:   }
934:   {
935:     PetscInt depths[4], dims[4], shift = 0, i, c;

937:     /* Cells (depth), Vertices (0), Faces (depth-1), Edges (1)
938:          Consider DM_POLYTOPE_FV_GHOST, DM_POLYTOPE_INTERIOR_GHOST and DM_POLYTOPE_UNKNOWN_CELL as cells
939:      */
940:     depths[0] = depth;
941:     depths[1] = 0;
942:     depths[2] = depth - 1;
943:     depths[3] = 1;
944:     dims[0]   = dim;
945:     dims[1]   = 0;
946:     dims[2]   = dim - 1;
947:     dims[3]   = 1;
948:     for (i = 0; i <= depth; ++i) {
949:       const PetscInt dep = depths[i];
950:       const PetscInt dim = dims[i];

952:       for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
953:         if (DMPolytopeTypeGetDim((DMPolytopeType)c) != dim && !(i == 0 && (c == DM_POLYTOPE_FV_GHOST || c == DM_POLYTOPE_INTERIOR_GHOST || c == DM_POLYTOPE_UNKNOWN_CELL))) continue;
954:         ctShift[c] = shift;
955:         shift += ctRecv[c];
956:       }
957:       depthShift[dep] = shift;
958:       shift += depthRecv[dep];
959:     }
960:     for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
961:       const PetscInt ctDim = DMPolytopeTypeGetDim((DMPolytopeType)c);

963:       if ((ctDim < 0 || ctDim > dim) && c != DM_POLYTOPE_FV_GHOST && c != DM_POLYTOPE_INTERIOR_GHOST && c != DM_POLYTOPE_UNKNOWN_CELL) {
964:         ctShift[c] = shift;
965:         shift += ctRecv[c];
966:       }
967:     }
968:   }
969:   /* Derive a new local permutation based on stratified indices */
970:   PetscCall(PetscMalloc1(nleaves, &ilocal));
971:   for (p = 0; p < nleaves; ++p) {
972:     const DMPolytopeType ct = (DMPolytopeType)remoteDepths[p].rank;

974:     ilocal[p] = ctShift[ct] + ctIdx[ct];
975:     ++ctIdx[ct];
976:   }
977:   PetscCall(PetscSFCreate(comm, migrationSF));
978:   PetscCall(PetscObjectSetName((PetscObject)*migrationSF, "Migration SF"));
979:   PetscCall(PetscSFSetGraph(*migrationSF, nroots, nleaves, ilocal, PETSC_OWN_POINTER, (PetscSFNode *)iremote, PETSC_COPY_VALUES));
980:   PetscCall(PetscFree2(pointDepths, remoteDepths));
981:   PetscCall(PetscFree6(depthRecv, depthShift, depthIdx, ctRecv, ctShift, ctIdx));
982:   PetscCall(PetscLogEventEnd(DMPLEX_PartStratSF, dm, 0, 0, 0));
983:   PetscFunctionReturn(PETSC_SUCCESS);
984: }

986: /*@
987:   DMPlexDistributeField - Distribute field data to match a given `PetscSF`, usually the `PetscSF` from mesh distribution

989:   Collective

991:   Input Parameters:
992: + dm              - The `DMPLEX` object
993: . pointSF         - The `PetscSF` describing the communication pattern
994: . originalSection - The `PetscSection` for existing data layout
995: - originalVec     - The existing data in a local vector

997:   Output Parameters:
998: + newSection - The `PetscSF` describing the new data layout
999: - newVec     - The new data in a local vector

1001:   Level: developer

1003: .seealso: `DMPLEX`, `DMPlexDistribute()`, `DMPlexDistributeFieldIS()`, `DMPlexDistributeData()`
1004: @*/
1005: PetscErrorCode DMPlexDistributeField(DM dm, PetscSF pointSF, PetscSection originalSection, Vec originalVec, PetscSection newSection, Vec newVec)
1006: {
1007:   PetscSF      fieldSF;
1008:   PetscInt    *remoteOffsets, fieldSize;
1009:   PetscScalar *originalValues, *newValues;

1011:   PetscFunctionBegin;
1012:   PetscCall(PetscLogEventBegin(DMPLEX_DistributeField, dm, 0, 0, 0));
1013:   PetscCall(PetscSFDistributeSection(pointSF, originalSection, &remoteOffsets, newSection));

1015:   PetscCall(PetscSectionGetStorageSize(newSection, &fieldSize));
1016:   PetscCall(VecSetSizes(newVec, fieldSize, PETSC_DETERMINE));
1017:   PetscCall(VecSetType(newVec, dm->vectype));

1019:   PetscCall(VecGetArray(originalVec, &originalValues));
1020:   PetscCall(VecGetArray(newVec, &newValues));
1021:   PetscCall(PetscSFCreateSectionSF(pointSF, originalSection, remoteOffsets, newSection, &fieldSF));
1022:   PetscCall(PetscFree(remoteOffsets));
1023:   PetscCall(PetscSFBcastBegin(fieldSF, MPIU_SCALAR, originalValues, newValues, MPI_REPLACE));
1024:   PetscCall(PetscSFBcastEnd(fieldSF, MPIU_SCALAR, originalValues, newValues, MPI_REPLACE));
1025:   PetscCall(PetscSFDestroy(&fieldSF));
1026:   PetscCall(VecRestoreArray(newVec, &newValues));
1027:   PetscCall(VecRestoreArray(originalVec, &originalValues));
1028:   PetscCall(PetscLogEventEnd(DMPLEX_DistributeField, dm, 0, 0, 0));
1029:   PetscFunctionReturn(PETSC_SUCCESS);
1030: }

1032: /*@
1033:   DMPlexDistributeFieldIS - Distribute field data to match a given `PetscSF`, usually the `PetscSF` from mesh distribution

1035:   Collective

1037:   Input Parameters:
1038: + dm              - The `DMPLEX` object
1039: . pointSF         - The `PetscSF` describing the communication pattern
1040: . originalSection - The `PetscSection` for existing data layout
1041: - originalIS      - The existing data

1043:   Output Parameters:
1044: + newSection - The `PetscSF` describing the new data layout
1045: - newIS      - The new data

1047:   Level: developer

1049: .seealso: `DMPLEX`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexDistributeData()`
1050: @*/
1051: PetscErrorCode DMPlexDistributeFieldIS(DM dm, PetscSF pointSF, PetscSection originalSection, IS originalIS, PetscSection newSection, IS *newIS)
1052: {
1053:   PetscSF         fieldSF;
1054:   PetscInt       *newValues, *remoteOffsets, fieldSize;
1055:   const PetscInt *originalValues;

1057:   PetscFunctionBegin;
1058:   PetscCall(PetscLogEventBegin(DMPLEX_DistributeField, dm, 0, 0, 0));
1059:   PetscCall(PetscSFDistributeSection(pointSF, originalSection, &remoteOffsets, newSection));

1061:   PetscCall(PetscSectionGetStorageSize(newSection, &fieldSize));
1062:   PetscCall(PetscMalloc1(fieldSize, &newValues));

1064:   PetscCall(ISGetIndices(originalIS, &originalValues));
1065:   PetscCall(PetscSFCreateSectionSF(pointSF, originalSection, remoteOffsets, newSection, &fieldSF));
1066:   PetscCall(PetscFree(remoteOffsets));
1067:   PetscCall(PetscSFBcastBegin(fieldSF, MPIU_INT, (PetscInt *)originalValues, newValues, MPI_REPLACE));
1068:   PetscCall(PetscSFBcastEnd(fieldSF, MPIU_INT, (PetscInt *)originalValues, newValues, MPI_REPLACE));
1069:   PetscCall(PetscSFDestroy(&fieldSF));
1070:   PetscCall(ISRestoreIndices(originalIS, &originalValues));
1071:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pointSF), fieldSize, newValues, PETSC_OWN_POINTER, newIS));
1072:   PetscCall(PetscLogEventEnd(DMPLEX_DistributeField, dm, 0, 0, 0));
1073:   PetscFunctionReturn(PETSC_SUCCESS);
1074: }

1076: /*@
1077:   DMPlexDistributeData - Distribute field data to match a given `PetscSF`, usually the `PetscSF` from mesh distribution

1079:   Collective

1081:   Input Parameters:
1082: + dm              - The `DMPLEX` object
1083: . pointSF         - The `PetscSF` describing the communication pattern
1084: . originalSection - The `PetscSection` for existing data layout
1085: . datatype        - The type of data
1086: - originalData    - The existing data

1088:   Output Parameters:
1089: + newSection - The `PetscSection` describing the new data layout
1090: - newData    - The new data

1092:   Level: developer

1094: .seealso: `DMPLEX`, `DMPlexDistribute()`, `DMPlexDistributeField()`
1095: @*/
1096: PetscErrorCode DMPlexDistributeData(DM dm, PetscSF pointSF, PetscSection originalSection, MPI_Datatype datatype, void *originalData, PetscSection newSection, void **newData)
1097: {
1098:   PetscSF     fieldSF;
1099:   PetscInt   *remoteOffsets, fieldSize;
1100:   PetscMPIInt dataSize;

1102:   PetscFunctionBegin;
1103:   PetscCall(PetscLogEventBegin(DMPLEX_DistributeData, dm, 0, 0, 0));
1104:   PetscCall(PetscSFDistributeSection(pointSF, originalSection, &remoteOffsets, newSection));

1106:   PetscCall(PetscSectionGetStorageSize(newSection, &fieldSize));
1107:   PetscCallMPI(MPI_Type_size(datatype, &dataSize));
1108:   PetscCall(PetscMalloc(fieldSize * dataSize, newData));

1110:   PetscCall(PetscSFCreateSectionSF(pointSF, originalSection, remoteOffsets, newSection, &fieldSF));
1111:   PetscCall(PetscFree(remoteOffsets));
1112:   PetscCall(PetscSFBcastBegin(fieldSF, datatype, originalData, *newData, MPI_REPLACE));
1113:   PetscCall(PetscSFBcastEnd(fieldSF, datatype, originalData, *newData, MPI_REPLACE));
1114:   PetscCall(PetscSFDestroy(&fieldSF));
1115:   PetscCall(PetscLogEventEnd(DMPLEX_DistributeData, dm, 0, 0, 0));
1116:   PetscFunctionReturn(PETSC_SUCCESS);
1117: }

1119: static PetscErrorCode DMPlexDistributeCones(DM dm, PetscSF migrationSF, ISLocalToGlobalMapping original, ISLocalToGlobalMapping renumbering, DM dmParallel)
1120: {
1121:   DM_Plex     *pmesh = (DM_Plex *)dmParallel->data;
1122:   MPI_Comm     comm;
1123:   PetscSF      coneSF;
1124:   PetscSection originalConeSection, newConeSection;
1125:   PetscInt    *remoteOffsets, *cones, *globCones, *newCones, newConesSize;
1126:   PetscBool    flg;

1128:   PetscFunctionBegin;
1131:   PetscCall(PetscLogEventBegin(DMPLEX_DistributeCones, dm, 0, 0, 0));
1132:   /* Distribute cone section */
1133:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1134:   PetscCall(DMPlexGetConeSection(dm, &originalConeSection));
1135:   PetscCall(DMPlexGetConeSection(dmParallel, &newConeSection));
1136:   PetscCall(PetscSFDistributeSection(migrationSF, originalConeSection, &remoteOffsets, newConeSection));
1137:   PetscCall(DMSetUp(dmParallel));
1138:   /* Communicate and renumber cones */
1139:   PetscCall(PetscSFCreateSectionSF(migrationSF, originalConeSection, remoteOffsets, newConeSection, &coneSF));
1140:   PetscCall(PetscFree(remoteOffsets));
1141:   PetscCall(DMPlexGetCones(dm, &cones));
1142:   if (original) {
1143:     PetscInt numCones;

1145:     PetscCall(PetscSectionGetStorageSize(originalConeSection, &numCones));
1146:     PetscCall(PetscMalloc1(numCones, &globCones));
1147:     PetscCall(ISLocalToGlobalMappingApplyBlock(original, numCones, cones, globCones));
1148:   } else {
1149:     globCones = cones;
1150:   }
1151:   PetscCall(DMPlexGetCones(dmParallel, &newCones));
1152:   PetscCall(PetscSFBcastBegin(coneSF, MPIU_INT, globCones, newCones, MPI_REPLACE));
1153:   PetscCall(PetscSFBcastEnd(coneSF, MPIU_INT, globCones, newCones, MPI_REPLACE));
1154:   if (original) PetscCall(PetscFree(globCones));
1155:   PetscCall(PetscSectionGetStorageSize(newConeSection, &newConesSize));
1156:   PetscCall(ISGlobalToLocalMappingApplyBlock(renumbering, IS_GTOLM_MASK, newConesSize, newCones, NULL, newCones));
1157:   if (PetscDefined(USE_DEBUG)) {
1158:     PetscInt  p;
1159:     PetscBool valid = PETSC_TRUE;
1160:     for (p = 0; p < newConesSize; ++p) {
1161:       if (newCones[p] < 0) {
1162:         valid = PETSC_FALSE;
1163:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Point %" PetscInt_FMT " not in overlap SF\n", PetscGlobalRank, p));
1164:       }
1165:     }
1166:     PetscCheck(valid, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid global to local map");
1167:   }
1168:   PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-cones_view", &flg));
1169:   if (flg) {
1170:     PetscCall(PetscPrintf(comm, "Serial Cone Section:\n"));
1171:     PetscCall(PetscSectionView(originalConeSection, PETSC_VIEWER_STDOUT_(comm)));
1172:     PetscCall(PetscPrintf(comm, "Parallel Cone Section:\n"));
1173:     PetscCall(PetscSectionView(newConeSection, PETSC_VIEWER_STDOUT_(comm)));
1174:     PetscCall(PetscSFView(coneSF, NULL));
1175:   }
1176:   PetscCall(DMPlexGetConeOrientations(dm, &cones));
1177:   PetscCall(DMPlexGetConeOrientations(dmParallel, &newCones));
1178:   PetscCall(PetscSFBcastBegin(coneSF, MPIU_INT, cones, newCones, MPI_REPLACE));
1179:   PetscCall(PetscSFBcastEnd(coneSF, MPIU_INT, cones, newCones, MPI_REPLACE));
1180:   PetscCall(PetscSFDestroy(&coneSF));
1181:   PetscCall(PetscLogEventEnd(DMPLEX_DistributeCones, dm, 0, 0, 0));
1182:   /* Create supports and stratify DMPlex */
1183:   {
1184:     PetscInt pStart, pEnd;

1186:     PetscCall(PetscSectionGetChart(pmesh->coneSection, &pStart, &pEnd));
1187:     PetscCall(PetscSectionSetChart(pmesh->supportSection, pStart, pEnd));
1188:   }
1189:   PetscCall(DMPlexSymmetrize(dmParallel));
1190:   PetscCall(DMPlexStratify(dmParallel));
1191:   {
1192:     PetscBool useCone, useClosure, useAnchors;

1194:     PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure));
1195:     PetscCall(DMSetBasicAdjacency(dmParallel, useCone, useClosure));
1196:     PetscCall(DMPlexGetAdjacencyUseAnchors(dm, &useAnchors));
1197:     PetscCall(DMPlexSetAdjacencyUseAnchors(dmParallel, useAnchors));
1198:   }
1199:   PetscFunctionReturn(PETSC_SUCCESS);
1200: }

1202: static PetscErrorCode DMPlexDistributeCoordinates(DM dm, PetscSF migrationSF, DM dmParallel)
1203: {
1204:   MPI_Comm         comm;
1205:   DM               cdm, cdmParallel;
1206:   PetscSection     originalCoordSection, newCoordSection;
1207:   Vec              originalCoordinates, newCoordinates;
1208:   PetscInt         bs;
1209:   const char      *name;
1210:   const PetscReal *maxCell, *Lstart, *L;

1212:   PetscFunctionBegin;

1216:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1217:   PetscCall(DMGetCoordinateDM(dm, &cdm));
1218:   PetscCall(DMGetCoordinateDM(dmParallel, &cdmParallel));
1219:   PetscCall(DMCopyDisc(cdm, cdmParallel));
1220:   PetscCall(DMGetCoordinateSection(dm, &originalCoordSection));
1221:   PetscCall(DMGetCoordinateSection(dmParallel, &newCoordSection));
1222:   PetscCall(DMGetCoordinatesLocal(dm, &originalCoordinates));
1223:   if (originalCoordinates) {
1224:     PetscCall(VecCreate(PETSC_COMM_SELF, &newCoordinates));
1225:     PetscCall(PetscObjectGetName((PetscObject)originalCoordinates, &name));
1226:     PetscCall(PetscObjectSetName((PetscObject)newCoordinates, name));

1228:     PetscCall(DMPlexDistributeField(dm, migrationSF, originalCoordSection, originalCoordinates, newCoordSection, newCoordinates));
1229:     PetscCall(DMSetCoordinatesLocal(dmParallel, newCoordinates));
1230:     PetscCall(VecGetBlockSize(originalCoordinates, &bs));
1231:     PetscCall(VecSetBlockSize(newCoordinates, bs));
1232:     PetscCall(VecDestroy(&newCoordinates));
1233:   }

1235:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1236:   PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
1237:   PetscCall(DMSetPeriodicity(dmParallel, maxCell, Lstart, L));
1238:   PetscCall(DMGetCellCoordinateDM(dm, &cdm));
1239:   if (cdm) {
1240:     PetscCall(DMClone(dmParallel, &cdmParallel));
1241:     PetscCall(DMSetCellCoordinateDM(dmParallel, cdmParallel));
1242:     PetscCall(DMCopyDisc(cdm, cdmParallel));
1243:     PetscCall(DMDestroy(&cdmParallel));
1244:     PetscCall(DMGetCellCoordinateSection(dm, &originalCoordSection));
1245:     PetscCall(DMGetCellCoordinatesLocal(dm, &originalCoordinates));
1246:     PetscCall(PetscSectionCreate(comm, &newCoordSection));
1247:     if (originalCoordinates) {
1248:       PetscCall(VecCreate(PETSC_COMM_SELF, &newCoordinates));
1249:       PetscCall(PetscObjectGetName((PetscObject)originalCoordinates, &name));
1250:       PetscCall(PetscObjectSetName((PetscObject)newCoordinates, name));

1252:       PetscCall(DMPlexDistributeField(dm, migrationSF, originalCoordSection, originalCoordinates, newCoordSection, newCoordinates));
1253:       PetscCall(VecGetBlockSize(originalCoordinates, &bs));
1254:       PetscCall(VecSetBlockSize(newCoordinates, bs));
1255:       PetscCall(DMSetCellCoordinateSection(dmParallel, bs, newCoordSection));
1256:       PetscCall(DMSetCellCoordinatesLocal(dmParallel, newCoordinates));
1257:       PetscCall(VecDestroy(&newCoordinates));
1258:     }
1259:     PetscCall(PetscSectionDestroy(&newCoordSection));
1260:   }
1261:   PetscFunctionReturn(PETSC_SUCCESS);
1262: }

1264: static PetscErrorCode DMPlexDistributeLabels(DM dm, PetscSF migrationSF, DM dmParallel)
1265: {
1266:   DM_Plex         *mesh = (DM_Plex *)dm->data;
1267:   MPI_Comm         comm;
1268:   DMLabel          depthLabel;
1269:   PetscMPIInt      rank;
1270:   PetscInt         depth, d, numLabels, numLocalLabels, l;
1271:   PetscBool        hasLabels  = PETSC_FALSE, lsendDepth, sendDepth;
1272:   PetscObjectState depthState = -1;

1274:   PetscFunctionBegin;

1278:   PetscCall(PetscLogEventBegin(DMPLEX_DistributeLabels, dm, 0, 0, 0));
1279:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1280:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

1282:   /* If the user has changed the depth label, communicate it instead */
1283:   PetscCall(DMPlexGetDepth(dm, &depth));
1284:   PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
1285:   if (depthLabel) PetscCall(PetscObjectStateGet((PetscObject)depthLabel, &depthState));
1286:   lsendDepth = mesh->depthState != depthState ? PETSC_TRUE : PETSC_FALSE;
1287:   PetscCallMPI(MPIU_Allreduce(&lsendDepth, &sendDepth, 1, MPI_C_BOOL, MPI_LOR, comm));
1288:   if (sendDepth) {
1289:     PetscCall(DMPlexGetDepthLabel(dmParallel, &dmParallel->depthLabel));
1290:     PetscCall(DMRemoveLabelBySelf(dmParallel, &dmParallel->depthLabel, PETSC_FALSE));
1291:   }
1292:   /* Everyone must have either the same number of labels, or none */
1293:   PetscCall(DMGetNumLabels(dm, &numLocalLabels));
1294:   numLabels = numLocalLabels;
1295:   PetscCallMPI(MPI_Bcast(&numLabels, 1, MPIU_INT, 0, comm));
1296:   if (numLabels == numLocalLabels) hasLabels = PETSC_TRUE;
1297:   for (l = 0; l < numLabels; ++l) {
1298:     DMLabel     label = NULL, labelNew = NULL;
1299:     PetscBool   isDepth, lisOutput     = PETSC_TRUE, isOutput;
1300:     const char *name = NULL;

1302:     if (hasLabels) {
1303:       PetscCall(DMGetLabelByNum(dm, l, &label));
1304:       /* Skip "depth" because it is recreated */
1305:       PetscCall(PetscObjectGetName((PetscObject)label, &name));
1306:       PetscCall(PetscStrcmp(name, "depth", &isDepth));
1307:     } else {
1308:       isDepth = PETSC_FALSE;
1309:     }
1310:     PetscCallMPI(MPI_Bcast(&isDepth, 1, MPI_C_BOOL, 0, comm));
1311:     if (isDepth && !sendDepth) continue;
1312:     PetscCall(DMLabelDistribute(label, migrationSF, &labelNew));
1313:     if (isDepth) {
1314:       /* Put in any missing strata which can occur if users are managing the depth label themselves */
1315:       PetscInt gdepth;

1317:       PetscCallMPI(MPIU_Allreduce(&depth, &gdepth, 1, MPIU_INT, MPI_MAX, comm));
1318:       PetscCheck(!(depth >= 0) || !(gdepth != depth), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent Plex depth %" PetscInt_FMT " != %" PetscInt_FMT, depth, gdepth);
1319:       for (d = 0; d <= gdepth; ++d) {
1320:         PetscBool has;

1322:         PetscCall(DMLabelHasStratum(labelNew, d, &has));
1323:         if (!has) PetscCall(DMLabelAddStratum(labelNew, d));
1324:       }
1325:     }
1326:     PetscCall(DMAddLabel(dmParallel, labelNew));
1327:     /* Put the output flag in the new label */
1328:     if (hasLabels) PetscCall(DMGetLabelOutput(dm, name, &lisOutput));
1329:     PetscCallMPI(MPIU_Allreduce(&lisOutput, &isOutput, 1, MPI_C_BOOL, MPI_LAND, comm));
1330:     PetscCall(PetscObjectGetName((PetscObject)labelNew, &name));
1331:     PetscCall(DMSetLabelOutput(dmParallel, name, isOutput));
1332:     PetscCall(DMLabelDestroy(&labelNew));
1333:   }
1334:   {
1335:     DMLabel ctLabel;

1337:     // Reset label for fast lookup
1338:     PetscCall(DMPlexGetCellTypeLabel(dmParallel, &ctLabel));
1339:     PetscCall(DMLabelMakeAllInvalid_Internal(ctLabel));
1340:   }
1341:   PetscCall(PetscLogEventEnd(DMPLEX_DistributeLabels, dm, 0, 0, 0));
1342:   PetscFunctionReturn(PETSC_SUCCESS);
1343: }

1345: static PetscErrorCode DMPlexDistributeSetupTree(DM dm, PetscSF migrationSF, ISLocalToGlobalMapping original, ISLocalToGlobalMapping renumbering, DM dmParallel)
1346: {
1347:   DM_Plex     *mesh  = (DM_Plex *)dm->data;
1348:   DM_Plex     *pmesh = (DM_Plex *)dmParallel->data;
1349:   MPI_Comm     comm;
1350:   DM           refTree;
1351:   PetscSection origParentSection, newParentSection;
1352:   PetscInt    *origParents, *origChildIDs;
1353:   PetscBool    flg;

1355:   PetscFunctionBegin;
1358:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));

1360:   /* Set up tree */
1361:   PetscCall(DMPlexGetReferenceTree(dm, &refTree));
1362:   PetscCall(DMPlexSetReferenceTree(dmParallel, refTree));
1363:   PetscCall(DMPlexGetTree(dm, &origParentSection, &origParents, &origChildIDs, NULL, NULL));
1364:   if (origParentSection) {
1365:     PetscInt  pStart, pEnd;
1366:     PetscInt *newParents, *newChildIDs, *globParents;
1367:     PetscInt *remoteOffsetsParents, newParentSize;
1368:     PetscSF   parentSF;

1370:     PetscCall(DMPlexGetChart(dmParallel, &pStart, &pEnd));
1371:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dmParallel), &newParentSection));
1372:     PetscCall(PetscSectionSetChart(newParentSection, pStart, pEnd));
1373:     PetscCall(PetscSFDistributeSection(migrationSF, origParentSection, &remoteOffsetsParents, newParentSection));
1374:     PetscCall(PetscSFCreateSectionSF(migrationSF, origParentSection, remoteOffsetsParents, newParentSection, &parentSF));
1375:     PetscCall(PetscFree(remoteOffsetsParents));
1376:     PetscCall(PetscSectionGetStorageSize(newParentSection, &newParentSize));
1377:     PetscCall(PetscMalloc2(newParentSize, &newParents, newParentSize, &newChildIDs));
1378:     if (original) {
1379:       PetscInt numParents;

1381:       PetscCall(PetscSectionGetStorageSize(origParentSection, &numParents));
1382:       PetscCall(PetscMalloc1(numParents, &globParents));
1383:       PetscCall(ISLocalToGlobalMappingApplyBlock(original, numParents, origParents, globParents));
1384:     } else {
1385:       globParents = origParents;
1386:     }
1387:     PetscCall(PetscSFBcastBegin(parentSF, MPIU_INT, globParents, newParents, MPI_REPLACE));
1388:     PetscCall(PetscSFBcastEnd(parentSF, MPIU_INT, globParents, newParents, MPI_REPLACE));
1389:     if (original) PetscCall(PetscFree(globParents));
1390:     PetscCall(PetscSFBcastBegin(parentSF, MPIU_INT, origChildIDs, newChildIDs, MPI_REPLACE));
1391:     PetscCall(PetscSFBcastEnd(parentSF, MPIU_INT, origChildIDs, newChildIDs, MPI_REPLACE));
1392:     PetscCall(ISGlobalToLocalMappingApplyBlock(renumbering, IS_GTOLM_MASK, newParentSize, newParents, NULL, newParents));
1393:     if (PetscDefined(USE_DEBUG)) {
1394:       PetscInt  p;
1395:       PetscBool valid = PETSC_TRUE;
1396:       for (p = 0; p < newParentSize; ++p) {
1397:         if (newParents[p] < 0) {
1398:           valid = PETSC_FALSE;
1399:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "Point %" PetscInt_FMT " not in overlap SF\n", p));
1400:         }
1401:       }
1402:       PetscCheck(valid, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid global to local map");
1403:     }
1404:     PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-parents_view", &flg));
1405:     if (flg) {
1406:       PetscCall(PetscPrintf(comm, "Serial Parent Section: \n"));
1407:       PetscCall(PetscSectionView(origParentSection, PETSC_VIEWER_STDOUT_(comm)));
1408:       PetscCall(PetscPrintf(comm, "Parallel Parent Section: \n"));
1409:       PetscCall(PetscSectionView(newParentSection, PETSC_VIEWER_STDOUT_(comm)));
1410:       PetscCall(PetscSFView(parentSF, NULL));
1411:     }
1412:     PetscCall(DMPlexSetTree(dmParallel, newParentSection, newParents, newChildIDs));
1413:     PetscCall(PetscSectionDestroy(&newParentSection));
1414:     PetscCall(PetscFree2(newParents, newChildIDs));
1415:     PetscCall(PetscSFDestroy(&parentSF));
1416:   }
1417:   pmesh->useAnchors = mesh->useAnchors;
1418:   PetscFunctionReturn(PETSC_SUCCESS);
1419: }

1421: /*@
1422:   DMPlexSetPartitionBalance - Should distribution of the `DM` attempt to balance the shared point partition?

1424:   Input Parameters:
1425: + dm  - The `DMPLEX` object
1426: - flg - Balance the partition?

1428:   Level: intermediate

1430: .seealso: `DMPLEX`, `DMPlexDistribute()`, `DMPlexGetPartitionBalance()`
1431: @*/
1432: PetscErrorCode DMPlexSetPartitionBalance(DM dm, PetscBool flg)
1433: {
1434:   DM_Plex *mesh = (DM_Plex *)dm->data;

1436:   PetscFunctionBegin;
1437:   mesh->partitionBalance = flg;
1438:   PetscFunctionReturn(PETSC_SUCCESS);
1439: }

1441: /*@
1442:   DMPlexGetPartitionBalance - Does distribution of the `DM` attempt to balance the shared point partition?

1444:   Input Parameter:
1445: . dm - The `DMPLEX` object

1447:   Output Parameter:
1448: . flg - Balance the partition?

1450:   Level: intermediate

1452: .seealso: `DMPLEX`, `DMPlexDistribute()`, `DMPlexSetPartitionBalance()`
1453: @*/
1454: PetscErrorCode DMPlexGetPartitionBalance(DM dm, PetscBool *flg)
1455: {
1456:   DM_Plex *mesh = (DM_Plex *)dm->data;

1458:   PetscFunctionBegin;
1460:   PetscAssertPointer(flg, 2);
1461:   *flg = mesh->partitionBalance;
1462:   PetscFunctionReturn(PETSC_SUCCESS);
1463: }

1465: typedef struct {
1466:   PetscInt vote, rank, index;
1467: } Petsc3Int;

1469: /* MaxLoc, but carry a third piece of information around */
1470: static void MPIAPI MaxLocCarry(void *in_, void *inout_, PetscMPIInt *len_, MPI_Datatype *dtype)
1471: {
1472:   Petsc3Int *a = (Petsc3Int *)inout_;
1473:   Petsc3Int *b = (Petsc3Int *)in_;
1474:   PetscInt   i, len = *len_;
1475:   for (i = 0; i < len; i++) {
1476:     if (a[i].vote < b[i].vote) {
1477:       a[i].vote  = b[i].vote;
1478:       a[i].rank  = b[i].rank;
1479:       a[i].index = b[i].index;
1480:     } else if (a[i].vote <= b[i].vote) {
1481:       if (a[i].rank >= b[i].rank) {
1482:         a[i].rank  = b[i].rank;
1483:         a[i].index = b[i].index;
1484:       }
1485:     }
1486:   }
1487: }

1489: /*@
1490:   DMPlexCreatePointSF - Build a point `PetscSF` from an `PetscSF` describing a point migration

1492:   Input Parameters:
1493: + dm          - The source `DMPLEX` object
1494: . migrationSF - The star forest that describes the parallel point remapping
1495: - ownership   - Flag causing a vote to determine point ownership

1497:   Output Parameter:
1498: . pointSF - The star forest describing the point overlap in the remapped `DM`

1500:   Level: developer

1502:   Note:
1503:   Output `pointSF` is guaranteed to have the array of local indices (leaves) sorted.

1505: .seealso: `DMPLEX`, `PetscSF`, `DM`, `DMPlexDistribute()`, `DMPlexDistributeOverlap()`
1506: @*/
1507: PetscErrorCode DMPlexCreatePointSF(DM dm, PetscSF migrationSF, PetscBool ownership, PetscSF *pointSF)
1508: {
1509:   PetscMPIInt        rank, size;
1510:   PetscInt           p, nroots, nleaves, idx, npointLeaves;
1511:   PetscInt          *pointLocal;
1512:   const PetscInt    *leaves;
1513:   const PetscSFNode *roots;
1514:   PetscSFNode       *rootNodes, *leafNodes, *pointRemote;
1515:   Vec                shifts;
1516:   const PetscInt     numShifts  = 13759;
1517:   const PetscScalar *shift      = NULL;
1518:   const PetscBool    shiftDebug = PETSC_FALSE;
1519:   PetscBool          balance;

1521:   PetscFunctionBegin;
1523:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1524:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
1525:   PetscCall(PetscLogEventBegin(DMPLEX_CreatePointSF, dm, 0, 0, 0));
1526:   PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), pointSF));
1527:   PetscCall(PetscSFSetFromOptions(*pointSF));
1528:   if (size == 1) PetscFunctionReturn(PETSC_SUCCESS);

1530:   PetscCall(DMPlexGetPartitionBalance(dm, &balance));
1531:   PetscCall(PetscSFGetGraph(migrationSF, &nroots, &nleaves, &leaves, &roots));
1532:   PetscCall(PetscMalloc2(nroots, &rootNodes, nleaves, &leafNodes));
1533:   if (ownership) {
1534:     MPI_Op       op;
1535:     MPI_Datatype datatype;
1536:     Petsc3Int   *rootVote = NULL, *leafVote = NULL;

1538:     /* If balancing, we compute a random cyclic shift of the rank for each remote point. That way, the max will evenly distribute among ranks. */
1539:     if (balance) {
1540:       PetscRandom r;

1542:       PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r));
1543:       PetscCall(PetscRandomSetInterval(r, 0, 2467 * size));
1544:       PetscCall(VecCreate(PETSC_COMM_SELF, &shifts));
1545:       PetscCall(VecSetSizes(shifts, numShifts, numShifts));
1546:       PetscCall(VecSetType(shifts, VECSTANDARD));
1547:       PetscCall(VecSetRandom(shifts, r));
1548:       PetscCall(PetscRandomDestroy(&r));
1549:       PetscCall(VecGetArrayRead(shifts, &shift));
1550:     }

1552:     PetscCall(PetscMalloc1(nroots, &rootVote));
1553:     PetscCall(PetscMalloc1(nleaves, &leafVote));
1554:     /* Point ownership vote: Process with highest rank owns shared points */
1555:     for (p = 0; p < nleaves; ++p) {
1556:       if (shiftDebug) {
1557:         PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Point %" PetscInt_FMT " RemotePoint %" PetscInt_FMT " Shift %" PetscInt_FMT " MyRank %" PetscInt_FMT "\n", rank, leaves ? leaves[p] : p, roots[p].index,
1558:                                           (PetscInt)PetscRealPart(shift[roots[p].index % numShifts]), (rank + (shift ? (PetscInt)PetscRealPart(shift[roots[p].index % numShifts]) : 0)) % size));
1559:       }
1560:       /* Either put in a bid or we know we own it */
1561:       leafVote[p].vote  = (rank + (shift ? (PetscInt)PetscRealPart(shift[roots[p].index % numShifts]) : 0)) % size;
1562:       leafVote[p].rank  = rank;
1563:       leafVote[p].index = p;
1564:     }
1565:     for (p = 0; p < nroots; p++) {
1566:       /* Root must not participate in the reduction, flag so that MAXLOC does not use */
1567:       rootVote[p].vote  = -3;
1568:       rootVote[p].rank  = -3;
1569:       rootVote[p].index = -3;
1570:     }
1571:     PetscCallMPI(MPI_Type_contiguous(3, MPIU_INT, &datatype));
1572:     PetscCallMPI(MPI_Type_commit(&datatype));
1573:     PetscCallMPI(MPI_Op_create(&MaxLocCarry, 1, &op));
1574:     PetscCall(PetscSFReduceBegin(migrationSF, datatype, leafVote, rootVote, op));
1575:     PetscCall(PetscSFReduceEnd(migrationSF, datatype, leafVote, rootVote, op));
1576:     PetscCallMPI(MPI_Op_free(&op));
1577:     PetscCallMPI(MPI_Type_free(&datatype));
1578:     for (p = 0; p < nroots; p++) {
1579:       rootNodes[p].rank  = rootVote[p].rank;
1580:       rootNodes[p].index = rootVote[p].index;
1581:     }
1582:     PetscCall(PetscFree(leafVote));
1583:     PetscCall(PetscFree(rootVote));
1584:   } else {
1585:     for (p = 0; p < nroots; p++) {
1586:       rootNodes[p].index = -1;
1587:       rootNodes[p].rank  = rank;
1588:     }
1589:     for (p = 0; p < nleaves; p++) {
1590:       /* Write new local id into old location */
1591:       if (roots[p].rank == rank) rootNodes[roots[p].index].index = leaves ? leaves[p] : p;
1592:     }
1593:   }
1594:   PetscCall(PetscSFBcastBegin(migrationSF, MPIU_SF_NODE, rootNodes, leafNodes, MPI_REPLACE));
1595:   PetscCall(PetscSFBcastEnd(migrationSF, MPIU_SF_NODE, rootNodes, leafNodes, MPI_REPLACE));

1597:   for (npointLeaves = 0, p = 0; p < nleaves; p++) {
1598:     if (leafNodes[p].rank != rank) npointLeaves++;
1599:   }
1600:   PetscCall(PetscMalloc1(npointLeaves, &pointLocal));
1601:   PetscCall(PetscMalloc1(npointLeaves, &pointRemote));
1602:   for (idx = 0, p = 0; p < nleaves; p++) {
1603:     if (leafNodes[p].rank != rank) {
1604:       /* Note that pointLocal is automatically sorted as it is sublist of 0, ..., nleaves-1 */
1605:       pointLocal[idx]  = p;
1606:       pointRemote[idx] = leafNodes[p];
1607:       idx++;
1608:     }
1609:   }
1610:   if (shift) {
1611:     PetscCall(VecRestoreArrayRead(shifts, &shift));
1612:     PetscCall(VecDestroy(&shifts));
1613:   }
1614:   if (shiftDebug) PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)dm), PETSC_STDOUT));
1615:   PetscCall(PetscSFSetGraph(*pointSF, nleaves, npointLeaves, pointLocal, PETSC_OWN_POINTER, pointRemote, PETSC_OWN_POINTER));
1616:   PetscCall(PetscFree2(rootNodes, leafNodes));
1617:   PetscCall(PetscLogEventEnd(DMPLEX_CreatePointSF, dm, 0, 0, 0));
1618:   if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, *pointSF, PETSC_FALSE));
1619:   PetscFunctionReturn(PETSC_SUCCESS);
1620: }

1622: /*@
1623:   DMPlexMigrate  - Migrates internal `DM` data over the supplied star forest

1625:   Collective

1627:   Input Parameters:
1628: + dm - The source `DMPLEX` object
1629: - sf - The star forest communication context describing the migration pattern

1631:   Output Parameter:
1632: . targetDM - The target `DMPLEX` object

1634:   Level: intermediate

1636: .seealso: `DMPLEX`, `PetscSF`, `DM`, `DMPlexDistribute()`, `DMPlexDistributeOverlap()`
1637: @*/
1638: PetscErrorCode DMPlexMigrate(DM dm, PetscSF sf, DM targetDM)
1639: {
1640:   MPI_Comm               comm;
1641:   PetscInt               dim, cdim, nroots;
1642:   PetscSF                sfPoint;
1643:   ISLocalToGlobalMapping ltogMigration;
1644:   ISLocalToGlobalMapping ltogOriginal = NULL;

1646:   PetscFunctionBegin;
1648:   PetscCall(PetscLogEventBegin(DMPLEX_Migrate, dm, 0, 0, 0));
1649:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1650:   PetscCall(DMGetDimension(dm, &dim));
1651:   PetscCall(DMSetDimension(targetDM, dim));
1652:   PetscCall(DMGetCoordinateDim(dm, &cdim));
1653:   PetscCall(DMSetCoordinateDim(targetDM, cdim));

1655:   /* Check for a one-to-all distribution pattern */
1656:   PetscCall(DMGetPointSF(dm, &sfPoint));
1657:   PetscCall(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL));
1658:   if (nroots >= 0) {
1659:     IS        isOriginal;
1660:     PetscInt  n, size, nleaves;
1661:     PetscInt *numbering_orig, *numbering_new;

1663:     /* Get the original point numbering */
1664:     PetscCall(DMPlexCreatePointNumbering(dm, &isOriginal));
1665:     PetscCall(ISLocalToGlobalMappingCreateIS(isOriginal, &ltogOriginal));
1666:     PetscCall(ISLocalToGlobalMappingGetSize(ltogOriginal, &size));
1667:     PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltogOriginal, (const PetscInt **)&numbering_orig));
1668:     /* Convert to positive global numbers */
1669:     for (n = 0; n < size; n++) {
1670:       if (numbering_orig[n] < 0) numbering_orig[n] = -(numbering_orig[n] + 1);
1671:     }
1672:     /* Derive the new local-to-global mapping from the old one */
1673:     PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, NULL, NULL));
1674:     PetscCall(PetscMalloc1(nleaves, &numbering_new));
1675:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, numbering_orig, numbering_new, MPI_REPLACE));
1676:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, numbering_orig, numbering_new, MPI_REPLACE));
1677:     PetscCall(ISLocalToGlobalMappingCreate(comm, 1, nleaves, numbering_new, PETSC_OWN_POINTER, &ltogMigration));
1678:     PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogOriginal, (const PetscInt **)&numbering_orig));
1679:     PetscCall(ISDestroy(&isOriginal));
1680:   } else {
1681:     /* One-to-all distribution pattern: We can derive LToG from SF */
1682:     PetscCall(ISLocalToGlobalMappingCreateSF(sf, 0, &ltogMigration));
1683:   }
1684:   PetscCall(PetscObjectSetName((PetscObject)ltogMigration, "Point renumbering for DM migration"));
1685:   PetscCall(ISLocalToGlobalMappingViewFromOptions(ltogMigration, (PetscObject)dm, "-partition_view"));
1686:   /* Migrate DM data to target DM */
1687:   PetscCall(DMPlexDistributeCones(dm, sf, ltogOriginal, ltogMigration, targetDM));
1688:   PetscCall(DMPlexDistributeLabels(dm, sf, targetDM));
1689:   PetscCall(DMPlexDistributeCoordinates(dm, sf, targetDM));
1690:   PetscCall(DMPlexDistributeSetupTree(dm, sf, ltogOriginal, ltogMigration, targetDM));
1691:   PetscCall(ISLocalToGlobalMappingDestroy(&ltogOriginal));
1692:   PetscCall(ISLocalToGlobalMappingDestroy(&ltogMigration));
1693:   PetscCall(PetscLogEventEnd(DMPLEX_Migrate, dm, 0, 0, 0));
1694:   PetscFunctionReturn(PETSC_SUCCESS);
1695: }

1697: /*@
1698:   DMPlexRemapMigrationSF - Rewrite the distribution SF to account for overlap

1700:   Collective

1702:   Input Parameters:
1703: + sfOverlap   - The `PetscSF` object just for the overlap
1704: - sfMigration - The original distribution `PetscSF` object

1706:   Output Parameters:
1707: . sfMigrationNew - A rewritten `PetscSF` object that incorporates the overlap

1709:   Level: developer

1711: .seealso: `DMPLEX`, `DM`, `DMPlexDistribute()`, `DMPlexDistributeOverlap()`, `DMPlexGetOverlap()`
1712: @*/
1713: PetscErrorCode DMPlexRemapMigrationSF(PetscSF sfOverlap, PetscSF sfMigration, PetscSF *sfMigrationNew)
1714: {
1715:   PetscSFNode       *newRemote, *permRemote = NULL;
1716:   const PetscInt    *oldLeaves;
1717:   const PetscSFNode *oldRemote;
1718:   PetscInt           nroots, nleaves, noldleaves;

1720:   PetscFunctionBegin;
1721:   PetscCall(PetscSFGetGraph(sfMigration, &nroots, &noldleaves, &oldLeaves, &oldRemote));
1722:   PetscCall(PetscSFGetGraph(sfOverlap, NULL, &nleaves, NULL, NULL));
1723:   PetscCall(PetscMalloc1(nleaves, &newRemote));
1724:   /* oldRemote: original root point mapping to original leaf point
1725:      newRemote: original leaf point mapping to overlapped leaf point */
1726:   if (oldLeaves) {
1727:     /* After stratification, the migration remotes may not be in root (canonical) order, so we reorder using the leaf numbering */
1728:     PetscCall(PetscMalloc1(noldleaves, &permRemote));
1729:     for (PetscInt l = 0; l < noldleaves; ++l) permRemote[oldLeaves[l]] = oldRemote[l];
1730:     oldRemote = permRemote;
1731:   }
1732:   PetscCall(PetscSFBcastBegin(sfOverlap, MPIU_SF_NODE, oldRemote, newRemote, MPI_REPLACE));
1733:   PetscCall(PetscSFBcastEnd(sfOverlap, MPIU_SF_NODE, oldRemote, newRemote, MPI_REPLACE));
1734:   PetscCall(PetscFree(permRemote));
1735:   PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)sfOverlap), sfMigrationNew));
1736:   PetscCall(PetscSFSetGraph(*sfMigrationNew, nroots, nleaves, NULL, PETSC_OWN_POINTER, newRemote, PETSC_OWN_POINTER));
1737:   PetscFunctionReturn(PETSC_SUCCESS);
1738: }

1740: /* For DG-like methods, the code below is equivalent (but faster) than calling
1741:    DMPlexCreateClosureIndex(dm,section) */
1742: static PetscErrorCode DMPlexCreateClosureIndex_CELL(DM dm, PetscSection section)
1743: {
1744:   PetscSection clSection;
1745:   IS           clPoints;
1746:   PetscInt     pStart, pEnd, point;
1747:   PetscInt    *closure, pos = 0;

1749:   PetscFunctionBegin;
1750:   if (!section) PetscCall(DMGetLocalSection(dm, &section));
1751:   PetscCall(DMPlexGetHeightStratum(dm, 0, &pStart, &pEnd));
1752:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &clSection));
1753:   PetscCall(PetscSectionSetChart(clSection, pStart, pEnd));
1754:   PetscCall(PetscMalloc1((2 * (pEnd - pStart)), &closure));
1755:   for (point = pStart; point < pEnd; point++) {
1756:     PetscCall(PetscSectionSetDof(clSection, point, 2));
1757:     closure[pos++] = point; /* point */
1758:     closure[pos++] = 0;     /* orientation */
1759:   }
1760:   PetscCall(PetscSectionSetUp(clSection));
1761:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 2 * (pEnd - pStart), closure, PETSC_OWN_POINTER, &clPoints));
1762:   PetscCall(PetscSectionSetClosureIndex(section, (PetscObject)dm, clSection, clPoints));
1763:   PetscCall(PetscSectionDestroy(&clSection));
1764:   PetscCall(ISDestroy(&clPoints));
1765:   PetscFunctionReturn(PETSC_SUCCESS);
1766: }

1768: static PetscErrorCode DMPlexDistribute_Multistage(DM dm, PetscInt overlap, PetscSF *sf, DM *dmParallel)
1769: {
1770:   MPI_Comm         comm = PetscObjectComm((PetscObject)dm);
1771:   PetscPartitioner part;
1772:   PetscBool        balance, printHeader;
1773:   PetscInt         nl = 0;

1775:   PetscFunctionBegin;
1776:   if (sf) *sf = NULL;
1777:   *dmParallel = NULL;

1779:   PetscCall(DMPlexGetPartitioner(dm, &part));
1780:   printHeader = part->printHeader;
1781:   PetscCall(DMPlexGetPartitionBalance(dm, &balance));
1782:   PetscCall(PetscPartitionerSetUp(part));
1783:   PetscCall(PetscLogEventBegin(DMPLEX_DistributeMultistage, dm, 0, 0, 0));
1784:   PetscCall(PetscPartitionerMultistageGetStages_Multistage(part, &nl, NULL));
1785:   for (PetscInt l = 0; l < nl; l++) {
1786:     PetscInt ovl = (l < nl - 1) ? 0 : overlap;
1787:     PetscSF  sfDist;
1788:     DM       dmDist;

1790:     PetscCall(DMPlexSetPartitionBalance(dm, balance));
1791:     PetscCall(DMViewFromOptions(dm, (PetscObject)part, "-petscpartitioner_multistage_dm_view"));
1792:     PetscCall(PetscPartitionerMultistageSetStage_Multistage(part, l, (PetscObject)dm));
1793:     PetscCall(DMPlexSetPartitioner(dm, part));
1794:     PetscCall(DMPlexDistribute(dm, ovl, &sfDist, &dmDist));
1795:     PetscCheck(dmDist, comm, PETSC_ERR_PLIB, "No distributed DM generated (stage %" PetscInt_FMT ")", l);
1796:     PetscCheck(sfDist, comm, PETSC_ERR_PLIB, "No SF generated (stage %" PetscInt_FMT ")", l);
1797:     part->printHeader = PETSC_FALSE;

1799:     /* Propagate cell weights to the next level (if any, and if not the final dm) */
1800:     if (part->usevwgt && dm->localSection && l != nl - 1) {
1801:       PetscSection oldSection, newSection;

1803:       PetscCall(DMGetLocalSection(dm, &oldSection));
1804:       PetscCall(DMGetLocalSection(dmDist, &newSection));
1805:       PetscCall(PetscSFDistributeSection(sfDist, oldSection, NULL, newSection));
1806:       PetscCall(DMPlexCreateClosureIndex_CELL(dmDist, newSection));
1807:     }
1808:     if (!sf) PetscCall(PetscSFDestroy(&sfDist));
1809:     if (l > 0) PetscCall(DMDestroy(&dm));

1811:     if (sf && *sf) {
1812:       PetscSF sfA = *sf, sfB = sfDist;
1813:       PetscCall(PetscSFCompose(sfA, sfB, &sfDist));
1814:       PetscCall(PetscSFDestroy(&sfA));
1815:       PetscCall(PetscSFDestroy(&sfB));
1816:     }

1818:     if (sf) *sf = sfDist;
1819:     dm = *dmParallel = dmDist;
1820:   }
1821:   PetscCall(PetscPartitionerMultistageSetStage_Multistage(part, 0, NULL)); /* reset */
1822:   PetscCall(PetscLogEventEnd(DMPLEX_DistributeMultistage, dm, 0, 0, 0));
1823:   part->printHeader = printHeader;
1824:   PetscFunctionReturn(PETSC_SUCCESS);
1825: }

1827: /*@
1828:   DMPlexDistribute - Distributes the mesh and any associated sections.

1830:   Collective

1832:   Input Parameters:
1833: + dm      - The original `DMPLEX` object
1834: - overlap - The overlap of partitions, 0 is the default

1836:   Output Parameters:
1837: + sf         - The `PetscSF` used for point distribution, or `NULL` if not needed
1838: - dmParallel - The distributed `DMPLEX` object

1840:   Level: intermediate

1842:   Note:
1843:   If the mesh was not distributed, the output `dmParallel` will be `NULL`.

1845:   The user can control the definition of adjacency for the mesh using `DMSetAdjacency()`. They should choose the combination appropriate for the function
1846:   representation on the mesh.

1848: .seealso: `DMPLEX`, `DM`, `DMPlexCreate()`, `DMSetAdjacency()`, `DMPlexGetOverlap()`
1849: @*/
1850: PetscErrorCode DMPlexDistribute(DM dm, PetscInt overlap, PeOp PetscSF *sf, DM *dmParallel)
1851: {
1852:   MPI_Comm         comm;
1853:   PetscPartitioner partitioner;
1854:   IS               cellPart;
1855:   PetscSection     cellPartSection;
1856:   DM               dmCoord;
1857:   DMLabel          lblPartition, lblMigration;
1858:   PetscSF          sfMigration, sfStratified, sfPoint;
1859:   PetscBool        flg, balance, isms;
1860:   PetscMPIInt      rank, size;

1862:   PetscFunctionBegin;
1865:   if (sf) PetscAssertPointer(sf, 3);
1866:   PetscAssertPointer(dmParallel, 4);

1868:   if (sf) *sf = NULL;
1869:   *dmParallel = NULL;
1870:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1871:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1872:   PetscCallMPI(MPI_Comm_size(comm, &size));
1873:   if (size == 1) PetscFunctionReturn(PETSC_SUCCESS);

1875:   /* Handle multistage partitioner */
1876:   PetscCall(DMPlexGetPartitioner(dm, &partitioner));
1877:   PetscCall(PetscObjectTypeCompare((PetscObject)partitioner, PETSCPARTITIONERMULTISTAGE, &isms));
1878:   if (isms) {
1879:     PetscObject stagedm;

1881:     PetscCall(PetscPartitionerMultistageGetStage_Multistage(partitioner, NULL, &stagedm));
1882:     if (!stagedm) { /* No stage dm present, start the multistage algorithm */
1883:       PetscCall(DMPlexDistribute_Multistage(dm, overlap, sf, dmParallel));
1884:       PetscFunctionReturn(PETSC_SUCCESS);
1885:     }
1886:   }
1887:   PetscCall(PetscLogEventBegin(DMPLEX_Distribute, dm, 0, 0, 0));
1888:   /* Create cell partition */
1889:   PetscCall(PetscLogEventBegin(DMPLEX_Partition, dm, 0, 0, 0));
1890:   PetscCall(PetscSectionCreate(comm, &cellPartSection));
1891:   PetscCall(PetscPartitionerDMPlexPartition(partitioner, dm, NULL, cellPartSection, &cellPart));
1892:   PetscCall(PetscLogEventBegin(DMPLEX_PartSelf, dm, 0, 0, 0));
1893:   {
1894:     /* Convert partition to DMLabel */
1895:     IS              is;
1896:     PetscHSetI      ht;
1897:     const PetscInt *points;
1898:     PetscInt       *iranks;
1899:     PetscInt        pStart, pEnd, proc, npoints, poff = 0, nranks;

1901:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Point Partition", &lblPartition));
1902:     /* Preallocate strata */
1903:     PetscCall(PetscHSetICreate(&ht));
1904:     PetscCall(PetscSectionGetChart(cellPartSection, &pStart, &pEnd));
1905:     for (proc = pStart; proc < pEnd; proc++) {
1906:       PetscCall(PetscSectionGetDof(cellPartSection, proc, &npoints));
1907:       if (npoints) PetscCall(PetscHSetIAdd(ht, proc));
1908:     }
1909:     PetscCall(PetscHSetIGetSize(ht, &nranks));
1910:     PetscCall(PetscMalloc1(nranks, &iranks));
1911:     PetscCall(PetscHSetIGetElems(ht, &poff, iranks));
1912:     PetscCall(PetscHSetIDestroy(&ht));
1913:     PetscCall(DMLabelAddStrata(lblPartition, nranks, iranks));
1914:     PetscCall(PetscFree(iranks));
1915:     /* Inline DMPlexPartitionLabelClosure() */
1916:     PetscCall(ISGetIndices(cellPart, &points));
1917:     PetscCall(PetscSectionGetChart(cellPartSection, &pStart, &pEnd));
1918:     for (proc = pStart; proc < pEnd; proc++) {
1919:       PetscCall(PetscSectionGetDof(cellPartSection, proc, &npoints));
1920:       if (!npoints) continue;
1921:       PetscCall(PetscSectionGetOffset(cellPartSection, proc, &poff));
1922:       PetscCall(DMPlexClosurePoints_Private(dm, npoints, points + poff, &is));
1923:       PetscCall(DMLabelSetStratumIS(lblPartition, proc, is));
1924:       PetscCall(ISDestroy(&is));
1925:     }
1926:     PetscCall(ISRestoreIndices(cellPart, &points));
1927:   }
1928:   PetscCall(PetscLogEventEnd(DMPLEX_PartSelf, dm, 0, 0, 0));

1930:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Point migration", &lblMigration));
1931:   PetscCall(DMPlexPartitionLabelInvert(dm, lblPartition, NULL, lblMigration));
1932:   PetscCall(DMPlexPartitionLabelCreateSF(dm, lblMigration, PETSC_TRUE, &sfMigration));
1933:   PetscCall(DMPlexStratifyMigrationSF(dm, sfMigration, &sfStratified));
1934:   PetscCall(PetscSFDestroy(&sfMigration));
1935:   sfMigration = sfStratified;
1936:   PetscCall(PetscSFSetUp(sfMigration));
1937:   PetscCall(PetscLogEventEnd(DMPLEX_Partition, dm, 0, 0, 0));
1938:   PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-partition_view", &flg));
1939:   if (flg) {
1940:     PetscCall(DMLabelView(lblPartition, PETSC_VIEWER_STDOUT_(comm)));
1941:     PetscCall(PetscSFView(sfMigration, PETSC_VIEWER_STDOUT_(comm)));
1942:   }

1944:   /* Create non-overlapping parallel DM and migrate internal data */
1945:   PetscCall(DMPlexCreate(comm, dmParallel));
1946:   PetscCall(PetscObjectSetName((PetscObject)*dmParallel, "Parallel Mesh"));
1947:   PetscCall(DMPlexMigrate(dm, sfMigration, *dmParallel));

1949:   /* Build the point SF without overlap */
1950:   PetscCall(DMPlexGetPartitionBalance(dm, &balance));
1951:   PetscCall(DMPlexSetPartitionBalance(*dmParallel, balance));
1952:   PetscCall(DMPlexCreatePointSF(*dmParallel, sfMigration, PETSC_TRUE, &sfPoint));
1953:   PetscCall(DMSetPointSF(*dmParallel, sfPoint));
1954:   PetscCall(DMPlexMigrateIsoperiodicFaceSF_Internal(dm, *dmParallel, sfMigration));
1955:   PetscCall(DMGetCoordinateDM(*dmParallel, &dmCoord));
1956:   if (dmCoord) PetscCall(DMSetPointSF(dmCoord, sfPoint));
1957:   if (flg) PetscCall(PetscSFView(sfPoint, NULL));

1959:   if (overlap > 0) {
1960:     DM      dmOverlap;
1961:     PetscSF sfOverlap, sfOverlapPoint;

1963:     /* Add the partition overlap to the distributed DM */
1964:     PetscCall(DMPlexDistributeOverlap(*dmParallel, overlap, &sfOverlap, &dmOverlap));
1965:     PetscCall(DMDestroy(dmParallel));
1966:     *dmParallel = dmOverlap;
1967:     if (flg) {
1968:       PetscCall(PetscPrintf(comm, "Overlap Migration SF:\n"));
1969:       PetscCall(PetscSFView(sfOverlap, NULL));
1970:     }

1972:     /* Re-map the migration SF to establish the full migration pattern */
1973:     PetscCall(DMPlexRemapMigrationSF(sfOverlap, sfMigration, &sfOverlapPoint));
1974:     PetscCall(PetscSFDestroy(&sfOverlap));
1975:     PetscCall(PetscSFDestroy(&sfMigration));
1976:     sfMigration = sfOverlapPoint;
1977:   }
1978:   /* Cleanup Partition */
1979:   PetscCall(DMLabelDestroy(&lblPartition));
1980:   PetscCall(DMLabelDestroy(&lblMigration));
1981:   PetscCall(PetscSectionDestroy(&cellPartSection));
1982:   PetscCall(ISDestroy(&cellPart));
1983:   PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_FALSE, *dmParallel));
1984:   // Create sfNatural, need discretization information
1985:   PetscCall(DMCopyDisc(dm, *dmParallel));
1986:   if (dm->localSection) {
1987:     PetscSection psection;
1988:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &psection));
1989:     PetscCall(PetscSFDistributeSection(sfMigration, dm->localSection, NULL, psection));
1990:     PetscCall(DMSetLocalSection(*dmParallel, psection));
1991:     PetscCall(PetscSectionDestroy(&psection));
1992:   }
1993:   if (dm->useNatural) {
1994:     PetscSection section;

1996:     PetscCall(DMSetUseNatural(*dmParallel, PETSC_TRUE));
1997:     PetscCall(DMGetLocalSection(dm, &section));

1999:     /* If DM has useNatural = PETSC_TRUE, but no sfNatural is set, the DM's Section is considered natural */
2000:     /* sfMigration and sfNatural are respectively the point and dofs SFs mapping to this natural DM */
2001:     /* Compose with a previous sfNatural if present */
2002:     if (dm->sfNatural) {
2003:       PetscCall(DMPlexMigrateGlobalToNaturalSF(dm, *dmParallel, dm->sfNatural, sfMigration, &(*dmParallel)->sfNatural));
2004:     } else {
2005:       PetscCall(DMPlexCreateGlobalToNaturalSF(*dmParallel, section, sfMigration, &(*dmParallel)->sfNatural));
2006:     }
2007:     /* Compose with a previous sfMigration if present */
2008:     if (dm->sfMigration) {
2009:       PetscSF naturalPointSF;

2011:       PetscCall(PetscSFCompose(dm->sfMigration, sfMigration, &naturalPointSF));
2012:       PetscCall(PetscSFDestroy(&sfMigration));
2013:       sfMigration = naturalPointSF;
2014:     }
2015:   }
2016:   /* Cleanup */
2017:   if (sf) {
2018:     *sf = sfMigration;
2019:   } else PetscCall(PetscSFDestroy(&sfMigration));
2020:   PetscCall(PetscSFDestroy(&sfPoint));
2021:   PetscCall(PetscLogEventEnd(DMPLEX_Distribute, dm, 0, 0, 0));
2022:   PetscFunctionReturn(PETSC_SUCCESS);
2023: }

2025: PetscErrorCode DMPlexDistributeOverlap_Internal(DM dm, PetscInt overlap, MPI_Comm newcomm, const char *ovlboundary, PetscSF *sf, DM *dmOverlap)
2026: {
2027:   DM_Plex     *mesh = (DM_Plex *)dm->data;
2028:   MPI_Comm     comm;
2029:   PetscMPIInt  size, rank;
2030:   PetscSection rootSection, leafSection;
2031:   IS           rootrank, leafrank;
2032:   DM           dmCoord;
2033:   DMLabel      lblOverlap;
2034:   PetscSF      sfOverlap, sfStratified, sfPoint;
2035:   PetscBool    clear_ovlboundary;

2037:   PetscFunctionBegin;
2038:   if (sf) *sf = NULL;
2039:   *dmOverlap = NULL;
2040:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2041:   PetscCallMPI(MPI_Comm_size(comm, &size));
2042:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
2043:   if (size == 1) PetscFunctionReturn(PETSC_SUCCESS);
2044:   {
2045:     // We need to get options for the _already_distributed mesh, so it must be done here
2046:     PetscInt    overlap;
2047:     const char *prefix;
2048:     char        oldPrefix[PETSC_MAX_PATH_LEN];

2050:     oldPrefix[0] = '\0';
2051:     PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
2052:     PetscCall(PetscStrncpy(oldPrefix, prefix, sizeof(oldPrefix)));
2053:     PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)dm, "dist_"));
2054:     PetscCall(DMPlexGetOverlap(dm, &overlap));
2055:     PetscObjectOptionsBegin((PetscObject)dm);
2056:     PetscCall(DMSetFromOptions_Overlap_Plex(dm, PetscOptionsObject, &overlap));
2057:     PetscOptionsEnd();
2058:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, oldPrefix[0] == '\0' ? NULL : oldPrefix));
2059:   }
2060:   if (ovlboundary) {
2061:     PetscBool flg;
2062:     PetscCall(DMHasLabel(dm, ovlboundary, &flg));
2063:     if (!flg) {
2064:       DMLabel label;

2066:       PetscCall(DMCreateLabel(dm, ovlboundary));
2067:       PetscCall(DMGetLabel(dm, ovlboundary, &label));
2068:       PetscCall(DMPlexMarkBoundaryFaces(dm, 1, label));
2069:       clear_ovlboundary = PETSC_TRUE;
2070:     }
2071:   }
2072:   PetscCall(PetscLogEventBegin(DMPLEX_DistributeOverlap, dm, 0, 0, 0));
2073:   /* Compute point overlap with neighbouring processes on the distributed DM */
2074:   PetscCall(PetscLogEventBegin(DMPLEX_Partition, dm, 0, 0, 0));
2075:   PetscCall(PetscSectionCreate(newcomm, &rootSection));
2076:   PetscCall(PetscSectionCreate(newcomm, &leafSection));
2077:   PetscCall(DMPlexDistributeOwnership(dm, rootSection, &rootrank, leafSection, &leafrank));
2078:   if (mesh->numOvLabels) PetscCall(DMPlexCreateOverlapLabelFromLabels(dm, mesh->numOvLabels, mesh->ovLabels, mesh->ovValues, mesh->numOvExLabels, mesh->ovExLabels, mesh->ovExValues, rootSection, rootrank, leafSection, leafrank, &lblOverlap));
2079:   else PetscCall(DMPlexCreateOverlapLabel(dm, overlap, rootSection, rootrank, leafSection, leafrank, &lblOverlap));
2080:   /* Convert overlap label to stratified migration SF */
2081:   PetscCall(DMPlexPartitionLabelCreateSF(dm, lblOverlap, PETSC_FALSE, &sfOverlap));
2082:   PetscCall(DMPlexStratifyMigrationSF(dm, sfOverlap, &sfStratified));
2083:   PetscCall(PetscSFDestroy(&sfOverlap));
2084:   sfOverlap = sfStratified;
2085:   PetscCall(PetscObjectSetName((PetscObject)sfOverlap, "Overlap SF"));
2086:   PetscCall(PetscSFSetFromOptions(sfOverlap));

2088:   PetscCall(PetscSectionDestroy(&rootSection));
2089:   PetscCall(PetscSectionDestroy(&leafSection));
2090:   PetscCall(ISDestroy(&rootrank));
2091:   PetscCall(ISDestroy(&leafrank));
2092:   PetscCall(PetscLogEventEnd(DMPLEX_Partition, dm, 0, 0, 0));

2094:   /* Build the overlapping DM */
2095:   PetscCall(DMPlexCreate(newcomm, dmOverlap));
2096:   PetscCall(PetscObjectSetName((PetscObject)*dmOverlap, "Parallel Mesh"));
2097:   PetscCall(DMPlexMigrate(dm, sfOverlap, *dmOverlap));
2098:   /* Store the overlap in the new DM */
2099:   PetscCall(DMPlexSetOverlap_Plex(*dmOverlap, dm, overlap));
2100:   /* Build the new point SF */
2101:   PetscCall(DMPlexCreatePointSF(*dmOverlap, sfOverlap, PETSC_FALSE, &sfPoint));
2102:   PetscCall(DMSetPointSF(*dmOverlap, sfPoint));
2103:   PetscCall(DMGetCoordinateDM(*dmOverlap, &dmCoord));
2104:   if (dmCoord) PetscCall(DMSetPointSF(dmCoord, sfPoint));
2105:   PetscCall(DMGetCellCoordinateDM(*dmOverlap, &dmCoord));
2106:   if (dmCoord) PetscCall(DMSetPointSF(dmCoord, sfPoint));
2107:   PetscCall(PetscSFDestroy(&sfPoint));
2108:   PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_FALSE, *dmOverlap));
2109:   /* TODO: labels stored inside the DS for regions should be handled here */
2110:   PetscCall(DMCopyDisc(dm, *dmOverlap));
2111:   /* Cleanup overlap partition */
2112:   PetscCall(DMLabelDestroy(&lblOverlap));
2113:   if (sf) *sf = sfOverlap;
2114:   else PetscCall(PetscSFDestroy(&sfOverlap));
2115:   if (ovlboundary) {
2116:     DMLabel   label;
2117:     PetscBool flg;

2119:     if (clear_ovlboundary) PetscCall(DMRemoveLabel(dm, ovlboundary, NULL));
2120:     PetscCall(DMHasLabel(*dmOverlap, ovlboundary, &flg));
2121:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing %s label", ovlboundary);
2122:     PetscCall(DMGetLabel(*dmOverlap, ovlboundary, &label));
2123:     PetscCall(DMPlexMarkBoundaryFaces_Internal(*dmOverlap, 1, 0, label, PETSC_TRUE));
2124:   }
2125:   PetscCall(PetscLogEventEnd(DMPLEX_DistributeOverlap, dm, 0, 0, 0));
2126:   PetscFunctionReturn(PETSC_SUCCESS);
2127: }

2129: /*@
2130:   DMPlexDistributeOverlap - Add partition overlap to a distributed non-overlapping `DM`.

2132:   Collective

2134:   Input Parameters:
2135: + dm      - The non-overlapping distributed `DMPLEX` object
2136: - overlap - The overlap of partitions (the same on all ranks)

2138:   Output Parameters:
2139: + sf        - The `PetscSF` used for point distribution, or pass `NULL` if not needed
2140: - dmOverlap - The overlapping distributed `DMPLEX` object

2142:   Options Database Keys:
2143: + -dm_plex_overlap_labels <name1,name2,...> - List of overlap label names
2144: . -dm_plex_overlap_values <int1,int2,...>   - List of overlap label values
2145: . -dm_plex_overlap_exclude_label <name>     - Label used to exclude points from overlap
2146: - -dm_plex_overlap_exclude_value <int>      - Label value used to exclude points from overlap

2148:   Level: advanced

2150:   Notes:
2151:   If the mesh was not distributed, the return value is `NULL`.

2153:   The user can control the definition of adjacency for the mesh using `DMSetAdjacency()`. They should choose the combination appropriate for the function
2154:   representation on the mesh.

2156: .seealso: `DMPLEX`, `PetscSF`, `DM`, `DMPlexCreate()`, `DMSetAdjacency()`, `DMPlexDistribute()`, `DMPlexCreateOverlapLabel()`, `DMPlexGetOverlap()`
2157: @*/
2158: PetscErrorCode DMPlexDistributeOverlap(DM dm, PetscInt overlap, PeOp PetscSF *sf, DM *dmOverlap)
2159: {
2160:   PetscFunctionBegin;
2163:   if (sf) PetscAssertPointer(sf, 3);
2164:   PetscAssertPointer(dmOverlap, 4);
2165:   PetscCall(DMPlexDistributeOverlap_Internal(dm, overlap, PetscObjectComm((PetscObject)dm), NULL, sf, dmOverlap));
2166:   PetscFunctionReturn(PETSC_SUCCESS);
2167: }

2169: PetscErrorCode DMPlexGetOverlap_Plex(DM dm, PetscInt *overlap)
2170: {
2171:   DM_Plex *mesh = (DM_Plex *)dm->data;

2173:   PetscFunctionBegin;
2174:   *overlap = mesh->overlap;
2175:   PetscFunctionReturn(PETSC_SUCCESS);
2176: }

2178: PetscErrorCode DMPlexSetOverlap_Plex(DM dm, DM dmSrc, PetscInt overlap)
2179: {
2180:   DM_Plex *mesh    = NULL;
2181:   DM_Plex *meshSrc = NULL;

2183:   PetscFunctionBegin;
2185:   mesh = (DM_Plex *)dm->data;
2186:   if (dmSrc) {
2188:     meshSrc       = (DM_Plex *)dmSrc->data;
2189:     mesh->overlap = meshSrc->overlap;
2190:   } else {
2191:     mesh->overlap = 0;
2192:   }
2193:   mesh->overlap += overlap;
2194:   PetscFunctionReturn(PETSC_SUCCESS);
2195: }

2197: /*@
2198:   DMPlexGetOverlap - Get the width of the cell overlap

2200:   Not Collective

2202:   Input Parameter:
2203: . dm - The `DM`

2205:   Output Parameter:
2206: . overlap - the width of the cell overlap

2208:   Level: intermediate

2210: .seealso: `DMPLEX`, `DMPlexSetOverlap()`, `DMPlexDistribute()`
2211: @*/
2212: PetscErrorCode DMPlexGetOverlap(DM dm, PetscInt *overlap)
2213: {
2214:   PetscFunctionBegin;
2216:   PetscAssertPointer(overlap, 2);
2217:   PetscUseMethod(dm, "DMPlexGetOverlap_C", (DM, PetscInt *), (dm, overlap));
2218:   PetscFunctionReturn(PETSC_SUCCESS);
2219: }

2221: /*@
2222:   DMPlexSetOverlap - Set the width of the cell overlap

2224:   Logically Collective

2226:   Input Parameters:
2227: + dm      - The `DM`
2228: . dmSrc   - The `DM` that produced this one, or `NULL`
2229: - overlap - the width of the cell overlap

2231:   Level: intermediate

2233:   Note:
2234:   The overlap from `dmSrc` is added to `dm`

2236: .seealso: `DMPLEX`, `DMPlexGetOverlap()`, `DMPlexDistribute()`
2237: @*/
2238: PetscErrorCode DMPlexSetOverlap(DM dm, DM dmSrc, PetscInt overlap)
2239: {
2240:   PetscFunctionBegin;
2243:   PetscTryMethod(dm, "DMPlexSetOverlap_C", (DM, DM, PetscInt), (dm, dmSrc, overlap));
2244:   PetscFunctionReturn(PETSC_SUCCESS);
2245: }

2247: PetscErrorCode DMPlexDistributeSetDefault_Plex(DM dm, PetscBool dist)
2248: {
2249:   DM_Plex *mesh = (DM_Plex *)dm->data;

2251:   PetscFunctionBegin;
2252:   mesh->distDefault = dist;
2253:   PetscFunctionReturn(PETSC_SUCCESS);
2254: }

2256: /*@
2257:   DMPlexDistributeSetDefault - Set flag indicating whether the `DM` should be distributed by default

2259:   Logically Collective

2261:   Input Parameters:
2262: + dm   - The `DM`
2263: - dist - Flag for distribution

2265:   Level: intermediate

2267: .seealso: `DMPLEX`, `DMPlexDistributeGetDefault()`, `DMPlexDistribute()`
2268: @*/
2269: PetscErrorCode DMPlexDistributeSetDefault(DM dm, PetscBool dist)
2270: {
2271:   PetscFunctionBegin;
2274:   PetscTryMethod(dm, "DMPlexDistributeSetDefault_C", (DM, PetscBool), (dm, dist));
2275:   PetscFunctionReturn(PETSC_SUCCESS);
2276: }

2278: PetscErrorCode DMPlexDistributeGetDefault_Plex(DM dm, PetscBool *dist)
2279: {
2280:   DM_Plex *mesh = (DM_Plex *)dm->data;

2282:   PetscFunctionBegin;
2283:   *dist = mesh->distDefault;
2284:   PetscFunctionReturn(PETSC_SUCCESS);
2285: }

2287: /*@
2288:   DMPlexDistributeGetDefault - Get flag indicating whether the `DM` should be distributed by default

2290:   Not Collective

2292:   Input Parameter:
2293: . dm - The `DM`

2295:   Output Parameter:
2296: . dist - Flag for distribution

2298:   Level: intermediate

2300: .seealso: `DMPLEX`, `DM`, `DMPlexDistributeSetDefault()`, `DMPlexDistribute()`
2301: @*/
2302: PetscErrorCode DMPlexDistributeGetDefault(DM dm, PetscBool *dist)
2303: {
2304:   PetscFunctionBegin;
2306:   PetscAssertPointer(dist, 2);
2307:   PetscUseMethod(dm, "DMPlexDistributeGetDefault_C", (DM, PetscBool *), (dm, dist));
2308:   PetscFunctionReturn(PETSC_SUCCESS);
2309: }

2311: /*@
2312:   DMPlexGetGatherDM - Get a copy of the `DMPLEX` that gathers all points on the
2313:   root process of the original's communicator.

2315:   Collective

2317:   Input Parameter:
2318: . dm - the original `DMPLEX` object

2320:   Output Parameters:
2321: + sf         - the `PetscSF` used for point distribution (optional)
2322: - gatherMesh - the gathered `DM` object, or `NULL`

2324:   Level: intermediate

2326: .seealso: `DMPLEX`, `DM`, `PetscSF`, `DMPlexDistribute()`, `DMPlexGetRedundantDM()`
2327: @*/
2328: PetscErrorCode DMPlexGetGatherDM(DM dm, PetscSF *sf, PeOp DM *gatherMesh)
2329: {
2330:   MPI_Comm         comm;
2331:   PetscMPIInt      size;
2332:   PetscPartitioner oldPart, gatherPart;

2334:   PetscFunctionBegin;
2336:   PetscAssertPointer(gatherMesh, 3);
2337:   *gatherMesh = NULL;
2338:   if (sf) *sf = NULL;
2339:   comm = PetscObjectComm((PetscObject)dm);
2340:   PetscCallMPI(MPI_Comm_size(comm, &size));
2341:   if (size == 1) PetscFunctionReturn(PETSC_SUCCESS);
2342:   PetscCall(DMPlexGetPartitioner(dm, &oldPart));
2343:   PetscCall(PetscObjectReference((PetscObject)oldPart));
2344:   PetscCall(PetscPartitionerCreate(comm, &gatherPart));
2345:   PetscCall(PetscPartitionerSetType(gatherPart, PETSCPARTITIONERGATHER));
2346:   PetscCall(DMPlexSetPartitioner(dm, gatherPart));
2347:   PetscCall(DMPlexDistribute(dm, 0, sf, gatherMesh));

2349:   PetscCall(DMPlexSetPartitioner(dm, oldPart));
2350:   PetscCall(PetscPartitionerDestroy(&gatherPart));
2351:   PetscCall(PetscPartitionerDestroy(&oldPart));
2352:   PetscFunctionReturn(PETSC_SUCCESS);
2353: }

2355: /*@
2356:   DMPlexGetRedundantDM - Get a copy of the `DMPLEX` that is completely copied on each process.

2358:   Collective

2360:   Input Parameter:
2361: . dm - the original `DMPLEX` object

2363:   Output Parameters:
2364: + sf            - the `PetscSF` used for point distribution (optional)
2365: - redundantMesh - the redundant `DM` object, or `NULL`

2367:   Level: intermediate

2369: .seealso: `DMPLEX`, `DMPlexDistribute()`, `DMPlexGetGatherDM()`
2370: @*/
2371: PetscErrorCode DMPlexGetRedundantDM(DM dm, PetscSF *sf, PeOp DM *redundantMesh)
2372: {
2373:   MPI_Comm     comm;
2374:   PetscMPIInt  size, rank;
2375:   PetscInt     pStart, pEnd, p;
2376:   PetscInt     numPoints = -1;
2377:   PetscSF      migrationSF, sfPoint, gatherSF;
2378:   DM           gatherDM, dmCoord;
2379:   PetscSFNode *points;

2381:   PetscFunctionBegin;
2383:   PetscAssertPointer(redundantMesh, 3);
2384:   *redundantMesh = NULL;
2385:   comm           = PetscObjectComm((PetscObject)dm);
2386:   PetscCallMPI(MPI_Comm_size(comm, &size));
2387:   if (size == 1) {
2388:     PetscCall(PetscObjectReference((PetscObject)dm));
2389:     *redundantMesh = dm;
2390:     if (sf) *sf = NULL;
2391:     PetscFunctionReturn(PETSC_SUCCESS);
2392:   }
2393:   PetscCall(DMPlexGetGatherDM(dm, &gatherSF, &gatherDM));
2394:   if (!gatherDM) PetscFunctionReturn(PETSC_SUCCESS);
2395:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
2396:   PetscCall(DMPlexGetChart(gatherDM, &pStart, &pEnd));
2397:   numPoints = pEnd - pStart;
2398:   PetscCallMPI(MPI_Bcast(&numPoints, 1, MPIU_INT, 0, comm));
2399:   PetscCall(PetscMalloc1(numPoints, &points));
2400:   PetscCall(PetscSFCreate(comm, &migrationSF));
2401:   for (p = 0; p < numPoints; p++) {
2402:     points[p].index = p;
2403:     points[p].rank  = 0;
2404:   }
2405:   PetscCall(PetscSFSetGraph(migrationSF, pEnd - pStart, numPoints, NULL, PETSC_OWN_POINTER, points, PETSC_OWN_POINTER));
2406:   PetscCall(DMPlexCreate(comm, redundantMesh));
2407:   PetscCall(PetscObjectSetName((PetscObject)*redundantMesh, "Redundant Mesh"));
2408:   PetscCall(DMPlexMigrate(gatherDM, migrationSF, *redundantMesh));
2409:   /* This is to express that all point are in overlap */
2410:   PetscCall(DMPlexSetOverlap_Plex(*redundantMesh, NULL, PETSC_INT_MAX));
2411:   PetscCall(DMPlexCreatePointSF(*redundantMesh, migrationSF, PETSC_FALSE, &sfPoint));
2412:   PetscCall(DMSetPointSF(*redundantMesh, sfPoint));
2413:   PetscCall(DMGetCoordinateDM(*redundantMesh, &dmCoord));
2414:   if (dmCoord) PetscCall(DMSetPointSF(dmCoord, sfPoint));
2415:   PetscCall(PetscSFDestroy(&sfPoint));
2416:   if (sf) {
2417:     PetscSF tsf;

2419:     PetscCall(PetscSFCompose(gatherSF, migrationSF, &tsf));
2420:     PetscCall(DMPlexStratifyMigrationSF(dm, tsf, sf));
2421:     PetscCall(PetscSFDestroy(&tsf));
2422:   }
2423:   PetscCall(PetscSFDestroy(&migrationSF));
2424:   PetscCall(PetscSFDestroy(&gatherSF));
2425:   PetscCall(DMDestroy(&gatherDM));
2426:   PetscCall(DMCopyDisc(dm, *redundantMesh));
2427:   PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_FALSE, *redundantMesh));
2428:   PetscFunctionReturn(PETSC_SUCCESS);
2429: }

2431: /*@
2432:   DMPlexIsDistributed - Find out whether this `DM` is distributed, i.e. more than one rank owns some points.

2434:   Collective

2436:   Input Parameter:
2437: . dm - The `DM` object

2439:   Output Parameter:
2440: . distributed - Flag whether the `DM` is distributed

2442:   Level: intermediate

2444:   Notes:
2445:   This currently finds out whether at least two ranks have any DAG points.
2446:   This involves `MPI_Allreduce()` with one integer.
2447:   The result is currently not stashed so every call to this routine involves this global communication.

2449: .seealso: `DMPLEX`, `DMPlexDistribute()`, `DMPlexGetOverlap()`, `DMPlexIsInterpolated()`
2450: @*/
2451: PetscErrorCode DMPlexIsDistributed(DM dm, PetscBool *distributed)
2452: {
2453:   PetscInt    pStart, pEnd, count;
2454:   MPI_Comm    comm;
2455:   PetscMPIInt size;

2457:   PetscFunctionBegin;
2459:   PetscAssertPointer(distributed, 2);
2460:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2461:   PetscCallMPI(MPI_Comm_size(comm, &size));
2462:   if (size == 1) {
2463:     *distributed = PETSC_FALSE;
2464:     PetscFunctionReturn(PETSC_SUCCESS);
2465:   }
2466:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
2467:   count = (pEnd - pStart) > 0 ? 1 : 0;
2468:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &count, 1, MPIU_INT, MPI_SUM, comm));
2469:   *distributed = count > 1 ? PETSC_TRUE : PETSC_FALSE;
2470:   PetscFunctionReturn(PETSC_SUCCESS);
2471: }

2473: /*@
2474:   DMPlexDistributionSetName - Set the name of the specific parallel distribution

2476:   Input Parameters:
2477: + dm   - The `DM`
2478: - name - The name of the specific parallel distribution

2480:   Level: developer

2482:   Note:
2483:   If distribution name is set when saving, `DMPlexTopologyView()` saves the plex's
2484:   parallel distribution (i.e., partition, ownership, and local ordering of points) under
2485:   this name. Conversely, if distribution name is set when loading, `DMPlexTopologyLoad()`
2486:   loads the parallel distribution stored in file under this name.

2488: .seealso: `DMPLEX`, `DMPlexDistributionGetName()`, `DMPlexTopologyView()`, `DMPlexTopologyLoad()`
2489: @*/
2490: PetscErrorCode DMPlexDistributionSetName(DM dm, const char name[])
2491: {
2492:   DM_Plex *mesh = (DM_Plex *)dm->data;

2494:   PetscFunctionBegin;
2496:   if (name) PetscAssertPointer(name, 2);
2497:   PetscCall(PetscFree(mesh->distributionName));
2498:   PetscCall(PetscStrallocpy(name, &mesh->distributionName));
2499:   PetscFunctionReturn(PETSC_SUCCESS);
2500: }

2502: /*@
2503:   DMPlexDistributionGetName - Retrieve the name of the specific parallel distribution

2505:   Input Parameter:
2506: . dm - The `DM`

2508:   Output Parameter:
2509: . name - The name of the specific parallel distribution

2511:   Level: developer

2513:   Note:
2514:   If distribution name is set when saving, `DMPlexTopologyView()` saves the plex's
2515:   parallel distribution (i.e., partition, ownership, and local ordering of points) under
2516:   this name. Conversely, if distribution name is set when loading, `DMPlexTopologyLoad()`
2517:   loads the parallel distribution stored in file under this name.

2519: .seealso: `DMPLEX`, `DMPlexDistributionSetName()`, `DMPlexTopologyView()`, `DMPlexTopologyLoad()`
2520: @*/
2521: PetscErrorCode DMPlexDistributionGetName(DM dm, const char *name[])
2522: {
2523:   DM_Plex *mesh = (DM_Plex *)dm->data;

2525:   PetscFunctionBegin;
2527:   PetscAssertPointer(name, 2);
2528:   *name = mesh->distributionName;
2529:   PetscFunctionReturn(PETSC_SUCCESS);
2530: }