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: [](chapter_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;

 30:   DMGetDimension(dm, &dim);
 31:   if (boundaryTypeX) *boundaryTypeX = stag->boundaryType[0];
 32:   if (boundaryTypeY && dim > 1) *boundaryTypeY = stag->boundaryType[1];
 33:   if (boundaryTypeZ && dim > 2) *boundaryTypeZ = stag->boundaryType[2];
 34:   return 0;
 35: }

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

 45:   DMGetDimension(dm, &dim);
 47:   arr[0] = arrX;
 48:   arr[1] = arrY;
 49:   arr[2] = arrZ;
 50:   DMGetCoordinateDM(dm, &dmCoord);
 52:   {
 53:     PetscBool isProduct;
 54:     DMType    dmType;
 55:     DMGetType(dmCoord, &dmType);
 56:     PetscStrcmp(DMPRODUCT, dmType, &isProduct);
 58:   }
 59:   for (s = 0; s < DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
 60:   checkDof = PETSC_FALSE;
 61:   for (d = 0; d < dim; ++d) {
 62:     DM        subDM;
 63:     DMType    dmType;
 64:     PetscBool isStag;
 65:     PetscInt  dof[DMSTAG_MAX_STRATA], subDim;
 66:     Vec       coord1d_local;

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

 71:     DMProductGetDM(dmCoord, d, &subDM);
 73:     DMGetDimension(subDM, &subDim);
 75:     DMGetType(subDM, &dmType);
 76:     PetscStrcmp(DMSTAG, dmType, &isStag);
 78:     DMStagGetDOF(subDM, &dof[0], &dof[1], &dof[2], &dof[3]);
 79:     if (!checkDof) {
 80:       for (s = 0; s < DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
 81:       checkDof = PETSC_TRUE;
 82:     } else {
 84:     }
 85:     DMGetCoordinatesLocal(subDM, &coord1d_local);
 86:     if (read) {
 87:       DMStagVecGetArrayRead(subDM, coord1d_local, arr[d]);
 88:     } else {
 89:       DMStagVecGetArray(subDM, coord1d_local, arr[d]);
 90:     }
 91:   }
 92:   return 0;
 93: }

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

 98:   Logically Collective

100:   Input Parameter:
101: . dm - the `DMSTAG` object

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

108:   Level: intermediate

110:   Notes:
111:   A high-level helper function to quickly extract local coordinate arrays.

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

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

123: .seealso: [](chapter_stag), `DMSTAG`, `DMPRODUCT`, `DMStagGetProductCoordinateArraysRead()`, `DMStagSetUniformCoordinates()`, `DMStagSetUniformCoordinatesProduct()`, `DMStagGetProductCoordinateLocationSlot()`
124: @*/
125: PetscErrorCode DMStagGetProductCoordinateArrays(DM dm, void *arrX, void *arrY, void *arrZ)
126: {
127:   DMStagGetProductCoordinateArrays_Private(dm, arrX, arrY, arrZ, PETSC_FALSE);
128:   return 0;
129: }

131: /*@C
132:   DMStagGetProductCoordinateArraysRead - extract product coordinate arrays, read-only

134:   Logically Collective

136:   Input Parameter:
137: . dm - the `DMSTAG` object

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

144:   Level: intermediate

146:   Note:
147:   See `DMStagGetProductCoordinateArrays()` for more information.

149: .seealso: [](chapter_stag), `DMSTAG`, `DMPRODUCT`, `DMStagGetProductCoordinateArrays()`, `DMStagSetUniformCoordinates()`, `DMStagSetUniformCoordinatesProduct()`, `DMStagGetProductCoordinateLocationSlot()`
150: @*/
151: PetscErrorCode DMStagGetProductCoordinateArraysRead(DM dm, void *arrX, void *arrY, void *arrZ)
152: {
153:   DMStagGetProductCoordinateArrays_Private(dm, arrX, arrY, arrZ, PETSC_TRUE);
154:   return 0;
155: }

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

160:   Not Collective

162:   Input Parameters:
163: + dm - the `DMSTAG` object
164: - loc - the grid location

166:   Output Parameter:
167: . slot - the index to use in local arrays

169:   Level: intermediate

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

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

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

183: .seealso: [](chapter_stag), `DMSTAG`, `DMPRODUCT`, `DMStagGetProductCoordinateArrays()`, `DMStagGetProductCoordinateArraysRead()`, `DMStagSetUniformCoordinates()`
184: @*/
185: PETSC_EXTERN PetscErrorCode DMStagGetProductCoordinateLocationSlot(DM dm, DMStagStencilLocation loc, PetscInt *slot)
186: {
187:   DM       dmCoord;
188:   PetscInt dim, dofCheck[DMSTAG_MAX_STRATA], s, d;

191:   DMGetDimension(dm, &dim);
192:   DMGetCoordinateDM(dm, &dmCoord);
194:   {
195:     PetscBool isProduct;
196:     DMType    dmType;
197:     DMGetType(dmCoord, &dmType);
198:     PetscStrcmp(DMPRODUCT, dmType, &isProduct);
200:   }
201:   for (s = 0; s < DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
202:   for (d = 0; d < dim; ++d) {
203:     DM        subDM;
204:     DMType    dmType;
205:     PetscBool isStag;
206:     PetscInt  dof[DMSTAG_MAX_STRATA], subDim;
207:     DMProductGetDM(dmCoord, d, &subDM);
209:     DMGetDimension(subDM, &subDim);
211:     DMGetType(subDM, &dmType);
212:     PetscStrcmp(DMSTAG, dmType, &isStag);
214:     DMStagGetDOF(subDM, &dof[0], &dof[1], &dof[2], &dof[3]);
215:     if (d == 0) {
216:       const PetscInt component = 0;
217:       for (s = 0; s < DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
218:       DMStagGetLocationSlot(subDM, loc, component, slot);
219:     } else {
221:     }
222:   }
223:   return 0;
224: }

226: /*@C
227:   DMStagGetCorners - return global element indices of the local region (excluding ghost points)

229:   Not Collective

231:   Input Parameter:
232: . dm - the `DMSTAG` object

234:   Output Parameters:
235: + x     - starting element index in first direction
236: . y     - starting element index in second direction
237: . z     - starting element index in third direction
238: . m     - element width in first direction
239: . n     - element width in second direction
240: . p     - element width in third direction
241: . nExtrax - number of extra partial elements in first direction
242: . nExtray - number of extra partial elements in second direction
243: - nExtraz - number of extra partial elements in third direction

245:   Level: beginner

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

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

254: .seealso: [](chapter_stag), `DMSTAG`, `DMStagGetGhostCorners()`, `DMDAGetCorners()`
255: @*/
256: PetscErrorCode DMStagGetCorners(DM dm, PetscInt *x, PetscInt *y, PetscInt *z, PetscInt *m, PetscInt *n, PetscInt *p, PetscInt *nExtrax, PetscInt *nExtray, PetscInt *nExtraz)
257: {
258:   const DM_Stag *const stag = (DM_Stag *)dm->data;

261:   if (x) *x = stag->start[0];
262:   if (y) *y = stag->start[1];
263:   if (z) *z = stag->start[2];
264:   if (m) *m = stag->n[0];
265:   if (n) *n = stag->n[1];
266:   if (p) *p = stag->n[2];
267:   if (nExtrax) *nExtrax = stag->boundaryType[0] != DM_BOUNDARY_PERIODIC && stag->lastRank[0] ? 1 : 0;
268:   if (nExtray) *nExtray = stag->boundaryType[1] != DM_BOUNDARY_PERIODIC && stag->lastRank[1] ? 1 : 0;
269:   if (nExtraz) *nExtraz = stag->boundaryType[2] != DM_BOUNDARY_PERIODIC && stag->lastRank[2] ? 1 : 0;
270:   return 0;
271: }

273: /*@C
274:   DMStagGetDOF - get number of DOF associated with each stratum of the grid

276:   Not Collective

278:   Input Parameter:
279: . dm - the `DMSTAG` object

281:   Output Parameters:
282: + dof0 - the number of points per 0-cell (vertex/node)
283: . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
284: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
285: - dof3 - the number of points per 3-cell (element in 3D)

287:   Level: beginner

289: .seealso: [](chapter_stag), `DMSTAG`, `DMStagGetCorners()`, `DMStagGetGhostCorners()`, `DMStagGetGlobalSizes()`, `DMStagGetStencilWidth()`, `DMStagGetBoundaryTypes()`, `DMStagGetLocationDOF()`, `DMDAGetDof()`
290: @*/
291: PetscErrorCode DMStagGetDOF(DM dm, PetscInt *dof0, PetscInt *dof1, PetscInt *dof2, PetscInt *dof3)
292: {
293:   const DM_Stag *const stag = (DM_Stag *)dm->data;

296:   if (dof0) *dof0 = stag->dof[0];
297:   if (dof1) *dof1 = stag->dof[1];
298:   if (dof2) *dof2 = stag->dof[2];
299:   if (dof3) *dof3 = stag->dof[3];
300:   return 0;
301: }

303: /*@C
304:   DMStagGetGhostCorners - return global element indices of the local region, including ghost points

306:   Not Collective

308:   Input Parameter:
309: . dm - the `DMSTAG` object

311:   Output Parameters:
312: + x - the starting element index in the first direction
313: . y - the starting element index in the second direction
314: . z - the starting element index in the third direction
315: . m - the element width in the first direction
316: . n - the element width in the second direction
317: - p - the element width in the third direction

319:   Level: beginner

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

324: .seealso: [](chapter_stag), `DMSTAG`, `DMStagGetCorners()`, `DMDAGetGhostCorners()`
325: @*/
326: PetscErrorCode DMStagGetGhostCorners(DM dm, PetscInt *x, PetscInt *y, PetscInt *z, PetscInt *m, PetscInt *n, PetscInt *p)
327: {
328:   const DM_Stag *const stag = (DM_Stag *)dm->data;

331:   if (x) *x = stag->startGhost[0];
332:   if (y) *y = stag->startGhost[1];
333:   if (z) *z = stag->startGhost[2];
334:   if (m) *m = stag->nGhost[0];
335:   if (n) *n = stag->nGhost[1];
336:   if (p) *p = stag->nGhost[2];
337:   return 0;
338: }

340: /*@C
341:   DMStagGetGlobalSizes - get global element counts

343:   Not Collective

345:   Input Parameter:
346: . dm - the `DMSTAG` object

348:   Output Parameters:
349: + M - global element counts in the x direction
350: . N - global element counts in the y direction
351: - P - global element counts in the z direction

353:   Level: beginner

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

358: .seealso: [](chapter_stag), `DMSTAG`, `DMStagGetLocalSizes()`, `DMDAGetInfo()`
359: @*/
360: PetscErrorCode DMStagGetGlobalSizes(DM dm, PetscInt *M, PetscInt *N, PetscInt *P)
361: {
362:   const DM_Stag *const stag = (DM_Stag *)dm->data;

365:   if (M) *M = stag->N[0];
366:   if (N) *N = stag->N[1];
367:   if (P) *P = stag->N[2];
368:   return 0;
369: }

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

374:   Not Collective

376:   Input Parameter:
377: . dm - the `DMSTAG` object

379:   Output Parameters:
380: + isFirstRank0 - whether this rank is first in the x direction
381: . isFirstRank1 - whether this rank is first in the y direction
382: - isFirstRank2 - whether this rank is first in the z direction

384:   Level: intermediate

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

389: .seealso: [](chapter_stag), `DMSTAG`, `DMStagGetIsLastRank()`
390: @*/
391: PetscErrorCode DMStagGetIsFirstRank(DM dm, PetscBool *isFirstRank0, PetscBool *isFirstRank1, PetscBool *isFirstRank2)
392: {
393:   const DM_Stag *const stag = (DM_Stag *)dm->data;

396:   if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
397:   if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
398:   if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
399:   return 0;
400: }

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

405:   Not Collective

407:   Input Parameter:
408: . dm - the `DMSTAG` object

410:   Output Parameters:
411: + isFirstRank0 - whether this rank is last in the x direction
412: . isFirstRank1 - whether this rank is last in the y direction
413: - isFirstRank2 - whether this rank is last in the z direction

415:   Level: intermediate

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

420: .seealso: [](chapter_stag), `DMSTAG`, `DMStagGetIsFirstRank()`
421: @*/
422: PetscErrorCode DMStagGetIsLastRank(DM dm, PetscBool *isLastRank0, PetscBool *isLastRank1, PetscBool *isLastRank2)
423: {
424:   const DM_Stag *const stag = (DM_Stag *)dm->data;

427:   if (isLastRank0) *isLastRank0 = stag->lastRank[0];
428:   if (isLastRank1) *isLastRank1 = stag->lastRank[1];
429:   if (isLastRank2) *isLastRank2 = stag->lastRank[2];
430:   return 0;
431: }

433: /*@C
434:   DMStagGetLocalSizes - get local elementwise sizes

436:   Not Collective

438:   Input Parameter:
439: . dm - the `DMSTAG` object

441:   Output Parameters:
442: + m - local element counts (excluding ghosts) in the x direction
443: . n - local element counts (excluding ghosts) in the y direction
444: - p - local element counts (excluding ghosts) in the z direction

446:   Level: beginner

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

451: .seealso: [](chapter_stag), `DMSTAG`, `DMStagGetGlobalSizes()`, `DMStagGetDOF()`, `DMStagGetNumRanks()`, `DMDAGetLocalInfo()`
452: @*/
453: PetscErrorCode DMStagGetLocalSizes(DM dm, PetscInt *m, PetscInt *n, PetscInt *p)
454: {
455:   const DM_Stag *const stag = (DM_Stag *)dm->data;

458:   if (m) *m = stag->n[0];
459:   if (n) *n = stag->n[1];
460:   if (p) *p = stag->n[2];
461:   return 0;
462: }

464: /*@C
465:   DMStagGetNumRanks - get number of ranks in each direction in the global grid decomposition

467:   Not Collective

469:   Input Parameter:
470: . dm - the `DMSTAG` object

472:   Output Parameters:
473: + nRanks0 - number of ranks in the x direction in the grid decomposition
474: . nRanks1 - number of ranks in the y direction in the grid decomposition
475: - nRanks2 - number of ranks in the z direction in the grid decomposition

477:  Level: intermediate

479: .seealso: [](chapter_stag), `DMSTAG`, `DMStagGetGlobalSizes()`, `DMStagGetLocalSize()`, `DMStagSetNumRanks()`, `DMDAGetInfo()`
480: @*/
481: PetscErrorCode DMStagGetNumRanks(DM dm, PetscInt *nRanks0, PetscInt *nRanks1, PetscInt *nRanks2)
482: {
483:   const DM_Stag *const stag = (DM_Stag *)dm->data;

486:   if (nRanks0) *nRanks0 = stag->nRanks[0];
487:   if (nRanks1) *nRanks1 = stag->nRanks[1];
488:   if (nRanks2) *nRanks2 = stag->nRanks[2];
489:   return 0;
490: }

492: /*@C
493:   DMStagGetEntries - get number of native entries in the global representation

495:   Not Collective

497:   Input Parameter:
498: . dm - the `DMSTAG` object

500:   Output Parameter:
501: . entries - number of rank-native entries in the global representation

503:   Level: developer

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

511: .seealso: [](chapter_stag), `DMSTAG`, `DMStagGetDOF()`, `DMStagGetEntriesLocal()`, `DMStagGetEntriesPerElement()`, `DMCreateLocalVector()`
512: @*/
513: PetscErrorCode DMStagGetEntries(DM dm, PetscInt *entries)
514: {
515:   const DM_Stag *const stag = (DM_Stag *)dm->data;

518:   if (entries) *entries = stag->entries;
519:   return 0;
520: }

522: /*@C
523:   DMStagGetEntriesLocal - get number of entries in the local representation

525:   Not Collective

527:   Input Parameter:
528: . dm - the `DMSTAG` object

530:   Output Parameter:
531: . entries - number of entries in the local representation

533:   Level: developer

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

541: .seealso: [](chapter_stag), DMSTAG, DMStagGetDOF(), DMStagGetEntries(), DMStagGetEntriesPerElement(), DMCreateLocalVector()
542: @*/
543: PetscErrorCode DMStagGetEntriesLocal(DM dm, PetscInt *entries)
544: {
545:   const DM_Stag *const stag = (DM_Stag *)dm->data;

548:   if (entries) *entries = stag->entriesGhost;
549:   return 0;
550: }

552: /*@C
553:   DMStagGetEntriesPerElement - get number of entries per element in the local representation

555:   Not Collective

557:   Input Parameter:
558: . dm - the `DMSTAG` object

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

563:   Level: developer

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

569: .seealso: [](chapter_stag), `DMSTAG`, `DMStagGetDOF()`
570: @*/
571: PetscErrorCode DMStagGetEntriesPerElement(DM dm, PetscInt *entriesPerElement)
572: {
573:   const DM_Stag *const stag = (DM_Stag *)dm->data;

576:   if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
577:   return 0;
578: }

580: /*@C
581:   DMStagGetStencilType - get elementwise ghost/halo stencil type

583:   Not Collective

585:   Input Parameter:
586: . dm - the `DMSTAG` object

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

591:   Level: beginner

593: .seealso: [](chapter_stag), `DMSTAG`, `DMStagSetStencilType()`, `DMStagGetStencilWidth`, `DMStagStencilType`
594: @*/
595: PetscErrorCode DMStagGetStencilType(DM dm, DMStagStencilType *stencilType)
596: {
597:   DM_Stag *const stag = (DM_Stag *)dm->data;

600:   *stencilType = stag->stencilType;
601:   return 0;
602: }

604: /*@C
605:   DMStagGetStencilWidth - get elementwise stencil width

607:   Not Collective

609:   Input Parameter:
610: . dm - the `DMSTAG` object

612:   Output Parameter:
613: . stencilWidth - stencil/halo/ghost width in elements

615:   Level: beginner

617: .seealso: [](chapter_stag), `DMSTAG`, `DMStagSetStencilWidth()`, `DMStagGetStencilType()`, `DMDAGetStencilType()`
618: @*/
619: PetscErrorCode DMStagGetStencilWidth(DM dm, PetscInt *stencilWidth)
620: {
621:   const DM_Stag *const stag = (DM_Stag *)dm->data;

624:   if (stencilWidth) *stencilWidth = stag->stencilWidth;
625:   return 0;
626: }

628: /*@C
629:   DMStagGetOwnershipRanges - get elements per rank in each direction

631:   Not Collective

633:   Input Parameter:
634: .     dm - the `DMSTAG` object

636:   Output Parameters:
637: +     lx - ownership along x direction (optional)
638: .     ly - ownership along y direction (optional)
639: -     lz - ownership along z direction (optional)

641:   Level: intermediate

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

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

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

651: .seealso: [](chapter_stag), `DMSTAG`, `DMStagSetGlobalSizes()`, `DMStagSetOwnershipRanges()`, `DMStagCreate1d()`, `DMStagCreate2d()`, `DMStagCreate3d()`, `DMDAGetOwnershipRanges()`
652: @*/
653: PetscErrorCode DMStagGetOwnershipRanges(DM dm, const PetscInt *lx[], const PetscInt *ly[], const PetscInt *lz[])
654: {
655:   const DM_Stag *const stag = (DM_Stag *)dm->data;

658:   if (lx) *lx = stag->l[0];
659:   if (ly) *ly = stag->l[1];
660:   if (lz) *lz = stag->l[2];
661:   return 0;
662: }

664: /*@C
665:   DMStagCreateCompatibleDMStag - create a compatible `DMSTAG` with different dof/stratum

667:   Collective

669:   Input Parameters:
670: + dm - the `DMSTAG` object
671: . dof0 - number of dof on the first stratum in the new `DMSTAG`
672: . dof1 - number of dof on the second stratum in the new `DMSTAG`
673: . dof2 - number of dof on the third stratum in the new `DMSTAG`
674: - dof3 - number of dof on the fourth stratum in the new `DMSTAG`

676:   Output Parameter:
677: . newdm - the new, compatible `DMSTAG`

679:   Level: intermediate

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

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

688: .seealso: [](chapter_stag), `DMSTAG`, `DMDACreateCompatibleDMDA()`, `DMGetCompatibility()`, `DMStagMigrateVec()`
689: @*/
690: PetscErrorCode DMStagCreateCompatibleDMStag(DM dm, PetscInt dof0, PetscInt dof1, PetscInt dof2, PetscInt dof3, DM *newdm)
691: {
693:   DMStagDuplicateWithoutSetup(dm, PetscObjectComm((PetscObject)dm), newdm);
694:   DMStagSetDOF(*newdm, dof0, dof1, dof2, dof3);
695:   DMSetUp(*newdm);
696:   return 0;
697: }

699: /*@C
700:   DMStagGetLocationSlot - get index to use in accessing raw local arrays

702:   Not Collective

704:   Input Parameters:
705: + dm - the `DMSTAG` object
706: . loc - location relative to an element
707: - c - component

709:   Output Parameter:
710: . slot - index to use

712:   Level: beginner

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

719: .seealso: [](chapter_stag), `DMSTAG`, `DMStagVecGetArray()`, `DMStagVecGetArrayRead()`, `DMStagGetDOF()`, `DMStagGetEntriesPerElement()`
720: @*/
721: PetscErrorCode DMStagGetLocationSlot(DM dm, DMStagStencilLocation loc, PetscInt c, PetscInt *slot)
722: {
723:   DM_Stag *const stag = (DM_Stag *)dm->data;

726:   if (PetscDefined(USE_DEBUG)) {
727:     PetscInt dof;
728:     DMStagGetLocationDOF(dm, loc, &dof);
731:   }
732:   *slot = stag->locationOffsets[loc] + c;
733:   return 0;
734: }

736: /*@C
737:   DMStagMigrateVec - transfer a vector associated with a `DMSTAG` to a vector associated with a compatible `DMSTAG`

739:   Collective

741:   Input Parameters:
742: + dm - the source `DMSTAG` object
743: . vec - the source vector, compatible with `dm`
744: . dmTo - the compatible destination `DMSTAG` object
745: - vecTo - the destination vector, compatible with `dmTo`

747:   Level: advanced

749:   Notes:
750:   Extra dof are ignored, and unfilled dof are zeroed.
751:   Currently only implemented to migrate global vectors to global vectors.
752:   For the definition of compatibility of `DM`s, see `DMGetCompatibility()`.

754: .seealso: [](chapter_stag), `DMSTAG`, `DMStagCreateCompatibleDMStag()`, `DMGetCompatibility()`, `DMStagVecSplitToDMDA()`
755: @*/
756: PetscErrorCode DMStagMigrateVec(DM dm, Vec vec, DM dmTo, Vec vecTo)
757: {
758:   DM_Stag *const     stag   = (DM_Stag *)dm->data;
759:   DM_Stag *const     stagTo = (DM_Stag *)dmTo->data;
760:   PetscInt           nLocalTo, nLocal, dim, i, j, k;
761:   PetscInt           start[DMSTAG_MAX_DIM], startGhost[DMSTAG_MAX_DIM], n[DMSTAG_MAX_DIM], nExtra[DMSTAG_MAX_DIM], offset[DMSTAG_MAX_DIM];
762:   Vec                vecToLocal, vecLocal;
763:   PetscBool          compatible, compatibleSet;
764:   const PetscScalar *arr;
765:   PetscScalar       *arrTo;
766:   const PetscInt     epe   = stag->entriesPerElement;
767:   const PetscInt     epeTo = stagTo->entriesPerElement;

773:   DMGetCompatibility(dm, dmTo, &compatible, &compatibleSet);
775:   DMGetDimension(dm, &dim);
776:   VecGetLocalSize(vec, &nLocal);
777:   VecGetLocalSize(vecTo, &nLocalTo);
779:   DMStagGetCorners(dm, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], &nExtra[0], &nExtra[1], &nExtra[2]);
780:   DMStagGetGhostCorners(dm, &startGhost[0], &startGhost[1], &startGhost[2], NULL, NULL, NULL);
781:   for (i = 0; i < DMSTAG_MAX_DIM; ++i) offset[i] = start[i] - startGhost[i];

783:   /* Proceed by transferring to a local vector, copying, and transferring back to a global vector */
784:   DMGetLocalVector(dm, &vecLocal);
785:   DMGetLocalVector(dmTo, &vecToLocal);
786:   DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal);
787:   DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal);
788:   VecGetArrayRead(vecLocal, &arr);
789:   VecGetArray(vecToLocal, &arrTo);
790:   /* Note that some superfluous copying of entries on partial dummy elements is done */
791:   if (dim == 1) {
792:     for (i = offset[0]; i < offset[0] + n[0] + nExtra[0]; ++i) {
793:       PetscInt d = 0, dTo = 0, b = 0, bTo = 0;
794:       PetscInt si;
795:       for (si = 0; si < 2; ++si) {
796:         b += stag->dof[si];
797:         bTo += stagTo->dof[si];
798:         for (; d < b && dTo < bTo; ++d, ++dTo) arrTo[i * epeTo + dTo] = arr[i * epe + d];
799:         for (; dTo < bTo; ++dTo) arrTo[i * epeTo + dTo] = 0.0;
800:         d = b;
801:       }
802:     }
803:   } else if (dim == 2) {
804:     const PetscInt epr   = stag->nGhost[0] * epe;
805:     const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
806:     for (j = offset[1]; j < offset[1] + n[1] + nExtra[1]; ++j) {
807:       for (i = offset[0]; i < offset[0] + n[0] + nExtra[0]; ++i) {
808:         const PetscInt base   = j * epr + i * epe;
809:         const PetscInt baseTo = j * eprTo + i * epeTo;
810:         PetscInt       d = 0, dTo = 0, b = 0, bTo = 0;
811:         const PetscInt s[4] = {0, 1, 1, 2}; /* Dimensions of points, in order */
812:         PetscInt       si;
813:         for (si = 0; si < 4; ++si) {
814:           b += stag->dof[s[si]];
815:           bTo += stagTo->dof[s[si]];
816:           for (; d < b && dTo < bTo; ++d, ++dTo) arrTo[baseTo + dTo] = arr[base + d];
817:           for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
818:           d = b;
819:         }
820:       }
821:     }
822:   } else if (dim == 3) {
823:     const PetscInt epr   = stag->nGhost[0] * epe;
824:     const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
825:     const PetscInt epl   = stag->nGhost[1] * epr;
826:     const PetscInt eplTo = stagTo->nGhost[1] * eprTo;
827:     for (k = offset[2]; k < offset[2] + n[2] + nExtra[2]; ++k) {
828:       for (j = offset[1]; j < offset[1] + n[1] + nExtra[1]; ++j) {
829:         for (i = offset[0]; i < offset[0] + n[0] + nExtra[0]; ++i) {
830:           PetscInt       d = 0, dTo = 0, b = 0, bTo = 0;
831:           const PetscInt base   = k * epl + j * epr + i * epe;
832:           const PetscInt baseTo = k * eplTo + j * eprTo + i * epeTo;
833:           const PetscInt s[8]   = {0, 1, 1, 2, 1, 2, 2, 3}; /* dimensions of points, in order */
834:           PetscInt       is;
835:           for (is = 0; is < 8; ++is) {
836:             b += stag->dof[s[is]];
837:             bTo += stagTo->dof[s[is]];
838:             for (; d < b && dTo < bTo; ++d, ++dTo) arrTo[baseTo + dTo] = arr[base + d];
839:             for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
840:             d = b;
841:           }
842:         }
843:       }
844:     }
845:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
846:   VecRestoreArrayRead(vecLocal, &arr);
847:   VecRestoreArray(vecToLocal, &arrTo);
848:   DMRestoreLocalVector(dm, &vecLocal);
849:   DMLocalToGlobalBegin(dmTo, vecToLocal, INSERT_VALUES, vecTo);
850:   DMLocalToGlobalEnd(dmTo, vecToLocal, INSERT_VALUES, vecTo);
851:   DMRestoreLocalVector(dmTo, &vecToLocal);
852:   return 0;
853: }

855: /*@C
856:   DMStagPopulateLocalToGlobalInjective - populate an internal 1-to-1 local-to-global map

858:   Collective

860:   Creates an internal object which explicitly maps a single local degree of
861:   freedom to each global degree of freedom. This is used, if populated,
862:   instead of SCATTER_REVERSE_LOCAL with the (1-to-many, in general)
863:   global-to-local map, when DMLocalToGlobal() is called with INSERT_VALUES.
864:   This allows usage, for example, even in the periodic, 1-rank case, where
865:   the inverse of the global-to-local map, even when restricted to on-rank
866:   communication, is non-injective. This is at the cost of storing an additional
867:   VecScatter object inside each `DMSTAG` object.

869:   Input Parameter:
870: . dm - the `DMSTAG` object

872:   Level: developer

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

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

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

885: .seealso: [](chapter_stag), `DMSTAG`, `DMLocalToGlobal()`, `VecScatter`
886: @*/
887: PetscErrorCode DMStagPopulateLocalToGlobalInjective(DM dm)
888: {
889:   PetscInt       dim;
890:   DM_Stag *const stag = (DM_Stag *)dm->data;

893:   if (stag->ltog_injective) return 0; /* Don't re-populate */
894:   DMGetDimension(dm, &dim);
895:   switch (dim) {
896:   case 1:
897:     DMStagPopulateLocalToGlobalInjective_1d(dm);
898:     break;
899:   case 2:
900:     DMStagPopulateLocalToGlobalInjective_2d(dm);
901:     break;
902:   case 3:
903:     DMStagPopulateLocalToGlobalInjective_3d(dm);
904:     break;
905:   default:
906:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
907:   }
908:   return 0;
909: }

911: static PetscErrorCode DMStagRestoreProductCoordinateArrays_Private(DM dm, void *arrX, void *arrY, void *arrZ, PetscBool read)
912: {
913:   PetscInt dim, d;
914:   void    *arr[DMSTAG_MAX_DIM];
915:   DM       dmCoord;

918:   DMGetDimension(dm, &dim);
920:   arr[0] = arrX;
921:   arr[1] = arrY;
922:   arr[2] = arrZ;
923:   DMGetCoordinateDM(dm, &dmCoord);
924:   for (d = 0; d < dim; ++d) {
925:     DM  subDM;
926:     Vec coord1d_local;

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

931:     DMProductGetDM(dmCoord, d, &subDM);
932:     DMGetCoordinatesLocal(subDM, &coord1d_local);
933:     if (read) {
934:       DMStagVecRestoreArrayRead(subDM, coord1d_local, arr[d]);
935:     } else {
936:       DMStagVecRestoreArray(subDM, coord1d_local, arr[d]);
937:     }
938:   }
939:   return 0;
940: }

942: /*@C
943:   DMStagRestoreProductCoordinateArrays - restore local array access

945:   Logically Collective

947:   Input Parameter:
948: . dm - the `DMSTAG` object

950:   Output Parameters:
951: + arrX - local 1D coordinate arrays for x direction
952: . arrY - local 1D coordinate arrays for y direction
953: - arrZ - local 1D coordinate arrays for z direction

955:   Level: intermediate

957:   Notes:
958:   This function does not automatically perform a local->global scatter to populate global coordinates from the local coordinates.
959:   Thus, it may be required to explicitly perform these operations in some situations, as in the following partial example:

961:   ```
962:   DMGetCoordinateDM(dm,&cdm);
963:   for (d=0; d<3; ++d) {
964:     DM  subdm;
965:     Vec coor,coor_local;

967:     DMProductGetDM(cdm,d,&subdm);
968:     DMGetCoordinates(subdm,&coor);
969:     DMGetCoordinatesLocal(subdm,&coor_local);
970:     DMLocalToGlobal(subdm,coor_local,INSERT_VALUES,coor);
971:     PetscPrintf(PETSC_COMM_WORLD,"Coordinates dim %" PetscInt_FMT ":\n",d);
972:     VecView(coor,PETSC_VIEWER_STDOUT_WORLD);
973:   }
974:    ```

976: .seealso: [](chapter_stag), `DMSTAG`, `DMStagGetProductCoordinateArrays()`, `DMStagGetProductCoordinateArraysRead()`
977: @*/
978: PetscErrorCode DMStagRestoreProductCoordinateArrays(DM dm, void *arrX, void *arrY, void *arrZ)
979: {
980:   DMStagRestoreProductCoordinateArrays_Private(dm, arrX, arrY, arrZ, PETSC_FALSE);
981:   return 0;
982: }

984: /*@C
985:   DMStagRestoreProductCoordinateArraysRead - restore local product array access, read-only

987:   Logically Collective

989:   Input Parameters:
990: + dm - the `DMSTAG` object
991: . arrX - local 1D coordinate arrays for x direction
992: . arrY - local 1D coordinate arrays for y direction
993: - arrZ - local 1D coordinate arrays for z direction

995:   Level: intermediate

997: .seealso: [](chapter_stag), `DMSTAG`, `DMStagGetProductCoordinateArrays()`, `DMStagGetProductCoordinateArraysRead()`
998: @*/
999: PetscErrorCode DMStagRestoreProductCoordinateArraysRead(DM dm, void *arrX, void *arrY, void *arrZ)
1000: {
1001:   DMStagRestoreProductCoordinateArrays_Private(dm, arrX, arrY, arrZ, PETSC_TRUE);
1002:   return 0;
1003: }

1005: /*@C
1006:   DMStagSetBoundaryTypes - set `DMSTAG` boundary types

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

1010:   Input Parameters:
1011: + dm - the `DMSTAG` object
1012: . boundaryTypeX - boundary type for x direction
1013: . boundaryTypeY - boundary type for y direction, not set for one dimensional problems
1014: - boundaryTypeZ - boundary type for z direction, not set for one and two dimensional problems

1016:   Level: advanced

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

1021: .seealso: [](chapter_stag), `DMSTAG`, `DMBoundaryType`, `DMStagCreate1d()`, `DMStagCreate2d()`, `DMStagCreate3d()`, `DMDASetBoundaryType()`
1022: @*/
1023: PetscErrorCode DMStagSetBoundaryTypes(DM dm, DMBoundaryType boundaryType0, DMBoundaryType boundaryType1, DMBoundaryType boundaryType2)
1024: {
1025:   DM_Stag *const stag = (DM_Stag *)dm->data;
1026:   PetscInt       dim;

1028:   DMGetDimension(dm, &dim);
1034:   stag->boundaryType[0] = boundaryType0;
1035:   if (dim > 1) stag->boundaryType[1] = boundaryType1;
1036:   if (dim > 2) stag->boundaryType[2] = boundaryType2;
1037:   return 0;
1038: }

1040: /*@C
1041:   DMStagSetCoordinateDMType - set DM type to store coordinates

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

1045:   Input Parameters:
1046: + dm - the `DMSTAG` object
1047: - dmtype - DMtype for coordinates, either `DMSTAG` or `DMPRODUCT`

1049:   Level: advanced

1051: .seealso: [](chapter_stag), `DMSTAG`, `DMPRODUCT`, `DMGetCoordinateDM()`, `DMStagSetUniformCoordinates()`, `DMStagSetUniformCoordinatesExplicit()`, `DMStagSetUniformCoordinatesProduct()`, `DMType`
1052: @*/
1053: PetscErrorCode DMStagSetCoordinateDMType(DM dm, DMType dmtype)
1054: {
1055:   DM_Stag *const stag = (DM_Stag *)dm->data;

1058:   PetscFree(stag->coordinateDMType);
1059:   PetscStrallocpy(dmtype, (char **)&stag->coordinateDMType);
1060:   return 0;
1061: }

1063: /*@C
1064:   DMStagSetDOF - set dof/stratum

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

1068:   Input Parameters:
1069: + dm - the `DMSTAG` object
1070: . dof0 - the number of points per 0-cell (vertex/node)
1071: . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
1072: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
1073: - dof3 - the number of points per 3-cell (element in 3D)

1075:   Level: advanced

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

1080: .seealso: [](chapter_stag), `DMSTAG`, `DMDASetDof()`
1081: @*/
1082: PetscErrorCode DMStagSetDOF(DM dm, PetscInt dof0, PetscInt dof1, PetscInt dof2, PetscInt dof3)
1083: {
1084:   DM_Stag *const stag = (DM_Stag *)dm->data;
1085:   PetscInt       dim;

1093:   DMGetDimension(dm, &dim);
1098:   stag->dof[0] = dof0;
1099:   stag->dof[1] = dof1;
1100:   if (dim > 1) stag->dof[2] = dof2;
1101:   if (dim > 2) stag->dof[3] = dof3;
1102:   return 0;
1103: }

1105: /*@C
1106:   DMStagSetNumRanks - set ranks in each direction in the global rank grid

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

1110:   Input Parameters:
1111: + dm - the `DMSTAG` object
1112: . nRanks0 - number of ranks in the x direction
1113: . nRanks1 - number of ranks in the y direction
1114: - nRanks2 - number of ranks in the z direction

1116:   Level: developer

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

1121: .seealso: [](chapter_stag), `DMSTAG`, `DMDASetNumProcs()`
1122: @*/
1123: PetscErrorCode DMStagSetNumRanks(DM dm, PetscInt nRanks0, PetscInt nRanks1, PetscInt nRanks2)
1124: {
1125:   DM_Stag *const stag = (DM_Stag *)dm->data;
1126:   PetscInt       dim;

1133:   DMGetDimension(dm, &dim);
1137:   if (nRanks0) stag->nRanks[0] = nRanks0;
1138:   if (dim > 1 && nRanks1) stag->nRanks[1] = nRanks1;
1139:   if (dim > 2 && nRanks2) stag->nRanks[2] = nRanks2;
1140:   return 0;
1141: }

1143: /*@C
1144:   DMStagSetStencilType - set elementwise ghost/halo stencil type

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

1148:   Input Parameters:
1149: + dm - the `DMSTAG` object
1150: - stencilType - the elementwise ghost stencil type: `DMSTAG_STENCIL_BOX`, `DMSTAG_STENCIL_STAR`, or `DMSTAG_STENCIL_NONE`

1152:   Level: beginner

1154: .seealso: [](chapter_stag), `DMSTAG`, `DMStagGetStencilType()`, `DMStagSetStencilWidth()`, `DMStagStencilType`
1155: @*/
1156: PetscErrorCode DMStagSetStencilType(DM dm, DMStagStencilType stencilType)
1157: {
1158:   DM_Stag *const stag = (DM_Stag *)dm->data;

1163:   stag->stencilType = stencilType;
1164:   return 0;
1165: }

1167: /*@C
1168:   DMStagSetStencilWidth - set elementwise stencil width

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

1172:   Input Parameters:
1173: + dm - the `DMSTAG` object
1174: - stencilWidth - stencil/halo/ghost width in elements

1176:   Level: beginner

1178:   Note:
1179:   The width value is not used when `DMSTAG_STENCIL_NONE` is specified.

1181: .seealso: [](chapter_stag), `DMSTAG`, `DMStagGetStencilWidth()`, `DMStagGetStencilType()`, `DMStagStencilType`
1182: @*/
1183: PetscErrorCode DMStagSetStencilWidth(DM dm, PetscInt stencilWidth)
1184: {
1185:   DM_Stag *const stag = (DM_Stag *)dm->data;

1191:   stag->stencilWidth = stencilWidth;
1192:   return 0;
1193: }

1195: /*@C
1196:   DMStagSetGlobalSizes - set global element counts in each direction

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

1200:   Input Parameters:
1201: + dm - the `DMSTAG` object
1202: . N0 - global elementwise size in the x direction
1203: . N1 - global elementwise size in the y direction
1204: - N2 - global elementwise size in the z direction

1206:   Level: advanced

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

1211: .seealso: [](chapter_stag), `DMSTAG`, `DMStagGetGlobalSizes()`, `DMDASetSizes()`
1212: @*/
1213: PetscErrorCode DMStagSetGlobalSizes(DM dm, PetscInt N0, PetscInt N1, PetscInt N2)
1214: {
1215:   DM_Stag *const stag = (DM_Stag *)dm->data;
1216:   PetscInt       dim;

1220:   DMGetDimension(dm, &dim);
1224:   if (N0) stag->N[0] = N0;
1225:   if (N1) stag->N[1] = N1;
1226:   if (N2) stag->N[2] = N2;
1227:   return 0;
1228: }

1230: /*@C
1231:   DMStagSetOwnershipRanges - set elements per rank in each direction

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

1235:   Input Parameters:
1236: + dm - the `DMSTAG` object
1237: . lx - element counts for each rank in the x direction
1238: . ly - element counts for each rank in the y direction
1239: - lz - element counts for each rank in the z direction

1241:   Level: developer

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

1246: .seealso: [](chapter_stag), `DMSTAG`, `DMStagSetGlobalSizes()`, `DMStagGetOwnershipRanges()`, `DMDASetOwnershipRanges()`
1247: @*/
1248: PetscErrorCode DMStagSetOwnershipRanges(DM dm, PetscInt const *lx, PetscInt const *ly, PetscInt const *lz)
1249: {
1250:   DM_Stag *const  stag = (DM_Stag *)dm->data;
1251:   const PetscInt *lin[3];
1252:   PetscInt        d, dim;

1256:   lin[0] = lx;
1257:   lin[1] = ly;
1258:   lin[2] = lz;
1259:   DMGetDimension(dm, &dim);
1260:   for (d = 0; d < dim; ++d) {
1261:     if (lin[d]) {
1263:       if (!stag->l[d]) PetscMalloc1(stag->nRanks[d], &stag->l[d]);
1264:       PetscArraycpy(stag->l[d], lin[d], stag->nRanks[d]);
1265:     }
1266:   }
1267:   return 0;
1268: }

1270: /*@C
1271:   DMStagSetUniformCoordinates - set `DMSTAG` coordinates to be a uniform grid

1273:   Collective

1275:   Input Parameters:
1276: + dm - the `DMSTAG` object
1277: . xmin - minimum global coordinate value in the x direction
1278: . xmax - maximum global coordinate values in the x direction
1279: . ymin - minimum global coordinate value in the y direction
1280: . ymax - maximum global coordinate value in the y direction
1281: . zmin - minimum global coordinate value in the z direction
1282: - zmax - maximum global coordinate value in the z direction

1284:   Level: advanced

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

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

1298: .seealso: [](chapter_stag), `DMSTAG`, `DMPRODUCT`, `DMStagSetUniformCoordinatesExplicit()`, `DMStagSetUniformCoordinatesProduct()`, `DMStagSetCoordinateDMType()`, `DMGetCoordinateDM()`, `DMGetCoordinates()`, `DMDASetUniformCoordinates()`, `DMBoundaryType`
1299: @*/
1300: PetscErrorCode DMStagSetUniformCoordinates(DM dm, PetscReal xmin, PetscReal xmax, PetscReal ymin, PetscReal ymax, PetscReal zmin, PetscReal zmax)
1301: {
1302:   DM_Stag *const stag = (DM_Stag *)dm->data;
1303:   PetscBool      flg_stag, flg_product;

1308:   PetscStrcmp(stag->coordinateDMType, DMSTAG, &flg_stag);
1309:   PetscStrcmp(stag->coordinateDMType, DMPRODUCT, &flg_product);
1310:   if (flg_stag) {
1311:     DMStagSetUniformCoordinatesExplicit(dm, xmin, xmax, ymin, ymax, zmin, zmax);
1312:   } else if (flg_product) {
1313:     DMStagSetUniformCoordinatesProduct(dm, xmin, xmax, ymin, ymax, zmin, zmax);
1314:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported DM Type %s", stag->coordinateDMType);
1315:   return 0;
1316: }

1318: /*@C
1319:   DMStagSetUniformCoordinatesExplicit - set `DMSTAG` coordinates to be a uniform grid, storing all values

1321:   Collective

1323:   Input Parameters:
1324: + dm - the `DMSTAG` object
1325: . xmin - minimum global coordinate value in the x direction
1326: . xmax - maximum global coordinate values in the x direction
1327: . ymin - minimum global coordinate value in the y direction
1328: . ymax - maximum global coordinate value in the y direction
1329: . zmin - minimum global coordinate value in the z direction
1330: - zmax - maximum global coordinate value in the z direction

1332:   Level: beginner

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

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

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

1343: .seealso: [](chapter_stag), `DMSTAG`, `DMStagSetUniformCoordinates()`, `DMStagSetUniformCoordinatesProduct()`, `DMStagSetCoordinateDMType()`
1344: @*/
1345: PetscErrorCode DMStagSetUniformCoordinatesExplicit(DM dm, PetscReal xmin, PetscReal xmax, PetscReal ymin, PetscReal ymax, PetscReal zmin, PetscReal zmax)
1346: {
1347:   DM_Stag *const stag = (DM_Stag *)dm->data;
1348:   PetscInt       dim;
1349:   PetscBool      flg;

1353:   PetscStrcmp(stag->coordinateDMType, DMSTAG, &flg);
1355:   DMStagSetCoordinateDMType(dm, DMSTAG);
1356:   DMGetDimension(dm, &dim);
1357:   switch (dim) {
1358:   case 1:
1359:     DMStagSetUniformCoordinatesExplicit_1d(dm, xmin, xmax);
1360:     break;
1361:   case 2:
1362:     DMStagSetUniformCoordinatesExplicit_2d(dm, xmin, xmax, ymin, ymax);
1363:     break;
1364:   case 3:
1365:     DMStagSetUniformCoordinatesExplicit_3d(dm, xmin, xmax, ymin, ymax, zmin, zmax);
1366:     break;
1367:   default:
1368:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
1369:   }
1370:   DMCoarsenHookRemove(dm, DMRestrictHook_Coordinates, NULL, NULL);
1371:   return 0;
1372: }

1374: /*@C
1375:   DMStagSetUniformCoordinatesProduct - create uniform coordinates, as a product of 1D arrays

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

1379:   Collective

1381:   Input Parameters:
1382: + dm - the `DMSTAG` object
1383: . xmin - minimum global coordinate value in the x direction
1384: . xmax - maximum global coordinate values in the x direction
1385: . ymin - minimum global coordinate value in the y direction
1386: . ymax - maximum global coordinate value in the y direction
1387: . zmin - minimum global coordinate value in the z direction
1388: - zmax - maximum global coordinate value in the z direction

1390:   Level: intermediate

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

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

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

1402: .seealso: [](chapter_stag), `DMSTAG`, `DMPRODUCT`, `DMStagSetUniformCoordinates()`, `DMStagSetUniformCoordinatesExplicit()`, `DMStagSetCoordinateDMType()`
1403: @*/
1404: PetscErrorCode DMStagSetUniformCoordinatesProduct(DM dm, PetscReal xmin, PetscReal xmax, PetscReal ymin, PetscReal ymax, PetscReal zmin, PetscReal zmax)
1405: {
1406:   DM_Stag *const stag = (DM_Stag *)dm->data;
1407:   DM             dmc;
1408:   PetscInt       dim, d, dof0, dof1;
1409:   PetscBool      flg;

1413:   PetscStrcmp(stag->coordinateDMType, DMPRODUCT, &flg);
1415:   DMStagSetCoordinateDMType(dm, DMPRODUCT);
1416:   DMGetCoordinateDM(dm, &dmc);
1417:   DMGetDimension(dm, &dim);

1419:   /* Create 1D sub-DMs, living on subcommunicators.
1420:      Always include both vertex and element dof, regardless of the active strata of the DMStag */
1421:   dof0 = 1;
1422:   dof1 = 1;

1424:   for (d = 0; d < dim; ++d) {
1425:     DM                subdm;
1426:     MPI_Comm          subcomm;
1427:     PetscMPIInt       color;
1428:     const PetscMPIInt key = 0; /* let existing rank break ties */

1430:     /* Choose colors based on position in the plane orthogonal to this dim, and split */
1431:     switch (d) {
1432:     case 0:
1433:       color = (dim > 1 ? stag->rank[1] : 0) + (dim > 2 ? stag->nRanks[1] * stag->rank[2] : 0);
1434:       break;
1435:     case 1:
1436:       color = stag->rank[0] + (dim > 2 ? stag->nRanks[0] * stag->rank[2] : 0);
1437:       break;
1438:     case 2:
1439:       color = stag->rank[0] + stag->nRanks[0] * stag->rank[1];
1440:       break;
1441:     default:
1442:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension index %" PetscInt_FMT, d);
1443:     }
1444:     MPI_Comm_split(PetscObjectComm((PetscObject)dm), color, key, &subcomm);

1446:     /* Create sub-DMs living on these new communicators (which are destroyed by DMProduct) */
1447:     DMStagCreate1d(subcomm, stag->boundaryType[d], stag->N[d], dof0, dof1, stag->stencilType, stag->stencilWidth, stag->l[d], &subdm);
1448:     /* Copy vector and matrix type information from parent DM */
1449:     DMSetVecType(subdm, dm->vectype);
1450:     DMSetMatType(subdm, dm->mattype);
1451:     DMSetUp(subdm);
1452:     switch (d) {
1453:     case 0:
1454:       DMStagSetUniformCoordinatesExplicit(subdm, xmin, xmax, 0, 0, 0, 0);
1455:       break;
1456:     case 1:
1457:       DMStagSetUniformCoordinatesExplicit(subdm, ymin, ymax, 0, 0, 0, 0);
1458:       break;
1459:     case 2:
1460:       DMStagSetUniformCoordinatesExplicit(subdm, zmin, zmax, 0, 0, 0, 0);
1461:       break;
1462:     default:
1463:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension index %" PetscInt_FMT, d);
1464:     }
1465:     DMProductSetDM(dmc, d, subdm);
1466:     DMProductSetDimensionIndex(dmc, d, 0);
1467:     DMDestroy(&subdm);
1468:     MPI_Comm_free(&subcomm);
1469:   }
1470:   DMCoarsenHookRemove(dm, DMRestrictHook_Coordinates, NULL, NULL);
1471:   return 0;
1472: }

1474: /*@C
1475:   DMStagVecGetArray - get access to local array

1477:   Logically Collective

1479:   Input Parameters:
1480: + dm - the `DMSTAG` object
1481: - vec - the `Vec` object

1483:   Output Parameter:
1484: . array - the array

1486:   Level: beginner

1488:   Note:
1489:   This function returns a (dim+1)-dimensional array for a dim-dimensional
1490:   `DMSTAG`.

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

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

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

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

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

1509: .seealso: [](chapter_stag), `DMSTAG`, `DMStagVecGetArrayRead()`, `DMStagGetLocationSlot()`, `DMGetLocalVector()`, `DMCreateLocalVector()`, `DMGetGlobalVector()`, `DMCreateGlobalVector()`, `DMDAVecGetArray()`, `DMDAVecGetArrayDOF()`
1510: @*/
1511: PetscErrorCode DMStagVecGetArray(DM dm, Vec vec, void *array)
1512: {
1513:   DM_Stag *const stag = (DM_Stag *)dm->data;
1514:   PetscInt       dim;
1515:   PetscInt       nLocal;

1519:   DMGetDimension(dm, &dim);
1520:   VecGetLocalSize(vec, &nLocal);
1522:   switch (dim) {
1523:   case 1:
1524:     VecGetArray2d(vec, stag->nGhost[0], stag->entriesPerElement, stag->startGhost[0], 0, (PetscScalar ***)array);
1525:     break;
1526:   case 2:
1527:     VecGetArray3d(vec, stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar ****)array);
1528:     break;
1529:   case 3:
1530:     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);
1531:     break;
1532:   default:
1533:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
1534:   }
1535:   return 0;
1536: }

1538: /*@C
1539:   DMStagVecGetArrayRead - get read-only access to a local array

1541:   Logically Collective

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

1545:   Input Parameters:
1546: + dm - the `DMSTAG` object
1547: - vec - the `Vec` object

1549:   Output Parameter:
1550: . array - the read-only array

1552:   Level: beginner

1554:   Note:
1555:   `DMStagVecRestoreArrayRead()` must be called, once finished with the array

1557: .seealso: [](chapter_stag), `DMSTAG`, `DMStagVecGetArrayRead()`, `DMStagGetLocationSlot()`, `DMGetLocalVector()`, `DMCreateLocalVector()`, `DMGetGlobalVector()`, `DMCreateGlobalVector()`, `DMDAVecGetArrayRead()`, `DMDAVecGetArrayDOFRead()`
1558: @*/
1559: PetscErrorCode DMStagVecGetArrayRead(DM dm, Vec vec, void *array)
1560: {
1561:   DM_Stag *const stag = (DM_Stag *)dm->data;
1562:   PetscInt       dim;
1563:   PetscInt       nLocal;

1567:   DMGetDimension(dm, &dim);
1568:   VecGetLocalSize(vec, &nLocal);
1570:   switch (dim) {
1571:   case 1:
1572:     VecGetArray2dRead(vec, stag->nGhost[0], stag->entriesPerElement, stag->startGhost[0], 0, (PetscScalar ***)array);
1573:     break;
1574:   case 2:
1575:     VecGetArray3dRead(vec, stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar ****)array);
1576:     break;
1577:   case 3:
1578:     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);
1579:     break;
1580:   default:
1581:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
1582:   }
1583:   return 0;
1584: }

1586: /*@C
1587:   DMStagVecRestoreArray - restore access to a raw array

1589:   Logically Collective

1591:   Input Parameters:
1592: + dm - the `DMSTAG` object
1593: - vec - the `Vec` object

1595:   Output Parameter:
1596: . array - the array

1598:   Level: beginner

1600: .seealso: [](chapter_stag), `DMSTAG`, `DMStagVecGetArray()`, `DMDAVecRestoreArray()`, `DMDAVecRestoreArrayDOF()`
1601: @*/
1602: PetscErrorCode DMStagVecRestoreArray(DM dm, Vec vec, void *array)
1603: {
1604:   DM_Stag *const stag = (DM_Stag *)dm->data;
1605:   PetscInt       dim;
1606:   PetscInt       nLocal;

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

1629: /*@C
1630:   DMStagVecRestoreArrayRead - restore read-only access to a raw array

1632:   Logically Collective

1634:   Input Parameters:
1635: + dm - the `DMSTAG` object
1636: - vec - the Vec object

1638:   Output Parameter:
1639: . array - the read-only array

1641:   Level: beginner

1643: .seealso: [](chapter_stag), `DMSTAG`, `DMStagVecGetArrayRead()`, `DMDAVecRestoreArrayRead()`, `DMDAVecRestoreArrayDOFRead()`
1644: @*/
1645: PetscErrorCode DMStagVecRestoreArrayRead(DM dm, Vec vec, void *array)
1646: {
1647:   DM_Stag *const stag = (DM_Stag *)dm->data;
1648:   PetscInt       dim;
1649:   PetscInt       nLocal;

1653:   DMGetDimension(dm, &dim);
1654:   VecGetLocalSize(vec, &nLocal);
1656:   switch (dim) {
1657:   case 1:
1658:     VecRestoreArray2dRead(vec, stag->nGhost[0], stag->entriesPerElement, stag->startGhost[0], 0, (PetscScalar ***)array);
1659:     break;
1660:   case 2:
1661:     VecRestoreArray3dRead(vec, stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar ****)array);
1662:     break;
1663:   case 3:
1664:     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);
1665:     break;
1666:   default:
1667:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
1668:   }
1669:   return 0;
1670: }