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: }