Actual source code: stagutils.c

      1: /* Additional functions in the DMStag API, which are not part of the general DM API. */
      2: #include <petsc/private/dmstagimpl.h>
      3: #include <petscdmproduct.h>

      5: PetscErrorCode DMRestrictHook_Coordinates(DM, DM, void *);

      7: /*@C
      8:   DMStagGetBoundaryTypes - get boundary types

     10:   Not Collective

     12:   Input Parameter:
     13: . dm - the `DMSTAG` object

     15:   Output Parameters:
     16: + boundaryTypeX - boundary type for x direction
     17: . boundaryTypeY - boundary type for y direction, not set for one dimensional problems
     18: - boundaryTypeZ - boundary type for z direction, not set for one and two dimensional problems

     20:   Level: intermediate

     22: .seealso: [](ch_stag), `DMSTAG`, `DMBoundaryType`
     23: @*/
     24: PetscErrorCode DMStagGetBoundaryTypes(DM dm, DMBoundaryType *boundaryTypeX, DMBoundaryType *boundaryTypeY, DMBoundaryType *boundaryTypeZ)
     25: {
     26:   const DM_Stag *const stag = (DM_Stag *)dm->data;
     27:   PetscInt             dim;

     29:   PetscFunctionBegin;
     31:   PetscCall(DMGetDimension(dm, &dim));
     32:   if (boundaryTypeX) *boundaryTypeX = stag->boundaryType[0];
     33:   if (boundaryTypeY && dim > 1) *boundaryTypeY = stag->boundaryType[1];
     34:   if (boundaryTypeZ && dim > 2) *boundaryTypeZ = stag->boundaryType[2];
     35:   PetscFunctionReturn(PETSC_SUCCESS);
     36: }

     38: static PetscErrorCode DMStagGetProductCoordinateArrays_Private(DM dm, void *arrX, void *arrY, void *arrZ, PetscBool read)
     39: {
     40:   PetscInt  dim, d, dofCheck[DMSTAG_MAX_STRATA], s;
     41:   DM        dmCoord;
     42:   void     *arr[DMSTAG_MAX_DIM];
     43:   PetscBool checkDof;

     45:   PetscFunctionBegin;
     47:   PetscCall(DMGetDimension(dm, &dim));
     48:   PetscCheck(dim <= DMSTAG_MAX_DIM, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for %" PetscInt_FMT " dimensions", dim);
     49:   arr[0] = arrX;
     50:   arr[1] = arrY;
     51:   arr[2] = arrZ;
     52:   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
     53:   PetscCheck(dmCoord, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM does not have a coordinate DM");
     54:   {
     55:     PetscBool isProduct;
     56:     DMType    dmType;
     57:     PetscCall(DMGetType(dmCoord, &dmType));
     58:     PetscCall(PetscStrcmp(DMPRODUCT, dmType, &isProduct));
     59:     PetscCheck(isProduct, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Coordinate DM is not of type DMPRODUCT");
     60:   }
     61:   for (s = 0; s < DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
     62:   checkDof = PETSC_FALSE;
     63:   for (d = 0; d < dim; ++d) {
     64:     DM        subDM;
     65:     DMType    dmType;
     66:     PetscBool isStag;
     67:     PetscInt  dof[DMSTAG_MAX_STRATA], subDim;
     68:     Vec       coord1d_local;

     70:     /* Ignore unrequested arrays */
     71:     if (!arr[d]) continue;

     73:     PetscCall(DMProductGetDM(dmCoord, d, &subDM));
     74:     PetscCheck(subDM, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Coordinate DM is missing sub DM %" PetscInt_FMT, d);
     75:     PetscCall(DMGetDimension(subDM, &subDim));
     76:     PetscCheck(subDim == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Coordinate sub-DM is not of dimension 1");
     77:     PetscCall(DMGetType(subDM, &dmType));
     78:     PetscCall(PetscStrcmp(DMSTAG, dmType, &isStag));
     79:     PetscCheck(isStag, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Coordinate sub-DM is not of type DMSTAG");
     80:     PetscCall(DMStagGetDOF(subDM, &dof[0], &dof[1], &dof[2], &dof[3]));
     81:     if (!checkDof) {
     82:       for (s = 0; s < DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
     83:       checkDof = PETSC_TRUE;
     84:     } else {
     85:       for (s = 0; s < DMSTAG_MAX_STRATA; ++s) PetscCheck(dofCheck[s] == dof[s], PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Coordinate sub-DMs have different dofs");
     86:     }
     87:     PetscCall(DMGetCoordinatesLocal(subDM, &coord1d_local));
     88:     if (read) {
     89:       PetscCall(DMStagVecGetArrayRead(subDM, coord1d_local, arr[d]));
     90:     } else {
     91:       PetscCall(DMStagVecGetArray(subDM, coord1d_local, arr[d]));
     92:     }
     93:   }
     94:   PetscFunctionReturn(PETSC_SUCCESS);
     95: }

     97: /*@C
     98:   DMStagGetProductCoordinateArrays - extract local product coordinate arrays, one per dimension

    100:   Logically Collective

    102:   Input Parameter:
    103: . dm - the `DMSTAG` object

    105:   Output Parameters:
    106: + arrX - local 1D coordinate arrays for x direction
    107: . arrY - local 1D coordinate arrays for y direction, not set for one dimensional problems
    108: - arrZ - local 1D coordinate arrays for z direction, not set for one and two dimensional problems

    110:   Level: intermediate

    112:   Notes:
    113:   A high-level helper function to quickly extract local coordinate arrays.

    115:   Note that 2-dimensional arrays are returned. See
    116:   `DMStagVecGetArray()`, which is called internally to produce these arrays
    117:   representing coordinates on elements and vertices (element boundaries)
    118:   for a 1-dimensional `DMSTAG` in each coordinate direction.

    120:   One should use `DMStagGetProductCoordinateLocationSlot()` to determine appropriate
    121:   indices for the second dimension in these returned arrays. This function
    122:   checks that the coordinate array is a suitable product of 1-dimensional
    123:   `DMSTAG` objects.

    125: .seealso: [](ch_stag), `DMSTAG`, `DMPRODUCT`, `DMStagGetProductCoordinateArraysRead()`, `DMStagSetUniformCoordinates()`, `DMStagSetUniformCoordinatesProduct()`, `DMStagGetProductCoordinateLocationSlot()`
    126: @*/
    127: PetscErrorCode DMStagGetProductCoordinateArrays(DM dm, void *arrX, void *arrY, void *arrZ)
    128: {
    129:   PetscFunctionBegin;
    130:   PetscCall(DMStagGetProductCoordinateArrays_Private(dm, arrX, arrY, arrZ, PETSC_FALSE));
    131:   PetscFunctionReturn(PETSC_SUCCESS);
    132: }

    134: /*@C
    135:   DMStagGetProductCoordinateArraysRead - extract product coordinate arrays, read-only

    137:   Logically Collective

    139:   Input Parameter:
    140: . dm - the `DMSTAG` object

    142:   Output Parameters:
    143: + arrX - local 1D coordinate arrays for x direction
    144: . arrY - local 1D coordinate arrays for y direction, not set for one dimensional problems
    145: - arrZ - local 1D coordinate arrays for z direction, not set for one and two dimensional problems

    147:   Level: intermediate

    149:   Note:
    150:   See `DMStagGetProductCoordinateArrays()` for more information.

    152: .seealso: [](ch_stag), `DMSTAG`, `DMPRODUCT`, `DMStagGetProductCoordinateArrays()`, `DMStagSetUniformCoordinates()`, `DMStagSetUniformCoordinatesProduct()`, `DMStagGetProductCoordinateLocationSlot()`
    153: @*/
    154: PetscErrorCode DMStagGetProductCoordinateArraysRead(DM dm, void *arrX, void *arrY, void *arrZ)
    155: {
    156:   PetscFunctionBegin;
    157:   PetscCall(DMStagGetProductCoordinateArrays_Private(dm, arrX, arrY, arrZ, PETSC_TRUE));
    158:   PetscFunctionReturn(PETSC_SUCCESS);
    159: }

    161: /*@C
    162:   DMStagGetProductCoordinateLocationSlot - get slot for use with local product coordinate arrays

    164:   Not Collective

    166:   Input Parameters:
    167: + dm  - the `DMSTAG` object
    168: - loc - the grid location

    170:   Output Parameter:
    171: . slot - the index to use in local arrays

    173:   Level: intermediate

    175:   Notes:
    176:   High-level helper function to get slot indices for 1D coordinate `DM`s,
    177:   for use with `DMStagGetProductCoordinateArrays()` and related functions.

    179:   For `loc`, one should use `DMSTAG_LEFT`, `DMSTAG_ELEMENT`, or `DMSTAG_RIGHT` for "previous", "center" and "next"
    180:   locations, respectively, in each dimension.
    181:   One can equivalently use `DMSTAG_DOWN` or `DMSTAG_BACK` in place of `DMSTAG_LEFT`,
    182:   and `DMSTAG_UP` or `DMSTACK_FRONT` in place of `DMSTAG_RIGHT`;

    184:   This function checks that the coordinates are actually set up so that using the
    185:   slots from any of the 1D coordinate sub-`DM`s are valid for all the 1D coordinate sub-`DM`s.

    187: .seealso: [](ch_stag), `DMSTAG`, `DMPRODUCT`, `DMStagGetProductCoordinateArrays()`, `DMStagGetProductCoordinateArraysRead()`, `DMStagSetUniformCoordinates()`
    188: @*/
    189: PETSC_EXTERN PetscErrorCode DMStagGetProductCoordinateLocationSlot(DM dm, DMStagStencilLocation loc, PetscInt *slot)
    190: {
    191:   DM       dmCoord;
    192:   PetscInt dim, dofCheck[DMSTAG_MAX_STRATA], s, d;

    194:   PetscFunctionBegin;
    196:   PetscCall(DMGetDimension(dm, &dim));
    197:   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
    198:   PetscCheck(dmCoord, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM does not have a coordinate DM");
    199:   {
    200:     PetscBool isProduct;
    201:     DMType    dmType;
    202:     PetscCall(DMGetType(dmCoord, &dmType));
    203:     PetscCall(PetscStrcmp(DMPRODUCT, dmType, &isProduct));
    204:     PetscCheck(isProduct, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Coordinate DM is not of type DMPRODUCT");
    205:   }
    206:   for (s = 0; s < DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
    207:   for (d = 0; d < dim; ++d) {
    208:     DM        subDM;
    209:     DMType    dmType;
    210:     PetscBool isStag;
    211:     PetscInt  dof[DMSTAG_MAX_STRATA], subDim;
    212:     PetscCall(DMProductGetDM(dmCoord, d, &subDM));
    213:     PetscCheck(subDM, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Coordinate DM is missing sub DM %" PetscInt_FMT, d);
    214:     PetscCall(DMGetDimension(subDM, &subDim));
    215:     PetscCheck(subDim == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Coordinate sub-DM is not of dimension 1");
    216:     PetscCall(DMGetType(subDM, &dmType));
    217:     PetscCall(PetscStrcmp(DMSTAG, dmType, &isStag));
    218:     PetscCheck(isStag, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Coordinate sub-DM is not of type DMSTAG");
    219:     PetscCall(DMStagGetDOF(subDM, &dof[0], &dof[1], &dof[2], &dof[3]));
    220:     if (d == 0) {
    221:       const PetscInt component = 0;
    222:       for (s = 0; s < DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
    223:       PetscCall(DMStagGetLocationSlot(subDM, loc, component, slot));
    224:     } else {
    225:       for (s = 0; s < DMSTAG_MAX_STRATA; ++s) PetscCheck(dofCheck[s] == dof[s], PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Coordinate sub-DMs have different dofs");
    226:     }
    227:   }
    228:   PetscFunctionReturn(PETSC_SUCCESS);
    229: }

    231: /*@
    232:   DMStagGetCorners - return global element indices of the local region (excluding ghost points)

    234:   Not Collective

    236:   Input Parameter:
    237: . dm - the `DMSTAG` object

    239:   Output Parameters:
    240: + x       - starting element index in first direction
    241: . y       - starting element index in second direction
    242: . z       - starting element index in third direction
    243: . m       - element width in first direction
    244: . n       - element width in second direction
    245: . p       - element width in third direction
    246: . nExtrax - number of extra partial elements in first direction
    247: . nExtray - number of extra partial elements in second direction
    248: - nExtraz - number of extra partial elements in third direction

    250:   Level: beginner

    252:   Notes:
    253:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to `NULL` in this case.

    255:   The number of extra partial elements is either 1 or 0.
    256:   The value is 1 on right, top, and front non-periodic domain ("physical") boundaries,
    257:   in the x, y, and z directions respectively, and otherwise 0.

    259: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetGhostCorners()`, `DMDAGetCorners()`
    260: @*/
    261: PetscErrorCode DMStagGetCorners(DM dm, PetscInt *x, PetscInt *y, PetscInt *z, PetscInt *m, PetscInt *n, PetscInt *p, PetscInt *nExtrax, PetscInt *nExtray, PetscInt *nExtraz)
    262: {
    263:   const DM_Stag *const stag = (DM_Stag *)dm->data;

    265:   PetscFunctionBegin;
    267:   if (x) *x = stag->start[0];
    268:   if (y) *y = stag->start[1];
    269:   if (z) *z = stag->start[2];
    270:   if (m) *m = stag->n[0];
    271:   if (n) *n = stag->n[1];
    272:   if (p) *p = stag->n[2];
    273:   if (nExtrax) *nExtrax = stag->boundaryType[0] != DM_BOUNDARY_PERIODIC && stag->lastRank[0] ? 1 : 0;
    274:   if (nExtray) *nExtray = stag->boundaryType[1] != DM_BOUNDARY_PERIODIC && stag->lastRank[1] ? 1 : 0;
    275:   if (nExtraz) *nExtraz = stag->boundaryType[2] != DM_BOUNDARY_PERIODIC && stag->lastRank[2] ? 1 : 0;
    276:   PetscFunctionReturn(PETSC_SUCCESS);
    277: }

    279: /*@
    280:   DMStagGetDOF - get number of DOF associated with each stratum of the grid

    282:   Not Collective

    284:   Input Parameter:
    285: . dm - the `DMSTAG` object

    287:   Output Parameters:
    288: + dof0 - the number of points per 0-cell (vertex/node)
    289: . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
    290: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
    291: - dof3 - the number of points per 3-cell (element in 3D)

    293:   Level: beginner

    295: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetCorners()`, `DMStagGetGhostCorners()`, `DMStagGetGlobalSizes()`, `DMStagGetStencilWidth()`, `DMStagGetBoundaryTypes()`, `DMStagGetLocationDOF()`, `DMDAGetDof()`
    296: @*/
    297: PetscErrorCode DMStagGetDOF(DM dm, PetscInt *dof0, PetscInt *dof1, PetscInt *dof2, PetscInt *dof3)
    298: {
    299:   const DM_Stag *const stag = (DM_Stag *)dm->data;

    301:   PetscFunctionBegin;
    303:   if (dof0) *dof0 = stag->dof[0];
    304:   if (dof1) *dof1 = stag->dof[1];
    305:   if (dof2) *dof2 = stag->dof[2];
    306:   if (dof3) *dof3 = stag->dof[3];
    307:   PetscFunctionReturn(PETSC_SUCCESS);
    308: }

    310: /*@
    311:   DMStagGetGhostCorners - return global element indices of the local region, including ghost points

    313:   Not Collective

    315:   Input Parameter:
    316: . dm - the `DMSTAG` object

    318:   Output Parameters:
    319: + x - the starting element index in the first direction
    320: . y - the starting element index in the second direction
    321: . z - the starting element index in the third direction
    322: . m - the element width in the first direction
    323: . n - the element width in the second direction
    324: - p - the element width in the third direction

    326:   Level: beginner

    328:   Note:
    329:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to `NULL` in this case.

    331: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetCorners()`, `DMDAGetGhostCorners()`
    332: @*/
    333: PetscErrorCode DMStagGetGhostCorners(DM dm, PetscInt *x, PetscInt *y, PetscInt *z, PetscInt *m, PetscInt *n, PetscInt *p)
    334: {
    335:   const DM_Stag *const stag = (DM_Stag *)dm->data;

    337:   PetscFunctionBegin;
    339:   if (x) *x = stag->startGhost[0];
    340:   if (y) *y = stag->startGhost[1];
    341:   if (z) *z = stag->startGhost[2];
    342:   if (m) *m = stag->nGhost[0];
    343:   if (n) *n = stag->nGhost[1];
    344:   if (p) *p = stag->nGhost[2];
    345:   PetscFunctionReturn(PETSC_SUCCESS);
    346: }

    348: /*@
    349:   DMStagGetGlobalSizes - get global element counts

    351:   Not Collective

    353:   Input Parameter:
    354: . dm - the `DMSTAG` object

    356:   Output Parameters:
    357: + M - global element counts in the x direction
    358: . N - global element counts in the y direction
    359: - P - global element counts in the z direction

    361:   Level: beginner

    363:   Note:
    364:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to `NULL` in this case.

    366: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetLocalSizes()`, `DMDAGetInfo()`
    367: @*/
    368: PetscErrorCode DMStagGetGlobalSizes(DM dm, PetscInt *M, PetscInt *N, PetscInt *P)
    369: {
    370:   const DM_Stag *const stag = (DM_Stag *)dm->data;

    372:   PetscFunctionBegin;
    374:   if (M) *M = stag->N[0];
    375:   if (N) *N = stag->N[1];
    376:   if (P) *P = stag->N[2];
    377:   PetscFunctionReturn(PETSC_SUCCESS);
    378: }

    380: /*@C
    381:   DMStagGetIsFirstRank - get boolean value for whether this rank is first in each direction in the rank grid

    383:   Not Collective

    385:   Input Parameter:
    386: . dm - the `DMSTAG` object

    388:   Output Parameters:
    389: + isFirstRank0 - whether this rank is first in the x direction
    390: . isFirstRank1 - whether this rank is first in the y direction
    391: - isFirstRank2 - whether this rank is first in the z direction

    393:   Level: intermediate

    395:   Note:
    396:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to `NULL` in this case.

    398: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetIsLastRank()`
    399: @*/
    400: PetscErrorCode DMStagGetIsFirstRank(DM dm, PetscBool *isFirstRank0, PetscBool *isFirstRank1, PetscBool *isFirstRank2)
    401: {
    402:   const DM_Stag *const stag = (DM_Stag *)dm->data;

    404:   PetscFunctionBegin;
    406:   if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
    407:   if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
    408:   if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
    409:   PetscFunctionReturn(PETSC_SUCCESS);
    410: }

    412: /*@C
    413:   DMStagGetIsLastRank - get boolean value for whether this rank is last in each direction in the rank grid

    415:   Not Collective

    417:   Input Parameter:
    418: . dm - the `DMSTAG` object

    420:   Output Parameters:
    421: + isLastRank0 - whether this rank is last in the x direction
    422: . isLastRank1 - whether this rank is last in the y direction
    423: - isLastRank2 - whether this rank is last in the z direction

    425:   Level: intermediate

    427:   Note:
    428:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to `NULL` in this case.

    430: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetIsFirstRank()`
    431: @*/
    432: PetscErrorCode DMStagGetIsLastRank(DM dm, PetscBool *isLastRank0, PetscBool *isLastRank1, PetscBool *isLastRank2)
    433: {
    434:   const DM_Stag *const stag = (DM_Stag *)dm->data;

    436:   PetscFunctionBegin;
    438:   if (isLastRank0) *isLastRank0 = stag->lastRank[0];
    439:   if (isLastRank1) *isLastRank1 = stag->lastRank[1];
    440:   if (isLastRank2) *isLastRank2 = stag->lastRank[2];
    441:   PetscFunctionReturn(PETSC_SUCCESS);
    442: }

    444: /*@
    445:   DMStagGetLocalSizes - get local elementwise sizes

    447:   Not Collective

    449:   Input Parameter:
    450: . dm - the `DMSTAG` object

    452:   Output Parameters:
    453: + m - local element counts (excluding ghosts) in the x direction
    454: . n - local element counts (excluding ghosts) in the y direction
    455: - p - local element counts (excluding ghosts) in the z direction

    457:   Level: beginner

    459:   Note:
    460:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to `NULL` in this case.

    462: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetGlobalSizes()`, `DMStagGetDOF()`, `DMStagGetNumRanks()`, `DMDAGetLocalInfo()`
    463: @*/
    464: PetscErrorCode DMStagGetLocalSizes(DM dm, PetscInt *m, PetscInt *n, PetscInt *p)
    465: {
    466:   const DM_Stag *const stag = (DM_Stag *)dm->data;

    468:   PetscFunctionBegin;
    470:   if (m) *m = stag->n[0];
    471:   if (n) *n = stag->n[1];
    472:   if (p) *p = stag->n[2];
    473:   PetscFunctionReturn(PETSC_SUCCESS);
    474: }

    476: /*@
    477:   DMStagGetNumRanks - get number of ranks in each direction in the global grid decomposition

    479:   Not Collective

    481:   Input Parameter:
    482: . dm - the `DMSTAG` object

    484:   Output Parameters:
    485: + nRanks0 - number of ranks in the x direction in the grid decomposition
    486: . nRanks1 - number of ranks in the y direction in the grid decomposition
    487: - nRanks2 - number of ranks in the z direction in the grid decomposition

    489:   Level: intermediate

    491: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetGlobalSizes()`, `DMStagGetLocalSize()`, `DMStagSetNumRanks()`, `DMDAGetInfo()`
    492: @*/
    493: PetscErrorCode DMStagGetNumRanks(DM dm, PetscInt *nRanks0, PetscInt *nRanks1, PetscInt *nRanks2)
    494: {
    495:   const DM_Stag *const stag = (DM_Stag *)dm->data;

    497:   PetscFunctionBegin;
    499:   if (nRanks0) *nRanks0 = stag->nRanks[0];
    500:   if (nRanks1) *nRanks1 = stag->nRanks[1];
    501:   if (nRanks2) *nRanks2 = stag->nRanks[2];
    502:   PetscFunctionReturn(PETSC_SUCCESS);
    503: }

    505: /*@
    506:   DMStagGetEntries - get number of native entries in the global representation

    508:   Not Collective

    510:   Input Parameter:
    511: . dm - the `DMSTAG` object

    513:   Output Parameter:
    514: . entries - number of rank-native entries in the global representation

    516:   Level: developer

    518:   Note:
    519:   This is the number of entries on this rank for a global vector associated with `dm`.
    520:   That is, it is value of `size` returned by `VecGetLocalSize(vec,&size)` when
    521:   `DMCreateGlobalVector(dm,&vec) is used to create a `Vec`. Users would typically
    522:   use these functions.

    524: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetDOF()`, `DMStagGetEntriesLocal()`, `DMStagGetEntriesPerElement()`, `DMCreateLocalVector()`
    525: @*/
    526: PetscErrorCode DMStagGetEntries(DM dm, PetscInt *entries)
    527: {
    528:   const DM_Stag *const stag = (DM_Stag *)dm->data;

    530:   PetscFunctionBegin;
    532:   if (entries) *entries = stag->entries;
    533:   PetscFunctionReturn(PETSC_SUCCESS);
    534: }

    536: /*@
    537:   DMStagGetEntriesLocal - get number of entries in the local representation

    539:   Not Collective

    541:   Input Parameter:
    542: . dm - the `DMSTAG` object

    544:   Output Parameter:
    545: . entries - number of entries in the local representation

    547:   Level: developer

    549:   Note:
    550:   This is the number of entries on this rank in the local representation.
    551:   That is, it is value of `size` returned by `VecGetSize(vec,&size)` or
    552:   `VecGetLocalSize(vec,&size)` when `DMCreateLocalVector(dm,&vec)` is used to
    553:   create a `Vec`. Users would typically use these functions.

    555: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetDOF()`, `DMStagGetEntries()`, `DMStagGetEntriesPerElement()`, `DMCreateLocalVector()`
    556: @*/
    557: PetscErrorCode DMStagGetEntriesLocal(DM dm, PetscInt *entries)
    558: {
    559:   const DM_Stag *const stag = (DM_Stag *)dm->data;

    561:   PetscFunctionBegin;
    563:   if (entries) *entries = stag->entriesGhost;
    564:   PetscFunctionReturn(PETSC_SUCCESS);
    565: }

    567: /*@
    568:   DMStagGetEntriesPerElement - get number of entries per element in the local representation

    570:   Not Collective

    572:   Input Parameter:
    573: . dm - the `DMSTAG` object

    575:   Output Parameter:
    576: . entriesPerElement - number of entries associated with each element in the local representation

    578:   Level: developer

    580:   Notes:
    581:   This is the natural block size for most local operations. In 1D it is equal to `dof0` $+$ `dof1`,
    582:   in 2D it is equal to `dof0` $+ 2$`dof1` $+$ `dof2`, and in 3D it is equal to `dof0` $+ 3$`dof1` $+ 3$`dof2` $+$ `dof3`

    584: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetDOF()`
    585: @*/
    586: PetscErrorCode DMStagGetEntriesPerElement(DM dm, PetscInt *entriesPerElement)
    587: {
    588:   const DM_Stag *const stag = (DM_Stag *)dm->data;

    590:   PetscFunctionBegin;
    592:   if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
    593:   PetscFunctionReturn(PETSC_SUCCESS);
    594: }

    596: /*@
    597:   DMStagGetStencilType - get elementwise ghost/halo stencil type

    599:   Not Collective

    601:   Input Parameter:
    602: . dm - the `DMSTAG` object

    604:   Output Parameter:
    605: . stencilType - the elementwise ghost stencil type: `DMSTAG_STENCIL_BOX`, `DMSTAG_STENCIL_STAR`, or `DMSTAG_STENCIL_NONE`

    607:   Level: beginner

    609: .seealso: [](ch_stag), `DMSTAG`, `DMStagSetStencilType()`, `DMStagGetStencilWidth`, `DMStagStencilType`
    610: @*/
    611: PetscErrorCode DMStagGetStencilType(DM dm, DMStagStencilType *stencilType)
    612: {
    613:   DM_Stag *const stag = (DM_Stag *)dm->data;

    615:   PetscFunctionBegin;
    617:   *stencilType = stag->stencilType;
    618:   PetscFunctionReturn(PETSC_SUCCESS);
    619: }

    621: /*@
    622:   DMStagGetStencilWidth - get elementwise stencil width

    624:   Not Collective

    626:   Input Parameter:
    627: . dm - the `DMSTAG` object

    629:   Output Parameter:
    630: . stencilWidth - stencil/halo/ghost width in elements

    632:   Level: beginner

    634: .seealso: [](ch_stag), `DMSTAG`, `DMStagSetStencilWidth()`, `DMStagGetStencilType()`, `DMDAGetStencilType()`
    635: @*/
    636: PetscErrorCode DMStagGetStencilWidth(DM dm, PetscInt *stencilWidth)
    637: {
    638:   const DM_Stag *const stag = (DM_Stag *)dm->data;

    640:   PetscFunctionBegin;
    642:   if (stencilWidth) *stencilWidth = stag->stencilWidth;
    643:   PetscFunctionReturn(PETSC_SUCCESS);
    644: }

    646: /*@C
    647:   DMStagGetOwnershipRanges - get elements per rank in each direction

    649:   Not Collective

    651:   Input Parameter:
    652: . dm - the `DMSTAG` object

    654:   Output Parameters:
    655: + lx - ownership along x direction (optional)
    656: . ly - ownership along y direction (optional)
    657: - lz - ownership along z direction (optional)

    659:   Level: intermediate

    661:   Notes:
    662:   These correspond to the optional final arguments passed to `DMStagCreate1d()`, `DMStagCreate2d()`, and `DMStagCreate3d()`.

    664:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to `NULL` in this case.

    666:   In C you should not free these arrays, nor change the values in them.
    667:   They will only have valid values while the `DMSTAG` they came from still exists (has not been destroyed).

    669: .seealso: [](ch_stag), `DMSTAG`, `DMStagSetGlobalSizes()`, `DMStagSetOwnershipRanges()`, `DMStagCreate1d()`, `DMStagCreate2d()`, `DMStagCreate3d()`, `DMDAGetOwnershipRanges()`
    670: @*/
    671: PetscErrorCode DMStagGetOwnershipRanges(DM dm, const PetscInt *lx[], const PetscInt *ly[], const PetscInt *lz[])
    672: {
    673:   const DM_Stag *const stag = (DM_Stag *)dm->data;

    675:   PetscFunctionBegin;
    677:   if (lx) *lx = stag->l[0];
    678:   if (ly) *ly = stag->l[1];
    679:   if (lz) *lz = stag->l[2];
    680:   PetscFunctionReturn(PETSC_SUCCESS);
    681: }

    683: /*@
    684:   DMStagCreateCompatibleDMStag - create a compatible `DMSTAG` with different dof/stratum

    686:   Collective

    688:   Input Parameters:
    689: + dm   - the `DMSTAG` object
    690: . dof0 - number of dof on the first stratum in the new `DMSTAG`
    691: . dof1 - number of dof on the second stratum in the new `DMSTAG`
    692: . dof2 - number of dof on the third stratum in the new `DMSTAG`
    693: - dof3 - number of dof on the fourth stratum in the new `DMSTAG`

    695:   Output Parameter:
    696: . newdm - the new, compatible `DMSTAG`

    698:   Level: intermediate

    700:   Notes:
    701:   DOF supplied for strata too big for the dimension are ignored; these may be set to `0`.
    702:   For example, for a 2-dimensional `DMSTAG`, `dof2` sets the number of dof per element,
    703:   and `dof3` is unused. For a 3-dimensional `DMSTAG`, `dof3` sets the number of DOF per element.

    705:   In contrast to `DMDACreateCompatibleDMDA()`, coordinates are not reused.

    707: .seealso: [](ch_stag), `DMSTAG`, `DMDACreateCompatibleDMDA()`, `DMGetCompatibility()`, `DMStagMigrateVec()`
    708: @*/
    709: PetscErrorCode DMStagCreateCompatibleDMStag(DM dm, PetscInt dof0, PetscInt dof1, PetscInt dof2, PetscInt dof3, DM *newdm)
    710: {
    711:   PetscFunctionBegin;
    713:   PetscCall(DMStagDuplicateWithoutSetup(dm, PetscObjectComm((PetscObject)dm), newdm));
    714:   PetscCall(DMStagSetDOF(*newdm, dof0, dof1, dof2, dof3));
    715:   PetscCall(DMSetUp(*newdm));
    716:   PetscFunctionReturn(PETSC_SUCCESS);
    717: }

    719: /*@
    720:   DMStagGetLocationSlot - get index to use in accessing raw local arrays

    722:   Not Collective

    724:   Input Parameters:
    725: + dm  - the `DMSTAG` object
    726: . loc - location relative to an element
    727: - c   - component

    729:   Output Parameter:
    730: . slot - index to use

    732:   Level: beginner

    734:   Notes:
    735:   Provides an appropriate index to use with `DMStagVecGetArray()` and friends.
    736:   This is required so that the user doesn't need to know about the ordering of
    737:   dof associated with each local element.

    739: .seealso: [](ch_stag), `DMSTAG`, `DMStagVecGetArray()`, `DMStagVecGetArrayRead()`, `DMStagGetDOF()`, `DMStagGetEntriesPerElement()`
    740: @*/
    741: PetscErrorCode DMStagGetLocationSlot(DM dm, DMStagStencilLocation loc, PetscInt c, PetscInt *slot)
    742: {
    743:   DM_Stag *const stag = (DM_Stag *)dm->data;

    745:   PetscFunctionBegin;
    747:   if (PetscDefined(USE_DEBUG)) {
    748:     PetscInt dof;
    749:     PetscCall(DMStagGetLocationDOF(dm, loc, &dof));
    750:     PetscCheck(dof >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Location %s has no dof attached", DMStagStencilLocations[loc]);
    751:     PetscCheck(c <= dof - 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Supplied component number (%" PetscInt_FMT ") for location %s is too big (maximum %" PetscInt_FMT ")", c, DMStagStencilLocations[loc], dof - 1);
    752:   }
    753:   *slot = stag->locationOffsets[loc] + c;
    754:   PetscFunctionReturn(PETSC_SUCCESS);
    755: }

    757: /*@
    758:   DMStagGetRefinementFactor - get refinement ratios in each direction

    760:   Not Collective

    762:   Input Parameter:
    763: . dm - the `DMSTAG` object

    765:   Output Parameters:
    766: + refine_x - ratio of fine grid to coarse in x-direction (2 by default)
    767: . refine_y - ratio of fine grid to coarse in y-direction (2 by default)
    768: - refine_z - ratio of fine grid to coarse in z-direction (2 by default)

    770:   Level: intermediate

    772: .seealso: [](ch_stag), `DMSTAG`, `DMRefine()`, `DMCoarsen()`, `DMStagSetRefinementFactor()`, `DMDASetRefinementFactor()`
    773: @*/
    774: PetscErrorCode DMStagGetRefinementFactor(DM dm, PetscInt *refine_x, PetscInt *refine_y, PetscInt *refine_z)
    775: {
    776:   DM_Stag *const stag = (DM_Stag *)dm->data;

    778:   PetscFunctionBegin;
    780:   if (refine_x) *refine_x = stag->refineFactor[0];
    781:   if (refine_y) *refine_y = stag->refineFactor[1];
    782:   if (refine_z) *refine_z = stag->refineFactor[2];
    783:   PetscFunctionReturn(PETSC_SUCCESS);
    784: }

    786: /*@
    787:   DMStagMigrateVec - transfer a vector associated with a `DMSTAG` to a vector associated with a compatible `DMSTAG`

    789:   Collective

    791:   Input Parameters:
    792: + dm    - the source `DMSTAG` object
    793: . vec   - the source vector, compatible with `dm`
    794: . dmTo  - the compatible destination `DMSTAG` object
    795: - vecTo - the destination vector, compatible with `dmTo`

    797:   Level: advanced

    799:   Notes:
    800:   Extra dof are ignored, and unfilled dof are zeroed.
    801:   Currently only implemented to migrate global vectors to global vectors.
    802:   For the definition of compatibility of `DM`s, see `DMGetCompatibility()`.

    804: .seealso: [](ch_stag), `DMSTAG`, `DMStagCreateCompatibleDMStag()`, `DMGetCompatibility()`, `DMStagVecSplitToDMDA()`
    805: @*/
    806: PetscErrorCode DMStagMigrateVec(DM dm, Vec vec, DM dmTo, Vec vecTo)
    807: {
    808:   DM_Stag *const     stag   = (DM_Stag *)dm->data;
    809:   DM_Stag *const     stagTo = (DM_Stag *)dmTo->data;
    810:   PetscInt           nLocalTo, nLocal, dim, i, j, k;
    811:   PetscInt           start[DMSTAG_MAX_DIM], startGhost[DMSTAG_MAX_DIM], n[DMSTAG_MAX_DIM], nExtra[DMSTAG_MAX_DIM], offset[DMSTAG_MAX_DIM];
    812:   Vec                vecToLocal, vecLocal;
    813:   PetscBool          compatible, compatibleSet;
    814:   const PetscScalar *arr;
    815:   PetscScalar       *arrTo;
    816:   const PetscInt     epe   = stag->entriesPerElement;
    817:   const PetscInt     epeTo = stagTo->entriesPerElement;

    819:   PetscFunctionBegin;
    824:   PetscCall(DMGetCompatibility(dm, dmTo, &compatible, &compatibleSet));
    825:   PetscCheck(compatibleSet && compatible, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_INCOMP, "DMStag objects must be shown to be compatible");
    826:   PetscCall(DMGetDimension(dm, &dim));
    827:   PetscCall(VecGetLocalSize(vec, &nLocal));
    828:   PetscCall(VecGetLocalSize(vecTo, &nLocalTo));
    829:   PetscCheck(nLocal == stag->entries && nLocalTo == stagTo->entries, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Vector migration only implemented for global vector to global vector.");
    830:   PetscCall(DMStagGetCorners(dm, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], &nExtra[0], &nExtra[1], &nExtra[2]));
    831:   PetscCall(DMStagGetGhostCorners(dm, &startGhost[0], &startGhost[1], &startGhost[2], NULL, NULL, NULL));
    832:   for (i = 0; i < DMSTAG_MAX_DIM; ++i) offset[i] = start[i] - startGhost[i];

    834:   /* Proceed by transferring to a local vector, copying, and transferring back to a global vector */
    835:   PetscCall(DMGetLocalVector(dm, &vecLocal));
    836:   PetscCall(DMGetLocalVector(dmTo, &vecToLocal));
    837:   PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal));
    838:   PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal));
    839:   PetscCall(VecGetArrayRead(vecLocal, &arr));
    840:   PetscCall(VecGetArray(vecToLocal, &arrTo));
    841:   /* Note that some superfluous copying of entries on partial dummy elements is done */
    842:   if (dim == 1) {
    843:     for (i = offset[0]; i < offset[0] + n[0] + nExtra[0]; ++i) {
    844:       PetscInt d = 0, dTo = 0, b = 0, bTo = 0;
    845:       PetscInt si;
    846:       for (si = 0; si < 2; ++si) {
    847:         b += stag->dof[si];
    848:         bTo += stagTo->dof[si];
    849:         for (; d < b && dTo < bTo; ++d, ++dTo) arrTo[i * epeTo + dTo] = arr[i * epe + d];
    850:         for (; dTo < bTo; ++dTo) arrTo[i * epeTo + dTo] = 0.0;
    851:         d = b;
    852:       }
    853:     }
    854:   } else if (dim == 2) {
    855:     const PetscInt epr   = stag->nGhost[0] * epe;
    856:     const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
    857:     for (j = offset[1]; j < offset[1] + n[1] + nExtra[1]; ++j) {
    858:       for (i = offset[0]; i < offset[0] + n[0] + nExtra[0]; ++i) {
    859:         const PetscInt base   = j * epr + i * epe;
    860:         const PetscInt baseTo = j * eprTo + i * epeTo;
    861:         PetscInt       d = 0, dTo = 0, b = 0, bTo = 0;
    862:         const PetscInt s[4] = {0, 1, 1, 2}; /* Dimensions of points, in order */
    863:         PetscInt       si;
    864:         for (si = 0; si < 4; ++si) {
    865:           b += stag->dof[s[si]];
    866:           bTo += stagTo->dof[s[si]];
    867:           for (; d < b && dTo < bTo; ++d, ++dTo) arrTo[baseTo + dTo] = arr[base + d];
    868:           for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
    869:           d = b;
    870:         }
    871:       }
    872:     }
    873:   } else if (dim == 3) {
    874:     const PetscInt epr   = stag->nGhost[0] * epe;
    875:     const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
    876:     const PetscInt epl   = stag->nGhost[1] * epr;
    877:     const PetscInt eplTo = stagTo->nGhost[1] * eprTo;
    878:     for (k = offset[2]; k < offset[2] + n[2] + nExtra[2]; ++k) {
    879:       for (j = offset[1]; j < offset[1] + n[1] + nExtra[1]; ++j) {
    880:         for (i = offset[0]; i < offset[0] + n[0] + nExtra[0]; ++i) {
    881:           PetscInt       d = 0, dTo = 0, b = 0, bTo = 0;
    882:           const PetscInt base   = k * epl + j * epr + i * epe;
    883:           const PetscInt baseTo = k * eplTo + j * eprTo + i * epeTo;
    884:           const PetscInt s[8]   = {0, 1, 1, 2, 1, 2, 2, 3}; /* dimensions of points, in order */
    885:           PetscInt       is;
    886:           for (is = 0; is < 8; ++is) {
    887:             b += stag->dof[s[is]];
    888:             bTo += stagTo->dof[s[is]];
    889:             for (; d < b && dTo < bTo; ++d, ++dTo) arrTo[baseTo + dTo] = arr[base + d];
    890:             for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
    891:             d = b;
    892:           }
    893:         }
    894:       }
    895:     }
    896:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
    897:   PetscCall(VecRestoreArrayRead(vecLocal, &arr));
    898:   PetscCall(VecRestoreArray(vecToLocal, &arrTo));
    899:   PetscCall(DMRestoreLocalVector(dm, &vecLocal));
    900:   PetscCall(DMLocalToGlobalBegin(dmTo, vecToLocal, INSERT_VALUES, vecTo));
    901:   PetscCall(DMLocalToGlobalEnd(dmTo, vecToLocal, INSERT_VALUES, vecTo));
    902:   PetscCall(DMRestoreLocalVector(dmTo, &vecToLocal));
    903:   PetscFunctionReturn(PETSC_SUCCESS);
    904: }

    906: /*@
    907:   DMStagPopulateLocalToGlobalInjective - populate an internal 1-to-1 local-to-global map

    909:   Collective

    911:   Creates an internal object which explicitly maps a single local degree of
    912:   freedom to each global degree of freedom. This is used, if populated,
    913:   instead of SCATTER_REVERSE_LOCAL with the (1-to-many, in general)
    914:   global-to-local map, when DMLocalToGlobal() is called with INSERT_VALUES.
    915:   This allows usage, for example, even in the periodic, 1-rank case, where
    916:   the inverse of the global-to-local map, even when restricted to on-rank
    917:   communication, is non-injective. This is at the cost of storing an additional
    918:   VecScatter object inside each `DMSTAG` object.

    920:   Input Parameter:
    921: . dm - the `DMSTAG` object

    923:   Level: developer

    925:   Notes:
    926:   In normal usage, library users shouldn't be concerned with this function,
    927:   as it is called during `DMSetUp()`, when required.

    929:   Returns immediately if the internal map is already populated.

    931:   Developer Notes:
    932:   This could, if desired, be moved up to a general `DM` routine. It would allow,
    933:   for example, `DMDA` to support `DMLocalToGlobal()` with `INSERT_VALUES`,
    934:   even in the single-rank periodic case.

    936: .seealso: [](ch_stag), `DMSTAG`, `DMLocalToGlobal()`, `VecScatter`
    937: @*/
    938: PetscErrorCode DMStagPopulateLocalToGlobalInjective(DM dm)
    939: {
    940:   PetscInt       dim;
    941:   DM_Stag *const stag = (DM_Stag *)dm->data;

    943:   PetscFunctionBegin;
    945:   if (stag->ltog_injective) PetscFunctionReturn(PETSC_SUCCESS); /* Don't re-populate */
    946:   PetscCall(DMGetDimension(dm, &dim));
    947:   switch (dim) {
    948:   case 1:
    949:     PetscCall(DMStagPopulateLocalToGlobalInjective_1d(dm));
    950:     break;
    951:   case 2:
    952:     PetscCall(DMStagPopulateLocalToGlobalInjective_2d(dm));
    953:     break;
    954:   case 3:
    955:     PetscCall(DMStagPopulateLocalToGlobalInjective_3d(dm));
    956:     break;
    957:   default:
    958:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
    959:   }
    960:   PetscFunctionReturn(PETSC_SUCCESS);
    961: }

    963: static PetscErrorCode DMStagRestoreProductCoordinateArrays_Private(DM dm, void *arrX, void *arrY, void *arrZ, PetscBool read)
    964: {
    965:   PetscInt dim, d;
    966:   void    *arr[DMSTAG_MAX_DIM];
    967:   DM       dmCoord;

    969:   PetscFunctionBegin;
    971:   PetscCall(DMGetDimension(dm, &dim));
    972:   PetscCheck(dim <= DMSTAG_MAX_DIM, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for %" PetscInt_FMT " dimensions", dim);
    973:   arr[0] = arrX;
    974:   arr[1] = arrY;
    975:   arr[2] = arrZ;
    976:   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
    977:   for (d = 0; d < dim; ++d) {
    978:     DM  subDM;
    979:     Vec coord1d_local;

    981:     /* Ignore unrequested arrays */
    982:     if (!arr[d]) continue;

    984:     PetscCall(DMProductGetDM(dmCoord, d, &subDM));
    985:     PetscCall(DMGetCoordinatesLocal(subDM, &coord1d_local));
    986:     if (read) {
    987:       PetscCall(DMStagVecRestoreArrayRead(subDM, coord1d_local, arr[d]));
    988:     } else {
    989:       PetscCall(DMStagVecRestoreArray(subDM, coord1d_local, arr[d]));
    990:     }
    991:   }
    992:   PetscFunctionReturn(PETSC_SUCCESS);
    993: }

    995: /*@C
    996:   DMStagRestoreProductCoordinateArrays - restore local array access

    998:   Logically Collective

    1000:   Input Parameter:
    1001: . dm - the `DMSTAG` object

    1003:   Output Parameters:
    1004: + arrX - local 1D coordinate arrays for x direction
    1005: . arrY - local 1D coordinate arrays for y direction
    1006: - arrZ - local 1D coordinate arrays for z direction

    1008:   Level: intermediate

    1010:   Notes:
    1011:   This function does not automatically perform a local->global scatter to populate global coordinates from the local coordinates.
    1012:   Thus, it may be required to explicitly perform these operations in some situations, as in the following partial example\:
    1013: .vb
    1014:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    1015:   for (PetscInt d = 0; d < 3; ++d) {
    1016:     DM  subdm;
    1017:     Vec coor, coor_local;

    1019:     PetscCall(DMProductGetDM(cdm, d, &subdm));
    1020:     PetscCall(DMGetCoordinates(subdm, &coor));
    1021:     PetscCall(DMGetCoordinatesLocal(subdm, &coor_local));
    1022:     PetscCall(DMLocalToGlobal(subdm, coor_local, INSERT_VALUES, coor));
    1023:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Coordinates dim %" PetscInt_FMT ":\n", d));
    1024:     PetscCall(VecView(coor, PETSC_VIEWER_STDOUT_WORLD));
    1025:   }
    1026: .ve

    1028: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetProductCoordinateArrays()`, `DMStagGetProductCoordinateArraysRead()`
    1029: @*/
    1030: PetscErrorCode DMStagRestoreProductCoordinateArrays(DM dm, void *arrX, void *arrY, void *arrZ)
    1031: {
    1032:   PetscFunctionBegin;
    1033:   PetscCall(DMStagRestoreProductCoordinateArrays_Private(dm, arrX, arrY, arrZ, PETSC_FALSE));
    1034:   PetscFunctionReturn(PETSC_SUCCESS);
    1035: }

    1037: /*@C
    1038:   DMStagRestoreProductCoordinateArraysRead - restore local product array access, read-only

    1040:   Logically Collective

    1042:   Input Parameters:
    1043: + dm   - the `DMSTAG` object
    1044: . arrX - local 1D coordinate arrays for x direction
    1045: . arrY - local 1D coordinate arrays for y direction
    1046: - arrZ - local 1D coordinate arrays for z direction

    1048:   Level: intermediate

    1050: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetProductCoordinateArrays()`, `DMStagGetProductCoordinateArraysRead()`
    1051: @*/
    1052: PetscErrorCode DMStagRestoreProductCoordinateArraysRead(DM dm, void *arrX, void *arrY, void *arrZ)
    1053: {
    1054:   PetscFunctionBegin;
    1055:   PetscCall(DMStagRestoreProductCoordinateArrays_Private(dm, arrX, arrY, arrZ, PETSC_TRUE));
    1056:   PetscFunctionReturn(PETSC_SUCCESS);
    1057: }

    1059: /*@
    1060:   DMStagSetBoundaryTypes - set `DMSTAG` boundary types

    1062:   Logically Collective; boundaryType0, boundaryType1, and boundaryType2 must contain common values

    1064:   Input Parameters:
    1065: + dm            - the `DMSTAG` object
    1066: . boundaryType2 - boundary type for x direction
    1067: . boundaryType1 - boundary type for y direction, not set for one dimensional problems
    1068: - boundaryType0 - boundary type for z direction, not set for one and two dimensional problems

    1070:   Level: advanced

    1072:   Note:
    1073:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

    1075: .seealso: [](ch_stag), `DMSTAG`, `DMBoundaryType`, `DMStagCreate1d()`, `DMStagCreate2d()`, `DMStagCreate3d()`, `DMDASetBoundaryType()`
    1076: @*/
    1077: PetscErrorCode DMStagSetBoundaryTypes(DM dm, DMBoundaryType boundaryType0, DMBoundaryType boundaryType1, DMBoundaryType boundaryType2)
    1078: {
    1079:   DM_Stag *const stag = (DM_Stag *)dm->data;
    1080:   PetscInt       dim;

    1082:   PetscFunctionBegin;
    1083:   PetscCall(DMGetDimension(dm, &dim));
    1088:   PetscCheck(!dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
    1089:   stag->boundaryType[0] = boundaryType0;
    1090:   if (dim > 1) stag->boundaryType[1] = boundaryType1;
    1091:   if (dim > 2) stag->boundaryType[2] = boundaryType2;
    1092:   PetscFunctionReturn(PETSC_SUCCESS);
    1093: }

    1095: /*@C
    1096:   DMStagSetCoordinateDMType - set DM type to store coordinates

    1098:   Logically Collective; `dmtype` must contain common value

    1100:   Input Parameters:
    1101: + dm     - the `DMSTAG` object
    1102: - dmtype - `DMtype` for coordinates, either `DMSTAG` or `DMPRODUCT`

    1104:   Level: advanced

    1106: .seealso: [](ch_stag), `DMSTAG`, `DMPRODUCT`, `DMGetCoordinateDM()`, `DMStagSetUniformCoordinates()`, `DMStagSetUniformCoordinatesExplicit()`, `DMStagSetUniformCoordinatesProduct()`, `DMType`
    1107: @*/
    1108: PetscErrorCode DMStagSetCoordinateDMType(DM dm, DMType dmtype)
    1109: {
    1110:   DM_Stag *const stag = (DM_Stag *)dm->data;

    1112:   PetscFunctionBegin;
    1114:   PetscCall(PetscFree(stag->coordinateDMType));
    1115:   PetscCall(PetscStrallocpy(dmtype, (char **)&stag->coordinateDMType));
    1116:   PetscFunctionReturn(PETSC_SUCCESS);
    1117: }

    1119: /*@
    1120:   DMStagSetDOF - set dof/stratum

    1122:   Logically Collective; `dof0`, `dof1`, `dof2`, and `dof3` must contain common values

    1124:   Input Parameters:
    1125: + dm   - the `DMSTAG` object
    1126: . dof0 - the number of points per 0-cell (vertex/node)
    1127: . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
    1128: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
    1129: - dof3 - the number of points per 3-cell (element in 3D)

    1131:   Level: advanced

    1133:   Note:
    1134:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

    1136: .seealso: [](ch_stag), `DMSTAG`, `DMDASetDof()`
    1137: @*/
    1138: PetscErrorCode DMStagSetDOF(DM dm, PetscInt dof0, PetscInt dof1, PetscInt dof2, PetscInt dof3)
    1139: {
    1140:   DM_Stag *const stag = (DM_Stag *)dm->data;
    1141:   PetscInt       dim;

    1143:   PetscFunctionBegin;
    1149:   PetscCheck(!dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
    1150:   PetscCall(DMGetDimension(dm, &dim));
    1151:   PetscCheck(dof0 >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "dof0 cannot be negative");
    1152:   PetscCheck(dof1 >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "dof1 cannot be negative");
    1153:   PetscCheck(dim <= 1 || dof2 >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "dof2 cannot be negative");
    1154:   PetscCheck(dim <= 2 || dof3 >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "dof3 cannot be negative");
    1155:   stag->dof[0] = dof0;
    1156:   stag->dof[1] = dof1;
    1157:   if (dim > 1) stag->dof[2] = dof2;
    1158:   if (dim > 2) stag->dof[3] = dof3;
    1159:   PetscFunctionReturn(PETSC_SUCCESS);
    1160: }

    1162: /*@
    1163:   DMStagSetNumRanks - set ranks in each direction in the global rank grid

    1165:   Logically Collective; `nRanks0`, `nRanks1`, and `nRanks2` must contain common values

    1167:   Input Parameters:
    1168: + dm      - the `DMSTAG` object
    1169: . nRanks0 - number of ranks in the x direction
    1170: . nRanks1 - number of ranks in the y direction
    1171: - nRanks2 - number of ranks in the z direction

    1173:   Level: developer

    1175:   Note:
    1176:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

    1178: .seealso: [](ch_stag), `DMSTAG`, `DMDASetNumProcs()`
    1179: @*/
    1180: PetscErrorCode DMStagSetNumRanks(DM dm, PetscInt nRanks0, PetscInt nRanks1, PetscInt nRanks2)
    1181: {
    1182:   DM_Stag *const stag = (DM_Stag *)dm->data;
    1183:   PetscInt       dim;

    1185:   PetscFunctionBegin;
    1190:   PetscCheck(!dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
    1191:   PetscCall(DMGetDimension(dm, &dim));
    1192:   PetscCheck(nRanks0 == PETSC_DECIDE || nRanks0 >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "number of ranks in X direction cannot be less than 1");
    1193:   PetscCheck(dim <= 1 || nRanks1 == PETSC_DECIDE || nRanks1 >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "number of ranks in Y direction cannot be less than 1");
    1194:   PetscCheck(dim <= 2 || nRanks2 == PETSC_DECIDE || nRanks2 >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "number of ranks in Z direction cannot be less than 1");
    1195:   if (nRanks0) PetscCall(PetscMPIIntCast(nRanks0, &stag->nRanks[0]));
    1196:   if (dim > 1 && nRanks1) PetscCall(PetscMPIIntCast(nRanks1, &stag->nRanks[1]));
    1197:   if (dim > 2 && nRanks2) PetscCall(PetscMPIIntCast(nRanks2, &stag->nRanks[2]));
    1198:   PetscFunctionReturn(PETSC_SUCCESS);
    1199: }

    1201: /*@
    1202:   DMStagSetStencilType - set elementwise ghost/halo stencil type

    1204:   Logically Collective; `stencilType` must contain common value

    1206:   Input Parameters:
    1207: + dm          - the `DMSTAG` object
    1208: - stencilType - the elementwise ghost stencil type: `DMSTAG_STENCIL_BOX`, `DMSTAG_STENCIL_STAR`, or `DMSTAG_STENCIL_NONE`

    1210:   Level: beginner

    1212: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetStencilType()`, `DMStagSetStencilWidth()`, `DMStagStencilType`
    1213: @*/
    1214: PetscErrorCode DMStagSetStencilType(DM dm, DMStagStencilType stencilType)
    1215: {
    1216:   DM_Stag *const stag = (DM_Stag *)dm->data;

    1218:   PetscFunctionBegin;
    1221:   PetscCheck(!dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
    1222:   stag->stencilType = stencilType;
    1223:   PetscFunctionReturn(PETSC_SUCCESS);
    1224: }

    1226: /*@
    1227:   DMStagSetStencilWidth - set elementwise stencil width

    1229:   Logically Collective; `stencilWidth` must contain common value

    1231:   Input Parameters:
    1232: + dm           - the `DMSTAG` object
    1233: - stencilWidth - stencil/halo/ghost width in elements

    1235:   Level: beginner

    1237:   Note:
    1238:   The width value is not used when `DMSTAG_STENCIL_NONE` is specified.

    1240: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetStencilWidth()`, `DMStagGetStencilType()`, `DMStagStencilType`
    1241: @*/
    1242: PetscErrorCode DMStagSetStencilWidth(DM dm, PetscInt stencilWidth)
    1243: {
    1244:   DM_Stag *const stag = (DM_Stag *)dm->data;

    1246:   PetscFunctionBegin;
    1249:   PetscCheck(!dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
    1250:   PetscCheck(stencilWidth >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Stencil width must be non-negative");
    1251:   stag->stencilWidth = stencilWidth;
    1252:   PetscFunctionReturn(PETSC_SUCCESS);
    1253: }

    1255: /*@
    1256:   DMStagSetGlobalSizes - set global element counts in each direction

    1258:   Logically Collective; `N0`, `N1`, and `N2` must contain common values

    1260:   Input Parameters:
    1261: + dm - the `DMSTAG` object
    1262: . N0 - global elementwise size in the x direction
    1263: . N1 - global elementwise size in the y direction
    1264: - N2 - global elementwise size in the z direction

    1266:   Level: advanced

    1268:   Note:
    1269:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

    1271: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetGlobalSizes()`, `DMDASetSizes()`
    1272: @*/
    1273: PetscErrorCode DMStagSetGlobalSizes(DM dm, PetscInt N0, PetscInt N1, PetscInt N2)
    1274: {
    1275:   DM_Stag *const stag = (DM_Stag *)dm->data;
    1276:   PetscInt       dim;

    1278:   PetscFunctionBegin;
    1280:   PetscCheck(!dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
    1281:   PetscCall(DMGetDimension(dm, &dim));
    1282:   PetscCheck(N0 >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Number of elements in X direction must be positive");
    1283:   PetscCheck(dim <= 1 || N1 >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Number of elements in Y direction must be positive");
    1284:   PetscCheck(dim <= 2 || N2 >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Number of elements in Z direction must be positive");
    1285:   if (N0) stag->N[0] = N0;
    1286:   if (N1) stag->N[1] = N1;
    1287:   if (N2) stag->N[2] = N2;
    1288:   PetscFunctionReturn(PETSC_SUCCESS);
    1289: }

    1291: /*@
    1292:   DMStagSetOwnershipRanges - set elements per rank in each direction

    1294:   Logically Collective; `lx`, `ly`, and `lz` must contain common values

    1296:   Input Parameters:
    1297: + dm - the `DMSTAG` object
    1298: . lx - element counts for each rank in the x direction, may be `NULL`
    1299: . ly - element counts for each rank in the y direction, may be `NULL`
    1300: - lz - element counts for each rank in the z direction, may be `NULL`

    1302:   Level: developer

    1304:   Note:
    1305:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to `NULL` in this case.

    1307: .seealso: [](ch_stag), `DMSTAG`, `DMStagSetGlobalSizes()`, `DMStagGetOwnershipRanges()`, `DMDASetOwnershipRanges()`
    1308: @*/
    1309: PetscErrorCode DMStagSetOwnershipRanges(DM dm, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
    1310: {
    1311:   DM_Stag *const  stag = (DM_Stag *)dm->data;
    1312:   const PetscInt *lin[3];
    1313:   PetscInt        d, dim;

    1315:   PetscFunctionBegin;
    1317:   PetscCheck(!dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
    1318:   lin[0] = lx;
    1319:   lin[1] = ly;
    1320:   lin[2] = lz;
    1321:   PetscCall(DMGetDimension(dm, &dim));
    1322:   for (d = 0; d < dim; ++d) {
    1323:     if (lin[d]) {
    1324:       PetscCheck(stag->nRanks[d] >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Cannot set ownership ranges before setting number of ranks");
    1325:       if (!stag->l[d]) PetscCall(PetscMalloc1(stag->nRanks[d], &stag->l[d]));
    1326:       PetscCall(PetscArraycpy(stag->l[d], lin[d], stag->nRanks[d]));
    1327:     }
    1328:   }
    1329:   PetscFunctionReturn(PETSC_SUCCESS);
    1330: }

    1332: /*@
    1333:   DMStagSetRefinementFactor - set refinement ratios in each direction

    1335:   Logically Collective

    1337:   Input Parameters:
    1338: + dm       - the `DMSTAG` object
    1339: . refine_x - ratio of fine grid to coarse in x-direction (2 by default)
    1340: . refine_y - ratio of fine grid to coarse in y-direction (2 by default)
    1341: - refine_z - ratio of fine grid to coarse in z-direction (2 by default)

    1343:   Level: intermediate

    1345:   Note:
    1346:   Pass `PETSC_IGNORE` to leave a value unchanged

    1348: .seealso: [](ch_stag), `DMSTAG`, `DMRefine()`, `DMCoarsen()`, `DMStagGetRefinementFactor()`, `DMDAGetRefinementFactor()`
    1349: @*/
    1350: PetscErrorCode DMStagSetRefinementFactor(DM dm, PetscInt refine_x, PetscInt refine_y, PetscInt refine_z)
    1351: {
    1352:   DM_Stag *const stag = (DM_Stag *)dm->data;

    1354:   PetscFunctionBegin;
    1356:   if (refine_x > 0) stag->refineFactor[0] = refine_x;
    1357:   if (refine_y > 0) stag->refineFactor[1] = refine_y;
    1358:   if (refine_z > 0) stag->refineFactor[2] = refine_z;
    1359:   PetscFunctionReturn(PETSC_SUCCESS);
    1360: }

    1362: /*@
    1363:   DMStagSetUniformCoordinates - set `DMSTAG` coordinates to be a uniform grid

    1365:   Collective

    1367:   Input Parameters:
    1368: + dm   - the `DMSTAG` object
    1369: . xmin - minimum global coordinate value in the x direction
    1370: . xmax - maximum global coordinate values in the x direction
    1371: . ymin - minimum global coordinate value in the y direction
    1372: . ymax - maximum global coordinate value in the y direction
    1373: . zmin - minimum global coordinate value in the z direction
    1374: - zmax - maximum global coordinate value in the z direction

    1376:   Level: advanced

    1378:   Notes:
    1379:   `DMSTAG` supports 2 different types of coordinate `DM`: `DMSTAG` and `DMPRODUCT`.
    1380:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

    1382:   Local coordinates are populated (using `DMSetCoordinatesLocal()`), linearly
    1383:   extrapolated to ghost cells, including those outside the physical domain.
    1384:   This is also done in case of periodic boundaries, meaning that the same
    1385:   global point may have different coordinates in different local representations,
    1386:   which are equivalent assuming a periodicity implied by the arguments to this function,
    1387:   i.e. two points are equivalent if their difference is a multiple of $($`xmax` $-$ `xmin` $)$
    1388:   in the x direction, $($ `ymax` $-$ `ymin` $)$ in the y direction, and $($ `zmax` $-$ `zmin` $)$ in the z direction.

    1390: .seealso: [](ch_stag), `DMSTAG`, `DMPRODUCT`, `DMStagSetUniformCoordinatesExplicit()`, `DMStagSetUniformCoordinatesProduct()`, `DMStagSetCoordinateDMType()`, `DMGetCoordinateDM()`, `DMGetCoordinates()`, `DMDASetUniformCoordinates()`, `DMBoundaryType`
    1391: @*/
    1392: PetscErrorCode DMStagSetUniformCoordinates(DM dm, PetscReal xmin, PetscReal xmax, PetscReal ymin, PetscReal ymax, PetscReal zmin, PetscReal zmax)
    1393: {
    1394:   DM_Stag *const stag = (DM_Stag *)dm->data;
    1395:   PetscBool      flg_stag, flg_product;

    1397:   PetscFunctionBegin;
    1399:   PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called after DMSetUp()");
    1400:   PetscCheck(stag->coordinateDMType, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "You must first call DMStagSetCoordinateDMType()");
    1401:   PetscCall(PetscStrcmp(stag->coordinateDMType, DMSTAG, &flg_stag));
    1402:   PetscCall(PetscStrcmp(stag->coordinateDMType, DMPRODUCT, &flg_product));
    1403:   if (flg_stag) {
    1404:     PetscCall(DMStagSetUniformCoordinatesExplicit(dm, xmin, xmax, ymin, ymax, zmin, zmax));
    1405:   } else if (flg_product) {
    1406:     PetscCall(DMStagSetUniformCoordinatesProduct(dm, xmin, xmax, ymin, ymax, zmin, zmax));
    1407:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported DM Type %s", stag->coordinateDMType);
    1408:   PetscFunctionReturn(PETSC_SUCCESS);
    1409: }

    1411: /*@
    1412:   DMStagSetUniformCoordinatesExplicit - set `DMSTAG` coordinates to be a uniform grid, storing all values

    1414:   Collective

    1416:   Input Parameters:
    1417: + dm   - the `DMSTAG` object
    1418: . xmin - minimum global coordinate value in the x direction
    1419: . xmax - maximum global coordinate value in the x direction
    1420: . ymin - minimum global coordinate value in the y direction
    1421: . ymax - maximum global coordinate value in the y direction
    1422: . zmin - minimum global coordinate value in the z direction
    1423: - zmax - maximum global coordinate value in the z direction

    1425:   Level: beginner

    1427:   Notes:
    1428:   `DMSTAG` supports 2 different types of coordinate `DM`: either another `DMSTAG`, or a `DMPRODUCT`.
    1429:   If the grid is orthogonal, using `DMPRODUCT` should be more efficient.

    1431:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

    1433:   See the manual page for `DMStagSetUniformCoordinates()` for information on how
    1434:   coordinates for dummy cells outside the physical domain boundary are populated.

    1436: .seealso: [](ch_stag), `DMSTAG`, `DMStagSetUniformCoordinates()`, `DMStagSetUniformCoordinatesProduct()`, `DMStagSetCoordinateDMType()`
    1437: @*/
    1438: PetscErrorCode DMStagSetUniformCoordinatesExplicit(DM dm, PetscReal xmin, PetscReal xmax, PetscReal ymin, PetscReal ymax, PetscReal zmin, PetscReal zmax)
    1439: {
    1440:   DM_Stag *const stag = (DM_Stag *)dm->data;
    1441:   PetscInt       dim;
    1442:   PetscBool      flg;

    1444:   PetscFunctionBegin;
    1446:   PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called after DMSetUp()");
    1447:   PetscCall(PetscStrcmp(stag->coordinateDMType, DMSTAG, &flg));
    1448:   PetscCheck(!stag->coordinateDMType || flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Refusing to change an already-set DM coordinate type");
    1449:   PetscCall(DMStagSetCoordinateDMType(dm, DMSTAG));
    1450:   PetscCall(DMGetDimension(dm, &dim));
    1451:   switch (dim) {
    1452:   case 1:
    1453:     PetscCall(DMStagSetUniformCoordinatesExplicit_1d(dm, xmin, xmax));
    1454:     break;
    1455:   case 2:
    1456:     PetscCall(DMStagSetUniformCoordinatesExplicit_2d(dm, xmin, xmax, ymin, ymax));
    1457:     break;
    1458:   case 3:
    1459:     PetscCall(DMStagSetUniformCoordinatesExplicit_3d(dm, xmin, xmax, ymin, ymax, zmin, zmax));
    1460:     break;
    1461:   default:
    1462:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
    1463:   }
    1464:   PetscCall(DMCoarsenHookRemove(dm, DMRestrictHook_Coordinates, NULL, NULL));
    1465:   PetscFunctionReturn(PETSC_SUCCESS);
    1466: }

    1468: /*@
    1469:   DMStagSetUniformCoordinatesProduct - create uniform coordinates, as a product of 1D arrays

    1471:   Set the coordinate `DM` to be a `DMPRODUCT` of 1D `DMSTAG` objects, each of which have a coordinate `DM` (also a 1d `DMSTAG`) holding uniform coordinates.

    1473:   Collective

    1475:   Input Parameters:
    1476: + dm   - the `DMSTAG` object
    1477: . xmin - minimum global coordinate value in the x direction
    1478: . xmax - maximum global coordinate value in the x direction
    1479: . ymin - minimum global coordinate value in the y direction
    1480: . ymax - maximum global coordinate value in the y direction
    1481: . zmin - minimum global coordinate value in the z direction
    1482: - zmax - maximum global coordinate value in the z direction

    1484:   Level: intermediate

    1486:   Notes:
    1487:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

    1489:   The per-dimension 1-dimensional `DMSTAG` objects that comprise the product
    1490:   always have active 0-cells (vertices, element boundaries) and 1-cells
    1491:   (element centers).

    1493:   See the manual page for `DMStagSetUniformCoordinates()` for information on how
    1494:   coordinates for dummy cells outside the physical domain boundary are populated.

    1496: .seealso: [](ch_stag), `DMSTAG`, `DMPRODUCT`, `DMStagSetUniformCoordinates()`, `DMStagSetUniformCoordinatesExplicit()`, `DMStagSetCoordinateDMType()`
    1497: @*/
    1498: PetscErrorCode DMStagSetUniformCoordinatesProduct(DM dm, PetscReal xmin, PetscReal xmax, PetscReal ymin, PetscReal ymax, PetscReal zmin, PetscReal zmax)
    1499: {
    1500:   DM_Stag *const stag = (DM_Stag *)dm->data;
    1501:   DM             dmc;
    1502:   PetscInt       dim, d, dof0, dof1;
    1503:   PetscBool      flg;

    1505:   PetscFunctionBegin;
    1507:   PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called after DMSetUp()");
    1508:   PetscCall(PetscStrcmp(stag->coordinateDMType, DMPRODUCT, &flg));
    1509:   PetscCheck(!stag->coordinateDMType || flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Refusing to change an already-set DM coordinate type");
    1510:   PetscCall(DMStagSetCoordinateDMType(dm, DMPRODUCT));
    1511:   PetscCall(DMGetCoordinateDM(dm, &dmc));
    1512:   PetscCall(DMGetDimension(dm, &dim));

    1514:   /* Create 1D sub-DMs, living on subcommunicators.
    1515:      Always include both vertex and element dof, regardless of the active strata of the DMStag */
    1516:   dof0 = 1;
    1517:   dof1 = 1;

    1519:   for (d = 0; d < dim; ++d) {
    1520:     DM                subdm;
    1521:     MPI_Comm          subcomm;
    1522:     PetscMPIInt       color;
    1523:     const PetscMPIInt key = 0; /* let existing rank break ties */

    1525:     /* Choose colors based on position in the plane orthogonal to this dim, and split */
    1526:     switch (d) {
    1527:     case 0:
    1528:       color = (dim > 1 ? stag->rank[1] : 0) + (dim > 2 ? stag->nRanks[1] * stag->rank[2] : 0);
    1529:       break;
    1530:     case 1:
    1531:       color = stag->rank[0] + (dim > 2 ? stag->nRanks[0] * stag->rank[2] : 0);
    1532:       break;
    1533:     case 2:
    1534:       color = stag->rank[0] + stag->nRanks[0] * stag->rank[1];
    1535:       break;
    1536:     default:
    1537:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension index %" PetscInt_FMT, d);
    1538:     }
    1539:     PetscCallMPI(MPI_Comm_split(PetscObjectComm((PetscObject)dm), color, key, &subcomm));

    1541:     /* Create sub-DMs living on these new communicators (which are destroyed by DMProduct) */
    1542:     PetscCall(DMStagCreate1d(subcomm, stag->boundaryType[d], stag->N[d], dof0, dof1, stag->stencilType, stag->stencilWidth, stag->l[d], &subdm));
    1543:     /* Copy vector and matrix type information from parent DM */
    1544:     PetscCall(DMSetVecType(subdm, dm->vectype));
    1545:     PetscCall(DMSetMatType(subdm, dm->mattype));
    1546:     PetscCall(DMSetUp(subdm));
    1547:     switch (d) {
    1548:     case 0:
    1549:       PetscCall(DMStagSetUniformCoordinatesExplicit(subdm, xmin, xmax, 0, 0, 0, 0));
    1550:       break;
    1551:     case 1:
    1552:       PetscCall(DMStagSetUniformCoordinatesExplicit(subdm, ymin, ymax, 0, 0, 0, 0));
    1553:       break;
    1554:     case 2:
    1555:       PetscCall(DMStagSetUniformCoordinatesExplicit(subdm, zmin, zmax, 0, 0, 0, 0));
    1556:       break;
    1557:     default:
    1558:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension index %" PetscInt_FMT, d);
    1559:     }
    1560:     PetscCall(DMProductSetDM(dmc, d, subdm));
    1561:     PetscCall(DMProductSetDimensionIndex(dmc, d, 0));
    1562:     PetscCall(DMDestroy(&subdm));
    1563:     PetscCallMPI(MPI_Comm_free(&subcomm));
    1564:   }
    1565:   PetscCall(DMCoarsenHookRemove(dm, DMRestrictHook_Coordinates, NULL, NULL));
    1566:   PetscFunctionReturn(PETSC_SUCCESS);
    1567: }

    1569: /*@C
    1570:   DMStagVecGetArray - get access to local array

    1572:   Logically Collective

    1574:   Input Parameters:
    1575: + dm  - the `DMSTAG` object
    1576: - vec - the `Vec` object

    1578:   Output Parameter:
    1579: . array - the array

    1581:   Level: beginner

    1583:   Note:
    1584:   This function returns a (dim+1)-dimensional array for a dim-dimensional
    1585:   `DMSTAG`.

    1587:   The first 1-3 dimensions indicate an element in the global
    1588:   numbering, using the standard C ordering.

    1590:   The final dimension in this array corresponds to a degree
    1591:   of freedom with respect to this element, for example corresponding to
    1592:   the element or one of its neighboring faces, edges, or vertices.

    1594:   For example, for a 3D `DMSTAG`, indexing is `array[k][j][i][idx]`, where `k` is the
    1595:   index in the z-direction, `j` is the index in the y-direction, and `i` is the
    1596:   index in the x-direction.

    1598:   `idx` is obtained with `DMStagGetLocationSlot()`, since the correct offset
    1599:   into the $(d+1)$-dimensional C array for a $d$-dimensional `DMSTAG` depends on the grid size and the number
    1600:   of DOF stored at each location.

    1602:   `DMStagVecRestoreArray()` must be called, once finished with the array

    1604: .seealso: [](ch_stag), `DMSTAG`, `DMStagVecGetArrayRead()`, `DMStagGetLocationSlot()`, `DMGetLocalVector()`, `DMCreateLocalVector()`, `DMGetGlobalVector()`, `DMCreateGlobalVector()`, `DMDAVecGetArray()`, `DMDAVecGetArrayDOF()`
    1605: @*/
    1606: PetscErrorCode DMStagVecGetArray(DM dm, Vec vec, void *array)
    1607: {
    1608:   DM_Stag *const stag = (DM_Stag *)dm->data;
    1609:   PetscInt       dim;
    1610:   PetscInt       nLocal;

    1612:   PetscFunctionBegin;
    1615:   PetscCall(DMGetDimension(dm, &dim));
    1616:   PetscCall(VecGetLocalSize(vec, &nLocal));
    1617:   PetscCheck(nLocal == stag->entriesGhost, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Vector local size %" PetscInt_FMT " is not compatible with DMStag local size %" PetscInt_FMT, nLocal, stag->entriesGhost);
    1618:   switch (dim) {
    1619:   case 1:
    1620:     PetscCall(VecGetArray2d(vec, stag->nGhost[0], stag->entriesPerElement, stag->startGhost[0], 0, (PetscScalar ***)array));
    1621:     break;
    1622:   case 2:
    1623:     PetscCall(VecGetArray3d(vec, stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar ****)array));
    1624:     break;
    1625:   case 3:
    1626:     PetscCall(VecGetArray4d(vec, stag->nGhost[2], stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[2], stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar *****)array));
    1627:     break;
    1628:   default:
    1629:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
    1630:   }
    1631:   PetscFunctionReturn(PETSC_SUCCESS);
    1632: }

    1634: /*@C
    1635:   DMStagVecGetArrayRead - get read-only access to a local array

    1637:   Logically Collective

    1639:   See the man page for `DMStagVecGetArray()` for more information.

    1641:   Input Parameters:
    1642: + dm  - the `DMSTAG` object
    1643: - vec - the `Vec` object

    1645:   Output Parameter:
    1646: . array - the read-only array

    1648:   Level: beginner

    1650:   Note:
    1651:   `DMStagVecRestoreArrayRead()` must be called, once finished with the array

    1653: .seealso: [](ch_stag), `DMSTAG`, `DMStagGetLocationSlot()`, `DMGetLocalVector()`, `DMCreateLocalVector()`, `DMGetGlobalVector()`, `DMCreateGlobalVector()`, `DMDAVecGetArrayRead()`, `DMDAVecGetArrayDOFRead()`
    1654: @*/
    1655: PetscErrorCode DMStagVecGetArrayRead(DM dm, Vec vec, void *array)
    1656: {
    1657:   DM_Stag *const stag = (DM_Stag *)dm->data;
    1658:   PetscInt       dim;
    1659:   PetscInt       nLocal;

    1661:   PetscFunctionBegin;
    1664:   PetscCall(DMGetDimension(dm, &dim));
    1665:   PetscCall(VecGetLocalSize(vec, &nLocal));
    1666:   PetscCheck(nLocal == stag->entriesGhost, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Vector local size %" PetscInt_FMT " is not compatible with DMSTAG local size %" PetscInt_FMT, nLocal, stag->entriesGhost);
    1667:   switch (dim) {
    1668:   case 1:
    1669:     PetscCall(VecGetArray2dRead(vec, stag->nGhost[0], stag->entriesPerElement, stag->startGhost[0], 0, (PetscScalar ***)array));
    1670:     break;
    1671:   case 2:
    1672:     PetscCall(VecGetArray3dRead(vec, stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar ****)array));
    1673:     break;
    1674:   case 3:
    1675:     PetscCall(VecGetArray4dRead(vec, stag->nGhost[2], stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[2], stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar *****)array));
    1676:     break;
    1677:   default:
    1678:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
    1679:   }
    1680:   PetscFunctionReturn(PETSC_SUCCESS);
    1681: }

    1683: /*@C
    1684:   DMStagVecRestoreArray - restore access to a raw array

    1686:   Logically Collective

    1688:   Input Parameters:
    1689: + dm  - the `DMSTAG` object
    1690: - vec - the `Vec` object

    1692:   Output Parameter:
    1693: . array - the array

    1695:   Level: beginner

    1697: .seealso: [](ch_stag), `DMSTAG`, `DMStagVecGetArray()`, `DMDAVecRestoreArray()`, `DMDAVecRestoreArrayDOF()`
    1698: @*/
    1699: PetscErrorCode DMStagVecRestoreArray(DM dm, Vec vec, void *array)
    1700: {
    1701:   DM_Stag *const stag = (DM_Stag *)dm->data;
    1702:   PetscInt       dim;
    1703:   PetscInt       nLocal;

    1705:   PetscFunctionBegin;
    1708:   PetscCall(DMGetDimension(dm, &dim));
    1709:   PetscCall(VecGetLocalSize(vec, &nLocal));
    1710:   PetscCheck(nLocal == stag->entriesGhost, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Vector local size %" PetscInt_FMT " is not compatible with DMSTAG local size %" PetscInt_FMT, nLocal, stag->entriesGhost);
    1711:   switch (dim) {
    1712:   case 1:
    1713:     PetscCall(VecRestoreArray2d(vec, stag->nGhost[0], stag->entriesPerElement, stag->startGhost[0], 0, (PetscScalar ***)array));
    1714:     break;
    1715:   case 2:
    1716:     PetscCall(VecRestoreArray3d(vec, stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar ****)array));
    1717:     break;
    1718:   case 3:
    1719:     PetscCall(VecRestoreArray4d(vec, stag->nGhost[2], stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[2], stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar *****)array));
    1720:     break;
    1721:   default:
    1722:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
    1723:   }
    1724:   PetscFunctionReturn(PETSC_SUCCESS);
    1725: }

    1727: /*@C
    1728:   DMStagVecRestoreArrayRead - restore read-only access to a raw array

    1730:   Logically Collective

    1732:   Input Parameters:
    1733: + dm  - the `DMSTAG` object
    1734: - vec - the Vec object

    1736:   Output Parameter:
    1737: . array - the read-only array

    1739:   Level: beginner

    1741: .seealso: [](ch_stag), `DMSTAG`, `DMStagVecGetArrayRead()`, `DMDAVecRestoreArrayRead()`, `DMDAVecRestoreArrayDOFRead()`
    1742: @*/
    1743: PetscErrorCode DMStagVecRestoreArrayRead(DM dm, Vec vec, void *array)
    1744: {
    1745:   DM_Stag *const stag = (DM_Stag *)dm->data;
    1746:   PetscInt       dim;
    1747:   PetscInt       nLocal;

    1749:   PetscFunctionBegin;
    1752:   PetscCall(DMGetDimension(dm, &dim));
    1753:   PetscCall(VecGetLocalSize(vec, &nLocal));
    1754:   PetscCheck(nLocal == stag->entriesGhost, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Vector local size %" PetscInt_FMT " is not compatible with DMSTAG local size %" PetscInt_FMT, nLocal, stag->entriesGhost);
    1755:   switch (dim) {
    1756:   case 1:
    1757:     PetscCall(VecRestoreArray2dRead(vec, stag->nGhost[0], stag->entriesPerElement, stag->startGhost[0], 0, (PetscScalar ***)array));
    1758:     break;
    1759:   case 2:
    1760:     PetscCall(VecRestoreArray3dRead(vec, stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar ****)array));
    1761:     break;
    1762:   case 3:
    1763:     PetscCall(VecRestoreArray4dRead(vec, stag->nGhost[2], stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[2], stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar *****)array));
    1764:     break;
    1765:   default:
    1766:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
    1767:   }
    1768:   PetscFunctionReturn(PETSC_SUCCESS);
    1769: }