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