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,boundaryTypeY,boundaryTypeZ - boundary types

 18:   Level: intermediate

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

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

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

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

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

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

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

 96:   Logically Collective

 98:   A high-level helper function to quickly extract local coordinate arrays.

100:   Note that 2-dimensional arrays are returned. See
101:   `DMStagVecGetArray()`, which is called internally to produce these arrays
102:   representing coordinates on elements and vertices (element boundaries)
103:   for a 1-dimensional DMStag in each coordinate direction.

105:   One should use `DMStagGetProductCoordinateLocationSlot()` to determine appropriate
106:   indices for the second dimension in these returned arrays. This function
107:   checks that the coordinate array is a suitable product of 1-dimensional
108:   DMStag objects.

110:   Input Parameter:
111: . dm - the DMStag object

113:   Output Parameters:
114: . arrX,arrY,arrZ - local 1D coordinate arrays

116:   Level: intermediate

118: .seealso: `DMSTAG`, `DMPRODUCT`, `DMStagGetProductCoordinateArraysRead()`, `DMStagSetUniformCoordinates()`, `DMStagSetUniformCoordinatesProduct()`, `DMStagGetProductCoordinateLocationSlot()`
119: @*/
120: PetscErrorCode DMStagGetProductCoordinateArrays(DM dm, void *arrX, void *arrY, void *arrZ)
121: {
122:   DMStagGetProductCoordinateArrays_Private(dm, arrX, arrY, arrZ, PETSC_FALSE);
123:   return 0;
124: }

126: /*@C
127:   DMStagGetProductCoordinateArraysRead - extract product coordinate arrays, read-only

129:   Logically Collective

131:   See the man page for `DMStagGetProductCoordinateArrays()` for more information.

133:   Input Parameter:
134: . dm - the DMStag object

136:   Output Parameters:
137: . arrX,arrY,arrZ - local 1D coordinate arrays

139:   Level: intermediate

141: .seealso: `DMSTAG`, `DMPRODUCT`, `DMStagGetProductCoordinateArrays()`, `DMStagSetUniformCoordinates()`, `DMStagSetUniformCoordinatesProduct()`, `DMStagGetProductCoordinateLocationSlot()`
142: @*/
143: PetscErrorCode DMStagGetProductCoordinateArraysRead(DM dm, void *arrX, void *arrY, void *arrZ)
144: {
145:   DMStagGetProductCoordinateArrays_Private(dm, arrX, arrY, arrZ, PETSC_TRUE);
146:   return 0;
147: }

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

152:   Not Collective

154:   High-level helper function to get slot indices for 1D coordinate DMs,
155:   for use with `DMStagGetProductCoordinateArrays()` and related functions.

157:   Input Parameters:
158: + dm - the DMStag object
159: - loc - the grid location

161:   Output Parameter:
162: . slot - the index to use in local arrays

164:   Notes:
165:   For `loc`, one should use `DMSTAG_LEFT`, `DMSTAG_ELEMENT`, or `DMSTAG_RIGHT` for "previous", "center" and "next"
166:   locations, respectively, in each dimension.
167:   One can equivalently use `DMSTAG_DOWN` or `DMSTAG_BACK` in place of `DMSTAG_LEFT`,
168:   and `DMSTAG_UP` or `DMSTACK_FRONT` in place of `DMSTAG_RIGHT`;

170:   This function checks that the coordinates are actually set up so that using the
171:   slots from any of the 1D coordinate sub-DMs are valid for all the 1D coordinate sub-DMs.

173:   Level: intermediate

175: .seealso: `DMSTAG`, `DMPRODUCT`, `DMStagGetProductCoordinateArrays()`, `DMStagGetProductCoordinateArraysRead()`, `DMStagSetUniformCoordinates()`
176: @*/
177: PETSC_EXTERN PetscErrorCode DMStagGetProductCoordinateLocationSlot(DM dm, DMStagStencilLocation loc, PetscInt *slot)
178: {
179:   DM       dmCoord;
180:   PetscInt dim, dofCheck[DMSTAG_MAX_STRATA], s, d;

183:   DMGetDimension(dm, &dim);
184:   DMGetCoordinateDM(dm, &dmCoord);
186:   {
187:     PetscBool isProduct;
188:     DMType    dmType;
189:     DMGetType(dmCoord, &dmType);
190:     PetscStrcmp(DMPRODUCT, dmType, &isProduct);
192:   }
193:   for (s = 0; s < DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
194:   for (d = 0; d < dim; ++d) {
195:     DM        subDM;
196:     DMType    dmType;
197:     PetscBool isStag;
198:     PetscInt  dof[DMSTAG_MAX_STRATA], subDim;
199:     DMProductGetDM(dmCoord, d, &subDM);
201:     DMGetDimension(subDM, &subDim);
203:     DMGetType(subDM, &dmType);
204:     PetscStrcmp(DMSTAG, dmType, &isStag);
206:     DMStagGetDOF(subDM, &dof[0], &dof[1], &dof[2], &dof[3]);
207:     if (d == 0) {
208:       const PetscInt component = 0;
209:       for (s = 0; s < DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
210:       DMStagGetLocationSlot(subDM, loc, component, slot);
211:     } else {
213:     }
214:   }
215:   return 0;
216: }

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

221:   Not Collective

223:   Input Parameter:
224: . dm - the DMStag object

226:   Output Parameters:
227: + x     - starting element index in first direction
228: . y     - starting element index in second direction
229: . z     - starting element index in third direction
230: . m     - element width in first direction
231: . n     - element width in second direction
232: . p     - element width in third direction
233: . nExtrax - number of extra partial elements in first direction
234: . nExtray - number of extra partial elements in second direction
235: - nExtraz - number of extra partial elements in third direction

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

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

244:   Level: beginner

246: .seealso: `DMSTAG`, `DMStagGetGhostCorners()`, `DMDAGetCorners()`
247: @*/
248: PetscErrorCode DMStagGetCorners(DM dm, PetscInt *x, PetscInt *y, PetscInt *z, PetscInt *m, PetscInt *n, PetscInt *p, PetscInt *nExtrax, PetscInt *nExtray, PetscInt *nExtraz)
249: {
250:   const DM_Stag *const stag = (DM_Stag *)dm->data;

253:   if (x) *x = stag->start[0];
254:   if (y) *y = stag->start[1];
255:   if (z) *z = stag->start[2];
256:   if (m) *m = stag->n[0];
257:   if (n) *n = stag->n[1];
258:   if (p) *p = stag->n[2];
259:   if (nExtrax) *nExtrax = stag->boundaryType[0] != DM_BOUNDARY_PERIODIC && stag->lastRank[0] ? 1 : 0;
260:   if (nExtray) *nExtray = stag->boundaryType[1] != DM_BOUNDARY_PERIODIC && stag->lastRank[1] ? 1 : 0;
261:   if (nExtraz) *nExtraz = stag->boundaryType[2] != DM_BOUNDARY_PERIODIC && stag->lastRank[2] ? 1 : 0;
262:   return 0;
263: }

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

268:   Not Collective

270:   Input Parameter:
271: . dm - the DMStag object

273:   Output Parameters:
274: + dof0 - the number of points per 0-cell (vertex/node)
275: . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
276: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
277: - dof3 - the number of points per 3-cell (element in 3D)

279:   Level: beginner

281: .seealso: `DMSTAG`, `DMStagGetCorners()`, `DMStagGetGhostCorners()`, `DMStagGetGlobalSizes()`, `DMStagGetStencilWidth()`, `DMStagGetBoundaryTypes()`, `DMStagGetLocationDOF()`, `DMDAGetDof()`
282: @*/
283: PetscErrorCode DMStagGetDOF(DM dm, PetscInt *dof0, PetscInt *dof1, PetscInt *dof2, PetscInt *dof3)
284: {
285:   const DM_Stag *const stag = (DM_Stag *)dm->data;

288:   if (dof0) *dof0 = stag->dof[0];
289:   if (dof1) *dof1 = stag->dof[1];
290:   if (dof2) *dof2 = stag->dof[2];
291:   if (dof3) *dof3 = stag->dof[3];
292:   return 0;
293: }

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

298:   Not Collective

300:   Input Parameter:
301: . dm - the DMStag object

303:   Output Parameters:
304: + x - the starting element index in the first direction
305: . y - the starting element index in the second direction
306: . z - the starting element index in the third direction
307: . m - the element width in the first direction
308: . n - the element width in the second direction
309: - p - the element width in the third direction

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

314:   Level: beginner

316: .seealso: `DMSTAG`, `DMStagGetCorners()`, `DMDAGetGhostCorners()`
317: @*/
318: PetscErrorCode DMStagGetGhostCorners(DM dm, PetscInt *x, PetscInt *y, PetscInt *z, PetscInt *m, PetscInt *n, PetscInt *p)
319: {
320:   const DM_Stag *const stag = (DM_Stag *)dm->data;

323:   if (x) *x = stag->startGhost[0];
324:   if (y) *y = stag->startGhost[1];
325:   if (z) *z = stag->startGhost[2];
326:   if (m) *m = stag->nGhost[0];
327:   if (n) *n = stag->nGhost[1];
328:   if (p) *p = stag->nGhost[2];
329:   return 0;
330: }

332: /*@C
333:   DMStagGetGlobalSizes - get global element counts

335:   Not Collective

337:   Input Parameter:
338: . dm - the DMStag object

340:   Output Parameters:
341: . M,N,P - global element counts in each direction

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

346:   Level: beginner

348: .seealso: `DMSTAG`, `DMStagGetLocalSizes()`, `DMDAGetInfo()`
349: @*/
350: PetscErrorCode DMStagGetGlobalSizes(DM dm, PetscInt *M, PetscInt *N, PetscInt *P)
351: {
352:   const DM_Stag *const stag = (DM_Stag *)dm->data;

355:   if (M) *M = stag->N[0];
356:   if (N) *N = stag->N[1];
357:   if (P) *P = stag->N[2];
358:   return 0;
359: }

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

364:   Not Collective

366:   Input Parameter:
367: . dm - the DMStag object

369:   Output Parameters:
370: . isFirstRank0,isFirstRank1,isFirstRank2 - whether this rank is first in each direction

372:   Level: intermediate

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

377: .seealso: `DMSTAG`, `DMStagGetIsLastRank()`
378: @*/
379: PetscErrorCode DMStagGetIsFirstRank(DM dm, PetscBool *isFirstRank0, PetscBool *isFirstRank1, PetscBool *isFirstRank2)
380: {
381:   const DM_Stag *const stag = (DM_Stag *)dm->data;

384:   if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
385:   if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
386:   if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
387:   return 0;
388: }

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

393:   Not Collective

395:   Input Parameter:
396: . dm - the DMStag object

398:   Output Parameters:
399: . isLastRank0,isLastRank1,isLastRank2 - whether this rank is last in each direction

401:   Level: intermediate

403:   Notes:
404:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to `NULL` in this case.
405:   Level: intermediate

407: .seealso: `DMSTAG`, `DMStagGetIsFirstRank()`
408: @*/
409: PetscErrorCode DMStagGetIsLastRank(DM dm, PetscBool *isLastRank0, PetscBool *isLastRank1, PetscBool *isLastRank2)
410: {
411:   const DM_Stag *const stag = (DM_Stag *)dm->data;

414:   if (isLastRank0) *isLastRank0 = stag->lastRank[0];
415:   if (isLastRank1) *isLastRank1 = stag->lastRank[1];
416:   if (isLastRank2) *isLastRank2 = stag->lastRank[2];
417:   return 0;
418: }

420: /*@C
421:   DMStagGetLocalSizes - get local elementwise sizes

423:   Not Collective

425:   Input Parameter:
426: . dm - the DMStag object

428:   Output Parameters:
429: . m,n,p - local element counts (excluding ghosts) in each direction

431:   Notes:
432:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to `NULL` in this case.
433:   Level: intermediate

435:   Level: beginner

437: .seealso: `DMSTAG`, `DMStagGetGlobalSizes()`, `DMStagGetDOF()`, `DMStagGetNumRanks()`, `DMDAGetLocalInfo()`
438: @*/
439: PetscErrorCode DMStagGetLocalSizes(DM dm, PetscInt *m, PetscInt *n, PetscInt *p)
440: {
441:   const DM_Stag *const stag = (DM_Stag *)dm->data;

444:   if (m) *m = stag->n[0];
445:   if (n) *n = stag->n[1];
446:   if (p) *p = stag->n[2];
447:   return 0;
448: }

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

453:   Not Collective

455:   Input Parameter:
456: . dm - the DMStag object

458:   Output Parameters:
459: . nRanks0,nRanks1,nRanks2 - number of ranks in each direction in the grid decomposition

461:   Notes:
462:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to `NULL` in this case.
463:   Level: intermediate

465:   Level: beginner

467: .seealso: `DMSTAG`, `DMStagGetGlobalSizes()`, `DMStagGetLocalSize()`, `DMStagSetNumRanks()`, `DMDAGetInfo()`
468: @*/
469: PetscErrorCode DMStagGetNumRanks(DM dm, PetscInt *nRanks0, PetscInt *nRanks1, PetscInt *nRanks2)
470: {
471:   const DM_Stag *const stag = (DM_Stag *)dm->data;

474:   if (nRanks0) *nRanks0 = stag->nRanks[0];
475:   if (nRanks1) *nRanks1 = stag->nRanks[1];
476:   if (nRanks2) *nRanks2 = stag->nRanks[2];
477:   return 0;
478: }

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

483:   Not Collective

485:   Input Parameter:
486: . dm - the DMStag object

488:   Output Parameters:
489: . entries - number of rank-native entries in the global representation

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

497:   Level: developer

499: .seealso: `DMSTAG`, `DMStagGetDOF()`, `DMStagGetEntriesLocal()`, `DMStagGetEntriesPerElement()`, `DMCreateLocalVector()`
500: @*/
501: PetscErrorCode DMStagGetEntries(DM dm, PetscInt *entries)
502: {
503:   const DM_Stag *const stag = (DM_Stag *)dm->data;

506:   if (entries) *entries = stag->entries;
507:   return 0;
508: }

510: /*@C
511:   DMStagGetEntriesLocal - get number of entries in the local representation

513:   Not Collective

515:   Input Parameter:
516: . dm - the DMStag object

518:   Output Parameters:
519: . entries - number of entries in the local representation

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

527:   Level: developer

529: .seealso: DMSTAG, DMStagGetDOF(), DMStagGetEntries(), DMStagGetEntriesPerElement(), DMCreateLocalVector()
530: @*/
531: PetscErrorCode DMStagGetEntriesLocal(DM dm, PetscInt *entries)
532: {
533:   const DM_Stag *const stag = (DM_Stag *)dm->data;

536:   if (entries) *entries = stag->entriesGhost;
537:   return 0;
538: }

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

543:   Not Collective

545:   Input Parameter:
546: . dm - the DMStag object

548:   Output Parameters:
549: . entriesPerElement - number of entries associated with each element in the local representation

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

555:   Level: developer

557: .seealso: `DMSTAG`, `DMStagGetDOF()`
558: @*/
559: PetscErrorCode DMStagGetEntriesPerElement(DM dm, PetscInt *entriesPerElement)
560: {
561:   const DM_Stag *const stag = (DM_Stag *)dm->data;

564:   if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
565:   return 0;
566: }

568: /*@C
569:   DMStagGetStencilType - get elementwise ghost/halo stencil type

571:   Not Collective

573:   Input Parameter:
574: . dm - the DMStag object

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

579:   Level: beginner

581: .seealso: `DMSTAG`, `DMStagSetStencilType()`, `DMStagGetStencilWidth`, `DMStagStencilType`
582: @*/
583: PetscErrorCode DMStagGetStencilType(DM dm, DMStagStencilType *stencilType)
584: {
585:   DM_Stag *const stag = (DM_Stag *)dm->data;

588:   *stencilType = stag->stencilType;
589:   return 0;
590: }

592: /*@C
593:   DMStagGetStencilWidth - get elementwise stencil width

595:   Not Collective

597:   Input Parameter:
598: . dm - the DMStag object

600:   Output Parameters:
601: . stencilWidth - stencil/halo/ghost width in elements

603:   Level: beginner

605: .seealso: `DMSTAG`, `DMStagSetStencilWidth()`, `DMStagGetStencilType()`, `DMDAGetStencilType()`
606: @*/
607: PetscErrorCode DMStagGetStencilWidth(DM dm, PetscInt *stencilWidth)
608: {
609:   const DM_Stag *const stag = (DM_Stag *)dm->data;

612:   if (stencilWidth) *stencilWidth = stag->stencilWidth;
613:   return 0;
614: }

616: /*@C
617:   DMStagGetOwnershipRanges - get elements per rank in each direction

619:   Not Collective

621:   Input Parameter:
622: .     dm - the DMStag object

624:   Output Parameters:
625: +     lx - ownership along x direction (optional)
626: .     ly - ownership along y direction (optional)
627: -     lz - ownership along z direction (optional)

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

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

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

637:   Level: intermediate

639: .seealso: `DMSTAG`, `DMStagSetGlobalSizes()`, `DMStagSetOwnershipRanges()`, `DMStagCreate1d()`, `DMStagCreate2d()`, `DMStagCreate3d()`, `DMDAGetOwnershipRanges()`
640: @*/
641: PetscErrorCode DMStagGetOwnershipRanges(DM dm, const PetscInt *lx[], const PetscInt *ly[], const PetscInt *lz[])
642: {
643:   const DM_Stag *const stag = (DM_Stag *)dm->data;

646:   if (lx) *lx = stag->l[0];
647:   if (ly) *ly = stag->l[1];
648:   if (lz) *lz = stag->l[2];
649:   return 0;
650: }

652: /*@C
653:   DMStagCreateCompatibleDMStag - create a compatible DMStag with different dof/stratum

655:   Collective

657:   Input Parameters:
658: + dm - the DMStag object
659: - dof0,dof1,dof2,dof3 - number of dof on each stratum in the new DMStag

661:   Output Parameters:
662: . newdm - the new, compatible DMStag

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

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

671:   Level: intermediate

673: .seealso: `DMSTAG`, `DMDACreateCompatibleDMDA()`, `DMGetCompatibility()`, `DMStagMigrateVec()`
674: @*/
675: PetscErrorCode DMStagCreateCompatibleDMStag(DM dm, PetscInt dof0, PetscInt dof1, PetscInt dof2, PetscInt dof3, DM *newdm)
676: {
678:   DMStagDuplicateWithoutSetup(dm, PetscObjectComm((PetscObject)dm), newdm);
679:   DMStagSetDOF(*newdm, dof0, dof1, dof2, dof3);
680:   DMSetUp(*newdm);
681:   return 0;
682: }

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

687:   Not Collective

689:   Input Parameters:
690: + dm - the DMStag object
691: . loc - location relative to an element
692: - c - component

694:   Output Parameter:
695: . slot - index to use

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

702:   Level: beginner

704: .seealso: `DMSTAG`, `DMStagVecGetArray()`, `DMStagVecGetArrayRead()`, `DMStagGetDOF()`, `DMStagGetEntriesPerElement()`
705: @*/
706: PetscErrorCode DMStagGetLocationSlot(DM dm, DMStagStencilLocation loc, PetscInt c, PetscInt *slot)
707: {
708:   DM_Stag *const stag = (DM_Stag *)dm->data;

711:   if (PetscDefined(USE_DEBUG)) {
712:     PetscInt dof;
713:     DMStagGetLocationDOF(dm, loc, &dof);
716:   }
717:   *slot = stag->locationOffsets[loc] + c;
718:   return 0;
719: }

721: /*@C
722:   DMStagMigrateVec - transfer a vector associated with a DMStag to a vector associated with a compatible DMStag

724:   Collective

726:   Input Parameters:
727: + dm - the source DMStag object
728: . vec - the source vector, compatible with `dm`
729: . dmTo - the compatible destination DMStag object
730: - vecTo - the destination vector, compatible with `dmTo`

732:   Notes:
733:   Extra dof are ignored, and unfilled dof are zeroed.
734:   Currently only implemented to migrate global vectors to global vectors.
735:   For the defintion of compatibility of `DM`s, see `DMGetCompatibility()`.

737:   Level: advanced

739: .seealso: `DMSTAG`, `DMStagCreateCompatibleDMStag()`, `DMGetCompatibility()`, `DMStagVecSplitToDMDA()`
740: @*/
741: PetscErrorCode DMStagMigrateVec(DM dm, Vec vec, DM dmTo, Vec vecTo)
742: {
743:   DM_Stag *const     stag   = (DM_Stag *)dm->data;
744:   DM_Stag *const     stagTo = (DM_Stag *)dmTo->data;
745:   PetscInt           nLocalTo, nLocal, dim, i, j, k;
746:   PetscInt           start[DMSTAG_MAX_DIM], startGhost[DMSTAG_MAX_DIM], n[DMSTAG_MAX_DIM], nExtra[DMSTAG_MAX_DIM], offset[DMSTAG_MAX_DIM];
747:   Vec                vecToLocal, vecLocal;
748:   PetscBool          compatible, compatibleSet;
749:   const PetscScalar *arr;
750:   PetscScalar       *arrTo;
751:   const PetscInt     epe   = stag->entriesPerElement;
752:   const PetscInt     epeTo = stagTo->entriesPerElement;

758:   DMGetCompatibility(dm, dmTo, &compatible, &compatibleSet);
760:   DMGetDimension(dm, &dim);
761:   VecGetLocalSize(vec, &nLocal);
762:   VecGetLocalSize(vecTo, &nLocalTo);
764:   DMStagGetCorners(dm, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], &nExtra[0], &nExtra[1], &nExtra[2]);
765:   DMStagGetGhostCorners(dm, &startGhost[0], &startGhost[1], &startGhost[2], NULL, NULL, NULL);
766:   for (i = 0; i < DMSTAG_MAX_DIM; ++i) offset[i] = start[i] - startGhost[i];

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

840: /*@C
841:   DMStagPopulateLocalToGlobalInjective - populate an internal 1-to-1 local-to-global map

843:   Collective

845:   Creates an internal object which explicitly maps a single local degree of
846:   freedom to each global degree of freedom. This is used, if populated,
847:   instead of SCATTER_REVERSE_LOCAL with the (1-to-many, in general)
848:   global-to-local map, when DMLocalToGlobal() is called with INSERT_VALUES.
849:   This allows usage, for example, even in the periodic, 1-rank case, where
850:   the inverse of the global-to-local map, even when restricted to on-rank
851:   communication, is non-injective. This is at the cost of storing an additional
852:   VecScatter object inside each DMStag object.

854:   Input Parameter:
855: . dm - the DMStag object

857:   Notes:
858:   In normal usage, library users shouldn't be concerned with this function,
859:   as it is called during DMSetUp(), when required.

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

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

868:   Level: developer

870: .seealso: `DMSTAG`, `DMLocalToGlobal()`, `VecScatter`
871: @*/
872: PetscErrorCode DMStagPopulateLocalToGlobalInjective(DM dm)
873: {
874:   PetscInt       dim;
875:   DM_Stag *const stag = (DM_Stag *)dm->data;

878:   if (stag->ltog_injective) return 0; /* Don't re-populate */
879:   DMGetDimension(dm, &dim);
880:   switch (dim) {
881:   case 1:
882:     DMStagPopulateLocalToGlobalInjective_1d(dm);
883:     break;
884:   case 2:
885:     DMStagPopulateLocalToGlobalInjective_2d(dm);
886:     break;
887:   case 3:
888:     DMStagPopulateLocalToGlobalInjective_3d(dm);
889:     break;
890:   default:
891:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
892:   }
893:   return 0;
894: }

896: static PetscErrorCode DMStagRestoreProductCoordinateArrays_Private(DM dm, void *arrX, void *arrY, void *arrZ, PetscBool read)
897: {
898:   PetscInt dim, d;
899:   void    *arr[DMSTAG_MAX_DIM];
900:   DM       dmCoord;

903:   DMGetDimension(dm, &dim);
905:   arr[0] = arrX;
906:   arr[1] = arrY;
907:   arr[2] = arrZ;
908:   DMGetCoordinateDM(dm, &dmCoord);
909:   for (d = 0; d < dim; ++d) {
910:     DM  subDM;
911:     Vec coord1d_local;

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

916:     DMProductGetDM(dmCoord, d, &subDM);
917:     DMGetCoordinatesLocal(subDM, &coord1d_local);
918:     if (read) {
919:       DMStagVecRestoreArrayRead(subDM, coord1d_local, arr[d]);
920:     } else {
921:       DMStagVecRestoreArray(subDM, coord1d_local, arr[d]);
922:     }
923:   }
924:   return 0;
925: }

927: /*@C
928:   DMStagRestoreProductCoordinateArrays - restore local array access

930:   Logically Collective

932:   Input Parameter:
933: . dm - the DMStag object

935:   Output Parameters:
936: . arrX,arrY,arrZ - local 1D coordinate arrays

938:   Level: intermediate

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

943:   ```
944:   DMGetCoordinateDM(dm,&cdm);
945:   for (d=0; d<3; ++d) {
946:     DM  subdm;
947:     Vec coor,coor_local;

949:     DMProductGetDM(cdm,d,&subdm);
950:     DMGetCoordinates(subdm,&coor);
951:     DMGetCoordinatesLocal(subdm,&coor_local);
952:     DMLocalToGlobal(subdm,coor_local,INSERT_VALUES,coor);
953:     PetscPrintf(PETSC_COMM_WORLD,"Coordinates dim %" PetscInt_FMT ":\n",d);
954:     VecView(coor,PETSC_VIEWER_STDOUT_WORLD);
955:   }
956:    ```

958: .seealso: `DMSTAG`, `DMStagGetProductCoordinateArrays()`, `DMStagGetProductCoordinateArraysRead()`
959: @*/
960: PetscErrorCode DMStagRestoreProductCoordinateArrays(DM dm, void *arrX, void *arrY, void *arrZ)
961: {
962:   DMStagRestoreProductCoordinateArrays_Private(dm, arrX, arrY, arrZ, PETSC_FALSE);
963:   return 0;
964: }

966: /*@C
967:   DMStagRestoreProductCoordinateArraysRead - restore local product array access, read-only

969:   Logically Collective

971:   Input Parameter:
972: . dm - the DMStag object

974:   Output Parameters:
975: . arrX,arrY,arrZ - local 1D coordinate arrays

977:   Level: intermediate

979: .seealso: `DMSTAG`, `DMStagGetProductCoordinateArrays()`, `DMStagGetProductCoordinateArraysRead()`
980: @*/
981: PetscErrorCode DMStagRestoreProductCoordinateArraysRead(DM dm, void *arrX, void *arrY, void *arrZ)
982: {
983:   DMStagRestoreProductCoordinateArrays_Private(dm, arrX, arrY, arrZ, PETSC_TRUE);
984:   return 0;
985: }

987: /*@C
988:   DMStagSetBoundaryTypes - set DMStag boundary types

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

992:   Input Parameters:
993: + dm - the DMStag object
994: - boundaryType0,boundaryType1,boundaryType2 - boundary types in each direction

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

999:   Level: advanced

1001: .seealso: `DMSTAG`, `DMBoundaryType`, `DMStagCreate1d()`, `DMStagCreate2d()`, `DMStagCreate3d()`, `DMDASetBoundaryType()`
1002: @*/
1003: PetscErrorCode DMStagSetBoundaryTypes(DM dm, DMBoundaryType boundaryType0, DMBoundaryType boundaryType1, DMBoundaryType boundaryType2)
1004: {
1005:   DM_Stag *const stag = (DM_Stag *)dm->data;
1006:   PetscInt       dim;

1008:   DMGetDimension(dm, &dim);
1014:   stag->boundaryType[0] = boundaryType0;
1015:   if (dim > 1) stag->boundaryType[1] = boundaryType1;
1016:   if (dim > 2) stag->boundaryType[2] = boundaryType2;
1017:   return 0;
1018: }

1020: /*@C
1021:   DMStagSetCoordinateDMType - set DM type to store coordinates

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

1025:   Input Parameters:
1026: + dm - the DMStag object
1027: - dmtype - DMtype for coordinates, either `DMSTAG` or `DMPRODUCT`

1029:   Level: advanced

1031: .seealso: `DMSTAG`, `DMPRODUCT`, `DMGetCoordinateDM()`, `DMStagSetUniformCoordinates()`, `DMStagSetUniformCoordinatesExplicit()`, `DMStagSetUniformCoordinatesProduct()`, `DMType`
1032: @*/
1033: PetscErrorCode DMStagSetCoordinateDMType(DM dm, DMType dmtype)
1034: {
1035:   DM_Stag *const stag = (DM_Stag *)dm->data;

1038:   PetscFree(stag->coordinateDMType);
1039:   PetscStrallocpy(dmtype, (char **)&stag->coordinateDMType);
1040:   return 0;
1041: }

1043: /*@C
1044:   DMStagSetDOF - set dof/stratum

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

1048:   Input Parameters:
1049: + dm - the DMStag object
1050: - dof0,dof1,dof2,dof3 - dof per stratum

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

1055:   Level: advanced

1057: .seealso: `DMSTAG`, `DMDASetDof()`
1058: @*/
1059: PetscErrorCode DMStagSetDOF(DM dm, PetscInt dof0, PetscInt dof1, PetscInt dof2, PetscInt dof3)
1060: {
1061:   DM_Stag *const stag = (DM_Stag *)dm->data;
1062:   PetscInt       dim;

1070:   DMGetDimension(dm, &dim);
1075:   stag->dof[0] = dof0;
1076:   stag->dof[1] = dof1;
1077:   if (dim > 1) stag->dof[2] = dof2;
1078:   if (dim > 2) stag->dof[3] = dof3;
1079:   return 0;
1080: }

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

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

1087:   Input Parameters:
1088: + dm - the DMStag object
1089: - nRanks0,nRanks1,nRanks2 - number of ranks in each direction

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

1094:   Level: developer

1096: .seealso: `DMSTAG`, `DMDASetNumProcs()`
1097: @*/
1098: PetscErrorCode DMStagSetNumRanks(DM dm, PetscInt nRanks0, PetscInt nRanks1, PetscInt nRanks2)
1099: {
1100:   DM_Stag *const stag = (DM_Stag *)dm->data;
1101:   PetscInt       dim;

1108:   DMGetDimension(dm, &dim);
1112:   if (nRanks0) stag->nRanks[0] = nRanks0;
1113:   if (dim > 1 && nRanks1) stag->nRanks[1] = nRanks1;
1114:   if (dim > 2 && nRanks2) stag->nRanks[2] = nRanks2;
1115:   return 0;
1116: }

1118: /*@C
1119:   DMStagSetStencilType - set elementwise ghost/halo stencil type

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

1123:   Input Parameters:
1124: + dm - the DMStag object
1125: - stencilType - the elementwise ghost stencil type: `DMSTAG_STENCIL_BOX`, `DMSTAG_STENCIL_STAR`, or `DMSTAG_STENCIL_NONE`

1127:   Level: beginner

1129: .seealso: `DMSTAG`, `DMStagGetStencilType()`, `DMStagSetStencilWidth()`, `DMStagStencilType`
1130: @*/
1131: PetscErrorCode DMStagSetStencilType(DM dm, DMStagStencilType stencilType)
1132: {
1133:   DM_Stag *const stag = (DM_Stag *)dm->data;

1138:   stag->stencilType = stencilType;
1139:   return 0;
1140: }

1142: /*@C
1143:   DMStagSetStencilWidth - set elementwise stencil width

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

1147:   Input Parameters:
1148: + dm - the DMStag object
1149: - stencilWidth - stencil/halo/ghost width in elements

1151:   Notes:
1152:   The width value is not used when `DMSTAG_STENCIL_NONE` is specified.

1154:   Level: beginner

1156: .seealso: `DMSTAG`, `DMStagGetStencilWidth()`, `DMStagGetStencilType()`, `DMStagStencilType`
1157: @*/
1158: PetscErrorCode DMStagSetStencilWidth(DM dm, PetscInt stencilWidth)
1159: {
1160:   DM_Stag *const stag = (DM_Stag *)dm->data;

1166:   stag->stencilWidth = stencilWidth;
1167:   return 0;
1168: }

1170: /*@C
1171:   DMStagSetGlobalSizes - set global element counts in each direction

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

1175:   Input Parameters:
1176: + dm - the DMStag object
1177: - N0,N1,N2 - global elementwise sizes

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

1182:   Level: advanced

1184: .seealso: `DMSTAG`, `DMStagGetGlobalSizes()`, `DMDASetSizes()`
1185: @*/
1186: PetscErrorCode DMStagSetGlobalSizes(DM dm, PetscInt N0, PetscInt N1, PetscInt N2)
1187: {
1188:   DM_Stag *const stag = (DM_Stag *)dm->data;
1189:   PetscInt       dim;

1193:   DMGetDimension(dm, &dim);
1197:   if (N0) stag->N[0] = N0;
1198:   if (N1) stag->N[1] = N1;
1199:   if (N2) stag->N[2] = N2;
1200:   return 0;
1201: }

1203: /*@C
1204:   DMStagSetOwnershipRanges - set elements per rank in each direction

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

1208:   Input Parameters:
1209: + dm - the DMStag object
1210: - lx,ly,lz - element counts for each rank in each direction

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

1215:   Level: developer

1217: .seealso: `DMSTAG`, `DMStagSetGlobalSizes()`, `DMStagGetOwnershipRanges()`, `DMDASetOwnershipRanges()`
1218: @*/
1219: PetscErrorCode DMStagSetOwnershipRanges(DM dm, PetscInt const *lx, PetscInt const *ly, PetscInt const *lz)
1220: {
1221:   DM_Stag *const  stag = (DM_Stag *)dm->data;
1222:   const PetscInt *lin[3];
1223:   PetscInt        d, dim;

1227:   lin[0] = lx;
1228:   lin[1] = ly;
1229:   lin[2] = lz;
1230:   DMGetDimension(dm, &dim);
1231:   for (d = 0; d < dim; ++d) {
1232:     if (lin[d]) {
1234:       if (!stag->l[d]) PetscMalloc1(stag->nRanks[d], &stag->l[d]);
1235:       PetscArraycpy(stag->l[d], lin[d], stag->nRanks[d]);
1236:     }
1237:   }
1238:   return 0;
1239: }

1241: /*@C
1242:   DMStagSetUniformCoordinates - set DMStag coordinates to be a uniform grid

1244:   Collective

1246:   Input Parameters:
1247: + dm - the DMStag object
1248: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

1250:   Notes:
1251:   DMStag supports 2 different types of coordinate DM: `DMSTAG` and `DMPRODUCT`.
1252:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

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

1262:   Level: advanced

1264: .seealso: `DMSTAG`, `DMPRODUCT`, `DMStagSetUniformCoordinatesExplicit()`, `DMStagSetUniformCoordinatesProduct()`, `DMStagSetCoordinateDMType()`, `DMGetCoordinateDM()`, `DMGetCoordinates()`, `DMDASetUniformCoordinates()`, `DMBoundaryType`
1265: @*/
1266: PetscErrorCode DMStagSetUniformCoordinates(DM dm, PetscReal xmin, PetscReal xmax, PetscReal ymin, PetscReal ymax, PetscReal zmin, PetscReal zmax)
1267: {
1268:   DM_Stag *const stag = (DM_Stag *)dm->data;
1269:   PetscBool      flg_stag, flg_product;

1274:   PetscStrcmp(stag->coordinateDMType, DMSTAG, &flg_stag);
1275:   PetscStrcmp(stag->coordinateDMType, DMPRODUCT, &flg_product);
1276:   if (flg_stag) {
1277:     DMStagSetUniformCoordinatesExplicit(dm, xmin, xmax, ymin, ymax, zmin, zmax);
1278:   } else if (flg_product) {
1279:     DMStagSetUniformCoordinatesProduct(dm, xmin, xmax, ymin, ymax, zmin, zmax);
1280:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported DM Type %s", stag->coordinateDMType);
1281:   return 0;
1282: }

1284: /*@C
1285:   DMStagSetUniformCoordinatesExplicit - set DMStag coordinates to be a uniform grid, storing all values

1287:   Collective

1289:   Input Parameters:
1290: + dm - the DMStag object
1291: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

1293:   Notes:
1294:   DMStag supports 2 different types of coordinate DM: either another DMStag, or a DMProduct.
1295:   If the grid is orthogonal, using DMProduct should be more efficient.

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

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

1302:   Level: beginner

1304: .seealso: `DMSTAG`, `DMStagSetUniformCoordinates()`, `DMStagSetUniformCoordinatesProduct()`, `DMStagSetCoordinateDMType()`
1305: @*/
1306: PetscErrorCode DMStagSetUniformCoordinatesExplicit(DM dm, PetscReal xmin, PetscReal xmax, PetscReal ymin, PetscReal ymax, PetscReal zmin, PetscReal zmax)
1307: {
1308:   DM_Stag *const stag = (DM_Stag *)dm->data;
1309:   PetscInt       dim;
1310:   PetscBool      flg;

1314:   PetscStrcmp(stag->coordinateDMType, DMSTAG, &flg);
1316:   DMStagSetCoordinateDMType(dm, DMSTAG);
1317:   DMGetDimension(dm, &dim);
1318:   switch (dim) {
1319:   case 1:
1320:     DMStagSetUniformCoordinatesExplicit_1d(dm, xmin, xmax);
1321:     break;
1322:   case 2:
1323:     DMStagSetUniformCoordinatesExplicit_2d(dm, xmin, xmax, ymin, ymax);
1324:     break;
1325:   case 3:
1326:     DMStagSetUniformCoordinatesExplicit_3d(dm, xmin, xmax, ymin, ymax, zmin, zmax);
1327:     break;
1328:   default:
1329:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
1330:   }
1331:   DMCoarsenHookRemove(dm, DMRestrictHook_Coordinates, NULL, NULL);
1332:   return 0;
1333: }

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

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

1340:   Collective

1342:   Input Parameters:
1343: + dm - the DMStag object
1344: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

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

1349:   The per-dimension 1-dimensional DMStag objects that comprise the product
1350:   always have active 0-cells (vertices, element boundaries) and 1-cells
1351:   (element centers).

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

1356:   Level: intermediate

1358: .seealso: `DMSTAG`, `DMPRODUCT`, `DMStagSetUniformCoordinates()`, `DMStagSetUniformCoordinatesExplicit()`, `DMStagSetCoordinateDMType()`
1359: @*/
1360: PetscErrorCode DMStagSetUniformCoordinatesProduct(DM dm, PetscReal xmin, PetscReal xmax, PetscReal ymin, PetscReal ymax, PetscReal zmin, PetscReal zmax)
1361: {
1362:   DM_Stag *const stag = (DM_Stag *)dm->data;
1363:   DM             dmc;
1364:   PetscInt       dim, d, dof0, dof1;
1365:   PetscBool      flg;

1369:   PetscStrcmp(stag->coordinateDMType, DMPRODUCT, &flg);
1371:   DMStagSetCoordinateDMType(dm, DMPRODUCT);
1372:   DMGetCoordinateDM(dm, &dmc);
1373:   DMGetDimension(dm, &dim);

1375:   /* Create 1D sub-DMs, living on subcommunicators.
1376:      Always include both vertex and element dof, regardless of the active strata of the DMStag */
1377:   dof0 = 1;
1378:   dof1 = 1;

1380:   for (d = 0; d < dim; ++d) {
1381:     DM                subdm;
1382:     MPI_Comm          subcomm;
1383:     PetscMPIInt       color;
1384:     const PetscMPIInt key = 0; /* let existing rank break ties */

1386:     /* Choose colors based on position in the plane orthogonal to this dim, and split */
1387:     switch (d) {
1388:     case 0:
1389:       color = (dim > 1 ? stag->rank[1] : 0) + (dim > 2 ? stag->nRanks[1] * stag->rank[2] : 0);
1390:       break;
1391:     case 1:
1392:       color = stag->rank[0] + (dim > 2 ? stag->nRanks[0] * stag->rank[2] : 0);
1393:       break;
1394:     case 2:
1395:       color = stag->rank[0] + stag->nRanks[0] * stag->rank[1];
1396:       break;
1397:     default:
1398:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension index %" PetscInt_FMT, d);
1399:     }
1400:     MPI_Comm_split(PetscObjectComm((PetscObject)dm), color, key, &subcomm);

1402:     /* Create sub-DMs living on these new communicators (which are destroyed by DMProduct) */
1403:     DMStagCreate1d(subcomm, stag->boundaryType[d], stag->N[d], dof0, dof1, stag->stencilType, stag->stencilWidth, stag->l[d], &subdm);
1404:     /* Copy vector and matrix type information from parent DM */
1405:     DMSetVecType(subdm, dm->vectype);
1406:     DMSetMatType(subdm, dm->mattype);
1407:     DMSetUp(subdm);
1408:     switch (d) {
1409:     case 0:
1410:       DMStagSetUniformCoordinatesExplicit(subdm, xmin, xmax, 0, 0, 0, 0);
1411:       break;
1412:     case 1:
1413:       DMStagSetUniformCoordinatesExplicit(subdm, ymin, ymax, 0, 0, 0, 0);
1414:       break;
1415:     case 2:
1416:       DMStagSetUniformCoordinatesExplicit(subdm, zmin, zmax, 0, 0, 0, 0);
1417:       break;
1418:     default:
1419:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension index %" PetscInt_FMT, d);
1420:     }
1421:     DMProductSetDM(dmc, d, subdm);
1422:     DMProductSetDimensionIndex(dmc, d, 0);
1423:     DMDestroy(&subdm);
1424:     MPI_Comm_free(&subcomm);
1425:   }
1426:   DMCoarsenHookRemove(dm, DMRestrictHook_Coordinates, NULL, NULL);
1427:   return 0;
1428: }

1430: /*@C
1431:   DMStagVecGetArray - get access to local array

1433:   Logically Collective

1435:   This function returns a (dim+1)-dimensional array for a dim-dimensional
1436:   DMStag.

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

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

1445:   For example, for a 3D DMStag, indexing is `array[k][j][i][idx]`, where `k` is the
1446:   index in the z-direction, `j` is the index in the y-direction, and `i` is the
1447:   index in the x-direction.

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

1453:   Input Parameters:
1454: + dm - the DMStag object
1455: - vec - the `Vec` object

1457:   Output Parameters:
1458: . array - the array

1460:   Notes:
1461:   `DMStagVecRestoreArray()` must be called, once finished with the array

1463:   Level: beginner

1465: .seealso: `DMSTAG`, `DMStagVecGetArrayRead()`, `DMStagGetLocationSlot()`, `DMGetLocalVector()`, `DMCreateLocalVector()`, `DMGetGlobalVector()`, `DMCreateGlobalVector()`, `DMDAVecGetArray()`, `DMDAVecGetArrayDOF()`
1466: @*/
1467: PetscErrorCode DMStagVecGetArray(DM dm, Vec vec, void *array)
1468: {
1469:   DM_Stag *const stag = (DM_Stag *)dm->data;
1470:   PetscInt       dim;
1471:   PetscInt       nLocal;

1475:   DMGetDimension(dm, &dim);
1476:   VecGetLocalSize(vec, &nLocal);
1478:   switch (dim) {
1479:   case 1:
1480:     VecGetArray2d(vec, stag->nGhost[0], stag->entriesPerElement, stag->startGhost[0], 0, (PetscScalar ***)array);
1481:     break;
1482:   case 2:
1483:     VecGetArray3d(vec, stag->nGhost[1], stag->nGhost[0], stag->entriesPerElement, stag->startGhost[1], stag->startGhost[0], 0, (PetscScalar ****)array);
1484:     break;
1485:   case 3:
1486:     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);
1487:     break;
1488:   default:
1489:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
1490:   }
1491:   return 0;
1492: }

1494: /*@C
1495:   DMStagVecGetArrayRead - get read-only access to a local array

1497:   Logically Collective

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

1501:   Input Parameters:
1502: + dm - the DMStag object
1503: - vec - the `Vec` object

1505:   Output Parameters:
1506: . array - the read-only array

1508:   Notes:
1509:   `DMStagVecRestoreArrayRead()` must be called, once finished with the array

1511:   Level: beginner

1513: .seealso: `DMSTAG`, `DMStagVecGetArrayRead()`, `DMStagGetLocationSlot()`, `DMGetLocalVector()`, `DMCreateLocalVector()`, `DMGetGlobalVector()`, `DMCreateGlobalVector()`, `DMDAVecGetArrayRead()`, `DMDAVecGetArrayDOFRead()`
1514: @*/
1515: PetscErrorCode DMStagVecGetArrayRead(DM dm, Vec vec, void *array)
1516: {
1517:   DM_Stag *const stag = (DM_Stag *)dm->data;
1518:   PetscInt       dim;
1519:   PetscInt       nLocal;

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

1542: /*@C
1543:   DMStagVecRestoreArray - restore access to a raw array

1545:   Logically Collective

1547:   Input Parameters:
1548: + dm - the DMStag object
1549: - vec - the `Vec` object

1551:   Output Parameters:
1552: . array - the array

1554:   Level: beginner

1556: .seealso: `DMSTAG`, `DMStagVecGetArray()`, `DMDAVecRestoreArray()`, `DMDAVecRestoreArrayDOF()`
1557: @*/
1558: PetscErrorCode DMStagVecRestoreArray(DM dm, Vec vec, void *array)
1559: {
1560:   DM_Stag *const stag = (DM_Stag *)dm->data;
1561:   PetscInt       dim;
1562:   PetscInt       nLocal;

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

1585: /*@C
1586:   DMStagVecRestoreArrayRead - restore read-only access to a raw array

1588:   Logically Collective

1590:   Input Parameters:
1591: + dm - the DMStag object
1592: - vec - the Vec object

1594:   Output Parameters:
1595: . array - the read-only array

1597:   Level: beginner

1599: .seealso: `DMSTAG`, `DMStagVecGetArrayRead()`, `DMDAVecRestoreArrayRead()`, `DMDAVecRestoreArrayDOFRead()`
1600: @*/
1601: PetscErrorCode DMStagVecRestoreArrayRead(DM dm, Vec vec, void *array)
1602: {
1603:   DM_Stag *const stag = (DM_Stag *)dm->data;
1604:   PetscInt       dim;
1605:   PetscInt       nLocal;

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