Actual source code: stagutils.c
petsc-3.11.4 2019-09-28
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: /*@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: arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
64: DMGetDimension(dm,&dim);
65: DMGetCoordinateDM(dm,&dmCoord);
66: if (!dmCoord) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM does not have a coordinate DM");
67: {
68: PetscBool isProduct;
69: DMType dmType;
70: DMGetType(dmCoord,&dmType);
71: PetscStrcmp(DMPRODUCT,dmType,&isProduct);
72: if (!isProduct) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is not of type DMPRODUCT");
73: }
74: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
75: for (d=0; d<dim; ++d) {
76: DM subDM;
77: DMType dmType;
78: PetscBool isStag;
79: PetscInt dof[DMSTAG_MAX_STRATA],subDim;
80: Vec coord1d;
81: DMProductGetDM(dmCoord,d,&subDM);
82: if (!subDM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is missing sub DM %D",d);
83: DMGetDimension(subDM,&subDim);
84: if (subDim != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of dimension 1");
85: DMGetType(subDM,&dmType);
86: PetscStrcmp(DMSTAG,dmType,&isStag);
87: if (!isStag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of type DMSTAG");
88: DMStagGetDOF(subDM,&dof[0],&dof[1],&dof[2],&dof[3]);
89: if (d == 0) {
90: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
91: } else {
92: for (s=0; s<DMSTAG_MAX_STRATA; ++s) {
93: if (dofCheck[s] != dof[s]) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DMs have different dofs");
94: }
95: }
96: DMGetCoordinatesLocal(subDM,&coord1d);
97: DMStagVecGetArrayDOFRead(subDM,coord1d,arr[d]);
98: }
99: return(0);
100: }
102: /*@C
103: DMStagGet1dCoordinateLocationSlot - get slot for use with local 1D coordinate arrays
105: High-level helper function to get slot ids for 1D coordinate DMs.
106: For use with DMStagGetIDCoordinateArraysDOFRead() and related functions.
108: Not Collective
110: Input Parameters:
111: + dm - the DMStag object
112: - loc - the grid location
114: Output Parameter:
115: . slot - the index to use in local arrays
117: Notes:
118: Checks that the coordinates are actually set up so that using the
119: slots from the first 1d coordinate sub-DM is valid for all the 1D coordinate sub-DMs.
121: Level: intermediate
123: .seealso: DMSTAG, DMPRODUCT, DMStagGet1dCoordinateArraysDOFRead(), DMStagSetUniformCoordinates()
124: @*/
125: PETSC_EXTERN PetscErrorCode DMStagGet1dCoordinateLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt *slot)
126: {
128: DM dmCoord;
129: PetscInt dim,dofCheck[DMSTAG_MAX_STRATA],s,d;
133: DMGetDimension(dm,&dim);
134: DMGetCoordinateDM(dm,&dmCoord);
135: if (!dmCoord) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM does not have a coordinate DM");
136: {
137: PetscBool isProduct;
138: DMType dmType;
139: DMGetType(dmCoord,&dmType);
140: PetscStrcmp(DMPRODUCT,dmType,&isProduct);
141: if (!isProduct) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is not of type DMPRODUCT");
142: }
143: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
144: for (d=0; d<dim; ++d) {
145: DM subDM;
146: DMType dmType;
147: PetscBool isStag;
148: PetscInt dof[DMSTAG_MAX_STRATA],subDim;
149: DMProductGetDM(dmCoord,d,&subDM);
150: if (!subDM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is missing sub DM %D",d);
151: DMGetDimension(subDM,&subDim);
152: if (subDim != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of dimension 1");
153: DMGetType(subDM,&dmType);
154: PetscStrcmp(DMSTAG,dmType,&isStag);
155: if (!isStag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of type DMSTAG");
156: DMStagGetDOF(subDM,&dof[0],&dof[1],&dof[2],&dof[3]);
157: if (d == 0) {
158: const PetscInt component = 0;
159: for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
160: DMStagGetLocationSlot(subDM,loc,component,slot);
161: } else {
162: for (s=0; s<DMSTAG_MAX_STRATA; ++s) {
163: if (dofCheck[s] != dof[s]) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DMs have different dofs");
164: }
165: }
166: }
167: return(0);
168: }
170: /*@C
171: DMStagGetCorners - return global element indices of the local region (excluding ghost points)
173: Not Collective
175: Input Parameter:
176: . dm - the DMStag object
178: Output Parameters:
179: + x,y,z - starting element indices in each direction
180: . m,n,p - element widths in each direction
181: - nExtrax,nExtray,nExtraz - number of extra partial elements in each direction. The number is 1 on the right, top, and front boundaries of the grid, otherwise 0.
183: Notes:
184: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
186: Level: beginner
188: .seealso: DMSTAG, DMStagGetGhostCorners()
189: @*/
190: PetscErrorCode DMStagGetCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *nExtrax,PetscInt *nExtray,PetscInt *nExtraz)
191: {
192: const DM_Stag * const stag = (DM_Stag*)dm->data;
196: if (x) *x = stag->start[0];
197: if (y) *y = stag->start[1];
198: if (z) *z = stag->start[2];
199: if (m) *m = stag->n[0];
200: if (n) *n = stag->n[1];
201: if (p) *p = stag->n[2];
202: if (nExtrax) *nExtrax = stag->lastRank[0] ? 1 : 0;
203: if (nExtray) *nExtray = stag->lastRank[1] ? 1 : 0;
204: if (nExtraz) *nExtraz = stag->lastRank[2] ? 1 : 0;
205: return(0);
206: }
208: /*@C
209: DMStagGetDOF - get number of DOF associated with each stratum of the grid
211: Not Collective
213: Input Parameter:
214: . dm - the DMStag object
216: Output Parameters:
217: + dof0 - the number of points per 0-cell (vertex/node)
218: . dof1 - the number of points per 1-cell (elementin 1D, edge in 2D and 3D)
219: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
220: - dof3 - the number of points per 3-cell (elementin 3D)
222: Level: beginner
224: .seealso: DMSTAG, DMStagGetCorners(), DMStagGetGhostCorners(), DMStagGetGlobalSizes(), DMStagGetStencilWidth(), DMStagGetBoundaryTypes()
225: @*/
226: PetscErrorCode DMStagGetDOF(DM dm,PetscInt *dof0,PetscInt *dof1,PetscInt *dof2,PetscInt *dof3)
227: {
228: const DM_Stag * const stag = (DM_Stag*)dm->data;
232: if (dof0) *dof0 = stag->dof[0];
233: if (dof1) *dof1 = stag->dof[1];
234: if (dof2) *dof2 = stag->dof[2];
235: if (dof3) *dof3 = stag->dof[3];
236: return(0);
237: }
239: /*@C
240: DMStagGetGhostCorners - return global element indices of the local region, including ghost points
242: Not Collective
244: Input Argument:
245: + dm - the DMStag object
247: Output Arguments:
248: + x,y,z - starting element indices in each direction
249: - m,n,p - element widths in each direction
251: Notes:
252: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
254: Level: beginner
256: .seealso: DMSTAG, DMStagGetCorners()
257: @*/
258: PetscErrorCode DMStagGetGhostCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
259: {
260: const DM_Stag * const stag = (DM_Stag*)dm->data;
264: if (x) *x = stag->startGhost[0];
265: if (y) *y = stag->startGhost[1];
266: if (z) *z = stag->startGhost[2];
267: if (m) *m = stag->nGhost[0];
268: if (n) *n = stag->nGhost[1];
269: if (p) *p = stag->nGhost[2];
270: return(0);
271: }
273: /*@C
274: DMStagGetGlobalSizes - get global element counts
276: Not Collective
278: Input Parameter:
279: . dm - the DMStag object
281: Output Parameters:
282: . M,N,P - global element counts in each direction
284: Notes:
285: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
287: Level: beginner
289: .seealso: DMSTAG, DMStagGetLocalSizes()
290: @*/
291: PetscErrorCode DMStagGetGlobalSizes(DM dm,PetscInt* M,PetscInt* N,PetscInt* P)
292: {
293: const DM_Stag * const stag = (DM_Stag*)dm->data;
297: if (M) *M = stag->N[0];
298: if (N) *N = stag->N[1];
299: if (P) *P = stag->N[2];
300: return(0);
301: }
303: /*@C
304: DMStagGetIsFirstRank - get boolean value for whether this rank is first in each direction in the rank grid
306: Not Collective
308: Input Parameter:
309: . dm - the DMStag object
311: Output Parameters:
312: . isFirstRank0,isFirstRank1,isFirstRank2 - whether this rank is first in each direction
314: Level: intermediate
316: Notes:
317: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
319: .seealso: DMSTAG, DMStagGetIsLastRank()
320: @*/
321: PetscErrorCode DMStagGetIsFirstRank(DM dm,PetscBool *isFirstRank0,PetscBool *isFirstRank1,PetscBool *isFirstRank2)
322: {
323: const DM_Stag * const stag = (DM_Stag*)dm->data;
327: if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
328: if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
329: if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
330: return(0);
331: }
333: /*@C
334: DMStagGetIsLastRank - get boolean value for whether this rank is last in each direction in the rank grid
336: Not Collective
338: Input Parameter:
339: . dm - the DMStag object
341: Output Parameters:
342: . isLastRank0,isLastRank1,isLastRank2 - whether this rank is last in each direction
344: Level: intermediate
346: Notes:
347: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
348: Level: intermediate
350: .seealso: DMSTAG, DMStagGetIsFirstRank()
351: @*/
352: PetscErrorCode DMStagGetIsLastRank(DM dm,PetscBool *isLastRank0,PetscBool *isLastRank1,PetscBool *isLastRank2)
353: {
354: const DM_Stag * const stag = (DM_Stag*)dm->data;
358: if (isLastRank0) *isLastRank0 = stag->lastRank[0];
359: if (isLastRank1) *isLastRank1 = stag->lastRank[1];
360: if (isLastRank2) *isLastRank2 = stag->lastRank[2];
361: return(0);
362: }
364: /*@C
365: DMStagGetLocalSizes - get local elementwise sizes
367: Not Collective
369: Input Parameter:
370: . dm - the DMStag object
372: Output Parameters:
373: . m,n,p - local element counts (excluding ghosts) in each direction
375: Notes:
376: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
377: Level: intermediate
379: Level: beginner
381: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetDOF(), DMStagGetNumRanks()
382: @*/
383: PetscErrorCode DMStagGetLocalSizes(DM dm,PetscInt* m,PetscInt* n,PetscInt* p)
384: {
385: const DM_Stag * const stag = (DM_Stag*)dm->data;
389: if (m) *m = stag->n[0];
390: if (n) *n = stag->n[1];
391: if (p) *p = stag->n[2];
392: return(0);
393: }
395: /*@C
396: DMStagGetNumRanks - get number of ranks in each direction in the global grid decomposition
398: Not Collective
400: Input Parameter:
401: . dm - the DMStag object
403: Output Parameters:
404: . nRanks0,nRanks1,nRanks2 - number of ranks in each direction in the grid decomposition
406: Notes:
407: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
408: Level: intermediate
410: Level: beginner
412: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetLocalSize(), DMStagSetNumRank()
413: @*/
414: PetscErrorCode DMStagGetNumRanks(DM dm,PetscInt *nRanks0,PetscInt *nRanks1,PetscInt *nRanks2)
415: {
416: const DM_Stag * const stag = (DM_Stag*)dm->data;
420: if (nRanks0) *nRanks0 = stag->nRanks[0];
421: if (nRanks1) *nRanks1 = stag->nRanks[1];
422: if (nRanks2) *nRanks2 = stag->nRanks[2];
423: return(0);
424: }
426: /*@C
427: DMStagGetEntriesPerElement - get number of entries per element in the local representation
429: Not Collective
431: Input Parameter:
432: . dm - the DMStag object
434: Output Parameters:
435: . entriesPerElement - number of entries associated with each element in the local representation
437: Notes:
438: This is the natural block size for most local operations. In 1D it is equal to dof0 + dof1,
439: in 2D it is equal to dof0 + 2*dof1 + dof2, and in 3D it is equal to dof0 + 3*dof1 + 3*dof2 + dof3
441: Level: developer
443: .seealso: DMSTAG, DMStagGetDOF()
444: @*/
445: PetscErrorCode DMStagGetEntriesPerElement(DM dm,PetscInt *entriesPerElement)
446: {
447: const DM_Stag * const stag = (DM_Stag*)dm->data;
451: if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
452: return(0);
453: }
455: /*@C
456: DMStagGetStencilWidth - get elementwise stencil width
458: Not Collective
460: Input Parameter:
461: . dm - the DMStag object
463: Output Parameters:
464: . stencilWidth - stencil/halo/ghost width in elements
466: Level: beginner
468: .seealso: DMSTAG
469: @*/
470: PetscErrorCode DMStagGetStencilWidth(DM dm,PetscInt *stencilWidth)
471: {
472: const DM_Stag * const stag = (DM_Stag*)dm->data;
476: if (stencilWidth) *stencilWidth = stag->stencilWidth;
477: return(0);
478: }
480: /*@C
481: DMStagCreateCompatibleDMStag - create a compatible DMStag with different dof/stratum
483: Collective
485: Input Parameters
486: + dm - the DMStag object
487: - dof0,dof1,dof2,dof3 - number of dof on each stratum in the new DMStag
489: Output Parameters:
490: . newdm - the new, compatible DMStag
492: Notes:
493: Dof supplied for strata too big for the dimension are ignored; these may be set to 0.
494: In contrast to DMDACreateCompatibleDMDA(), coordinates are not reused.
496: Level: intermediate
498: .seealso: DMSTAG, DMDACreateCompatibleDMDA(), DMGetCompatibility(), DMStagMigrateVec()
499: @*/
500: PetscErrorCode DMStagCreateCompatibleDMStag(DM dm,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DM *newdm)
501: {
502: PetscErrorCode ierr;
503: const DM_Stag * const stag = (DM_Stag*)dm->data;
504: PetscInt dim;
508: DMGetDimension(dm,&dim);
509: switch (dim) {
510: case 1:
511: DMStagCreate1d(PetscObjectComm((PetscObject)dm),stag->boundaryType[0],stag->N[0],dof0,dof1,stag->stencilType,stag->stencilWidth,NULL,newdm);
512: break;
513: case 2:
514: 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);
515: break;
516: case 3:
517: 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);
518: break;
519: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
520: }
521: DMSetUp(*newdm);
522: return(0);
523: }
525: /*@C
526: DMStagGetLocationSlot - get index to use in accessing raw local arrays
528: Not Collective
530: Input Parameters:
531: + dm - the DMStag object
532: . loc - location relative to an element
533: - c - component
535: Output Parameter:
536: . slot - index to use
538: Notes:
539: Provides an appropriate index to use with DMStagVecGetArrayDOF() and friends.
540: This is required so that the user doesn't need to know about the ordering of
541: dof associated with each local element.
543: Level: beginner
545: .seealso: DMSTAG, DMStagVecGetArrayDOF(), DMStagVecGetArrayDOFRead(), DMStagGetDOF(), DStagMGetEntriesPerElement()
546: @*/
547: PetscErrorCode DMStagGetLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt c,PetscInt *slot)
548: {
549: DM_Stag * const stag = (DM_Stag*)dm->data;
553: #if defined(PETSC_USE_DEBUG)
554: {
556: PetscInt dof;
557: DMStagGetLocationDOF(dm,loc,&dof);
558: if (dof < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Location %s has no dof attached",DMStagStencilLocations[loc]);
559: 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);
560: }
561: #endif
562: *slot = stag->locationOffsets[loc] + c;
563: return(0);
564: }
566: /*@C
567: DMStagMigrateVec - transfer a vector associated with a DMStag to a vector associated with a compatible DMStag
569: Collective
571: Input Parameters:
572: + dm - the source DMStag object
573: . vec - the destination vector, compatible with dm
574: . dmTo - the compatible destination DMStag object
575: - vecTo - the destination vector, compatible with dmTo
577: Notes:
578: Extra dof are ignored, and unfilled dof are zeroed.
579: Currently only implemented to migrate global vectors to global vectors.
581: Level: advanced
583: .seealso: DMSTAG, DMStagCreateCompatibleDMStag(), DMGetCompatibility(), DMStagVecSplitToDMDA()
584: @*/
585: PetscErrorCode DMStagMigrateVec(DM dm,Vec vec,DM dmTo,Vec vecTo)
586: {
587: PetscErrorCode ierr;
588: DM_Stag * const stag = (DM_Stag*)dm->data;
589: DM_Stag * const stagTo = (DM_Stag*)dmTo->data;
590: PetscInt nLocalTo,nLocal,dim,i,j,k;
591: PetscInt start[DMSTAG_MAX_DIM],startGhost[DMSTAG_MAX_DIM],n[DMSTAG_MAX_DIM],nExtra[DMSTAG_MAX_DIM],offset[DMSTAG_MAX_DIM];
592: Vec vecToLocal,vecLocal;
593: PetscBool compatible,compatibleSet;
594: const PetscScalar *arr;
595: PetscScalar *arrTo;
596: const PetscInt epe = stag->entriesPerElement;
597: const PetscInt epeTo = stagTo->entriesPerElement;
604: DMGetCompatibility(dm,dmTo,&compatible,&compatibleSet);
605: if (!compatibleSet || !compatible) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"DMStag objects must be shown to be compatible");
606: DMGetDimension(dm,&dim);
607: VecGetLocalSize(vec,&nLocal);
608: VecGetLocalSize(vecTo,&nLocalTo);
609: if (nLocal != stag->entries|| nLocalTo !=stagTo->entries) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Vector migration only implemented for global vector to global vector.");
610: DMStagGetCorners(dm,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
611: DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],NULL,NULL,NULL);
612: for (i=0; i<DMSTAG_MAX_DIM; ++i) offset[i] = start[i]-startGhost[i];
614: /* Proceed by transferring to a local vector, copying, and transferring back to a global vector */
615: DMGetLocalVector(dm,&vecLocal);
616: DMGetLocalVector(dmTo,&vecToLocal);
617: DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal);
618: DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal);
619: VecGetArrayRead(vecLocal,&arr);
620: VecGetArray(vecToLocal,&arrTo);
621: /* Note that some superfluous copying of entries on partial dummy elements is done */
622: if (dim == 1) {
623: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
624: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
625: PetscInt si;
626: for (si=0; si<2; ++si) {
627: b += stag->dof[si];
628: bTo += stagTo->dof[si];
629: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[i*epeTo + dTo] = arr[i*epe + d];
630: for (; dTo < bTo ; ++dTo) arrTo[i*epeTo + dTo] = 0.0;
631: d = b;
632: }
633: }
634: } else if (dim == 2) {
635: const PetscInt epr = stag->nGhost[0] * epe;
636: const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
637: for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
638: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
639: const PetscInt base = j*epr + i*epe;
640: const PetscInt baseTo = j*eprTo + i*epeTo;
641: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
642: const PetscInt s[4] = {0,1,1,2}; /* Dimensions of points, in order */
643: PetscInt si;
644: for (si=0; si<4; ++si) {
645: b += stag->dof[s[si]];
646: bTo += stagTo->dof[s[si]];
647: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
648: for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
649: d = b;
650: }
651: }
652: }
653: } else if (dim == 3) {
654: const PetscInt epr = stag->nGhost[0] * epe;
655: const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
656: const PetscInt epl = stag->nGhost[1] * epr;
657: const PetscInt eplTo = stagTo->nGhost[1] * eprTo;
658: for (k=offset[2]; k<offset[2] + n[2] + nExtra[2]; ++k) {
659: for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
660: for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
661: PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
662: const PetscInt base = k*epl + j*epr + i*epe;
663: const PetscInt baseTo = k*eplTo + j*eprTo + i*epeTo;
664: const PetscInt s[8] = {0,1,1,2,1,2,2,3}; /* dimensions of points, in order */
665: PetscInt is;
666: for (is=0; is<8; ++is) {
667: b += stag->dof[s[is]];
668: bTo += stagTo->dof[s[is]];
669: for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
670: for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
671: d = b;
672: }
673: }
674: }
675: }
676: } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
677: VecRestoreArrayRead(vecLocal,&arr);
678: VecRestoreArray(vecToLocal,&arrTo);
679: DMRestoreLocalVector(dm,&vecLocal);
680: DMLocalToGlobalBegin(dmTo,vecToLocal,INSERT_VALUES,vecTo);
681: DMLocalToGlobalEnd(dmTo,vecToLocal,INSERT_VALUES,vecTo);
682: DMRestoreLocalVector(dmTo,&vecToLocal);
683: return(0);
684: }
686: /*@C
687: DMStagRestore1dCoordinateArraysDOFRead - restore local array access
689: Logically Collective
691: Input Parameter:
692: . dm - the DMStag object
694: Output Parameters:
695: . arrX,arrY,arrX - local 1D coordinate arrays
697: Level: intermediate
699: .seealso: DMSTAG, DMStagGet1dCoordinateArraysDOFRead()
700: @*/
701: PetscErrorCode DMStagRestore1dCoordinateArraysDOFRead(DM dm,void *arrX,void *arrY,void *arrZ)
702: {
703: PetscErrorCode ierr;
704: PetscInt dim,d;
705: void* arr[DMSTAG_MAX_DIM];
706: DM dmCoord;
710: DMGetDimension(dm,&dim);
711: arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
712: DMGetCoordinateDM(dm,&dmCoord);
713: for (d=0; d<dim; ++d) {
714: DM subDM;
715: Vec coord1d;
716: DMProductGetDM(dmCoord,d,&subDM);
717: DMGetCoordinatesLocal(subDM,&coord1d);
718: DMStagVecRestoreArrayDOFRead(subDM,coord1d,arr[d]);
719: }
720: return(0);
721: }
723: /*@C
724: DMStagSetBoundaryTypes - set DMStag boundary types
726: Logically Collective
728: Input Parameters:
729: + dm - the DMStag object
730: - boundaryType0,boundaryType1,boundaryType2 - boundary types in each direction
732: Notes:
733: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
735: Level: advanced
737: .seealso: DMSTAG, DMBoundaryType, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d()
738: @*/
739: PetscErrorCode DMStagSetBoundaryTypes(DM dm,DMBoundaryType boundaryType0,DMBoundaryType boundaryType1,DMBoundaryType boundaryType2)
740: {
741: PetscErrorCode ierr;
742: DM_Stag * const stag = (DM_Stag*)dm->data;
743: PetscInt dim;
750: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
751: DMGetDimension(dm,&dim);
752: if (boundaryType0 ) stag->boundaryType[0] = boundaryType0;
753: if (boundaryType1 && dim > 1) stag->boundaryType[1] = boundaryType1;
754: if (boundaryType2 && dim > 2) stag->boundaryType[2] = boundaryType2;
755: return(0);
756: }
758: /*@C
759: DMStagSetCoordinateDMType - set DM type to store coordinates
761: Logically Collective
763: Input Parameters:
764: + dm - the DMStag object
765: - dmtype - DMtype for coordinates, either DMSTAG or DMPRODUCT
767: Level: advanced
769: .seealso: DMSTAG, DMPRODUCT, DMGetCoordinateDM(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMType
770: @*/
771: PetscErrorCode DMStagSetCoordinateDMType(DM dm,DMType dmtype)
772: {
773: DM_Stag * const stag = (DM_Stag*)dm->data;
777: stag->coordinateDMType = dmtype;
778: return(0);
779: }
781: /*@C
782: DMStagSetDOF - set dof/stratum
784: Logically Collective
786: Input Parameters:
787: + dm - the DMStag object
788: - dof0,dof1,dof2,dof3 - dof per stratum
790: Notes:
791: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
793: Level: advanced
795: .seealso: DMSTAG
796: @*/
797: PetscErrorCode DMStagSetDOF(DM dm,PetscInt dof0, PetscInt dof1,PetscInt dof2,PetscInt dof3)
798: {
799: PetscErrorCode ierr;
800: DM_Stag * const stag = (DM_Stag*)dm->data;
801: PetscInt dim;
809: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
810: DMGetDimension(dm,&dim);
811: if (dof0 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof0 cannot be negative");
812: if (dof1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof1 cannot be negative");
813: if (dim > 1 && dof2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof2 cannot be negative");
814: if (dim > 2 && dof3 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof3 cannot be negative");
815: stag->dof[0] = dof0;
816: stag->dof[1] = dof1;
817: if (dim > 1) stag->dof[2] = dof2;
818: if (dim > 2) stag->dof[3] = dof3;
819: return(0);
820: }
822: /*@C
823: DMStagSetNumRanks - set ranks in each direction in the global rank grid
825: Logically Collective
827: Input Parameters:
828: + dm - the DMStag object
829: - nRanks0,nRanks1,nRanks2 - number of ranks in each direction
831: Notes:
832: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
834: Level: developer
836: .seealso: DMSTAG
837: @*/
838: PetscErrorCode DMStagSetNumRanks(DM dm,PetscInt nRanks0,PetscInt nRanks1,PetscInt nRanks2)
839: {
840: DM_Stag * const stag = (DM_Stag*)dm->data;
847: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
848: if (nRanks0) stag->nRanks[0] = nRanks0;
849: if (nRanks1) stag->nRanks[1] = nRanks1;
850: if (nRanks2) stag->nRanks[2] = nRanks2;
851: return(0);
852: }
854: /*@C
855: DMStagSetGhostType - set elementwise ghost/halo stencil type
857: Logically Collective
859: Input Parameters:
860: + dm - the DMStag object
861: - stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE
863: Level: beginner
865: .seealso: DMSTAG
866: @*/
867: PetscErrorCode DMStagSetGhostType(DM dm,DMStagStencilType stencilType)
868: {
869: DM_Stag * const stag = (DM_Stag*)dm->data;
874: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
875: stag->stencilType = stencilType;
876: return(0);
877: }
879: /*@C
880: DMStagSetGlobalSizes -
882: Logically Collective
884: Input Parameters:
885: + dm - the DMStag object
886: - N0,N1,N2 - global elementwise sizes
888: Notes:
889: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
891: Level: advanced
893: .seealso: DMSTAG, DMStagGetGlobalSizes()
894: @*/
895: PetscErrorCode DMStagSetGlobalSizes(DM dm,PetscInt N0,PetscInt N1,PetscInt N2)
896: {
897: DM_Stag * const stag = (DM_Stag*)dm->data;
901: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
902: if (N0) stag->N[0] = N0;
903: if (N1) stag->N[1] = N1;
904: if (N2) stag->N[2] = N2;
905: return(0);
906: }
908: /*@C
909: DMStagSetOwnershipRanges - set elements per rank in each direction
911: Logically Collective
913: Input Parameters:
914: + dm - the DMStag object
915: - lx,ly,lz - element counts for each rank in each direction
917: Notes:
918: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
920: Level: developer
922: .seealso: DMSTAG, DMStagSetGlobalSizes
923: @*/
924: PetscErrorCode DMStagSetOwnershipRanges(DM dm,PetscInt const *lx,PetscInt const *ly,PetscInt const *lz)
925: {
926: PetscErrorCode ierr;
927: DM_Stag * const stag = (DM_Stag*)dm->data;
928: const PetscInt *lin[3];
929: PetscInt d;
933: if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
934: lin[0] = lx; lin[1] = ly; lin[2] = lz;
935: for (d=0; d<3; ++d) {
936: if (lin[d]) {
937: if (stag->nRanks[d] < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of ranks");
938: if (!stag->l[d]) {
939: PetscMalloc1(stag->nRanks[d], &stag->l[d]);
940: }
941: PetscMemcpy(stag->l[d], lin[d], stag->nRanks[d]*sizeof(PetscInt));
942: }
943: }
944: return(0);
945: }
947: /*@C
948: DMStagSetUniformCoordinates - set DMStag coordinates to be a uniform grid
950: Collective
952: Input Parameters:
953: + dm - the DMStag object
954: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
956: Notes:
957: DMStag supports 2 different types of coordinate DM: DMSTAG and DMPRODUCT.
958: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
960: Level: advanced
962: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType(), DMGetCoordinateDM(), DMGetCoordinates()
963: @*/
964: PetscErrorCode DMStagSetUniformCoordinates(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
965: {
966: PetscErrorCode ierr;
967: DM_Stag * const stag = (DM_Stag*)dm->data;
968: PetscBool flg;
972: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
973: if (!stag->coordinateDMType) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"You must first call DMStagSetCoordinateDMType()");
974: PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg);
975: if (flg) {
976: DMStagSetUniformCoordinatesExplicit(dm,xmin,xmax,ymin,ymax,zmin,zmax);
977: } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported DM Type %s",stag->coordinateDMType);
978: return(0);
979: }
981: /*@C
982: DMStagSetUniformCoordinatesExplicit - set DMStag coordinates to be a uniform grid, storing all values
984: Collective
986: Input Parameters:
987: + dm - the DMStag object
988: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
990: Notes:
991: DMStag supports 2 different types of coordinate DM: either another DMStag, or a DMProduct.
992: If the grid is orthogonal, using DMProduct should be more efficient.
993: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
995: Level: beginner
997: .seealso: DMSTAG, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType()
998: @*/
999: PetscErrorCode DMStagSetUniformCoordinatesExplicit(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1000: {
1001: PetscErrorCode ierr;
1002: DM_Stag * const stag = (DM_Stag*)dm->data;
1003: PetscInt dim;
1004: PetscBool flg;
1008: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1009: PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg);
1010: if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1011: DMStagSetCoordinateDMType(dm,DMSTAG);
1012: DMGetDimension(dm,&dim);
1013: switch (dim) {
1014: case 1: DMStagSetUniformCoordinatesExplicit_1d(dm,xmin,xmax); break;
1015: case 2: DMStagSetUniformCoordinatesExplicit_2d(dm,xmin,xmax,ymin,ymax); break;
1016: case 3: DMStagSetUniformCoordinatesExplicit_3d(dm,xmin,xmax,ymin,ymax,zmin,zmax); break;
1017: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1018: }
1019: return(0);
1020: }
1022: /*@C
1023: DMStagSetUniformCoordinatesProduct - create uniform coordinates, as a product of 1D arrays
1025: 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
1026: coordinates.
1028: Collective
1030: Input Parameters:
1031: + dm - the DMStag object
1032: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1034: Notes:
1035: Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1037: Level: intermediate
1039: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetCoordinateDMType()
1040: @*/
1041: PetscErrorCode DMStagSetUniformCoordinatesProduct(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1042: {
1043: PetscErrorCode ierr;
1044: DM_Stag * const stag = (DM_Stag*)dm->data;
1045: DM dmc;
1046: PetscInt dim,d,dof0,dof1;
1047: PetscBool flg;
1051: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1052: PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg);
1053: if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1054: DMStagSetCoordinateDMType(dm,DMPRODUCT);
1055: DMGetCoordinateDM(dm,&dmc);
1056: DMGetDimension(dm,&dim);
1058: /* Create 1D sub-DMs, living on subcommunicators */
1060: dof0 = 0;
1061: for (d=0; d<dim; ++d) { /* Include point dof in the sub-DMs if any non-element strata are active in the original DMStag */
1062: if (stag->dof[d]) {
1063: dof0 = 1;
1064: break;
1065: }
1066: }
1067: dof1 = stag->dof[dim] ? 1 : 0; /* Include element dof in the sub-DMs if the elements are active in the original DMStag */
1069: for (d=0; d<dim; ++d) {
1070: DM subdm;
1071: MPI_Comm subcomm;
1072: PetscMPIInt color;
1073: const PetscMPIInt key = 0; /* let existing rank break ties */
1075: /* Choose colors based on position in the plane orthogonal to this dim, and split */
1076: switch (d) {
1077: case 0: color = (dim > 1 ? stag->rank[1] : 0) + (dim > 2 ? stag->nRanks[1]*stag->rank[2] : 0); break;
1078: case 1: color = stag->rank[0] + (dim > 2 ? stag->nRanks[0]*stag->rank[2] : 0); break;
1079: case 2: color = stag->rank[0] + stag->nRanks[0]*stag->rank[1] ; break;
1080: default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1081: }
1082: MPI_Comm_split(PetscObjectComm((PetscObject)dm),color,key,&subcomm);
1084: /* Create sub-DMs living on these new communicators (which are destroyed by DMProduct) */
1085: DMStagCreate1d(subcomm,stag->boundaryType[d],stag->N[d],dof0,dof1,stag->stencilType,stag->stencilWidth,stag->l[d],&subdm);
1086: DMSetUp(subdm);
1087: switch (d) {
1088: case 0:
1089: DMStagSetUniformCoordinatesExplicit(subdm,xmin,xmax,0,0,0,0);
1090: break;
1091: case 1:
1092: DMStagSetUniformCoordinatesExplicit(subdm,ymin,ymax,0,0,0,0);
1093: break;
1094: case 2:
1095: DMStagSetUniformCoordinatesExplicit(subdm,zmin,zmax,0,0,0,0);
1096: break;
1097: default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1098: }
1099: DMProductSetDM(dmc,d,subdm);
1100: DMProductSetDimensionIndex(dmc,d,0);
1101: DMDestroy(&subdm);
1102: MPI_Comm_free(&subcomm);
1103: }
1104: return(0);
1105: }
1107: /*@C
1108: DMStagVecGetArrayDOF - get access to raw local array
1110: Logically Collective
1112: Input Parameters:
1113: + dm - the DMStag object
1114: - vec - the Vec object
1116: Output Parameters:
1117: . array - the array
1119: Notes:
1120: Indexing is array[k][j][i][idx].
1121: Obtain idx with DMStagGetLocationSlot().
1123: Level: beginner
1125: .seealso: DMSTAG, DMStagVecGetArrayDOFRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector()
1126: @*/
1127: PetscErrorCode DMStagVecGetArrayDOF(DM dm,Vec vec,void *array)
1128: {
1129: PetscErrorCode ierr;
1130: DM_Stag * const stag = (DM_Stag*)dm->data;
1131: PetscInt dim;
1132: PetscInt nLocal;
1137: DMGetDimension(dm,&dim);
1138: VecGetLocalSize(vec,&nLocal);
1139: 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);
1140: switch (dim) {
1141: case 1:
1142: VecGetArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1143: break;
1144: case 2:
1145: VecGetArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1146: break;
1147: case 3:
1148: 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);
1149: break;
1150: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1151: }
1152: return(0);
1153: }
1155: /*@C
1156: DMStagVecGetArrayDOFRead - get read-only access to raw local array
1158: Logically Collective
1160: Input Parameters:
1161: + dm - the DMStag object
1162: - vec - the Vec object
1164: Output Parameters:
1165: . array - read-only the array
1167: Notes:
1168: Indexing is array[k][j][i][idx].
1169: Obtain idx with DMStagGetLocationSlot()
1171: Level: beginner
1173: .seealso: DMSTAG, DMStagVecGetArrayDOFRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector()
1174: @*/
1175: PetscErrorCode DMStagVecGetArrayDOFRead(DM dm,Vec vec,void *array)
1176: {
1177: PetscErrorCode ierr;
1178: DM_Stag * const stag = (DM_Stag*)dm->data;
1179: PetscInt dim;
1180: PetscInt nLocal;
1185: DMGetDimension(dm,&dim);
1186: VecGetLocalSize(vec,&nLocal);
1187: 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);
1188: switch (dim) {
1189: case 1:
1190: VecGetArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1191: break;
1192: case 2:
1193: VecGetArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1194: break;
1195: case 3:
1196: 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);
1197: break;
1198: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1199: }
1200: return(0);
1201: }
1203: /*@C
1204: DMStagVecRestoreArrayDOF - restore read-only access to a raw array
1206: Logically Collective
1208: Input Parameters:
1209: + dm - the DMStag object
1210: - vec - the Vec object
1212: Output Parameters:
1213: . array - the array
1215: Level: beginner
1217: .seealso: DMSTAG, DMStagVecGetArrayDOF()
1218: @*/
1219: PetscErrorCode DMStagVecRestoreArrayDOF(DM dm,Vec vec,void *array)
1220: {
1221: PetscErrorCode ierr;
1222: DM_Stag * const stag = (DM_Stag*)dm->data;
1223: PetscInt dim;
1224: PetscInt nLocal;
1229: DMGetDimension(dm,&dim);
1230: VecGetLocalSize(vec,&nLocal);
1231: 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);
1232: switch (dim) {
1233: case 1:
1234: VecRestoreArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1235: break;
1236: case 2:
1237: VecRestoreArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1238: break;
1239: case 3:
1240: 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);
1241: break;
1242: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1243: }
1244: return(0);
1245: }
1247: /*@C
1248: DMStagVecRestoreArrayDOFRead - restore read-only access to a raw array
1250: Logically Collective
1252: Input Parameters:
1253: + dm - the DMStag object
1254: - vec - the Vec object
1256: Output Parameters:
1257: . array - the read-only array
1259: Level: beginner
1261: .seealso: DMSTAG, DMStagVecGetArrayDOFRead()
1262: @*/
1263: PetscErrorCode DMStagVecRestoreArrayDOFRead(DM dm,Vec vec,void *array)
1264: {
1265: PetscErrorCode ierr;
1266: DM_Stag * const stag = (DM_Stag*)dm->data;
1267: PetscInt dim;
1268: PetscInt nLocal;
1273: DMGetDimension(dm,&dim);
1274: VecGetLocalSize(vec,&nLocal);
1275: 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);
1276: switch (dim) {
1277: case 1:
1278: VecRestoreArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1279: break;
1280: case 2:
1281: VecRestoreArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1282: break;
1283: case 3:
1284: 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);
1285: break;
1286: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1287: }
1288: return(0);
1289: }