Actual source code: plexdistribute.c

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

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

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

     12:   Level: advanced

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

     17:      Any setting here overrides other configuration of DMPlex adjacency determination.

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

     27:   mesh->useradjacency = user;
     28:   mesh->useradjacencyctx = ctx;
     29:   return(0);
     30: }

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

     35:   Input Parameter:
     36: . dm      - The DM object

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

     42:   Level: advanced

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

     52:   if (user) *user = mesh->useradjacency;
     53:   if (ctx) *ctx = mesh->useradjacencyctx;
     54:   return(0);
     55: }

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

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

     64:   Level: intermediate

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

     74:   mesh->useAnchors = useAnchors;
     75:   return(0);
     76: }

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

     81:   Input Parameter:
     82: . dm      - The DM object

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

     87:   Level: intermediate

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

     98:   *useAnchors = mesh->useAnchors;
     99:   return(0);
    100: }

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

    109:   DMPlexGetConeSize(dm, p, &coneSize);
    110:   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:     DMPlexGetSupportSize(dm, point, &supportSize);
    117:     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:       if (numAdj > maxAdjSize) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize);
    123:     }
    124:   }
    125:   *adjSize = numAdj;
    126:   return(0);
    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;
    133:   PetscErrorCode  ierr;

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

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

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

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

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

    182: PetscErrorCode DMPlexGetAdjacency_Internal(DM dm, PetscInt p, PetscBool useCone, PetscBool useTransitiveClosure, PetscBool useAnchors, PetscInt *adjSize, PetscInt *adj[])
    183: {
    184:   static PetscInt asiz = 0;
    185:   PetscInt maxAnchors = 1;
    186:   PetscInt aStart = -1, aEnd = -1;
    187:   PetscInt maxAdjSize;
    188:   PetscSection aSec = NULL;
    189:   IS aIS = NULL;
    190:   const PetscInt *anchors;
    191:   DM_Plex *mesh = (DM_Plex *)dm->data;
    192:   PetscErrorCode  ierr;

    195:   if (useAnchors) {
    196:     DMPlexGetAnchors(dm,&aSec,&aIS);
    197:     if (aSec) {
    198:       PetscSectionGetMaxDof(aSec,&maxAnchors);
    199:       maxAnchors = PetscMax(1,maxAnchors);
    200:       PetscSectionGetChart(aSec,&aStart,&aEnd);
    201:       ISGetIndices(aIS,&anchors);
    202:     }
    203:   }
    204:   if (!*adj) {
    205:     PetscInt depth, coneSeries, supportSeries, maxC, maxS, pStart, pEnd;

    207:     DMPlexGetChart(dm, &pStart,&pEnd);
    208:     DMPlexGetDepth(dm, &depth);
    209:     depth = PetscMax(depth, -depth);
    210:     DMPlexGetMaxSizes(dm, &maxC, &maxS);
    211:     coneSeries    = (maxC > 1) ? ((PetscPowInt(maxC,depth+1)-1)/(maxC-1)) : depth+1;
    212:     supportSeries = (maxS > 1) ? ((PetscPowInt(maxS,depth+1)-1)/(maxS-1)) : depth+1;
    213:     asiz  = PetscMax(PetscPowInt(maxS,depth)*coneSeries,PetscPowInt(maxC,depth)*supportSeries);
    214:     asiz *= maxAnchors;
    215:     asiz  = PetscMin(asiz,pEnd-pStart);
    216:     PetscMalloc1(asiz,adj);
    217:   }
    218:   if (*adjSize < 0) *adjSize = asiz;
    219:   maxAdjSize = *adjSize;
    220:   if (mesh->useradjacency) {
    221:     (*mesh->useradjacency)(dm, p, adjSize, *adj, mesh->useradjacencyctx);
    222:   } else if (useTransitiveClosure) {
    223:     DMPlexGetAdjacency_Transitive_Internal(dm, p, useCone, adjSize, *adj);
    224:   } else if (useCone) {
    225:     DMPlexGetAdjacency_Cone_Internal(dm, p, adjSize, *adj);
    226:   } else {
    227:     DMPlexGetAdjacency_Support_Internal(dm, p, adjSize, *adj);
    228:   }
    229:   if (useAnchors && aSec) {
    230:     PetscInt origSize = *adjSize;
    231:     PetscInt numAdj = origSize;
    232:     PetscInt i = 0, j;
    233:     PetscInt *orig = *adj;

    235:     while (i < origSize) {
    236:       PetscInt p = orig[i];
    237:       PetscInt aDof = 0;

    239:       if (p >= aStart && p < aEnd) {
    240:         PetscSectionGetDof(aSec,p,&aDof);
    241:       }
    242:       if (aDof) {
    243:         PetscInt aOff;
    244:         PetscInt s, q;

    246:         for (j = i + 1; j < numAdj; j++) {
    247:           orig[j - 1] = orig[j];
    248:         }
    249:         origSize--;
    250:         numAdj--;
    251:         PetscSectionGetOffset(aSec,p,&aOff);
    252:         for (s = 0; s < aDof; ++s) {
    253:           for (q = 0; q < numAdj || ((void)(orig[numAdj++] = anchors[aOff+s]),0); ++q) {
    254:             if (anchors[aOff+s] == orig[q]) break;
    255:           }
    256:           if (numAdj > maxAdjSize) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize);
    257:         }
    258:       }
    259:       else {
    260:         i++;
    261:       }
    262:     }
    263:     *adjSize = numAdj;
    264:     ISRestoreIndices(aIS,&anchors);
    265:   }
    266:   return(0);
    267: }

    269: /*@
    270:   DMPlexGetAdjacency - Return all points adjacent to the given point

    272:   Input Parameters:
    273: + dm - The DM object
    274: . p  - The point
    275: . adjSize - The maximum size of adj if it is non-NULL, or PETSC_DETERMINE
    276: - adj - Either NULL so that the array is allocated, or an existing array with size adjSize

    278:   Output Parameters:
    279: + adjSize - The number of adjacent points
    280: - adj - The adjacent points

    282:   Level: advanced

    284:   Notes:
    285:     The user must PetscFree the adj array if it was not passed in.

    287: .seealso: DMSetAdjacency(), DMPlexDistribute(), DMCreateMatrix(), DMPlexPreallocateOperator()
    288: @*/
    289: PetscErrorCode DMPlexGetAdjacency(DM dm, PetscInt p, PetscInt *adjSize, PetscInt *adj[])
    290: {
    291:   PetscBool      useCone, useClosure, useAnchors;

    298:   DMGetBasicAdjacency(dm, &useCone, &useClosure);
    299:   DMPlexGetAdjacencyUseAnchors(dm, &useAnchors);
    300:   DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, adjSize, adj);
    301:   return(0);
    302: }

    304: /*@
    305:   DMPlexCreateTwoSidedProcessSF - Create an SF which just has process connectivity

    307:   Collective on dm

    309:   Input Parameters:
    310: + dm      - The DM
    311: - sfPoint - The PetscSF which encodes point connectivity

    313:   Output Parameters:
    314: + processRanks - A list of process neighbors, or NULL
    315: - sfProcess    - An SF encoding the two-sided process connectivity, or NULL

    317:   Level: developer

    319: .seealso: PetscSFCreate(), DMPlexCreateProcessSF()
    320: @*/
    321: PetscErrorCode DMPlexCreateTwoSidedProcessSF(DM dm, PetscSF sfPoint, PetscSection rootRankSection, IS rootRanks, PetscSection leafRankSection, IS leafRanks, IS *processRanks, PetscSF *sfProcess)
    322: {
    323:   const PetscSFNode *remotePoints;
    324:   PetscInt          *localPointsNew;
    325:   PetscSFNode       *remotePointsNew;
    326:   const PetscInt    *nranks;
    327:   PetscInt          *ranksNew;
    328:   PetscBT            neighbors;
    329:   PetscInt           pStart, pEnd, p, numLeaves, l, numNeighbors, n;
    330:   PetscMPIInt        size, proc, rank;
    331:   PetscErrorCode     ierr;

    338:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
    339:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
    340:   PetscSFGetGraph(sfPoint, NULL, &numLeaves, NULL, &remotePoints);
    341:   PetscBTCreate(size, &neighbors);
    342:   PetscBTMemzero(size, neighbors);
    343:   /* Compute root-to-leaf process connectivity */
    344:   PetscSectionGetChart(rootRankSection, &pStart, &pEnd);
    345:   ISGetIndices(rootRanks, &nranks);
    346:   for (p = pStart; p < pEnd; ++p) {
    347:     PetscInt ndof, noff, n;

    349:     PetscSectionGetDof(rootRankSection, p, &ndof);
    350:     PetscSectionGetOffset(rootRankSection, p, &noff);
    351:     for (n = 0; n < ndof; ++n) {PetscBTSet(neighbors, nranks[noff+n]);}
    352:   }
    353:   ISRestoreIndices(rootRanks, &nranks);
    354:   /* Compute leaf-to-neighbor process connectivity */
    355:   PetscSectionGetChart(leafRankSection, &pStart, &pEnd);
    356:   ISGetIndices(leafRanks, &nranks);
    357:   for (p = pStart; p < pEnd; ++p) {
    358:     PetscInt ndof, noff, n;

    360:     PetscSectionGetDof(leafRankSection, p, &ndof);
    361:     PetscSectionGetOffset(leafRankSection, p, &noff);
    362:     for (n = 0; n < ndof; ++n) {PetscBTSet(neighbors, nranks[noff+n]);}
    363:   }
    364:   ISRestoreIndices(leafRanks, &nranks);
    365:   /* Compute leaf-to-root process connectivity */
    366:   for (l = 0; l < numLeaves; ++l) {PetscBTSet(neighbors, remotePoints[l].rank);}
    367:   /* Calculate edges */
    368:   PetscBTClear(neighbors, rank);
    369:   for (proc = 0, numNeighbors = 0; proc < size; ++proc) {if (PetscBTLookup(neighbors, proc)) ++numNeighbors;}
    370:   PetscMalloc1(numNeighbors, &ranksNew);
    371:   PetscMalloc1(numNeighbors, &localPointsNew);
    372:   PetscMalloc1(numNeighbors, &remotePointsNew);
    373:   for (proc = 0, n = 0; proc < size; ++proc) {
    374:     if (PetscBTLookup(neighbors, proc)) {
    375:       ranksNew[n]              = proc;
    376:       localPointsNew[n]        = proc;
    377:       remotePointsNew[n].index = rank;
    378:       remotePointsNew[n].rank  = proc;
    379:       ++n;
    380:     }
    381:   }
    382:   PetscBTDestroy(&neighbors);
    383:   if (processRanks) {ISCreateGeneral(PetscObjectComm((PetscObject)dm), numNeighbors, ranksNew, PETSC_OWN_POINTER, processRanks);}
    384:   else              {PetscFree(ranksNew);}
    385:   if (sfProcess) {
    386:     PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);
    387:     PetscObjectSetName((PetscObject) *sfProcess, "Two-Sided Process SF");
    388:     PetscSFSetFromOptions(*sfProcess);
    389:     PetscSFSetGraph(*sfProcess, size, numNeighbors, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
    390:   }
    391:   return(0);
    392: }

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

    397:   Collective on dm

    399:   Input Parameter:
    400: . dm - The DM

    402:   Output Parameters:
    403: + rootSection - The number of leaves for a given root point
    404: . rootrank    - The rank of each edge into the root point
    405: . leafSection - The number of processes sharing a given leaf point
    406: - leafrank    - The rank of each process sharing a leaf point

    408:   Level: developer

    410: .seealso: DMPlexCreateOverlapLabel(), DMPlexDistribute(), DMPlexDistributeOverlap()
    411: @*/
    412: PetscErrorCode DMPlexDistributeOwnership(DM dm, PetscSection rootSection, IS *rootrank, PetscSection leafSection, IS *leafrank)
    413: {
    414:   MPI_Comm        comm;
    415:   PetscSF         sfPoint;
    416:   const PetscInt *rootdegree;
    417:   PetscInt       *myrank, *remoterank;
    418:   PetscInt        pStart, pEnd, p, nedges;
    419:   PetscMPIInt     rank;
    420:   PetscErrorCode  ierr;

    423:   PetscObjectGetComm((PetscObject) dm, &comm);
    424:   MPI_Comm_rank(comm, &rank);
    425:   DMPlexGetChart(dm, &pStart, &pEnd);
    426:   DMGetPointSF(dm, &sfPoint);
    427:   /* Compute number of leaves for each root */
    428:   PetscObjectSetName((PetscObject) rootSection, "Root Section");
    429:   PetscSectionSetChart(rootSection, pStart, pEnd);
    430:   PetscSFComputeDegreeBegin(sfPoint, &rootdegree);
    431:   PetscSFComputeDegreeEnd(sfPoint, &rootdegree);
    432:   for (p = pStart; p < pEnd; ++p) {PetscSectionSetDof(rootSection, p, rootdegree[p-pStart]);}
    433:   PetscSectionSetUp(rootSection);
    434:   /* Gather rank of each leaf to root */
    435:   PetscSectionGetStorageSize(rootSection, &nedges);
    436:   PetscMalloc1(pEnd-pStart, &myrank);
    437:   PetscMalloc1(nedges,  &remoterank);
    438:   for (p = 0; p < pEnd-pStart; ++p) myrank[p] = rank;
    439:   PetscSFGatherBegin(sfPoint, MPIU_INT, myrank, remoterank);
    440:   PetscSFGatherEnd(sfPoint, MPIU_INT, myrank, remoterank);
    441:   PetscFree(myrank);
    442:   ISCreateGeneral(comm, nedges, remoterank, PETSC_OWN_POINTER, rootrank);
    443:   /* Distribute remote ranks to leaves */
    444:   PetscObjectSetName((PetscObject) leafSection, "Leaf Section");
    445:   DMPlexDistributeFieldIS(dm, sfPoint, rootSection, *rootrank, leafSection, leafrank);
    446:   return(0);
    447: }

    449: /*@C
    450:   DMPlexCreateOverlapLabel - Compute owner information for shared points. This basically gets two-sided for an SF.

    452:   Collective on dm

    454:   Input Parameters:
    455: + dm          - The DM
    456: . levels      - Number of overlap levels
    457: . rootSection - The number of leaves for a given root point
    458: . rootrank    - The rank of each edge into the root point
    459: . leafSection - The number of processes sharing a given leaf point
    460: - leafrank    - The rank of each process sharing a leaf point

    462:   Output Parameter:
    463: . ovLabel     - DMLabel containing remote overlap contributions as point/rank pairings

    465:   Level: developer

    467: .seealso: DMPlexDistributeOwnership(), DMPlexDistribute()
    468: @*/
    469: PetscErrorCode DMPlexCreateOverlapLabel(DM dm, PetscInt levels, PetscSection rootSection, IS rootrank, PetscSection leafSection, IS leafrank, DMLabel *ovLabel)
    470: {
    471:   MPI_Comm           comm;
    472:   DMLabel            ovAdjByRank; /* A DMLabel containing all points adjacent to shared points, separated by rank (value in label) */
    473:   PetscSF            sfPoint;
    474:   const PetscSFNode *remote;
    475:   const PetscInt    *local;
    476:   const PetscInt    *nrank, *rrank;
    477:   PetscInt          *adj = NULL;
    478:   PetscInt           pStart, pEnd, p, sStart, sEnd, nleaves, l;
    479:   PetscMPIInt        rank, size;
    480:   PetscBool          flg;
    481:   PetscErrorCode     ierr;

    484:   *ovLabel = NULL;
    485:   PetscObjectGetComm((PetscObject) dm, &comm);
    486:   MPI_Comm_size(comm, &size);
    487:   MPI_Comm_rank(comm, &rank);
    488:   if (size == 1) return(0);
    489:   DMGetPointSF(dm, &sfPoint);
    490:   DMPlexGetChart(dm, &pStart, &pEnd);
    491:   if (!levels) {
    492:     /* Add owned points */
    493:     DMLabelCreate(PETSC_COMM_SELF, "Overlap label", ovLabel);
    494:     for (p = pStart; p < pEnd; ++p) {DMLabelSetValue(*ovLabel, p, rank);}
    495:     return(0);
    496:   }
    497:   PetscSectionGetChart(leafSection, &sStart, &sEnd);
    498:   PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote);
    499:   DMLabelCreate(PETSC_COMM_SELF, "Overlap adjacency", &ovAdjByRank);
    500:   /* Handle leaves: shared with the root point */
    501:   for (l = 0; l < nleaves; ++l) {
    502:     PetscInt adjSize = PETSC_DETERMINE, a;

    504:     DMPlexGetAdjacency(dm, local ? local[l] : l, &adjSize, &adj);
    505:     for (a = 0; a < adjSize; ++a) {DMLabelSetValue(ovAdjByRank, adj[a], remote[l].rank);}
    506:   }
    507:   ISGetIndices(rootrank, &rrank);
    508:   ISGetIndices(leafrank, &nrank);
    509:   /* Handle roots */
    510:   for (p = pStart; p < pEnd; ++p) {
    511:     PetscInt adjSize = PETSC_DETERMINE, neighbors = 0, noff, n, a;

    513:     if ((p >= sStart) && (p < sEnd)) {
    514:       /* Some leaves share a root with other leaves on different processes */
    515:       PetscSectionGetDof(leafSection, p, &neighbors);
    516:       if (neighbors) {
    517:         PetscSectionGetOffset(leafSection, p, &noff);
    518:         DMPlexGetAdjacency(dm, p, &adjSize, &adj);
    519:         for (n = 0; n < neighbors; ++n) {
    520:           const PetscInt remoteRank = nrank[noff+n];

    522:           if (remoteRank == rank) continue;
    523:           for (a = 0; a < adjSize; ++a) {DMLabelSetValue(ovAdjByRank, adj[a], remoteRank);}
    524:         }
    525:       }
    526:     }
    527:     /* Roots are shared with leaves */
    528:     PetscSectionGetDof(rootSection, p, &neighbors);
    529:     if (!neighbors) continue;
    530:     PetscSectionGetOffset(rootSection, p, &noff);
    531:     DMPlexGetAdjacency(dm, p, &adjSize, &adj);
    532:     for (n = 0; n < neighbors; ++n) {
    533:       const PetscInt remoteRank = rrank[noff+n];

    535:       if (remoteRank == rank) continue;
    536:       for (a = 0; a < adjSize; ++a) {DMLabelSetValue(ovAdjByRank, adj[a], remoteRank);}
    537:     }
    538:   }
    539:   PetscFree(adj);
    540:   ISRestoreIndices(rootrank, &rrank);
    541:   ISRestoreIndices(leafrank, &nrank);
    542:   /* Add additional overlap levels */
    543:   for (l = 1; l < levels; l++) {
    544:     /* Propagate point donations over SF to capture remote connections */
    545:     DMPlexPartitionLabelPropagate(dm, ovAdjByRank);
    546:     /* Add next level of point donations to the label */
    547:     DMPlexPartitionLabelAdjacency(dm, ovAdjByRank);
    548:   }
    549:   /* We require the closure in the overlap */
    550:   DMPlexPartitionLabelClosure(dm, ovAdjByRank);
    551:   PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-overlap_view", &flg);
    552:   if (flg) {
    553:     PetscViewer viewer;
    554:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm), &viewer);
    555:     DMLabelView(ovAdjByRank, viewer);
    556:   }
    557:   /* Invert sender to receiver label */
    558:   DMLabelCreate(PETSC_COMM_SELF, "Overlap label", ovLabel);
    559:   DMPlexPartitionLabelInvert(dm, ovAdjByRank, NULL, *ovLabel);
    560:   /* Add owned points, except for shared local points */
    561:   for (p = pStart; p < pEnd; ++p) {DMLabelSetValue(*ovLabel, p, rank);}
    562:   for (l = 0; l < nleaves; ++l) {
    563:     DMLabelClearValue(*ovLabel, local[l], rank);
    564:     DMLabelSetValue(*ovLabel, remote[l].index, remote[l].rank);
    565:   }
    566:   /* Clean up */
    567:   DMLabelDestroy(&ovAdjByRank);
    568:   return(0);
    569: }

    571: /*@C
    572:   DMPlexCreateOverlapMigrationSF - Create an SF describing the new mesh distribution to make the overlap described by the input SF

    574:   Collective on dm

    576:   Input Parameters:
    577: + dm          - The DM
    578: - overlapSF   - The SF mapping ghost points in overlap to owner points on other processes

    580:   Output Parameter:
    581: . migrationSF - An SF that maps original points in old locations to points in new locations

    583:   Level: developer

    585: .seealso: DMPlexCreateOverlapLabel(), DMPlexDistributeOverlap(), DMPlexDistribute()
    586: @*/
    587: PetscErrorCode DMPlexCreateOverlapMigrationSF(DM dm, PetscSF overlapSF, PetscSF *migrationSF)
    588: {
    589:   MPI_Comm           comm;
    590:   PetscMPIInt        rank, size;
    591:   PetscInt           d, dim, p, pStart, pEnd, nroots, nleaves, newLeaves, point, numSharedPoints;
    592:   PetscInt          *pointDepths, *remoteDepths, *ilocal;
    593:   PetscInt          *depthRecv, *depthShift, *depthIdx;
    594:   PetscSFNode       *iremote;
    595:   PetscSF            pointSF;
    596:   const PetscInt    *sharedLocal;
    597:   const PetscSFNode *overlapRemote, *sharedRemote;
    598:   PetscErrorCode     ierr;

    602:   PetscObjectGetComm((PetscObject)dm, &comm);
    603:   MPI_Comm_rank(comm, &rank);
    604:   MPI_Comm_size(comm, &size);
    605:   DMGetDimension(dm, &dim);

    607:   /* Before building the migration SF we need to know the new stratum offsets */
    608:   PetscSFGetGraph(overlapSF, &nroots, &nleaves, NULL, &overlapRemote);
    609:   PetscMalloc2(nroots, &pointDepths, nleaves, &remoteDepths);
    610:   for (d=0; d<dim+1; d++) {
    611:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
    612:     for (p=pStart; p<pEnd; p++) pointDepths[p] = d;
    613:   }
    614:   for (p=0; p<nleaves; p++) remoteDepths[p] = -1;
    615:   PetscSFBcastBegin(overlapSF, MPIU_INT, pointDepths, remoteDepths,MPI_REPLACE);
    616:   PetscSFBcastEnd(overlapSF, MPIU_INT, pointDepths, remoteDepths,MPI_REPLACE);

    618:   /* Count received points in each stratum and compute the internal strata shift */
    619:   PetscMalloc3(dim+1, &depthRecv, dim+1, &depthShift, dim+1, &depthIdx);
    620:   for (d=0; d<dim+1; d++) depthRecv[d]=0;
    621:   for (p=0; p<nleaves; p++) depthRecv[remoteDepths[p]]++;
    622:   depthShift[dim] = 0;
    623:   for (d=0; d<dim; d++) depthShift[d] = depthRecv[dim];
    624:   for (d=1; d<dim; d++) depthShift[d] += depthRecv[0];
    625:   for (d=dim-2; d>0; d--) depthShift[d] += depthRecv[d+1];
    626:   for (d=0; d<dim+1; d++) {
    627:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
    628:     depthIdx[d] = pStart + depthShift[d];
    629:   }

    631:   /* Form the overlap SF build an SF that describes the full overlap migration SF */
    632:   DMPlexGetChart(dm, &pStart, &pEnd);
    633:   newLeaves = pEnd - pStart + nleaves;
    634:   PetscMalloc1(newLeaves, &ilocal);
    635:   PetscMalloc1(newLeaves, &iremote);
    636:   /* First map local points to themselves */
    637:   for (d=0; d<dim+1; d++) {
    638:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
    639:     for (p=pStart; p<pEnd; p++) {
    640:       point = p + depthShift[d];
    641:       ilocal[point] = point;
    642:       iremote[point].index = p;
    643:       iremote[point].rank = rank;
    644:       depthIdx[d]++;
    645:     }
    646:   }

    648:   /* Add in the remote roots for currently shared points */
    649:   DMGetPointSF(dm, &pointSF);
    650:   PetscSFGetGraph(pointSF, NULL, &numSharedPoints, &sharedLocal, &sharedRemote);
    651:   for (d=0; d<dim+1; d++) {
    652:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
    653:     for (p=0; p<numSharedPoints; p++) {
    654:       if (pStart <= sharedLocal[p] && sharedLocal[p] < pEnd) {
    655:         point = sharedLocal[p] + depthShift[d];
    656:         iremote[point].index = sharedRemote[p].index;
    657:         iremote[point].rank = sharedRemote[p].rank;
    658:       }
    659:     }
    660:   }

    662:   /* Now add the incoming overlap points */
    663:   for (p=0; p<nleaves; p++) {
    664:     point = depthIdx[remoteDepths[p]];
    665:     ilocal[point] = point;
    666:     iremote[point].index = overlapRemote[p].index;
    667:     iremote[point].rank = overlapRemote[p].rank;
    668:     depthIdx[remoteDepths[p]]++;
    669:   }
    670:   PetscFree2(pointDepths,remoteDepths);

    672:   PetscSFCreate(comm, migrationSF);
    673:   PetscObjectSetName((PetscObject) *migrationSF, "Overlap Migration SF");
    674:   PetscSFSetFromOptions(*migrationSF);
    675:   DMPlexGetChart(dm, &pStart, &pEnd);
    676:   PetscSFSetGraph(*migrationSF, pEnd-pStart, newLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);

    678:   PetscFree3(depthRecv, depthShift, depthIdx);
    679:   return(0);
    680: }

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

    685:   Input Parameters:
    686: + dm          - The DM
    687: - sf          - A star forest with non-ordered leaves, usually defining a DM point migration

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

    692:   Note:
    693:   This lexicographically sorts by (depth, cellType)

    695:   Level: developer

    697: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexDistributeOverlap()
    698: @*/
    699: PetscErrorCode DMPlexStratifyMigrationSF(DM dm, PetscSF sf, PetscSF *migrationSF)
    700: {
    701:   MPI_Comm           comm;
    702:   PetscMPIInt        rank, size;
    703:   PetscInt           d, ldepth, depth, dim, p, pStart, pEnd, nroots, nleaves;
    704:   PetscSFNode       *pointDepths, *remoteDepths;
    705:   PetscInt          *ilocal;
    706:   PetscInt          *depthRecv, *depthShift, *depthIdx;
    707:   PetscInt          *ctRecv,    *ctShift,    *ctIdx;
    708:   const PetscSFNode *iremote;
    709:   PetscErrorCode     ierr;

    713:   PetscObjectGetComm((PetscObject) dm, &comm);
    714:   MPI_Comm_rank(comm, &rank);
    715:   MPI_Comm_size(comm, &size);
    716:   DMPlexGetDepth(dm, &ldepth);
    717:   DMGetDimension(dm, &dim);
    718:   MPIU_Allreduce(&ldepth, &depth, 1, MPIU_INT, MPI_MAX, comm);
    719:   if ((ldepth >= 0) && (depth != ldepth)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent Plex depth %d != %d", ldepth, depth);
    720:   PetscLogEventBegin(DMPLEX_PartStratSF,dm,0,0,0);

    722:   /* Before building the migration SF we need to know the new stratum offsets */
    723:   PetscSFGetGraph(sf, &nroots, &nleaves, NULL, &iremote);
    724:   PetscMalloc2(nroots, &pointDepths, nleaves, &remoteDepths);
    725:   for (d = 0; d < depth+1; ++d) {
    726:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
    727:     for (p = pStart; p < pEnd; ++p) {
    728:       DMPolytopeType ct;

    730:       DMPlexGetCellType(dm, p, &ct);
    731:       pointDepths[p].index = d;
    732:       pointDepths[p].rank  = ct;
    733:     }
    734:   }
    735:   for (p = 0; p < nleaves; ++p) {remoteDepths[p].index = -1; remoteDepths[p].rank = -1;}
    736:   PetscSFBcastBegin(sf, MPIU_2INT, pointDepths, remoteDepths,MPI_REPLACE);
    737:   PetscSFBcastEnd(sf, MPIU_2INT, pointDepths, remoteDepths,MPI_REPLACE);
    738:   /* Count received points in each stratum and compute the internal strata shift */
    739:   PetscCalloc6(depth+1, &depthRecv, depth+1, &depthShift, depth+1, &depthIdx, DM_NUM_POLYTOPES, &ctRecv, DM_NUM_POLYTOPES, &ctShift, DM_NUM_POLYTOPES, &ctIdx);
    740:   for (p = 0; p < nleaves; ++p) {
    741:     if (remoteDepths[p].rank < 0) {
    742:       ++depthRecv[remoteDepths[p].index];
    743:     } else {
    744:       ++ctRecv[remoteDepths[p].rank];
    745:     }
    746:   }
    747:   {
    748:     PetscInt depths[4], dims[4], shift = 0, i, c;

    750:     /* Cells (depth), Vertices (0), Faces (depth-1), Edges (1)
    751:          Consider DM_POLYTOPE_FV_GHOST and DM_POLYTOPE_INTERIOR_GHOST as cells
    752:      */
    753:     depths[0] = depth; depths[1] = 0; depths[2] = depth-1; depths[3] = 1;
    754:     dims[0]   = dim;   dims[1]   = 0; dims[2]   = dim-1;   dims[3]   = 1;
    755:     for (i = 0; i <= depth; ++i) {
    756:       const PetscInt dep = depths[i];
    757:       const PetscInt dim = dims[i];

    759:       for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
    760:         if (DMPolytopeTypeGetDim((DMPolytopeType) c) != dim && !(i == 0 && (c == DM_POLYTOPE_FV_GHOST || c == DM_POLYTOPE_INTERIOR_GHOST))) continue;
    761:         ctShift[c] = shift;
    762:         shift     += ctRecv[c];
    763:       }
    764:       depthShift[dep] = shift;
    765:       shift          += depthRecv[dep];
    766:     }
    767:     for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
    768:       const PetscInt ctDim = DMPolytopeTypeGetDim((DMPolytopeType) c);

    770:       if ((ctDim < 0 || ctDim > dim) && (c != DM_POLYTOPE_FV_GHOST && c != DM_POLYTOPE_INTERIOR_GHOST)) {
    771:         ctShift[c] = shift;
    772:         shift     += ctRecv[c];
    773:       }
    774:     }
    775:   }
    776:   /* Derive a new local permutation based on stratified indices */
    777:   PetscMalloc1(nleaves, &ilocal);
    778:   for (p = 0; p < nleaves; ++p) {
    779:     const PetscInt       dep = remoteDepths[p].index;
    780:     const DMPolytopeType ct  = (DMPolytopeType) remoteDepths[p].rank;

    782:     if ((PetscInt) ct < 0) {
    783:       ilocal[p] = depthShift[dep] + depthIdx[dep];
    784:       ++depthIdx[dep];
    785:     } else {
    786:       ilocal[p] = ctShift[ct] + ctIdx[ct];
    787:       ++ctIdx[ct];
    788:     }
    789:   }
    790:   PetscSFCreate(comm, migrationSF);
    791:   PetscObjectSetName((PetscObject) *migrationSF, "Migration SF");
    792:   PetscSFSetGraph(*migrationSF, nroots, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_COPY_VALUES);
    793:   PetscFree2(pointDepths,remoteDepths);
    794:   PetscFree6(depthRecv, depthShift, depthIdx, ctRecv, ctShift, ctIdx);
    795:   PetscLogEventEnd(DMPLEX_PartStratSF,dm,0,0,0);
    796:   return(0);
    797: }

    799: /*@
    800:   DMPlexDistributeField - Distribute field data to match a given PetscSF, usually the SF from mesh distribution

    802:   Collective on dm

    804:   Input Parameters:
    805: + dm - The DMPlex object
    806: . pointSF - The PetscSF describing the communication pattern
    807: . originalSection - The PetscSection for existing data layout
    808: - originalVec - The existing data in a local vector

    810:   Output Parameters:
    811: + newSection - The PetscSF describing the new data layout
    812: - newVec - The new data in a local vector

    814:   Level: developer

    816: .seealso: DMPlexDistribute(), DMPlexDistributeFieldIS(), DMPlexDistributeData()
    817: @*/
    818: PetscErrorCode DMPlexDistributeField(DM dm, PetscSF pointSF, PetscSection originalSection, Vec originalVec, PetscSection newSection, Vec newVec)
    819: {
    820:   PetscSF        fieldSF;
    821:   PetscInt      *remoteOffsets, fieldSize;
    822:   PetscScalar   *originalValues, *newValues;

    826:   PetscLogEventBegin(DMPLEX_DistributeField,dm,0,0,0);
    827:   PetscSFDistributeSection(pointSF, originalSection, &remoteOffsets, newSection);

    829:   PetscSectionGetStorageSize(newSection, &fieldSize);
    830:   VecSetSizes(newVec, fieldSize, PETSC_DETERMINE);
    831:   VecSetType(newVec,dm->vectype);

    833:   VecGetArray(originalVec, &originalValues);
    834:   VecGetArray(newVec, &newValues);
    835:   PetscSFCreateSectionSF(pointSF, originalSection, remoteOffsets, newSection, &fieldSF);
    836:   PetscFree(remoteOffsets);
    837:   PetscSFBcastBegin(fieldSF, MPIU_SCALAR, originalValues, newValues,MPI_REPLACE);
    838:   PetscSFBcastEnd(fieldSF, MPIU_SCALAR, originalValues, newValues,MPI_REPLACE);
    839:   PetscSFDestroy(&fieldSF);
    840:   VecRestoreArray(newVec, &newValues);
    841:   VecRestoreArray(originalVec, &originalValues);
    842:   PetscLogEventEnd(DMPLEX_DistributeField,dm,0,0,0);
    843:   return(0);
    844: }

    846: /*@
    847:   DMPlexDistributeFieldIS - Distribute field data to match a given PetscSF, usually the SF from mesh distribution

    849:   Collective on dm

    851:   Input Parameters:
    852: + dm - The DMPlex object
    853: . pointSF - The PetscSF describing the communication pattern
    854: . originalSection - The PetscSection for existing data layout
    855: - originalIS - The existing data

    857:   Output Parameters:
    858: + newSection - The PetscSF describing the new data layout
    859: - newIS - The new data

    861:   Level: developer

    863: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexDistributeData()
    864: @*/
    865: PetscErrorCode DMPlexDistributeFieldIS(DM dm, PetscSF pointSF, PetscSection originalSection, IS originalIS, PetscSection newSection, IS *newIS)
    866: {
    867:   PetscSF         fieldSF;
    868:   PetscInt       *newValues, *remoteOffsets, fieldSize;
    869:   const PetscInt *originalValues;
    870:   PetscErrorCode  ierr;

    873:   PetscLogEventBegin(DMPLEX_DistributeField,dm,0,0,0);
    874:   PetscSFDistributeSection(pointSF, originalSection, &remoteOffsets, newSection);

    876:   PetscSectionGetStorageSize(newSection, &fieldSize);
    877:   PetscMalloc1(fieldSize, &newValues);

    879:   ISGetIndices(originalIS, &originalValues);
    880:   PetscSFCreateSectionSF(pointSF, originalSection, remoteOffsets, newSection, &fieldSF);
    881:   PetscFree(remoteOffsets);
    882:   PetscSFBcastBegin(fieldSF, MPIU_INT, (PetscInt *) originalValues, newValues,MPI_REPLACE);
    883:   PetscSFBcastEnd(fieldSF, MPIU_INT, (PetscInt *) originalValues, newValues,MPI_REPLACE);
    884:   PetscSFDestroy(&fieldSF);
    885:   ISRestoreIndices(originalIS, &originalValues);
    886:   ISCreateGeneral(PetscObjectComm((PetscObject) pointSF), fieldSize, newValues, PETSC_OWN_POINTER, newIS);
    887:   PetscLogEventEnd(DMPLEX_DistributeField,dm,0,0,0);
    888:   return(0);
    889: }

    891: /*@
    892:   DMPlexDistributeData - Distribute field data to match a given PetscSF, usually the SF from mesh distribution

    894:   Collective on dm

    896:   Input Parameters:
    897: + dm - The DMPlex object
    898: . pointSF - The PetscSF describing the communication pattern
    899: . originalSection - The PetscSection for existing data layout
    900: . datatype - The type of data
    901: - originalData - The existing data

    903:   Output Parameters:
    904: + newSection - The PetscSection describing the new data layout
    905: - newData - The new data

    907:   Level: developer

    909: .seealso: DMPlexDistribute(), DMPlexDistributeField()
    910: @*/
    911: PetscErrorCode DMPlexDistributeData(DM dm, PetscSF pointSF, PetscSection originalSection, MPI_Datatype datatype, void *originalData, PetscSection newSection, void **newData)
    912: {
    913:   PetscSF        fieldSF;
    914:   PetscInt      *remoteOffsets, fieldSize;
    915:   PetscMPIInt    dataSize;

    919:   PetscLogEventBegin(DMPLEX_DistributeData,dm,0,0,0);
    920:   PetscSFDistributeSection(pointSF, originalSection, &remoteOffsets, newSection);

    922:   PetscSectionGetStorageSize(newSection, &fieldSize);
    923:   MPI_Type_size(datatype, &dataSize);
    924:   PetscMalloc(fieldSize * dataSize, newData);

    926:   PetscSFCreateSectionSF(pointSF, originalSection, remoteOffsets, newSection, &fieldSF);
    927:   PetscFree(remoteOffsets);
    928:   PetscSFBcastBegin(fieldSF, datatype, originalData, *newData,MPI_REPLACE);
    929:   PetscSFBcastEnd(fieldSF, datatype, originalData, *newData,MPI_REPLACE);
    930:   PetscSFDestroy(&fieldSF);
    931:   PetscLogEventEnd(DMPLEX_DistributeData,dm,0,0,0);
    932:   return(0);
    933: }

    935: static PetscErrorCode DMPlexDistributeCones(DM dm, PetscSF migrationSF, ISLocalToGlobalMapping original, ISLocalToGlobalMapping renumbering, DM dmParallel)
    936: {
    937:   DM_Plex               *pmesh = (DM_Plex*) (dmParallel)->data;
    938:   MPI_Comm               comm;
    939:   PetscSF                coneSF;
    940:   PetscSection           originalConeSection, newConeSection;
    941:   PetscInt              *remoteOffsets, *cones, *globCones, *newCones, newConesSize;
    942:   PetscBool              flg;
    943:   PetscErrorCode         ierr;

    948:   PetscLogEventBegin(DMPLEX_DistributeCones,dm,0,0,0);
    949:   /* Distribute cone section */
    950:   PetscObjectGetComm((PetscObject)dm, &comm);
    951:   DMPlexGetConeSection(dm, &originalConeSection);
    952:   DMPlexGetConeSection(dmParallel, &newConeSection);
    953:   PetscSFDistributeSection(migrationSF, originalConeSection, &remoteOffsets, newConeSection);
    954:   DMSetUp(dmParallel);
    955:   {
    956:     PetscInt pStart, pEnd, p;

    958:     PetscSectionGetChart(newConeSection, &pStart, &pEnd);
    959:     for (p = pStart; p < pEnd; ++p) {
    960:       PetscInt coneSize;
    961:       PetscSectionGetDof(newConeSection, p, &coneSize);
    962:       pmesh->maxConeSize = PetscMax(pmesh->maxConeSize, coneSize);
    963:     }
    964:   }
    965:   /* Communicate and renumber cones */
    966:   PetscSFCreateSectionSF(migrationSF, originalConeSection, remoteOffsets, newConeSection, &coneSF);
    967:   PetscFree(remoteOffsets);
    968:   DMPlexGetCones(dm, &cones);
    969:   if (original) {
    970:     PetscInt numCones;

    972:     PetscSectionGetStorageSize(originalConeSection,&numCones);
    973:     PetscMalloc1(numCones,&globCones);
    974:     ISLocalToGlobalMappingApplyBlock(original, numCones, cones, globCones);
    975:   } else {
    976:     globCones = cones;
    977:   }
    978:   DMPlexGetCones(dmParallel, &newCones);
    979:   PetscSFBcastBegin(coneSF, MPIU_INT, globCones, newCones,MPI_REPLACE);
    980:   PetscSFBcastEnd(coneSF, MPIU_INT, globCones, newCones,MPI_REPLACE);
    981:   if (original) {
    982:     PetscFree(globCones);
    983:   }
    984:   PetscSectionGetStorageSize(newConeSection, &newConesSize);
    985:   ISGlobalToLocalMappingApplyBlock(renumbering, IS_GTOLM_MASK, newConesSize, newCones, NULL, newCones);
    986:   if (PetscDefined(USE_DEBUG)) {
    987:     PetscInt  p;
    988:     PetscBool valid = PETSC_TRUE;
    989:     for (p = 0; p < newConesSize; ++p) {
    990:       if (newCones[p] < 0) {valid = PETSC_FALSE; PetscPrintf(PETSC_COMM_SELF, "[%d] Point %D not in overlap SF\n", PetscGlobalRank,p);}
    991:     }
    992:     if (!valid) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid global to local map");
    993:   }
    994:   PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-cones_view", &flg);
    995:   if (flg) {
    996:     PetscPrintf(comm, "Serial Cone Section:\n");
    997:     PetscSectionView(originalConeSection, PETSC_VIEWER_STDOUT_(comm));
    998:     PetscPrintf(comm, "Parallel Cone Section:\n");
    999:     PetscSectionView(newConeSection, PETSC_VIEWER_STDOUT_(comm));
    1000:     PetscSFView(coneSF, NULL);
    1001:   }
    1002:   DMPlexGetConeOrientations(dm, &cones);
    1003:   DMPlexGetConeOrientations(dmParallel, &newCones);
    1004:   PetscSFBcastBegin(coneSF, MPIU_INT, cones, newCones,MPI_REPLACE);
    1005:   PetscSFBcastEnd(coneSF, MPIU_INT, cones, newCones,MPI_REPLACE);
    1006:   PetscSFDestroy(&coneSF);
    1007:   PetscLogEventEnd(DMPLEX_DistributeCones,dm,0,0,0);
    1008:   /* Create supports and stratify DMPlex */
    1009:   {
    1010:     PetscInt pStart, pEnd;

    1012:     PetscSectionGetChart(pmesh->coneSection, &pStart, &pEnd);
    1013:     PetscSectionSetChart(pmesh->supportSection, pStart, pEnd);
    1014:   }
    1015:   DMPlexSymmetrize(dmParallel);
    1016:   DMPlexStratify(dmParallel);
    1017:   {
    1018:     PetscBool useCone, useClosure, useAnchors;

    1020:     DMGetBasicAdjacency(dm, &useCone, &useClosure);
    1021:     DMSetBasicAdjacency(dmParallel, useCone, useClosure);
    1022:     DMPlexGetAdjacencyUseAnchors(dm, &useAnchors);
    1023:     DMPlexSetAdjacencyUseAnchors(dmParallel, useAnchors);
    1024:   }
    1025:   return(0);
    1026: }

    1028: static PetscErrorCode DMPlexDistributeCoordinates(DM dm, PetscSF migrationSF, DM dmParallel)
    1029: {
    1030:   MPI_Comm         comm;
    1031:   PetscSection     originalCoordSection, newCoordSection;
    1032:   Vec              originalCoordinates, newCoordinates;
    1033:   PetscInt         bs;
    1034:   PetscBool        isper;
    1035:   const char      *name;
    1036:   const PetscReal *maxCell, *L;
    1037:   const DMBoundaryType *bd;
    1038:   PetscErrorCode   ierr;

    1044:   PetscObjectGetComm((PetscObject)dm, &comm);
    1045:   DMGetCoordinateSection(dm, &originalCoordSection);
    1046:   DMGetCoordinateSection(dmParallel, &newCoordSection);
    1047:   DMGetCoordinatesLocal(dm, &originalCoordinates);
    1048:   if (originalCoordinates) {
    1049:     VecCreate(PETSC_COMM_SELF, &newCoordinates);
    1050:     PetscObjectGetName((PetscObject) originalCoordinates, &name);
    1051:     PetscObjectSetName((PetscObject) newCoordinates, name);

    1053:     DMPlexDistributeField(dm, migrationSF, originalCoordSection, originalCoordinates, newCoordSection, newCoordinates);
    1054:     DMSetCoordinatesLocal(dmParallel, newCoordinates);
    1055:     VecGetBlockSize(originalCoordinates, &bs);
    1056:     VecSetBlockSize(newCoordinates, bs);
    1057:     VecDestroy(&newCoordinates);
    1058:   }
    1059:   DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
    1060:   DMSetPeriodicity(dmParallel, isper, maxCell, L, bd);
    1061:   return(0);
    1062: }

    1064: static PetscErrorCode DMPlexDistributeLabels(DM dm, PetscSF migrationSF, DM dmParallel)
    1065: {
    1066:   DM_Plex         *mesh = (DM_Plex*) dm->data;
    1067:   MPI_Comm         comm;
    1068:   DMLabel          depthLabel;
    1069:   PetscMPIInt      rank;
    1070:   PetscInt         depth, d, numLabels, numLocalLabels, l;
    1071:   PetscBool        hasLabels = PETSC_FALSE, lsendDepth, sendDepth;
    1072:   PetscObjectState depthState = -1;
    1073:   PetscErrorCode   ierr;

    1079:   PetscLogEventBegin(DMPLEX_DistributeLabels,dm,0,0,0);
    1080:   PetscObjectGetComm((PetscObject)dm, &comm);
    1081:   MPI_Comm_rank(comm, &rank);

    1083:   /* If the user has changed the depth label, communicate it instead */
    1084:   DMPlexGetDepth(dm, &depth);
    1085:   DMPlexGetDepthLabel(dm, &depthLabel);
    1086:   if (depthLabel) {PetscObjectStateGet((PetscObject) depthLabel, &depthState);}
    1087:   lsendDepth = mesh->depthState != depthState ? PETSC_TRUE : PETSC_FALSE;
    1088:   MPIU_Allreduce(&lsendDepth, &sendDepth, 1, MPIU_BOOL, MPI_LOR, comm);
    1089:   if (sendDepth) {
    1090:     DMPlexGetDepthLabel(dmParallel, &dmParallel->depthLabel);
    1091:     DMRemoveLabelBySelf(dmParallel, &dmParallel->depthLabel, PETSC_FALSE);
    1092:   }
    1093:   /* Everyone must have either the same number of labels, or none */
    1094:   DMGetNumLabels(dm, &numLocalLabels);
    1095:   numLabels = numLocalLabels;
    1096:   MPI_Bcast(&numLabels, 1, MPIU_INT, 0, comm);
    1097:   if (numLabels == numLocalLabels) hasLabels = PETSC_TRUE;
    1098:   for (l = 0; l < numLabels; ++l) {
    1099:     DMLabel     label = NULL, labelNew = NULL;
    1100:     PetscBool   isDepth, lisOutput = PETSC_TRUE, isOutput;
    1101:     const char *name = NULL;

    1103:     if (hasLabels) {
    1104:       DMGetLabelByNum(dm, l, &label);
    1105:       /* Skip "depth" because it is recreated */
    1106:       PetscObjectGetName((PetscObject) label, &name);
    1107:       PetscStrcmp(name, "depth", &isDepth);
    1108:     }
    1109:     MPI_Bcast(&isDepth, 1, MPIU_BOOL, 0, comm);
    1110:     if (isDepth && !sendDepth) continue;
    1111:     DMLabelDistribute(label, migrationSF, &labelNew);
    1112:     if (isDepth) {
    1113:       /* Put in any missing strata which can occur if users are managing the depth label themselves */
    1114:       PetscInt gdepth;

    1116:       MPIU_Allreduce(&depth, &gdepth, 1, MPIU_INT, MPI_MAX, comm);
    1117:       if ((depth >= 0) && (gdepth != depth)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent Plex depth %d != %d", depth, gdepth);
    1118:       for (d = 0; d <= gdepth; ++d) {
    1119:         PetscBool has;

    1121:         DMLabelHasStratum(labelNew, d, &has);
    1122:         if (!has) {DMLabelAddStratum(labelNew, d);}
    1123:       }
    1124:     }
    1125:     DMAddLabel(dmParallel, labelNew);
    1126:     /* Put the output flag in the new label */
    1127:     if (hasLabels) {DMGetLabelOutput(dm, name, &lisOutput);}
    1128:     MPIU_Allreduce(&lisOutput, &isOutput, 1, MPIU_BOOL, MPI_LAND, comm);
    1129:     PetscObjectGetName((PetscObject) labelNew, &name);
    1130:     DMSetLabelOutput(dmParallel, name, isOutput);
    1131:     DMLabelDestroy(&labelNew);
    1132:   }
    1133:   PetscLogEventEnd(DMPLEX_DistributeLabels,dm,0,0,0);
    1134:   return(0);
    1135: }

    1137: static PetscErrorCode DMPlexDistributeSetupTree(DM dm, PetscSF migrationSF, ISLocalToGlobalMapping original, ISLocalToGlobalMapping renumbering, DM dmParallel)
    1138: {
    1139:   DM_Plex        *mesh  = (DM_Plex*) dm->data;
    1140:   DM_Plex        *pmesh = (DM_Plex*) (dmParallel)->data;
    1141:   MPI_Comm        comm;
    1142:   DM              refTree;
    1143:   PetscSection    origParentSection, newParentSection;
    1144:   PetscInt        *origParents, *origChildIDs;
    1145:   PetscBool       flg;
    1146:   PetscErrorCode  ierr;

    1151:   PetscObjectGetComm((PetscObject)dm, &comm);

    1153:   /* Set up tree */
    1154:   DMPlexGetReferenceTree(dm,&refTree);
    1155:   DMPlexSetReferenceTree(dmParallel,refTree);
    1156:   DMPlexGetTree(dm,&origParentSection,&origParents,&origChildIDs,NULL,NULL);
    1157:   if (origParentSection) {
    1158:     PetscInt        pStart, pEnd;
    1159:     PetscInt        *newParents, *newChildIDs, *globParents;
    1160:     PetscInt        *remoteOffsetsParents, newParentSize;
    1161:     PetscSF         parentSF;

    1163:     DMPlexGetChart(dmParallel, &pStart, &pEnd);
    1164:     PetscSectionCreate(PetscObjectComm((PetscObject)dmParallel),&newParentSection);
    1165:     PetscSectionSetChart(newParentSection,pStart,pEnd);
    1166:     PetscSFDistributeSection(migrationSF, origParentSection, &remoteOffsetsParents, newParentSection);
    1167:     PetscSFCreateSectionSF(migrationSF, origParentSection, remoteOffsetsParents, newParentSection, &parentSF);
    1168:     PetscFree(remoteOffsetsParents);
    1169:     PetscSectionGetStorageSize(newParentSection,&newParentSize);
    1170:     PetscMalloc2(newParentSize,&newParents,newParentSize,&newChildIDs);
    1171:     if (original) {
    1172:       PetscInt numParents;

    1174:       PetscSectionGetStorageSize(origParentSection,&numParents);
    1175:       PetscMalloc1(numParents,&globParents);
    1176:       ISLocalToGlobalMappingApplyBlock(original, numParents, origParents, globParents);
    1177:     }
    1178:     else {
    1179:       globParents = origParents;
    1180:     }
    1181:     PetscSFBcastBegin(parentSF, MPIU_INT, globParents, newParents,MPI_REPLACE);
    1182:     PetscSFBcastEnd(parentSF, MPIU_INT, globParents, newParents,MPI_REPLACE);
    1183:     if (original) {
    1184:       PetscFree(globParents);
    1185:     }
    1186:     PetscSFBcastBegin(parentSF, MPIU_INT, origChildIDs, newChildIDs,MPI_REPLACE);
    1187:     PetscSFBcastEnd(parentSF, MPIU_INT, origChildIDs, newChildIDs,MPI_REPLACE);
    1188:     ISGlobalToLocalMappingApplyBlock(renumbering,IS_GTOLM_MASK, newParentSize, newParents, NULL, newParents);
    1189:     if (PetscDefined(USE_DEBUG)) {
    1190:       PetscInt  p;
    1191:       PetscBool valid = PETSC_TRUE;
    1192:       for (p = 0; p < newParentSize; ++p) {
    1193:         if (newParents[p] < 0) {valid = PETSC_FALSE; PetscPrintf(PETSC_COMM_SELF, "Point %d not in overlap SF\n", p);}
    1194:       }
    1195:       if (!valid) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid global to local map");
    1196:     }
    1197:     PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-parents_view", &flg);
    1198:     if (flg) {
    1199:       PetscPrintf(comm, "Serial Parent Section: \n");
    1200:       PetscSectionView(origParentSection, PETSC_VIEWER_STDOUT_(comm));
    1201:       PetscPrintf(comm, "Parallel Parent Section: \n");
    1202:       PetscSectionView(newParentSection, PETSC_VIEWER_STDOUT_(comm));
    1203:       PetscSFView(parentSF, NULL);
    1204:     }
    1205:     DMPlexSetTree(dmParallel,newParentSection,newParents,newChildIDs);
    1206:     PetscSectionDestroy(&newParentSection);
    1207:     PetscFree2(newParents,newChildIDs);
    1208:     PetscSFDestroy(&parentSF);
    1209:   }
    1210:   pmesh->useAnchors = mesh->useAnchors;
    1211:   return(0);
    1212: }

    1214: PETSC_UNUSED static PetscErrorCode DMPlexDistributeSF(DM dm, PetscSF migrationSF, DM dmParallel)
    1215: {
    1216:   PetscMPIInt            rank, size;
    1217:   MPI_Comm               comm;
    1218:   PetscErrorCode         ierr;

    1224:   /* Create point SF for parallel mesh */
    1225:   PetscLogEventBegin(DMPLEX_DistributeSF,dm,0,0,0);
    1226:   PetscObjectGetComm((PetscObject)dm, &comm);
    1227:   MPI_Comm_rank(comm, &rank);
    1228:   MPI_Comm_size(comm, &size);
    1229:   {
    1230:     const PetscInt *leaves;
    1231:     PetscSFNode    *remotePoints, *rowners, *lowners;
    1232:     PetscInt        numRoots, numLeaves, numGhostPoints = 0, p, gp, *ghostPoints;
    1233:     PetscInt        pStart, pEnd;

    1235:     DMPlexGetChart(dmParallel, &pStart, &pEnd);
    1236:     PetscSFGetGraph(migrationSF, &numRoots, &numLeaves, &leaves, NULL);
    1237:     PetscMalloc2(numRoots,&rowners,numLeaves,&lowners);
    1238:     for (p=0; p<numRoots; p++) {
    1239:       rowners[p].rank  = -1;
    1240:       rowners[p].index = -1;
    1241:     }
    1242:     PetscSFBcastBegin(migrationSF, MPIU_2INT, rowners, lowners,MPI_REPLACE);
    1243:     PetscSFBcastEnd(migrationSF, MPIU_2INT, rowners, lowners,MPI_REPLACE);
    1244:     for (p = 0; p < numLeaves; ++p) {
    1245:       if (lowners[p].rank < 0 || lowners[p].rank == rank) { /* Either put in a bid or we know we own it */
    1246:         lowners[p].rank  = rank;
    1247:         lowners[p].index = leaves ? leaves[p] : p;
    1248:       } else if (lowners[p].rank >= 0) { /* Point already claimed so flag so that MAXLOC does not listen to us */
    1249:         lowners[p].rank  = -2;
    1250:         lowners[p].index = -2;
    1251:       }
    1252:     }
    1253:     for (p=0; p<numRoots; p++) { /* Root must not participate in the rediction, flag so that MAXLOC does not use */
    1254:       rowners[p].rank  = -3;
    1255:       rowners[p].index = -3;
    1256:     }
    1257:     PetscSFReduceBegin(migrationSF, MPIU_2INT, lowners, rowners, MPI_MAXLOC);
    1258:     PetscSFReduceEnd(migrationSF, MPIU_2INT, lowners, rowners, MPI_MAXLOC);
    1259:     PetscSFBcastBegin(migrationSF, MPIU_2INT, rowners, lowners,MPI_REPLACE);
    1260:     PetscSFBcastEnd(migrationSF, MPIU_2INT, rowners, lowners,MPI_REPLACE);
    1261:     for (p = 0; p < numLeaves; ++p) {
    1262:       if (lowners[p].rank < 0 || lowners[p].index < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cell partition corrupt: point not claimed");
    1263:       if (lowners[p].rank != rank) ++numGhostPoints;
    1264:     }
    1265:     PetscMalloc1(numGhostPoints, &ghostPoints);
    1266:     PetscMalloc1(numGhostPoints, &remotePoints);
    1267:     for (p = 0, gp = 0; p < numLeaves; ++p) {
    1268:       if (lowners[p].rank != rank) {
    1269:         ghostPoints[gp]        = leaves ? leaves[p] : p;
    1270:         remotePoints[gp].rank  = lowners[p].rank;
    1271:         remotePoints[gp].index = lowners[p].index;
    1272:         ++gp;
    1273:       }
    1274:     }
    1275:     PetscFree2(rowners,lowners);
    1276:     PetscSFSetGraph((dmParallel)->sf, pEnd - pStart, numGhostPoints, ghostPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
    1277:     PetscSFSetFromOptions((dmParallel)->sf);
    1278:   }
    1279:   {
    1280:     PetscBool useCone, useClosure, useAnchors;

    1282:     DMGetBasicAdjacency(dm, &useCone, &useClosure);
    1283:     DMSetBasicAdjacency(dmParallel, useCone, useClosure);
    1284:     DMPlexGetAdjacencyUseAnchors(dm, &useAnchors);
    1285:     DMPlexSetAdjacencyUseAnchors(dmParallel, useAnchors);
    1286:   }
    1287:   PetscLogEventEnd(DMPLEX_DistributeSF,dm,0,0,0);
    1288:   return(0);
    1289: }

    1291: /*@
    1292:   DMPlexSetPartitionBalance - Should distribution of the DM attempt to balance the shared point partition?

    1294:   Input Parameters:
    1295: + dm - The DMPlex object
    1296: - flg - Balance the partition?

    1298:   Level: intermediate

    1300: .seealso: DMPlexDistribute(), DMPlexGetPartitionBalance()
    1301: @*/
    1302: PetscErrorCode DMPlexSetPartitionBalance(DM dm, PetscBool flg)
    1303: {
    1304:   DM_Plex *mesh = (DM_Plex *)dm->data;

    1307:   mesh->partitionBalance = flg;
    1308:   return(0);
    1309: }

    1311: /*@
    1312:   DMPlexGetPartitionBalance - Does distribution of the DM attempt to balance the shared point partition?

    1314:   Input Parameter:
    1315: . dm - The DMPlex object

    1317:   Output Parameter:
    1318: . flg - Balance the partition?

    1320:   Level: intermediate

    1322: .seealso: DMPlexDistribute(), DMPlexSetPartitionBalance()
    1323: @*/
    1324: PetscErrorCode DMPlexGetPartitionBalance(DM dm, PetscBool *flg)
    1325: {
    1326:   DM_Plex *mesh = (DM_Plex *)dm->data;

    1331:   *flg = mesh->partitionBalance;
    1332:   return(0);
    1333: }

    1335: typedef struct {
    1336:   PetscInt vote, rank, index;
    1337: } Petsc3Int;

    1339: /* MaxLoc, but carry a third piece of information around */
    1340: static void MPIAPI MaxLocCarry(void *in_, void *inout_, PetscMPIInt *len_, MPI_Datatype *dtype)
    1341: {
    1342:   Petsc3Int *a = (Petsc3Int *)inout_;
    1343:   Petsc3Int *b = (Petsc3Int *)in_;
    1344:   PetscInt i, len = *len_;
    1345:   for (i = 0; i < len; i++) {
    1346:     if (a[i].vote < b[i].vote) {
    1347:       a[i].vote = b[i].vote;
    1348:       a[i].rank = b[i].rank;
    1349:       a[i].index = b[i].index;
    1350:     } else if (a[i].vote <= b[i].vote) {
    1351:       if (a[i].rank >= b[i].rank) {
    1352:         a[i].rank = b[i].rank;
    1353:         a[i].index = b[i].index;
    1354:       }
    1355:     }
    1356:   }
    1357: }

    1359: /*@C
    1360:   DMPlexCreatePointSF - Build a point SF from an SF describing a point migration

    1362:   Input Parameters:
    1363: + dm          - The source DMPlex object
    1364: . migrationSF - The star forest that describes the parallel point remapping
    1365: . ownership   - Flag causing a vote to determine point ownership

    1367:   Output Parameter:
    1368: - pointSF     - The star forest describing the point overlap in the remapped DM

    1370:   Notes:
    1371:   Output pointSF is guaranteed to have the array of local indices (leaves) sorted.

    1373:   Level: developer

    1375: .seealso: DMPlexDistribute(), DMPlexDistributeOverlap()
    1376: @*/
    1377: PetscErrorCode DMPlexCreatePointSF(DM dm, PetscSF migrationSF, PetscBool ownership, PetscSF *pointSF)
    1378: {
    1379:   PetscMPIInt        rank, size;
    1380:   PetscInt           p, nroots, nleaves, idx, npointLeaves;
    1381:   PetscInt          *pointLocal;
    1382:   const PetscInt    *leaves;
    1383:   const PetscSFNode *roots;
    1384:   PetscSFNode       *rootNodes, *leafNodes, *pointRemote;
    1385:   Vec                shifts;
    1386:   const PetscInt     numShifts = 13759;
    1387:   const PetscScalar *shift = NULL;
    1388:   const PetscBool    shiftDebug = PETSC_FALSE;
    1389:   PetscBool          balance;
    1390:   PetscErrorCode     ierr;

    1394:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
    1395:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
    1396:   PetscLogEventBegin(DMPLEX_CreatePointSF,dm,0,0,0);

    1398:   DMPlexGetPartitionBalance(dm, &balance);
    1399:   PetscSFGetGraph(migrationSF, &nroots, &nleaves, &leaves, &roots);
    1400:   PetscMalloc2(nroots, &rootNodes, nleaves, &leafNodes);
    1401:   if (ownership) {
    1402:     MPI_Op       op;
    1403:     MPI_Datatype datatype;
    1404:     Petsc3Int   *rootVote = NULL, *leafVote = NULL;
    1405:     /* If balancing, we compute a random cyclic shift of the rank for each remote point. That way, the max will evenly distribute among ranks. */
    1406:     if (balance) {
    1407:       PetscRandom r;

    1409:       PetscRandomCreate(PETSC_COMM_SELF, &r);
    1410:       PetscRandomSetInterval(r, 0, 2467*size);
    1411:       VecCreate(PETSC_COMM_SELF, &shifts);
    1412:       VecSetSizes(shifts, numShifts, numShifts);
    1413:       VecSetType(shifts, VECSTANDARD);
    1414:       VecSetRandom(shifts, r);
    1415:       PetscRandomDestroy(&r);
    1416:       VecGetArrayRead(shifts, &shift);
    1417:     }

    1419:     PetscMalloc1(nroots, &rootVote);
    1420:     PetscMalloc1(nleaves, &leafVote);
    1421:     /* Point ownership vote: Process with highest rank owns shared points */
    1422:     for (p = 0; p < nleaves; ++p) {
    1423:       if (shiftDebug) {
    1424:         PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Point %D RemotePoint %D Shift %D MyRank %D\n", rank, leaves ? leaves[p] : p, roots[p].index, (PetscInt) PetscRealPart(shift[roots[p].index%numShifts]), (rank + (shift ? (PetscInt) PetscRealPart(shift[roots[p].index%numShifts]) : 0))%size);
    1425:       }
    1426:       /* Either put in a bid or we know we own it */
    1427:       leafVote[p].vote  = (rank + (shift ? (PetscInt) PetscRealPart(shift[roots[p].index%numShifts]) : 0))%size;
    1428:       leafVote[p].rank = rank;
    1429:       leafVote[p].index = p;
    1430:     }
    1431:     for (p = 0; p < nroots; p++) {
    1432:       /* Root must not participate in the reduction, flag so that MAXLOC does not use */
    1433:       rootVote[p].vote  = -3;
    1434:       rootVote[p].rank  = -3;
    1435:       rootVote[p].index = -3;
    1436:     }
    1437:     MPI_Type_contiguous(3, MPIU_INT, &datatype);
    1438:     MPI_Type_commit(&datatype);
    1439:     MPI_Op_create(&MaxLocCarry, 1, &op);
    1440:     PetscSFReduceBegin(migrationSF, datatype, leafVote, rootVote, op);
    1441:     PetscSFReduceEnd(migrationSF, datatype, leafVote, rootVote, op);
    1442:     MPI_Op_free(&op);
    1443:     MPI_Type_free(&datatype);
    1444:     for (p = 0; p < nroots; p++) {
    1445:       rootNodes[p].rank = rootVote[p].rank;
    1446:       rootNodes[p].index = rootVote[p].index;
    1447:     }
    1448:     PetscFree(leafVote);
    1449:     PetscFree(rootVote);
    1450:   } else {
    1451:     for (p = 0; p < nroots; p++) {
    1452:       rootNodes[p].index = -1;
    1453:       rootNodes[p].rank = rank;
    1454:     }
    1455:     for (p = 0; p < nleaves; p++) {
    1456:       /* Write new local id into old location */
    1457:       if (roots[p].rank == rank) {
    1458:         rootNodes[roots[p].index].index = leaves ? leaves[p] : p;
    1459:       }
    1460:     }
    1461:   }
    1462:   PetscSFBcastBegin(migrationSF, MPIU_2INT, rootNodes, leafNodes,MPI_REPLACE);
    1463:   PetscSFBcastEnd(migrationSF, MPIU_2INT, rootNodes, leafNodes,MPI_REPLACE);

    1465:   for (npointLeaves = 0, p = 0; p < nleaves; p++) {
    1466:     if (leafNodes[p].rank != rank) npointLeaves++;
    1467:   }
    1468:   PetscMalloc1(npointLeaves, &pointLocal);
    1469:   PetscMalloc1(npointLeaves, &pointRemote);
    1470:   for (idx = 0, p = 0; p < nleaves; p++) {
    1471:     if (leafNodes[p].rank != rank) {
    1472:       /* Note that pointLocal is automatically sorted as it is sublist of 0, ..., nleaves-1 */
    1473:       pointLocal[idx] = p;
    1474:       pointRemote[idx] = leafNodes[p];
    1475:       idx++;
    1476:     }
    1477:   }
    1478:   if (shift) {
    1479:     VecRestoreArrayRead(shifts, &shift);
    1480:     VecDestroy(&shifts);
    1481:   }
    1482:   if (shiftDebug) {PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), PETSC_STDOUT);}
    1483:   PetscSFCreate(PetscObjectComm((PetscObject) dm), pointSF);
    1484:   PetscSFSetFromOptions(*pointSF);
    1485:   PetscSFSetGraph(*pointSF, nleaves, npointLeaves, pointLocal, PETSC_OWN_POINTER, pointRemote, PETSC_OWN_POINTER);
    1486:   PetscFree2(rootNodes, leafNodes);
    1487:   PetscLogEventEnd(DMPLEX_CreatePointSF,dm,0,0,0);
    1488:   return(0);
    1489: }

    1491: /*@C
    1492:   DMPlexMigrate  - Migrates internal DM data over the supplied star forest

    1494:   Collective on dm

    1496:   Input Parameters:
    1497: + dm       - The source DMPlex object
    1498: . sf       - The star forest communication context describing the migration pattern

    1500:   Output Parameter:
    1501: - targetDM - The target DMPlex object

    1503:   Level: intermediate

    1505: .seealso: DMPlexDistribute(), DMPlexDistributeOverlap()
    1506: @*/
    1507: PetscErrorCode DMPlexMigrate(DM dm, PetscSF sf, DM targetDM)
    1508: {
    1509:   MPI_Comm               comm;
    1510:   PetscInt               dim, cdim, nroots;
    1511:   PetscSF                sfPoint;
    1512:   ISLocalToGlobalMapping ltogMigration;
    1513:   ISLocalToGlobalMapping ltogOriginal = NULL;
    1514:   PetscBool              flg;
    1515:   PetscErrorCode         ierr;

    1519:   PetscLogEventBegin(DMPLEX_Migrate, dm, 0, 0, 0);
    1520:   PetscObjectGetComm((PetscObject) dm, &comm);
    1521:   DMGetDimension(dm, &dim);
    1522:   DMSetDimension(targetDM, dim);
    1523:   DMGetCoordinateDim(dm, &cdim);
    1524:   DMSetCoordinateDim(targetDM, cdim);

    1526:   /* Check for a one-to-all distribution pattern */
    1527:   DMGetPointSF(dm, &sfPoint);
    1528:   PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
    1529:   if (nroots >= 0) {
    1530:     IS        isOriginal;
    1531:     PetscInt  n, size, nleaves;
    1532:     PetscInt  *numbering_orig, *numbering_new;

    1534:     /* Get the original point numbering */
    1535:     DMPlexCreatePointNumbering(dm, &isOriginal);
    1536:     ISLocalToGlobalMappingCreateIS(isOriginal, &ltogOriginal);
    1537:     ISLocalToGlobalMappingGetSize(ltogOriginal, &size);
    1538:     ISLocalToGlobalMappingGetBlockIndices(ltogOriginal, (const PetscInt**)&numbering_orig);
    1539:     /* Convert to positive global numbers */
    1540:     for (n=0; n<size; n++) {if (numbering_orig[n] < 0) numbering_orig[n] = -(numbering_orig[n]+1);}
    1541:     /* Derive the new local-to-global mapping from the old one */
    1542:     PetscSFGetGraph(sf, NULL, &nleaves, NULL, NULL);
    1543:     PetscMalloc1(nleaves, &numbering_new);
    1544:     PetscSFBcastBegin(sf, MPIU_INT, numbering_orig, numbering_new,MPI_REPLACE);
    1545:     PetscSFBcastEnd(sf, MPIU_INT, numbering_orig, numbering_new,MPI_REPLACE);
    1546:     ISLocalToGlobalMappingCreate(comm, 1, nleaves, numbering_new, PETSC_OWN_POINTER, &ltogMigration);
    1547:     ISLocalToGlobalMappingRestoreIndices(ltogOriginal, (const PetscInt**)&numbering_orig);
    1548:     ISDestroy(&isOriginal);
    1549:   } else {
    1550:     /* One-to-all distribution pattern: We can derive LToG from SF */
    1551:     ISLocalToGlobalMappingCreateSF(sf, 0, &ltogMigration);
    1552:   }
    1553:   PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-partition_view", &flg);
    1554:   if (flg) {
    1555:     PetscPrintf(comm, "Point renumbering for DM migration:\n");
    1556:     ISLocalToGlobalMappingView(ltogMigration, NULL);
    1557:   }
    1558:   /* Migrate DM data to target DM */
    1559:   DMPlexDistributeCones(dm, sf, ltogOriginal, ltogMigration, targetDM);
    1560:   DMPlexDistributeLabels(dm, sf, targetDM);
    1561:   DMPlexDistributeCoordinates(dm, sf, targetDM);
    1562:   DMPlexDistributeSetupTree(dm, sf, ltogOriginal, ltogMigration, targetDM);
    1563:   ISLocalToGlobalMappingDestroy(&ltogOriginal);
    1564:   ISLocalToGlobalMappingDestroy(&ltogMigration);
    1565:   PetscLogEventEnd(DMPLEX_Migrate, dm, 0, 0, 0);
    1566:   return(0);
    1567: }

    1569: /*@C
    1570:   DMPlexDistribute - Distributes the mesh and any associated sections.

    1572:   Collective on dm

    1574:   Input Parameters:
    1575: + dm  - The original DMPlex object
    1576: - overlap - The overlap of partitions, 0 is the default

    1578:   Output Parameters:
    1579: + sf - The PetscSF used for point distribution, or NULL if not needed
    1580: - dmParallel - The distributed DMPlex object

    1582:   Note: If the mesh was not distributed, the output dmParallel will be NULL.

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

    1587:   Level: intermediate

    1589: .seealso: DMPlexCreate(), DMSetAdjacency(), DMPlexGetOverlap()
    1590: @*/
    1591: PetscErrorCode DMPlexDistribute(DM dm, PetscInt overlap, PetscSF *sf, DM *dmParallel)
    1592: {
    1593:   MPI_Comm               comm;
    1594:   PetscPartitioner       partitioner;
    1595:   IS                     cellPart;
    1596:   PetscSection           cellPartSection;
    1597:   DM                     dmCoord;
    1598:   DMLabel                lblPartition, lblMigration;
    1599:   PetscSF                sfMigration, sfStratified, sfPoint;
    1600:   PetscBool              flg, balance;
    1601:   PetscMPIInt            rank, size;
    1602:   PetscErrorCode         ierr;

    1610:   if (sf) *sf = NULL;
    1611:   *dmParallel = NULL;
    1612:   PetscObjectGetComm((PetscObject)dm,&comm);
    1613:   MPI_Comm_rank(comm, &rank);
    1614:   MPI_Comm_size(comm, &size);
    1615:   if (size == 1) return(0);

    1617:   PetscLogEventBegin(DMPLEX_Distribute,dm,0,0,0);
    1618:   /* Create cell partition */
    1619:   PetscLogEventBegin(DMPLEX_Partition,dm,0,0,0);
    1620:   PetscSectionCreate(comm, &cellPartSection);
    1621:   DMPlexGetPartitioner(dm, &partitioner);
    1622:   PetscPartitionerDMPlexPartition(partitioner, dm, NULL, cellPartSection, &cellPart);
    1623:   PetscLogEventBegin(DMPLEX_PartSelf,dm,0,0,0);
    1624:   {
    1625:     /* Convert partition to DMLabel */
    1626:     IS             is;
    1627:     PetscHSetI     ht;
    1628:     const PetscInt *points;
    1629:     PetscInt       *iranks;
    1630:     PetscInt       pStart, pEnd, proc, npoints, poff = 0, nranks;

    1632:     DMLabelCreate(PETSC_COMM_SELF, "Point Partition", &lblPartition);
    1633:     /* Preallocate strata */
    1634:     PetscHSetICreate(&ht);
    1635:     PetscSectionGetChart(cellPartSection, &pStart, &pEnd);
    1636:     for (proc = pStart; proc < pEnd; proc++) {
    1637:       PetscSectionGetDof(cellPartSection, proc, &npoints);
    1638:       if (npoints) {PetscHSetIAdd(ht, proc);}
    1639:     }
    1640:     PetscHSetIGetSize(ht, &nranks);
    1641:     PetscMalloc1(nranks, &iranks);
    1642:     PetscHSetIGetElems(ht, &poff, iranks);
    1643:     PetscHSetIDestroy(&ht);
    1644:     DMLabelAddStrata(lblPartition, nranks, iranks);
    1645:     PetscFree(iranks);
    1646:     /* Inline DMPlexPartitionLabelClosure() */
    1647:     ISGetIndices(cellPart, &points);
    1648:     PetscSectionGetChart(cellPartSection, &pStart, &pEnd);
    1649:     for (proc = pStart; proc < pEnd; proc++) {
    1650:       PetscSectionGetDof(cellPartSection, proc, &npoints);
    1651:       if (!npoints) continue;
    1652:       PetscSectionGetOffset(cellPartSection, proc, &poff);
    1653:       DMPlexClosurePoints_Private(dm, npoints, points+poff, &is);
    1654:       DMLabelSetStratumIS(lblPartition, proc, is);
    1655:       ISDestroy(&is);
    1656:     }
    1657:     ISRestoreIndices(cellPart, &points);
    1658:   }
    1659:   PetscLogEventEnd(DMPLEX_PartSelf,dm,0,0,0);

    1661:   DMLabelCreate(PETSC_COMM_SELF, "Point migration", &lblMigration);
    1662:   DMPlexPartitionLabelInvert(dm, lblPartition, NULL, lblMigration);
    1663:   DMPlexPartitionLabelCreateSF(dm, lblMigration, &sfMigration);
    1664:   DMPlexStratifyMigrationSF(dm, sfMigration, &sfStratified);
    1665:   PetscSFDestroy(&sfMigration);
    1666:   sfMigration = sfStratified;
    1667:   PetscSFSetUp(sfMigration);
    1668:   PetscLogEventEnd(DMPLEX_Partition,dm,0,0,0);
    1669:   PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-partition_view", &flg);
    1670:   if (flg) {
    1671:     DMLabelView(lblPartition, PETSC_VIEWER_STDOUT_(comm));
    1672:     PetscSFView(sfMigration, PETSC_VIEWER_STDOUT_(comm));
    1673:   }

    1675:   /* Create non-overlapping parallel DM and migrate internal data */
    1676:   DMPlexCreate(comm, dmParallel);
    1677:   PetscObjectSetName((PetscObject) *dmParallel, "Parallel Mesh");
    1678:   DMPlexMigrate(dm, sfMigration, *dmParallel);

    1680:   /* Build the point SF without overlap */
    1681:   DMPlexGetPartitionBalance(dm, &balance);
    1682:   DMPlexSetPartitionBalance(*dmParallel, balance);
    1683:   DMPlexCreatePointSF(*dmParallel, sfMigration, PETSC_TRUE, &sfPoint);
    1684:   DMSetPointSF(*dmParallel, sfPoint);
    1685:   DMGetCoordinateDM(*dmParallel, &dmCoord);
    1686:   if (dmCoord) {DMSetPointSF(dmCoord, sfPoint);}
    1687:   if (flg) {PetscSFView(sfPoint, NULL);}

    1689:   if (overlap > 0) {
    1690:     DM                 dmOverlap;
    1691:     PetscInt           nroots, nleaves, noldleaves, l;
    1692:     const PetscInt    *oldLeaves;
    1693:     PetscSFNode       *newRemote, *permRemote;
    1694:     const PetscSFNode *oldRemote;
    1695:     PetscSF            sfOverlap, sfOverlapPoint;

    1697:     /* Add the partition overlap to the distributed DM */
    1698:     DMPlexDistributeOverlap(*dmParallel, overlap, &sfOverlap, &dmOverlap);
    1699:     DMDestroy(dmParallel);
    1700:     *dmParallel = dmOverlap;
    1701:     if (flg) {
    1702:       PetscPrintf(comm, "Overlap Migration SF:\n");
    1703:       PetscSFView(sfOverlap, NULL);
    1704:     }

    1706:     /* Re-map the migration SF to establish the full migration pattern */
    1707:     PetscSFGetGraph(sfMigration, &nroots, &noldleaves, &oldLeaves, &oldRemote);
    1708:     PetscSFGetGraph(sfOverlap, NULL, &nleaves, NULL, NULL);
    1709:     PetscMalloc1(nleaves, &newRemote);
    1710:     /* oldRemote: original root point mapping to original leaf point
    1711:        newRemote: original leaf point mapping to overlapped leaf point */
    1712:     if (oldLeaves) {
    1713:       /* After stratification, the migration remotes may not be in root (canonical) order, so we reorder using the leaf numbering */
    1714:       PetscMalloc1(noldleaves, &permRemote);
    1715:       for (l = 0; l < noldleaves; ++l) permRemote[oldLeaves[l]] = oldRemote[l];
    1716:       oldRemote = permRemote;
    1717:     }
    1718:     PetscSFBcastBegin(sfOverlap, MPIU_2INT, oldRemote, newRemote,MPI_REPLACE);
    1719:     PetscSFBcastEnd(sfOverlap, MPIU_2INT, oldRemote, newRemote,MPI_REPLACE);
    1720:     if (oldLeaves) {PetscFree(oldRemote);}
    1721:     PetscSFCreate(comm, &sfOverlapPoint);
    1722:     PetscSFSetGraph(sfOverlapPoint, nroots, nleaves, NULL, PETSC_OWN_POINTER, newRemote, PETSC_OWN_POINTER);
    1723:     PetscSFDestroy(&sfOverlap);
    1724:     PetscSFDestroy(&sfMigration);
    1725:     sfMigration = sfOverlapPoint;
    1726:   }
    1727:   /* Cleanup Partition */
    1728:   DMLabelDestroy(&lblPartition);
    1729:   DMLabelDestroy(&lblMigration);
    1730:   PetscSectionDestroy(&cellPartSection);
    1731:   ISDestroy(&cellPart);
    1732:   /* Copy BC */
    1733:   DMCopyBoundary(dm, *dmParallel);
    1734:   /* Create sfNatural */
    1735:   if (dm->useNatural) {
    1736:     PetscSection section;

    1738:     DMGetLocalSection(dm, &section);
    1739:     DMPlexCreateGlobalToNaturalSF(*dmParallel, section, sfMigration, &(*dmParallel)->sfNatural);
    1740:     DMSetUseNatural(*dmParallel, PETSC_TRUE);
    1741:   }
    1742:   /* Cleanup */
    1743:   if (sf) {*sf = sfMigration;}
    1744:   else    {PetscSFDestroy(&sfMigration);}
    1745:   PetscSFDestroy(&sfPoint);
    1746:   PetscLogEventEnd(DMPLEX_Distribute,dm,0,0,0);
    1747:   return(0);
    1748: }

    1750: /*@C
    1751:   DMPlexDistributeOverlap - Add partition overlap to a distributed non-overlapping DM.

    1753:   Collective on dm

    1755:   Input Parameters:
    1756: + dm  - The non-overlapping distributed DMPlex object
    1757: - overlap - The overlap of partitions (the same on all ranks)

    1759:   Output Parameters:
    1760: + sf - The PetscSF used for point distribution
    1761: - dmOverlap - The overlapping distributed DMPlex object, or NULL

    1763:   Notes:
    1764:   If the mesh was not distributed, the return value is NULL.

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

    1769:   Level: advanced

    1771: .seealso: DMPlexCreate(), DMSetAdjacency(), DMPlexDistribute(), DMPlexCreateOverlapLabel(), DMPlexGetOverlap()
    1772: @*/
    1773: PetscErrorCode DMPlexDistributeOverlap(DM dm, PetscInt overlap, PetscSF *sf, DM *dmOverlap)
    1774: {
    1775:   MPI_Comm               comm;
    1776:   PetscMPIInt            size, rank;
    1777:   PetscSection           rootSection, leafSection;
    1778:   IS                     rootrank, leafrank;
    1779:   DM                     dmCoord;
    1780:   DMLabel                lblOverlap;
    1781:   PetscSF                sfOverlap, sfStratified, sfPoint;
    1782:   PetscErrorCode         ierr;

    1790:   if (sf) *sf = NULL;
    1791:   *dmOverlap  = NULL;
    1792:   PetscObjectGetComm((PetscObject)dm,&comm);
    1793:   MPI_Comm_size(comm, &size);
    1794:   MPI_Comm_rank(comm, &rank);
    1795:   if (size == 1) return(0);

    1797:   PetscLogEventBegin(DMPLEX_DistributeOverlap, dm, 0, 0, 0);
    1798:   /* Compute point overlap with neighbouring processes on the distributed DM */
    1799:   PetscLogEventBegin(DMPLEX_Partition,dm,0,0,0);
    1800:   PetscSectionCreate(comm, &rootSection);
    1801:   PetscSectionCreate(comm, &leafSection);
    1802:   DMPlexDistributeOwnership(dm, rootSection, &rootrank, leafSection, &leafrank);
    1803:   DMPlexCreateOverlapLabel(dm, overlap, rootSection, rootrank, leafSection, leafrank, &lblOverlap);
    1804:   /* Convert overlap label to stratified migration SF */
    1805:   DMPlexPartitionLabelCreateSF(dm, lblOverlap, &sfOverlap);
    1806:   DMPlexStratifyMigrationSF(dm, sfOverlap, &sfStratified);
    1807:   PetscSFDestroy(&sfOverlap);
    1808:   sfOverlap = sfStratified;
    1809:   PetscObjectSetName((PetscObject) sfOverlap, "Overlap SF");
    1810:   PetscSFSetFromOptions(sfOverlap);

    1812:   PetscSectionDestroy(&rootSection);
    1813:   PetscSectionDestroy(&leafSection);
    1814:   ISDestroy(&rootrank);
    1815:   ISDestroy(&leafrank);
    1816:   PetscLogEventEnd(DMPLEX_Partition,dm,0,0,0);

    1818:   /* Build the overlapping DM */
    1819:   DMPlexCreate(comm, dmOverlap);
    1820:   PetscObjectSetName((PetscObject) *dmOverlap, "Parallel Mesh");
    1821:   DMPlexMigrate(dm, sfOverlap, *dmOverlap);
    1822:   /* Store the overlap in the new DM */
    1823:   ((DM_Plex*)(*dmOverlap)->data)->overlap = overlap + ((DM_Plex*)dm->data)->overlap;
    1824:   /* Build the new point SF */
    1825:   DMPlexCreatePointSF(*dmOverlap, sfOverlap, PETSC_FALSE, &sfPoint);
    1826:   DMSetPointSF(*dmOverlap, sfPoint);
    1827:   DMGetCoordinateDM(*dmOverlap, &dmCoord);
    1828:   if (dmCoord) {DMSetPointSF(dmCoord, sfPoint);}
    1829:   PetscSFDestroy(&sfPoint);
    1830:   /* Cleanup overlap partition */
    1831:   DMLabelDestroy(&lblOverlap);
    1832:   if (sf) *sf = sfOverlap;
    1833:   else    {PetscSFDestroy(&sfOverlap);}
    1834:   PetscLogEventEnd(DMPLEX_DistributeOverlap, dm, 0, 0, 0);
    1835:   return(0);
    1836: }

    1838: PetscErrorCode DMPlexGetOverlap_Plex(DM dm, PetscInt *overlap)
    1839: {
    1840:   DM_Plex        *mesh  = (DM_Plex*) dm->data;

    1843:   *overlap = mesh->overlap;
    1844:   return(0);
    1845: }

    1847: /*@
    1848:   DMPlexGetOverlap - Get the DMPlex partition overlap.

    1850:   Not collective

    1852:   Input Parameter:
    1853: . dm - The DM

    1855:   Output Parameter:
    1856: . overlap - The overlap of this DM

    1858:   Level: intermediate

    1860: .seealso: DMPlexDistribute(), DMPlexDistributeOverlap(), DMPlexCreateOverlapLabel()
    1861: @*/
    1862: PetscErrorCode DMPlexGetOverlap(DM dm, PetscInt *overlap)
    1863: {
    1864:   PetscErrorCode     ierr;

    1868:   PetscUseMethod(dm,"DMPlexGetOverlap_C",(DM,PetscInt*),(dm,overlap));
    1869:   return(0);
    1870: }

    1873: /*@C
    1874:   DMPlexGetGatherDM - Get a copy of the DMPlex that gathers all points on the
    1875:   root process of the original's communicator.

    1877:   Collective on dm

    1879:   Input Parameter:
    1880: . dm - the original DMPlex object

    1882:   Output Parameters:
    1883: + sf - the PetscSF used for point distribution (optional)
    1884: - gatherMesh - the gathered DM object, or NULL

    1886:   Level: intermediate

    1888: .seealso: DMPlexDistribute(), DMPlexGetRedundantDM()
    1889: @*/
    1890: PetscErrorCode DMPlexGetGatherDM(DM dm, PetscSF *sf, DM *gatherMesh)
    1891: {
    1892:   MPI_Comm       comm;
    1893:   PetscMPIInt    size;
    1894:   PetscPartitioner oldPart, gatherPart;

    1900:   *gatherMesh = NULL;
    1901:   if (sf) *sf = NULL;
    1902:   comm = PetscObjectComm((PetscObject)dm);
    1903:   MPI_Comm_size(comm,&size);
    1904:   if (size == 1) return(0);
    1905:   DMPlexGetPartitioner(dm,&oldPart);
    1906:   PetscObjectReference((PetscObject)oldPart);
    1907:   PetscPartitionerCreate(comm,&gatherPart);
    1908:   PetscPartitionerSetType(gatherPart,PETSCPARTITIONERGATHER);
    1909:   DMPlexSetPartitioner(dm,gatherPart);
    1910:   DMPlexDistribute(dm,0,sf,gatherMesh);

    1912:   DMPlexSetPartitioner(dm,oldPart);
    1913:   PetscPartitionerDestroy(&gatherPart);
    1914:   PetscPartitionerDestroy(&oldPart);
    1915:   return(0);
    1916: }

    1918: /*@C
    1919:   DMPlexGetRedundantDM - Get a copy of the DMPlex that is completely copied on each process.

    1921:   Collective on dm

    1923:   Input Parameter:
    1924: . dm - the original DMPlex object

    1926:   Output Parameters:
    1927: + sf - the PetscSF used for point distribution (optional)
    1928: - redundantMesh - the redundant DM object, or NULL

    1930:   Level: intermediate

    1932: .seealso: DMPlexDistribute(), DMPlexGetGatherDM()
    1933: @*/
    1934: PetscErrorCode DMPlexGetRedundantDM(DM dm, PetscSF *sf, DM *redundantMesh)
    1935: {
    1936:   MPI_Comm       comm;
    1937:   PetscMPIInt    size, rank;
    1938:   PetscInt       pStart, pEnd, p;
    1939:   PetscInt       numPoints = -1;
    1940:   PetscSF        migrationSF, sfPoint, gatherSF;
    1941:   DM             gatherDM, dmCoord;
    1942:   PetscSFNode    *points;

    1948:   *redundantMesh = NULL;
    1949:   comm = PetscObjectComm((PetscObject)dm);
    1950:   MPI_Comm_size(comm,&size);
    1951:   if (size == 1) {
    1952:     PetscObjectReference((PetscObject) dm);
    1953:     *redundantMesh = dm;
    1954:     if (sf) *sf = NULL;
    1955:     return(0);
    1956:   }
    1957:   DMPlexGetGatherDM(dm,&gatherSF,&gatherDM);
    1958:   if (!gatherDM) return(0);
    1959:   MPI_Comm_rank(comm,&rank);
    1960:   DMPlexGetChart(gatherDM,&pStart,&pEnd);
    1961:   numPoints = pEnd - pStart;
    1962:   MPI_Bcast(&numPoints,1,MPIU_INT,0,comm);
    1963:   PetscMalloc1(numPoints,&points);
    1964:   PetscSFCreate(comm,&migrationSF);
    1965:   for (p = 0; p < numPoints; p++) {
    1966:     points[p].index = p;
    1967:     points[p].rank  = 0;
    1968:   }
    1969:   PetscSFSetGraph(migrationSF,pEnd-pStart,numPoints,NULL,PETSC_OWN_POINTER,points,PETSC_OWN_POINTER);
    1970:   DMPlexCreate(comm, redundantMesh);
    1971:   PetscObjectSetName((PetscObject) *redundantMesh, "Redundant Mesh");
    1972:   DMPlexMigrate(gatherDM, migrationSF, *redundantMesh);
    1973:   DMPlexCreatePointSF(*redundantMesh, migrationSF, PETSC_FALSE, &sfPoint);
    1974:   DMSetPointSF(*redundantMesh, sfPoint);
    1975:   DMGetCoordinateDM(*redundantMesh, &dmCoord);
    1976:   if (dmCoord) {DMSetPointSF(dmCoord, sfPoint);}
    1977:   PetscSFDestroy(&sfPoint);
    1978:   if (sf) {
    1979:     PetscSF tsf;

    1981:     PetscSFCompose(gatherSF,migrationSF,&tsf);
    1982:     DMPlexStratifyMigrationSF(dm, tsf, sf);
    1983:     PetscSFDestroy(&tsf);
    1984:   }
    1985:   PetscSFDestroy(&migrationSF);
    1986:   PetscSFDestroy(&gatherSF);
    1987:   DMDestroy(&gatherDM);
    1988:   return(0);
    1989: }

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

    1994:   Collective

    1996:   Input Parameter:
    1997: . dm      - The DM object

    1999:   Output Parameter:
    2000: . distributed - Flag whether the DM is distributed

    2002:   Level: intermediate

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

    2009: .seealso: DMPlexDistribute(), DMPlexGetOverlap(), DMPlexIsInterpolated()
    2010: @*/
    2011: PetscErrorCode DMPlexIsDistributed(DM dm, PetscBool *distributed)
    2012: {
    2013:   PetscInt          pStart, pEnd, count;
    2014:   MPI_Comm          comm;
    2015:   PetscErrorCode    ierr;

    2020:   PetscObjectGetComm((PetscObject)dm,&comm);
    2021:   DMPlexGetChart(dm, &pStart, &pEnd);
    2022:   count = !!(pEnd - pStart);
    2023:   MPI_Allreduce(MPI_IN_PLACE, &count, 1, MPIU_INT, MPI_SUM, comm);
    2024:   *distributed = count > 1 ? PETSC_TRUE : PETSC_FALSE;
    2025:   return(0);
    2026: }