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 DMStagGetProductCoordinateSlot() 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: Checks that the coordinates are actually set up so that using the
173: slots from the first 1d coordinate sub-DM is valid for all the 1D coordinate sub-DMs.
175: Level: intermediate
177: .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead(), DMStagSetUniformCoordinates()
178: @*/
179: PETSC_EXTERN PetscErrorCode DMStagGetProductCoordinateLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt *slot)
180: {
182: DM dmCoord;
183: PetscInt dim,dofCheck[DMSTAG_MAX_STRATA],s,d;
187: DMGetDimension(dm,&dim);
188: DMGetCoordinateDM(dm,&dmCoord);
189: if (!dmCoord) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM does not have a coordinate DM");
190: {
191: PetscBool isProduct;
192: DMType dmType;
193: DMGetType(dmCoord,&dmType);
194: PetscStrcmp(DMPRODUCT,dmType,&isProduct);
195: if (!isProduct) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is not of type DMPRODUCT");
196: }
197: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
198: for (d=0; d<dim; ++d) {
199: DM subDM;
200: DMType dmType;
201: PetscBool isStag;
202: PetscInt dof[DMSTAG_MAX_STRATA],subDim;
203: DMProductGetDM(dmCoord,d,&subDM);
204: if (!subDM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is missing sub DM %D",d);
205: DMGetDimension(subDM,&subDim);
206: if (subDim != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of dimension 1");
207: DMGetType(subDM,&dmType);
208: PetscStrcmp(DMSTAG,dmType,&isStag);
209: if (!isStag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of type DMSTAG");
210: DMStagGetDOF(subDM,&dof[0],&dof[1],&dof[2],&dof[3]);
211: if (d == 0) {
212: const PetscInt component = 0;
213: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
214: DMStagGetLocationSlot(subDM,loc,component,slot);
215: } else {
216: for (s=0; s<DMSTAG_MAX_STRATA; ++s) {
217: if (dofCheck[s] != dof[s]) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DMs have different dofs");
218: }
219: }
220: }
221: return(0);
222: }
224: /*@C
225: DMStagGetCorners - return global element indices of the local region (excluding ghost points)
227: Not Collective
229: Input Parameter:
230: . dm - the DMStag object
232: Output Parameters:
233: + x,y,z - starting element indices in each direction
234: . m,n,p - element widths in each direction
235: - nExtrax,nExtray,nExtraz - number of extra partial elements in each direction.
237: Notes:
238: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
240: The number of extra partial elements is either 1 or 0.
241: The value is 1 on right, top, and front non-periodic domain ("physical") boundaries,
242: in the x, y, and z directions respectively, and otherwise 0.
244: Level: beginner
246: .seealso: DMSTAG, DMStagGetGhostCorners(), DMDAGetCorners()
247: @*/
248: PetscErrorCode DMStagGetCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *nExtrax,PetscInt *nExtray,PetscInt *nExtraz)
249: {
250: const DM_Stag * const stag = (DM_Stag*)dm->data;
254: if (x) *x = stag->start[0];
255: if (y) *y = stag->start[1];
256: if (z) *z = stag->start[2];
257: if (m) *m = stag->n[0];
258: if (n) *n = stag->n[1];
259: if (p) *p = stag->n[2];
260: if (nExtrax) *nExtrax = stag->boundaryType[0] != DM_BOUNDARY_PERIODIC && stag->lastRank[0] ? 1 : 0;
261: if (nExtray) *nExtray = stag->boundaryType[1] != DM_BOUNDARY_PERIODIC && stag->lastRank[1] ? 1 : 0;
262: if (nExtraz) *nExtraz = stag->boundaryType[2] != DM_BOUNDARY_PERIODIC && stag->lastRank[2] ? 1 : 0;
263: return(0);
264: }
266: /*@C
267: DMStagGetDOF - get number of DOF associated with each stratum of the grid
269: Not Collective
271: Input Parameter:
272: . dm - the DMStag object
274: Output Parameters:
275: + dof0 - the number of points per 0-cell (vertex/node)
276: . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
277: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
278: - dof3 - the number of points per 3-cell (element in 3D)
280: Level: beginner
282: .seealso: DMSTAG, DMStagGetCorners(), DMStagGetGhostCorners(), DMStagGetGlobalSizes(), DMStagGetStencilWidth(), DMStagGetBoundaryTypes(), DMStagGetLocationDof(), DMDAGetDof()
283: @*/
284: PetscErrorCode DMStagGetDOF(DM dm,PetscInt *dof0,PetscInt *dof1,PetscInt *dof2,PetscInt *dof3)
285: {
286: const DM_Stag * const stag = (DM_Stag*)dm->data;
290: if (dof0) *dof0 = stag->dof[0];
291: if (dof1) *dof1 = stag->dof[1];
292: if (dof2) *dof2 = stag->dof[2];
293: if (dof3) *dof3 = stag->dof[3];
294: return(0);
295: }
297: /*@C
298: DMStagGetGhostCorners - return global element indices of the local region, including ghost points
300: Not Collective
302: Input Argument:
303: . dm - the DMStag object
305: Output Arguments:
306: + x,y,z - starting element indices in each direction
307: - m,n,p - element widths in each direction
309: Notes:
310: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
312: Level: beginner
314: .seealso: DMSTAG, DMStagGetCorners(), DMDAGetGhostCorners()
315: @*/
316: PetscErrorCode DMStagGetGhostCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
317: {
318: 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;
355: if (M) *M = stag->N[0];
356: if (N) *N = stag->N[1];
357: if (P) *P = stag->N[2];
358: return(0);
359: }
361: /*@C
362: DMStagGetIsFirstRank - get boolean value for whether this rank is first in each direction in the rank grid
364: Not Collective
366: Input Parameter:
367: . dm - the DMStag object
369: Output Parameters:
370: . isFirstRank0,isFirstRank1,isFirstRank2 - whether this rank is first in each direction
372: Level: intermediate
374: Notes:
375: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
377: .seealso: DMSTAG, DMStagGetIsLastRank()
378: @*/
379: PetscErrorCode DMStagGetIsFirstRank(DM dm,PetscBool *isFirstRank0,PetscBool *isFirstRank1,PetscBool *isFirstRank2)
380: {
381: const DM_Stag * const stag = (DM_Stag*)dm->data;
385: if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
386: if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
387: if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
388: return(0);
389: }
391: /*@C
392: DMStagGetIsLastRank - get boolean value for whether this rank is last in each direction in the rank grid
394: Not Collective
396: Input Parameter:
397: . dm - the DMStag object
399: Output Parameters:
400: . isLastRank0,isLastRank1,isLastRank2 - whether this rank is last in each direction
402: Level: intermediate
404: Notes:
405: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
406: Level: intermediate
408: .seealso: DMSTAG, DMStagGetIsFirstRank()
409: @*/
410: PetscErrorCode DMStagGetIsLastRank(DM dm,PetscBool *isLastRank0,PetscBool *isLastRank1,PetscBool *isLastRank2)
411: {
412: const DM_Stag * const stag = (DM_Stag*)dm->data;
416: if (isLastRank0) *isLastRank0 = stag->lastRank[0];
417: if (isLastRank1) *isLastRank1 = stag->lastRank[1];
418: if (isLastRank2) *isLastRank2 = stag->lastRank[2];
419: return(0);
420: }
422: /*@C
423: DMStagGetLocalSizes - get local elementwise sizes
425: Not Collective
427: Input Parameter:
428: . dm - the DMStag object
430: Output Parameters:
431: . m,n,p - local element counts (excluding ghosts) in each direction
433: Notes:
434: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
435: Level: intermediate
437: Level: beginner
439: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetDOF(), DMStagGetNumRanks(), DMDAGetLocalInfo()
440: @*/
441: PetscErrorCode DMStagGetLocalSizes(DM dm,PetscInt* m,PetscInt* n,PetscInt* p)
442: {
443: const DM_Stag * const stag = (DM_Stag*)dm->data;
447: if (m) *m = stag->n[0];
448: if (n) *n = stag->n[1];
449: if (p) *p = stag->n[2];
450: return(0);
451: }
453: /*@C
454: DMStagGetNumRanks - get number of ranks in each direction in the global grid decomposition
456: Not Collective
458: Input Parameter:
459: . dm - the DMStag object
461: Output Parameters:
462: . nRanks0,nRanks1,nRanks2 - number of ranks in each direction in the grid decomposition
464: Notes:
465: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
466: Level: intermediate
468: Level: beginner
470: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetLocalSize(), DMStagSetNumRanks(), DMDAGetInfo()
471: @*/
472: PetscErrorCode DMStagGetNumRanks(DM dm,PetscInt *nRanks0,PetscInt *nRanks1,PetscInt *nRanks2)
473: {
474: const DM_Stag * const stag = (DM_Stag*)dm->data;
478: if (nRanks0) *nRanks0 = stag->nRanks[0];
479: if (nRanks1) *nRanks1 = stag->nRanks[1];
480: if (nRanks2) *nRanks2 = stag->nRanks[2];
481: return(0);
482: }
484: /*@C
485: DMStagGetEntries - get number of native entries in the global representation
487: Not Collective
489: Input Parameter:
490: . dm - the DMStag object
492: Output Parameters:
493: . entries - number of rank-native entries in the global representation
495: Note:
496: This is the number of entries on this rank for a global vector associated with dm.
498: Level: developer
500: .seealso: DMSTAG, DMStagGetDOF(), DMStagGetEntriesPerElement(), DMCreateLocalVector()
501: @*/
502: PetscErrorCode DMStagGetEntries(DM dm,PetscInt *entries)
503: {
504: const DM_Stag * const stag = (DM_Stag*)dm->data;
508: if (entries) *entries = stag->entries;
509: return(0);
510: }
512: /*@C
513: DMStagGetEntriesPerElement - get number of entries per element in the local representation
515: Not Collective
517: Input Parameter:
518: . dm - the DMStag object
520: Output Parameters:
521: . entriesPerElement - number of entries associated with each element in the local representation
523: Notes:
524: This is the natural block size for most local operations. In 1D it is equal to dof0 + dof1,
525: in 2D it is equal to dof0 + 2*dof1 + dof2, and in 3D it is equal to dof0 + 3*dof1 + 3*dof2 + dof3
527: Level: developer
529: .seealso: DMSTAG, DMStagGetDOF()
530: @*/
531: PetscErrorCode DMStagGetEntriesPerElement(DM dm,PetscInt *entriesPerElement)
532: {
533: const DM_Stag * const stag = (DM_Stag*)dm->data;
537: if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
538: return(0);
539: }
541: /*@C
542: DMStagGetStencilType - get elementwise ghost/halo stencil type
544: Not Collective
546: Input Parameter:
547: . dm - the DMStag object
549: Output Parameter:
550: . stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE
552: Level: beginner
554: .seealso: DMSTAG, DMStagSetStencilType(), DMStagStencilType, DMDAGetInfo()
555: @*/
556: PetscErrorCode DMStagGetStencilType(DM dm,DMStagStencilType *stencilType)
557: {
558: DM_Stag * const stag = (DM_Stag*)dm->data;
562: *stencilType = stag->stencilType;
563: return(0);
564: }
566: /*@C
567: DMStagGetStencilWidth - get elementwise stencil width
569: Not Collective
571: Input Parameter:
572: . dm - the DMStag object
574: Output Parameters:
575: . stencilWidth - stencil/halo/ghost width in elements
577: Level: beginner
579: .seealso: DMSTAG, DMDAGetStencilWidth(), DMDAGetInfo()
580: @*/
581: PetscErrorCode DMStagGetStencilWidth(DM dm,PetscInt *stencilWidth)
582: {
583: const DM_Stag * const stag = (DM_Stag*)dm->data;
587: if (stencilWidth) *stencilWidth = stag->stencilWidth;
588: return(0);
589: }
591: /*@C
592: DMStagGetOwnershipRanges - get elements per rank in each direction
594: Not Collective
596: Input Parameter:
597: . dm - the DMStag object
599: Output Parameters:
600: + lx - ownership along x direction (optional)
601: . ly - ownership along y direction (optional)
602: - lz - ownership along z direction (optional)
604: Notes:
605: These correspond to the optional final arguments passed to DMStagCreate1d(), DMStagCreate2d(), and DMStagCreate3d().
607: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
609: In C you should not free these arrays, nor change the values in them.
610: They will only have valid values while the DMStag they came from still exists (has not been destroyed).
612: Level: intermediate
614: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagSetOwnershipRanges(), DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDAGetOwnershipRanges()
615: @*/
616: PetscErrorCode DMStagGetOwnershipRanges(DM dm,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
617: {
618: const DM_Stag * const stag = (DM_Stag*)dm->data;
622: if (lx) *lx = stag->l[0];
623: if (ly) *ly = stag->l[1];
624: if (lz) *lz = stag->l[2];
625: return(0);
626: }
628: /*@C
629: DMStagCreateCompatibleDMStag - create a compatible DMStag with different dof/stratum
631: Collective
633: Input Parameters:
634: + dm - the DMStag object
635: - dof0,dof1,dof2,dof3 - number of dof on each stratum in the new DMStag
637: Output Parameters:
638: . newdm - the new, compatible DMStag
640: Notes:
641: Dof supplied for strata too big for the dimension are ignored; these may be set to 0.
642: For example, for a 2-dimensional DMStag, dof2 sets the number of dof per element,
643: and dof3 is unused. For a 3-dimensional DMStag, dof3 sets the number of dof per element.
645: In contrast to DMDACreateCompatibleDMDA(), coordinates are not reused.
647: Level: intermediate
649: .seealso: DMSTAG, DMDACreateCompatibleDMDA(), DMGetCompatibility(), DMStagMigrateVec()
650: @*/
651: PetscErrorCode DMStagCreateCompatibleDMStag(DM dm,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DM *newdm)
652: {
653: PetscErrorCode ierr;
657: DMStagDuplicateWithoutSetup(dm,PetscObjectComm((PetscObject)dm),newdm);
658: DMStagSetDOF(*newdm,dof0,dof1,dof2,dof3);
659: DMSetUp(*newdm);
660: return(0);
661: }
663: /*@C
664: DMStagGetLocationSlot - get index to use in accessing raw local arrays
666: Not Collective
668: Input Parameters:
669: + dm - the DMStag object
670: . loc - location relative to an element
671: - c - component
673: Output Parameter:
674: . slot - index to use
676: Notes:
677: Provides an appropriate index to use with DMStagVecGetArray() and friends.
678: This is required so that the user doesn't need to know about the ordering of
679: dof associated with each local element.
681: Level: beginner
683: .seealso: DMSTAG, DMStagVecGetArray(), DMStagVecGetArrayRead(), DMStagGetDOF(), DMStagGetEntriesPerElement()
684: @*/
685: PetscErrorCode DMStagGetLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt c,PetscInt *slot)
686: {
687: DM_Stag * const stag = (DM_Stag*)dm->data;
691: if (PetscDefined(USE_DEBUG)) {
693: PetscInt dof;
694: DMStagGetLocationDOF(dm,loc,&dof);
695: if (dof < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Location %s has no dof attached",DMStagStencilLocations[loc]);
696: 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);
697: }
698: *slot = stag->locationOffsets[loc] + c;
699: return(0);
700: }
702: /*@C
703: DMStagMigrateVec - transfer a vector associated with a DMStag to a vector associated with a compatible DMStag
705: Collective
707: Input Parameters:
708: + dm - the source DMStag object
709: . vec - the source vector, compatible with dm
710: . dmTo - the compatible destination DMStag object
711: - vecTo - the destination vector, compatible with dmTo
713: Notes:
714: Extra dof are ignored, and unfilled dof are zeroed.
715: Currently only implemented to migrate global vectors to global vectors.
717: Level: advanced
719: .seealso: DMSTAG, DMStagCreateCompatibleDMStag(), DMGetCompatibility(), DMStagVecSplitToDMDA()
720: @*/
721: PetscErrorCode DMStagMigrateVec(DM dm,Vec vec,DM dmTo,Vec vecTo)
722: {
723: PetscErrorCode ierr;
724: DM_Stag * const stag = (DM_Stag*)dm->data;
725: DM_Stag * const stagTo = (DM_Stag*)dmTo->data;
726: PetscInt nLocalTo,nLocal,dim,i,j,k;
727: PetscInt start[DMSTAG_MAX_DIM],startGhost[DMSTAG_MAX_DIM],n[DMSTAG_MAX_DIM],nExtra[DMSTAG_MAX_DIM],offset[DMSTAG_MAX_DIM];
728: Vec vecToLocal,vecLocal;
729: PetscBool compatible,compatibleSet;
730: const PetscScalar *arr;
731: PetscScalar *arrTo;
732: const PetscInt epe = stag->entriesPerElement;
733: const PetscInt epeTo = stagTo->entriesPerElement;
740: DMGetCompatibility(dm,dmTo,&compatible,&compatibleSet);
741: if (!compatibleSet || !compatible) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"DMStag objects must be shown to be compatible");
742: DMGetDimension(dm,&dim);
743: VecGetLocalSize(vec,&nLocal);
744: VecGetLocalSize(vecTo,&nLocalTo);
745: if (nLocal != stag->entries|| nLocalTo !=stagTo->entries) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Vector migration only implemented for global vector to global vector.");
746: DMStagGetCorners(dm,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
747: DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],NULL,NULL,NULL);
748: for (i=0; i<DMSTAG_MAX_DIM; ++i) offset[i] = start[i]-startGhost[i];
750: /* Proceed by transferring to a local vector, copying, and transferring back to a global vector */
751: DMGetLocalVector(dm,&vecLocal);
752: DMGetLocalVector(dmTo,&vecToLocal);
753: DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal);
754: DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal);
755: VecGetArrayRead(vecLocal,&arr);
756: VecGetArray(vecToLocal,&arrTo);
757: /* Note that some superfluous copying of entries on partial dummy elements is done */
758: if (dim == 1) {
759: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
760: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
761: PetscInt si;
762: for (si=0; si<2; ++si) {
763: b += stag->dof[si];
764: bTo += stagTo->dof[si];
765: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[i*epeTo + dTo] = arr[i*epe + d];
766: for (; dTo < bTo ; ++dTo) arrTo[i*epeTo + dTo] = 0.0;
767: d = b;
768: }
769: }
770: } else if (dim == 2) {
771: const PetscInt epr = stag->nGhost[0] * epe;
772: const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
773: for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
774: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
775: const PetscInt base = j*epr + i*epe;
776: const PetscInt baseTo = j*eprTo + i*epeTo;
777: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
778: const PetscInt s[4] = {0,1,1,2}; /* Dimensions of points, in order */
779: PetscInt si;
780: for (si=0; si<4; ++si) {
781: b += stag->dof[s[si]];
782: bTo += stagTo->dof[s[si]];
783: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
784: for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
785: d = b;
786: }
787: }
788: }
789: } else if (dim == 3) {
790: const PetscInt epr = stag->nGhost[0] * epe;
791: const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
792: const PetscInt epl = stag->nGhost[1] * epr;
793: const PetscInt eplTo = stagTo->nGhost[1] * eprTo;
794: for (k=offset[2]; k<offset[2] + n[2] + nExtra[2]; ++k) {
795: for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
796: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
797: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
798: const PetscInt base = k*epl + j*epr + i*epe;
799: const PetscInt baseTo = k*eplTo + j*eprTo + i*epeTo;
800: const PetscInt s[8] = {0,1,1,2,1,2,2,3}; /* dimensions of points, in order */
801: PetscInt is;
802: for (is=0; is<8; ++is) {
803: b += stag->dof[s[is]];
804: bTo += stagTo->dof[s[is]];
805: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
806: for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
807: d = b;
808: }
809: }
810: }
811: }
812: } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
813: VecRestoreArrayRead(vecLocal,&arr);
814: VecRestoreArray(vecToLocal,&arrTo);
815: DMRestoreLocalVector(dm,&vecLocal);
816: DMLocalToGlobalBegin(dmTo,vecToLocal,INSERT_VALUES,vecTo);
817: DMLocalToGlobalEnd(dmTo,vecToLocal,INSERT_VALUES,vecTo);
818: DMRestoreLocalVector(dmTo,&vecToLocal);
819: return(0);
820: }
822: /*@C
823: DMStagPopulateLocalToGlobalInjective - populate an internal 1-to-1 local-to-global map
825: Collective
827: Creates an internal object which explicitly maps a single local degree of
828: freedom to each global degree of freedom. This is used, if populated,
829: instead of SCATTER_REVERSE_LOCAL with the (1-to-many, in general)
830: global-to-local map, when DMLocalToGlobal() is called with INSERT_VALUES.
831: This allows usage, for example, even in the periodic, 1-rank case, where
832: the inverse of the global-to-local map, even when restricted to on-rank
833: communication, is non-injective. This is at the cost of storing an additional
834: VecScatter object inside each DMStag object.
836: Input Parameter:
837: . dm - the DMStag object
839: Notes:
840: In normal usage, library users shouldn't be concerned with this function,
841: as it is called during DMSetUp(), when required.
843: Returns immediately if the internal map is already populated.
845: Developer Notes:
846: This could, if desired, be moved up to a general DM routine. It would allow,
847: for example, DMDA to support DMLocalToGlobal() with INSERT_VALUES,
848: even in the single-rank periodic case.
850: Level: developer
852: .seealso: DMSTAG, DMLocalToGlobal(), VecScatter
853: @*/
854: PetscErrorCode DMStagPopulateLocalToGlobalInjective(DM dm)
855: {
856: PetscErrorCode ierr;
857: PetscInt dim;
858: DM_Stag * const stag = (DM_Stag*)dm->data;
862: if (stag->ltog_injective) return(0); /* Don't re-populate */
863: DMGetDimension(dm,&dim);
864: switch (dim) {
865: case 1: DMStagPopulateLocalToGlobalInjective_1d(dm); break;
866: case 2: DMStagPopulateLocalToGlobalInjective_2d(dm); break;
867: case 3: DMStagPopulateLocalToGlobalInjective_3d(dm); break;
868: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
869: }
870: return(0);
871: }
873: static PetscErrorCode DMStagRestoreProductCoordinateArrays_Private(DM dm,void *arrX,void *arrY,void *arrZ,PetscBool read)
874: {
875: PetscErrorCode ierr;
876: PetscInt dim,d;
877: void* arr[DMSTAG_MAX_DIM];
878: DM dmCoord;
882: DMGetDimension(dm,&dim);
883: if (dim > DMSTAG_MAX_DIM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for %D dimensions",dim);
884: arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
885: DMGetCoordinateDM(dm,&dmCoord);
886: for (d=0; d<dim; ++d) {
887: DM subDM;
888: Vec coord1d_local;
890: /* Ignore unrequested arrays */
891: if (!arr[d]) continue;
893: DMProductGetDM(dmCoord,d,&subDM);
894: DMGetCoordinatesLocal(subDM,&coord1d_local);
895: if (read) {
896: DMStagVecRestoreArrayRead(subDM,coord1d_local,arr[d]);
897: } else {
898: DMStagVecRestoreArray(subDM,coord1d_local,arr[d]);
899: }
900: }
901: return(0);
902: }
904: /*@C
905: DMStagRestoreProductCoordinateArrays - restore local array access
907: Logically Collective
909: Input Parameter:
910: . dm - the DMStag object
912: Output Parameters:
913: . arrX,arrY,arrZ - local 1D coordinate arrays
915: Level: intermediate
917: Notes:
918: 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:
920: $ DMGetCoordinateDM(dm,&cdm);
921: $ for (d=0; d<3; ++d) {
922: $ DM subdm;
923: $ Vec coor,coor_local;
925: $ DMProductGetDM(cdm,d,&subdm);
926: $ DMGetCoordinates(subdm,&coor);
927: $ DMGetCoordinatesLocal(subdm,&coor_local);
928: $ DMLocalToGlobal(subdm,coor_local,INSERT_VALUES,coor);
929: $ PetscPrintf(PETSC_COMM_WORLD,"Coordinates dim %D:\n",d);
930: $ VecView(coor,PETSC_VIEWER_STDOUT_WORLD);
931: $ }
933: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
934: @*/
935: PetscErrorCode DMStagRestoreProductCoordinateArrays(DM dm,void *arrX,void *arrY,void *arrZ)
936: {
937: PetscErrorCode ierr;
940: DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_FALSE);
941: return(0);
942: }
944: /*@C
945: DMStagRestoreProductCoordinateArraysRead - restore local product array access, read-only
947: Logically Collective
949: Input Parameter:
950: . dm - the DMStag object
952: Output Parameters:
953: . arrX,arrY,arrZ - local 1D coordinate arrays
955: Level: intermediate
957: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
958: @*/
959: PetscErrorCode DMStagRestoreProductCoordinateArraysRead(DM dm,void *arrX,void *arrY,void *arrZ)
960: {
961: PetscErrorCode ierr;
964: DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_TRUE);
965: return(0);
966: }
968: /*@C
969: DMStagSetBoundaryTypes - set DMStag boundary types
971: Logically Collective; boundaryType0, boundaryType1, and boundaryType2 must contain common values
973: Input Parameters:
974: + dm - the DMStag object
975: - boundaryType0,boundaryType1,boundaryType2 - boundary types in each direction
977: Notes:
978: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
980: Level: advanced
982: .seealso: DMSTAG, DMBoundaryType, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDASetBoundaryType()
983: @*/
984: PetscErrorCode DMStagSetBoundaryTypes(DM dm,DMBoundaryType boundaryType0,DMBoundaryType boundaryType1,DMBoundaryType boundaryType2)
985: {
986: PetscErrorCode ierr;
987: DM_Stag * const stag = (DM_Stag*)dm->data;
988: PetscInt dim;
991: DMGetDimension(dm,&dim);
996: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
997: stag->boundaryType[0] = boundaryType0;
998: if (dim > 1) stag->boundaryType[1] = boundaryType1;
999: if (dim > 2) stag->boundaryType[2] = boundaryType2;
1000: return(0);
1001: }
1003: /*@C
1004: DMStagSetCoordinateDMType - set DM type to store coordinates
1006: Logically Collective; dmtype must contain common value
1008: Input Parameters:
1009: + dm - the DMStag object
1010: - dmtype - DMtype for coordinates, either DMSTAG or DMPRODUCT
1012: Level: advanced
1014: .seealso: DMSTAG, DMPRODUCT, DMGetCoordinateDM(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMType
1015: @*/
1016: PetscErrorCode DMStagSetCoordinateDMType(DM dm,DMType dmtype)
1017: {
1018: PetscErrorCode ierr;
1019: DM_Stag * const stag = (DM_Stag*)dm->data;
1023: PetscFree(stag->coordinateDMType);
1024: PetscStrallocpy(dmtype,(char**)&stag->coordinateDMType);
1025: return(0);
1026: }
1028: /*@C
1029: DMStagSetDOF - set dof/stratum
1031: Logically Collective; dof0, dof1, dof2, and dof3 must contain common values
1033: Input Parameters:
1034: + dm - the DMStag object
1035: - dof0,dof1,dof2,dof3 - dof per stratum
1037: Notes:
1038: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1040: Level: advanced
1042: .seealso: DMSTAG, DMDASetDof()
1043: @*/
1044: PetscErrorCode DMStagSetDOF(DM dm,PetscInt dof0, PetscInt dof1,PetscInt dof2,PetscInt dof3)
1045: {
1046: PetscErrorCode ierr;
1047: DM_Stag * const stag = (DM_Stag*)dm->data;
1048: PetscInt dim;
1056: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1057: DMGetDimension(dm,&dim);
1058: if (dof0 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof0 cannot be negative");
1059: if (dof1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof1 cannot be negative");
1060: if (dim > 1 && dof2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof2 cannot be negative");
1061: if (dim > 2 && dof3 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof3 cannot be negative");
1062: stag->dof[0] = dof0;
1063: stag->dof[1] = dof1;
1064: if (dim > 1) stag->dof[2] = dof2;
1065: if (dim > 2) stag->dof[3] = dof3;
1066: return(0);
1067: }
1069: /*@C
1070: DMStagSetNumRanks - set ranks in each direction in the global rank grid
1072: Logically Collective; nRanks0, nRanks1, and nRanks2 must contain common values
1074: Input Parameters:
1075: + dm - the DMStag object
1076: - nRanks0,nRanks1,nRanks2 - number of ranks in each direction
1078: Notes:
1079: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1081: Level: developer
1083: .seealso: DMSTAG, DMDASetNumProcs()
1084: @*/
1085: PetscErrorCode DMStagSetNumRanks(DM dm,PetscInt nRanks0,PetscInt nRanks1,PetscInt nRanks2)
1086: {
1087: PetscErrorCode ierr;
1088: DM_Stag * const stag = (DM_Stag*)dm->data;
1089: PetscInt dim;
1096: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1097: DMGetDimension(dm,&dim);
1098: 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");
1099: 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");
1100: 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");
1101: if (nRanks0) stag->nRanks[0] = nRanks0;
1102: if (dim > 1 && nRanks1) stag->nRanks[1] = nRanks1;
1103: if (dim > 2 && nRanks2) stag->nRanks[2] = nRanks2;
1104: return(0);
1105: }
1107: /*@C
1108: DMStagSetStencilType - set elementwise ghost/halo stencil type
1110: Logically Collective; stencilType must contain common value
1112: Input Parameters:
1113: + dm - the DMStag object
1114: - stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE
1116: Level: beginner
1118: .seealso: DMSTAG, DMStagGetStencilType(), DMStagStencilType, DMDASetStencilType()
1119: @*/
1120: PetscErrorCode DMStagSetStencilType(DM dm,DMStagStencilType stencilType)
1121: {
1122: DM_Stag * const stag = (DM_Stag*)dm->data;
1127: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1128: stag->stencilType = stencilType;
1129: return(0);
1130: }
1132: /*@C
1133: DMStagSetStencilWidth - set elementwise stencil width
1135: Logically Collective; stencilWidth must contain common value
1137: Input Parameters:
1138: + dm - the DMStag object
1139: - stencilWidth - stencil/halo/ghost width in elements
1141: Level: beginner
1143: .seealso: DMSTAG, DMDASetStencilWidth()
1144: @*/
1145: PetscErrorCode DMStagSetStencilWidth(DM dm,PetscInt stencilWidth)
1146: {
1147: DM_Stag * const stag = (DM_Stag*)dm->data;
1152: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1153: if (stencilWidth < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Stencil width must be non-negative");
1154: stag->stencilWidth = stencilWidth;
1155: return(0);
1156: }
1158: /*@C
1159: DMStagSetGlobalSizes - set global element counts in each direction
1161: Logically Collective; N0, N1, and N2 must contain common values
1163: Input Parameters:
1164: + dm - the DMStag object
1165: - N0,N1,N2 - global elementwise sizes
1167: Notes:
1168: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1170: Level: advanced
1172: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMDASetSizes()
1173: @*/
1174: PetscErrorCode DMStagSetGlobalSizes(DM dm,PetscInt N0,PetscInt N1,PetscInt N2)
1175: {
1176: PetscErrorCode ierr;
1177: DM_Stag * const stag = (DM_Stag*)dm->data;
1178: PetscInt dim;
1182: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1183: DMGetDimension(dm,&dim);
1184: if (N0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in X direction must be positive");
1185: if (dim > 1 && N1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Y direction must be positive");
1186: if (dim > 2 && N2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Z direction must be positive");
1187: if (N0) stag->N[0] = N0;
1188: if (N1) stag->N[1] = N1;
1189: if (N2) stag->N[2] = N2;
1190: return(0);
1191: }
1193: /*@C
1194: DMStagSetOwnershipRanges - set elements per rank in each direction
1196: Logically Collective; lx, ly, and lz must contain common values
1198: Input Parameters:
1199: + dm - the DMStag object
1200: - lx,ly,lz - element counts for each rank in each direction
1202: Notes:
1203: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
1205: Level: developer
1207: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagGetOwnershipRanges(), DMDASetOwnershipRanges()
1208: @*/
1209: PetscErrorCode DMStagSetOwnershipRanges(DM dm,PetscInt const *lx,PetscInt const *ly,PetscInt const *lz)
1210: {
1211: PetscErrorCode ierr;
1212: DM_Stag * const stag = (DM_Stag*)dm->data;
1213: const PetscInt *lin[3];
1214: PetscInt d,dim;
1218: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1219: lin[0] = lx; lin[1] = ly; lin[2] = lz;
1220: DMGetDimension(dm,&dim);
1221: for (d=0; d<dim; ++d) {
1222: if (lin[d]) {
1223: if (stag->nRanks[d] < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of ranks");
1224: if (!stag->l[d]) {
1225: PetscMalloc1(stag->nRanks[d], &stag->l[d]);
1226: }
1227: PetscArraycpy(stag->l[d], lin[d], stag->nRanks[d]);
1228: }
1229: }
1230: return(0);
1231: }
1233: /*@C
1234: DMStagSetUniformCoordinates - set DMStag coordinates to be a uniform grid
1236: Collective
1238: Input Parameters:
1239: + dm - the DMStag object
1240: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1242: Notes:
1243: DMStag supports 2 different types of coordinate DM: DMSTAG and DMPRODUCT.
1244: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1246: Level: advanced
1248: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType(), DMGetCoordinateDM(), DMGetCoordinates(), DMDASetUniformCoordinates()
1249: @*/
1250: PetscErrorCode DMStagSetUniformCoordinates(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1251: {
1252: PetscErrorCode ierr;
1253: DM_Stag * const stag = (DM_Stag*)dm->data;
1254: PetscBool flg_stag,flg_product;
1258: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1259: if (!stag->coordinateDMType) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"You must first call DMStagSetCoordinateDMType()");
1260: PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg_stag);
1261: PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg_product);
1262: if (flg_stag) {
1263: DMStagSetUniformCoordinatesExplicit(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1264: } else if (flg_product) {
1265: DMStagSetUniformCoordinatesProduct(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1266: } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported DM Type %s",stag->coordinateDMType);
1267: return(0);
1268: }
1270: /*@C
1271: DMStagSetUniformCoordinatesExplicit - set DMStag coordinates to be a uniform grid, storing all values
1273: Collective
1275: Input Parameters:
1276: + dm - the DMStag object
1277: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1279: Notes:
1280: DMStag supports 2 different types of coordinate DM: either another DMStag, or a DMProduct.
1281: If the grid is orthogonal, using DMProduct should be more efficient.
1282: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1284: Level: beginner
1286: .seealso: DMSTAG, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType()
1287: @*/
1288: PetscErrorCode DMStagSetUniformCoordinatesExplicit(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1289: {
1290: PetscErrorCode ierr;
1291: DM_Stag * const stag = (DM_Stag*)dm->data;
1292: PetscInt dim;
1293: PetscBool flg;
1297: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1298: PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg);
1299: if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1300: DMStagSetCoordinateDMType(dm,DMSTAG);
1301: DMGetDimension(dm,&dim);
1302: switch (dim) {
1303: case 1: DMStagSetUniformCoordinatesExplicit_1d(dm,xmin,xmax); break;
1304: case 2: DMStagSetUniformCoordinatesExplicit_2d(dm,xmin,xmax,ymin,ymax); break;
1305: case 3: DMStagSetUniformCoordinatesExplicit_3d(dm,xmin,xmax,ymin,ymax,zmin,zmax); break;
1306: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1307: }
1308: return(0);
1309: }
1311: /*@C
1312: DMStagSetUniformCoordinatesProduct - create uniform coordinates, as a product of 1D arrays
1314: 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.
1316: Collective
1318: Input Parameters:
1319: + dm - the DMStag object
1320: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1322: Notes:
1323: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1325: The per-dimension 1-dimensional DMStag objects that comprise the product
1326: always have active 0-cells (vertices, element boundaries) and 1-cells
1327: (element centers).
1329: Level: intermediate
1331: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetCoordinateDMType()
1332: @*/
1333: PetscErrorCode DMStagSetUniformCoordinatesProduct(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1334: {
1335: PetscErrorCode ierr;
1336: DM_Stag * const stag = (DM_Stag*)dm->data;
1337: DM dmc;
1338: PetscInt dim,d,dof0,dof1;
1339: PetscBool flg;
1343: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1344: PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg);
1345: if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1346: DMStagSetCoordinateDMType(dm,DMPRODUCT);
1347: DMGetCoordinateDM(dm,&dmc);
1348: DMGetDimension(dm,&dim);
1350: /* Create 1D sub-DMs, living on subcommunicators.
1351: Always include both vertex and element dof, regardless of the active strata of the DMStag */
1352: dof0 = 1;
1353: dof1 = 1;
1355: for (d=0; d<dim; ++d) {
1356: DM subdm;
1357: MPI_Comm subcomm;
1358: PetscMPIInt color;
1359: const PetscMPIInt key = 0; /* let existing rank break ties */
1361: /* Choose colors based on position in the plane orthogonal to this dim, and split */
1362: switch (d) {
1363: case 0: color = (dim > 1 ? stag->rank[1] : 0) + (dim > 2 ? stag->nRanks[1]*stag->rank[2] : 0); break;
1364: case 1: color = stag->rank[0] + (dim > 2 ? stag->nRanks[0]*stag->rank[2] : 0); break;
1365: case 2: color = stag->rank[0] + stag->nRanks[0]*stag->rank[1] ; break;
1366: default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1367: }
1368: MPI_Comm_split(PetscObjectComm((PetscObject)dm),color,key,&subcomm);
1370: /* Create sub-DMs living on these new communicators (which are destroyed by DMProduct) */
1371: DMStagCreate1d(subcomm,stag->boundaryType[d],stag->N[d],dof0,dof1,stag->stencilType,stag->stencilWidth,stag->l[d],&subdm);
1372: DMSetUp(subdm);
1373: switch (d) {
1374: case 0:
1375: DMStagSetUniformCoordinatesExplicit(subdm,xmin,xmax,0,0,0,0);
1376: break;
1377: case 1:
1378: DMStagSetUniformCoordinatesExplicit(subdm,ymin,ymax,0,0,0,0);
1379: break;
1380: case 2:
1381: DMStagSetUniformCoordinatesExplicit(subdm,zmin,zmax,0,0,0,0);
1382: break;
1383: default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1384: }
1385: DMProductSetDM(dmc,d,subdm);
1386: DMProductSetDimensionIndex(dmc,d,0);
1387: DMDestroy(&subdm);
1388: MPI_Comm_free(&subcomm);
1389: }
1390: return(0);
1391: }
1393: /*@C
1394: DMStagVecGetArray - get access to local array
1396: Logically Collective
1398: This function returns a (dim+1)-dimensional array for a dim-dimensional
1399: DMStag.
1401: The first 1-3 dimensions indicate an element in the global
1402: numbering, using the standard C ordering.
1404: The final dimension in this array corresponds to a degree
1405: of freedom with respect to this element, for example corresponding to
1406: the element or one of its neighboring faces, edges, or vertices.
1408: For example, for a 3D DMStag, indexing is array[k][j][i][idx], where k is the
1409: index in the z-direction, j is the index in the y-direction, and i is the
1410: index in the x-direction.
1412: "idx" is obtained with DMStagGetLocationSlot(), since the correct offset
1413: into the (dim+1)-dimensional C array depends on the grid size and the number
1414: of dof stored at each location.
1416: Input Parameters:
1417: + dm - the DMStag object
1418: - vec - the Vec object
1420: Output Parameters:
1421: . array - the array
1423: Notes:
1424: DMStagVecRestoreArray() must be called, once finished with the array
1426: Level: beginner
1428: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArray(), DMDAVecGetArrayDOF()
1429: @*/
1430: PetscErrorCode DMStagVecGetArray(DM dm,Vec vec,void *array)
1431: {
1432: PetscErrorCode ierr;
1433: DM_Stag * const stag = (DM_Stag*)dm->data;
1434: PetscInt dim;
1435: PetscInt nLocal;
1440: DMGetDimension(dm,&dim);
1441: VecGetLocalSize(vec,&nLocal);
1442: 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);
1443: switch (dim) {
1444: case 1:
1445: VecGetArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1446: break;
1447: case 2:
1448: VecGetArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1449: break;
1450: case 3:
1451: 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);
1452: break;
1453: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1454: }
1455: return(0);
1456: }
1458: /*@C
1459: DMStagVecGetArrayRead - get read-only access to a local array
1461: Logically Collective
1463: See the man page for DMStagVecGetArray() for more information.
1465: Input Parameters:
1466: + dm - the DMStag object
1467: - vec - the Vec object
1469: Output Parameters:
1470: . array - the read-only array
1472: Notes:
1473: DMStagVecRestoreArrayRead() must be called, once finished with the array
1475: Level: beginner
1477: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArrayRead(), DMDAVecGetArrayDOFRead()
1478: @*/
1479: PetscErrorCode DMStagVecGetArrayRead(DM dm,Vec vec,void *array)
1480: {
1481: PetscErrorCode ierr;
1482: DM_Stag * const stag = (DM_Stag*)dm->data;
1483: PetscInt dim;
1484: PetscInt nLocal;
1489: DMGetDimension(dm,&dim);
1490: VecGetLocalSize(vec,&nLocal);
1491: 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);
1492: switch (dim) {
1493: case 1:
1494: VecGetArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1495: break;
1496: case 2:
1497: VecGetArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1498: break;
1499: case 3:
1500: 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);
1501: break;
1502: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1503: }
1504: return(0);
1505: }
1507: /*@C
1508: DMStagVecRestoreArray - restore access to a raw array
1510: Logically Collective
1512: Input Parameters:
1513: + dm - the DMStag object
1514: - vec - the Vec object
1516: Output Parameters:
1517: . array - the array
1519: Level: beginner
1521: .seealso: DMSTAG, DMStagVecGetArray(), DMDAVecRestoreArray(), DMDAVecRestoreArrayDOF()
1522: @*/
1523: PetscErrorCode DMStagVecRestoreArray(DM dm,Vec vec,void *array)
1524: {
1525: PetscErrorCode ierr;
1526: DM_Stag * const stag = (DM_Stag*)dm->data;
1527: PetscInt dim;
1528: PetscInt nLocal;
1533: DMGetDimension(dm,&dim);
1534: VecGetLocalSize(vec,&nLocal);
1535: 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);
1536: switch (dim) {
1537: case 1:
1538: VecRestoreArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1539: break;
1540: case 2:
1541: VecRestoreArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1542: break;
1543: case 3:
1544: 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);
1545: break;
1546: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1547: }
1548: return(0);
1549: }
1551: /*@C
1552: DMStagVecRestoreArrayRead - restore read-only access to a raw array
1554: Logically Collective
1556: Input Parameters:
1557: + dm - the DMStag object
1558: - vec - the Vec object
1560: Output Parameters:
1561: . array - the read-only array
1563: Level: beginner
1565: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMDAVecRestoreArrayRead(), DMDAVecRestoreArrayDOFRead()
1566: @*/
1567: PetscErrorCode DMStagVecRestoreArrayRead(DM dm,Vec vec,void *array)
1568: {
1569: PetscErrorCode ierr;
1570: DM_Stag * const stag = (DM_Stag*)dm->data;
1571: PetscInt dim;
1572: PetscInt nLocal;
1577: DMGetDimension(dm,&dim);
1578: VecGetLocalSize(vec,&nLocal);
1579: 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);
1580: switch (dim) {
1581: case 1:
1582: VecRestoreArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1583: break;
1584: case 2:
1585: VecRestoreArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1586: break;
1587: case 3:
1588: 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);
1589: break;
1590: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1591: }
1592: return(0);
1593: }