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 - starting element index in first direction
234: . y - starting element index in second direction
235: . z - starting element index in third direction
236: . m - element width in first direction
237: . n - element width in second direction
238: . p - element width in third direction
239: . nExtrax - number of extra partial elements in first direction
240: . nExtray - number of extra partial elements in second direction
241: - nExtraz - number of extra partial elements in third direction
243: Notes:
244: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
246: The number of extra partial elements is either 1 or 0.
247: The value is 1 on right, top, and front non-periodic domain ("physical") boundaries,
248: in the x, y, and z directions respectively, and otherwise 0.
250: Level: beginner
252: .seealso: DMSTAG, DMStagGetGhostCorners(), DMDAGetCorners()
253: @*/
254: PetscErrorCode DMStagGetCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *nExtrax,PetscInt *nExtray,PetscInt *nExtraz)
255: {
256: const DM_Stag * const stag = (DM_Stag*)dm->data;
260: if (x) *x = stag->start[0];
261: if (y) *y = stag->start[1];
262: if (z) *z = stag->start[2];
263: if (m) *m = stag->n[0];
264: if (n) *n = stag->n[1];
265: if (p) *p = stag->n[2];
266: if (nExtrax) *nExtrax = stag->boundaryType[0] != DM_BOUNDARY_PERIODIC && stag->lastRank[0] ? 1 : 0;
267: if (nExtray) *nExtray = stag->boundaryType[1] != DM_BOUNDARY_PERIODIC && stag->lastRank[1] ? 1 : 0;
268: if (nExtraz) *nExtraz = stag->boundaryType[2] != DM_BOUNDARY_PERIODIC && stag->lastRank[2] ? 1 : 0;
269: return(0);
270: }
272: /*@C
273: DMStagGetDOF - get number of DOF associated with each stratum of the grid
275: Not Collective
277: Input Parameter:
278: . dm - the DMStag object
280: Output Parameters:
281: + dof0 - the number of points per 0-cell (vertex/node)
282: . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
283: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
284: - dof3 - the number of points per 3-cell (element in 3D)
286: Level: beginner
288: .seealso: DMSTAG, DMStagGetCorners(), DMStagGetGhostCorners(), DMStagGetGlobalSizes(), DMStagGetStencilWidth(), DMStagGetBoundaryTypes(), DMStagGetLocationDof(), DMDAGetDof()
289: @*/
290: PetscErrorCode DMStagGetDOF(DM dm,PetscInt *dof0,PetscInt *dof1,PetscInt *dof2,PetscInt *dof3)
291: {
292: const DM_Stag * const stag = (DM_Stag*)dm->data;
296: if (dof0) *dof0 = stag->dof[0];
297: if (dof1) *dof1 = stag->dof[1];
298: if (dof2) *dof2 = stag->dof[2];
299: if (dof3) *dof3 = stag->dof[3];
300: return(0);
301: }
303: /*@C
304: DMStagGetGhostCorners - return global element indices of the local region, including ghost points
306: Not Collective
308: Input Parameter:
309: . dm - the DMStag object
311: Output Parameters:
312: + x - the starting element index in the first direction
313: . y - the starting element index in the second direction
314: . z - the starting element index in the third direction
315: . m - the element width in the first direction
316: . n - the element width in the second direction
317: - p - the element width in the third direction
319: Notes:
320: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
322: Level: beginner
324: .seealso: DMSTAG, DMStagGetCorners(), DMDAGetGhostCorners()
325: @*/
326: PetscErrorCode DMStagGetGhostCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
327: {
328: const DM_Stag * const stag = (DM_Stag*)dm->data;
332: if (x) *x = stag->startGhost[0];
333: if (y) *y = stag->startGhost[1];
334: if (z) *z = stag->startGhost[2];
335: if (m) *m = stag->nGhost[0];
336: if (n) *n = stag->nGhost[1];
337: if (p) *p = stag->nGhost[2];
338: return(0);
339: }
341: /*@C
342: DMStagGetGlobalSizes - get global element counts
344: Not Collective
346: Input Parameter:
347: . dm - the DMStag object
349: Output Parameters:
350: . M,N,P - global element counts in each direction
352: Notes:
353: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
355: Level: beginner
357: .seealso: DMSTAG, DMStagGetLocalSizes(), DMDAGetInfo()
358: @*/
359: PetscErrorCode DMStagGetGlobalSizes(DM dm,PetscInt* M,PetscInt* N,PetscInt* P)
360: {
361: const DM_Stag * const stag = (DM_Stag*)dm->data;
365: if (M) *M = stag->N[0];
366: if (N) *N = stag->N[1];
367: if (P) *P = stag->N[2];
368: return(0);
369: }
371: /*@C
372: DMStagGetIsFirstRank - get boolean value for whether this rank is first in each direction in the rank grid
374: Not Collective
376: Input Parameter:
377: . dm - the DMStag object
379: Output Parameters:
380: . isFirstRank0,isFirstRank1,isFirstRank2 - whether this rank is first in each direction
382: Level: intermediate
384: Notes:
385: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
387: .seealso: DMSTAG, DMStagGetIsLastRank()
388: @*/
389: PetscErrorCode DMStagGetIsFirstRank(DM dm,PetscBool *isFirstRank0,PetscBool *isFirstRank1,PetscBool *isFirstRank2)
390: {
391: const DM_Stag * const stag = (DM_Stag*)dm->data;
395: if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
396: if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
397: if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
398: return(0);
399: }
401: /*@C
402: DMStagGetIsLastRank - get boolean value for whether this rank is last in each direction in the rank grid
404: Not Collective
406: Input Parameter:
407: . dm - the DMStag object
409: Output Parameters:
410: . isLastRank0,isLastRank1,isLastRank2 - whether this rank is last in each direction
412: Level: intermediate
414: Notes:
415: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
416: Level: intermediate
418: .seealso: DMSTAG, DMStagGetIsFirstRank()
419: @*/
420: PetscErrorCode DMStagGetIsLastRank(DM dm,PetscBool *isLastRank0,PetscBool *isLastRank1,PetscBool *isLastRank2)
421: {
422: const DM_Stag * const stag = (DM_Stag*)dm->data;
426: if (isLastRank0) *isLastRank0 = stag->lastRank[0];
427: if (isLastRank1) *isLastRank1 = stag->lastRank[1];
428: if (isLastRank2) *isLastRank2 = stag->lastRank[2];
429: return(0);
430: }
432: /*@C
433: DMStagGetLocalSizes - get local elementwise sizes
435: Not Collective
437: Input Parameter:
438: . dm - the DMStag object
440: Output Parameters:
441: . m,n,p - local element counts (excluding ghosts) in each direction
443: Notes:
444: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
445: Level: intermediate
447: Level: beginner
449: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetDOF(), DMStagGetNumRanks(), DMDAGetLocalInfo()
450: @*/
451: PetscErrorCode DMStagGetLocalSizes(DM dm,PetscInt* m,PetscInt* n,PetscInt* p)
452: {
453: const DM_Stag * const stag = (DM_Stag*)dm->data;
457: if (m) *m = stag->n[0];
458: if (n) *n = stag->n[1];
459: if (p) *p = stag->n[2];
460: return(0);
461: }
463: /*@C
464: DMStagGetNumRanks - get number of ranks in each direction in the global grid decomposition
466: Not Collective
468: Input Parameter:
469: . dm - the DMStag object
471: Output Parameters:
472: . nRanks0,nRanks1,nRanks2 - number of ranks in each direction in the grid decomposition
474: Notes:
475: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
476: Level: intermediate
478: Level: beginner
480: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetLocalSize(), DMStagSetNumRanks(), DMDAGetInfo()
481: @*/
482: PetscErrorCode DMStagGetNumRanks(DM dm,PetscInt *nRanks0,PetscInt *nRanks1,PetscInt *nRanks2)
483: {
484: const DM_Stag * const stag = (DM_Stag*)dm->data;
488: if (nRanks0) *nRanks0 = stag->nRanks[0];
489: if (nRanks1) *nRanks1 = stag->nRanks[1];
490: if (nRanks2) *nRanks2 = stag->nRanks[2];
491: return(0);
492: }
494: /*@C
495: DMStagGetEntries - get number of native entries in the global representation
497: Not Collective
499: Input Parameter:
500: . dm - the DMStag object
502: Output Parameters:
503: . entries - number of rank-native entries in the global representation
505: Note:
506: This is the number of entries on this rank for a global vector associated with dm.
508: Level: developer
510: .seealso: DMSTAG, DMStagGetDOF(), DMStagGetEntriesPerElement(), DMCreateLocalVector()
511: @*/
512: PetscErrorCode DMStagGetEntries(DM dm,PetscInt *entries)
513: {
514: const DM_Stag * const stag = (DM_Stag*)dm->data;
518: if (entries) *entries = stag->entries;
519: return(0);
520: }
522: /*@C
523: DMStagGetEntriesPerElement - get number of entries per element in the local representation
525: Not Collective
527: Input Parameter:
528: . dm - the DMStag object
530: Output Parameters:
531: . entriesPerElement - number of entries associated with each element in the local representation
533: Notes:
534: This is the natural block size for most local operations. In 1D it is equal to dof0 + dof1,
535: in 2D it is equal to dof0 + 2*dof1 + dof2, and in 3D it is equal to dof0 + 3*dof1 + 3*dof2 + dof3
537: Level: developer
539: .seealso: DMSTAG, DMStagGetDOF()
540: @*/
541: PetscErrorCode DMStagGetEntriesPerElement(DM dm,PetscInt *entriesPerElement)
542: {
543: const DM_Stag * const stag = (DM_Stag*)dm->data;
547: if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
548: return(0);
549: }
551: /*@C
552: DMStagGetStencilType - get elementwise ghost/halo stencil type
554: Not Collective
556: Input Parameter:
557: . dm - the DMStag object
559: Output Parameter:
560: . stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE
562: Level: beginner
564: .seealso: DMSTAG, DMStagSetStencilType(), DMStagStencilType, DMDAGetInfo()
565: @*/
566: PetscErrorCode DMStagGetStencilType(DM dm,DMStagStencilType *stencilType)
567: {
568: DM_Stag * const stag = (DM_Stag*)dm->data;
572: *stencilType = stag->stencilType;
573: return(0);
574: }
576: /*@C
577: DMStagGetStencilWidth - get elementwise stencil width
579: Not Collective
581: Input Parameter:
582: . dm - the DMStag object
584: Output Parameters:
585: . stencilWidth - stencil/halo/ghost width in elements
587: Level: beginner
589: .seealso: DMSTAG, DMDAGetStencilWidth(), DMDAGetInfo()
590: @*/
591: PetscErrorCode DMStagGetStencilWidth(DM dm,PetscInt *stencilWidth)
592: {
593: const DM_Stag * const stag = (DM_Stag*)dm->data;
597: if (stencilWidth) *stencilWidth = stag->stencilWidth;
598: return(0);
599: }
601: /*@C
602: DMStagGetOwnershipRanges - get elements per rank in each direction
604: Not Collective
606: Input Parameter:
607: . dm - the DMStag object
609: Output Parameters:
610: + lx - ownership along x direction (optional)
611: . ly - ownership along y direction (optional)
612: - lz - ownership along z direction (optional)
614: Notes:
615: These correspond to the optional final arguments passed to DMStagCreate1d(), DMStagCreate2d(), and DMStagCreate3d().
617: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
619: In C you should not free these arrays, nor change the values in them.
620: They will only have valid values while the DMStag they came from still exists (has not been destroyed).
622: Level: intermediate
624: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagSetOwnershipRanges(), DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDAGetOwnershipRanges()
625: @*/
626: PetscErrorCode DMStagGetOwnershipRanges(DM dm,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
627: {
628: const DM_Stag * const stag = (DM_Stag*)dm->data;
632: if (lx) *lx = stag->l[0];
633: if (ly) *ly = stag->l[1];
634: if (lz) *lz = stag->l[2];
635: return(0);
636: }
638: /*@C
639: DMStagCreateCompatibleDMStag - create a compatible DMStag with different dof/stratum
641: Collective
643: Input Parameters:
644: + dm - the DMStag object
645: - dof0,dof1,dof2,dof3 - number of dof on each stratum in the new DMStag
647: Output Parameters:
648: . newdm - the new, compatible DMStag
650: Notes:
651: Dof supplied for strata too big for the dimension are ignored; these may be set to 0.
652: For example, for a 2-dimensional DMStag, dof2 sets the number of dof per element,
653: and dof3 is unused. For a 3-dimensional DMStag, dof3 sets the number of dof per element.
655: In contrast to DMDACreateCompatibleDMDA(), coordinates are not reused.
657: Level: intermediate
659: .seealso: DMSTAG, DMDACreateCompatibleDMDA(), DMGetCompatibility(), DMStagMigrateVec()
660: @*/
661: PetscErrorCode DMStagCreateCompatibleDMStag(DM dm,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DM *newdm)
662: {
663: PetscErrorCode ierr;
667: DMStagDuplicateWithoutSetup(dm,PetscObjectComm((PetscObject)dm),newdm);
668: DMStagSetDOF(*newdm,dof0,dof1,dof2,dof3);
669: DMSetUp(*newdm);
670: return(0);
671: }
673: /*@C
674: DMStagGetLocationSlot - get index to use in accessing raw local arrays
676: Not Collective
678: Input Parameters:
679: + dm - the DMStag object
680: . loc - location relative to an element
681: - c - component
683: Output Parameter:
684: . slot - index to use
686: Notes:
687: Provides an appropriate index to use with DMStagVecGetArray() and friends.
688: This is required so that the user doesn't need to know about the ordering of
689: dof associated with each local element.
691: Level: beginner
693: .seealso: DMSTAG, DMStagVecGetArray(), DMStagVecGetArrayRead(), DMStagGetDOF(), DMStagGetEntriesPerElement()
694: @*/
695: PetscErrorCode DMStagGetLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt c,PetscInt *slot)
696: {
697: DM_Stag * const stag = (DM_Stag*)dm->data;
701: if (PetscDefined(USE_DEBUG)) {
703: PetscInt dof;
704: DMStagGetLocationDOF(dm,loc,&dof);
705: if (dof < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Location %s has no dof attached",DMStagStencilLocations[loc]);
706: 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);
707: }
708: *slot = stag->locationOffsets[loc] + c;
709: return(0);
710: }
712: /*@C
713: DMStagMigrateVec - transfer a vector associated with a DMStag to a vector associated with a compatible DMStag
715: Collective
717: Input Parameters:
718: + dm - the source DMStag object
719: . vec - the source vector, compatible with dm
720: . dmTo - the compatible destination DMStag object
721: - vecTo - the destination vector, compatible with dmTo
723: Notes:
724: Extra dof are ignored, and unfilled dof are zeroed.
725: Currently only implemented to migrate global vectors to global vectors.
727: Level: advanced
729: .seealso: DMSTAG, DMStagCreateCompatibleDMStag(), DMGetCompatibility(), DMStagVecSplitToDMDA()
730: @*/
731: PetscErrorCode DMStagMigrateVec(DM dm,Vec vec,DM dmTo,Vec vecTo)
732: {
733: PetscErrorCode ierr;
734: DM_Stag * const stag = (DM_Stag*)dm->data;
735: DM_Stag * const stagTo = (DM_Stag*)dmTo->data;
736: PetscInt nLocalTo,nLocal,dim,i,j,k;
737: PetscInt start[DMSTAG_MAX_DIM],startGhost[DMSTAG_MAX_DIM],n[DMSTAG_MAX_DIM],nExtra[DMSTAG_MAX_DIM],offset[DMSTAG_MAX_DIM];
738: Vec vecToLocal,vecLocal;
739: PetscBool compatible,compatibleSet;
740: const PetscScalar *arr;
741: PetscScalar *arrTo;
742: const PetscInt epe = stag->entriesPerElement;
743: const PetscInt epeTo = stagTo->entriesPerElement;
750: DMGetCompatibility(dm,dmTo,&compatible,&compatibleSet);
751: if (!compatibleSet || !compatible) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"DMStag objects must be shown to be compatible");
752: DMGetDimension(dm,&dim);
753: VecGetLocalSize(vec,&nLocal);
754: VecGetLocalSize(vecTo,&nLocalTo);
755: if (nLocal != stag->entries|| nLocalTo !=stagTo->entries) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Vector migration only implemented for global vector to global vector.");
756: DMStagGetCorners(dm,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
757: DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],NULL,NULL,NULL);
758: for (i=0; i<DMSTAG_MAX_DIM; ++i) offset[i] = start[i]-startGhost[i];
760: /* Proceed by transferring to a local vector, copying, and transferring back to a global vector */
761: DMGetLocalVector(dm,&vecLocal);
762: DMGetLocalVector(dmTo,&vecToLocal);
763: DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal);
764: DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal);
765: VecGetArrayRead(vecLocal,&arr);
766: VecGetArray(vecToLocal,&arrTo);
767: /* Note that some superfluous copying of entries on partial dummy elements is done */
768: if (dim == 1) {
769: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
770: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
771: PetscInt si;
772: for (si=0; si<2; ++si) {
773: b += stag->dof[si];
774: bTo += stagTo->dof[si];
775: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[i*epeTo + dTo] = arr[i*epe + d];
776: for (; dTo < bTo ; ++dTo) arrTo[i*epeTo + dTo] = 0.0;
777: d = b;
778: }
779: }
780: } else if (dim == 2) {
781: const PetscInt epr = stag->nGhost[0] * epe;
782: const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
783: for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
784: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
785: const PetscInt base = j*epr + i*epe;
786: const PetscInt baseTo = j*eprTo + i*epeTo;
787: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
788: const PetscInt s[4] = {0,1,1,2}; /* Dimensions of points, in order */
789: PetscInt si;
790: for (si=0; si<4; ++si) {
791: b += stag->dof[s[si]];
792: bTo += stagTo->dof[s[si]];
793: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
794: for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
795: d = b;
796: }
797: }
798: }
799: } else if (dim == 3) {
800: const PetscInt epr = stag->nGhost[0] * epe;
801: const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
802: const PetscInt epl = stag->nGhost[1] * epr;
803: const PetscInt eplTo = stagTo->nGhost[1] * eprTo;
804: for (k=offset[2]; k<offset[2] + n[2] + nExtra[2]; ++k) {
805: for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
806: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
807: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
808: const PetscInt base = k*epl + j*epr + i*epe;
809: const PetscInt baseTo = k*eplTo + j*eprTo + i*epeTo;
810: const PetscInt s[8] = {0,1,1,2,1,2,2,3}; /* dimensions of points, in order */
811: PetscInt is;
812: for (is=0; is<8; ++is) {
813: b += stag->dof[s[is]];
814: bTo += stagTo->dof[s[is]];
815: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
816: for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
817: d = b;
818: }
819: }
820: }
821: }
822: } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
823: VecRestoreArrayRead(vecLocal,&arr);
824: VecRestoreArray(vecToLocal,&arrTo);
825: DMRestoreLocalVector(dm,&vecLocal);
826: DMLocalToGlobalBegin(dmTo,vecToLocal,INSERT_VALUES,vecTo);
827: DMLocalToGlobalEnd(dmTo,vecToLocal,INSERT_VALUES,vecTo);
828: DMRestoreLocalVector(dmTo,&vecToLocal);
829: return(0);
830: }
832: /*@C
833: DMStagPopulateLocalToGlobalInjective - populate an internal 1-to-1 local-to-global map
835: Collective
837: Creates an internal object which explicitly maps a single local degree of
838: freedom to each global degree of freedom. This is used, if populated,
839: instead of SCATTER_REVERSE_LOCAL with the (1-to-many, in general)
840: global-to-local map, when DMLocalToGlobal() is called with INSERT_VALUES.
841: This allows usage, for example, even in the periodic, 1-rank case, where
842: the inverse of the global-to-local map, even when restricted to on-rank
843: communication, is non-injective. This is at the cost of storing an additional
844: VecScatter object inside each DMStag object.
846: Input Parameter:
847: . dm - the DMStag object
849: Notes:
850: In normal usage, library users shouldn't be concerned with this function,
851: as it is called during DMSetUp(), when required.
853: Returns immediately if the internal map is already populated.
855: Developer Notes:
856: This could, if desired, be moved up to a general DM routine. It would allow,
857: for example, DMDA to support DMLocalToGlobal() with INSERT_VALUES,
858: even in the single-rank periodic case.
860: Level: developer
862: .seealso: DMSTAG, DMLocalToGlobal(), VecScatter
863: @*/
864: PetscErrorCode DMStagPopulateLocalToGlobalInjective(DM dm)
865: {
866: PetscErrorCode ierr;
867: PetscInt dim;
868: DM_Stag * const stag = (DM_Stag*)dm->data;
872: if (stag->ltog_injective) return(0); /* Don't re-populate */
873: DMGetDimension(dm,&dim);
874: switch (dim) {
875: case 1: DMStagPopulateLocalToGlobalInjective_1d(dm); break;
876: case 2: DMStagPopulateLocalToGlobalInjective_2d(dm); break;
877: case 3: DMStagPopulateLocalToGlobalInjective_3d(dm); break;
878: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
879: }
880: return(0);
881: }
883: static PetscErrorCode DMStagRestoreProductCoordinateArrays_Private(DM dm,void *arrX,void *arrY,void *arrZ,PetscBool read)
884: {
885: PetscErrorCode ierr;
886: PetscInt dim,d;
887: void* arr[DMSTAG_MAX_DIM];
888: DM dmCoord;
892: DMGetDimension(dm,&dim);
893: if (dim > DMSTAG_MAX_DIM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for %D dimensions",dim);
894: arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
895: DMGetCoordinateDM(dm,&dmCoord);
896: for (d=0; d<dim; ++d) {
897: DM subDM;
898: Vec coord1d_local;
900: /* Ignore unrequested arrays */
901: if (!arr[d]) continue;
903: DMProductGetDM(dmCoord,d,&subDM);
904: DMGetCoordinatesLocal(subDM,&coord1d_local);
905: if (read) {
906: DMStagVecRestoreArrayRead(subDM,coord1d_local,arr[d]);
907: } else {
908: DMStagVecRestoreArray(subDM,coord1d_local,arr[d]);
909: }
910: }
911: return(0);
912: }
914: /*@C
915: DMStagRestoreProductCoordinateArrays - restore local array access
917: Logically Collective
919: Input Parameter:
920: . dm - the DMStag object
922: Output Parameters:
923: . arrX,arrY,arrZ - local 1D coordinate arrays
925: Level: intermediate
927: Notes:
928: 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:
930: $ DMGetCoordinateDM(dm,&cdm);
931: $ for (d=0; d<3; ++d) {
932: $ DM subdm;
933: $ Vec coor,coor_local;
935: $ DMProductGetDM(cdm,d,&subdm);
936: $ DMGetCoordinates(subdm,&coor);
937: $ DMGetCoordinatesLocal(subdm,&coor_local);
938: $ DMLocalToGlobal(subdm,coor_local,INSERT_VALUES,coor);
939: $ PetscPrintf(PETSC_COMM_WORLD,"Coordinates dim %D:\n",d);
940: $ VecView(coor,PETSC_VIEWER_STDOUT_WORLD);
941: $ }
943: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
944: @*/
945: PetscErrorCode DMStagRestoreProductCoordinateArrays(DM dm,void *arrX,void *arrY,void *arrZ)
946: {
947: PetscErrorCode ierr;
950: DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_FALSE);
951: return(0);
952: }
954: /*@C
955: DMStagRestoreProductCoordinateArraysRead - restore local product array access, read-only
957: Logically Collective
959: Input Parameter:
960: . dm - the DMStag object
962: Output Parameters:
963: . arrX,arrY,arrZ - local 1D coordinate arrays
965: Level: intermediate
967: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
968: @*/
969: PetscErrorCode DMStagRestoreProductCoordinateArraysRead(DM dm,void *arrX,void *arrY,void *arrZ)
970: {
971: PetscErrorCode ierr;
974: DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_TRUE);
975: return(0);
976: }
978: /*@C
979: DMStagSetBoundaryTypes - set DMStag boundary types
981: Logically Collective; boundaryType0, boundaryType1, and boundaryType2 must contain common values
983: Input Parameters:
984: + dm - the DMStag object
985: - boundaryType0,boundaryType1,boundaryType2 - boundary types in each direction
987: Notes:
988: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
990: Level: advanced
992: .seealso: DMSTAG, DMBoundaryType, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDASetBoundaryType()
993: @*/
994: PetscErrorCode DMStagSetBoundaryTypes(DM dm,DMBoundaryType boundaryType0,DMBoundaryType boundaryType1,DMBoundaryType boundaryType2)
995: {
996: PetscErrorCode ierr;
997: DM_Stag * const stag = (DM_Stag*)dm->data;
998: PetscInt dim;
1001: DMGetDimension(dm,&dim);
1006: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1007: stag->boundaryType[0] = boundaryType0;
1008: if (dim > 1) stag->boundaryType[1] = boundaryType1;
1009: if (dim > 2) stag->boundaryType[2] = boundaryType2;
1010: return(0);
1011: }
1013: /*@C
1014: DMStagSetCoordinateDMType - set DM type to store coordinates
1016: Logically Collective; dmtype must contain common value
1018: Input Parameters:
1019: + dm - the DMStag object
1020: - dmtype - DMtype for coordinates, either DMSTAG or DMPRODUCT
1022: Level: advanced
1024: .seealso: DMSTAG, DMPRODUCT, DMGetCoordinateDM(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMType
1025: @*/
1026: PetscErrorCode DMStagSetCoordinateDMType(DM dm,DMType dmtype)
1027: {
1028: PetscErrorCode ierr;
1029: DM_Stag * const stag = (DM_Stag*)dm->data;
1033: PetscFree(stag->coordinateDMType);
1034: PetscStrallocpy(dmtype,(char**)&stag->coordinateDMType);
1035: return(0);
1036: }
1038: /*@C
1039: DMStagSetDOF - set dof/stratum
1041: Logically Collective; dof0, dof1, dof2, and dof3 must contain common values
1043: Input Parameters:
1044: + dm - the DMStag object
1045: - dof0,dof1,dof2,dof3 - dof per stratum
1047: Notes:
1048: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1050: Level: advanced
1052: .seealso: DMSTAG, DMDASetDof()
1053: @*/
1054: PetscErrorCode DMStagSetDOF(DM dm,PetscInt dof0, PetscInt dof1,PetscInt dof2,PetscInt dof3)
1055: {
1056: PetscErrorCode ierr;
1057: DM_Stag * const stag = (DM_Stag*)dm->data;
1058: PetscInt dim;
1066: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1067: DMGetDimension(dm,&dim);
1068: if (dof0 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof0 cannot be negative");
1069: if (dof1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof1 cannot be negative");
1070: if (dim > 1 && dof2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof2 cannot be negative");
1071: if (dim > 2 && dof3 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof3 cannot be negative");
1072: stag->dof[0] = dof0;
1073: stag->dof[1] = dof1;
1074: if (dim > 1) stag->dof[2] = dof2;
1075: if (dim > 2) stag->dof[3] = dof3;
1076: return(0);
1077: }
1079: /*@C
1080: DMStagSetNumRanks - set ranks in each direction in the global rank grid
1082: Logically Collective; nRanks0, nRanks1, and nRanks2 must contain common values
1084: Input Parameters:
1085: + dm - the DMStag object
1086: - nRanks0,nRanks1,nRanks2 - number of ranks in each direction
1088: Notes:
1089: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1091: Level: developer
1093: .seealso: DMSTAG, DMDASetNumProcs()
1094: @*/
1095: PetscErrorCode DMStagSetNumRanks(DM dm,PetscInt nRanks0,PetscInt nRanks1,PetscInt nRanks2)
1096: {
1097: PetscErrorCode ierr;
1098: DM_Stag * const stag = (DM_Stag*)dm->data;
1099: PetscInt dim;
1106: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1107: DMGetDimension(dm,&dim);
1108: 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");
1109: 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");
1110: 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");
1111: if (nRanks0) stag->nRanks[0] = nRanks0;
1112: if (dim > 1 && nRanks1) stag->nRanks[1] = nRanks1;
1113: if (dim > 2 && nRanks2) stag->nRanks[2] = nRanks2;
1114: return(0);
1115: }
1117: /*@C
1118: DMStagSetStencilType - set elementwise ghost/halo stencil type
1120: Logically Collective; stencilType must contain common value
1122: Input Parameters:
1123: + dm - the DMStag object
1124: - stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE
1126: Level: beginner
1128: .seealso: DMSTAG, DMStagGetStencilType(), DMStagStencilType, DMDASetStencilType()
1129: @*/
1130: PetscErrorCode DMStagSetStencilType(DM dm,DMStagStencilType stencilType)
1131: {
1132: DM_Stag * const stag = (DM_Stag*)dm->data;
1137: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1138: stag->stencilType = stencilType;
1139: return(0);
1140: }
1142: /*@C
1143: DMStagSetStencilWidth - set elementwise stencil width
1145: Logically Collective; stencilWidth must contain common value
1147: Input Parameters:
1148: + dm - the DMStag object
1149: - stencilWidth - stencil/halo/ghost width in elements
1151: Level: beginner
1153: .seealso: DMSTAG, DMDASetStencilWidth()
1154: @*/
1155: PetscErrorCode DMStagSetStencilWidth(DM dm,PetscInt stencilWidth)
1156: {
1157: DM_Stag * const stag = (DM_Stag*)dm->data;
1162: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1163: if (stencilWidth < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Stencil width must be non-negative");
1164: stag->stencilWidth = stencilWidth;
1165: return(0);
1166: }
1168: /*@C
1169: DMStagSetGlobalSizes - set global element counts in each direction
1171: Logically Collective; N0, N1, and N2 must contain common values
1173: Input Parameters:
1174: + dm - the DMStag object
1175: - N0,N1,N2 - global elementwise sizes
1177: Notes:
1178: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1180: Level: advanced
1182: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMDASetSizes()
1183: @*/
1184: PetscErrorCode DMStagSetGlobalSizes(DM dm,PetscInt N0,PetscInt N1,PetscInt N2)
1185: {
1186: PetscErrorCode ierr;
1187: DM_Stag * const stag = (DM_Stag*)dm->data;
1188: PetscInt dim;
1192: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1193: DMGetDimension(dm,&dim);
1194: if (N0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in X direction must be positive");
1195: if (dim > 1 && N1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Y direction must be positive");
1196: if (dim > 2 && N2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Z direction must be positive");
1197: if (N0) stag->N[0] = N0;
1198: if (N1) stag->N[1] = N1;
1199: if (N2) stag->N[2] = N2;
1200: return(0);
1201: }
1203: /*@C
1204: DMStagSetOwnershipRanges - set elements per rank in each direction
1206: Logically Collective; lx, ly, and lz must contain common values
1208: Input Parameters:
1209: + dm - the DMStag object
1210: - lx,ly,lz - element counts for each rank in each direction
1212: Notes:
1213: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
1215: Level: developer
1217: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagGetOwnershipRanges(), DMDASetOwnershipRanges()
1218: @*/
1219: PetscErrorCode DMStagSetOwnershipRanges(DM dm,PetscInt const *lx,PetscInt const *ly,PetscInt const *lz)
1220: {
1221: PetscErrorCode ierr;
1222: DM_Stag * const stag = (DM_Stag*)dm->data;
1223: const PetscInt *lin[3];
1224: PetscInt d,dim;
1228: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1229: lin[0] = lx; lin[1] = ly; lin[2] = lz;
1230: DMGetDimension(dm,&dim);
1231: for (d=0; d<dim; ++d) {
1232: if (lin[d]) {
1233: if (stag->nRanks[d] < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of ranks");
1234: if (!stag->l[d]) {
1235: PetscMalloc1(stag->nRanks[d], &stag->l[d]);
1236: }
1237: PetscArraycpy(stag->l[d], lin[d], stag->nRanks[d]);
1238: }
1239: }
1240: return(0);
1241: }
1243: /*@C
1244: DMStagSetUniformCoordinates - set DMStag coordinates to be a uniform grid
1246: Collective
1248: Input Parameters:
1249: + dm - the DMStag object
1250: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1252: Notes:
1253: DMStag supports 2 different types of coordinate DM: DMSTAG and DMPRODUCT.
1254: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1256: Local coordinates are populated (using DMSetCoordinatesLocal()), linearly
1257: extrapolated to ghost cells, including those outside the physical domain.
1258: This is also done in case of periodic boundaries, meaning that the same
1259: global point may have different coordinates in different local representations,
1260: which are equivalent assuming a periodicity implied by the arguments to this function,
1261: i.e. two points are equivalent if their difference is a multiple of (xmax - xmin)
1262: in the x direction, (ymax - ymin) in the y direction, and (zmax - zmin) in the z direction.
1264: Level: advanced
1266: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType(), DMGetCoordinateDM(), DMGetCoordinates(), DMDASetUniformCoordinates(), DMBoundaryType
1267: @*/
1268: PetscErrorCode DMStagSetUniformCoordinates(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1269: {
1270: PetscErrorCode ierr;
1271: DM_Stag * const stag = (DM_Stag*)dm->data;
1272: PetscBool flg_stag,flg_product;
1276: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1277: if (!stag->coordinateDMType) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"You must first call DMStagSetCoordinateDMType()");
1278: PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg_stag);
1279: PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg_product);
1280: if (flg_stag) {
1281: DMStagSetUniformCoordinatesExplicit(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1282: } else if (flg_product) {
1283: DMStagSetUniformCoordinatesProduct(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1284: } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported DM Type %s",stag->coordinateDMType);
1285: return(0);
1286: }
1288: /*@C
1289: DMStagSetUniformCoordinatesExplicit - set DMStag coordinates to be a uniform grid, storing all values
1291: Collective
1293: Input Parameters:
1294: + dm - the DMStag object
1295: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1297: Notes:
1298: DMStag supports 2 different types of coordinate DM: either another DMStag, or a DMProduct.
1299: If the grid is orthogonal, using DMProduct should be more efficient.
1301: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1303: See the manual page for DMStagSetUniformCoordinates() for information on how
1304: coordinates for dummy cells outside the physical domain boundary are populated.
1306: Level: beginner
1308: .seealso: DMSTAG, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType()
1309: @*/
1310: PetscErrorCode DMStagSetUniformCoordinatesExplicit(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1311: {
1312: PetscErrorCode ierr;
1313: DM_Stag * const stag = (DM_Stag*)dm->data;
1314: PetscInt dim;
1315: PetscBool flg;
1319: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1320: PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg);
1321: if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1322: DMStagSetCoordinateDMType(dm,DMSTAG);
1323: DMGetDimension(dm,&dim);
1324: switch (dim) {
1325: case 1: DMStagSetUniformCoordinatesExplicit_1d(dm,xmin,xmax); break;
1326: case 2: DMStagSetUniformCoordinatesExplicit_2d(dm,xmin,xmax,ymin,ymax); break;
1327: case 3: DMStagSetUniformCoordinatesExplicit_3d(dm,xmin,xmax,ymin,ymax,zmin,zmax); break;
1328: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1329: }
1330: return(0);
1331: }
1333: /*@C
1334: DMStagSetUniformCoordinatesProduct - create uniform coordinates, as a product of 1D arrays
1336: 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.
1338: Collective
1340: Input Parameters:
1341: + dm - the DMStag object
1342: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1344: Notes:
1345: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1347: The per-dimension 1-dimensional DMStag objects that comprise the product
1348: always have active 0-cells (vertices, element boundaries) and 1-cells
1349: (element centers).
1351: See the manual page for DMStagSetUniformCoordinates() for information on how
1352: coordinates for dummy cells outside the physical domain boundary are populated.
1354: Level: intermediate
1356: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetCoordinateDMType()
1357: @*/
1358: PetscErrorCode DMStagSetUniformCoordinatesProduct(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1359: {
1360: PetscErrorCode ierr;
1361: DM_Stag * const stag = (DM_Stag*)dm->data;
1362: DM dmc;
1363: PetscInt dim,d,dof0,dof1;
1364: PetscBool flg;
1368: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1369: PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg);
1370: if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1371: DMStagSetCoordinateDMType(dm,DMPRODUCT);
1372: DMGetCoordinateDM(dm,&dmc);
1373: DMGetDimension(dm,&dim);
1375: /* Create 1D sub-DMs, living on subcommunicators.
1376: Always include both vertex and element dof, regardless of the active strata of the DMStag */
1377: dof0 = 1;
1378: dof1 = 1;
1380: for (d=0; d<dim; ++d) {
1381: DM subdm;
1382: MPI_Comm subcomm;
1383: PetscMPIInt color;
1384: const PetscMPIInt key = 0; /* let existing rank break ties */
1386: /* Choose colors based on position in the plane orthogonal to this dim, and split */
1387: switch (d) {
1388: case 0: color = (dim > 1 ? stag->rank[1] : 0) + (dim > 2 ? stag->nRanks[1]*stag->rank[2] : 0); break;
1389: case 1: color = stag->rank[0] + (dim > 2 ? stag->nRanks[0]*stag->rank[2] : 0); break;
1390: case 2: color = stag->rank[0] + stag->nRanks[0]*stag->rank[1] ; break;
1391: default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1392: }
1393: MPI_Comm_split(PetscObjectComm((PetscObject)dm),color,key,&subcomm);
1395: /* Create sub-DMs living on these new communicators (which are destroyed by DMProduct) */
1396: DMStagCreate1d(subcomm,stag->boundaryType[d],stag->N[d],dof0,dof1,stag->stencilType,stag->stencilWidth,stag->l[d],&subdm);
1397: DMSetUp(subdm);
1398: switch (d) {
1399: case 0:
1400: DMStagSetUniformCoordinatesExplicit(subdm,xmin,xmax,0,0,0,0);
1401: break;
1402: case 1:
1403: DMStagSetUniformCoordinatesExplicit(subdm,ymin,ymax,0,0,0,0);
1404: break;
1405: case 2:
1406: DMStagSetUniformCoordinatesExplicit(subdm,zmin,zmax,0,0,0,0);
1407: break;
1408: default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1409: }
1410: DMProductSetDM(dmc,d,subdm);
1411: DMProductSetDimensionIndex(dmc,d,0);
1412: DMDestroy(&subdm);
1413: MPI_Comm_free(&subcomm);
1414: }
1415: return(0);
1416: }
1418: /*@C
1419: DMStagVecGetArray - get access to local array
1421: Logically Collective
1423: This function returns a (dim+1)-dimensional array for a dim-dimensional
1424: DMStag.
1426: The first 1-3 dimensions indicate an element in the global
1427: numbering, using the standard C ordering.
1429: The final dimension in this array corresponds to a degree
1430: of freedom with respect to this element, for example corresponding to
1431: the element or one of its neighboring faces, edges, or vertices.
1433: For example, for a 3D DMStag, indexing is array[k][j][i][idx], where k is the
1434: index in the z-direction, j is the index in the y-direction, and i is the
1435: index in the x-direction.
1437: "idx" is obtained with DMStagGetLocationSlot(), since the correct offset
1438: into the (dim+1)-dimensional C array depends on the grid size and the number
1439: of dof stored at each location.
1441: Input Parameters:
1442: + dm - the DMStag object
1443: - vec - the Vec object
1445: Output Parameters:
1446: . array - the array
1448: Notes:
1449: DMStagVecRestoreArray() must be called, once finished with the array
1451: Level: beginner
1453: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArray(), DMDAVecGetArrayDOF()
1454: @*/
1455: PetscErrorCode DMStagVecGetArray(DM dm,Vec vec,void *array)
1456: {
1457: PetscErrorCode ierr;
1458: DM_Stag * const stag = (DM_Stag*)dm->data;
1459: PetscInt dim;
1460: PetscInt nLocal;
1465: DMGetDimension(dm,&dim);
1466: VecGetLocalSize(vec,&nLocal);
1467: 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);
1468: switch (dim) {
1469: case 1:
1470: VecGetArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1471: break;
1472: case 2:
1473: VecGetArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1474: break;
1475: case 3:
1476: 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);
1477: break;
1478: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1479: }
1480: return(0);
1481: }
1483: /*@C
1484: DMStagVecGetArrayRead - get read-only access to a local array
1486: Logically Collective
1488: See the man page for DMStagVecGetArray() for more information.
1490: Input Parameters:
1491: + dm - the DMStag object
1492: - vec - the Vec object
1494: Output Parameters:
1495: . array - the read-only array
1497: Notes:
1498: DMStagVecRestoreArrayRead() must be called, once finished with the array
1500: Level: beginner
1502: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArrayRead(), DMDAVecGetArrayDOFRead()
1503: @*/
1504: PetscErrorCode DMStagVecGetArrayRead(DM dm,Vec vec,void *array)
1505: {
1506: PetscErrorCode ierr;
1507: DM_Stag * const stag = (DM_Stag*)dm->data;
1508: PetscInt dim;
1509: PetscInt nLocal;
1514: DMGetDimension(dm,&dim);
1515: VecGetLocalSize(vec,&nLocal);
1516: 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);
1517: switch (dim) {
1518: case 1:
1519: VecGetArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1520: break;
1521: case 2:
1522: VecGetArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1523: break;
1524: case 3:
1525: 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);
1526: break;
1527: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1528: }
1529: return(0);
1530: }
1532: /*@C
1533: DMStagVecRestoreArray - restore access to a raw array
1535: Logically Collective
1537: Input Parameters:
1538: + dm - the DMStag object
1539: - vec - the Vec object
1541: Output Parameters:
1542: . array - the array
1544: Level: beginner
1546: .seealso: DMSTAG, DMStagVecGetArray(), DMDAVecRestoreArray(), DMDAVecRestoreArrayDOF()
1547: @*/
1548: PetscErrorCode DMStagVecRestoreArray(DM dm,Vec vec,void *array)
1549: {
1550: PetscErrorCode ierr;
1551: DM_Stag * const stag = (DM_Stag*)dm->data;
1552: PetscInt dim;
1553: PetscInt nLocal;
1558: DMGetDimension(dm,&dim);
1559: VecGetLocalSize(vec,&nLocal);
1560: 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);
1561: switch (dim) {
1562: case 1:
1563: VecRestoreArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1564: break;
1565: case 2:
1566: VecRestoreArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1567: break;
1568: case 3:
1569: 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);
1570: break;
1571: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1572: }
1573: return(0);
1574: }
1576: /*@C
1577: DMStagVecRestoreArrayRead - restore read-only access to a raw array
1579: Logically Collective
1581: Input Parameters:
1582: + dm - the DMStag object
1583: - vec - the Vec object
1585: Output Parameters:
1586: . array - the read-only array
1588: Level: beginner
1590: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMDAVecRestoreArrayRead(), DMDAVecRestoreArrayDOFRead()
1591: @*/
1592: PetscErrorCode DMStagVecRestoreArrayRead(DM dm,Vec vec,void *array)
1593: {
1594: PetscErrorCode ierr;
1595: DM_Stag * const stag = (DM_Stag*)dm->data;
1596: PetscInt dim;
1597: PetscInt nLocal;
1602: DMGetDimension(dm,&dim);
1603: VecGetLocalSize(vec,&nLocal);
1604: 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);
1605: switch (dim) {
1606: case 1:
1607: VecRestoreArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1608: break;
1609: case 2:
1610: VecRestoreArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1611: break;
1612: case 3:
1613: 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);
1614: break;
1615: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1616: }
1617: return(0);
1618: }