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: PetscErrorCode ierr;
22: const DM_Stag * const stag = (DM_Stag*)dm->data;
23: PetscInt dim;
27: DMGetDimension(dm,&dim);
28: if (boundaryTypeX) *boundaryTypeX = stag->boundaryType[0];
29: if (boundaryTypeY && dim > 1) *boundaryTypeY = stag->boundaryType[1];
30: if (boundaryTypeZ && dim > 2) *boundaryTypeZ = stag->boundaryType[2];
31: return(0);
32: }
34: static PetscErrorCode DMStagGetProductCoordinateArrays_Private(DM dm,void* arrX,void* arrY,void* arrZ,PetscBool read)
35: {
37: PetscInt dim,d,dofCheck[DMSTAG_MAX_STRATA],s;
38: DM dmCoord;
39: void* arr[DMSTAG_MAX_DIM];
40: PetscBool checkDof;
44: DMGetDimension(dm,&dim);
45: if (dim > DMSTAG_MAX_DIM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for %D dimensions",dim);
46: arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
47: DMGetCoordinateDM(dm,&dmCoord);
48: if (!dmCoord) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM does not have a coordinate DM");
49: {
50: PetscBool isProduct;
51: DMType dmType;
52: DMGetType(dmCoord,&dmType);
53: PetscStrcmp(DMPRODUCT,dmType,&isProduct);
54: if (!isProduct) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is not of type DMPRODUCT");
55: }
56: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
57: checkDof = PETSC_FALSE;
58: for (d=0; d<dim; ++d) {
59: DM subDM;
60: DMType dmType;
61: PetscBool isStag;
62: PetscInt dof[DMSTAG_MAX_STRATA],subDim;
63: Vec coord1d_local;
65: /* Ignore unrequested arrays */
66: if (!arr[d]) continue;
68: DMProductGetDM(dmCoord,d,&subDM);
69: if (!subDM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is missing sub DM %D",d);
70: DMGetDimension(subDM,&subDim);
71: if (subDim != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of dimension 1");
72: DMGetType(subDM,&dmType);
73: PetscStrcmp(DMSTAG,dmType,&isStag);
74: if (!isStag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of type DMSTAG");
75: DMStagGetDOF(subDM,&dof[0],&dof[1],&dof[2],&dof[3]);
76: if (!checkDof) {
77: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
78: checkDof = PETSC_TRUE;
79: } else {
80: for (s=0; s<DMSTAG_MAX_STRATA; ++s) {
81: if (dofCheck[s] != dof[s]) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DMs have different dofs");
82: }
83: }
84: DMGetCoordinatesLocal(subDM,&coord1d_local);
85: if (read) {
86: DMStagVecGetArrayRead(subDM,coord1d_local,arr[d]);
87: } else {
88: DMStagVecGetArray(subDM,coord1d_local,arr[d]);
89: }
90: }
91: return(0);
92: }
94: /*@C
95: DMStagGetProductCoordinateArrays - extract local product coordinate arrays, one per dimension
97: Logically Collective
99: A high-level helper function to quickly extract local coordinate arrays.
101: Note that 2-dimensional arrays are returned. See
102: DMStagVecGetArray(), which is called internally to produce these arrays
103: representing coordinates on elements and vertices (element boundaries)
104: for a 1-dimensional DMStag in each coordinate direction.
106: One should use DMStagGetProductCoordinateLocationSlot() to determine appropriate
107: indices for the second dimension in these returned arrays. This function
108: checks that the coordinate array is a suitable product of 1-dimensional
109: DMStag objects.
111: Input Parameter:
112: . dm - the DMStag object
114: Output Parameters:
115: . arrX,arrY,arrZ - local 1D coordinate arrays
117: Level: intermediate
119: .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArraysRead(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagGetProductCoordinateLocationSlot()
120: @*/
121: PetscErrorCode DMStagGetProductCoordinateArrays(DM dm,void* arrX,void* arrY,void* arrZ)
122: {
126: DMStagGetProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_FALSE);
127: return(0);
128: }
130: /*@C
131: DMStagGetProductCoordinateArraysRead - extract product coordinate arrays, read-only
133: Logically Collective
135: See the man page for DMStagGetProductCoordinateArrays() for more information.
137: Input Parameter:
138: . dm - the DMStag object
140: Output Parameters:
141: . arrX,arrY,arrZ - local 1D coordinate arrays
143: Level: intermediate
145: .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArrays(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagGetProductCoordinateLocationSlot()
146: @*/
147: PetscErrorCode DMStagGetProductCoordinateArraysRead(DM dm,void* arrX,void* arrY,void* arrZ)
148: {
152: DMStagGetProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_TRUE);
153: return(0);
154: }
156: /*@C
157: DMStagGetProductCoordinateLocationSlot - get slot for use with local product coordinate arrays
159: Not Collective
161: High-level helper function to get slot indices for 1D coordinate DMs,
162: for use with DMStagGetProductCoordinateArrays() and related functions.
164: Input Parameters:
165: + dm - the DMStag object
166: - loc - the grid location
168: Output Parameter:
169: . slot - the index to use in local arrays
171: Notes:
172: For loc, one should use DMSTAG_LEFT, DMSTAG_ELEMENT, or DMSTAG_RIGHT for "previous", "center" and "next"
173: locations, respectively, in each dimension.
174: One can equivalently use DMSTAG_DOWN or DMSTAG_BACK in place of DMSTAG_LEFT,
175: and DMSTAG_UP or DMSTACK_FRONT in place of DMSTAG_RIGHT;
177: This function checks that the coordinates are actually set up so that using the
178: slots from any of the 1D coordinate sub-DMs are valid for all the 1D coordinate sub-DMs.
180: Level: intermediate
182: .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead(), DMStagSetUniformCoordinates()
183: @*/
184: PETSC_EXTERN PetscErrorCode DMStagGetProductCoordinateLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt *slot)
185: {
187: DM dmCoord;
188: PetscInt dim,dofCheck[DMSTAG_MAX_STRATA],s,d;
192: DMGetDimension(dm,&dim);
193: DMGetCoordinateDM(dm,&dmCoord);
194: if (!dmCoord) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM does not have a coordinate DM");
195: {
196: PetscBool isProduct;
197: DMType dmType;
198: DMGetType(dmCoord,&dmType);
199: PetscStrcmp(DMPRODUCT,dmType,&isProduct);
200: if (!isProduct) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is not of type DMPRODUCT");
201: }
202: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
203: for (d=0; d<dim; ++d) {
204: DM subDM;
205: DMType dmType;
206: PetscBool isStag;
207: PetscInt dof[DMSTAG_MAX_STRATA],subDim;
208: DMProductGetDM(dmCoord,d,&subDM);
209: if (!subDM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is missing sub DM %D",d);
210: DMGetDimension(subDM,&subDim);
211: if (subDim != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of dimension 1");
212: DMGetType(subDM,&dmType);
213: PetscStrcmp(DMSTAG,dmType,&isStag);
214: if (!isStag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of type DMSTAG");
215: DMStagGetDOF(subDM,&dof[0],&dof[1],&dof[2],&dof[3]);
216: if (d == 0) {
217: const PetscInt component = 0;
218: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
219: DMStagGetLocationSlot(subDM,loc,component,slot);
220: } else {
221: for (s=0; s<DMSTAG_MAX_STRATA; ++s) {
222: if (dofCheck[s] != dof[s]) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DMs have different dofs");
223: }
224: }
225: }
226: return(0);
227: }
229: /*@C
230: DMStagGetCorners - return global element indices of the local region (excluding ghost points)
232: Not Collective
234: Input Parameter:
235: . dm - the DMStag object
237: Output Parameters:
238: + x - starting element index in first direction
239: . y - starting element index in second direction
240: . z - starting element index in third direction
241: . m - element width in first direction
242: . n - element width in second direction
243: . p - element width in third direction
244: . nExtrax - number of extra partial elements in first direction
245: . nExtray - number of extra partial elements in second direction
246: - nExtraz - number of extra partial elements in third direction
248: Notes:
249: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
251: The number of extra partial elements is either 1 or 0.
252: The value is 1 on right, top, and front non-periodic domain ("physical") boundaries,
253: in the x, y, and z directions respectively, and otherwise 0.
255: Level: beginner
257: .seealso: DMSTAG, DMStagGetGhostCorners(), DMDAGetCorners()
258: @*/
259: PetscErrorCode DMStagGetCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *nExtrax,PetscInt *nExtray,PetscInt *nExtraz)
260: {
261: const DM_Stag * const stag = (DM_Stag*)dm->data;
265: if (x) *x = stag->start[0];
266: if (y) *y = stag->start[1];
267: if (z) *z = stag->start[2];
268: if (m) *m = stag->n[0];
269: if (n) *n = stag->n[1];
270: if (p) *p = stag->n[2];
271: if (nExtrax) *nExtrax = stag->boundaryType[0] != DM_BOUNDARY_PERIODIC && stag->lastRank[0] ? 1 : 0;
272: if (nExtray) *nExtray = stag->boundaryType[1] != DM_BOUNDARY_PERIODIC && stag->lastRank[1] ? 1 : 0;
273: if (nExtraz) *nExtraz = stag->boundaryType[2] != DM_BOUNDARY_PERIODIC && stag->lastRank[2] ? 1 : 0;
274: return(0);
275: }
277: /*@C
278: DMStagGetDOF - get number of DOF associated with each stratum of the grid
280: Not Collective
282: Input Parameter:
283: . dm - the DMStag object
285: Output Parameters:
286: + dof0 - the number of points per 0-cell (vertex/node)
287: . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
288: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
289: - dof3 - the number of points per 3-cell (element in 3D)
291: Level: beginner
293: .seealso: DMSTAG, DMStagGetCorners(), DMStagGetGhostCorners(), DMStagGetGlobalSizes(), DMStagGetStencilWidth(), DMStagGetBoundaryTypes(), DMStagGetLocationDOF(), DMDAGetDof()
294: @*/
295: PetscErrorCode DMStagGetDOF(DM dm,PetscInt *dof0,PetscInt *dof1,PetscInt *dof2,PetscInt *dof3)
296: {
297: const DM_Stag * const stag = (DM_Stag*)dm->data;
301: if (dof0) *dof0 = stag->dof[0];
302: if (dof1) *dof1 = stag->dof[1];
303: if (dof2) *dof2 = stag->dof[2];
304: if (dof3) *dof3 = stag->dof[3];
305: return(0);
306: }
308: /*@C
309: DMStagGetGhostCorners - return global element indices of the local region, including ghost points
311: Not Collective
313: Input Parameter:
314: . dm - the DMStag object
316: Output Parameters:
317: + x - the starting element index in the first direction
318: . y - the starting element index in the second direction
319: . z - the starting element index in the third direction
320: . m - the element width in the first direction
321: . n - the element width in the second direction
322: - p - the element width in the third direction
324: Notes:
325: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
327: Level: beginner
329: .seealso: DMSTAG, DMStagGetCorners(), DMDAGetGhostCorners()
330: @*/
331: PetscErrorCode DMStagGetGhostCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
332: {
333: const DM_Stag * const stag = (DM_Stag*)dm->data;
337: if (x) *x = stag->startGhost[0];
338: if (y) *y = stag->startGhost[1];
339: if (z) *z = stag->startGhost[2];
340: if (m) *m = stag->nGhost[0];
341: if (n) *n = stag->nGhost[1];
342: if (p) *p = stag->nGhost[2];
343: return(0);
344: }
346: /*@C
347: DMStagGetGlobalSizes - get global element counts
349: Not Collective
351: Input Parameter:
352: . dm - the DMStag object
354: Output Parameters:
355: . M,N,P - global element counts in each direction
357: Notes:
358: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
360: Level: beginner
362: .seealso: DMSTAG, DMStagGetLocalSizes(), DMDAGetInfo()
363: @*/
364: PetscErrorCode DMStagGetGlobalSizes(DM dm,PetscInt* M,PetscInt* N,PetscInt* P)
365: {
366: const DM_Stag * const stag = (DM_Stag*)dm->data;
370: if (M) *M = stag->N[0];
371: if (N) *N = stag->N[1];
372: if (P) *P = stag->N[2];
373: return(0);
374: }
376: /*@C
377: DMStagGetIsFirstRank - get boolean value for whether this rank is first in each direction in the rank grid
379: Not Collective
381: Input Parameter:
382: . dm - the DMStag object
384: Output Parameters:
385: . isFirstRank0,isFirstRank1,isFirstRank2 - whether this rank is first in each direction
387: Level: intermediate
389: Notes:
390: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
392: .seealso: DMSTAG, DMStagGetIsLastRank()
393: @*/
394: PetscErrorCode DMStagGetIsFirstRank(DM dm,PetscBool *isFirstRank0,PetscBool *isFirstRank1,PetscBool *isFirstRank2)
395: {
396: const DM_Stag * const stag = (DM_Stag*)dm->data;
400: if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
401: if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
402: if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
403: return(0);
404: }
406: /*@C
407: DMStagGetIsLastRank - get boolean value for whether this rank is last in each direction in the rank grid
409: Not Collective
411: Input Parameter:
412: . dm - the DMStag object
414: Output Parameters:
415: . isLastRank0,isLastRank1,isLastRank2 - whether this rank is last in each direction
417: Level: intermediate
419: Notes:
420: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
421: Level: intermediate
423: .seealso: DMSTAG, DMStagGetIsFirstRank()
424: @*/
425: PetscErrorCode DMStagGetIsLastRank(DM dm,PetscBool *isLastRank0,PetscBool *isLastRank1,PetscBool *isLastRank2)
426: {
427: const DM_Stag * const stag = (DM_Stag*)dm->data;
431: if (isLastRank0) *isLastRank0 = stag->lastRank[0];
432: if (isLastRank1) *isLastRank1 = stag->lastRank[1];
433: if (isLastRank2) *isLastRank2 = stag->lastRank[2];
434: return(0);
435: }
437: /*@C
438: DMStagGetLocalSizes - get local elementwise sizes
440: Not Collective
442: Input Parameter:
443: . dm - the DMStag object
445: Output Parameters:
446: . m,n,p - local element counts (excluding ghosts) in each direction
448: Notes:
449: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
450: Level: intermediate
452: Level: beginner
454: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetDOF(), DMStagGetNumRanks(), DMDAGetLocalInfo()
455: @*/
456: PetscErrorCode DMStagGetLocalSizes(DM dm,PetscInt* m,PetscInt* n,PetscInt* p)
457: {
458: const DM_Stag * const stag = (DM_Stag*)dm->data;
462: if (m) *m = stag->n[0];
463: if (n) *n = stag->n[1];
464: if (p) *p = stag->n[2];
465: return(0);
466: }
468: /*@C
469: DMStagGetNumRanks - get number of ranks in each direction in the global grid decomposition
471: Not Collective
473: Input Parameter:
474: . dm - the DMStag object
476: Output Parameters:
477: . nRanks0,nRanks1,nRanks2 - number of ranks in each direction in the grid decomposition
479: Notes:
480: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
481: Level: intermediate
483: Level: beginner
485: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetLocalSize(), DMStagSetNumRanks(), DMDAGetInfo()
486: @*/
487: PetscErrorCode DMStagGetNumRanks(DM dm,PetscInt *nRanks0,PetscInt *nRanks1,PetscInt *nRanks2)
488: {
489: const DM_Stag * const stag = (DM_Stag*)dm->data;
493: if (nRanks0) *nRanks0 = stag->nRanks[0];
494: if (nRanks1) *nRanks1 = stag->nRanks[1];
495: if (nRanks2) *nRanks2 = stag->nRanks[2];
496: return(0);
497: }
499: /*@C
500: DMStagGetEntries - get number of native entries in the global representation
502: Not Collective
504: Input Parameter:
505: . dm - the DMStag object
507: Output Parameters:
508: . entries - number of rank-native entries in the global representation
510: Note:
511: This is the number of entries on this rank for a global vector associated with dm.
513: Level: developer
515: .seealso: DMSTAG, DMStagGetDOF(), DMStagGetEntriesPerElement(), DMCreateLocalVector()
516: @*/
517: PetscErrorCode DMStagGetEntries(DM dm,PetscInt *entries)
518: {
519: const DM_Stag * const stag = (DM_Stag*)dm->data;
523: if (entries) *entries = stag->entries;
524: return(0);
525: }
527: /*@C
528: DMStagGetEntriesPerElement - get number of entries per element in the local representation
530: Not Collective
532: Input Parameter:
533: . dm - the DMStag object
535: Output Parameters:
536: . entriesPerElement - number of entries associated with each element in the local representation
538: Notes:
539: This is the natural block size for most local operations. In 1D it is equal to dof0 + dof1,
540: in 2D it is equal to dof0 + 2*dof1 + dof2, and in 3D it is equal to dof0 + 3*dof1 + 3*dof2 + dof3
542: Level: developer
544: .seealso: DMSTAG, DMStagGetDOF()
545: @*/
546: PetscErrorCode DMStagGetEntriesPerElement(DM dm,PetscInt *entriesPerElement)
547: {
548: const DM_Stag * const stag = (DM_Stag*)dm->data;
552: if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
553: return(0);
554: }
556: /*@C
557: DMStagGetStencilType - get elementwise ghost/halo stencil type
559: Not Collective
561: Input Parameter:
562: . dm - the DMStag object
564: Output Parameter:
565: . stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE
567: Level: beginner
569: .seealso: DMSTAG, DMStagSetStencilType(), DMStagGetStencilWidth, DMStagStencilType
570: @*/
571: PetscErrorCode DMStagGetStencilType(DM dm,DMStagStencilType *stencilType)
572: {
573: DM_Stag * const stag = (DM_Stag*)dm->data;
577: *stencilType = stag->stencilType;
578: return(0);
579: }
581: /*@C
582: DMStagGetStencilWidth - get elementwise stencil width
584: Not Collective
586: Input Parameter:
587: . dm - the DMStag object
589: Output Parameters:
590: . stencilWidth - stencil/halo/ghost width in elements
592: Level: beginner
594: .seealso: DMSTAG, DMStagSetStencilWidth(), DMStagGetStencilType(), DMDAGetStencilType()
595: @*/
596: PetscErrorCode DMStagGetStencilWidth(DM dm,PetscInt *stencilWidth)
597: {
598: const DM_Stag * const stag = (DM_Stag*)dm->data;
602: if (stencilWidth) *stencilWidth = stag->stencilWidth;
603: return(0);
604: }
606: /*@C
607: DMStagGetOwnershipRanges - get elements per rank in each direction
609: Not Collective
611: Input Parameter:
612: . dm - the DMStag object
614: Output Parameters:
615: + lx - ownership along x direction (optional)
616: . ly - ownership along y direction (optional)
617: - lz - ownership along z direction (optional)
619: Notes:
620: These correspond to the optional final arguments passed to DMStagCreate1d(), DMStagCreate2d(), and DMStagCreate3d().
622: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
624: In C you should not free these arrays, nor change the values in them.
625: They will only have valid values while the DMStag they came from still exists (has not been destroyed).
627: Level: intermediate
629: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagSetOwnershipRanges(), DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDAGetOwnershipRanges()
630: @*/
631: PetscErrorCode DMStagGetOwnershipRanges(DM dm,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
632: {
633: const DM_Stag * const stag = (DM_Stag*)dm->data;
637: if (lx) *lx = stag->l[0];
638: if (ly) *ly = stag->l[1];
639: if (lz) *lz = stag->l[2];
640: return(0);
641: }
643: /*@C
644: DMStagCreateCompatibleDMStag - create a compatible DMStag with different dof/stratum
646: Collective
648: Input Parameters:
649: + dm - the DMStag object
650: - dof0,dof1,dof2,dof3 - number of dof on each stratum in the new DMStag
652: Output Parameters:
653: . newdm - the new, compatible DMStag
655: Notes:
656: Dof supplied for strata too big for the dimension are ignored; these may be set to 0.
657: For example, for a 2-dimensional DMStag, dof2 sets the number of dof per element,
658: and dof3 is unused. For a 3-dimensional DMStag, dof3 sets the number of dof per element.
660: In contrast to DMDACreateCompatibleDMDA(), coordinates are not reused.
662: Level: intermediate
664: .seealso: DMSTAG, DMDACreateCompatibleDMDA(), DMGetCompatibility(), DMStagMigrateVec()
665: @*/
666: PetscErrorCode DMStagCreateCompatibleDMStag(DM dm,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DM *newdm)
667: {
668: PetscErrorCode ierr;
672: DMStagDuplicateWithoutSetup(dm,PetscObjectComm((PetscObject)dm),newdm);
673: DMStagSetDOF(*newdm,dof0,dof1,dof2,dof3);
674: DMSetUp(*newdm);
675: return(0);
676: }
678: /*@C
679: DMStagGetLocationSlot - get index to use in accessing raw local arrays
681: Not Collective
683: Input Parameters:
684: + dm - the DMStag object
685: . loc - location relative to an element
686: - c - component
688: Output Parameter:
689: . slot - index to use
691: Notes:
692: Provides an appropriate index to use with DMStagVecGetArray() and friends.
693: This is required so that the user doesn't need to know about the ordering of
694: dof associated with each local element.
696: Level: beginner
698: .seealso: DMSTAG, DMStagVecGetArray(), DMStagVecGetArrayRead(), DMStagGetDOF(), DMStagGetEntriesPerElement()
699: @*/
700: PetscErrorCode DMStagGetLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt c,PetscInt *slot)
701: {
702: DM_Stag * const stag = (DM_Stag*)dm->data;
706: if (PetscDefined(USE_DEBUG)) {
708: PetscInt dof;
709: DMStagGetLocationDOF(dm,loc,&dof);
710: if (dof < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Location %s has no dof attached",DMStagStencilLocations[loc]);
711: if (c > dof-1) SETERRQ3(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Supplied component number (%D) for location %s is too big (maximum %D)",c,DMStagStencilLocations[loc],dof-1);
712: }
713: *slot = stag->locationOffsets[loc] + c;
714: return(0);
715: }
717: /*@C
718: DMStagMigrateVec - transfer a vector associated with a DMStag to a vector associated with a compatible DMStag
720: Collective
722: Input Parameters:
723: + dm - the source DMStag object
724: . vec - the source vector, compatible with dm
725: . dmTo - the compatible destination DMStag object
726: - vecTo - the destination vector, compatible with dmTo
728: Notes:
729: Extra dof are ignored, and unfilled dof are zeroed.
730: Currently only implemented to migrate global vectors to global vectors.
732: Level: advanced
734: .seealso: DMSTAG, DMStagCreateCompatibleDMStag(), DMGetCompatibility(), DMStagVecSplitToDMDA()
735: @*/
736: PetscErrorCode DMStagMigrateVec(DM dm,Vec vec,DM dmTo,Vec vecTo)
737: {
738: PetscErrorCode ierr;
739: DM_Stag * const stag = (DM_Stag*)dm->data;
740: DM_Stag * const stagTo = (DM_Stag*)dmTo->data;
741: PetscInt nLocalTo,nLocal,dim,i,j,k;
742: PetscInt start[DMSTAG_MAX_DIM],startGhost[DMSTAG_MAX_DIM],n[DMSTAG_MAX_DIM],nExtra[DMSTAG_MAX_DIM],offset[DMSTAG_MAX_DIM];
743: Vec vecToLocal,vecLocal;
744: PetscBool compatible,compatibleSet;
745: const PetscScalar *arr;
746: PetscScalar *arrTo;
747: const PetscInt epe = stag->entriesPerElement;
748: const PetscInt epeTo = stagTo->entriesPerElement;
755: DMGetCompatibility(dm,dmTo,&compatible,&compatibleSet);
756: if (!compatibleSet || !compatible) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"DMStag objects must be shown to be compatible");
757: DMGetDimension(dm,&dim);
758: VecGetLocalSize(vec,&nLocal);
759: VecGetLocalSize(vecTo,&nLocalTo);
760: if (nLocal != stag->entries|| nLocalTo !=stagTo->entries) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Vector migration only implemented for global vector to global vector.");
761: DMStagGetCorners(dm,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
762: DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],NULL,NULL,NULL);
763: for (i=0; i<DMSTAG_MAX_DIM; ++i) offset[i] = start[i]-startGhost[i];
765: /* Proceed by transferring to a local vector, copying, and transferring back to a global vector */
766: DMGetLocalVector(dm,&vecLocal);
767: DMGetLocalVector(dmTo,&vecToLocal);
768: DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal);
769: DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal);
770: VecGetArrayRead(vecLocal,&arr);
771: VecGetArray(vecToLocal,&arrTo);
772: /* Note that some superfluous copying of entries on partial dummy elements is done */
773: if (dim == 1) {
774: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
775: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
776: PetscInt si;
777: for (si=0; si<2; ++si) {
778: b += stag->dof[si];
779: bTo += stagTo->dof[si];
780: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[i*epeTo + dTo] = arr[i*epe + d];
781: for (; dTo < bTo ; ++dTo) arrTo[i*epeTo + dTo] = 0.0;
782: d = b;
783: }
784: }
785: } else if (dim == 2) {
786: const PetscInt epr = stag->nGhost[0] * epe;
787: const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
788: for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
789: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
790: const PetscInt base = j*epr + i*epe;
791: const PetscInt baseTo = j*eprTo + i*epeTo;
792: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
793: const PetscInt s[4] = {0,1,1,2}; /* Dimensions of points, in order */
794: PetscInt si;
795: for (si=0; si<4; ++si) {
796: b += stag->dof[s[si]];
797: bTo += stagTo->dof[s[si]];
798: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
799: for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
800: d = b;
801: }
802: }
803: }
804: } else if (dim == 3) {
805: const PetscInt epr = stag->nGhost[0] * epe;
806: const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
807: const PetscInt epl = stag->nGhost[1] * epr;
808: const PetscInt eplTo = stagTo->nGhost[1] * eprTo;
809: for (k=offset[2]; k<offset[2] + n[2] + nExtra[2]; ++k) {
810: for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
811: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
812: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
813: const PetscInt base = k*epl + j*epr + i*epe;
814: const PetscInt baseTo = k*eplTo + j*eprTo + i*epeTo;
815: const PetscInt s[8] = {0,1,1,2,1,2,2,3}; /* dimensions of points, in order */
816: PetscInt is;
817: for (is=0; is<8; ++is) {
818: b += stag->dof[s[is]];
819: bTo += stagTo->dof[s[is]];
820: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
821: for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
822: d = b;
823: }
824: }
825: }
826: }
827: } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
828: VecRestoreArrayRead(vecLocal,&arr);
829: VecRestoreArray(vecToLocal,&arrTo);
830: DMRestoreLocalVector(dm,&vecLocal);
831: DMLocalToGlobalBegin(dmTo,vecToLocal,INSERT_VALUES,vecTo);
832: DMLocalToGlobalEnd(dmTo,vecToLocal,INSERT_VALUES,vecTo);
833: DMRestoreLocalVector(dmTo,&vecToLocal);
834: return(0);
835: }
837: /*@C
838: DMStagPopulateLocalToGlobalInjective - populate an internal 1-to-1 local-to-global map
840: Collective
842: Creates an internal object which explicitly maps a single local degree of
843: freedom to each global degree of freedom. This is used, if populated,
844: instead of SCATTER_REVERSE_LOCAL with the (1-to-many, in general)
845: global-to-local map, when DMLocalToGlobal() is called with INSERT_VALUES.
846: This allows usage, for example, even in the periodic, 1-rank case, where
847: the inverse of the global-to-local map, even when restricted to on-rank
848: communication, is non-injective. This is at the cost of storing an additional
849: VecScatter object inside each DMStag object.
851: Input Parameter:
852: . dm - the DMStag object
854: Notes:
855: In normal usage, library users shouldn't be concerned with this function,
856: as it is called during DMSetUp(), when required.
858: Returns immediately if the internal map is already populated.
860: Developer Notes:
861: This could, if desired, be moved up to a general DM routine. It would allow,
862: for example, DMDA to support DMLocalToGlobal() with INSERT_VALUES,
863: even in the single-rank periodic case.
865: Level: developer
867: .seealso: DMSTAG, DMLocalToGlobal(), VecScatter
868: @*/
869: PetscErrorCode DMStagPopulateLocalToGlobalInjective(DM dm)
870: {
871: PetscErrorCode ierr;
872: PetscInt dim;
873: DM_Stag * const stag = (DM_Stag*)dm->data;
877: if (stag->ltog_injective) return(0); /* Don't re-populate */
878: DMGetDimension(dm,&dim);
879: switch (dim) {
880: case 1: DMStagPopulateLocalToGlobalInjective_1d(dm); break;
881: case 2: DMStagPopulateLocalToGlobalInjective_2d(dm); break;
882: case 3: DMStagPopulateLocalToGlobalInjective_3d(dm); break;
883: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
884: }
885: return(0);
886: }
888: static PetscErrorCode DMStagRestoreProductCoordinateArrays_Private(DM dm,void *arrX,void *arrY,void *arrZ,PetscBool read)
889: {
890: PetscErrorCode ierr;
891: PetscInt dim,d;
892: void* arr[DMSTAG_MAX_DIM];
893: DM dmCoord;
897: DMGetDimension(dm,&dim);
898: if (dim > DMSTAG_MAX_DIM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for %D dimensions",dim);
899: arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
900: DMGetCoordinateDM(dm,&dmCoord);
901: for (d=0; d<dim; ++d) {
902: DM subDM;
903: Vec coord1d_local;
905: /* Ignore unrequested arrays */
906: if (!arr[d]) continue;
908: DMProductGetDM(dmCoord,d,&subDM);
909: DMGetCoordinatesLocal(subDM,&coord1d_local);
910: if (read) {
911: DMStagVecRestoreArrayRead(subDM,coord1d_local,arr[d]);
912: } else {
913: DMStagVecRestoreArray(subDM,coord1d_local,arr[d]);
914: }
915: }
916: return(0);
917: }
919: /*@C
920: DMStagRestoreProductCoordinateArrays - restore local array access
922: Logically Collective
924: Input Parameter:
925: . dm - the DMStag object
927: Output Parameters:
928: . arrX,arrY,arrZ - local 1D coordinate arrays
930: Level: intermediate
932: Notes:
933: 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:
935: $ DMGetCoordinateDM(dm,&cdm);
936: $ for (d=0; d<3; ++d) {
937: $ DM subdm;
938: $ Vec coor,coor_local;
940: $ DMProductGetDM(cdm,d,&subdm);
941: $ DMGetCoordinates(subdm,&coor);
942: $ DMGetCoordinatesLocal(subdm,&coor_local);
943: $ DMLocalToGlobal(subdm,coor_local,INSERT_VALUES,coor);
944: $ PetscPrintf(PETSC_COMM_WORLD,"Coordinates dim %D:\n",d);
945: $ VecView(coor,PETSC_VIEWER_STDOUT_WORLD);
946: $ }
948: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
949: @*/
950: PetscErrorCode DMStagRestoreProductCoordinateArrays(DM dm,void *arrX,void *arrY,void *arrZ)
951: {
952: PetscErrorCode ierr;
955: DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_FALSE);
956: return(0);
957: }
959: /*@C
960: DMStagRestoreProductCoordinateArraysRead - restore local product array access, read-only
962: Logically Collective
964: Input Parameter:
965: . dm - the DMStag object
967: Output Parameters:
968: . arrX,arrY,arrZ - local 1D coordinate arrays
970: Level: intermediate
972: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
973: @*/
974: PetscErrorCode DMStagRestoreProductCoordinateArraysRead(DM dm,void *arrX,void *arrY,void *arrZ)
975: {
976: PetscErrorCode ierr;
979: DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_TRUE);
980: return(0);
981: }
983: /*@C
984: DMStagSetBoundaryTypes - set DMStag boundary types
986: Logically Collective; boundaryType0, boundaryType1, and boundaryType2 must contain common values
988: Input Parameters:
989: + dm - the DMStag object
990: - boundaryType0,boundaryType1,boundaryType2 - boundary types in each direction
992: Notes:
993: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
995: Level: advanced
997: .seealso: DMSTAG, DMBoundaryType, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDASetBoundaryType()
998: @*/
999: PetscErrorCode DMStagSetBoundaryTypes(DM dm,DMBoundaryType boundaryType0,DMBoundaryType boundaryType1,DMBoundaryType boundaryType2)
1000: {
1001: PetscErrorCode ierr;
1002: DM_Stag * const stag = (DM_Stag*)dm->data;
1003: PetscInt dim;
1006: DMGetDimension(dm,&dim);
1011: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1012: stag->boundaryType[0] = boundaryType0;
1013: if (dim > 1) stag->boundaryType[1] = boundaryType1;
1014: if (dim > 2) stag->boundaryType[2] = boundaryType2;
1015: return(0);
1016: }
1018: /*@C
1019: DMStagSetCoordinateDMType - set DM type to store coordinates
1021: Logically Collective; dmtype must contain common value
1023: Input Parameters:
1024: + dm - the DMStag object
1025: - dmtype - DMtype for coordinates, either DMSTAG or DMPRODUCT
1027: Level: advanced
1029: .seealso: DMSTAG, DMPRODUCT, DMGetCoordinateDM(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMType
1030: @*/
1031: PetscErrorCode DMStagSetCoordinateDMType(DM dm,DMType dmtype)
1032: {
1033: PetscErrorCode ierr;
1034: 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: PetscErrorCode ierr;
1062: DM_Stag * const stag = (DM_Stag*)dm->data;
1063: PetscInt dim;
1071: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1072: DMGetDimension(dm,&dim);
1073: if (dof0 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof0 cannot be negative");
1074: if (dof1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof1 cannot be negative");
1075: if (dim > 1 && dof2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof2 cannot be negative");
1076: if (dim > 2 && dof3 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof3 cannot be negative");
1077: stag->dof[0] = dof0;
1078: stag->dof[1] = dof1;
1079: if (dim > 1) stag->dof[2] = dof2;
1080: if (dim > 2) stag->dof[3] = dof3;
1081: return(0);
1082: }
1084: /*@C
1085: DMStagSetNumRanks - set ranks in each direction in the global rank grid
1087: Logically Collective; nRanks0, nRanks1, and nRanks2 must contain common values
1089: Input Parameters:
1090: + dm - the DMStag object
1091: - nRanks0,nRanks1,nRanks2 - number of ranks in each direction
1093: Notes:
1094: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1096: Level: developer
1098: .seealso: DMSTAG, DMDASetNumProcs()
1099: @*/
1100: PetscErrorCode DMStagSetNumRanks(DM dm,PetscInt nRanks0,PetscInt nRanks1,PetscInt nRanks2)
1101: {
1102: PetscErrorCode ierr;
1103: DM_Stag * const stag = (DM_Stag*)dm->data;
1104: PetscInt dim;
1111: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1112: DMGetDimension(dm,&dim);
1113: if (nRanks0 != PETSC_DECIDE && nRanks0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in X direction cannot be less than 1");
1114: if (dim > 1 && nRanks1 != PETSC_DECIDE && nRanks1 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in Y direction cannot be less than 1");
1115: if (dim > 2 && nRanks2 != PETSC_DECIDE && nRanks2 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in Z direction cannot be less than 1");
1116: if (nRanks0) stag->nRanks[0] = nRanks0;
1117: if (dim > 1 && nRanks1) stag->nRanks[1] = nRanks1;
1118: if (dim > 2 && nRanks2) stag->nRanks[2] = nRanks2;
1119: return(0);
1120: }
1122: /*@C
1123: DMStagSetStencilType - set elementwise ghost/halo stencil type
1125: Logically Collective; stencilType must contain common value
1127: Input Parameters:
1128: + dm - the DMStag object
1129: - stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE
1131: Level: beginner
1133: .seealso: DMSTAG, DMStagGetStencilType(), DMStagSetStencilWidth(), DMStagStencilType
1134: @*/
1135: PetscErrorCode DMStagSetStencilType(DM dm,DMStagStencilType stencilType)
1136: {
1137: DM_Stag * const stag = (DM_Stag*)dm->data;
1142: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1143: stag->stencilType = stencilType;
1144: return(0);
1145: }
1147: /*@C
1148: DMStagSetStencilWidth - set elementwise stencil width
1150: Logically Collective; stencilWidth must contain common value
1152: Input Parameters:
1153: + dm - the DMStag object
1154: - stencilWidth - stencil/halo/ghost width in elements
1156: Notes:
1157: The width value is not used when DMSTAG_STENCIL_NONE is specified.
1159: Level: beginner
1161: .seealso: DMSTAG, DMStagGetStencilWidth(), DMStagGetStencilType(), DMStagStencilType
1162: @*/
1163: PetscErrorCode DMStagSetStencilWidth(DM dm,PetscInt stencilWidth)
1164: {
1165: DM_Stag * const stag = (DM_Stag*)dm->data;
1170: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1171: if (stencilWidth < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Stencil width must be non-negative");
1172: stag->stencilWidth = stencilWidth;
1173: return(0);
1174: }
1176: /*@C
1177: DMStagSetGlobalSizes - set global element counts in each direction
1179: Logically Collective; N0, N1, and N2 must contain common values
1181: Input Parameters:
1182: + dm - the DMStag object
1183: - N0,N1,N2 - global elementwise sizes
1185: Notes:
1186: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1188: Level: advanced
1190: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMDASetSizes()
1191: @*/
1192: PetscErrorCode DMStagSetGlobalSizes(DM dm,PetscInt N0,PetscInt N1,PetscInt N2)
1193: {
1194: PetscErrorCode ierr;
1195: DM_Stag * const stag = (DM_Stag*)dm->data;
1196: PetscInt dim;
1200: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1201: DMGetDimension(dm,&dim);
1202: if (N0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in X direction must be positive");
1203: if (dim > 1 && N1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Y direction must be positive");
1204: if (dim > 2 && N2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Z direction must be positive");
1205: if (N0) stag->N[0] = N0;
1206: if (N1) stag->N[1] = N1;
1207: if (N2) stag->N[2] = N2;
1208: return(0);
1209: }
1211: /*@C
1212: DMStagSetOwnershipRanges - set elements per rank in each direction
1214: Logically Collective; lx, ly, and lz must contain common values
1216: Input Parameters:
1217: + dm - the DMStag object
1218: - lx,ly,lz - element counts for each rank in each direction
1220: Notes:
1221: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
1223: Level: developer
1225: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagGetOwnershipRanges(), DMDASetOwnershipRanges()
1226: @*/
1227: PetscErrorCode DMStagSetOwnershipRanges(DM dm,PetscInt const *lx,PetscInt const *ly,PetscInt const *lz)
1228: {
1229: PetscErrorCode ierr;
1230: DM_Stag * const stag = (DM_Stag*)dm->data;
1231: const PetscInt *lin[3];
1232: PetscInt d,dim;
1236: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1237: lin[0] = lx; lin[1] = ly; lin[2] = lz;
1238: DMGetDimension(dm,&dim);
1239: for (d=0; d<dim; ++d) {
1240: if (lin[d]) {
1241: if (stag->nRanks[d] < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of ranks");
1242: if (!stag->l[d]) {
1243: PetscMalloc1(stag->nRanks[d], &stag->l[d]);
1244: }
1245: PetscArraycpy(stag->l[d], lin[d], stag->nRanks[d]);
1246: }
1247: }
1248: return(0);
1249: }
1251: /*@C
1252: DMStagSetUniformCoordinates - set DMStag coordinates to be a uniform grid
1254: Collective
1256: Input Parameters:
1257: + dm - the DMStag object
1258: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1260: Notes:
1261: DMStag supports 2 different types of coordinate DM: DMSTAG and DMPRODUCT.
1262: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1264: Local coordinates are populated (using DMSetCoordinatesLocal()), linearly
1265: extrapolated to ghost cells, including those outside the physical domain.
1266: This is also done in case of periodic boundaries, meaning that the same
1267: global point may have different coordinates in different local representations,
1268: which are equivalent assuming a periodicity implied by the arguments to this function,
1269: i.e. two points are equivalent if their difference is a multiple of (xmax - xmin)
1270: in the x direction, (ymax - ymin) in the y direction, and (zmax - zmin) in the z direction.
1272: Level: advanced
1274: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType(), DMGetCoordinateDM(), DMGetCoordinates(), DMDASetUniformCoordinates(), DMBoundaryType
1275: @*/
1276: PetscErrorCode DMStagSetUniformCoordinates(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1277: {
1278: PetscErrorCode ierr;
1279: DM_Stag * const stag = (DM_Stag*)dm->data;
1280: PetscBool flg_stag,flg_product;
1284: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1285: if (!stag->coordinateDMType) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"You must first call DMStagSetCoordinateDMType()");
1286: PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg_stag);
1287: PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg_product);
1288: if (flg_stag) {
1289: DMStagSetUniformCoordinatesExplicit(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1290: } else if (flg_product) {
1291: DMStagSetUniformCoordinatesProduct(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1292: } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported DM Type %s",stag->coordinateDMType);
1293: return(0);
1294: }
1296: /*@C
1297: DMStagSetUniformCoordinatesExplicit - set DMStag coordinates to be a uniform grid, storing all values
1299: Collective
1301: Input Parameters:
1302: + dm - the DMStag object
1303: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1305: Notes:
1306: DMStag supports 2 different types of coordinate DM: either another DMStag, or a DMProduct.
1307: If the grid is orthogonal, using DMProduct should be more efficient.
1309: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1311: See the manual page for DMStagSetUniformCoordinates() for information on how
1312: coordinates for dummy cells outside the physical domain boundary are populated.
1314: Level: beginner
1316: .seealso: DMSTAG, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType()
1317: @*/
1318: PetscErrorCode DMStagSetUniformCoordinatesExplicit(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1319: {
1320: PetscErrorCode ierr;
1321: DM_Stag * const stag = (DM_Stag*)dm->data;
1322: PetscInt dim;
1323: PetscBool flg;
1327: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1328: PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg);
1329: if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1330: DMStagSetCoordinateDMType(dm,DMSTAG);
1331: DMGetDimension(dm,&dim);
1332: switch (dim) {
1333: case 1: DMStagSetUniformCoordinatesExplicit_1d(dm,xmin,xmax); break;
1334: case 2: DMStagSetUniformCoordinatesExplicit_2d(dm,xmin,xmax,ymin,ymax); break;
1335: case 3: DMStagSetUniformCoordinatesExplicit_3d(dm,xmin,xmax,ymin,ymax,zmin,zmax); break;
1336: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1337: }
1338: return(0);
1339: }
1341: /*@C
1342: DMStagSetUniformCoordinatesProduct - create uniform coordinates, as a product of 1D arrays
1344: 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.
1346: Collective
1348: Input Parameters:
1349: + dm - the DMStag object
1350: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1352: Notes:
1353: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1355: The per-dimension 1-dimensional DMStag objects that comprise the product
1356: always have active 0-cells (vertices, element boundaries) and 1-cells
1357: (element centers).
1359: See the manual page for DMStagSetUniformCoordinates() for information on how
1360: coordinates for dummy cells outside the physical domain boundary are populated.
1362: Level: intermediate
1364: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetCoordinateDMType()
1365: @*/
1366: PetscErrorCode DMStagSetUniformCoordinatesProduct(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1367: {
1368: PetscErrorCode ierr;
1369: DM_Stag * const stag = (DM_Stag*)dm->data;
1370: DM dmc;
1371: PetscInt dim,d,dof0,dof1;
1372: PetscBool flg;
1376: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1377: PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg);
1378: if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1379: DMStagSetCoordinateDMType(dm,DMPRODUCT);
1380: DMGetCoordinateDM(dm,&dmc);
1381: DMGetDimension(dm,&dim);
1383: /* Create 1D sub-DMs, living on subcommunicators.
1384: Always include both vertex and element dof, regardless of the active strata of the DMStag */
1385: dof0 = 1;
1386: dof1 = 1;
1388: for (d=0; d<dim; ++d) {
1389: DM subdm;
1390: MPI_Comm subcomm;
1391: PetscMPIInt color;
1392: const PetscMPIInt key = 0; /* let existing rank break ties */
1394: /* Choose colors based on position in the plane orthogonal to this dim, and split */
1395: switch (d) {
1396: case 0: color = (dim > 1 ? stag->rank[1] : 0) + (dim > 2 ? stag->nRanks[1]*stag->rank[2] : 0); break;
1397: case 1: color = stag->rank[0] + (dim > 2 ? stag->nRanks[0]*stag->rank[2] : 0); break;
1398: case 2: color = stag->rank[0] + stag->nRanks[0]*stag->rank[1] ; break;
1399: default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1400: }
1401: MPI_Comm_split(PetscObjectComm((PetscObject)dm),color,key,&subcomm);
1403: /* Create sub-DMs living on these new communicators (which are destroyed by DMProduct) */
1404: DMStagCreate1d(subcomm,stag->boundaryType[d],stag->N[d],dof0,dof1,stag->stencilType,stag->stencilWidth,stag->l[d],&subdm);
1405: DMSetUp(subdm);
1406: switch (d) {
1407: case 0:
1408: DMStagSetUniformCoordinatesExplicit(subdm,xmin,xmax,0,0,0,0);
1409: break;
1410: case 1:
1411: DMStagSetUniformCoordinatesExplicit(subdm,ymin,ymax,0,0,0,0);
1412: break;
1413: case 2:
1414: DMStagSetUniformCoordinatesExplicit(subdm,zmin,zmax,0,0,0,0);
1415: break;
1416: default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1417: }
1418: DMProductSetDM(dmc,d,subdm);
1419: DMProductSetDimensionIndex(dmc,d,0);
1420: DMDestroy(&subdm);
1421: MPI_Comm_free(&subcomm);
1422: }
1423: return(0);
1424: }
1426: /*@C
1427: DMStagVecGetArray - get access to local array
1429: Logically Collective
1431: This function returns a (dim+1)-dimensional array for a dim-dimensional
1432: DMStag.
1434: The first 1-3 dimensions indicate an element in the global
1435: numbering, using the standard C ordering.
1437: The final dimension in this array corresponds to a degree
1438: of freedom with respect to this element, for example corresponding to
1439: the element or one of its neighboring faces, edges, or vertices.
1441: For example, for a 3D DMStag, indexing is array[k][j][i][idx], where k is the
1442: index in the z-direction, j is the index in the y-direction, and i is the
1443: index in the x-direction.
1445: "idx" is obtained with DMStagGetLocationSlot(), since the correct offset
1446: into the (dim+1)-dimensional C array depends on the grid size and the number
1447: of dof stored at each location.
1449: Input Parameters:
1450: + dm - the DMStag object
1451: - vec - the Vec object
1453: Output Parameters:
1454: . array - the array
1456: Notes:
1457: DMStagVecRestoreArray() must be called, once finished with the array
1459: Level: beginner
1461: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArray(), DMDAVecGetArrayDOF()
1462: @*/
1463: PetscErrorCode DMStagVecGetArray(DM dm,Vec vec,void *array)
1464: {
1465: PetscErrorCode ierr;
1466: DM_Stag * const stag = (DM_Stag*)dm->data;
1467: PetscInt dim;
1468: PetscInt nLocal;
1473: DMGetDimension(dm,&dim);
1474: VecGetLocalSize(vec,&nLocal);
1475: if (nLocal != stag->entriesGhost) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMStag local size %D\n",nLocal,stag->entriesGhost);
1476: switch (dim) {
1477: case 1:
1478: VecGetArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1479: break;
1480: case 2:
1481: VecGetArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1482: break;
1483: case 3:
1484: 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);
1485: break;
1486: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1487: }
1488: return(0);
1489: }
1491: /*@C
1492: DMStagVecGetArrayRead - get read-only access to a local array
1494: Logically Collective
1496: See the man page for DMStagVecGetArray() for more information.
1498: Input Parameters:
1499: + dm - the DMStag object
1500: - vec - the Vec object
1502: Output Parameters:
1503: . array - the read-only array
1505: Notes:
1506: DMStagVecRestoreArrayRead() must be called, once finished with the array
1508: Level: beginner
1510: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArrayRead(), DMDAVecGetArrayDOFRead()
1511: @*/
1512: PetscErrorCode DMStagVecGetArrayRead(DM dm,Vec vec,void *array)
1513: {
1514: PetscErrorCode ierr;
1515: DM_Stag * const stag = (DM_Stag*)dm->data;
1516: PetscInt dim;
1517: PetscInt nLocal;
1522: DMGetDimension(dm,&dim);
1523: VecGetLocalSize(vec,&nLocal);
1524: if (nLocal != stag->entriesGhost) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMStag local size %D\n",nLocal,stag->entriesGhost);
1525: switch (dim) {
1526: case 1:
1527: VecGetArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1528: break;
1529: case 2:
1530: VecGetArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1531: break;
1532: case 3:
1533: 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);
1534: break;
1535: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1536: }
1537: return(0);
1538: }
1540: /*@C
1541: DMStagVecRestoreArray - restore access to a raw array
1543: Logically Collective
1545: Input Parameters:
1546: + dm - the DMStag object
1547: - vec - the Vec object
1549: Output Parameters:
1550: . array - the array
1552: Level: beginner
1554: .seealso: DMSTAG, DMStagVecGetArray(), DMDAVecRestoreArray(), DMDAVecRestoreArrayDOF()
1555: @*/
1556: PetscErrorCode DMStagVecRestoreArray(DM dm,Vec vec,void *array)
1557: {
1558: PetscErrorCode ierr;
1559: DM_Stag * const stag = (DM_Stag*)dm->data;
1560: PetscInt dim;
1561: PetscInt nLocal;
1566: DMGetDimension(dm,&dim);
1567: VecGetLocalSize(vec,&nLocal);
1568: if (nLocal != stag->entriesGhost) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMStag local size %D\n",nLocal,stag->entriesGhost);
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: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1580: }
1581: return(0);
1582: }
1584: /*@C
1585: DMStagVecRestoreArrayRead - restore read-only access to a raw array
1587: Logically Collective
1589: Input Parameters:
1590: + dm - the DMStag object
1591: - vec - the Vec object
1593: Output Parameters:
1594: . array - the read-only array
1596: Level: beginner
1598: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMDAVecRestoreArrayRead(), DMDAVecRestoreArrayDOFRead()
1599: @*/
1600: PetscErrorCode DMStagVecRestoreArrayRead(DM dm,Vec vec,void *array)
1601: {
1602: PetscErrorCode ierr;
1603: DM_Stag * const stag = (DM_Stag*)dm->data;
1604: PetscInt dim;
1605: PetscInt nLocal;
1610: DMGetDimension(dm,&dim);
1611: VecGetLocalSize(vec,&nLocal);
1612: if (nLocal != stag->entriesGhost) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMStag local size %D\n",nLocal,stag->entriesGhost);
1613: switch (dim) {
1614: case 1:
1615: VecRestoreArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1616: break;
1617: case 2:
1618: VecRestoreArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1619: break;
1620: case 3:
1621: 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);
1622: break;
1623: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1624: }
1625: return(0);
1626: }