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>
4: /*@C
5: DMStagGetBoundaryTypes - get boundary types
7: Not Collective
9: Input Parameter:
10: . dm - the DMStag object
12: Output Parameters:
13: . boundaryTypeX,boundaryTypeY,boundaryTypeZ - boundary types
15: Level: intermediate
17: .seealso: DMSTAG
18: @*/
19: PetscErrorCode DMStagGetBoundaryTypes(DM dm,DMBoundaryType *boundaryTypeX,DMBoundaryType *boundaryTypeY,DMBoundaryType *boundaryTypeZ)
20: {
21: const DM_Stag * const stag = (DM_Stag*)dm->data;
22: PetscInt dim;
25: DMGetDimension(dm,&dim);
26: if (boundaryTypeX) *boundaryTypeX = stag->boundaryType[0];
27: if (boundaryTypeY && dim > 1) *boundaryTypeY = stag->boundaryType[1];
28: if (boundaryTypeZ && dim > 2) *boundaryTypeZ = stag->boundaryType[2];
29: return 0;
30: }
32: static PetscErrorCode DMStagGetProductCoordinateArrays_Private(DM dm,void* arrX,void* arrY,void* arrZ,PetscBool read)
33: {
34: PetscInt dim,d,dofCheck[DMSTAG_MAX_STRATA],s;
35: DM dmCoord;
36: void* arr[DMSTAG_MAX_DIM];
37: PetscBool checkDof;
40: DMGetDimension(dm,&dim);
42: arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
43: DMGetCoordinateDM(dm,&dmCoord);
45: {
46: PetscBool isProduct;
47: DMType dmType;
48: DMGetType(dmCoord,&dmType);
49: PetscStrcmp(DMPRODUCT,dmType,&isProduct);
51: }
52: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
53: checkDof = PETSC_FALSE;
54: for (d=0; d<dim; ++d) {
55: DM subDM;
56: DMType dmType;
57: PetscBool isStag;
58: PetscInt dof[DMSTAG_MAX_STRATA],subDim;
59: Vec coord1d_local;
61: /* Ignore unrequested arrays */
62: if (!arr[d]) continue;
64: DMProductGetDM(dmCoord,d,&subDM);
66: DMGetDimension(subDM,&subDim);
68: DMGetType(subDM,&dmType);
69: PetscStrcmp(DMSTAG,dmType,&isStag);
71: DMStagGetDOF(subDM,&dof[0],&dof[1],&dof[2],&dof[3]);
72: if (!checkDof) {
73: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
74: checkDof = PETSC_TRUE;
75: } else {
76: for (s=0; s<DMSTAG_MAX_STRATA; ++s) {
78: }
79: }
80: DMGetCoordinatesLocal(subDM,&coord1d_local);
81: if (read) {
82: DMStagVecGetArrayRead(subDM,coord1d_local,arr[d]);
83: } else {
84: DMStagVecGetArray(subDM,coord1d_local,arr[d]);
85: }
86: }
87: return 0;
88: }
90: /*@C
91: DMStagGetProductCoordinateArrays - extract local product coordinate arrays, one per dimension
93: Logically Collective
95: A high-level helper function to quickly extract local coordinate arrays.
97: Note that 2-dimensional arrays are returned. See
98: DMStagVecGetArray(), which is called internally to produce these arrays
99: representing coordinates on elements and vertices (element boundaries)
100: for a 1-dimensional DMStag in each coordinate direction.
102: One should use DMStagGetProductCoordinateLocationSlot() to determine appropriate
103: indices for the second dimension in these returned arrays. This function
104: checks that the coordinate array is a suitable product of 1-dimensional
105: DMStag objects.
107: Input Parameter:
108: . dm - the DMStag object
110: Output Parameters:
111: . arrX,arrY,arrZ - local 1D coordinate arrays
113: Level: intermediate
115: .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArraysRead(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagGetProductCoordinateLocationSlot()
116: @*/
117: PetscErrorCode DMStagGetProductCoordinateArrays(DM dm,void* arrX,void* arrY,void* arrZ)
118: {
119: DMStagGetProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_FALSE);
120: return 0;
121: }
123: /*@C
124: DMStagGetProductCoordinateArraysRead - extract product coordinate arrays, read-only
126: Logically Collective
128: See the man page for DMStagGetProductCoordinateArrays() for more information.
130: Input Parameter:
131: . dm - the DMStag object
133: Output Parameters:
134: . arrX,arrY,arrZ - local 1D coordinate arrays
136: Level: intermediate
138: .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArrays(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagGetProductCoordinateLocationSlot()
139: @*/
140: PetscErrorCode DMStagGetProductCoordinateArraysRead(DM dm,void* arrX,void* arrY,void* arrZ)
141: {
142: DMStagGetProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_TRUE);
143: return 0;
144: }
146: /*@C
147: DMStagGetProductCoordinateLocationSlot - get slot for use with local product coordinate arrays
149: Not Collective
151: High-level helper function to get slot indices for 1D coordinate DMs,
152: for use with DMStagGetProductCoordinateArrays() and related functions.
154: Input Parameters:
155: + dm - the DMStag object
156: - loc - the grid location
158: Output Parameter:
159: . slot - the index to use in local arrays
161: Notes:
162: For loc, one should use DMSTAG_LEFT, DMSTAG_ELEMENT, or DMSTAG_RIGHT for "previous", "center" and "next"
163: locations, respectively, in each dimension.
164: One can equivalently use DMSTAG_DOWN or DMSTAG_BACK in place of DMSTAG_LEFT,
165: and DMSTAG_UP or DMSTACK_FRONT in place of DMSTAG_RIGHT;
167: This function checks that the coordinates are actually set up so that using the
168: slots from any of the 1D coordinate sub-DMs are valid for all the 1D coordinate sub-DMs.
170: Level: intermediate
172: .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead(), DMStagSetUniformCoordinates()
173: @*/
174: PETSC_EXTERN PetscErrorCode DMStagGetProductCoordinateLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt *slot)
175: {
176: DM dmCoord;
177: PetscInt dim,dofCheck[DMSTAG_MAX_STRATA],s,d;
180: DMGetDimension(dm,&dim);
181: DMGetCoordinateDM(dm,&dmCoord);
183: {
184: PetscBool isProduct;
185: DMType dmType;
186: DMGetType(dmCoord,&dmType);
187: PetscStrcmp(DMPRODUCT,dmType,&isProduct);
189: }
190: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
191: for (d=0; d<dim; ++d) {
192: DM subDM;
193: DMType dmType;
194: PetscBool isStag;
195: PetscInt dof[DMSTAG_MAX_STRATA],subDim;
196: DMProductGetDM(dmCoord,d,&subDM);
198: DMGetDimension(subDM,&subDim);
200: DMGetType(subDM,&dmType);
201: PetscStrcmp(DMSTAG,dmType,&isStag);
203: DMStagGetDOF(subDM,&dof[0],&dof[1],&dof[2],&dof[3]);
204: if (d == 0) {
205: const PetscInt component = 0;
206: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
207: DMStagGetLocationSlot(subDM,loc,component,slot);
208: } else {
209: for (s=0; s<DMSTAG_MAX_STRATA; ++s) {
211: }
212: }
213: }
214: return 0;
215: }
217: /*@C
218: DMStagGetCorners - return global element indices of the local region (excluding ghost points)
220: Not Collective
222: Input Parameter:
223: . dm - the DMStag object
225: Output Parameters:
226: + x - starting element index in first direction
227: . y - starting element index in second direction
228: . z - starting element index in third direction
229: . m - element width in first direction
230: . n - element width in second direction
231: . p - element width in third direction
232: . nExtrax - number of extra partial elements in first direction
233: . nExtray - number of extra partial elements in second direction
234: - nExtraz - number of extra partial elements in third direction
236: Notes:
237: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
239: The number of extra partial elements is either 1 or 0.
240: The value is 1 on right, top, and front non-periodic domain ("physical") boundaries,
241: in the x, y, and z directions respectively, and otherwise 0.
243: Level: beginner
245: .seealso: DMSTAG, DMStagGetGhostCorners(), DMDAGetCorners()
246: @*/
247: PetscErrorCode DMStagGetCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *nExtrax,PetscInt *nExtray,PetscInt *nExtraz)
248: {
249: const DM_Stag * const stag = (DM_Stag*)dm->data;
252: if (x) *x = stag->start[0];
253: if (y) *y = stag->start[1];
254: if (z) *z = stag->start[2];
255: if (m) *m = stag->n[0];
256: if (n) *n = stag->n[1];
257: if (p) *p = stag->n[2];
258: if (nExtrax) *nExtrax = stag->boundaryType[0] != DM_BOUNDARY_PERIODIC && stag->lastRank[0] ? 1 : 0;
259: if (nExtray) *nExtray = stag->boundaryType[1] != DM_BOUNDARY_PERIODIC && stag->lastRank[1] ? 1 : 0;
260: if (nExtraz) *nExtraz = stag->boundaryType[2] != DM_BOUNDARY_PERIODIC && stag->lastRank[2] ? 1 : 0;
261: return 0;
262: }
264: /*@C
265: DMStagGetDOF - get number of DOF associated with each stratum of the grid
267: Not Collective
269: Input Parameter:
270: . dm - the DMStag object
272: Output Parameters:
273: + dof0 - the number of points per 0-cell (vertex/node)
274: . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
275: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
276: - dof3 - the number of points per 3-cell (element in 3D)
278: Level: beginner
280: .seealso: DMSTAG, DMStagGetCorners(), DMStagGetGhostCorners(), DMStagGetGlobalSizes(), DMStagGetStencilWidth(), DMStagGetBoundaryTypes(), DMStagGetLocationDOF(), DMDAGetDof()
281: @*/
282: PetscErrorCode DMStagGetDOF(DM dm,PetscInt *dof0,PetscInt *dof1,PetscInt *dof2,PetscInt *dof3)
283: {
284: const DM_Stag * const stag = (DM_Stag*)dm->data;
287: if (dof0) *dof0 = stag->dof[0];
288: if (dof1) *dof1 = stag->dof[1];
289: if (dof2) *dof2 = stag->dof[2];
290: if (dof3) *dof3 = stag->dof[3];
291: return 0;
292: }
294: /*@C
295: DMStagGetGhostCorners - return global element indices of the local region, including ghost points
297: Not Collective
299: Input Parameter:
300: . dm - the DMStag object
302: Output Parameters:
303: + x - the starting element index in the first direction
304: . y - the starting element index in the second direction
305: . z - the starting element index in the third direction
306: . m - the element width in the first direction
307: . n - the element width in the second direction
308: - p - the element width in the third direction
310: Notes:
311: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
313: Level: beginner
315: .seealso: DMSTAG, DMStagGetCorners(), DMDAGetGhostCorners()
316: @*/
317: PetscErrorCode DMStagGetGhostCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
318: {
319: const DM_Stag * const stag = (DM_Stag*)dm->data;
322: if (x) *x = stag->startGhost[0];
323: if (y) *y = stag->startGhost[1];
324: if (z) *z = stag->startGhost[2];
325: if (m) *m = stag->nGhost[0];
326: if (n) *n = stag->nGhost[1];
327: if (p) *p = stag->nGhost[2];
328: return 0;
329: }
331: /*@C
332: DMStagGetGlobalSizes - get global element counts
334: Not Collective
336: Input Parameter:
337: . dm - the DMStag object
339: Output Parameters:
340: . M,N,P - global element counts in each direction
342: Notes:
343: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
345: Level: beginner
347: .seealso: DMSTAG, DMStagGetLocalSizes(), DMDAGetInfo()
348: @*/
349: PetscErrorCode DMStagGetGlobalSizes(DM dm,PetscInt* M,PetscInt* N,PetscInt* P)
350: {
351: const DM_Stag * const stag = (DM_Stag*)dm->data;
354: if (M) *M = stag->N[0];
355: if (N) *N = stag->N[1];
356: if (P) *P = stag->N[2];
357: return 0;
358: }
360: /*@C
361: DMStagGetIsFirstRank - get boolean value for whether this rank is first in each direction in the rank grid
363: Not Collective
365: Input Parameter:
366: . dm - the DMStag object
368: Output Parameters:
369: . isFirstRank0,isFirstRank1,isFirstRank2 - whether this rank is first in each direction
371: Level: intermediate
373: Notes:
374: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
376: .seealso: DMSTAG, DMStagGetIsLastRank()
377: @*/
378: PetscErrorCode DMStagGetIsFirstRank(DM dm,PetscBool *isFirstRank0,PetscBool *isFirstRank1,PetscBool *isFirstRank2)
379: {
380: const DM_Stag * const stag = (DM_Stag*)dm->data;
383: if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
384: if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
385: if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
386: return 0;
387: }
389: /*@C
390: DMStagGetIsLastRank - get boolean value for whether this rank is last in each direction in the rank grid
392: Not Collective
394: Input Parameter:
395: . dm - the DMStag object
397: Output Parameters:
398: . isLastRank0,isLastRank1,isLastRank2 - whether this rank is last in each direction
400: Level: intermediate
402: Notes:
403: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
404: Level: intermediate
406: .seealso: DMSTAG, DMStagGetIsFirstRank()
407: @*/
408: PetscErrorCode DMStagGetIsLastRank(DM dm,PetscBool *isLastRank0,PetscBool *isLastRank1,PetscBool *isLastRank2)
409: {
410: const DM_Stag * const stag = (DM_Stag*)dm->data;
413: if (isLastRank0) *isLastRank0 = stag->lastRank[0];
414: if (isLastRank1) *isLastRank1 = stag->lastRank[1];
415: if (isLastRank2) *isLastRank2 = stag->lastRank[2];
416: return 0;
417: }
419: /*@C
420: DMStagGetLocalSizes - get local elementwise sizes
422: Not Collective
424: Input Parameter:
425: . dm - the DMStag object
427: Output Parameters:
428: . m,n,p - local element counts (excluding ghosts) in each direction
430: Notes:
431: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
432: Level: intermediate
434: Level: beginner
436: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetDOF(), DMStagGetNumRanks(), DMDAGetLocalInfo()
437: @*/
438: PetscErrorCode DMStagGetLocalSizes(DM dm,PetscInt* m,PetscInt* n,PetscInt* p)
439: {
440: const DM_Stag * const stag = (DM_Stag*)dm->data;
443: if (m) *m = stag->n[0];
444: if (n) *n = stag->n[1];
445: if (p) *p = stag->n[2];
446: return 0;
447: }
449: /*@C
450: DMStagGetNumRanks - get number of ranks in each direction in the global grid decomposition
452: Not Collective
454: Input Parameter:
455: . dm - the DMStag object
457: Output Parameters:
458: . nRanks0,nRanks1,nRanks2 - number of ranks in each direction in the grid decomposition
460: Notes:
461: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
462: Level: intermediate
464: Level: beginner
466: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetLocalSize(), DMStagSetNumRanks(), DMDAGetInfo()
467: @*/
468: PetscErrorCode DMStagGetNumRanks(DM dm,PetscInt *nRanks0,PetscInt *nRanks1,PetscInt *nRanks2)
469: {
470: const DM_Stag * const stag = (DM_Stag*)dm->data;
473: if (nRanks0) *nRanks0 = stag->nRanks[0];
474: if (nRanks1) *nRanks1 = stag->nRanks[1];
475: if (nRanks2) *nRanks2 = stag->nRanks[2];
476: return 0;
477: }
479: /*@C
480: DMStagGetEntries - get number of native entries in the global representation
482: Not Collective
484: Input Parameter:
485: . dm - the DMStag object
487: Output Parameters:
488: . entries - number of rank-native entries in the global representation
490: Note:
491: This is the number of entries on this rank for a global vector associated with dm.
493: Level: developer
495: .seealso: DMSTAG, DMStagGetDOF(), DMStagGetEntriesPerElement(), DMCreateLocalVector()
496: @*/
497: PetscErrorCode DMStagGetEntries(DM dm,PetscInt *entries)
498: {
499: const DM_Stag * const stag = (DM_Stag*)dm->data;
502: if (entries) *entries = stag->entries;
503: return 0;
504: }
506: /*@C
507: DMStagGetEntriesPerElement - get number of entries per element in the local representation
509: Not Collective
511: Input Parameter:
512: . dm - the DMStag object
514: Output Parameters:
515: . entriesPerElement - number of entries associated with each element in the local representation
517: Notes:
518: This is the natural block size for most local operations. In 1D it is equal to dof0 + dof1,
519: in 2D it is equal to dof0 + 2*dof1 + dof2, and in 3D it is equal to dof0 + 3*dof1 + 3*dof2 + dof3
521: Level: developer
523: .seealso: DMSTAG, DMStagGetDOF()
524: @*/
525: PetscErrorCode DMStagGetEntriesPerElement(DM dm,PetscInt *entriesPerElement)
526: {
527: const DM_Stag * const stag = (DM_Stag*)dm->data;
530: if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
531: return 0;
532: }
534: /*@C
535: DMStagGetStencilType - get elementwise ghost/halo stencil type
537: Not Collective
539: Input Parameter:
540: . dm - the DMStag object
542: Output Parameter:
543: . stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE
545: Level: beginner
547: .seealso: DMSTAG, DMStagSetStencilType(), DMStagGetStencilWidth, DMStagStencilType
548: @*/
549: PetscErrorCode DMStagGetStencilType(DM dm,DMStagStencilType *stencilType)
550: {
551: DM_Stag * const stag = (DM_Stag*)dm->data;
554: *stencilType = stag->stencilType;
555: return 0;
556: }
558: /*@C
559: DMStagGetStencilWidth - get elementwise stencil width
561: Not Collective
563: Input Parameter:
564: . dm - the DMStag object
566: Output Parameters:
567: . stencilWidth - stencil/halo/ghost width in elements
569: Level: beginner
571: .seealso: DMSTAG, DMStagSetStencilWidth(), DMStagGetStencilType(), DMDAGetStencilType()
572: @*/
573: PetscErrorCode DMStagGetStencilWidth(DM dm,PetscInt *stencilWidth)
574: {
575: const DM_Stag * const stag = (DM_Stag*)dm->data;
578: if (stencilWidth) *stencilWidth = stag->stencilWidth;
579: return 0;
580: }
582: /*@C
583: DMStagGetOwnershipRanges - get elements per rank in each direction
585: Not Collective
587: Input Parameter:
588: . dm - the DMStag object
590: Output Parameters:
591: + lx - ownership along x direction (optional)
592: . ly - ownership along y direction (optional)
593: - lz - ownership along z direction (optional)
595: Notes:
596: These correspond to the optional final arguments passed to DMStagCreate1d(), DMStagCreate2d(), and DMStagCreate3d().
598: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
600: In C you should not free these arrays, nor change the values in them.
601: They will only have valid values while the DMStag they came from still exists (has not been destroyed).
603: Level: intermediate
605: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagSetOwnershipRanges(), DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDAGetOwnershipRanges()
606: @*/
607: PetscErrorCode DMStagGetOwnershipRanges(DM dm,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
608: {
609: const DM_Stag * const stag = (DM_Stag*)dm->data;
612: if (lx) *lx = stag->l[0];
613: if (ly) *ly = stag->l[1];
614: if (lz) *lz = stag->l[2];
615: return 0;
616: }
618: /*@C
619: DMStagCreateCompatibleDMStag - create a compatible DMStag with different dof/stratum
621: Collective
623: Input Parameters:
624: + dm - the DMStag object
625: - dof0,dof1,dof2,dof3 - number of dof on each stratum in the new DMStag
627: Output Parameters:
628: . newdm - the new, compatible DMStag
630: Notes:
631: Dof supplied for strata too big for the dimension are ignored; these may be set to 0.
632: For example, for a 2-dimensional DMStag, dof2 sets the number of dof per element,
633: and dof3 is unused. For a 3-dimensional DMStag, dof3 sets the number of dof per element.
635: In contrast to DMDACreateCompatibleDMDA(), coordinates are not reused.
637: Level: intermediate
639: .seealso: DMSTAG, DMDACreateCompatibleDMDA(), DMGetCompatibility(), DMStagMigrateVec()
640: @*/
641: PetscErrorCode DMStagCreateCompatibleDMStag(DM dm,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DM *newdm)
642: {
644: DMStagDuplicateWithoutSetup(dm,PetscObjectComm((PetscObject)dm),newdm);
645: DMStagSetDOF(*newdm,dof0,dof1,dof2,dof3);
646: DMSetUp(*newdm);
647: return 0;
648: }
650: /*@C
651: DMStagGetLocationSlot - get index to use in accessing raw local arrays
653: Not Collective
655: Input Parameters:
656: + dm - the DMStag object
657: . loc - location relative to an element
658: - c - component
660: Output Parameter:
661: . slot - index to use
663: Notes:
664: Provides an appropriate index to use with DMStagVecGetArray() and friends.
665: This is required so that the user doesn't need to know about the ordering of
666: dof associated with each local element.
668: Level: beginner
670: .seealso: DMSTAG, DMStagVecGetArray(), DMStagVecGetArrayRead(), DMStagGetDOF(), DMStagGetEntriesPerElement()
671: @*/
672: PetscErrorCode DMStagGetLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt c,PetscInt *slot)
673: {
674: DM_Stag * const stag = (DM_Stag*)dm->data;
677: if (PetscDefined(USE_DEBUG)) {
678: PetscInt dof;
679: DMStagGetLocationDOF(dm,loc,&dof);
682: }
683: *slot = stag->locationOffsets[loc] + c;
684: return 0;
685: }
687: /*@C
688: DMStagMigrateVec - transfer a vector associated with a DMStag to a vector associated with a compatible DMStag
690: Collective
692: Input Parameters:
693: + dm - the source DMStag object
694: . vec - the source vector, compatible with dm
695: . dmTo - the compatible destination DMStag object
696: - vecTo - the destination vector, compatible with dmTo
698: Notes:
699: Extra dof are ignored, and unfilled dof are zeroed.
700: Currently only implemented to migrate global vectors to global vectors.
702: Level: advanced
704: .seealso: DMSTAG, DMStagCreateCompatibleDMStag(), DMGetCompatibility(), DMStagVecSplitToDMDA()
705: @*/
706: PetscErrorCode DMStagMigrateVec(DM dm,Vec vec,DM dmTo,Vec vecTo)
707: {
708: DM_Stag * const stag = (DM_Stag*)dm->data;
709: DM_Stag * const stagTo = (DM_Stag*)dmTo->data;
710: PetscInt nLocalTo,nLocal,dim,i,j,k;
711: PetscInt start[DMSTAG_MAX_DIM],startGhost[DMSTAG_MAX_DIM],n[DMSTAG_MAX_DIM],nExtra[DMSTAG_MAX_DIM],offset[DMSTAG_MAX_DIM];
712: Vec vecToLocal,vecLocal;
713: PetscBool compatible,compatibleSet;
714: const PetscScalar *arr;
715: PetscScalar *arrTo;
716: const PetscInt epe = stag->entriesPerElement;
717: const PetscInt epeTo = stagTo->entriesPerElement;
723: DMGetCompatibility(dm,dmTo,&compatible,&compatibleSet);
725: DMGetDimension(dm,&dim);
726: VecGetLocalSize(vec,&nLocal);
727: VecGetLocalSize(vecTo,&nLocalTo);
729: DMStagGetCorners(dm,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
730: DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],NULL,NULL,NULL);
731: for (i=0; i<DMSTAG_MAX_DIM; ++i) offset[i] = start[i]-startGhost[i];
733: /* Proceed by transferring to a local vector, copying, and transferring back to a global vector */
734: DMGetLocalVector(dm,&vecLocal);
735: DMGetLocalVector(dmTo,&vecToLocal);
736: DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal);
737: DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal);
738: VecGetArrayRead(vecLocal,&arr);
739: VecGetArray(vecToLocal,&arrTo);
740: /* Note that some superfluous copying of entries on partial dummy elements is done */
741: if (dim == 1) {
742: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
743: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
744: PetscInt si;
745: for (si=0; si<2; ++si) {
746: b += stag->dof[si];
747: bTo += stagTo->dof[si];
748: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[i*epeTo + dTo] = arr[i*epe + d];
749: for (; dTo < bTo ; ++dTo) arrTo[i*epeTo + dTo] = 0.0;
750: d = b;
751: }
752: }
753: } else if (dim == 2) {
754: const PetscInt epr = stag->nGhost[0] * epe;
755: const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
756: for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
757: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
758: const PetscInt base = j*epr + i*epe;
759: const PetscInt baseTo = j*eprTo + i*epeTo;
760: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
761: const PetscInt s[4] = {0,1,1,2}; /* Dimensions of points, in order */
762: PetscInt si;
763: for (si=0; si<4; ++si) {
764: b += stag->dof[s[si]];
765: bTo += stagTo->dof[s[si]];
766: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
767: for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
768: d = b;
769: }
770: }
771: }
772: } else if (dim == 3) {
773: const PetscInt epr = stag->nGhost[0] * epe;
774: const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
775: const PetscInt epl = stag->nGhost[1] * epr;
776: const PetscInt eplTo = stagTo->nGhost[1] * eprTo;
777: for (k=offset[2]; k<offset[2] + n[2] + nExtra[2]; ++k) {
778: for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
779: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
780: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
781: const PetscInt base = k*epl + j*epr + i*epe;
782: const PetscInt baseTo = k*eplTo + j*eprTo + i*epeTo;
783: const PetscInt s[8] = {0,1,1,2,1,2,2,3}; /* dimensions of points, in order */
784: PetscInt is;
785: for (is=0; is<8; ++is) {
786: b += stag->dof[s[is]];
787: bTo += stagTo->dof[s[is]];
788: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
789: for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
790: d = b;
791: }
792: }
793: }
794: }
795: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
796: VecRestoreArrayRead(vecLocal,&arr);
797: VecRestoreArray(vecToLocal,&arrTo);
798: DMRestoreLocalVector(dm,&vecLocal);
799: DMLocalToGlobalBegin(dmTo,vecToLocal,INSERT_VALUES,vecTo);
800: DMLocalToGlobalEnd(dmTo,vecToLocal,INSERT_VALUES,vecTo);
801: DMRestoreLocalVector(dmTo,&vecToLocal);
802: return 0;
803: }
805: /*@C
806: DMStagPopulateLocalToGlobalInjective - populate an internal 1-to-1 local-to-global map
808: Collective
810: Creates an internal object which explicitly maps a single local degree of
811: freedom to each global degree of freedom. This is used, if populated,
812: instead of SCATTER_REVERSE_LOCAL with the (1-to-many, in general)
813: global-to-local map, when DMLocalToGlobal() is called with INSERT_VALUES.
814: This allows usage, for example, even in the periodic, 1-rank case, where
815: the inverse of the global-to-local map, even when restricted to on-rank
816: communication, is non-injective. This is at the cost of storing an additional
817: VecScatter object inside each DMStag object.
819: Input Parameter:
820: . dm - the DMStag object
822: Notes:
823: In normal usage, library users shouldn't be concerned with this function,
824: as it is called during DMSetUp(), when required.
826: Returns immediately if the internal map is already populated.
828: Developer Notes:
829: This could, if desired, be moved up to a general DM routine. It would allow,
830: for example, DMDA to support DMLocalToGlobal() with INSERT_VALUES,
831: even in the single-rank periodic case.
833: Level: developer
835: .seealso: DMSTAG, DMLocalToGlobal(), VecScatter
836: @*/
837: PetscErrorCode DMStagPopulateLocalToGlobalInjective(DM dm)
838: {
839: PetscInt dim;
840: DM_Stag * const stag = (DM_Stag*)dm->data;
843: if (stag->ltog_injective) return 0; /* Don't re-populate */
844: DMGetDimension(dm,&dim);
845: switch (dim) {
846: case 1: DMStagPopulateLocalToGlobalInjective_1d(dm); break;
847: case 2: DMStagPopulateLocalToGlobalInjective_2d(dm); break;
848: case 3: DMStagPopulateLocalToGlobalInjective_3d(dm); break;
849: default: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
850: }
851: return 0;
852: }
854: static PetscErrorCode DMStagRestoreProductCoordinateArrays_Private(DM dm,void *arrX,void *arrY,void *arrZ,PetscBool read)
855: {
856: PetscInt dim,d;
857: void* arr[DMSTAG_MAX_DIM];
858: DM dmCoord;
861: DMGetDimension(dm,&dim);
863: arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
864: DMGetCoordinateDM(dm,&dmCoord);
865: for (d=0; d<dim; ++d) {
866: DM subDM;
867: Vec coord1d_local;
869: /* Ignore unrequested arrays */
870: if (!arr[d]) continue;
872: DMProductGetDM(dmCoord,d,&subDM);
873: DMGetCoordinatesLocal(subDM,&coord1d_local);
874: if (read) {
875: DMStagVecRestoreArrayRead(subDM,coord1d_local,arr[d]);
876: } else {
877: DMStagVecRestoreArray(subDM,coord1d_local,arr[d]);
878: }
879: }
880: return 0;
881: }
883: /*@C
884: DMStagRestoreProductCoordinateArrays - restore local array access
886: Logically Collective
888: Input Parameter:
889: . dm - the DMStag object
891: Output Parameters:
892: . arrX,arrY,arrZ - local 1D coordinate arrays
894: Level: intermediate
896: Notes:
897: 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:
899: $ DMGetCoordinateDM(dm,&cdm);
900: $ for (d=0; d<3; ++d) {
901: $ DM subdm;
902: $ Vec coor,coor_local;
904: $ DMProductGetDM(cdm,d,&subdm);
905: $ DMGetCoordinates(subdm,&coor);
906: $ DMGetCoordinatesLocal(subdm,&coor_local);
907: $ DMLocalToGlobal(subdm,coor_local,INSERT_VALUES,coor);
908: $ PetscPrintf(PETSC_COMM_WORLD,"Coordinates dim %D:\n",d);
909: $ VecView(coor,PETSC_VIEWER_STDOUT_WORLD);
910: $ }
912: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
913: @*/
914: PetscErrorCode DMStagRestoreProductCoordinateArrays(DM dm,void *arrX,void *arrY,void *arrZ)
915: {
916: DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_FALSE);
917: return 0;
918: }
920: /*@C
921: DMStagRestoreProductCoordinateArraysRead - restore local product array access, read-only
923: Logically Collective
925: Input Parameter:
926: . dm - the DMStag object
928: Output Parameters:
929: . arrX,arrY,arrZ - local 1D coordinate arrays
931: Level: intermediate
933: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
934: @*/
935: PetscErrorCode DMStagRestoreProductCoordinateArraysRead(DM dm,void *arrX,void *arrY,void *arrZ)
936: {
937: DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_TRUE);
938: return 0;
939: }
941: /*@C
942: DMStagSetBoundaryTypes - set DMStag boundary types
944: Logically Collective; boundaryType0, boundaryType1, and boundaryType2 must contain common values
946: Input Parameters:
947: + dm - the DMStag object
948: - boundaryType0,boundaryType1,boundaryType2 - boundary types in each direction
950: Notes:
951: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
953: Level: advanced
955: .seealso: DMSTAG, DMBoundaryType, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDASetBoundaryType()
956: @*/
957: PetscErrorCode DMStagSetBoundaryTypes(DM dm,DMBoundaryType boundaryType0,DMBoundaryType boundaryType1,DMBoundaryType boundaryType2)
958: {
959: DM_Stag * const stag = (DM_Stag*)dm->data;
960: PetscInt dim;
962: DMGetDimension(dm,&dim);
968: stag->boundaryType[0] = boundaryType0;
969: if (dim > 1) stag->boundaryType[1] = boundaryType1;
970: if (dim > 2) stag->boundaryType[2] = boundaryType2;
971: return 0;
972: }
974: /*@C
975: DMStagSetCoordinateDMType - set DM type to store coordinates
977: Logically Collective; dmtype must contain common value
979: Input Parameters:
980: + dm - the DMStag object
981: - dmtype - DMtype for coordinates, either DMSTAG or DMPRODUCT
983: Level: advanced
985: .seealso: DMSTAG, DMPRODUCT, DMGetCoordinateDM(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMType
986: @*/
987: PetscErrorCode DMStagSetCoordinateDMType(DM dm,DMType dmtype)
988: {
989: DM_Stag * const stag = (DM_Stag*)dm->data;
992: PetscFree(stag->coordinateDMType);
993: PetscStrallocpy(dmtype,(char**)&stag->coordinateDMType);
994: return 0;
995: }
997: /*@C
998: DMStagSetDOF - set dof/stratum
1000: Logically Collective; dof0, dof1, dof2, and dof3 must contain common values
1002: Input Parameters:
1003: + dm - the DMStag object
1004: - dof0,dof1,dof2,dof3 - dof per stratum
1006: Notes:
1007: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1009: Level: advanced
1011: .seealso: DMSTAG, DMDASetDof()
1012: @*/
1013: PetscErrorCode DMStagSetDOF(DM dm,PetscInt dof0, PetscInt dof1,PetscInt dof2,PetscInt dof3)
1014: {
1015: DM_Stag * const stag = (DM_Stag*)dm->data;
1016: PetscInt dim;
1024: DMGetDimension(dm,&dim);
1029: stag->dof[0] = dof0;
1030: stag->dof[1] = dof1;
1031: if (dim > 1) stag->dof[2] = dof2;
1032: if (dim > 2) stag->dof[3] = dof3;
1033: return 0;
1034: }
1036: /*@C
1037: DMStagSetNumRanks - set ranks in each direction in the global rank grid
1039: Logically Collective; nRanks0, nRanks1, and nRanks2 must contain common values
1041: Input Parameters:
1042: + dm - the DMStag object
1043: - nRanks0,nRanks1,nRanks2 - number of ranks in each direction
1045: Notes:
1046: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1048: Level: developer
1050: .seealso: DMSTAG, DMDASetNumProcs()
1051: @*/
1052: PetscErrorCode DMStagSetNumRanks(DM dm,PetscInt nRanks0,PetscInt nRanks1,PetscInt nRanks2)
1053: {
1054: DM_Stag * const stag = (DM_Stag*)dm->data;
1055: PetscInt dim;
1062: DMGetDimension(dm,&dim);
1066: if (nRanks0) stag->nRanks[0] = nRanks0;
1067: if (dim > 1 && nRanks1) stag->nRanks[1] = nRanks1;
1068: if (dim > 2 && nRanks2) stag->nRanks[2] = nRanks2;
1069: return 0;
1070: }
1072: /*@C
1073: DMStagSetStencilType - set elementwise ghost/halo stencil type
1075: Logically Collective; stencilType must contain common value
1077: Input Parameters:
1078: + dm - the DMStag object
1079: - stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE
1081: Level: beginner
1083: .seealso: DMSTAG, DMStagGetStencilType(), DMStagSetStencilWidth(), DMStagStencilType
1084: @*/
1085: PetscErrorCode DMStagSetStencilType(DM dm,DMStagStencilType stencilType)
1086: {
1087: DM_Stag * const stag = (DM_Stag*)dm->data;
1092: stag->stencilType = stencilType;
1093: return 0;
1094: }
1096: /*@C
1097: DMStagSetStencilWidth - set elementwise stencil width
1099: Logically Collective; stencilWidth must contain common value
1101: Input Parameters:
1102: + dm - the DMStag object
1103: - stencilWidth - stencil/halo/ghost width in elements
1105: Notes:
1106: The width value is not used when DMSTAG_STENCIL_NONE is specified.
1108: Level: beginner
1110: .seealso: DMSTAG, DMStagGetStencilWidth(), DMStagGetStencilType(), DMStagStencilType
1111: @*/
1112: PetscErrorCode DMStagSetStencilWidth(DM dm,PetscInt stencilWidth)
1113: {
1114: DM_Stag * const stag = (DM_Stag*)dm->data;
1120: stag->stencilWidth = stencilWidth;
1121: return 0;
1122: }
1124: /*@C
1125: DMStagSetGlobalSizes - set global element counts in each direction
1127: Logically Collective; N0, N1, and N2 must contain common values
1129: Input Parameters:
1130: + dm - the DMStag object
1131: - N0,N1,N2 - global elementwise sizes
1133: Notes:
1134: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1136: Level: advanced
1138: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMDASetSizes()
1139: @*/
1140: PetscErrorCode DMStagSetGlobalSizes(DM dm,PetscInt N0,PetscInt N1,PetscInt N2)
1141: {
1142: DM_Stag * const stag = (DM_Stag*)dm->data;
1143: PetscInt dim;
1147: DMGetDimension(dm,&dim);
1151: if (N0) stag->N[0] = N0;
1152: if (N1) stag->N[1] = N1;
1153: if (N2) stag->N[2] = N2;
1154: return 0;
1155: }
1157: /*@C
1158: DMStagSetOwnershipRanges - set elements per rank in each direction
1160: Logically Collective; lx, ly, and lz must contain common values
1162: Input Parameters:
1163: + dm - the DMStag object
1164: - lx,ly,lz - element counts for each rank in each direction
1166: Notes:
1167: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
1169: Level: developer
1171: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagGetOwnershipRanges(), DMDASetOwnershipRanges()
1172: @*/
1173: PetscErrorCode DMStagSetOwnershipRanges(DM dm,PetscInt const *lx,PetscInt const *ly,PetscInt const *lz)
1174: {
1175: DM_Stag * const stag = (DM_Stag*)dm->data;
1176: const PetscInt *lin[3];
1177: PetscInt d,dim;
1181: lin[0] = lx; lin[1] = ly; lin[2] = lz;
1182: DMGetDimension(dm,&dim);
1183: for (d=0; d<dim; ++d) {
1184: if (lin[d]) {
1186: if (!stag->l[d]) {
1187: PetscMalloc1(stag->nRanks[d], &stag->l[d]);
1188: }
1189: PetscArraycpy(stag->l[d], lin[d], stag->nRanks[d]);
1190: }
1191: }
1192: return 0;
1193: }
1195: /*@C
1196: DMStagSetUniformCoordinates - set DMStag coordinates to be a uniform grid
1198: Collective
1200: Input Parameters:
1201: + dm - the DMStag object
1202: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1204: Notes:
1205: DMStag supports 2 different types of coordinate DM: DMSTAG and DMPRODUCT.
1206: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1208: Local coordinates are populated (using DMSetCoordinatesLocal()), linearly
1209: extrapolated to ghost cells, including those outside the physical domain.
1210: This is also done in case of periodic boundaries, meaning that the same
1211: global point may have different coordinates in different local representations,
1212: which are equivalent assuming a periodicity implied by the arguments to this function,
1213: i.e. two points are equivalent if their difference is a multiple of (xmax - xmin)
1214: in the x direction, (ymax - ymin) in the y direction, and (zmax - zmin) in the z direction.
1216: Level: advanced
1218: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType(), DMGetCoordinateDM(), DMGetCoordinates(), DMDASetUniformCoordinates(), DMBoundaryType
1219: @*/
1220: PetscErrorCode DMStagSetUniformCoordinates(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1221: {
1222: DM_Stag * const stag = (DM_Stag*)dm->data;
1223: PetscBool flg_stag,flg_product;
1228: PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg_stag);
1229: PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg_product);
1230: if (flg_stag) {
1231: DMStagSetUniformCoordinatesExplicit(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1232: } else if (flg_product) {
1233: DMStagSetUniformCoordinatesProduct(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1234: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported DM Type %s",stag->coordinateDMType);
1235: return 0;
1236: }
1238: /*@C
1239: DMStagSetUniformCoordinatesExplicit - set DMStag coordinates to be a uniform grid, storing all values
1241: Collective
1243: Input Parameters:
1244: + dm - the DMStag object
1245: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1247: Notes:
1248: DMStag supports 2 different types of coordinate DM: either another DMStag, or a DMProduct.
1249: If the grid is orthogonal, using DMProduct should be more efficient.
1251: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1253: See the manual page for DMStagSetUniformCoordinates() for information on how
1254: coordinates for dummy cells outside the physical domain boundary are populated.
1256: Level: beginner
1258: .seealso: DMSTAG, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType()
1259: @*/
1260: PetscErrorCode DMStagSetUniformCoordinatesExplicit(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1261: {
1262: DM_Stag * const stag = (DM_Stag*)dm->data;
1263: PetscInt dim;
1264: PetscBool flg;
1268: PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg);
1270: DMStagSetCoordinateDMType(dm,DMSTAG);
1271: DMGetDimension(dm,&dim);
1272: switch (dim) {
1273: case 1: DMStagSetUniformCoordinatesExplicit_1d(dm,xmin,xmax); break;
1274: case 2: DMStagSetUniformCoordinatesExplicit_2d(dm,xmin,xmax,ymin,ymax); break;
1275: case 3: DMStagSetUniformCoordinatesExplicit_3d(dm,xmin,xmax,ymin,ymax,zmin,zmax); break;
1276: default: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1277: }
1278: return 0;
1279: }
1281: /*@C
1282: DMStagSetUniformCoordinatesProduct - create uniform coordinates, as a product of 1D arrays
1284: 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.
1286: Collective
1288: Input Parameters:
1289: + dm - the DMStag object
1290: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1292: Notes:
1293: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1295: The per-dimension 1-dimensional DMStag objects that comprise the product
1296: always have active 0-cells (vertices, element boundaries) and 1-cells
1297: (element centers).
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: intermediate
1304: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetCoordinateDMType()
1305: @*/
1306: PetscErrorCode DMStagSetUniformCoordinatesProduct(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: DM dmc;
1310: PetscInt dim,d,dof0,dof1;
1311: PetscBool flg;
1315: PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg);
1317: DMStagSetCoordinateDMType(dm,DMPRODUCT);
1318: DMGetCoordinateDM(dm,&dmc);
1319: DMGetDimension(dm,&dim);
1321: /* Create 1D sub-DMs, living on subcommunicators.
1322: Always include both vertex and element dof, regardless of the active strata of the DMStag */
1323: dof0 = 1;
1324: dof1 = 1;
1326: for (d=0; d<dim; ++d) {
1327: DM subdm;
1328: MPI_Comm subcomm;
1329: PetscMPIInt color;
1330: const PetscMPIInt key = 0; /* let existing rank break ties */
1332: /* Choose colors based on position in the plane orthogonal to this dim, and split */
1333: switch (d) {
1334: case 0: color = (dim > 1 ? stag->rank[1] : 0) + (dim > 2 ? stag->nRanks[1]*stag->rank[2] : 0); break;
1335: case 1: color = stag->rank[0] + (dim > 2 ? stag->nRanks[0]*stag->rank[2] : 0); break;
1336: case 2: color = stag->rank[0] + stag->nRanks[0]*stag->rank[1] ; break;
1337: default: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1338: }
1339: MPI_Comm_split(PetscObjectComm((PetscObject)dm),color,key,&subcomm);
1341: /* Create sub-DMs living on these new communicators (which are destroyed by DMProduct) */
1342: DMStagCreate1d(subcomm,stag->boundaryType[d],stag->N[d],dof0,dof1,stag->stencilType,stag->stencilWidth,stag->l[d],&subdm);
1343: DMSetUp(subdm);
1344: switch (d) {
1345: case 0:
1346: DMStagSetUniformCoordinatesExplicit(subdm,xmin,xmax,0,0,0,0);
1347: break;
1348: case 1:
1349: DMStagSetUniformCoordinatesExplicit(subdm,ymin,ymax,0,0,0,0);
1350: break;
1351: case 2:
1352: DMStagSetUniformCoordinatesExplicit(subdm,zmin,zmax,0,0,0,0);
1353: break;
1354: default: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1355: }
1356: DMProductSetDM(dmc,d,subdm);
1357: DMProductSetDimensionIndex(dmc,d,0);
1358: DMDestroy(&subdm);
1359: MPI_Comm_free(&subcomm);
1360: }
1361: return 0;
1362: }
1364: /*@C
1365: DMStagVecGetArray - get access to local array
1367: Logically Collective
1369: This function returns a (dim+1)-dimensional array for a dim-dimensional
1370: DMStag.
1372: The first 1-3 dimensions indicate an element in the global
1373: numbering, using the standard C ordering.
1375: The final dimension in this array corresponds to a degree
1376: of freedom with respect to this element, for example corresponding to
1377: the element or one of its neighboring faces, edges, or vertices.
1379: For example, for a 3D DMStag, indexing is array[k][j][i][idx], where k is the
1380: index in the z-direction, j is the index in the y-direction, and i is the
1381: index in the x-direction.
1383: "idx" is obtained with DMStagGetLocationSlot(), since the correct offset
1384: into the (dim+1)-dimensional C array depends on the grid size and the number
1385: of dof stored at each location.
1387: Input Parameters:
1388: + dm - the DMStag object
1389: - vec - the Vec object
1391: Output Parameters:
1392: . array - the array
1394: Notes:
1395: DMStagVecRestoreArray() must be called, once finished with the array
1397: Level: beginner
1399: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArray(), DMDAVecGetArrayDOF()
1400: @*/
1401: PetscErrorCode DMStagVecGetArray(DM dm,Vec vec,void *array)
1402: {
1403: DM_Stag * const stag = (DM_Stag*)dm->data;
1404: PetscInt dim;
1405: PetscInt nLocal;
1409: DMGetDimension(dm,&dim);
1410: VecGetLocalSize(vec,&nLocal);
1412: switch (dim) {
1413: case 1:
1414: VecGetArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1415: break;
1416: case 2:
1417: VecGetArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1418: break;
1419: case 3:
1420: 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);
1421: break;
1422: default: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1423: }
1424: return 0;
1425: }
1427: /*@C
1428: DMStagVecGetArrayRead - get read-only access to a local array
1430: Logically Collective
1432: See the man page for DMStagVecGetArray() for more information.
1434: Input Parameters:
1435: + dm - the DMStag object
1436: - vec - the Vec object
1438: Output Parameters:
1439: . array - the read-only array
1441: Notes:
1442: DMStagVecRestoreArrayRead() must be called, once finished with the array
1444: Level: beginner
1446: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArrayRead(), DMDAVecGetArrayDOFRead()
1447: @*/
1448: PetscErrorCode DMStagVecGetArrayRead(DM dm,Vec vec,void *array)
1449: {
1450: DM_Stag * const stag = (DM_Stag*)dm->data;
1451: PetscInt dim;
1452: PetscInt nLocal;
1456: DMGetDimension(dm,&dim);
1457: VecGetLocalSize(vec,&nLocal);
1459: switch (dim) {
1460: case 1:
1461: VecGetArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1462: break;
1463: case 2:
1464: VecGetArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1465: break;
1466: case 3:
1467: 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);
1468: break;
1469: default: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1470: }
1471: return 0;
1472: }
1474: /*@C
1475: DMStagVecRestoreArray - restore access to a raw array
1477: Logically Collective
1479: Input Parameters:
1480: + dm - the DMStag object
1481: - vec - the Vec object
1483: Output Parameters:
1484: . array - the array
1486: Level: beginner
1488: .seealso: DMSTAG, DMStagVecGetArray(), DMDAVecRestoreArray(), DMDAVecRestoreArrayDOF()
1489: @*/
1490: PetscErrorCode DMStagVecRestoreArray(DM dm,Vec vec,void *array)
1491: {
1492: DM_Stag * const stag = (DM_Stag*)dm->data;
1493: PetscInt dim;
1494: PetscInt nLocal;
1498: DMGetDimension(dm,&dim);
1499: VecGetLocalSize(vec,&nLocal);
1501: switch (dim) {
1502: case 1:
1503: VecRestoreArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1504: break;
1505: case 2:
1506: VecRestoreArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1507: break;
1508: case 3:
1509: 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);
1510: break;
1511: default: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1512: }
1513: return 0;
1514: }
1516: /*@C
1517: DMStagVecRestoreArrayRead - restore read-only access to a raw array
1519: Logically Collective
1521: Input Parameters:
1522: + dm - the DMStag object
1523: - vec - the Vec object
1525: Output Parameters:
1526: . array - the read-only array
1528: Level: beginner
1530: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMDAVecRestoreArrayRead(), DMDAVecRestoreArrayDOFRead()
1531: @*/
1532: PetscErrorCode DMStagVecRestoreArrayRead(DM dm,Vec vec,void *array)
1533: {
1534: DM_Stag * const stag = (DM_Stag*)dm->data;
1535: PetscInt dim;
1536: PetscInt nLocal;
1540: DMGetDimension(dm,&dim);
1541: VecGetLocalSize(vec,&nLocal);
1543: switch (dim) {
1544: case 1:
1545: VecRestoreArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1546: break;
1547: case 2:
1548: VecRestoreArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1549: break;
1550: case 3:
1551: 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);
1552: break;
1553: default: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1554: }
1555: return 0;
1556: }