Actual source code: plexnatural.c

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

      3: /*@
      4:   DMPlexSetMigrationSF - Sets the `PetscSF` for migrating from a parent `DM` into this `DM`

      6:   Logically Collective

      8:   Input Parameters:
      9: + dm          - The `DM`
     10: - migrationSF - The `PetscSF`

     12:   Level: intermediate

     14:   Note:
     15:   It is necessary to call this in order to have `DMCreateSubDM()` or `DMCreateSuperDM()` build the Global-To-Natural map

     17: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateMigrationSF()`, `DMPlexGetMigrationSF()`
     18: @*/
     19: PetscErrorCode DMPlexSetMigrationSF(DM dm, PetscSF migrationSF)
     20: {
     21:   PetscFunctionBegin;
     24:   PetscCall(PetscObjectReference((PetscObject)migrationSF));
     25:   PetscCall(PetscSFDestroy(&dm->sfMigration));
     26:   dm->sfMigration = migrationSF;
     27:   PetscFunctionReturn(PETSC_SUCCESS);
     28: }

     30: /*@
     31:   DMPlexGetMigrationSF - Gets the `PetscSF` for migrating from a parent `DM` into this `DM`

     33:   Note Collective

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

     38:   Output Parameter:
     39: . migrationSF - The `PetscSF`

     41:   Level: intermediate

     43: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateMigrationSF()`, `DMPlexSetMigrationSF`
     44: @*/
     45: PetscErrorCode DMPlexGetMigrationSF(DM dm, PetscSF *migrationSF)
     46: {
     47:   PetscFunctionBegin;
     48:   *migrationSF = dm->sfMigration;
     49:   PetscFunctionReturn(PETSC_SUCCESS);
     50: }

     52: /*@
     53:   DMPlexSetGlobalToNaturalSF - Sets the `PetscSF` for mapping Global `Vec` to the Natural `Vec`

     55:   Input Parameters:
     56: + dm        - The `DM`
     57: - naturalSF - The `PetscSF`

     59:   Level: intermediate

     61: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateGlobalToNaturalSF()`, `DMPlexGetGlobaltoNaturalSF()`
     62: @*/
     63: PetscErrorCode DMPlexSetGlobalToNaturalSF(DM dm, PetscSF naturalSF)
     64: {
     65:   PetscFunctionBegin;
     66:   dm->sfNatural = naturalSF;
     67:   PetscCall(PetscObjectReference((PetscObject)naturalSF));
     68:   dm->useNatural = naturalSF ? PETSC_TRUE : PETSC_FALSE;
     69:   PetscFunctionReturn(PETSC_SUCCESS);
     70: }

     72: /*@
     73:   DMPlexGetGlobalToNaturalSF - Gets the `PetscSF` for mapping Global `Vec` to the Natural `Vec`

     75:   Input Parameter:
     76: . dm - The `DM`

     78:   Output Parameter:
     79: . naturalSF - The `PetscSF`

     81:   Level: intermediate

     83: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateGlobalToNaturalSF()`, `DMPlexSetGlobaltoNaturalSF`
     84: @*/
     85: PetscErrorCode DMPlexGetGlobalToNaturalSF(DM dm, PetscSF *naturalSF)
     86: {
     87:   PetscFunctionBegin;
     88:   *naturalSF = dm->sfNatural;
     89:   PetscFunctionReturn(PETSC_SUCCESS);
     90: }

     92: /*@
     93:   DMPlexCreateGlobalToNaturalSF - Creates the `PetscSF` for mapping Global `Vec` to the Natural `Vec`

     95:   Input Parameters:
     96: + dm          - The redistributed `DM`
     97: . section     - The local `PetscSection` describing the `Vec` before the mesh was distributed, or `NULL` if not available
     98: - sfMigration - The `PetscSF` used to distribute the mesh, or `NULL` if it cannot be computed

    100:   Output Parameter:
    101: . sfNatural - `PetscSF` for mapping the `Vec` in PETSc ordering to the canonical ordering

    103:   Level: intermediate

    105:   Note:
    106:   This is not typically called by the user.

    108: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `PetscSection`, `DMPlexDistribute()`, `DMPlexDistributeField()`
    109:  @*/
    110: PetscErrorCode DMPlexCreateGlobalToNaturalSF(DM dm, PetscSection section, PetscSF sfMigration, PetscSF *sfNatural)
    111: {
    112:   MPI_Comm     comm;
    113:   PetscSF      sf, sfEmbed, sfField;
    114:   PetscSection gSection, sectionDist, gLocSection;
    115:   PetscInt    *spoints, *remoteOffsets;
    116:   PetscInt     ssize, pStart, pEnd, p, localSize, maxStorageSize;
    117:   PetscBool    destroyFlag = PETSC_FALSE, debug = PETSC_FALSE;

    119:   PetscFunctionBegin;
    120:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    121:   if (!sfMigration) {
    122:     /* If sfMigration is missing, sfNatural cannot be computed and is set to NULL */
    123:     *sfNatural = NULL;
    124:     PetscFunctionReturn(PETSC_SUCCESS);
    125:   } else if (!section) {
    126:     /* If the sequential section is not provided (NULL), it is reconstructed from the parallel section */
    127:     PetscSF      sfMigrationInv;
    128:     PetscSection localSection;

    130:     PetscCall(DMGetLocalSection(dm, &localSection));
    131:     PetscCall(PetscSFCreateInverseSF(sfMigration, &sfMigrationInv));
    132:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section));
    133:     PetscCall(PetscSFDistributeSection(sfMigrationInv, localSection, NULL, section));
    134:     PetscCall(PetscSFDestroy(&sfMigrationInv));
    135:     destroyFlag = PETSC_TRUE;
    136:   }
    137:   if (debug) PetscCall(PetscSFView(sfMigration, NULL));
    138:   /* Create a new section from distributing the original section */
    139:   PetscCall(PetscSectionCreate(comm, &sectionDist));
    140:   PetscCall(PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist));
    141:   PetscCall(PetscObjectSetName((PetscObject)sectionDist, "Migrated Section"));
    142:   if (debug) PetscCall(PetscSectionView(sectionDist, NULL));
    143:   PetscCall(DMSetLocalSection(dm, sectionDist));
    144:   /* If a sequential section is provided but no dof is affected, sfNatural cannot be computed and is set to NULL */
    145:   PetscCall(PetscSectionGetStorageSize(sectionDist, &localSize));
    146:   PetscCallMPI(MPIU_Allreduce(&localSize, &maxStorageSize, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
    147:   if (maxStorageSize) {
    148:     const PetscInt *leaves;
    149:     PetscInt       *sortleaves, *indices;
    150:     PetscInt        Nl;

    152:     /* Get a pruned version of migration SF */
    153:     PetscCall(DMGetGlobalSection(dm, &gSection));
    154:     if (debug) PetscCall(PetscSectionView(gSection, NULL));
    155:     PetscCall(PetscSFGetGraph(sfMigration, NULL, &Nl, &leaves, NULL));
    156:     PetscCall(PetscSectionGetChart(gSection, &pStart, &pEnd));
    157:     for (p = pStart, ssize = 0; p < pEnd; ++p) {
    158:       PetscInt dof, off;

    160:       PetscCall(PetscSectionGetDof(gSection, p, &dof));
    161:       PetscCall(PetscSectionGetOffset(gSection, p, &off));
    162:       if ((dof > 0) && (off >= 0)) ++ssize;
    163:     }
    164:     PetscCall(PetscMalloc3(ssize, &spoints, Nl, &sortleaves, Nl, &indices));
    165:     for (p = 0; p < Nl; ++p) {
    166:       sortleaves[p] = leaves ? leaves[p] : p;
    167:       indices[p]    = p;
    168:     }
    169:     PetscCall(PetscSortIntWithArray(Nl, sortleaves, indices));
    170:     for (p = pStart, ssize = 0; p < pEnd; ++p) {
    171:       PetscInt dof, off, loc;

    173:       PetscCall(PetscSectionGetDof(gSection, p, &dof));
    174:       PetscCall(PetscSectionGetOffset(gSection, p, &off));
    175:       if ((dof > 0) && (off >= 0)) {
    176:         PetscCall(PetscFindInt(p, Nl, sortleaves, &loc));
    177:         PetscCheck(loc >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " with nonzero dof is not a leaf of the migration SF", p);
    178:         spoints[ssize++] = indices[loc];
    179:       }
    180:     }
    181:     PetscCall(PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed));
    182:     PetscCall(PetscObjectSetName((PetscObject)sfEmbed, "Embedded SF"));
    183:     PetscCall(PetscFree3(spoints, sortleaves, indices));
    184:     if (debug) PetscCall(PetscSFView(sfEmbed, NULL));
    185:     /* Create the SF associated with this section
    186:          Roots are natural dofs, leaves are global dofs */
    187:     PetscCall(DMGetPointSF(dm, &sf));
    188:     PetscCall(PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, &gLocSection));
    189:     PetscCall(PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField));
    190:     PetscCall(PetscSFDestroy(&sfEmbed));
    191:     PetscCall(PetscSectionDestroy(&gLocSection));
    192:     PetscCall(PetscObjectSetName((PetscObject)sfField, "Natural-to-Global SF"));
    193:     if (debug) PetscCall(PetscSFView(sfField, NULL));
    194:     /* Invert the field SF
    195:          Roots are global dofs, leaves are natural dofs */
    196:     PetscCall(PetscSFCreateInverseSF(sfField, sfNatural));
    197:     PetscCall(PetscObjectSetName((PetscObject)*sfNatural, "Global-to-Natural SF"));
    198:     PetscCall(PetscObjectViewFromOptions((PetscObject)*sfNatural, NULL, "-globaltonatural_sf_view"));
    199:     /* Clean up */
    200:     PetscCall(PetscSFDestroy(&sfField));
    201:   } else {
    202:     *sfNatural = NULL;
    203:   }
    204:   PetscCall(PetscSectionDestroy(&sectionDist));
    205:   PetscCall(PetscFree(remoteOffsets));
    206:   if (destroyFlag) PetscCall(PetscSectionDestroy(&section));
    207:   PetscFunctionReturn(PETSC_SUCCESS);
    208: }

    210: /*@
    211:   DMPlexMigrateGlobalToNaturalSF - Migrates the input `sfNatural` based on sfMigration

    213:   Input Parameters:
    214: + dmOld        - The original `DM`
    215: . dmNew        - The `DM` to be migrated to
    216: . sfNaturalOld - The sfNatural for the `dmOld`
    217: - sfMigration  - The `PetscSF` used to distribute the mesh, or `NULL` if it cannot be computed

    219:   Output Parameter:
    220: . sfNaturalNew - `PetscSF` for mapping the `Vec` in PETSc ordering to the canonical ordering

    222:   Level: intermediate

    224:   Notes:
    225:   `sfNaturalOld` maps from the old Global section (roots) to the natural Vec layout (leaves, may or may not be described by a PetscSection).
    226:   `DMPlexMigrateGlobalToNaturalSF` creates an SF to map from the old global section to the new global section (generated from `sfMigration`).
    227:   That SF is then composed with the `sfNaturalOld` to generate `sfNaturalNew`.
    228:   This also distributes and sets the local section for `dmNew`.

    230:   This is not typically called by the user.

    232: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `PetscSection`, `DMPlexDistribute()`, `DMPlexDistributeField()`
    233:  @*/
    234: PetscErrorCode DMPlexMigrateGlobalToNaturalSF(DM dmOld, DM dmNew, PetscSF sfNaturalOld, PetscSF sfMigration, PetscSF *sfNaturalNew)
    235: {
    236:   MPI_Comm     comm;
    237:   PetscSection oldGlobalSection;
    238:   PetscSection newGlobalSection;
    239:   PetscInt    *remoteOffsets;
    240:   PetscBool    debug = PETSC_FALSE;

    242:   PetscFunctionBegin;
    243:   PetscCall(PetscObjectGetComm((PetscObject)dmNew, &comm));
    244:   if (!sfMigration) {
    245:     /* If sfMigration is missing, sfNatural cannot be computed and is set to NULL */
    246:     *sfNaturalNew = NULL;
    247:     PetscFunctionReturn(PETSC_SUCCESS);
    248:   }
    249:   if (debug) PetscCall(PetscSFView(sfMigration, NULL));

    251:   { // Create oldGlobalSection and newGlobalSection *with* localOffsets
    252:     PetscSection oldLocalSection, newLocalSection;
    253:     PetscSF      pointSF;

    255:     PetscCall(DMGetLocalSection(dmOld, &oldLocalSection));
    256:     PetscCall(DMGetPointSF(dmOld, &pointSF));
    257:     PetscCall(PetscSectionCreateGlobalSection(oldLocalSection, pointSF, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, &oldGlobalSection));

    259:     PetscCall(PetscSectionCreate(comm, &newLocalSection));
    260:     PetscCall(PetscSFDistributeSection(sfMigration, oldLocalSection, NULL, newLocalSection));
    261:     PetscCall(DMSetLocalSection(dmNew, newLocalSection));

    263:     PetscCall(PetscSectionCreate(comm, &newGlobalSection));
    264:     PetscCall(DMGetPointSF(dmNew, &pointSF));
    265:     PetscCall(PetscSectionCreateGlobalSection(newLocalSection, pointSF, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, &newGlobalSection));

    267:     PetscCall(PetscObjectSetName((PetscObject)oldLocalSection, "Old Local Section"));
    268:     if (debug) PetscCall(PetscSectionView(oldLocalSection, NULL));
    269:     PetscCall(PetscObjectSetName((PetscObject)oldGlobalSection, "Old Global Section"));
    270:     if (debug) PetscCall(PetscSectionView(oldGlobalSection, NULL));
    271:     PetscCall(PetscObjectSetName((PetscObject)newLocalSection, "New Local Section"));
    272:     if (debug) PetscCall(PetscSectionView(newLocalSection, NULL));
    273:     PetscCall(PetscObjectSetName((PetscObject)newGlobalSection, "New Global Section"));
    274:     if (debug) PetscCall(PetscSectionView(newGlobalSection, NULL));
    275:     PetscCall(PetscSectionDestroy(&newLocalSection));
    276:   }

    278:   { // Create remoteOffsets array, mapping the oldGlobalSection offsets to the local points (according to sfMigration)
    279:     PetscInt lpStart, lpEnd, rpStart, rpEnd;

    281:     PetscCall(PetscSectionGetChart(oldGlobalSection, &rpStart, &rpEnd));
    282:     PetscCall(PetscSectionGetChart(newGlobalSection, &lpStart, &lpEnd));

    284:     // in `PetscSFDistributeSection` (where this is taken from), it possibly makes a new embedded SF. Should possibly do that here?
    285:     PetscCall(PetscMalloc1(lpEnd - lpStart, &remoteOffsets));
    286:     PetscCall(PetscSFBcastBegin(sfMigration, MPIU_INT, PetscSafePointerPlusOffset(oldGlobalSection->atlasOff, -rpStart), PetscSafePointerPlusOffset(remoteOffsets, -lpStart), MPI_REPLACE));
    287:     PetscCall(PetscSFBcastEnd(sfMigration, MPIU_INT, PetscSafePointerPlusOffset(oldGlobalSection->atlasOff, -rpStart), PetscSafePointerPlusOffset(remoteOffsets, -lpStart), MPI_REPLACE));
    288:     if (debug) {
    289:       PetscViewer viewer;

    291:       PetscCall(PetscPrintf(comm, "Remote Offsets:\n"));
    292:       PetscCall(PetscViewerASCIIGetStdout(comm, &viewer));
    293:       PetscCall(PetscIntView(lpEnd - lpStart, remoteOffsets, viewer));
    294:     }
    295:   }

    297:   { // Create SF from oldGlobalSection to newGlobalSection and compose with sfNaturalOld
    298:     PetscSF oldglob_to_newglob_sf, newglob_to_oldglob_sf;

    300:     PetscCall(PetscSFCreateSectionSF(sfMigration, oldGlobalSection, remoteOffsets, newGlobalSection, &oldglob_to_newglob_sf));
    301:     PetscCall(PetscObjectSetName((PetscObject)oldglob_to_newglob_sf, "OldGlobal-to-NewGlobal SF"));
    302:     if (debug) PetscCall(PetscSFView(oldglob_to_newglob_sf, NULL));

    304:     PetscCall(PetscSFCreateInverseSF(oldglob_to_newglob_sf, &newglob_to_oldglob_sf));
    305:     PetscCall(PetscObjectSetName((PetscObject)newglob_to_oldglob_sf, "NewGlobal-to-OldGlobal SF"));
    306:     PetscCall(PetscObjectViewFromOptions((PetscObject)newglob_to_oldglob_sf, (PetscObject)dmOld, "-natural_migrate_sf_view"));
    307:     PetscCall(PetscSFCompose(newglob_to_oldglob_sf, sfNaturalOld, sfNaturalNew));

    309:     PetscCall(PetscSFDestroy(&oldglob_to_newglob_sf));
    310:     PetscCall(PetscSFDestroy(&newglob_to_oldglob_sf));
    311:   }

    313:   PetscCall(PetscSectionDestroy(&oldGlobalSection));
    314:   PetscCall(PetscSectionDestroy(&newGlobalSection));
    315:   PetscCall(PetscFree(remoteOffsets));
    316:   PetscFunctionReturn(PETSC_SUCCESS);
    317: }

    319: /*@
    320:   DMPlexGlobalToNaturalBegin - Rearranges a global `Vec` in the natural order.

    322:   Collective

    324:   Input Parameters:
    325: + dm - The distributed `DMPLEX`
    326: - gv - The global `Vec`

    328:   Output Parameter:
    329: . nv - `Vec` in the canonical ordering distributed over all processors associated with `gv`

    331:   Level: intermediate

    333:   Note:
    334:   The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.

    336: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()`
    337: @*/
    338: PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
    339: {
    340:   const PetscScalar *inarray;
    341:   PetscScalar       *outarray;
    342:   MPI_Comm           comm;
    343:   PetscMPIInt        size;

    345:   PetscFunctionBegin;
    346:   PetscCall(PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin, dm, 0, 0, 0));
    347:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    348:   PetscCallMPI(MPI_Comm_size(comm, &size));
    349:   if (dm->sfNatural) {
    350:     if (PetscDefined(USE_DEBUG)) {
    351:       PetscSection gs;
    352:       PetscInt     Nl, n;

    354:       PetscCall(PetscSFGetGraph(dm->sfNatural, NULL, &Nl, NULL, NULL));
    355:       PetscCall(VecGetLocalSize(nv, &n));
    356:       PetscCheck(n == Nl, comm, PETSC_ERR_ARG_INCOMP, "Natural vector local size %" PetscInt_FMT " != %" PetscInt_FMT " local size of natural section", n, Nl);

    358:       PetscCall(DMGetGlobalSection(dm, &gs));
    359:       PetscCall(PetscSectionGetConstrainedStorageSize(gs, &Nl));
    360:       PetscCall(VecGetLocalSize(gv, &n));
    361:       PetscCheck(n == Nl, comm, PETSC_ERR_ARG_INCOMP, "Global vector local size %" PetscInt_FMT " != %" PetscInt_FMT " local size of global section", n, Nl);
    362:     }
    363:     PetscCall(VecGetArray(nv, &outarray));
    364:     PetscCall(VecGetArrayRead(gv, &inarray));
    365:     PetscCall(PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_REPLACE));
    366:     PetscCall(VecRestoreArrayRead(gv, &inarray));
    367:     PetscCall(VecRestoreArray(nv, &outarray));
    368:   } else if (size == 1) {
    369:     PetscCall(VecCopy(gv, nv));
    370:   } else {
    371:     PetscCheck(!dm->useNatural, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "DM global to natural SF not present. If DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.");
    372:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
    373:   }
    374:   PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin, dm, 0, 0, 0));
    375:   PetscFunctionReturn(PETSC_SUCCESS);
    376: }

    378: /*@
    379:   DMPlexGlobalToNaturalEnd - Rearranges a global `Vec` in the natural order.

    381:   Collective

    383:   Input Parameters:
    384: + dm - The distributed `DMPLEX`
    385: - gv - The global `Vec`

    387:   Output Parameter:
    388: . nv - The natural `Vec`

    390:   Level: intermediate

    392:   Note:
    393:   The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.

    395: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
    396:  @*/
    397: PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
    398: {
    399:   const PetscScalar *inarray;
    400:   PetscScalar       *outarray;
    401:   PetscMPIInt        size;

    403:   PetscFunctionBegin;
    404:   PetscCall(PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd, dm, 0, 0, 0));
    405:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
    406:   if (dm->sfNatural) {
    407:     PetscCall(VecGetArrayRead(gv, &inarray));
    408:     PetscCall(VecGetArray(nv, &outarray));
    409:     PetscCall(PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_REPLACE));
    410:     PetscCall(VecRestoreArrayRead(gv, &inarray));
    411:     PetscCall(VecRestoreArray(nv, &outarray));
    412:   } else if (size == 1) {
    413:   } else {
    414:     PetscCheck(!dm->useNatural, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "DM global to natural SF not present. If DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.");
    415:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
    416:   }
    417:   PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd, dm, 0, 0, 0));
    418:   PetscFunctionReturn(PETSC_SUCCESS);
    419: }

    421: /*@
    422:   DMPlexNaturalToGlobalBegin - Rearranges a `Vec` in the natural order to the Global order.

    424:   Collective

    426:   Input Parameters:
    427: + dm - The distributed `DMPLEX`
    428: - nv - The natural `Vec`

    430:   Output Parameter:
    431: . gv - The global `Vec`

    433:   Level: intermediate

    435:   Note:
    436:   The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.

    438: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexGlobalToNaturalEnd()`
    439: @*/
    440: PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
    441: {
    442:   const PetscScalar *inarray;
    443:   PetscScalar       *outarray;
    444:   PetscMPIInt        size;

    446:   PetscFunctionBegin;
    447:   PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin, dm, 0, 0, 0));
    448:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
    449:   if (dm->sfNatural) {
    450:     /* We only have access to the SF that goes from Global to Natural.
    451:        Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM.
    452:        Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */
    453:     PetscCall(VecZeroEntries(gv));
    454:     PetscCall(VecGetArray(gv, &outarray));
    455:     PetscCall(VecGetArrayRead(nv, &inarray));
    456:     PetscCall(PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_SUM));
    457:     PetscCall(VecRestoreArrayRead(nv, &inarray));
    458:     PetscCall(VecRestoreArray(gv, &outarray));
    459:   } else if (size == 1) {
    460:     PetscCall(VecCopy(nv, gv));
    461:   } else {
    462:     PetscCheck(!dm->useNatural, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "DM global to natural SF not present. If DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.");
    463:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
    464:   }
    465:   PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin, dm, 0, 0, 0));
    466:   PetscFunctionReturn(PETSC_SUCCESS);
    467: }

    469: /*@
    470:   DMPlexNaturalToGlobalEnd - Rearranges a `Vec` in the natural order to the Global order.

    472:   Collective

    474:   Input Parameters:
    475: + dm - The distributed `DMPLEX`
    476: - nv - The natural `Vec`

    478:   Output Parameter:
    479: . gv - The global `Vec`

    481:   Level: intermediate

    483:   Note:
    484:   The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.

    486: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
    487:  @*/
    488: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
    489: {
    490:   const PetscScalar *inarray;
    491:   PetscScalar       *outarray;
    492:   PetscMPIInt        size;

    494:   PetscFunctionBegin;
    495:   PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0));
    496:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
    497:   if (dm->sfNatural) {
    498:     PetscCall(VecGetArrayRead(nv, &inarray));
    499:     PetscCall(VecGetArray(gv, &outarray));
    500:     PetscCall(PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *)inarray, outarray, MPI_SUM));
    501:     PetscCall(VecRestoreArrayRead(nv, &inarray));
    502:     PetscCall(VecRestoreArray(gv, &outarray));
    503:   } else if (size == 1) {
    504:   } else {
    505:     PetscCheck(!dm->useNatural, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "DM global to natural SF not present. If DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.");
    506:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
    507:   }
    508:   PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd, dm, 0, 0, 0));
    509:   PetscFunctionReturn(PETSC_SUCCESS);
    510: }

    512: /*@
    513:   DMPlexCreateNaturalVector - Provide a `Vec` capable of holding the natural ordering and distribution.

    515:   Collective

    517:   Input Parameter:
    518: . dm - The distributed `DMPLEX`

    520:   Output Parameter:
    521: . nv - The natural `Vec`

    523:   Level: intermediate

    525:   Note:
    526:   The user must call `DMSetUseNatural`(dm, `PETSC_TRUE`) before `DMPlexDistribute()`.

    528: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `Vec`, `DMPlexDistribute()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()`
    529:  @*/
    530: PetscErrorCode DMPlexCreateNaturalVector(DM dm, Vec *nv)
    531: {
    532:   PetscMPIInt size;

    534:   PetscFunctionBegin;
    535:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
    536:   if (dm->sfNatural) {
    537:     PetscInt nleaves, bs, maxbs;
    538:     Vec      v;

    540:     /*
    541:       Setting the natural vector block size.
    542:       We can't get it from a global vector because of constraints, and the block size in the local vector
    543:       may be inconsistent across processes, typically when some local vectors have size 0, their block size is set to 1
    544:     */
    545:     PetscCall(DMGetLocalVector(dm, &v));
    546:     PetscCall(VecGetBlockSize(v, &bs));
    547:     PetscCallMPI(MPIU_Allreduce(&bs, &maxbs, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
    548:     if (bs == 1 && maxbs > 1) bs = maxbs;
    549:     PetscCall(DMRestoreLocalVector(dm, &v));

    551:     PetscCall(PetscSFGetGraph(dm->sfNatural, NULL, &nleaves, NULL, NULL));
    552:     PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), nv));
    553:     PetscCall(VecSetSizes(*nv, nleaves, PETSC_DETERMINE));
    554:     PetscCall(VecSetBlockSize(*nv, bs));
    555:     PetscCall(VecSetType(*nv, dm->vectype));
    556:     PetscCall(VecSetDM(*nv, dm));
    557:   } else if (size == 1) {
    558:     PetscCall(DMCreateLocalVector(dm, nv));
    559:   } else {
    560:     PetscCheck(!dm->useNatural, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "DM global to natural SF not present. If DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.");
    561:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created. You must call DMSetUseNatural() before DMPlexDistribute().");
    562:   }
    563:   PetscFunctionReturn(PETSC_SUCCESS);
    564: }