Actual source code: stagutils.c

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

  7:   Not Collective

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

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

 15:   Level: intermediate

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

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

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

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

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

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

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

 97:   Logically Collective

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

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

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

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

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

117:   Level: intermediate

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

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

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

133:   Logically Collective

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

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

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

143:   Level: intermediate

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

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

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

159:   Not Collective

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

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

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

171:   Notes:
172:   Checks that the coordinates are actually set up so that using the
173:   slots from the first 1d coordinate sub-DM is valid for all the 1D coordinate sub-DMs.

175:   Level: intermediate

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

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

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

227:   Not Collective

229:   Input Parameter:
230: . dm - the DMStag object

232:   Output Parameters:
233: + x     - starting element index in first direction
234: . y     - starting element index in second direction
235: . z     - starting element index in third direction
236: . m     - element width in first direction
237: . n     - element width in second direction
238: . p     - element width in third direction
239: . nExtrax - number of extra partial elements in first direction
240: . nExtray - number of extra partial elements in second direction
241: - nExtraz - number of extra partial elements in third direction

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

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

250:   Level: beginner

252: .seealso: DMSTAG, DMStagGetGhostCorners(), DMDAGetCorners()
253: @*/
254: PetscErrorCode DMStagGetCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *nExtrax,PetscInt *nExtray,PetscInt *nExtraz)
255: {
256:   const DM_Stag * const stag = (DM_Stag*)dm->data;

260:   if (x) *x = stag->start[0];
261:   if (y) *y = stag->start[1];
262:   if (z) *z = stag->start[2];
263:   if (m) *m = stag->n[0];
264:   if (n) *n = stag->n[1];
265:   if (p) *p = stag->n[2];
266:   if (nExtrax) *nExtrax = stag->boundaryType[0] != DM_BOUNDARY_PERIODIC && stag->lastRank[0] ? 1 : 0;
267:   if (nExtray) *nExtray = stag->boundaryType[1] != DM_BOUNDARY_PERIODIC && stag->lastRank[1] ? 1 : 0;
268:   if (nExtraz) *nExtraz = stag->boundaryType[2] != DM_BOUNDARY_PERIODIC && stag->lastRank[2] ? 1 : 0;
269:   return(0);
270: }

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

275:   Not Collective

277:   Input Parameter:
278: . dm - the DMStag object

280:   Output Parameters:
281: + dof0 - the number of points per 0-cell (vertex/node)
282: . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
283: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
284: - dof3 - the number of points per 3-cell (element in 3D)

286:   Level: beginner

288: .seealso: DMSTAG, DMStagGetCorners(), DMStagGetGhostCorners(), DMStagGetGlobalSizes(), DMStagGetStencilWidth(), DMStagGetBoundaryTypes(), DMStagGetLocationDof(), DMDAGetDof()
289: @*/
290: PetscErrorCode DMStagGetDOF(DM dm,PetscInt *dof0,PetscInt *dof1,PetscInt *dof2,PetscInt *dof3)
291: {
292:   const DM_Stag * const stag = (DM_Stag*)dm->data;

296:   if (dof0) *dof0 = stag->dof[0];
297:   if (dof1) *dof1 = stag->dof[1];
298:   if (dof2) *dof2 = stag->dof[2];
299:   if (dof3) *dof3 = stag->dof[3];
300:   return(0);
301: }

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

306:   Not Collective

308:   Input Parameter:
309: . dm - the DMStag object

311:   Output Parameters:
312: + x - the starting element index in the first direction
313: . y - the starting element index in the second direction
314: . z - the starting element index in the third direction
315: . m - the element width in the first direction
316: . n - the element width in the second direction
317: - p - the element width in the third direction

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

322:   Level: beginner

324: .seealso: DMSTAG, DMStagGetCorners(), DMDAGetGhostCorners()
325: @*/
326: PetscErrorCode DMStagGetGhostCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
327: {
328:   const DM_Stag * const stag = (DM_Stag*)dm->data;

332:   if (x) *x = stag->startGhost[0];
333:   if (y) *y = stag->startGhost[1];
334:   if (z) *z = stag->startGhost[2];
335:   if (m) *m = stag->nGhost[0];
336:   if (n) *n = stag->nGhost[1];
337:   if (p) *p = stag->nGhost[2];
338:   return(0);
339: }

341: /*@C
342:   DMStagGetGlobalSizes - get global element counts

344:   Not Collective

346:   Input Parameter:
347: . dm - the DMStag object

349:   Output Parameters:
350: . M,N,P - global element counts in each direction

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

355:   Level: beginner

357: .seealso: DMSTAG, DMStagGetLocalSizes(), DMDAGetInfo()
358: @*/
359: PetscErrorCode DMStagGetGlobalSizes(DM dm,PetscInt* M,PetscInt* N,PetscInt* P)
360: {
361:   const DM_Stag * const stag = (DM_Stag*)dm->data;

365:   if (M) *M = stag->N[0];
366:   if (N) *N = stag->N[1];
367:   if (P) *P = stag->N[2];
368:   return(0);
369: }

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

374:   Not Collective

376:   Input Parameter:
377: . dm - the DMStag object

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

382:   Level: intermediate

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

387: .seealso: DMSTAG, DMStagGetIsLastRank()
388: @*/
389: PetscErrorCode DMStagGetIsFirstRank(DM dm,PetscBool *isFirstRank0,PetscBool *isFirstRank1,PetscBool *isFirstRank2)
390: {
391:   const DM_Stag * const stag = (DM_Stag*)dm->data;

395:   if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
396:   if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
397:   if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
398:   return(0);
399: }

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

404:   Not Collective

406:   Input Parameter:
407: . dm - the DMStag object

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

412:   Level: intermediate

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

418: .seealso: DMSTAG, DMStagGetIsFirstRank()
419: @*/
420: PetscErrorCode DMStagGetIsLastRank(DM dm,PetscBool *isLastRank0,PetscBool *isLastRank1,PetscBool *isLastRank2)
421: {
422:   const DM_Stag * const stag = (DM_Stag*)dm->data;

426:   if (isLastRank0) *isLastRank0 = stag->lastRank[0];
427:   if (isLastRank1) *isLastRank1 = stag->lastRank[1];
428:   if (isLastRank2) *isLastRank2 = stag->lastRank[2];
429:   return(0);
430: }

432: /*@C
433:   DMStagGetLocalSizes - get local elementwise sizes

435:   Not Collective

437:   Input Parameter:
438: . dm - the DMStag object

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

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

447:   Level: beginner

449: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetDOF(), DMStagGetNumRanks(), DMDAGetLocalInfo()
450: @*/
451: PetscErrorCode DMStagGetLocalSizes(DM dm,PetscInt* m,PetscInt* n,PetscInt* p)
452: {
453:   const DM_Stag * const stag = (DM_Stag*)dm->data;

457:   if (m) *m = stag->n[0];
458:   if (n) *n = stag->n[1];
459:   if (p) *p = stag->n[2];
460:   return(0);
461: }

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

466:   Not Collective

468:   Input Parameter:
469: . dm - the DMStag object

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

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

478:   Level: beginner

480: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetLocalSize(), DMStagSetNumRanks(), DMDAGetInfo()
481: @*/
482: PetscErrorCode DMStagGetNumRanks(DM dm,PetscInt *nRanks0,PetscInt *nRanks1,PetscInt *nRanks2)
483: {
484:   const DM_Stag * const stag = (DM_Stag*)dm->data;

488:   if (nRanks0) *nRanks0 = stag->nRanks[0];
489:   if (nRanks1) *nRanks1 = stag->nRanks[1];
490:   if (nRanks2) *nRanks2 = stag->nRanks[2];
491:   return(0);
492: }

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

497:   Not Collective

499:   Input Parameter:
500: . dm - the DMStag object

502:   Output Parameters:
503: . entries - number of rank-native entries in the global representation

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

508:   Level: developer

510: .seealso: DMSTAG, DMStagGetDOF(), DMStagGetEntriesPerElement(), DMCreateLocalVector()
511: @*/
512: PetscErrorCode DMStagGetEntries(DM dm,PetscInt *entries)
513: {
514:   const DM_Stag * const stag = (DM_Stag*)dm->data;

518:   if (entries) *entries = stag->entries;
519:   return(0);
520: }

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

525:   Not Collective

527:   Input Parameter:
528: . dm - the DMStag object

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

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

537:   Level: developer

539: .seealso: DMSTAG, DMStagGetDOF()
540: @*/
541: PetscErrorCode DMStagGetEntriesPerElement(DM dm,PetscInt *entriesPerElement)
542: {
543:   const DM_Stag * const stag = (DM_Stag*)dm->data;

547:   if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
548:   return(0);
549: }

551: /*@C
552:   DMStagGetStencilType - get elementwise ghost/halo stencil type

554:   Not Collective

556:   Input Parameter:
557: . dm - the DMStag object

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

562:   Level: beginner

564: .seealso: DMSTAG, DMStagSetStencilType(), DMStagStencilType, DMDAGetInfo()
565: @*/
566: PetscErrorCode DMStagGetStencilType(DM dm,DMStagStencilType *stencilType)
567: {
568:   DM_Stag * const stag = (DM_Stag*)dm->data;

572:   *stencilType = stag->stencilType;
573:   return(0);
574: }

576: /*@C
577:   DMStagGetStencilWidth - get elementwise stencil width

579:   Not Collective

581:   Input Parameter:
582: . dm - the DMStag object

584:   Output Parameters:
585: . stencilWidth - stencil/halo/ghost width in elements

587:   Level: beginner

589: .seealso: DMSTAG, DMDAGetStencilWidth(), DMDAGetInfo()
590: @*/
591: PetscErrorCode DMStagGetStencilWidth(DM dm,PetscInt *stencilWidth)
592: {
593:   const DM_Stag * const stag = (DM_Stag*)dm->data;

597:   if (stencilWidth) *stencilWidth = stag->stencilWidth;
598:   return(0);
599: }

601: /*@C
602:   DMStagGetOwnershipRanges - get elements per rank in each direction

604:   Not Collective

606:   Input Parameter:
607: .     dm - the DMStag object

609:   Output Parameters:
610: +     lx - ownership along x direction (optional)
611: .     ly - ownership along y direction (optional)
612: -     lz - ownership along z direction (optional)

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

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

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

622:   Level: intermediate

624: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagSetOwnershipRanges(), DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDAGetOwnershipRanges()
625: @*/
626: PetscErrorCode DMStagGetOwnershipRanges(DM dm,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
627: {
628:   const DM_Stag * const stag = (DM_Stag*)dm->data;

632:   if (lx) *lx = stag->l[0];
633:   if (ly) *ly = stag->l[1];
634:   if (lz) *lz = stag->l[2];
635:   return(0);
636: }

638: /*@C
639:   DMStagCreateCompatibleDMStag - create a compatible DMStag with different dof/stratum

641:   Collective

643:   Input Parameters:
644: + dm - the DMStag object
645: - dof0,dof1,dof2,dof3 - number of dof on each stratum in the new DMStag

647:   Output Parameters:
648: . newdm - the new, compatible DMStag

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

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

657:   Level: intermediate

659: .seealso: DMSTAG, DMDACreateCompatibleDMDA(), DMGetCompatibility(), DMStagMigrateVec()
660: @*/
661: PetscErrorCode DMStagCreateCompatibleDMStag(DM dm,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DM *newdm)
662: {
663:   PetscErrorCode  ierr;

667:   DMStagDuplicateWithoutSetup(dm,PetscObjectComm((PetscObject)dm),newdm);
668:   DMStagSetDOF(*newdm,dof0,dof1,dof2,dof3);
669:   DMSetUp(*newdm);
670:   return(0);
671: }

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

676:   Not Collective

678:   Input Parameters:
679: + dm - the DMStag object
680: . loc - location relative to an element
681: - c - component

683:   Output Parameter:
684: . slot - index to use

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

691:   Level: beginner

693: .seealso: DMSTAG, DMStagVecGetArray(), DMStagVecGetArrayRead(), DMStagGetDOF(), DMStagGetEntriesPerElement()
694: @*/
695: PetscErrorCode DMStagGetLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt c,PetscInt *slot)
696: {
697:   DM_Stag * const stag = (DM_Stag*)dm->data;

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

712: /*@C
713:   DMStagMigrateVec - transfer a vector associated with a DMStag to a vector associated with a compatible DMStag

715:   Collective

717:   Input Parameters:
718: + dm - the source DMStag object
719: . vec - the source vector, compatible with dm
720: . dmTo - the compatible destination DMStag object
721: - vecTo - the destination vector, compatible with dmTo

723:   Notes:
724:   Extra dof are ignored, and unfilled dof are zeroed.
725:   Currently only implemented to migrate global vectors to global vectors.

727:   Level: advanced

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

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

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

832: /*@C
833:   DMStagPopulateLocalToGlobalInjective - populate an internal 1-to-1 local-to-global map

835:   Collective

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

846:   Input Parameter:
847: . dm - the DMStag object

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

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

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

860:   Level: developer

862: .seealso: DMSTAG, DMLocalToGlobal(), VecScatter
863: @*/
864: PetscErrorCode DMStagPopulateLocalToGlobalInjective(DM dm)
865: {
866:   PetscErrorCode  ierr;
867:   PetscInt        dim;
868:   DM_Stag * const stag  = (DM_Stag*)dm->data;

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

883: static PetscErrorCode DMStagRestoreProductCoordinateArrays_Private(DM dm,void *arrX,void *arrY,void *arrZ,PetscBool read)
884: {
885:   PetscErrorCode  ierr;
886:   PetscInt        dim,d;
887:   void*           arr[DMSTAG_MAX_DIM];
888:   DM              dmCoord;

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

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

903:     DMProductGetDM(dmCoord,d,&subDM);
904:     DMGetCoordinatesLocal(subDM,&coord1d_local);
905:     if (read) {
906:       DMStagVecRestoreArrayRead(subDM,coord1d_local,arr[d]);
907:     } else {
908:       DMStagVecRestoreArray(subDM,coord1d_local,arr[d]);
909:     }
910:   }
911:   return(0);
912: }

914: /*@C
915:   DMStagRestoreProductCoordinateArrays - restore local array access

917:   Logically Collective

919:   Input Parameter:
920: . dm - the DMStag object

922:   Output Parameters:
923: . arrX,arrY,arrZ - local 1D coordinate arrays

925:   Level: intermediate

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

930: $   DMGetCoordinateDM(dm,&cdm);
931: $   for (d=0; d<3; ++d) {
932: $     DM  subdm;
933: $     Vec coor,coor_local;

935: $     DMProductGetDM(cdm,d,&subdm);
936: $     DMGetCoordinates(subdm,&coor);
937: $     DMGetCoordinatesLocal(subdm,&coor_local);
938: $     DMLocalToGlobal(subdm,coor_local,INSERT_VALUES,coor);
939: $     PetscPrintf(PETSC_COMM_WORLD,"Coordinates dim %D:\n",d);
940: $     VecView(coor,PETSC_VIEWER_STDOUT_WORLD);
941: $   }

943: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
944: @*/
945: PetscErrorCode DMStagRestoreProductCoordinateArrays(DM dm,void *arrX,void *arrY,void *arrZ)
946: {
947:   PetscErrorCode  ierr;

950:   DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_FALSE);
951:   return(0);
952: }

954: /*@C
955:   DMStagRestoreProductCoordinateArraysRead - restore local product array access, read-only

957:   Logically Collective

959:   Input Parameter:
960: . dm - the DMStag object

962:   Output Parameters:
963: . arrX,arrY,arrZ - local 1D coordinate arrays

965:   Level: intermediate

967: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
968: @*/
969: PetscErrorCode DMStagRestoreProductCoordinateArraysRead(DM dm,void *arrX,void *arrY,void *arrZ)
970: {
971:   PetscErrorCode  ierr;

974:   DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_TRUE);
975:   return(0);
976: }

978: /*@C
979:   DMStagSetBoundaryTypes - set DMStag boundary types

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

983:   Input Parameters:
984: + dm - the DMStag object
985: - boundaryType0,boundaryType1,boundaryType2 - boundary types in each direction

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

990:   Level: advanced

992: .seealso: DMSTAG, DMBoundaryType, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDASetBoundaryType()
993: @*/
994: PetscErrorCode DMStagSetBoundaryTypes(DM dm,DMBoundaryType boundaryType0,DMBoundaryType boundaryType1,DMBoundaryType boundaryType2)
995: {
996:   PetscErrorCode  ierr;
997:   DM_Stag * const stag  = (DM_Stag*)dm->data;
998:   PetscInt        dim;

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

1013: /*@C
1014:   DMStagSetCoordinateDMType - set DM type to store coordinates

1016:   Logically Collective; dmtype must contain common value

1018:   Input Parameters:
1019: + dm - the DMStag object
1020: - dmtype - DMtype for coordinates, either DMSTAG or DMPRODUCT

1022:   Level: advanced

1024: .seealso: DMSTAG, DMPRODUCT, DMGetCoordinateDM(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMType
1025: @*/
1026: PetscErrorCode DMStagSetCoordinateDMType(DM dm,DMType dmtype)
1027: {
1028:   PetscErrorCode  ierr;
1029:   DM_Stag * const stag = (DM_Stag*)dm->data;

1033:   PetscFree(stag->coordinateDMType);
1034:   PetscStrallocpy(dmtype,(char**)&stag->coordinateDMType);
1035:   return(0);
1036: }

1038: /*@C
1039:   DMStagSetDOF - set dof/stratum

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

1043:   Input Parameters:
1044: + dm - the DMStag object
1045: - dof0,dof1,dof2,dof3 - dof per stratum

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

1050:   Level: advanced

1052: .seealso: DMSTAG, DMDASetDof()
1053: @*/
1054: PetscErrorCode DMStagSetDOF(DM dm,PetscInt dof0, PetscInt dof1,PetscInt dof2,PetscInt dof3)
1055: {
1056:   PetscErrorCode  ierr;
1057:   DM_Stag * const stag = (DM_Stag*)dm->data;
1058:   PetscInt        dim;

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

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

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

1084:   Input Parameters:
1085: + dm - the DMStag object
1086: - nRanks0,nRanks1,nRanks2 - number of ranks in each direction

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

1091:   Level: developer

1093: .seealso: DMSTAG, DMDASetNumProcs()
1094: @*/
1095: PetscErrorCode DMStagSetNumRanks(DM dm,PetscInt nRanks0,PetscInt nRanks1,PetscInt nRanks2)
1096: {
1097:   PetscErrorCode  ierr;
1098:   DM_Stag * const stag = (DM_Stag*)dm->data;
1099:   PetscInt        dim;

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

1117: /*@C
1118:   DMStagSetStencilType - set elementwise ghost/halo stencil type

1120:   Logically Collective; stencilType must contain common value

1122:   Input Parameters:
1123: + dm - the DMStag object
1124: - stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE

1126:   Level: beginner

1128: .seealso: DMSTAG, DMStagGetStencilType(), DMStagStencilType, DMDASetStencilType()
1129: @*/
1130: PetscErrorCode DMStagSetStencilType(DM dm,DMStagStencilType stencilType)
1131: {
1132:   DM_Stag * const stag = (DM_Stag*)dm->data;

1137:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1138:   stag->stencilType = stencilType;
1139:   return(0);
1140: }

1142: /*@C
1143:   DMStagSetStencilWidth - set elementwise stencil width

1145:   Logically Collective; stencilWidth must contain common value

1147:   Input Parameters:
1148: + dm - the DMStag object
1149: - stencilWidth - stencil/halo/ghost width in elements

1151:   Level: beginner

1153: .seealso: DMSTAG, DMDASetStencilWidth()
1154: @*/
1155: PetscErrorCode DMStagSetStencilWidth(DM dm,PetscInt stencilWidth)
1156: {
1157:   DM_Stag * const stag = (DM_Stag*)dm->data;

1162:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1163:   if (stencilWidth < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Stencil width must be non-negative");
1164:   stag->stencilWidth = stencilWidth;
1165:   return(0);
1166: }

1168: /*@C
1169:   DMStagSetGlobalSizes - set global element counts in each direction

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

1173:   Input Parameters:
1174: + dm - the DMStag object
1175: - N0,N1,N2 - global elementwise sizes

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

1180:   Level: advanced

1182: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMDASetSizes()
1183: @*/
1184: PetscErrorCode DMStagSetGlobalSizes(DM dm,PetscInt N0,PetscInt N1,PetscInt N2)
1185: {
1186:   PetscErrorCode  ierr;
1187:   DM_Stag * const stag = (DM_Stag*)dm->data;
1188:   PetscInt        dim;

1192:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1193:   DMGetDimension(dm,&dim);
1194:   if (N0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in X direction must be positive");
1195:   if (dim > 1 && N1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Y direction must be positive");
1196:   if (dim > 2 && N2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Z direction must be positive");
1197:   if (N0) stag->N[0] = N0;
1198:   if (N1) stag->N[1] = N1;
1199:   if (N2) stag->N[2] = N2;
1200:   return(0);
1201: }

1203: /*@C
1204:   DMStagSetOwnershipRanges - set elements per rank in each direction

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

1208:   Input Parameters:
1209: + dm - the DMStag object
1210: - lx,ly,lz - element counts for each rank in each direction

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

1215:   Level: developer

1217: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagGetOwnershipRanges(), DMDASetOwnershipRanges()
1218: @*/
1219: PetscErrorCode DMStagSetOwnershipRanges(DM dm,PetscInt const *lx,PetscInt const *ly,PetscInt const *lz)
1220: {
1221:   PetscErrorCode  ierr;
1222:   DM_Stag * const stag = (DM_Stag*)dm->data;
1223:   const PetscInt  *lin[3];
1224:   PetscInt        d,dim;

1228:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1229:   lin[0] = lx; lin[1] = ly; lin[2] = lz;
1230:   DMGetDimension(dm,&dim);
1231:   for (d=0; d<dim; ++d) {
1232:     if (lin[d]) {
1233:       if (stag->nRanks[d] < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of ranks");
1234:       if (!stag->l[d]) {
1235:         PetscMalloc1(stag->nRanks[d], &stag->l[d]);
1236:       }
1237:       PetscArraycpy(stag->l[d], lin[d], stag->nRanks[d]);
1238:     }
1239:   }
1240:   return(0);
1241: }

1243: /*@C
1244:   DMStagSetUniformCoordinates - set DMStag coordinates to be a uniform grid

1246:   Collective

1248:   Input Parameters:
1249: + dm - the DMStag object
1250: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

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

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

1264:   Level: advanced

1266: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType(), DMGetCoordinateDM(), DMGetCoordinates(), DMDASetUniformCoordinates(), DMBoundaryType
1267: @*/
1268: PetscErrorCode DMStagSetUniformCoordinates(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1269: {
1270:   PetscErrorCode  ierr;
1271:   DM_Stag * const stag = (DM_Stag*)dm->data;
1272:   PetscBool       flg_stag,flg_product;

1276:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1277:   if (!stag->coordinateDMType) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"You must first call DMStagSetCoordinateDMType()");
1278:   PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg_stag);
1279:   PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg_product);
1280:   if (flg_stag) {
1281:     DMStagSetUniformCoordinatesExplicit(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1282:   } else if (flg_product) {
1283:     DMStagSetUniformCoordinatesProduct(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1284:   } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported DM Type %s",stag->coordinateDMType);
1285:   return(0);
1286: }

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

1291:   Collective

1293:   Input Parameters:
1294: + dm - the DMStag object
1295: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

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

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

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

1306:   Level: beginner

1308: .seealso: DMSTAG, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType()
1309: @*/
1310: PetscErrorCode DMStagSetUniformCoordinatesExplicit(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1311: {
1312:   PetscErrorCode  ierr;
1313:   DM_Stag * const stag = (DM_Stag*)dm->data;
1314:   PetscInt        dim;
1315:   PetscBool       flg;

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

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

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

1338:   Collective

1340:   Input Parameters:
1341: + dm - the DMStag object
1342: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

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

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

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

1354:   Level: intermediate

1356: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetCoordinateDMType()
1357: @*/
1358: PetscErrorCode DMStagSetUniformCoordinatesProduct(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1359: {
1360:   PetscErrorCode  ierr;
1361:   DM_Stag * const stag = (DM_Stag*)dm->data;
1362:   DM              dmc;
1363:   PetscInt        dim,d,dof0,dof1;
1364:   PetscBool       flg;

1368:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1369:   PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg);
1370:   if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1371:   DMStagSetCoordinateDMType(dm,DMPRODUCT);
1372:   DMGetCoordinateDM(dm,&dmc);
1373:   DMGetDimension(dm,&dim);

1375:   /* Create 1D sub-DMs, living on subcommunicators.
1376:      Always include both vertex and element dof, regardless of the active strata of the DMStag */
1377:   dof0 = 1;
1378:   dof1 = 1;

1380:   for (d=0; d<dim; ++d) {
1381:     DM                subdm;
1382:     MPI_Comm          subcomm;
1383:     PetscMPIInt       color;
1384:     const PetscMPIInt key = 0; /* let existing rank break ties */

1386:     /* Choose colors based on position in the plane orthogonal to this dim, and split */
1387:     switch (d) {
1388:       case 0: color = (dim > 1 ? stag->rank[1] : 0)  + (dim > 2 ? stag->nRanks[1]*stag->rank[2] : 0); break;
1389:       case 1: color =            stag->rank[0]       + (dim > 2 ? stag->nRanks[0]*stag->rank[2] : 0); break;
1390:       case 2: color =            stag->rank[0]       +            stag->nRanks[0]*stag->rank[1]     ; break;
1391:       default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1392:     }
1393:     MPI_Comm_split(PetscObjectComm((PetscObject)dm),color,key,&subcomm);

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

1418: /*@C
1419:   DMStagVecGetArray - get access to local array

1421:   Logically Collective

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

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

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

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

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

1441:   Input Parameters:
1442: + dm - the DMStag object
1443: - vec - the Vec object

1445:   Output Parameters:
1446: . array - the array

1448:   Notes:
1449:   DMStagVecRestoreArray() must be called, once finished with the array

1451:   Level: beginner

1453: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArray(), DMDAVecGetArrayDOF()
1454: @*/
1455: PetscErrorCode DMStagVecGetArray(DM dm,Vec vec,void *array)
1456: {
1457:   PetscErrorCode  ierr;
1458:   DM_Stag * const stag = (DM_Stag*)dm->data;
1459:   PetscInt        dim;
1460:   PetscInt        nLocal;

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

1483: /*@C
1484:   DMStagVecGetArrayRead - get read-only access to a local array

1486:   Logically Collective

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

1490:   Input Parameters:
1491: + dm - the DMStag object
1492: - vec - the Vec object

1494:   Output Parameters:
1495: . array - the read-only array

1497:   Notes:
1498:   DMStagVecRestoreArrayRead() must be called, once finished with the array

1500:   Level: beginner

1502: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArrayRead(), DMDAVecGetArrayDOFRead()
1503: @*/
1504: PetscErrorCode DMStagVecGetArrayRead(DM dm,Vec vec,void *array)
1505: {
1506:   PetscErrorCode  ierr;
1507:   DM_Stag * const stag = (DM_Stag*)dm->data;
1508:   PetscInt        dim;
1509:   PetscInt        nLocal;

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

1532: /*@C
1533:   DMStagVecRestoreArray - restore access to a raw array

1535:   Logically Collective

1537:   Input Parameters:
1538: + dm - the DMStag object
1539: - vec - the Vec object

1541:   Output Parameters:
1542: . array - the array

1544:   Level: beginner

1546: .seealso: DMSTAG, DMStagVecGetArray(), DMDAVecRestoreArray(), DMDAVecRestoreArrayDOF()
1547: @*/
1548: PetscErrorCode DMStagVecRestoreArray(DM dm,Vec vec,void *array)
1549: {
1550:   PetscErrorCode  ierr;
1551:   DM_Stag * const stag = (DM_Stag*)dm->data;
1552:   PetscInt        dim;
1553:   PetscInt        nLocal;

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

1576: /*@C
1577:   DMStagVecRestoreArrayRead - restore read-only access to a raw array

1579:   Logically Collective

1581:   Input Parameters:
1582: + dm - the DMStag object
1583: - vec - the Vec object

1585:   Output Parameters:
1586: . array - the read-only array

1588:   Level: beginner

1590: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMDAVecRestoreArrayRead(), DMDAVecRestoreArrayDOFRead()
1591: @*/
1592: PetscErrorCode DMStagVecRestoreArrayRead(DM dm,Vec vec,void *array)
1593: {
1594:   PetscErrorCode  ierr;
1595:   DM_Stag * const stag = (DM_Stag*)dm->data;
1596:   PetscInt        dim;
1597:   PetscInt        nLocal;

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