Actual source code: stagutils.c
petsc-3.12.5 2020-03-29
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, DMDAGetBoundaryTypes()
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: /*@C
35: DMStagGet1dCoordinateArraysDOFRead - extract 1D coordinate arrays
37: Logically Collective
39: A high-level helper function to quickly extract raw 1D local coordinate arrays.
40: Checks that the coordinate DM is a DMProduct or 1D DMStags, with the same number of dof.
41: Check on the number of dof and dimension ensures that the elementwise data
42: is the same for each, so the same indexing can be used on the arrays.
44: Input Parameter:
45: . dm - the DMStag object
47: Output Parameters:
48: . arrX,arrY,arrX - local 1D coordinate arrays
50: Level: intermediate
52: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagGet1dCoordinateLocationSlot()
53: @*/
54: PetscErrorCode DMStagGet1dCoordinateArraysDOFRead(DM dm,void* arrX,void* arrY,void* arrZ)
55: {
57: PetscInt dim,d,dofCheck[DMSTAG_MAX_STRATA],s;
58: DM dmCoord;
59: void* arr[DMSTAG_MAX_DIM];
63: DMGetDimension(dm,&dim);
64: if (dim > DMSTAG_MAX_DIM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for %D dimensions",dim);
65: arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
66: DMGetCoordinateDM(dm,&dmCoord);
67: if (!dmCoord) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM does not have a coordinate DM");
68: {
69: PetscBool isProduct;
70: DMType dmType;
71: DMGetType(dmCoord,&dmType);
72: PetscStrcmp(DMPRODUCT,dmType,&isProduct);
73: if (!isProduct) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is not of type DMPRODUCT");
74: }
75: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
76: for (d=0; d<dim; ++d) {
77: DM subDM;
78: DMType dmType;
79: PetscBool isStag;
80: PetscInt dof[DMSTAG_MAX_STRATA],subDim;
81: Vec coord1d;
82: DMProductGetDM(dmCoord,d,&subDM);
83: if (!subDM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is missing sub DM %D",d);
84: DMGetDimension(subDM,&subDim);
85: if (subDim != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of dimension 1");
86: DMGetType(subDM,&dmType);
87: PetscStrcmp(DMSTAG,dmType,&isStag);
88: if (!isStag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of type DMSTAG");
89: DMStagGetDOF(subDM,&dof[0],&dof[1],&dof[2],&dof[3]);
90: if (d == 0) {
91: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
92: } else {
93: for (s=0; s<DMSTAG_MAX_STRATA; ++s) {
94: if (dofCheck[s] != dof[s]) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DMs have different dofs");
95: }
96: }
97: DMGetCoordinatesLocal(subDM,&coord1d);
98: DMStagVecGetArrayDOFRead(subDM,coord1d,arr[d]);
99: }
100: return(0);
101: }
103: /*@C
104: DMStagGet1dCoordinateLocationSlot - get slot for use with local 1D coordinate arrays
106: High-level helper function to get slot ids for 1D coordinate DMs.
107: For use with DMStagGetIDCoordinateArraysDOFRead() and related functions.
109: Not Collective
111: Input Parameters:
112: + dm - the DMStag object
113: - loc - the grid location
115: Output Parameter:
116: . slot - the index to use in local arrays
118: Notes:
119: Checks that the coordinates are actually set up so that using the
120: slots from the first 1d coordinate sub-DM is valid for all the 1D coordinate sub-DMs.
122: Level: intermediate
124: .seealso: DMSTAG, DMPRODUCT, DMStagGet1dCoordinateArraysDOFRead(), DMStagSetUniformCoordinates()
125: @*/
126: PETSC_EXTERN PetscErrorCode DMStagGet1dCoordinateLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt *slot)
127: {
129: DM dmCoord;
130: PetscInt dim,dofCheck[DMSTAG_MAX_STRATA],s,d;
134: DMGetDimension(dm,&dim);
135: DMGetCoordinateDM(dm,&dmCoord);
136: if (!dmCoord) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM does not have a coordinate DM");
137: {
138: PetscBool isProduct;
139: DMType dmType;
140: DMGetType(dmCoord,&dmType);
141: PetscStrcmp(DMPRODUCT,dmType,&isProduct);
142: if (!isProduct) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is not of type DMPRODUCT");
143: }
144: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
145: for (d=0; d<dim; ++d) {
146: DM subDM;
147: DMType dmType;
148: PetscBool isStag;
149: PetscInt dof[DMSTAG_MAX_STRATA],subDim;
150: DMProductGetDM(dmCoord,d,&subDM);
151: if (!subDM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is missing sub DM %D",d);
152: DMGetDimension(subDM,&subDim);
153: if (subDim != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of dimension 1");
154: DMGetType(subDM,&dmType);
155: PetscStrcmp(DMSTAG,dmType,&isStag);
156: if (!isStag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of type DMSTAG");
157: DMStagGetDOF(subDM,&dof[0],&dof[1],&dof[2],&dof[3]);
158: if (d == 0) {
159: const PetscInt component = 0;
160: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
161: DMStagGetLocationSlot(subDM,loc,component,slot);
162: } else {
163: for (s=0; s<DMSTAG_MAX_STRATA; ++s) {
164: if (dofCheck[s] != dof[s]) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DMs have different dofs");
165: }
166: }
167: }
168: return(0);
169: }
171: /*@C
172: DMStagGetCorners - return global element indices of the local region (excluding ghost points)
174: Not Collective
176: Input Parameter:
177: . dm - the DMStag object
179: Output Parameters:
180: + x,y,z - starting element indices in each direction
181: . m,n,p - element widths in each direction
182: - nExtrax,nExtray,nExtraz - number of extra partial elements in each direction.
184: Notes:
185: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
187: The number of extra partial elements is either 1 or 0.
188: The value is 1 on right, top, and front non-periodic domain ("physical") boundaries,
189: in the x, y, and z directions respectively, and otherwise 0.
191: Level: beginner
193: .seealso: DMSTAG, DMStagGetGhostCorners(), DMDAGetCorners()
194: @*/
195: PetscErrorCode DMStagGetCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *nExtrax,PetscInt *nExtray,PetscInt *nExtraz)
196: {
197: const DM_Stag * const stag = (DM_Stag*)dm->data;
201: if (x) *x = stag->start[0];
202: if (y) *y = stag->start[1];
203: if (z) *z = stag->start[2];
204: if (m) *m = stag->n[0];
205: if (n) *n = stag->n[1];
206: if (p) *p = stag->n[2];
207: if (nExtrax) *nExtrax = stag->boundaryType[0] != DM_BOUNDARY_PERIODIC && stag->lastRank[0] ? 1 : 0;
208: if (nExtray) *nExtray = stag->boundaryType[1] != DM_BOUNDARY_PERIODIC && stag->lastRank[1] ? 1 : 0;
209: if (nExtraz) *nExtraz = stag->boundaryType[2] != DM_BOUNDARY_PERIODIC && stag->lastRank[2] ? 1 : 0;
210: return(0);
211: }
213: /*@C
214: DMStagGetDOF - get number of DOF associated with each stratum of the grid
216: Not Collective
218: Input Parameter:
219: . dm - the DMStag object
221: Output Parameters:
222: + dof0 - the number of points per 0-cell (vertex/node)
223: . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
224: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
225: - dof3 - the number of points per 3-cell (element in 3D)
227: Level: beginner
229: .seealso: DMSTAG, DMStagGetCorners(), DMStagGetGhostCorners(), DMStagGetGlobalSizes(), DMStagGetStencilWidth(), DMStagGetBoundaryTypes(), DMStagGetLocationDof(), DMDAGetDof()
230: @*/
231: PetscErrorCode DMStagGetDOF(DM dm,PetscInt *dof0,PetscInt *dof1,PetscInt *dof2,PetscInt *dof3)
232: {
233: const DM_Stag * const stag = (DM_Stag*)dm->data;
237: if (dof0) *dof0 = stag->dof[0];
238: if (dof1) *dof1 = stag->dof[1];
239: if (dof2) *dof2 = stag->dof[2];
240: if (dof3) *dof3 = stag->dof[3];
241: return(0);
242: }
244: /*@C
245: DMStagGetGhostCorners - return global element indices of the local region, including ghost points
247: Not Collective
249: Input Argument:
250: . dm - the DMStag object
252: Output Arguments:
253: + x,y,z - starting element indices in each direction
254: - m,n,p - element widths in each direction
256: Notes:
257: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
259: Level: beginner
261: .seealso: DMSTAG, DMStagGetCorners(), DMDAGetGhostCorners()
262: @*/
263: PetscErrorCode DMStagGetGhostCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
264: {
265: const DM_Stag * const stag = (DM_Stag*)dm->data;
269: if (x) *x = stag->startGhost[0];
270: if (y) *y = stag->startGhost[1];
271: if (z) *z = stag->startGhost[2];
272: if (m) *m = stag->nGhost[0];
273: if (n) *n = stag->nGhost[1];
274: if (p) *p = stag->nGhost[2];
275: return(0);
276: }
278: /*@C
279: DMStagGetGlobalSizes - get global element counts
281: Not Collective
283: Input Parameter:
284: . dm - the DMStag object
286: Output Parameters:
287: . M,N,P - global element counts in each direction
289: Notes:
290: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
292: Level: beginner
294: .seealso: DMSTAG, DMStagGetLocalSizes(), DMDAGetInfo()
295: @*/
296: PetscErrorCode DMStagGetGlobalSizes(DM dm,PetscInt* M,PetscInt* N,PetscInt* P)
297: {
298: const DM_Stag * const stag = (DM_Stag*)dm->data;
302: if (M) *M = stag->N[0];
303: if (N) *N = stag->N[1];
304: if (P) *P = stag->N[2];
305: return(0);
306: }
308: /*@C
309: DMStagGetIsFirstRank - get boolean value for whether this rank is first in each direction in the rank grid
311: Not Collective
313: Input Parameter:
314: . dm - the DMStag object
316: Output Parameters:
317: . isFirstRank0,isFirstRank1,isFirstRank2 - whether this rank is first in each direction
319: Level: intermediate
321: Notes:
322: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
324: .seealso: DMSTAG, DMStagGetIsLastRank()
325: @*/
326: PetscErrorCode DMStagGetIsFirstRank(DM dm,PetscBool *isFirstRank0,PetscBool *isFirstRank1,PetscBool *isFirstRank2)
327: {
328: const DM_Stag * const stag = (DM_Stag*)dm->data;
332: if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
333: if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
334: if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
335: return(0);
336: }
338: /*@C
339: DMStagGetIsLastRank - get boolean value for whether this rank is last in each direction in the rank grid
341: Not Collective
343: Input Parameter:
344: . dm - the DMStag object
346: Output Parameters:
347: . isLastRank0,isLastRank1,isLastRank2 - whether this rank is last in each direction
349: Level: intermediate
351: Notes:
352: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
353: Level: intermediate
355: .seealso: DMSTAG, DMStagGetIsFirstRank()
356: @*/
357: PetscErrorCode DMStagGetIsLastRank(DM dm,PetscBool *isLastRank0,PetscBool *isLastRank1,PetscBool *isLastRank2)
358: {
359: const DM_Stag * const stag = (DM_Stag*)dm->data;
363: if (isLastRank0) *isLastRank0 = stag->lastRank[0];
364: if (isLastRank1) *isLastRank1 = stag->lastRank[1];
365: if (isLastRank2) *isLastRank2 = stag->lastRank[2];
366: return(0);
367: }
369: /*@C
370: DMStagGetLocalSizes - get local elementwise sizes
372: Not Collective
374: Input Parameter:
375: . dm - the DMStag object
377: Output Parameters:
378: . m,n,p - local element counts (excluding ghosts) in each direction
380: Notes:
381: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
382: Level: intermediate
384: Level: beginner
386: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetDOF(), DMStagGetNumRanks(), DMDAGetLocalInfo()
387: @*/
388: PetscErrorCode DMStagGetLocalSizes(DM dm,PetscInt* m,PetscInt* n,PetscInt* p)
389: {
390: const DM_Stag * const stag = (DM_Stag*)dm->data;
394: if (m) *m = stag->n[0];
395: if (n) *n = stag->n[1];
396: if (p) *p = stag->n[2];
397: return(0);
398: }
400: /*@C
401: DMStagGetNumRanks - get number of ranks in each direction in the global grid decomposition
403: Not Collective
405: Input Parameter:
406: . dm - the DMStag object
408: Output Parameters:
409: . nRanks0,nRanks1,nRanks2 - number of ranks in each direction in the grid decomposition
411: Notes:
412: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
413: Level: intermediate
415: Level: beginner
417: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetLocalSize(), DMStagSetNumRank(), DMDAGetInfo()
418: @*/
419: PetscErrorCode DMStagGetNumRanks(DM dm,PetscInt *nRanks0,PetscInt *nRanks1,PetscInt *nRanks2)
420: {
421: const DM_Stag * const stag = (DM_Stag*)dm->data;
425: if (nRanks0) *nRanks0 = stag->nRanks[0];
426: if (nRanks1) *nRanks1 = stag->nRanks[1];
427: if (nRanks2) *nRanks2 = stag->nRanks[2];
428: return(0);
429: }
431: /*@C
432: DMStagGetEntriesPerElement - get number of entries per element in the local representation
434: Not Collective
436: Input Parameter:
437: . dm - the DMStag object
439: Output Parameters:
440: . entriesPerElement - number of entries associated with each element in the local representation
442: Notes:
443: This is the natural block size for most local operations. In 1D it is equal to dof0 + dof1,
444: in 2D it is equal to dof0 + 2*dof1 + dof2, and in 3D it is equal to dof0 + 3*dof1 + 3*dof2 + dof3
446: Level: developer
448: .seealso: DMSTAG, DMStagGetDOF()
449: @*/
450: PetscErrorCode DMStagGetEntriesPerElement(DM dm,PetscInt *entriesPerElement)
451: {
452: const DM_Stag * const stag = (DM_Stag*)dm->data;
456: if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
457: return(0);
458: }
460: /*@C
461: DMStagGetStencilType - get elementwise ghost/halo stencil type
463: Not Collective
465: Input Parameter:
466: . dm - the DMStag object
468: Output Parameter:
469: . stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE
471: Level: beginner
473: .seealso: DMSTAG, DMStagSetStencilType(), DMStagStencilType, DMDAGetInfo()
474: @*/
475: PetscErrorCode DMStagGetStencilType(DM dm,DMStagStencilType *stencilType)
476: {
477: DM_Stag * const stag = (DM_Stag*)dm->data;
481: *stencilType = stag->stencilType;
482: return(0);
483: }
485: /*@C
486: DMStagGetStencilWidth - get elementwise stencil width
488: Not Collective
490: Input Parameter:
491: . dm - the DMStag object
493: Output Parameters:
494: . stencilWidth - stencil/halo/ghost width in elements
496: Level: beginner
498: .seealso: DMSTAG, DMDAGetStencilWidth(), DMDAGetInfo()
499: @*/
500: PetscErrorCode DMStagGetStencilWidth(DM dm,PetscInt *stencilWidth)
501: {
502: const DM_Stag * const stag = (DM_Stag*)dm->data;
506: if (stencilWidth) *stencilWidth = stag->stencilWidth;
507: return(0);
508: }
510: /*@C
511: DMStagGetOwnershipRanges - get elements per rank in each direction
513: Not Collective
515: Input Parameter:
516: . dm - the DMStag object
518: Output Parameters:
519: + lx - ownership along x direction (optional)
520: . ly - ownership along y direction (optional)
521: - lz - ownership along z direction (optional)
523: Notes:
524: These correspond to the optional final arguments passed to DMStagCreate1d(), DMStagCreate2d(), and DMStagCreate3d().
526: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
528: In C you should not free these arrays, nor change the values in them.
529: They will only have valid values while the DMStag they came from still exists (has not been destroyed).
531: Level: intermediate
533: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagSetOwnershipRanges(), DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDAGetOwnershipRanges()
534: @*/
535: PetscErrorCode DMStagGetOwnershipRanges(DM dm,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
536: {
537: const DM_Stag * const stag = (DM_Stag*)dm->data;
541: if (lx) *lx = stag->l[0];
542: if (ly) *ly = stag->l[1];
543: if (lz) *lz = stag->l[2];
544: return(0);
545: }
547: /*@C
548: DMStagCreateCompatibleDMStag - create a compatible DMStag with different dof/stratum
550: Collective
552: Input Parameters
553: + dm - the DMStag object
554: - dof0,dof1,dof2,dof3 - number of dof on each stratum in the new DMStag
556: Output Parameters:
557: . newdm - the new, compatible DMStag
559: Notes:
560: Dof supplied for strata too big for the dimension are ignored; these may be set to 0.
561: In contrast to DMDACreateCompatibleDMDA(), coordinates are not reused.
563: Level: intermediate
565: .seealso: DMSTAG, DMDACreateCompatibleDMDA(), DMGetCompatibility(), DMStagMigrateVec()
566: @*/
567: PetscErrorCode DMStagCreateCompatibleDMStag(DM dm,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DM *newdm)
568: {
569: PetscErrorCode ierr;
570: const DM_Stag * const stag = (DM_Stag*)dm->data;
571: PetscInt dim;
575: DMGetDimension(dm,&dim);
576: switch (dim) {
577: case 1:
578: DMStagCreate1d(PetscObjectComm((PetscObject)dm),stag->boundaryType[0],stag->N[0],dof0,dof1,stag->stencilType,stag->stencilWidth,NULL,newdm);
579: break;
580: case 2:
581: DMStagCreate2d(PetscObjectComm((PetscObject)dm),stag->boundaryType[0],stag->boundaryType[1],stag->N[0],stag->N[1],stag->nRanks[0],stag->nRanks[1],dof0,dof1,dof2,stag->stencilType,stag->stencilWidth,NULL,NULL,newdm);
582: break;
583: case 3:
584: DMStagCreate3d(PetscObjectComm((PetscObject)dm),stag->boundaryType[0],stag->boundaryType[1],stag->boundaryType[2],stag->N[0],stag->N[1],stag->N[2],stag->nRanks[0],stag->nRanks[1],stag->nRanks[2],dof0,dof1,dof2,dof3,stag->stencilType,stag->stencilWidth,NULL,NULL,NULL,newdm);
585: break;
586: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
587: }
588: DMSetUp(*newdm);
589: return(0);
590: }
592: /*@C
593: DMStagGetLocationSlot - get index to use in accessing raw local arrays
595: Not Collective
597: Input Parameters:
598: + dm - the DMStag object
599: . loc - location relative to an element
600: - c - component
602: Output Parameter:
603: . slot - index to use
605: Notes:
606: Provides an appropriate index to use with DMStagVecGetArrayDOF() and friends.
607: This is required so that the user doesn't need to know about the ordering of
608: dof associated with each local element.
610: Level: beginner
612: .seealso: DMSTAG, DMStagVecGetArrayDOF(), DMStagVecGetArrayDOFRead(), DMStagGetDOF(), DMStagGetEntriesPerElement()
613: @*/
614: PetscErrorCode DMStagGetLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt c,PetscInt *slot)
615: {
616: DM_Stag * const stag = (DM_Stag*)dm->data;
620: #if defined(PETSC_USE_DEBUG)
621: {
623: PetscInt dof;
624: DMStagGetLocationDOF(dm,loc,&dof);
625: if (dof < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Location %s has no dof attached",DMStagStencilLocations[loc]);
626: 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);
627: }
628: #endif
629: *slot = stag->locationOffsets[loc] + c;
630: return(0);
631: }
633: /*@C
634: DMStagMigrateVec - transfer a vector associated with a DMStag to a vector associated with a compatible DMStag
636: Collective
638: Input Parameters:
639: + dm - the source DMStag object
640: . vec - the source vector, compatible with dm
641: . dmTo - the compatible destination DMStag object
642: - vecTo - the destination vector, compatible with dmTo
644: Notes:
645: Extra dof are ignored, and unfilled dof are zeroed.
646: Currently only implemented to migrate global vectors to global vectors.
648: Level: advanced
650: .seealso: DMSTAG, DMStagCreateCompatibleDMStag(), DMGetCompatibility(), DMStagVecSplitToDMDA()
651: @*/
652: PetscErrorCode DMStagMigrateVec(DM dm,Vec vec,DM dmTo,Vec vecTo)
653: {
654: PetscErrorCode ierr;
655: DM_Stag * const stag = (DM_Stag*)dm->data;
656: DM_Stag * const stagTo = (DM_Stag*)dmTo->data;
657: PetscInt nLocalTo,nLocal,dim,i,j,k;
658: PetscInt start[DMSTAG_MAX_DIM],startGhost[DMSTAG_MAX_DIM],n[DMSTAG_MAX_DIM],nExtra[DMSTAG_MAX_DIM],offset[DMSTAG_MAX_DIM];
659: Vec vecToLocal,vecLocal;
660: PetscBool compatible,compatibleSet;
661: const PetscScalar *arr;
662: PetscScalar *arrTo;
663: const PetscInt epe = stag->entriesPerElement;
664: const PetscInt epeTo = stagTo->entriesPerElement;
671: DMGetCompatibility(dm,dmTo,&compatible,&compatibleSet);
672: if (!compatibleSet || !compatible) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"DMStag objects must be shown to be compatible");
673: DMGetDimension(dm,&dim);
674: VecGetLocalSize(vec,&nLocal);
675: VecGetLocalSize(vecTo,&nLocalTo);
676: if (nLocal != stag->entries|| nLocalTo !=stagTo->entries) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Vector migration only implemented for global vector to global vector.");
677: DMStagGetCorners(dm,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
678: DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],NULL,NULL,NULL);
679: for (i=0; i<DMSTAG_MAX_DIM; ++i) offset[i] = start[i]-startGhost[i];
681: /* Proceed by transferring to a local vector, copying, and transferring back to a global vector */
682: DMGetLocalVector(dm,&vecLocal);
683: DMGetLocalVector(dmTo,&vecToLocal);
684: DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal);
685: DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal);
686: VecGetArrayRead(vecLocal,&arr);
687: VecGetArray(vecToLocal,&arrTo);
688: /* Note that some superfluous copying of entries on partial dummy elements is done */
689: if (dim == 1) {
690: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
691: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
692: PetscInt si;
693: for (si=0; si<2; ++si) {
694: b += stag->dof[si];
695: bTo += stagTo->dof[si];
696: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[i*epeTo + dTo] = arr[i*epe + d];
697: for (; dTo < bTo ; ++dTo) arrTo[i*epeTo + dTo] = 0.0;
698: d = b;
699: }
700: }
701: } else if (dim == 2) {
702: const PetscInt epr = stag->nGhost[0] * epe;
703: const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
704: for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
705: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
706: const PetscInt base = j*epr + i*epe;
707: const PetscInt baseTo = j*eprTo + i*epeTo;
708: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
709: const PetscInt s[4] = {0,1,1,2}; /* Dimensions of points, in order */
710: PetscInt si;
711: for (si=0; si<4; ++si) {
712: b += stag->dof[s[si]];
713: bTo += stagTo->dof[s[si]];
714: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
715: for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
716: d = b;
717: }
718: }
719: }
720: } else if (dim == 3) {
721: const PetscInt epr = stag->nGhost[0] * epe;
722: const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
723: const PetscInt epl = stag->nGhost[1] * epr;
724: const PetscInt eplTo = stagTo->nGhost[1] * eprTo;
725: for (k=offset[2]; k<offset[2] + n[2] + nExtra[2]; ++k) {
726: for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
727: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
728: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
729: const PetscInt base = k*epl + j*epr + i*epe;
730: const PetscInt baseTo = k*eplTo + j*eprTo + i*epeTo;
731: const PetscInt s[8] = {0,1,1,2,1,2,2,3}; /* dimensions of points, in order */
732: PetscInt is;
733: for (is=0; is<8; ++is) {
734: b += stag->dof[s[is]];
735: bTo += stagTo->dof[s[is]];
736: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
737: for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
738: d = b;
739: }
740: }
741: }
742: }
743: } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
744: VecRestoreArrayRead(vecLocal,&arr);
745: VecRestoreArray(vecToLocal,&arrTo);
746: DMRestoreLocalVector(dm,&vecLocal);
747: DMLocalToGlobalBegin(dmTo,vecToLocal,INSERT_VALUES,vecTo);
748: DMLocalToGlobalEnd(dmTo,vecToLocal,INSERT_VALUES,vecTo);
749: DMRestoreLocalVector(dmTo,&vecToLocal);
750: return(0);
751: }
753: /*@C
754: DMStagPopulateLocalToGlobalInjective - populate an internal 1-to-1 local-to-global map
756: Collective
758: Creates an internal object which explicitly maps a single local degree of
759: freedom to each global degree of freedom. This is used, if populated,
760: instead of SCATTER_REVERSE_LOCAL with the (1-to-many, in general)
761: global-to-local map, when DMLocalToGlobal() is called with INSERT_VALUES.
762: This allows usage, for example, even in the periodic, 1-rank case, where
763: the inverse of the global-to-local map, even when restricted to on-rank
764: communication, is non-injective. This is at the cost of storing an additional
765: VecScatter object inside each DMStag object.
767: Input Parameter:
768: . dm - the DMStag object
770: Notes:
771: In normal usage, library users shouldn't be concerned with this function,
772: as it is called during DMSetUp(), when required.
774: Returns immediately if the internal map is already populated.
776: Developer Notes:
777: This could, if desired, be moved up to a general DM routine. It would allow,
778: for example, DMDA to support DMLocalToGlobal() with INSERT_VALUES,
779: even in the single-rank periodic case.
781: Level: developer
783: .seealso: DMSTAG, DMLocalToGlobal(), VecScatter
784: @*/
785: PetscErrorCode DMStagPopulateLocalToGlobalInjective(DM dm)
786: {
787: PetscErrorCode ierr;
788: PetscInt dim;
789: DM_Stag * const stag = (DM_Stag*)dm->data;
793: if (stag->ltog_injective) return(0); /* Don't re-populate */
794: DMGetDimension(dm,&dim);
795: switch (dim) {
796: case 1: DMStagPopulateLocalToGlobalInjective_1d(dm); break;
797: case 2: DMStagPopulateLocalToGlobalInjective_2d(dm); break;
798: case 3: DMStagPopulateLocalToGlobalInjective_3d(dm); break;
799: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
800: }
801: return(0);
802: }
804: /*@C
805: DMStagRestore1dCoordinateArraysDOFRead - restore local array access
807: Logically Collective
809: Input Parameter:
810: . dm - the DMStag object
812: Output Parameters:
813: . arrX,arrY,arrX - local 1D coordinate arrays
815: Level: intermediate
817: .seealso: DMSTAG, DMStagGet1dCoordinateArraysDOFRead()
818: @*/
819: PetscErrorCode DMStagRestore1dCoordinateArraysDOFRead(DM dm,void *arrX,void *arrY,void *arrZ)
820: {
821: PetscErrorCode ierr;
822: PetscInt dim,d;
823: void* arr[DMSTAG_MAX_DIM];
824: DM dmCoord;
828: DMGetDimension(dm,&dim);
829: if (dim > DMSTAG_MAX_DIM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for %D dimensions",dim);
830: arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
831: DMGetCoordinateDM(dm,&dmCoord);
832: for (d=0; d<dim; ++d) {
833: DM subDM;
834: Vec coord1d;
835: DMProductGetDM(dmCoord,d,&subDM);
836: DMGetCoordinatesLocal(subDM,&coord1d);
837: DMStagVecRestoreArrayDOFRead(subDM,coord1d,arr[d]);
838: }
839: return(0);
840: }
842: /*@C
843: DMStagSetBoundaryTypes - set DMStag boundary types
845: Logically Collective; boundaryType0, boundaryType1, and boundaryType2 must contain common values
847: Input Parameters:
848: + dm - the DMStag object
849: - boundaryType0,boundaryType1,boundaryType2 - boundary types in each direction
851: Notes:
852: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
854: Level: advanced
856: .seealso: DMSTAG, DMBoundaryType, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDASetBoundaryType()
857: @*/
858: PetscErrorCode DMStagSetBoundaryTypes(DM dm,DMBoundaryType boundaryType0,DMBoundaryType boundaryType1,DMBoundaryType boundaryType2)
859: {
860: PetscErrorCode ierr;
861: DM_Stag * const stag = (DM_Stag*)dm->data;
862: PetscInt dim;
869: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
870: DMGetDimension(dm,&dim);
871: if (boundaryType0 ) stag->boundaryType[0] = boundaryType0;
872: if (boundaryType1 && dim > 1) stag->boundaryType[1] = boundaryType1;
873: if (boundaryType2 && dim > 2) stag->boundaryType[2] = boundaryType2;
874: return(0);
875: }
877: /*@C
878: DMStagSetCoordinateDMType - set DM type to store coordinates
880: Logically Collective; dmtype must contain common value
882: Input Parameters:
883: + dm - the DMStag object
884: - dmtype - DMtype for coordinates, either DMSTAG or DMPRODUCT
886: Level: advanced
888: .seealso: DMSTAG, DMPRODUCT, DMGetCoordinateDM(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMType
889: @*/
890: PetscErrorCode DMStagSetCoordinateDMType(DM dm,DMType dmtype)
891: {
892: PetscErrorCode ierr;
893: DM_Stag * const stag = (DM_Stag*)dm->data;
897: PetscFree(stag->coordinateDMType);
898: PetscStrallocpy(dmtype,(char**)&stag->coordinateDMType);
899: return(0);
900: }
902: /*@C
903: DMStagSetDOF - set dof/stratum
905: Logically Collective; dof0, dof1, dof2, and dof3 must contain common values
907: Input Parameters:
908: + dm - the DMStag object
909: - dof0,dof1,dof2,dof3 - dof per stratum
911: Notes:
912: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
914: Level: advanced
916: .seealso: DMSTAG, DMDASetDof()
917: @*/
918: PetscErrorCode DMStagSetDOF(DM dm,PetscInt dof0, PetscInt dof1,PetscInt dof2,PetscInt dof3)
919: {
920: PetscErrorCode ierr;
921: DM_Stag * const stag = (DM_Stag*)dm->data;
922: PetscInt dim;
930: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
931: DMGetDimension(dm,&dim);
932: if (dof0 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof0 cannot be negative");
933: if (dof1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof1 cannot be negative");
934: if (dim > 1 && dof2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof2 cannot be negative");
935: if (dim > 2 && dof3 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof3 cannot be negative");
936: stag->dof[0] = dof0;
937: stag->dof[1] = dof1;
938: if (dim > 1) stag->dof[2] = dof2;
939: if (dim > 2) stag->dof[3] = dof3;
940: return(0);
941: }
943: /*@C
944: DMStagSetNumRanks - set ranks in each direction in the global rank grid
946: Logically Collective; nRanks0, nRanks1, and nRanks2 must contain common values
948: Input Parameters:
949: + dm - the DMStag object
950: - nRanks0,nRanks1,nRanks2 - number of ranks in each direction
952: Notes:
953: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
955: Level: developer
957: .seealso: DMSTAG, DMDASetNumProcs()
958: @*/
959: PetscErrorCode DMStagSetNumRanks(DM dm,PetscInt nRanks0,PetscInt nRanks1,PetscInt nRanks2)
960: {
961: PetscErrorCode ierr;
962: DM_Stag * const stag = (DM_Stag*)dm->data;
963: PetscInt dim;
970: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
971: DMGetDimension(dm,&dim);
972: if (nRanks0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in X direction cannot be less than 1");
973: if (dim > 1 && nRanks1 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in Y direction cannot be less than 1");
974: if (dim > 2 && nRanks2 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in Z direction cannot be less than 1");
975: if (nRanks0) stag->nRanks[0] = nRanks0;
976: if (nRanks1) stag->nRanks[1] = nRanks1;
977: if (nRanks2) stag->nRanks[2] = nRanks2;
978: return(0);
979: }
981: /*@C
982: DMStagSetStencilType - set elementwise ghost/halo stencil type
984: Logically Collective; stencilType must contain common value
986: Input Parameters:
987: + dm - the DMStag object
988: - stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE
990: Level: beginner
992: .seealso: DMSTAG, DMStagGetStencilType(), DMStagStencilType, DMDASetStencilType()
993: @*/
994: PetscErrorCode DMStagSetStencilType(DM dm,DMStagStencilType stencilType)
995: {
996: DM_Stag * const stag = (DM_Stag*)dm->data;
1001: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1002: stag->stencilType = stencilType;
1003: return(0);
1004: }
1006: /*@C
1007: DMStagSetStencilWidth - set elementwise stencil width
1009: Logically Collective; stencilWidth must contain common value
1011: Input Parameters:
1012: + dm - the DMStag object
1013: - stencilWidth - stencil/halo/ghost width in elements
1015: Level: beginner
1017: .seealso: DMSTAG, DMDASetStencilWidth()
1018: @*/
1019: PetscErrorCode DMStagSetStencilWidth(DM dm,PetscInt stencilWidth)
1020: {
1021: DM_Stag * const stag = (DM_Stag*)dm->data;
1026: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1027: if (stencilWidth < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Stencil width must be non-negative");
1028: stag->stencilWidth = stencilWidth;
1029: return(0);
1030: }
1032: /*@C
1033: DMStagSetGlobalSizes - set global element counts in each direction
1035: Logically Collective; N0, N1, and N2 must contain common values
1037: Input Parameters:
1038: + dm - the DMStag object
1039: - N0,N1,N2 - global elementwise sizes
1041: Notes:
1042: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1044: Level: advanced
1046: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMDASetSizes()
1047: @*/
1048: PetscErrorCode DMStagSetGlobalSizes(DM dm,PetscInt N0,PetscInt N1,PetscInt N2)
1049: {
1050: PetscErrorCode ierr;
1051: DM_Stag * const stag = (DM_Stag*)dm->data;
1052: 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 (N0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in X direction must be positive");
1059: if (dim > 1 && N1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Y direction must be positive");
1060: if (dim > 2 && N2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Z direction must be positive");
1061: if (N0) stag->N[0] = N0;
1062: if (N1) stag->N[1] = N1;
1063: if (N2) stag->N[2] = N2;
1064: return(0);
1065: }
1067: /*@C
1068: DMStagSetOwnershipRanges - set elements per rank in each direction
1070: Logically Collective; lx, ly, and lz must contain common values
1072: Input Parameters:
1073: + dm - the DMStag object
1074: - lx,ly,lz - element counts for each rank in each direction
1076: Notes:
1077: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
1079: Level: developer
1081: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagGetOwnershipRanges(), DMDASetOwnershipRanges()
1082: @*/
1083: PetscErrorCode DMStagSetOwnershipRanges(DM dm,PetscInt const *lx,PetscInt const *ly,PetscInt const *lz)
1084: {
1085: PetscErrorCode ierr;
1086: DM_Stag * const stag = (DM_Stag*)dm->data;
1087: const PetscInt *lin[3];
1088: PetscInt d;
1092: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1093: lin[0] = lx; lin[1] = ly; lin[2] = lz;
1094: for (d=0; d<3; ++d) {
1095: if (lin[d]) {
1096: if (stag->nRanks[d] < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of ranks");
1097: if (!stag->l[d]) {
1098: PetscMalloc1(stag->nRanks[d], &stag->l[d]);
1099: }
1100: PetscArraycpy(stag->l[d], lin[d], stag->nRanks[d]);
1101: }
1102: }
1103: return(0);
1104: }
1106: /*@C
1107: DMStagSetUniformCoordinates - set DMStag coordinates to be a uniform grid
1109: Collective
1111: Input Parameters:
1112: + dm - the DMStag object
1113: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1115: Notes:
1116: DMStag supports 2 different types of coordinate DM: DMSTAG and DMPRODUCT.
1117: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1119: Level: advanced
1121: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType(), DMGetCoordinateDM(), DMGetCoordinates(), DMDASetUniformCoordinates()
1122: @*/
1123: PetscErrorCode DMStagSetUniformCoordinates(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1124: {
1125: PetscErrorCode ierr;
1126: DM_Stag * const stag = (DM_Stag*)dm->data;
1127: PetscBool flg_stag,flg_product;
1131: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1132: if (!stag->coordinateDMType) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"You must first call DMStagSetCoordinateDMType()");
1133: PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg_stag);
1134: PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg_product);
1135: if (flg_stag) {
1136: DMStagSetUniformCoordinatesExplicit(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1137: } else if (flg_product) {
1138: DMStagSetUniformCoordinatesProduct(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1139: } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported DM Type %s",stag->coordinateDMType);
1140: return(0);
1141: }
1143: /*@C
1144: DMStagSetUniformCoordinatesExplicit - set DMStag coordinates to be a uniform grid, storing all values
1146: Collective
1148: Input Parameters:
1149: + dm - the DMStag object
1150: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1152: Notes:
1153: DMStag supports 2 different types of coordinate DM: either another DMStag, or a DMProduct.
1154: If the grid is orthogonal, using DMProduct should be more efficient.
1155: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1157: Level: beginner
1159: .seealso: DMSTAG, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType()
1160: @*/
1161: PetscErrorCode DMStagSetUniformCoordinatesExplicit(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1162: {
1163: PetscErrorCode ierr;
1164: DM_Stag * const stag = (DM_Stag*)dm->data;
1165: PetscInt dim;
1166: PetscBool flg;
1170: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1171: PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg);
1172: if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1173: DMStagSetCoordinateDMType(dm,DMSTAG);
1174: DMGetDimension(dm,&dim);
1175: switch (dim) {
1176: case 1: DMStagSetUniformCoordinatesExplicit_1d(dm,xmin,xmax); break;
1177: case 2: DMStagSetUniformCoordinatesExplicit_2d(dm,xmin,xmax,ymin,ymax); break;
1178: case 3: DMStagSetUniformCoordinatesExplicit_3d(dm,xmin,xmax,ymin,ymax,zmin,zmax); break;
1179: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1180: }
1181: return(0);
1182: }
1184: /*@C
1185: DMStagSetUniformCoordinatesProduct - create uniform coordinates, as a product of 1D arrays
1187: 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
1188: coordinates.
1190: Collective
1192: Input Parameters:
1193: + dm - the DMStag object
1194: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1196: Notes:
1197: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1199: Level: intermediate
1201: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetCoordinateDMType()
1202: @*/
1203: PetscErrorCode DMStagSetUniformCoordinatesProduct(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1204: {
1205: PetscErrorCode ierr;
1206: DM_Stag * const stag = (DM_Stag*)dm->data;
1207: DM dmc;
1208: PetscInt dim,d,dof0,dof1;
1209: PetscBool flg;
1213: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1214: PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg);
1215: if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1216: DMStagSetCoordinateDMType(dm,DMPRODUCT);
1217: DMGetCoordinateDM(dm,&dmc);
1218: DMGetDimension(dm,&dim);
1220: /* Create 1D sub-DMs, living on subcommunicators */
1222: dof0 = 0;
1223: for (d=0; d<dim; ++d) { /* Include point dof in the sub-DMs if any non-element strata are active in the original DMStag */
1224: if (stag->dof[d]) {
1225: dof0 = 1;
1226: break;
1227: }
1228: }
1229: dof1 = stag->dof[dim] ? 1 : 0; /* Include element dof in the sub-DMs if the elements are active in the original DMStag */
1231: for (d=0; d<dim; ++d) {
1232: DM subdm;
1233: MPI_Comm subcomm;
1234: PetscMPIInt color;
1235: const PetscMPIInt key = 0; /* let existing rank break ties */
1237: /* Choose colors based on position in the plane orthogonal to this dim, and split */
1238: switch (d) {
1239: case 0: color = (dim > 1 ? stag->rank[1] : 0) + (dim > 2 ? stag->nRanks[1]*stag->rank[2] : 0); break;
1240: case 1: color = stag->rank[0] + (dim > 2 ? stag->nRanks[0]*stag->rank[2] : 0); break;
1241: case 2: color = stag->rank[0] + stag->nRanks[0]*stag->rank[1] ; break;
1242: default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1243: }
1244: MPI_Comm_split(PetscObjectComm((PetscObject)dm),color,key,&subcomm);
1246: /* Create sub-DMs living on these new communicators (which are destroyed by DMProduct) */
1247: DMStagCreate1d(subcomm,stag->boundaryType[d],stag->N[d],dof0,dof1,stag->stencilType,stag->stencilWidth,stag->l[d],&subdm);
1248: DMSetUp(subdm);
1249: switch (d) {
1250: case 0:
1251: DMStagSetUniformCoordinatesExplicit(subdm,xmin,xmax,0,0,0,0);
1252: break;
1253: case 1:
1254: DMStagSetUniformCoordinatesExplicit(subdm,ymin,ymax,0,0,0,0);
1255: break;
1256: case 2:
1257: DMStagSetUniformCoordinatesExplicit(subdm,zmin,zmax,0,0,0,0);
1258: break;
1259: default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1260: }
1261: DMProductSetDM(dmc,d,subdm);
1262: DMProductSetDimensionIndex(dmc,d,0);
1263: DMDestroy(&subdm);
1264: MPI_Comm_free(&subcomm);
1265: }
1266: return(0);
1267: }
1269: /*@C
1270: DMStagVecGetArrayDOF - get access to raw local array
1272: Logically Collective
1274: Input Parameters:
1275: + dm - the DMStag object
1276: - vec - the Vec object
1278: Output Parameters:
1279: . array - the array
1281: Notes:
1282: Indexing is array[k][j][i][idx].
1283: Obtain idx with DMStagGetLocationSlot().
1285: Level: beginner
1287: .seealso: DMSTAG, DMStagVecGetArrayDOFRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArrayDOF()
1288: @*/
1289: PetscErrorCode DMStagVecGetArrayDOF(DM dm,Vec vec,void *array)
1290: {
1291: PetscErrorCode ierr;
1292: DM_Stag * const stag = (DM_Stag*)dm->data;
1293: PetscInt dim;
1294: PetscInt nLocal;
1299: DMGetDimension(dm,&dim);
1300: VecGetLocalSize(vec,&nLocal);
1301: 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);
1302: switch (dim) {
1303: case 1:
1304: VecGetArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1305: break;
1306: case 2:
1307: VecGetArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1308: break;
1309: case 3:
1310: 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);
1311: break;
1312: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1313: }
1314: return(0);
1315: }
1317: /*@C
1318: DMStagVecGetArrayDOFRead - get read-only access to raw local array
1320: Logically Collective
1322: Input Parameters:
1323: + dm - the DMStag object
1324: - vec - the Vec object
1326: Output Parameters:
1327: . array - read-only the array
1329: Notes:
1330: Indexing is array[k][j][i][idx].
1331: Obtain idx with DMStagGetLocationSlot()
1333: Level: beginner
1335: .seealso: DMSTAG, DMStagVecGetArrayDOFRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArrayDOFRead()
1336: @*/
1337: PetscErrorCode DMStagVecGetArrayDOFRead(DM dm,Vec vec,void *array)
1338: {
1339: PetscErrorCode ierr;
1340: DM_Stag * const stag = (DM_Stag*)dm->data;
1341: PetscInt dim;
1342: PetscInt nLocal;
1347: DMGetDimension(dm,&dim);
1348: VecGetLocalSize(vec,&nLocal);
1349: 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);
1350: switch (dim) {
1351: case 1:
1352: VecGetArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1353: break;
1354: case 2:
1355: VecGetArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1356: break;
1357: case 3:
1358: 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);
1359: break;
1360: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1361: }
1362: return(0);
1363: }
1365: /*@C
1366: DMStagVecRestoreArrayDOF - restore read-only access to a raw array
1368: Logically Collective
1370: Input Parameters:
1371: + dm - the DMStag object
1372: - vec - the Vec object
1374: Output Parameters:
1375: . array - the array
1377: Level: beginner
1379: .seealso: DMSTAG, DMStagVecGetArrayDOF(), DMDAVecRestoreArrayDOFRead()
1380: @*/
1381: PetscErrorCode DMStagVecRestoreArrayDOF(DM dm,Vec vec,void *array)
1382: {
1383: PetscErrorCode ierr;
1384: DM_Stag * const stag = (DM_Stag*)dm->data;
1385: PetscInt dim;
1386: PetscInt nLocal;
1391: DMGetDimension(dm,&dim);
1392: VecGetLocalSize(vec,&nLocal);
1393: 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);
1394: switch (dim) {
1395: case 1:
1396: VecRestoreArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1397: break;
1398: case 2:
1399: VecRestoreArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1400: break;
1401: case 3:
1402: 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);
1403: break;
1404: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1405: }
1406: return(0);
1407: }
1409: /*@C
1410: DMStagVecRestoreArrayDOFRead - restore read-only access to a raw array
1412: Logically Collective
1414: Input Parameters:
1415: + dm - the DMStag object
1416: - vec - the Vec object
1418: Output Parameters:
1419: . array - the read-only array
1421: Level: beginner
1423: .seealso: DMSTAG, DMStagVecGetArrayDOFRead(), DMDAVecRestoreArrayDOFRead()
1424: @*/
1425: PetscErrorCode DMStagVecRestoreArrayDOFRead(DM dm,Vec vec,void *array)
1426: {
1427: PetscErrorCode ierr;
1428: DM_Stag * const stag = (DM_Stag*)dm->data;
1429: PetscInt dim;
1430: PetscInt nLocal;
1435: DMGetDimension(dm,&dim);
1436: VecGetLocalSize(vec,&nLocal);
1437: 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);
1438: switch (dim) {
1439: case 1:
1440: VecRestoreArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1441: break;
1442: case 2:
1443: VecRestoreArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1444: break;
1445: case 3:
1446: 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);
1447: break;
1448: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1449: }
1450: return(0);
1451: }