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:       for (PetscInt si = 0; si < 2; ++si) {
846:         b += stag->dof[si];
847:         bTo += stagTo->dof[si];
848:         for (; d < b && dTo < bTo; ++d, ++dTo) arrTo[i * epeTo + dTo] = arr[i * epe + d];
849:         for (; dTo < bTo; ++dTo) arrTo[i * epeTo + dTo] = 0.0;
850:         d = b;
851:       }
852:     }
853:   } else if (dim == 2) {
854:     const PetscInt epr   = stag->nGhost[0] * epe;
855:     const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
856:     for (j = offset[1]; j < offset[1] + n[1] + nExtra[1]; ++j) {
857:       for (i = offset[0]; i < offset[0] + n[0] + nExtra[0]; ++i) {
858:         const PetscInt base   = j * epr + i * epe;
859:         const PetscInt baseTo = j * eprTo + i * epeTo;
860:         PetscInt       d = 0, dTo = 0, b = 0, bTo = 0;
861:         const PetscInt s[4] = {0, 1, 1, 2}; /* Dimensions of points, in order */
862:         for (PetscInt si = 0; si < 4; ++si) {
863:           b += stag->dof[s[si]];
864:           bTo += stagTo->dof[s[si]];
865:           for (; d < b && dTo < bTo; ++d, ++dTo) arrTo[baseTo + dTo] = arr[base + d];
866:           for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
867:           d = b;
868:         }
869:       }
870:     }
871:   } else if (dim == 3) {
872:     const PetscInt epr   = stag->nGhost[0] * epe;
873:     const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
874:     const PetscInt epl   = stag->nGhost[1] * epr;
875:     const PetscInt eplTo = stagTo->nGhost[1] * eprTo;
876:     for (k = offset[2]; k < offset[2] + n[2] + nExtra[2]; ++k) {
877:       for (j = offset[1]; j < offset[1] + n[1] + nExtra[1]; ++j) {
878:         for (i = offset[0]; i < offset[0] + n[0] + nExtra[0]; ++i) {
879:           PetscInt       d = 0, dTo = 0, b = 0, bTo = 0;
880:           const PetscInt base   = k * epl + j * epr + i * epe;
881:           const PetscInt baseTo = k * eplTo + j * eprTo + i * epeTo;
882:           const PetscInt s[8]   = {0, 1, 1, 2, 1, 2, 2, 3}; /* dimensions of points, in order */
883:           for (PetscInt is = 0; is < 8; ++is) {
884:             b += stag->dof[s[is]];
885:             bTo += stagTo->dof[s[is]];
886:             for (; d < b && dTo < bTo; ++d, ++dTo) arrTo[baseTo + dTo] = arr[base + d];
887:             for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
888:             d = b;
889:           }
890:         }
891:       }
892:     }
893:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
894:   PetscCall(VecRestoreArrayRead(vecLocal, &arr));
895:   PetscCall(VecRestoreArray(vecToLocal, &arrTo));
896:   PetscCall(DMRestoreLocalVector(dm, &vecLocal));
897:   PetscCall(DMLocalToGlobalBegin(dmTo, vecToLocal, INSERT_VALUES, vecTo));
898:   PetscCall(DMLocalToGlobalEnd(dmTo, vecToLocal, INSERT_VALUES, vecTo));
899:   PetscCall(DMRestoreLocalVector(dmTo, &vecToLocal));
900:   PetscFunctionReturn(PETSC_SUCCESS);
901: }

903: /*@
904:   DMStagPopulateLocalToGlobalInjective - populate an internal 1-to-1 local-to-global map

906:   Collective

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

917:   Input Parameter:
918: . dm - the `DMSTAG` object

920:   Level: developer

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

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

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

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

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

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

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

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

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

992: /*@C
993:   DMStagRestoreProductCoordinateArrays - restore local array access

995:   Logically Collective

997:   Input Parameter:
998: . dm - the `DMSTAG` object

1000:   Output Parameters:
1001: + arrX - local 1D coordinate arrays for x direction
1002: . arrY - local 1D coordinate arrays for y direction
1003: - arrZ - local 1D coordinate arrays for z direction

1005:   Level: intermediate

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

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

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

1034: /*@C
1035:   DMStagRestoreProductCoordinateArraysRead - restore local product array access, read-only

1037:   Logically Collective

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

1045:   Level: intermediate

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

1056: /*@
1057:   DMStagSetBoundaryTypes - set `DMSTAG` boundary types

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

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

1067:   Level: advanced

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

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

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

1092: /*@C
1093:   DMStagSetCoordinateDMType - set DM type to store coordinates

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

1097:   Input Parameters:
1098: + dm     - the `DMSTAG` object
1099: - dmtype - `DMtype` for coordinates, either `DMSTAG` or `DMPRODUCT`

1101:   Level: advanced

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

1109:   PetscFunctionBegin;
1111:   PetscCall(PetscFree(stag->coordinateDMType));
1112:   PetscCall(PetscStrallocpy(dmtype, (char **)&stag->coordinateDMType));
1113:   PetscFunctionReturn(PETSC_SUCCESS);
1114: }

1116: /*@
1117:   DMStagSetDOF - set dof/stratum

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

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

1128:   Level: advanced

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

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

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

1159: /*@
1160:   DMStagSetNumRanks - set ranks in each direction in the global rank grid

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

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

1170:   Level: developer

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

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

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

1198: /*@
1199:   DMStagSetStencilType - set elementwise ghost/halo stencil type

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

1203:   Input Parameters:
1204: + dm          - the `DMSTAG` object
1205: - stencilType - the elementwise ghost stencil type: `DMSTAG_STENCIL_BOX`, `DMSTAG_STENCIL_STAR`, or `DMSTAG_STENCIL_NONE`

1207:   Level: beginner

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

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

1223: /*@
1224:   DMStagSetStencilWidth - set elementwise stencil width

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

1228:   Input Parameters:
1229: + dm           - the `DMSTAG` object
1230: - stencilWidth - stencil/halo/ghost width in elements

1232:   Level: beginner

1234:   Note:
1235:   The width value is not used when `DMSTAG_STENCIL_NONE` is specified.

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

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

1252: /*@
1253:   DMStagSetGlobalSizes - set global element counts in each direction

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

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

1263:   Level: advanced

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

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

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

1288: /*@
1289:   DMStagSetOwnershipRanges - set elements per rank in each direction

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

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

1299:   Level: developer

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

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

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

1329: /*@
1330:   DMStagSetRefinementFactor - set refinement ratios in each direction

1332:   Logically Collective

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

1340:   Level: intermediate

1342:   Note:
1343:   Pass `PETSC_IGNORE` to leave a value unchanged

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

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

1359: /*@
1360:   DMStagSetUniformCoordinates - set `DMSTAG` coordinates to be a uniform grid

1362:   Collective

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

1373:   Level: advanced

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

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

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

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

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

1411:   Collective

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

1422:   Level: beginner

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

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

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

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

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

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

1468:   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.

1470:   Collective

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

1481:   Level: intermediate

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

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

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

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

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

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

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

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

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

1566: /*@C
1567:   DMStagVecGetArray - get access to local array

1569:   Logically Collective

1571:   Input Parameters:
1572: + dm  - the `DMSTAG` object
1573: - vec - the `Vec` object

1575:   Output Parameter:
1576: . array - the array

1578:   Level: beginner

1580:   Note:
1581:   This function returns a (dim+1)-dimensional array for a dim-dimensional
1582:   `DMSTAG`.

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

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

1591:   For example, for a 3D `DMSTAG`, indexing is `array[k][j][i][slot]`, where `k` is the
1592:   index in the z-direction, `j` is the index in the y-direction, and `i` is the
1593:   index in the x-direction.

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

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

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

1609:   PetscFunctionBegin;
1612:   PetscCall(DMGetDimension(dm, &dim));
1613:   PetscCall(VecGetLocalSize(vec, &nLocal));
1614:   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);
1615:   switch (dim) {
1616:   case 1:
1617:     PetscCall(VecGetArray2d(vec, stag->nGhost[0], stag->entriesPerElement, stag->startGhost[0], 0, (PetscScalar ***)array));
1618:     break;
1619:   case 2:
1620:     PetscCall(VecGetArray3d(vec, stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar ****)array));
1621:     break;
1622:   case 3:
1623:     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));
1624:     break;
1625:   default:
1626:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
1627:   }
1628:   PetscFunctionReturn(PETSC_SUCCESS);
1629: }

1631: /*@C
1632:   DMStagVecGetArrayRead - get read-only access to a local array

1634:   Logically Collective

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

1638:   Input Parameters:
1639: + dm  - the `DMSTAG` object
1640: - vec - the `Vec` object

1642:   Output Parameter:
1643: . array - the read-only array

1645:   Level: beginner

1647:   Note:
1648:   `DMStagVecRestoreArrayRead()` must be called, once finished with the array

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

1658:   PetscFunctionBegin;
1661:   PetscCall(DMGetDimension(dm, &dim));
1662:   PetscCall(VecGetLocalSize(vec, &nLocal));
1663:   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);
1664:   switch (dim) {
1665:   case 1:
1666:     PetscCall(VecGetArray2dRead(vec, stag->nGhost[0], stag->entriesPerElement, stag->startGhost[0], 0, (PetscScalar ***)array));
1667:     break;
1668:   case 2:
1669:     PetscCall(VecGetArray3dRead(vec, stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar ****)array));
1670:     break;
1671:   case 3:
1672:     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));
1673:     break;
1674:   default:
1675:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
1676:   }
1677:   PetscFunctionReturn(PETSC_SUCCESS);
1678: }

1680: /*@C
1681:   DMStagVecRestoreArray - restore access to a raw array

1683:   Logically Collective

1685:   Input Parameters:
1686: + dm  - the `DMSTAG` object
1687: - vec - the `Vec` object

1689:   Output Parameter:
1690: . array - the array

1692:   Level: beginner

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

1702:   PetscFunctionBegin;
1705:   PetscCall(DMGetDimension(dm, &dim));
1706:   PetscCall(VecGetLocalSize(vec, &nLocal));
1707:   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);
1708:   switch (dim) {
1709:   case 1:
1710:     PetscCall(VecRestoreArray2d(vec, stag->nGhost[0], stag->entriesPerElement, stag->startGhost[0], 0, (PetscScalar ***)array));
1711:     break;
1712:   case 2:
1713:     PetscCall(VecRestoreArray3d(vec, stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar ****)array));
1714:     break;
1715:   case 3:
1716:     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));
1717:     break;
1718:   default:
1719:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
1720:   }
1721:   PetscFunctionReturn(PETSC_SUCCESS);
1722: }

1724: /*@C
1725:   DMStagVecRestoreArrayRead - restore read-only access to a raw array

1727:   Logically Collective

1729:   Input Parameters:
1730: + dm  - the `DMSTAG` object
1731: - vec - the Vec object

1733:   Output Parameter:
1734: . array - the read-only array

1736:   Level: beginner

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

1746:   PetscFunctionBegin;
1749:   PetscCall(DMGetDimension(dm, &dim));
1750:   PetscCall(VecGetLocalSize(vec, &nLocal));
1751:   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);
1752:   switch (dim) {
1753:   case 1:
1754:     PetscCall(VecRestoreArray2dRead(vec, stag->nGhost[0], stag->entriesPerElement, stag->startGhost[0], 0, (PetscScalar ***)array));
1755:     break;
1756:   case 2:
1757:     PetscCall(VecRestoreArray3dRead(vec, stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar ****)array));
1758:     break;
1759:   case 3:
1760:     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));
1761:     break;
1762:   default:
1763:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
1764:   }
1765:   PetscFunctionReturn(PETSC_SUCCESS);
1766: }