Actual source code: stagutils.c

  1: /* Additional functions in the DMStag API, which are not part of the general DM API. */
  2: #include <petsc/private/dmstagimpl.h>
  3: #include <petscdmproduct.h>
  4: /*@C
  5:   DMStagGetBoundaryTypes - get boundary types

  7:   Not Collective

  9:   Input Parameter:
 10: . dm - the DMStag object

 12:   Output Parameters:
 13: . boundaryTypeX,boundaryTypeY,boundaryTypeZ - boundary types

 15:   Level: intermediate

 17: .seealso: DMSTAG
 18: @*/
 19: PetscErrorCode DMStagGetBoundaryTypes(DM dm,DMBoundaryType *boundaryTypeX,DMBoundaryType *boundaryTypeY,DMBoundaryType *boundaryTypeZ)
 20: {
 21:   PetscErrorCode        ierr;
 22:   const DM_Stag * const stag  = (DM_Stag*)dm->data;
 23:   PetscInt              dim;

 27:   DMGetDimension(dm,&dim);
 28:   if (boundaryTypeX) *boundaryTypeX = stag->boundaryType[0];
 29:   if (boundaryTypeY && dim > 1) *boundaryTypeY = stag->boundaryType[1];
 30:   if (boundaryTypeZ && dim > 2) *boundaryTypeZ = stag->boundaryType[2];
 31:   return(0);
 32: }

 34: static PetscErrorCode DMStagGetProductCoordinateArrays_Private(DM dm,void* arrX,void* arrY,void* arrZ,PetscBool read)
 35: {
 37:   PetscInt       dim,d,dofCheck[DMSTAG_MAX_STRATA],s;
 38:   DM             dmCoord;
 39:   void*          arr[DMSTAG_MAX_DIM];
 40:   PetscBool      checkDof;

 44:   DMGetDimension(dm,&dim);
 45:   if (dim > DMSTAG_MAX_DIM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for %D dimensions",dim);
 46:   arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
 47:   DMGetCoordinateDM(dm,&dmCoord);
 48:   if (!dmCoord) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM does not have a coordinate DM");
 49:   {
 50:     PetscBool isProduct;
 51:     DMType    dmType;
 52:     DMGetType(dmCoord,&dmType);
 53:     PetscStrcmp(DMPRODUCT,dmType,&isProduct);
 54:     if (!isProduct) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is not of type DMPRODUCT");
 55:   }
 56:   for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
 57:   checkDof = PETSC_FALSE;
 58:   for (d=0; d<dim; ++d) {
 59:     DM        subDM;
 60:     DMType    dmType;
 61:     PetscBool isStag;
 62:     PetscInt  dof[DMSTAG_MAX_STRATA],subDim;
 63:     Vec       coord1d_local;

 65:     /* Ignore unrequested arrays */
 66:     if (!arr[d]) continue;

 68:     DMProductGetDM(dmCoord,d,&subDM);
 69:     if (!subDM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is missing sub DM %D",d);
 70:     DMGetDimension(subDM,&subDim);
 71:     if (subDim != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of dimension 1");
 72:     DMGetType(subDM,&dmType);
 73:     PetscStrcmp(DMSTAG,dmType,&isStag);
 74:     if (!isStag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of type DMSTAG");
 75:     DMStagGetDOF(subDM,&dof[0],&dof[1],&dof[2],&dof[3]);
 76:     if (!checkDof) {
 77:       for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
 78:       checkDof = PETSC_TRUE;
 79:     } else {
 80:       for (s=0; s<DMSTAG_MAX_STRATA; ++s) {
 81:         if (dofCheck[s] != dof[s]) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DMs have different dofs");
 82:       }
 83:     }
 84:     DMGetCoordinatesLocal(subDM,&coord1d_local);
 85:     if (read) {
 86:       DMStagVecGetArrayRead(subDM,coord1d_local,arr[d]);
 87:     } else {
 88:       DMStagVecGetArray(subDM,coord1d_local,arr[d]);
 89:     }
 90:   }
 91:   return(0);
 92: }

 94: /*@C
 95:   DMStagGetProductCoordinateArrays - extract local product coordinate arrays, one per dimension

 97:   Logically Collective

 99:   A high-level helper function to quickly extract local coordinate arrays.

101:   Note that 2-dimensional arrays are returned. See
102:   DMStagVecGetArray(), which is called internally to produce these arrays
103:   representing coordinates on elements and vertices (element boundaries)
104:   for a 1-dimensional DMStag in each coordinate direction.

106:   One should use DMStagGetProductCoordinateLocationSlot() to determine appropriate
107:   indices for the second dimension in these returned arrays. This function
108:   checks that the coordinate array is a suitable product of 1-dimensional
109:   DMStag objects.

111:   Input Parameter:
112: . dm - the DMStag object

114:   Output Parameters:
115: . arrX,arrY,arrZ - local 1D coordinate arrays

117:   Level: intermediate

119: .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArraysRead(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagGetProductCoordinateLocationSlot()
120: @*/
121: PetscErrorCode DMStagGetProductCoordinateArrays(DM dm,void* arrX,void* arrY,void* arrZ)
122: {

126:   DMStagGetProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_FALSE);
127:   return(0);
128: }

130: /*@C
131:   DMStagGetProductCoordinateArraysRead - extract product coordinate arrays, read-only

133:   Logically Collective

135:   See the man page for DMStagGetProductCoordinateArrays() for more information.

137:   Input Parameter:
138: . dm - the DMStag object

140:   Output Parameters:
141: . arrX,arrY,arrZ - local 1D coordinate arrays

143:   Level: intermediate

145: .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArrays(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagGetProductCoordinateLocationSlot()
146: @*/
147: PetscErrorCode DMStagGetProductCoordinateArraysRead(DM dm,void* arrX,void* arrY,void* arrZ)
148: {

152:   DMStagGetProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_TRUE);
153:   return(0);
154: }

156: /*@C
157:   DMStagGetProductCoordinateLocationSlot - get slot for use with local product coordinate arrays

159:   Not Collective

161:   High-level helper function to get slot indices for 1D coordinate DMs,
162:   for use with DMStagGetProductCoordinateArrays() and related functions.

164:   Input Parameters:
165: + dm - the DMStag object
166: - loc - the grid location

168:   Output Parameter:
169: . slot - the index to use in local arrays

171:   Notes:
172:   For loc, one should use DMSTAG_LEFT, DMSTAG_ELEMENT, or DMSTAG_RIGHT for "previous", "center" and "next"
173:   locations, respectively, in each dimension.
174:   One can equivalently use DMSTAG_DOWN or DMSTAG_BACK in place of DMSTAG_LEFT,
175:   and DMSTAG_UP or DMSTACK_FRONT in place of DMSTAG_RIGHT;

177:   This function checks that the coordinates are actually set up so that using the
178:   slots from any of the 1D coordinate sub-DMs are valid for all the 1D coordinate sub-DMs.

180:   Level: intermediate

182: .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead(), DMStagSetUniformCoordinates()
183: @*/
184: PETSC_EXTERN PetscErrorCode DMStagGetProductCoordinateLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt *slot)
185: {
187:   DM             dmCoord;
188:   PetscInt       dim,dofCheck[DMSTAG_MAX_STRATA],s,d;

192:   DMGetDimension(dm,&dim);
193:   DMGetCoordinateDM(dm,&dmCoord);
194:   if (!dmCoord) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM does not have a coordinate DM");
195:   {
196:     PetscBool isProduct;
197:     DMType    dmType;
198:     DMGetType(dmCoord,&dmType);
199:     PetscStrcmp(DMPRODUCT,dmType,&isProduct);
200:     if (!isProduct) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is not of type DMPRODUCT");
201:   }
202:   for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
203:   for (d=0; d<dim; ++d) {
204:     DM        subDM;
205:     DMType    dmType;
206:     PetscBool isStag;
207:     PetscInt  dof[DMSTAG_MAX_STRATA],subDim;
208:     DMProductGetDM(dmCoord,d,&subDM);
209:     if (!subDM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is missing sub DM %D",d);
210:     DMGetDimension(subDM,&subDim);
211:     if (subDim != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of dimension 1");
212:     DMGetType(subDM,&dmType);
213:     PetscStrcmp(DMSTAG,dmType,&isStag);
214:     if (!isStag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of type DMSTAG");
215:     DMStagGetDOF(subDM,&dof[0],&dof[1],&dof[2],&dof[3]);
216:     if (d == 0) {
217:       const PetscInt component = 0;
218:       for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
219:       DMStagGetLocationSlot(subDM,loc,component,slot);
220:     } else {
221:       for (s=0; s<DMSTAG_MAX_STRATA; ++s) {
222:         if (dofCheck[s] != dof[s]) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DMs have different dofs");
223:       }
224:     }
225:   }
226:   return(0);
227: }

229: /*@C
230:   DMStagGetCorners - return global element indices of the local region (excluding ghost points)

232:   Not Collective

234:   Input Parameter:
235: . dm - the DMStag object

237:   Output Parameters:
238: + x     - starting element index in first direction
239: . y     - starting element index in second direction
240: . z     - starting element index in third direction
241: . m     - element width in first direction
242: . n     - element width in second direction
243: . p     - element width in third direction
244: . nExtrax - number of extra partial elements in first direction
245: . nExtray - number of extra partial elements in second direction
246: - nExtraz - number of extra partial elements in third direction

248:   Notes:
249:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.

251:   The number of extra partial elements is either 1 or 0.
252:   The value is 1 on right, top, and front non-periodic domain ("physical") boundaries,
253:   in the x, y, and z directions respectively, and otherwise 0.

255:   Level: beginner

257: .seealso: DMSTAG, DMStagGetGhostCorners(), DMDAGetCorners()
258: @*/
259: PetscErrorCode DMStagGetCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *nExtrax,PetscInt *nExtray,PetscInt *nExtraz)
260: {
261:   const DM_Stag * const stag = (DM_Stag*)dm->data;

265:   if (x) *x = stag->start[0];
266:   if (y) *y = stag->start[1];
267:   if (z) *z = stag->start[2];
268:   if (m) *m = stag->n[0];
269:   if (n) *n = stag->n[1];
270:   if (p) *p = stag->n[2];
271:   if (nExtrax) *nExtrax = stag->boundaryType[0] != DM_BOUNDARY_PERIODIC && stag->lastRank[0] ? 1 : 0;
272:   if (nExtray) *nExtray = stag->boundaryType[1] != DM_BOUNDARY_PERIODIC && stag->lastRank[1] ? 1 : 0;
273:   if (nExtraz) *nExtraz = stag->boundaryType[2] != DM_BOUNDARY_PERIODIC && stag->lastRank[2] ? 1 : 0;
274:   return(0);
275: }

277: /*@C
278:   DMStagGetDOF - get number of DOF associated with each stratum of the grid

280:   Not Collective

282:   Input Parameter:
283: . dm - the DMStag object

285:   Output Parameters:
286: + dof0 - the number of points per 0-cell (vertex/node)
287: . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
288: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
289: - dof3 - the number of points per 3-cell (element in 3D)

291:   Level: beginner

293: .seealso: DMSTAG, DMStagGetCorners(), DMStagGetGhostCorners(), DMStagGetGlobalSizes(), DMStagGetStencilWidth(), DMStagGetBoundaryTypes(), DMStagGetLocationDOF(), DMDAGetDof()
294: @*/
295: PetscErrorCode DMStagGetDOF(DM dm,PetscInt *dof0,PetscInt *dof1,PetscInt *dof2,PetscInt *dof3)
296: {
297:   const DM_Stag * const stag = (DM_Stag*)dm->data;

301:   if (dof0) *dof0 = stag->dof[0];
302:   if (dof1) *dof1 = stag->dof[1];
303:   if (dof2) *dof2 = stag->dof[2];
304:   if (dof3) *dof3 = stag->dof[3];
305:   return(0);
306: }

308: /*@C
309:   DMStagGetGhostCorners - return global element indices of the local region, including ghost points

311:   Not Collective

313:   Input Parameter:
314: . dm - the DMStag object

316:   Output Parameters:
317: + x - the starting element index in the first direction
318: . y - the starting element index in the second direction
319: . z - the starting element index in the third direction
320: . m - the element width in the first direction
321: . n - the element width in the second direction
322: - p - the element width in the third direction

324:   Notes:
325:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.

327:   Level: beginner

329: .seealso: DMSTAG, DMStagGetCorners(), DMDAGetGhostCorners()
330: @*/
331: PetscErrorCode DMStagGetGhostCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
332: {
333:   const DM_Stag * const stag = (DM_Stag*)dm->data;

337:   if (x) *x = stag->startGhost[0];
338:   if (y) *y = stag->startGhost[1];
339:   if (z) *z = stag->startGhost[2];
340:   if (m) *m = stag->nGhost[0];
341:   if (n) *n = stag->nGhost[1];
342:   if (p) *p = stag->nGhost[2];
343:   return(0);
344: }

346: /*@C
347:   DMStagGetGlobalSizes - get global element counts

349:   Not Collective

351:   Input Parameter:
352: . dm - the DMStag object

354:   Output Parameters:
355: . M,N,P - global element counts in each direction

357:   Notes:
358:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.

360:   Level: beginner

362: .seealso: DMSTAG, DMStagGetLocalSizes(), DMDAGetInfo()
363: @*/
364: PetscErrorCode DMStagGetGlobalSizes(DM dm,PetscInt* M,PetscInt* N,PetscInt* P)
365: {
366:   const DM_Stag * const stag = (DM_Stag*)dm->data;

370:   if (M) *M = stag->N[0];
371:   if (N) *N = stag->N[1];
372:   if (P) *P = stag->N[2];
373:   return(0);
374: }

376: /*@C
377:   DMStagGetIsFirstRank - get boolean value for whether this rank is first in each direction in the rank grid

379:   Not Collective

381:   Input Parameter:
382: . dm - the DMStag object

384:   Output Parameters:
385: . isFirstRank0,isFirstRank1,isFirstRank2 - whether this rank is first in each direction

387:   Level: intermediate

389:   Notes:
390:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.

392: .seealso: DMSTAG, DMStagGetIsLastRank()
393: @*/
394: PetscErrorCode DMStagGetIsFirstRank(DM dm,PetscBool *isFirstRank0,PetscBool *isFirstRank1,PetscBool *isFirstRank2)
395: {
396:   const DM_Stag * const stag = (DM_Stag*)dm->data;

400:   if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
401:   if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
402:   if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
403:   return(0);
404: }

406: /*@C
407:   DMStagGetIsLastRank - get boolean value for whether this rank is last in each direction in the rank grid

409:   Not Collective

411:   Input Parameter:
412: . dm - the DMStag object

414:   Output Parameters:
415: . isLastRank0,isLastRank1,isLastRank2 - whether this rank is last in each direction

417:   Level: intermediate

419:   Notes:
420:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
421:   Level: intermediate

423: .seealso: DMSTAG, DMStagGetIsFirstRank()
424: @*/
425: PetscErrorCode DMStagGetIsLastRank(DM dm,PetscBool *isLastRank0,PetscBool *isLastRank1,PetscBool *isLastRank2)
426: {
427:   const DM_Stag * const stag = (DM_Stag*)dm->data;

431:   if (isLastRank0) *isLastRank0 = stag->lastRank[0];
432:   if (isLastRank1) *isLastRank1 = stag->lastRank[1];
433:   if (isLastRank2) *isLastRank2 = stag->lastRank[2];
434:   return(0);
435: }

437: /*@C
438:   DMStagGetLocalSizes - get local elementwise sizes

440:   Not Collective

442:   Input Parameter:
443: . dm - the DMStag object

445:   Output Parameters:
446: . m,n,p - local element counts (excluding ghosts) in each direction

448:   Notes:
449:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
450:   Level: intermediate

452:   Level: beginner

454: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetDOF(), DMStagGetNumRanks(), DMDAGetLocalInfo()
455: @*/
456: PetscErrorCode DMStagGetLocalSizes(DM dm,PetscInt* m,PetscInt* n,PetscInt* p)
457: {
458:   const DM_Stag * const stag = (DM_Stag*)dm->data;

462:   if (m) *m = stag->n[0];
463:   if (n) *n = stag->n[1];
464:   if (p) *p = stag->n[2];
465:   return(0);
466: }

468: /*@C
469:   DMStagGetNumRanks - get number of ranks in each direction in the global grid decomposition

471:   Not Collective

473:   Input Parameter:
474: . dm - the DMStag object

476:   Output Parameters:
477: . nRanks0,nRanks1,nRanks2 - number of ranks in each direction in the grid decomposition

479:   Notes:
480:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
481:   Level: intermediate

483:   Level: beginner

485: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetLocalSize(), DMStagSetNumRanks(), DMDAGetInfo()
486: @*/
487: PetscErrorCode DMStagGetNumRanks(DM dm,PetscInt *nRanks0,PetscInt *nRanks1,PetscInt *nRanks2)
488: {
489:   const DM_Stag * const stag = (DM_Stag*)dm->data;

493:   if (nRanks0) *nRanks0 = stag->nRanks[0];
494:   if (nRanks1) *nRanks1 = stag->nRanks[1];
495:   if (nRanks2) *nRanks2 = stag->nRanks[2];
496:   return(0);
497: }

499: /*@C
500:   DMStagGetEntries - get number of native entries in the global representation

502:   Not Collective

504:   Input Parameter:
505: . dm - the DMStag object

507:   Output Parameters:
508: . entries - number of rank-native entries in the global representation

510:   Note:
511:   This is the number of entries on this rank for a global vector associated with dm.

513:   Level: developer

515: .seealso: DMSTAG, DMStagGetDOF(), DMStagGetEntriesPerElement(), DMCreateLocalVector()
516: @*/
517: PetscErrorCode DMStagGetEntries(DM dm,PetscInt *entries)
518: {
519:   const DM_Stag * const stag = (DM_Stag*)dm->data;

523:   if (entries) *entries = stag->entries;
524:   return(0);
525: }

527: /*@C
528:   DMStagGetEntriesPerElement - get number of entries per element in the local representation

530:   Not Collective

532:   Input Parameter:
533: . dm - the DMStag object

535:   Output Parameters:
536: . entriesPerElement - number of entries associated with each element in the local representation

538:   Notes:
539:   This is the natural block size for most local operations. In 1D it is equal to dof0 + dof1,
540:   in 2D it is equal to dof0 + 2*dof1 + dof2, and in 3D it is equal to dof0 + 3*dof1 + 3*dof2 + dof3

542:   Level: developer

544: .seealso: DMSTAG, DMStagGetDOF()
545: @*/
546: PetscErrorCode DMStagGetEntriesPerElement(DM dm,PetscInt *entriesPerElement)
547: {
548:   const DM_Stag * const stag = (DM_Stag*)dm->data;

552:   if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
553:   return(0);
554: }

556: /*@C
557:   DMStagGetStencilType - get elementwise ghost/halo stencil type

559:   Not Collective

561:   Input Parameter:
562: . dm - the DMStag object

564:   Output Parameter:
565: . stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE

567:   Level: beginner

569: .seealso: DMSTAG, DMStagSetStencilType(), DMStagGetStencilWidth, DMStagStencilType
570: @*/
571: PetscErrorCode DMStagGetStencilType(DM dm,DMStagStencilType *stencilType)
572: {
573:   DM_Stag * const stag = (DM_Stag*)dm->data;

577:   *stencilType = stag->stencilType;
578:   return(0);
579: }

581: /*@C
582:   DMStagGetStencilWidth - get elementwise stencil width

584:   Not Collective

586:   Input Parameter:
587: . dm - the DMStag object

589:   Output Parameters:
590: . stencilWidth - stencil/halo/ghost width in elements

592:   Level: beginner

594: .seealso: DMSTAG, DMStagSetStencilWidth(), DMStagGetStencilType(), DMDAGetStencilType()
595: @*/
596: PetscErrorCode DMStagGetStencilWidth(DM dm,PetscInt *stencilWidth)
597: {
598:   const DM_Stag * const stag = (DM_Stag*)dm->data;

602:   if (stencilWidth) *stencilWidth = stag->stencilWidth;
603:   return(0);
604: }

606: /*@C
607:   DMStagGetOwnershipRanges - get elements per rank in each direction

609:   Not Collective

611:   Input Parameter:
612: .     dm - the DMStag object

614:   Output Parameters:
615: +     lx - ownership along x direction (optional)
616: .     ly - ownership along y direction (optional)
617: -     lz - ownership along z direction (optional)

619:   Notes:
620:   These correspond to the optional final arguments passed to DMStagCreate1d(), DMStagCreate2d(), and DMStagCreate3d().

622:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.

624:   In C you should not free these arrays, nor change the values in them.
625:   They will only have valid values while the DMStag they came from still exists (has not been destroyed).

627:   Level: intermediate

629: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagSetOwnershipRanges(), DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDAGetOwnershipRanges()
630: @*/
631: PetscErrorCode DMStagGetOwnershipRanges(DM dm,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
632: {
633:   const DM_Stag * const stag = (DM_Stag*)dm->data;

637:   if (lx) *lx = stag->l[0];
638:   if (ly) *ly = stag->l[1];
639:   if (lz) *lz = stag->l[2];
640:   return(0);
641: }

643: /*@C
644:   DMStagCreateCompatibleDMStag - create a compatible DMStag with different dof/stratum

646:   Collective

648:   Input Parameters:
649: + dm - the DMStag object
650: - dof0,dof1,dof2,dof3 - number of dof on each stratum in the new DMStag

652:   Output Parameters:
653: . newdm - the new, compatible DMStag

655:   Notes:
656:   Dof supplied for strata too big for the dimension are ignored; these may be set to 0.
657:   For example, for a 2-dimensional DMStag, dof2 sets the number of dof per element,
658:   and dof3 is unused. For a 3-dimensional DMStag, dof3 sets the number of dof per element.

660:   In contrast to DMDACreateCompatibleDMDA(), coordinates are not reused.

662:   Level: intermediate

664: .seealso: DMSTAG, DMDACreateCompatibleDMDA(), DMGetCompatibility(), DMStagMigrateVec()
665: @*/
666: PetscErrorCode DMStagCreateCompatibleDMStag(DM dm,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DM *newdm)
667: {
668:   PetscErrorCode  ierr;

672:   DMStagDuplicateWithoutSetup(dm,PetscObjectComm((PetscObject)dm),newdm);
673:   DMStagSetDOF(*newdm,dof0,dof1,dof2,dof3);
674:   DMSetUp(*newdm);
675:   return(0);
676: }

678: /*@C
679:   DMStagGetLocationSlot - get index to use in accessing raw local arrays

681:   Not Collective

683:   Input Parameters:
684: + dm - the DMStag object
685: . loc - location relative to an element
686: - c - component

688:   Output Parameter:
689: . slot - index to use

691:   Notes:
692:   Provides an appropriate index to use with DMStagVecGetArray() and friends.
693:   This is required so that the user doesn't need to know about the ordering of
694:   dof associated with each local element.

696:   Level: beginner

698: .seealso: DMSTAG, DMStagVecGetArray(), DMStagVecGetArrayRead(), DMStagGetDOF(), DMStagGetEntriesPerElement()
699: @*/
700: PetscErrorCode DMStagGetLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt c,PetscInt *slot)
701: {
702:   DM_Stag * const stag = (DM_Stag*)dm->data;

706:   if (PetscDefined(USE_DEBUG)) {
708:     PetscInt       dof;
709:     DMStagGetLocationDOF(dm,loc,&dof);
710:     if (dof < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Location %s has no dof attached",DMStagStencilLocations[loc]);
711:     if (c > dof-1) SETERRQ3(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Supplied component number (%D) for location %s is too big (maximum %D)",c,DMStagStencilLocations[loc],dof-1);
712:   }
713:   *slot = stag->locationOffsets[loc] + c;
714:   return(0);
715: }

717: /*@C
718:   DMStagMigrateVec - transfer a vector associated with a DMStag to a vector associated with a compatible DMStag

720:   Collective

722:   Input Parameters:
723: + dm - the source DMStag object
724: . vec - the source vector, compatible with dm
725: . dmTo - the compatible destination DMStag object
726: - vecTo - the destination vector, compatible with dmTo

728:   Notes:
729:   Extra dof are ignored, and unfilled dof are zeroed.
730:   Currently only implemented to migrate global vectors to global vectors.

732:   Level: advanced

734: .seealso: DMSTAG, DMStagCreateCompatibleDMStag(), DMGetCompatibility(), DMStagVecSplitToDMDA()
735: @*/
736: PetscErrorCode DMStagMigrateVec(DM dm,Vec vec,DM dmTo,Vec vecTo)
737: {
738:   PetscErrorCode    ierr;
739:   DM_Stag * const   stag = (DM_Stag*)dm->data;
740:   DM_Stag * const   stagTo = (DM_Stag*)dmTo->data;
741:   PetscInt          nLocalTo,nLocal,dim,i,j,k;
742:   PetscInt          start[DMSTAG_MAX_DIM],startGhost[DMSTAG_MAX_DIM],n[DMSTAG_MAX_DIM],nExtra[DMSTAG_MAX_DIM],offset[DMSTAG_MAX_DIM];
743:   Vec               vecToLocal,vecLocal;
744:   PetscBool         compatible,compatibleSet;
745:   const PetscScalar *arr;
746:   PetscScalar       *arrTo;
747:   const PetscInt    epe   = stag->entriesPerElement;
748:   const PetscInt    epeTo = stagTo->entriesPerElement;

755:   DMGetCompatibility(dm,dmTo,&compatible,&compatibleSet);
756:   if (!compatibleSet || !compatible) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"DMStag objects must be shown to be compatible");
757:   DMGetDimension(dm,&dim);
758:   VecGetLocalSize(vec,&nLocal);
759:   VecGetLocalSize(vecTo,&nLocalTo);
760:   if (nLocal != stag->entries|| nLocalTo !=stagTo->entries) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Vector migration only implemented for global vector to global vector.");
761:   DMStagGetCorners(dm,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
762:   DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],NULL,NULL,NULL);
763:   for (i=0; i<DMSTAG_MAX_DIM; ++i) offset[i] = start[i]-startGhost[i];

765:   /* Proceed by transferring to a local vector, copying, and transferring back to a global vector */
766:   DMGetLocalVector(dm,&vecLocal);
767:   DMGetLocalVector(dmTo,&vecToLocal);
768:   DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal);
769:   DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal);
770:   VecGetArrayRead(vecLocal,&arr);
771:   VecGetArray(vecToLocal,&arrTo);
772:   /* Note that some superfluous copying of entries on partial dummy elements is done */
773:   if (dim == 1) {
774:     for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
775:       PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
776:       PetscInt si;
777:       for (si=0; si<2; ++si) {
778:         b   += stag->dof[si];
779:         bTo += stagTo->dof[si];
780:         for (; d < b && dTo < bTo; ++d,++dTo) arrTo[i*epeTo + dTo] = arr[i*epe + d];
781:         for (; dTo < bTo         ;     ++dTo) arrTo[i*epeTo + dTo] = 0.0;
782:         d = b;
783:       }
784:     }
785:   } else if (dim == 2) {
786:     const PetscInt epr   = stag->nGhost[0] * epe;
787:     const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
788:     for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
789:       for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
790:         const PetscInt base   = j*epr   + i*epe;
791:         const PetscInt baseTo = j*eprTo + i*epeTo;
792:         PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
793:         const PetscInt s[4] = {0,1,1,2}; /* Dimensions of points, in order */
794:         PetscInt si;
795:         for (si=0; si<4; ++si) {
796:             b   += stag->dof[s[si]];
797:             bTo += stagTo->dof[s[si]];
798:             for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
799:             for (;          dTo < bTo;     ++dTo) arrTo[baseTo + dTo] = 0.0;
800:             d = b;
801:         }
802:       }
803:     }
804:   } else if (dim == 3) {
805:     const PetscInt epr   = stag->nGhost[0]   * epe;
806:     const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
807:     const PetscInt epl   = stag->nGhost[1]   * epr;
808:     const PetscInt eplTo = stagTo->nGhost[1] * eprTo;
809:     for (k=offset[2]; k<offset[2] + n[2] + nExtra[2]; ++k) {
810:       for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
811:         for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
812:           PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
813:           const PetscInt base   = k*epl   + j*epr   + i*epe;
814:           const PetscInt baseTo = k*eplTo + j*eprTo + i*epeTo;
815:           const PetscInt s[8] = {0,1,1,2,1,2,2,3}; /* dimensions of points, in order */
816:           PetscInt is;
817:           for (is=0; is<8; ++is) {
818:             b   += stag->dof[s[is]];
819:             bTo += stagTo->dof[s[is]];
820:             for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
821:             for (;          dTo < bTo;     ++dTo) arrTo[baseTo + dTo] = 0.0;
822:             d = b;
823:           }
824:         }
825:       }
826:     }
827:   } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
828:   VecRestoreArrayRead(vecLocal,&arr);
829:   VecRestoreArray(vecToLocal,&arrTo);
830:   DMRestoreLocalVector(dm,&vecLocal);
831:   DMLocalToGlobalBegin(dmTo,vecToLocal,INSERT_VALUES,vecTo);
832:   DMLocalToGlobalEnd(dmTo,vecToLocal,INSERT_VALUES,vecTo);
833:   DMRestoreLocalVector(dmTo,&vecToLocal);
834:   return(0);
835: }

837: /*@C
838:   DMStagPopulateLocalToGlobalInjective - populate an internal 1-to-1 local-to-global map

840:   Collective

842:   Creates an internal object which explicitly maps a single local degree of
843:   freedom to each global degree of freedom. This is used, if populated,
844:   instead of SCATTER_REVERSE_LOCAL with the (1-to-many, in general)
845:   global-to-local map, when DMLocalToGlobal() is called with INSERT_VALUES.
846:   This allows usage, for example, even in the periodic, 1-rank case, where
847:   the inverse of the global-to-local map, even when restricted to on-rank
848:   communication, is non-injective. This is at the cost of storing an additional
849:   VecScatter object inside each DMStag object.

851:   Input Parameter:
852: . dm - the DMStag object

854:   Notes:
855:   In normal usage, library users shouldn't be concerned with this function,
856:   as it is called during DMSetUp(), when required.

858:   Returns immediately if the internal map is already populated.

860:   Developer Notes:
861:   This could, if desired, be moved up to a general DM routine. It would allow,
862:   for example, DMDA to support DMLocalToGlobal() with INSERT_VALUES,
863:   even in the single-rank periodic case.

865:   Level: developer

867: .seealso: DMSTAG, DMLocalToGlobal(), VecScatter
868: @*/
869: PetscErrorCode DMStagPopulateLocalToGlobalInjective(DM dm)
870: {
871:   PetscErrorCode  ierr;
872:   PetscInt        dim;
873:   DM_Stag * const stag  = (DM_Stag*)dm->data;

877:   if (stag->ltog_injective) return(0); /* Don't re-populate */
878:   DMGetDimension(dm,&dim);
879:   switch (dim) {
880:     case 1: DMStagPopulateLocalToGlobalInjective_1d(dm); break;
881:     case 2: DMStagPopulateLocalToGlobalInjective_2d(dm); break;
882:     case 3: DMStagPopulateLocalToGlobalInjective_3d(dm); break;
883:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
884:   }
885:   return(0);
886: }

888: static PetscErrorCode DMStagRestoreProductCoordinateArrays_Private(DM dm,void *arrX,void *arrY,void *arrZ,PetscBool read)
889: {
890:   PetscErrorCode  ierr;
891:   PetscInt        dim,d;
892:   void*           arr[DMSTAG_MAX_DIM];
893:   DM              dmCoord;

897:   DMGetDimension(dm,&dim);
898:   if (dim > DMSTAG_MAX_DIM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for %D dimensions",dim);
899:   arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
900:   DMGetCoordinateDM(dm,&dmCoord);
901:   for (d=0; d<dim; ++d) {
902:     DM  subDM;
903:     Vec coord1d_local;

905:     /* Ignore unrequested arrays */
906:     if (!arr[d]) continue;

908:     DMProductGetDM(dmCoord,d,&subDM);
909:     DMGetCoordinatesLocal(subDM,&coord1d_local);
910:     if (read) {
911:       DMStagVecRestoreArrayRead(subDM,coord1d_local,arr[d]);
912:     } else {
913:       DMStagVecRestoreArray(subDM,coord1d_local,arr[d]);
914:     }
915:   }
916:   return(0);
917: }

919: /*@C
920:   DMStagRestoreProductCoordinateArrays - restore local array access

922:   Logically Collective

924:   Input Parameter:
925: . dm - the DMStag object

927:   Output Parameters:
928: . arrX,arrY,arrZ - local 1D coordinate arrays

930:   Level: intermediate

932:   Notes:
933:   This function does not automatically perform a local->global scatter to populate global coordinates from the local coordinates. Thus, it may be required to explicitly perform these operations in some situations, as in the following partial example:

935: $   DMGetCoordinateDM(dm,&cdm);
936: $   for (d=0; d<3; ++d) {
937: $     DM  subdm;
938: $     Vec coor,coor_local;

940: $     DMProductGetDM(cdm,d,&subdm);
941: $     DMGetCoordinates(subdm,&coor);
942: $     DMGetCoordinatesLocal(subdm,&coor_local);
943: $     DMLocalToGlobal(subdm,coor_local,INSERT_VALUES,coor);
944: $     PetscPrintf(PETSC_COMM_WORLD,"Coordinates dim %D:\n",d);
945: $     VecView(coor,PETSC_VIEWER_STDOUT_WORLD);
946: $   }

948: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
949: @*/
950: PetscErrorCode DMStagRestoreProductCoordinateArrays(DM dm,void *arrX,void *arrY,void *arrZ)
951: {
952:   PetscErrorCode  ierr;

955:   DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_FALSE);
956:   return(0);
957: }

959: /*@C
960:   DMStagRestoreProductCoordinateArraysRead - restore local product array access, read-only

962:   Logically Collective

964:   Input Parameter:
965: . dm - the DMStag object

967:   Output Parameters:
968: . arrX,arrY,arrZ - local 1D coordinate arrays

970:   Level: intermediate

972: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
973: @*/
974: PetscErrorCode DMStagRestoreProductCoordinateArraysRead(DM dm,void *arrX,void *arrY,void *arrZ)
975: {
976:   PetscErrorCode  ierr;

979:   DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_TRUE);
980:   return(0);
981: }

983: /*@C
984:   DMStagSetBoundaryTypes - set DMStag boundary types

986:   Logically Collective; boundaryType0, boundaryType1, and boundaryType2 must contain common values

988:   Input Parameters:
989: + dm - the DMStag object
990: - boundaryType0,boundaryType1,boundaryType2 - boundary types in each direction

992:   Notes:
993:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

995:   Level: advanced

997: .seealso: DMSTAG, DMBoundaryType, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDASetBoundaryType()
998: @*/
999: PetscErrorCode DMStagSetBoundaryTypes(DM dm,DMBoundaryType boundaryType0,DMBoundaryType boundaryType1,DMBoundaryType boundaryType2)
1000: {
1001:   PetscErrorCode  ierr;
1002:   DM_Stag * const stag  = (DM_Stag*)dm->data;
1003:   PetscInt        dim;

1006:   DMGetDimension(dm,&dim);
1011:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1012:                stag->boundaryType[0] = boundaryType0;
1013:   if (dim > 1) stag->boundaryType[1] = boundaryType1;
1014:   if (dim > 2) stag->boundaryType[2] = boundaryType2;
1015:   return(0);
1016: }

1018: /*@C
1019:   DMStagSetCoordinateDMType - set DM type to store coordinates

1021:   Logically Collective; dmtype must contain common value

1023:   Input Parameters:
1024: + dm - the DMStag object
1025: - dmtype - DMtype for coordinates, either DMSTAG or DMPRODUCT

1027:   Level: advanced

1029: .seealso: DMSTAG, DMPRODUCT, DMGetCoordinateDM(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMType
1030: @*/
1031: PetscErrorCode DMStagSetCoordinateDMType(DM dm,DMType dmtype)
1032: {
1033:   PetscErrorCode  ierr;
1034:   DM_Stag * const stag = (DM_Stag*)dm->data;

1038:   PetscFree(stag->coordinateDMType);
1039:   PetscStrallocpy(dmtype,(char**)&stag->coordinateDMType);
1040:   return(0);
1041: }

1043: /*@C
1044:   DMStagSetDOF - set dof/stratum

1046:   Logically Collective; dof0, dof1, dof2, and dof3 must contain common values

1048:   Input Parameters:
1049: + dm - the DMStag object
1050: - dof0,dof1,dof2,dof3 - dof per stratum

1052:   Notes:
1053:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

1055:   Level: advanced

1057: .seealso: DMSTAG, DMDASetDof()
1058: @*/
1059: PetscErrorCode DMStagSetDOF(DM dm,PetscInt dof0, PetscInt dof1,PetscInt dof2,PetscInt dof3)
1060: {
1061:   PetscErrorCode  ierr;
1062:   DM_Stag * const stag = (DM_Stag*)dm->data;
1063:   PetscInt        dim;

1071:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1072:   DMGetDimension(dm,&dim);
1073:   if (dof0 < 0)            SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof0 cannot be negative");
1074:   if (dof1 < 0)            SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof1 cannot be negative");
1075:   if (dim > 1 && dof2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof2 cannot be negative");
1076:   if (dim > 2 && dof3 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof3 cannot be negative");
1077:                stag->dof[0] = dof0;
1078:                stag->dof[1] = dof1;
1079:   if (dim > 1) stag->dof[2] = dof2;
1080:   if (dim > 2) stag->dof[3] = dof3;
1081:   return(0);
1082: }

1084: /*@C
1085:   DMStagSetNumRanks - set ranks in each direction in the global rank grid

1087:   Logically Collective; nRanks0, nRanks1, and nRanks2 must contain common values

1089:   Input Parameters:
1090: + dm - the DMStag object
1091: - nRanks0,nRanks1,nRanks2 - number of ranks in each direction

1093:   Notes:
1094:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

1096:   Level: developer

1098: .seealso: DMSTAG, DMDASetNumProcs()
1099: @*/
1100: PetscErrorCode DMStagSetNumRanks(DM dm,PetscInt nRanks0,PetscInt nRanks1,PetscInt nRanks2)
1101: {
1102:   PetscErrorCode  ierr;
1103:   DM_Stag * const stag = (DM_Stag*)dm->data;
1104:   PetscInt        dim;

1111:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1112:   DMGetDimension(dm,&dim);
1113:   if (nRanks0 != PETSC_DECIDE && nRanks0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in X direction cannot be less than 1");
1114:   if (dim > 1 && nRanks1 != PETSC_DECIDE && nRanks1 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in Y direction cannot be less than 1");
1115:   if (dim > 2 && nRanks2 != PETSC_DECIDE && nRanks2 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in Z direction cannot be less than 1");
1116:   if (nRanks0) stag->nRanks[0] = nRanks0;
1117:   if (dim > 1 && nRanks1) stag->nRanks[1] = nRanks1;
1118:   if (dim > 2 && nRanks2) stag->nRanks[2] = nRanks2;
1119:   return(0);
1120: }

1122: /*@C
1123:   DMStagSetStencilType - set elementwise ghost/halo stencil type

1125:   Logically Collective; stencilType must contain common value

1127:   Input Parameters:
1128: + dm - the DMStag object
1129: - stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE

1131:   Level: beginner

1133: .seealso: DMSTAG, DMStagGetStencilType(), DMStagSetStencilWidth(), DMStagStencilType
1134: @*/
1135: PetscErrorCode DMStagSetStencilType(DM dm,DMStagStencilType stencilType)
1136: {
1137:   DM_Stag * const stag = (DM_Stag*)dm->data;

1142:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1143:   stag->stencilType = stencilType;
1144:   return(0);
1145: }

1147: /*@C
1148:   DMStagSetStencilWidth - set elementwise stencil width

1150:   Logically Collective; stencilWidth must contain common value

1152:   Input Parameters:
1153: + dm - the DMStag object
1154: - stencilWidth - stencil/halo/ghost width in elements

1156:   Notes:
1157:   The width value is not used when DMSTAG_STENCIL_NONE is specified.

1159:   Level: beginner

1161: .seealso: DMSTAG, DMStagGetStencilWidth(), DMStagGetStencilType(), DMStagStencilType
1162: @*/
1163: PetscErrorCode DMStagSetStencilWidth(DM dm,PetscInt stencilWidth)
1164: {
1165:   DM_Stag * const stag = (DM_Stag*)dm->data;

1170:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1171:   if (stencilWidth < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Stencil width must be non-negative");
1172:   stag->stencilWidth = stencilWidth;
1173:   return(0);
1174: }

1176: /*@C
1177:   DMStagSetGlobalSizes - set global element counts in each direction

1179:   Logically Collective; N0, N1, and N2 must contain common values

1181:   Input Parameters:
1182: + dm - the DMStag object
1183: - N0,N1,N2 - global elementwise sizes

1185:   Notes:
1186:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

1188:   Level: advanced

1190: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMDASetSizes()
1191: @*/
1192: PetscErrorCode DMStagSetGlobalSizes(DM dm,PetscInt N0,PetscInt N1,PetscInt N2)
1193: {
1194:   PetscErrorCode  ierr;
1195:   DM_Stag * const stag = (DM_Stag*)dm->data;
1196:   PetscInt        dim;

1200:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1201:   DMGetDimension(dm,&dim);
1202:   if (N0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in X direction must be positive");
1203:   if (dim > 1 && N1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Y direction must be positive");
1204:   if (dim > 2 && N2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Z direction must be positive");
1205:   if (N0) stag->N[0] = N0;
1206:   if (N1) stag->N[1] = N1;
1207:   if (N2) stag->N[2] = N2;
1208:   return(0);
1209: }

1211: /*@C
1212:   DMStagSetOwnershipRanges - set elements per rank in each direction

1214:   Logically Collective; lx, ly, and lz must contain common values

1216:   Input Parameters:
1217: + dm - the DMStag object
1218: - lx,ly,lz - element counts for each rank in each direction

1220:   Notes:
1221:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.

1223:   Level: developer

1225: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagGetOwnershipRanges(), DMDASetOwnershipRanges()
1226: @*/
1227: PetscErrorCode DMStagSetOwnershipRanges(DM dm,PetscInt const *lx,PetscInt const *ly,PetscInt const *lz)
1228: {
1229:   PetscErrorCode  ierr;
1230:   DM_Stag * const stag = (DM_Stag*)dm->data;
1231:   const PetscInt  *lin[3];
1232:   PetscInt        d,dim;

1236:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1237:   lin[0] = lx; lin[1] = ly; lin[2] = lz;
1238:   DMGetDimension(dm,&dim);
1239:   for (d=0; d<dim; ++d) {
1240:     if (lin[d]) {
1241:       if (stag->nRanks[d] < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of ranks");
1242:       if (!stag->l[d]) {
1243:         PetscMalloc1(stag->nRanks[d], &stag->l[d]);
1244:       }
1245:       PetscArraycpy(stag->l[d], lin[d], stag->nRanks[d]);
1246:     }
1247:   }
1248:   return(0);
1249: }

1251: /*@C
1252:   DMStagSetUniformCoordinates - set DMStag coordinates to be a uniform grid

1254:   Collective

1256:   Input Parameters:
1257: + dm - the DMStag object
1258: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

1260:   Notes:
1261:   DMStag supports 2 different types of coordinate DM: DMSTAG and DMPRODUCT.
1262:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

1264:   Local coordinates are populated (using DMSetCoordinatesLocal()), linearly
1265:   extrapolated to ghost cells, including those outside the physical domain.
1266:   This is also done in case of periodic boundaries, meaning that the same
1267:   global point may have different coordinates in different local representations,
1268:   which are equivalent assuming a periodicity implied by the arguments to this function,
1269:   i.e. two points are equivalent if their difference is a multiple of (xmax - xmin)
1270:   in the x direction, (ymax - ymin) in the y direction, and (zmax - zmin) in the z direction.

1272:   Level: advanced

1274: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType(), DMGetCoordinateDM(), DMGetCoordinates(), DMDASetUniformCoordinates(), DMBoundaryType
1275: @*/
1276: PetscErrorCode DMStagSetUniformCoordinates(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1277: {
1278:   PetscErrorCode  ierr;
1279:   DM_Stag * const stag = (DM_Stag*)dm->data;
1280:   PetscBool       flg_stag,flg_product;

1284:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1285:   if (!stag->coordinateDMType) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"You must first call DMStagSetCoordinateDMType()");
1286:   PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg_stag);
1287:   PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg_product);
1288:   if (flg_stag) {
1289:     DMStagSetUniformCoordinatesExplicit(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1290:   } else if (flg_product) {
1291:     DMStagSetUniformCoordinatesProduct(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1292:   } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported DM Type %s",stag->coordinateDMType);
1293:   return(0);
1294: }

1296: /*@C
1297:   DMStagSetUniformCoordinatesExplicit - set DMStag coordinates to be a uniform grid, storing all values

1299:   Collective

1301:   Input Parameters:
1302: + dm - the DMStag object
1303: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

1305:   Notes:
1306:   DMStag supports 2 different types of coordinate DM: either another DMStag, or a DMProduct.
1307:   If the grid is orthogonal, using DMProduct should be more efficient.

1309:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

1311:   See the manual page for DMStagSetUniformCoordinates() for information on how
1312:   coordinates for dummy cells outside the physical domain boundary are populated.

1314:   Level: beginner

1316: .seealso: DMSTAG, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType()
1317: @*/
1318: PetscErrorCode DMStagSetUniformCoordinatesExplicit(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1319: {
1320:   PetscErrorCode  ierr;
1321:   DM_Stag * const stag = (DM_Stag*)dm->data;
1322:   PetscInt        dim;
1323:   PetscBool       flg;

1327:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1328:   PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg);
1329:   if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1330:   DMStagSetCoordinateDMType(dm,DMSTAG);
1331:   DMGetDimension(dm,&dim);
1332:   switch (dim) {
1333:     case 1: DMStagSetUniformCoordinatesExplicit_1d(dm,xmin,xmax);                     break;
1334:     case 2: DMStagSetUniformCoordinatesExplicit_2d(dm,xmin,xmax,ymin,ymax);           break;
1335:     case 3: DMStagSetUniformCoordinatesExplicit_3d(dm,xmin,xmax,ymin,ymax,zmin,zmax); break;
1336:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1337:   }
1338:   return(0);
1339: }

1341: /*@C
1342:   DMStagSetUniformCoordinatesProduct - create uniform coordinates, as a product of 1D arrays

1344:   Set the coordinate DM to be a DMProduct of 1D DMStag objects, each of which have a coordinate DM (also a 1d DMStag) holding uniform coordinates.

1346:   Collective

1348:   Input Parameters:
1349: + dm - the DMStag object
1350: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

1352:   Notes:
1353:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

1355:   The per-dimension 1-dimensional DMStag objects that comprise the product
1356:   always have active 0-cells (vertices, element boundaries) and 1-cells
1357:   (element centers).

1359:   See the manual page for DMStagSetUniformCoordinates() for information on how
1360:   coordinates for dummy cells outside the physical domain boundary are populated.

1362:   Level: intermediate

1364: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetCoordinateDMType()
1365: @*/
1366: PetscErrorCode DMStagSetUniformCoordinatesProduct(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1367: {
1368:   PetscErrorCode  ierr;
1369:   DM_Stag * const stag = (DM_Stag*)dm->data;
1370:   DM              dmc;
1371:   PetscInt        dim,d,dof0,dof1;
1372:   PetscBool       flg;

1376:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1377:   PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg);
1378:   if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1379:   DMStagSetCoordinateDMType(dm,DMPRODUCT);
1380:   DMGetCoordinateDM(dm,&dmc);
1381:   DMGetDimension(dm,&dim);

1383:   /* Create 1D sub-DMs, living on subcommunicators.
1384:      Always include both vertex and element dof, regardless of the active strata of the DMStag */
1385:   dof0 = 1;
1386:   dof1 = 1;

1388:   for (d=0; d<dim; ++d) {
1389:     DM                subdm;
1390:     MPI_Comm          subcomm;
1391:     PetscMPIInt       color;
1392:     const PetscMPIInt key = 0; /* let existing rank break ties */

1394:     /* Choose colors based on position in the plane orthogonal to this dim, and split */
1395:     switch (d) {
1396:       case 0: color = (dim > 1 ? stag->rank[1] : 0)  + (dim > 2 ? stag->nRanks[1]*stag->rank[2] : 0); break;
1397:       case 1: color =            stag->rank[0]       + (dim > 2 ? stag->nRanks[0]*stag->rank[2] : 0); break;
1398:       case 2: color =            stag->rank[0]       +            stag->nRanks[0]*stag->rank[1]     ; break;
1399:       default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1400:     }
1401:     MPI_Comm_split(PetscObjectComm((PetscObject)dm),color,key,&subcomm);

1403:     /* Create sub-DMs living on these new communicators (which are destroyed by DMProduct) */
1404:     DMStagCreate1d(subcomm,stag->boundaryType[d],stag->N[d],dof0,dof1,stag->stencilType,stag->stencilWidth,stag->l[d],&subdm);
1405:     DMSetUp(subdm);
1406:     switch (d) {
1407:       case 0:
1408:         DMStagSetUniformCoordinatesExplicit(subdm,xmin,xmax,0,0,0,0);
1409:         break;
1410:       case 1:
1411:         DMStagSetUniformCoordinatesExplicit(subdm,ymin,ymax,0,0,0,0);
1412:         break;
1413:       case 2:
1414:         DMStagSetUniformCoordinatesExplicit(subdm,zmin,zmax,0,0,0,0);
1415:         break;
1416:       default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1417:     }
1418:     DMProductSetDM(dmc,d,subdm);
1419:     DMProductSetDimensionIndex(dmc,d,0);
1420:     DMDestroy(&subdm);
1421:     MPI_Comm_free(&subcomm);
1422:   }
1423:   return(0);
1424: }

1426: /*@C
1427:   DMStagVecGetArray - get access to local array

1429:   Logically Collective

1431:   This function returns a (dim+1)-dimensional array for a dim-dimensional
1432:   DMStag.

1434:   The first 1-3 dimensions indicate an element in the global
1435:   numbering, using the standard C ordering.

1437:   The final dimension in this array corresponds to a degree
1438:   of freedom with respect to this element, for example corresponding to
1439:   the element or one of its neighboring faces, edges, or vertices.

1441:   For example, for a 3D DMStag, indexing is array[k][j][i][idx], where k is the
1442:   index in the z-direction, j is the index in the y-direction, and i is the
1443:   index in the x-direction.

1445:   "idx" is obtained with DMStagGetLocationSlot(), since the correct offset
1446:   into the (dim+1)-dimensional C array depends on the grid size and the number
1447:   of dof stored at each location.

1449:   Input Parameters:
1450: + dm - the DMStag object
1451: - vec - the Vec object

1453:   Output Parameters:
1454: . array - the array

1456:   Notes:
1457:   DMStagVecRestoreArray() must be called, once finished with the array

1459:   Level: beginner

1461: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArray(), DMDAVecGetArrayDOF()
1462: @*/
1463: PetscErrorCode DMStagVecGetArray(DM dm,Vec vec,void *array)
1464: {
1465:   PetscErrorCode  ierr;
1466:   DM_Stag * const stag = (DM_Stag*)dm->data;
1467:   PetscInt        dim;
1468:   PetscInt        nLocal;

1473:   DMGetDimension(dm,&dim);
1474:   VecGetLocalSize(vec,&nLocal);
1475:   if (nLocal != stag->entriesGhost) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMStag local size %D\n",nLocal,stag->entriesGhost);
1476:   switch (dim) {
1477:     case 1:
1478:       VecGetArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1479:       break;
1480:     case 2:
1481:       VecGetArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1482:       break;
1483:     case 3:
1484:       VecGetArray4d(vec,stag->nGhost[2],stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[2],stag->startGhost[1],stag->startGhost[0],0,(PetscScalar*****)array);
1485:       break;
1486:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1487:   }
1488:   return(0);
1489: }

1491: /*@C
1492:   DMStagVecGetArrayRead - get read-only access to a local array

1494:   Logically Collective

1496:   See the man page for DMStagVecGetArray() for more information.

1498:   Input Parameters:
1499: + dm - the DMStag object
1500: - vec - the Vec object

1502:   Output Parameters:
1503: . array - the read-only array

1505:   Notes:
1506:   DMStagVecRestoreArrayRead() must be called, once finished with the array

1508:   Level: beginner

1510: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArrayRead(), DMDAVecGetArrayDOFRead()
1511: @*/
1512: PetscErrorCode DMStagVecGetArrayRead(DM dm,Vec vec,void *array)
1513: {
1514:   PetscErrorCode  ierr;
1515:   DM_Stag * const stag = (DM_Stag*)dm->data;
1516:   PetscInt        dim;
1517:   PetscInt        nLocal;

1522:   DMGetDimension(dm,&dim);
1523:   VecGetLocalSize(vec,&nLocal);
1524:   if (nLocal != stag->entriesGhost) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMStag local size %D\n",nLocal,stag->entriesGhost);
1525:   switch (dim) {
1526:     case 1:
1527:       VecGetArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1528:       break;
1529:     case 2:
1530:       VecGetArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1531:       break;
1532:     case 3:
1533:       VecGetArray4dRead(vec,stag->nGhost[2],stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[2],stag->startGhost[1],stag->startGhost[0],0,(PetscScalar*****)array);
1534:       break;
1535:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1536:   }
1537:   return(0);
1538: }

1540: /*@C
1541:   DMStagVecRestoreArray - restore access to a raw array

1543:   Logically Collective

1545:   Input Parameters:
1546: + dm - the DMStag object
1547: - vec - the Vec object

1549:   Output Parameters:
1550: . array - the array

1552:   Level: beginner

1554: .seealso: DMSTAG, DMStagVecGetArray(), DMDAVecRestoreArray(), DMDAVecRestoreArrayDOF()
1555: @*/
1556: PetscErrorCode DMStagVecRestoreArray(DM dm,Vec vec,void *array)
1557: {
1558:   PetscErrorCode  ierr;
1559:   DM_Stag * const stag = (DM_Stag*)dm->data;
1560:   PetscInt        dim;
1561:   PetscInt        nLocal;

1566:   DMGetDimension(dm,&dim);
1567:   VecGetLocalSize(vec,&nLocal);
1568:   if (nLocal != stag->entriesGhost) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMStag local size %D\n",nLocal,stag->entriesGhost);
1569:   switch (dim) {
1570:     case 1:
1571:       VecRestoreArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1572:       break;
1573:     case 2:
1574:       VecRestoreArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1575:       break;
1576:     case 3:
1577:       VecRestoreArray4d(vec,stag->nGhost[2],stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[2],stag->startGhost[1],stag->startGhost[0],0,(PetscScalar*****)array);
1578:       break;
1579:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1580:   }
1581:   return(0);
1582: }

1584: /*@C
1585:   DMStagVecRestoreArrayRead - restore read-only access to a raw array

1587:   Logically Collective

1589:   Input Parameters:
1590: + dm - the DMStag object
1591: - vec - the Vec object

1593:   Output Parameters:
1594: . array - the read-only array

1596:   Level: beginner

1598: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMDAVecRestoreArrayRead(), DMDAVecRestoreArrayDOFRead()
1599: @*/
1600: PetscErrorCode DMStagVecRestoreArrayRead(DM dm,Vec vec,void *array)
1601: {
1602:   PetscErrorCode  ierr;
1603:   DM_Stag * const stag = (DM_Stag*)dm->data;
1604:   PetscInt        dim;
1605:   PetscInt        nLocal;

1610:   DMGetDimension(dm,&dim);
1611:   VecGetLocalSize(vec,&nLocal);
1612:   if (nLocal != stag->entriesGhost) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMStag local size %D\n",nLocal,stag->entriesGhost);
1613:   switch (dim) {
1614:     case 1:
1615:       VecRestoreArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1616:       break;
1617:     case 2:
1618:       VecRestoreArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1619:       break;
1620:     case 3:
1621:       VecRestoreArray4dRead(vec,stag->nGhost[2],stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[2],stag->startGhost[1],stag->startGhost[0],0,(PetscScalar*****)array);
1622:       break;
1623:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1624:   }
1625:   return(0);
1626: }