Actual source code: stagutils.c

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

  7:   Not Collective

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

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

 15:   Level: intermediate

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

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

 34: 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,y,z - starting element indices in each direction
234: . m,n,p - element widths in each direction
235: - nExtrax,nExtray,nExtraz - number of extra partial elements in each direction.

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

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

244:   Level: beginner

246: .seealso: DMSTAG, DMStagGetGhostCorners(), DMDAGetCorners()
247: @*/
248: PetscErrorCode DMStagGetCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *nExtrax,PetscInt *nExtray,PetscInt *nExtraz)
249: {
250:   const DM_Stag * const stag = (DM_Stag*)dm->data;

254:   if (x) *x = stag->start[0];
255:   if (y) *y = stag->start[1];
256:   if (z) *z = stag->start[2];
257:   if (m) *m = stag->n[0];
258:   if (n) *n = stag->n[1];
259:   if (p) *p = stag->n[2];
260:   if (nExtrax) *nExtrax = stag->boundaryType[0] != DM_BOUNDARY_PERIODIC && stag->lastRank[0] ? 1 : 0;
261:   if (nExtray) *nExtray = stag->boundaryType[1] != DM_BOUNDARY_PERIODIC && stag->lastRank[1] ? 1 : 0;
262:   if (nExtraz) *nExtraz = stag->boundaryType[2] != DM_BOUNDARY_PERIODIC && stag->lastRank[2] ? 1 : 0;
263:   return(0);
264: }

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

269:   Not Collective

271:   Input Parameter:
272: . dm - the DMStag object

274:   Output Parameters:
275: + dof0 - the number of points per 0-cell (vertex/node)
276: . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
277: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
278: - dof3 - the number of points per 3-cell (element in 3D)

280:   Level: beginner

282: .seealso: DMSTAG, DMStagGetCorners(), DMStagGetGhostCorners(), DMStagGetGlobalSizes(), DMStagGetStencilWidth(), DMStagGetBoundaryTypes(), DMStagGetLocationDof(), DMDAGetDof()
283: @*/
284: PetscErrorCode DMStagGetDOF(DM dm,PetscInt *dof0,PetscInt *dof1,PetscInt *dof2,PetscInt *dof3)
285: {
286:   const DM_Stag * const stag = (DM_Stag*)dm->data;

290:   if (dof0) *dof0 = stag->dof[0];
291:   if (dof1) *dof1 = stag->dof[1];
292:   if (dof2) *dof2 = stag->dof[2];
293:   if (dof3) *dof3 = stag->dof[3];
294:   return(0);
295: }

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

300:   Not Collective

302:   Input Argument:
303: . dm - the DMStag object

305:   Output Arguments:
306: + x,y,z - starting element indices in each direction
307: - m,n,p - element widths in each direction

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

312:   Level: beginner

314: .seealso: DMSTAG, DMStagGetCorners(), DMDAGetGhostCorners()
315: @*/
316: PetscErrorCode DMStagGetGhostCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
317: {
318:   const DM_Stag * const stag = (DM_Stag*)dm->data;

322:   if (x) *x = stag->startGhost[0];
323:   if (y) *y = stag->startGhost[1];
324:   if (z) *z = stag->startGhost[2];
325:   if (m) *m = stag->nGhost[0];
326:   if (n) *n = stag->nGhost[1];
327:   if (p) *p = stag->nGhost[2];
328:   return(0);
329: }

331: /*@C
332:   DMStagGetGlobalSizes - get global element counts

334:   Not Collective

336:   Input Parameter:
337: . dm - the DMStag object

339:   Output Parameters:
340: . M,N,P - global element counts in each direction

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

345:   Level: beginner

347: .seealso: DMSTAG, DMStagGetLocalSizes(), DMDAGetInfo()
348: @*/
349: PetscErrorCode DMStagGetGlobalSizes(DM dm,PetscInt* M,PetscInt* N,PetscInt* P)
350: {
351:   const DM_Stag * const stag = (DM_Stag*)dm->data;

355:   if (M) *M = stag->N[0];
356:   if (N) *N = stag->N[1];
357:   if (P) *P = stag->N[2];
358:   return(0);
359: }

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

364:   Not Collective

366:   Input Parameter:
367: . dm - the DMStag object

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

372:   Level: intermediate

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

377: .seealso: DMSTAG, DMStagGetIsLastRank()
378: @*/
379: PetscErrorCode DMStagGetIsFirstRank(DM dm,PetscBool *isFirstRank0,PetscBool *isFirstRank1,PetscBool *isFirstRank2)
380: {
381:   const DM_Stag * const stag = (DM_Stag*)dm->data;

385:   if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
386:   if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
387:   if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
388:   return(0);
389: }

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

394:   Not Collective

396:   Input Parameter:
397: . dm - the DMStag object

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

402:   Level: intermediate

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

408: .seealso: DMSTAG, DMStagGetIsFirstRank()
409: @*/
410: PetscErrorCode DMStagGetIsLastRank(DM dm,PetscBool *isLastRank0,PetscBool *isLastRank1,PetscBool *isLastRank2)
411: {
412:   const DM_Stag * const stag = (DM_Stag*)dm->data;

416:   if (isLastRank0) *isLastRank0 = stag->lastRank[0];
417:   if (isLastRank1) *isLastRank1 = stag->lastRank[1];
418:   if (isLastRank2) *isLastRank2 = stag->lastRank[2];
419:   return(0);
420: }

422: /*@C
423:   DMStagGetLocalSizes - get local elementwise sizes

425:   Not Collective

427:   Input Parameter:
428: . dm - the DMStag object

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

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

437:   Level: beginner

439: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetDOF(), DMStagGetNumRanks(), DMDAGetLocalInfo()
440: @*/
441: PetscErrorCode DMStagGetLocalSizes(DM dm,PetscInt* m,PetscInt* n,PetscInt* p)
442: {
443:   const DM_Stag * const stag = (DM_Stag*)dm->data;

447:   if (m) *m = stag->n[0];
448:   if (n) *n = stag->n[1];
449:   if (p) *p = stag->n[2];
450:   return(0);
451: }

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

456:   Not Collective

458:   Input Parameter:
459: . dm - the DMStag object

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

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

468:   Level: beginner

470: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetLocalSize(), DMStagSetNumRanks(), DMDAGetInfo()
471: @*/
472: PetscErrorCode DMStagGetNumRanks(DM dm,PetscInt *nRanks0,PetscInt *nRanks1,PetscInt *nRanks2)
473: {
474:   const DM_Stag * const stag = (DM_Stag*)dm->data;

478:   if (nRanks0) *nRanks0 = stag->nRanks[0];
479:   if (nRanks1) *nRanks1 = stag->nRanks[1];
480:   if (nRanks2) *nRanks2 = stag->nRanks[2];
481:   return(0);
482: }

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

487:   Not Collective

489:   Input Parameter:
490: . dm - the DMStag object

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

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

499:   Level: developer

501: .seealso: DMSTAG, DMStagGetDOF()
502: @*/
503: PetscErrorCode DMStagGetEntriesPerElement(DM dm,PetscInt *entriesPerElement)
504: {
505:   const DM_Stag * const stag = (DM_Stag*)dm->data;

509:   if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
510:   return(0);
511: }

513: /*@C
514:   DMStagGetStencilType - get elementwise ghost/halo stencil type

516:   Not Collective

518:   Input Parameter:
519: . dm - the DMStag object

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

524:   Level: beginner

526: .seealso: DMSTAG, DMStagSetStencilType(), DMStagStencilType, DMDAGetInfo()
527: @*/
528: PetscErrorCode DMStagGetStencilType(DM dm,DMStagStencilType *stencilType)
529: {
530:   DM_Stag * const stag = (DM_Stag*)dm->data;

534:   *stencilType = stag->stencilType;
535:   return(0);
536: }

538: /*@C
539:   DMStagGetStencilWidth - get elementwise stencil width

541:   Not Collective

543:   Input Parameter:
544: . dm - the DMStag object

546:   Output Parameters:
547: . stencilWidth - stencil/halo/ghost width in elements

549:   Level: beginner

551: .seealso: DMSTAG, DMDAGetStencilWidth(), DMDAGetInfo()
552: @*/
553: PetscErrorCode DMStagGetStencilWidth(DM dm,PetscInt *stencilWidth)
554: {
555:   const DM_Stag * const stag = (DM_Stag*)dm->data;

559:   if (stencilWidth) *stencilWidth = stag->stencilWidth;
560:   return(0);
561: }

563: /*@C
564:   DMStagGetOwnershipRanges - get elements per rank in each direction

566:   Not Collective

568:   Input Parameter:
569: .     dm - the DMStag object

571:   Output Parameters:
572: +     lx - ownership along x direction (optional)
573: .     ly - ownership along y direction (optional)
574: -     lz - ownership along z direction (optional)

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

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

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

584:   Level: intermediate

586: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagSetOwnershipRanges(), DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDAGetOwnershipRanges()
587: @*/
588: PetscErrorCode DMStagGetOwnershipRanges(DM dm,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
589: {
590:   const DM_Stag * const stag = (DM_Stag*)dm->data;

594:   if (lx) *lx = stag->l[0];
595:   if (ly) *ly = stag->l[1];
596:   if (lz) *lz = stag->l[2];
597:   return(0);
598: }

600: /*@C
601:   DMStagCreateCompatibleDMStag - create a compatible DMStag with different dof/stratum

603:   Collective

605:   Input Parameters:
606: + dm - the DMStag object
607: - dof0,dof1,dof2,dof3 - number of dof on each stratum in the new DMStag

609:   Output Parameters:
610: . newdm - the new, compatible DMStag

612:   Notes:
613:   Dof supplied for strata too big for the dimension are ignored; these may be set to 0.
614:   In contrast to DMDACreateCompatibleDMDA(), coordinates are not reused.

616:   Level: intermediate

618: .seealso: DMSTAG, DMDACreateCompatibleDMDA(), DMGetCompatibility(), DMStagMigrateVec()
619: @*/
620: PetscErrorCode DMStagCreateCompatibleDMStag(DM dm,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DM *newdm)
621: {
622:   PetscErrorCode  ierr;

626:   DMStagDuplicateWithoutSetup(dm,PetscObjectComm((PetscObject)dm),newdm);
627:   DMStagSetDOF(*newdm,dof0,dof1,dof2,dof3);
628:   DMSetUp(*newdm);
629:   return(0);
630: }

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

635:   Not Collective

637:   Input Parameters:
638: + dm - the DMStag object
639: . loc - location relative to an element
640: - c - component

642:   Output Parameter:
643: . slot - index to use

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

650:   Level: beginner

652: .seealso: DMSTAG, DMStagVecGetArray(), DMStagVecGetArrayRead(), DMStagGetDOF(), DMStagGetEntriesPerElement()
653: @*/
654: PetscErrorCode DMStagGetLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt c,PetscInt *slot)
655: {
656:   DM_Stag * const stag = (DM_Stag*)dm->data;

660: #if defined(PETSC_USE_DEBUG)
661:   {
663:     PetscInt       dof;
664:     DMStagGetLocationDOF(dm,loc,&dof);
665:     if (dof < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Location %s has no dof attached",DMStagStencilLocations[loc]);
666:     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);
667:   }
668: #endif
669:   *slot = stag->locationOffsets[loc] + c;
670:   return(0);
671: }

673: /*@C
674:   DMStagMigrateVec - transfer a vector associated with a DMStag to a vector associated with a compatible DMStag

676:   Collective

678:   Input Parameters:
679: + dm - the source DMStag object
680: . vec - the source vector, compatible with dm
681: . dmTo - the compatible destination DMStag object
682: - vecTo - the destination vector, compatible with dmTo

684:   Notes:
685:   Extra dof are ignored, and unfilled dof are zeroed.
686:   Currently only implemented to migrate global vectors to global vectors.

688:   Level: advanced

690: .seealso: DMSTAG, DMStagCreateCompatibleDMStag(), DMGetCompatibility(), DMStagVecSplitToDMDA()
691: @*/
692: PetscErrorCode DMStagMigrateVec(DM dm,Vec vec,DM dmTo,Vec vecTo)
693: {
694:   PetscErrorCode    ierr;
695:   DM_Stag * const   stag = (DM_Stag*)dm->data;
696:   DM_Stag * const   stagTo = (DM_Stag*)dmTo->data;
697:   PetscInt          nLocalTo,nLocal,dim,i,j,k;
698:   PetscInt          start[DMSTAG_MAX_DIM],startGhost[DMSTAG_MAX_DIM],n[DMSTAG_MAX_DIM],nExtra[DMSTAG_MAX_DIM],offset[DMSTAG_MAX_DIM];
699:   Vec               vecToLocal,vecLocal;
700:   PetscBool         compatible,compatibleSet;
701:   const PetscScalar *arr;
702:   PetscScalar       *arrTo;
703:   const PetscInt    epe   = stag->entriesPerElement;
704:   const PetscInt    epeTo = stagTo->entriesPerElement;

711:   DMGetCompatibility(dm,dmTo,&compatible,&compatibleSet);
712:   if (!compatibleSet || !compatible) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"DMStag objects must be shown to be compatible");
713:   DMGetDimension(dm,&dim);
714:   VecGetLocalSize(vec,&nLocal);
715:   VecGetLocalSize(vecTo,&nLocalTo);
716:   if (nLocal != stag->entries|| nLocalTo !=stagTo->entries) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Vector migration only implemented for global vector to global vector.");
717:   DMStagGetCorners(dm,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
718:   DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],NULL,NULL,NULL);
719:   for (i=0; i<DMSTAG_MAX_DIM; ++i) offset[i] = start[i]-startGhost[i];

721:   /* Proceed by transferring to a local vector, copying, and transferring back to a global vector */
722:   DMGetLocalVector(dm,&vecLocal);
723:   DMGetLocalVector(dmTo,&vecToLocal);
724:   DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal);
725:   DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal);
726:   VecGetArrayRead(vecLocal,&arr);
727:   VecGetArray(vecToLocal,&arrTo);
728:   /* Note that some superfluous copying of entries on partial dummy elements is done */
729:   if (dim == 1) {
730:     for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
731:       PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
732:       PetscInt si;
733:       for (si=0; si<2; ++si) {
734:         b   += stag->dof[si];
735:         bTo += stagTo->dof[si];
736:         for (; d < b && dTo < bTo; ++d,++dTo) arrTo[i*epeTo + dTo] = arr[i*epe + d];
737:         for (; dTo < bTo         ;     ++dTo) arrTo[i*epeTo + dTo] = 0.0;
738:         d = b;
739:       }
740:     }
741:   } else if (dim == 2) {
742:     const PetscInt epr   = stag->nGhost[0] * epe;
743:     const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
744:     for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
745:       for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
746:         const PetscInt base   = j*epr   + i*epe;
747:         const PetscInt baseTo = j*eprTo + i*epeTo;
748:         PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
749:         const PetscInt s[4] = {0,1,1,2}; /* Dimensions of points, in order */
750:         PetscInt si;
751:         for (si=0; si<4; ++si) {
752:             b   += stag->dof[s[si]];
753:             bTo += stagTo->dof[s[si]];
754:             for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
755:             for (;          dTo < bTo;     ++dTo) arrTo[baseTo + dTo] = 0.0;
756:             d = b;
757:         }
758:       }
759:     }
760:   } else if (dim == 3) {
761:     const PetscInt epr   = stag->nGhost[0]   * epe;
762:     const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
763:     const PetscInt epl   = stag->nGhost[1]   * epr;
764:     const PetscInt eplTo = stagTo->nGhost[1] * eprTo;
765:     for (k=offset[2]; k<offset[2] + n[2] + nExtra[2]; ++k) {
766:       for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
767:         for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
768:           PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
769:           const PetscInt base   = k*epl   + j*epr   + i*epe;
770:           const PetscInt baseTo = k*eplTo + j*eprTo + i*epeTo;
771:           const PetscInt s[8] = {0,1,1,2,1,2,2,3}; /* dimensions of points, in order */
772:           PetscInt is;
773:           for (is=0; is<8; ++is) {
774:             b   += stag->dof[s[is]];
775:             bTo += stagTo->dof[s[is]];
776:             for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
777:             for (;          dTo < bTo;     ++dTo) arrTo[baseTo + dTo] = 0.0;
778:             d = b;
779:           }
780:         }
781:       }
782:     }
783:   } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
784:   VecRestoreArrayRead(vecLocal,&arr);
785:   VecRestoreArray(vecToLocal,&arrTo);
786:   DMRestoreLocalVector(dm,&vecLocal);
787:   DMLocalToGlobalBegin(dmTo,vecToLocal,INSERT_VALUES,vecTo);
788:   DMLocalToGlobalEnd(dmTo,vecToLocal,INSERT_VALUES,vecTo);
789:   DMRestoreLocalVector(dmTo,&vecToLocal);
790:   return(0);
791: }

793: /*@C
794:   DMStagPopulateLocalToGlobalInjective - populate an internal 1-to-1 local-to-global map

796:   Collective

798:   Creates an internal object which explicitly maps a single local degree of
799:   freedom to each global degree of freedom. This is used, if populated,
800:   instead of SCATTER_REVERSE_LOCAL with the (1-to-many, in general)
801:   global-to-local map, when DMLocalToGlobal() is called with INSERT_VALUES.
802:   This allows usage, for example, even in the periodic, 1-rank case, where
803:   the inverse of the global-to-local map, even when restricted to on-rank
804:   communication, is non-injective. This is at the cost of storing an additional
805:   VecScatter object inside each DMStag object.

807:   Input Parameter:
808: . dm - the DMStag object

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

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

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

821:   Level: developer

823: .seealso: DMSTAG, DMLocalToGlobal(), VecScatter
824: @*/
825: PetscErrorCode DMStagPopulateLocalToGlobalInjective(DM dm)
826: {
827:   PetscErrorCode  ierr;
828:   PetscInt        dim;
829:   DM_Stag * const stag  = (DM_Stag*)dm->data;

833:   if (stag->ltog_injective) return(0); /* Don't re-populate */
834:   DMGetDimension(dm,&dim);
835:   switch (dim) {
836:     case 1: DMStagPopulateLocalToGlobalInjective_1d(dm); break;
837:     case 2: DMStagPopulateLocalToGlobalInjective_2d(dm); break;
838:     case 3: DMStagPopulateLocalToGlobalInjective_3d(dm); break;
839:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
840:   }
841:   return(0);
842: }

844: static PetscErrorCode DMStagRestoreProductCoordinateArrays_Private(DM dm,void *arrX,void *arrY,void *arrZ,PetscBool read)
845: {
846:   PetscErrorCode  ierr;
847:   PetscInt        dim,d;
848:   void*           arr[DMSTAG_MAX_DIM];
849:   DM              dmCoord;

853:   DMGetDimension(dm,&dim);
854:   if (dim > DMSTAG_MAX_DIM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for %D dimensions",dim);
855:   arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
856:   DMGetCoordinateDM(dm,&dmCoord);
857:   for (d=0; d<dim; ++d) {
858:     DM  subDM;
859:     Vec coord1d_local;

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

864:     DMProductGetDM(dmCoord,d,&subDM);
865:     DMGetCoordinatesLocal(subDM,&coord1d_local);
866:     if (read) {
867:       DMStagVecRestoreArrayRead(subDM,coord1d_local,arr[d]);
868:     } else {
869:       DMStagVecRestoreArray(subDM,coord1d_local,arr[d]);
870:     }
871:   }
872:   return(0);
873: }

875: /*@C
876:   DMStagRestoreProductCoordinateArrays - restore local array access

878:   Logically Collective

880:   Input Parameter:
881: . dm - the DMStag object

883:   Output Parameters:
884: . arrX,arrY,arrZ - local 1D coordinate arrays

886:   Level: intermediate

888:   Notes:
889:   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:

891: $   DMGetCoordinateDM(dm,&cdm);
892: $   for (d=0; d<3; ++d) {
893: $     DM  subdm;
894: $     Vec coor,coor_local;

896: $     DMProductGetDM(cdm,d,&subdm);
897: $     DMGetCoordinates(subdm,&coor);
898: $     DMGetCoordinatesLocal(subdm,&coor_local);
899: $     DMLocalToGlobal(subdm,coor_local,INSERT_VALUES,coor);
900: $     PetscPrintf(PETSC_COMM_WORLD,"Coordinates dim %D:\n",d);
901: $     VecView(coor,PETSC_VIEWER_STDOUT_WORLD);
902: $   }

904: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
905: @*/
906: PetscErrorCode DMStagRestoreProductCoordinateArrays(DM dm,void *arrX,void *arrY,void *arrZ)
907: {
908:   PetscErrorCode  ierr;

911:   DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_FALSE);
912:   return(0);
913: }

915: /*@C
916:   DMStagRestoreProductCoordinateArraysRead - restore local product array access, read-only

918:   Logically Collective

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

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

926:   Level: intermediate

928: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
929: @*/
930: PetscErrorCode DMStagRestoreProductCoordinateArraysRead(DM dm,void *arrX,void *arrY,void *arrZ)
931: {
932:   PetscErrorCode  ierr;

935:   DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_TRUE);
936:   return(0);
937: }

939: /*@C
940:   DMStagSetBoundaryTypes - set DMStag boundary types

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

944:   Input Parameters:
945: + dm - the DMStag object
946: - boundaryType0,boundaryType1,boundaryType2 - boundary types in each direction

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

951:   Level: advanced

953: .seealso: DMSTAG, DMBoundaryType, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDASetBoundaryType()
954: @*/
955: PetscErrorCode DMStagSetBoundaryTypes(DM dm,DMBoundaryType boundaryType0,DMBoundaryType boundaryType1,DMBoundaryType boundaryType2)
956: {
957:   PetscErrorCode  ierr;
958:   DM_Stag * const stag  = (DM_Stag*)dm->data;
959:   PetscInt        dim;

962:   DMGetDimension(dm,&dim);
967:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
968:                stag->boundaryType[0] = boundaryType0;
969:   if (dim > 1) stag->boundaryType[1] = boundaryType1;
970:   if (dim > 2) stag->boundaryType[2] = boundaryType2;
971:   return(0);
972: }

974: /*@C
975:   DMStagSetCoordinateDMType - set DM type to store coordinates

977:   Logically Collective; dmtype must contain common value

979:   Input Parameters:
980: + dm - the DMStag object
981: - dmtype - DMtype for coordinates, either DMSTAG or DMPRODUCT

983:   Level: advanced

985: .seealso: DMSTAG, DMPRODUCT, DMGetCoordinateDM(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMType
986: @*/
987: PetscErrorCode DMStagSetCoordinateDMType(DM dm,DMType dmtype)
988: {
989:   PetscErrorCode  ierr;
990:   DM_Stag * const stag = (DM_Stag*)dm->data;

994:   PetscFree(stag->coordinateDMType);
995:   PetscStrallocpy(dmtype,(char**)&stag->coordinateDMType);
996:   return(0);
997: }

999: /*@C
1000:   DMStagSetDOF - set dof/stratum

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

1004:   Input Parameters:
1005: + dm - the DMStag object
1006: - dof0,dof1,dof2,dof3 - dof per stratum

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

1011:   Level: advanced

1013: .seealso: DMSTAG, DMDASetDof()
1014: @*/
1015: PetscErrorCode DMStagSetDOF(DM dm,PetscInt dof0, PetscInt dof1,PetscInt dof2,PetscInt dof3)
1016: {
1017:   PetscErrorCode  ierr;
1018:   DM_Stag * const stag = (DM_Stag*)dm->data;
1019:   PetscInt        dim;

1027:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1028:   DMGetDimension(dm,&dim);
1029:   if (dof0 < 0)            SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof0 cannot be negative");
1030:   if (dof1 < 0)            SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof1 cannot be negative");
1031:   if (dim > 1 && dof2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof2 cannot be negative");
1032:   if (dim > 2 && dof3 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof3 cannot be negative");
1033:                stag->dof[0] = dof0;
1034:                stag->dof[1] = dof1;
1035:   if (dim > 1) stag->dof[2] = dof2;
1036:   if (dim > 2) stag->dof[3] = dof3;
1037:   return(0);
1038: }

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

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

1045:   Input Parameters:
1046: + dm - the DMStag object
1047: - nRanks0,nRanks1,nRanks2 - number of ranks in each direction

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

1052:   Level: developer

1054: .seealso: DMSTAG, DMDASetNumProcs()
1055: @*/
1056: PetscErrorCode DMStagSetNumRanks(DM dm,PetscInt nRanks0,PetscInt nRanks1,PetscInt nRanks2)
1057: {
1058:   PetscErrorCode  ierr;
1059:   DM_Stag * const stag = (DM_Stag*)dm->data;
1060:   PetscInt        dim;

1067:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1068:   DMGetDimension(dm,&dim);
1069:   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");
1070:   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");
1071:   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");
1072:   if (nRanks0) stag->nRanks[0] = nRanks0;
1073:   if (dim > 1 && nRanks1) stag->nRanks[1] = nRanks1;
1074:   if (dim > 2 && nRanks2) stag->nRanks[2] = nRanks2;
1075:   return(0);
1076: }

1078: /*@C
1079:   DMStagSetStencilType - set elementwise ghost/halo stencil type

1081:   Logically Collective; stencilType must contain common value

1083:   Input Parameters:
1084: + dm - the DMStag object
1085: - stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE

1087:   Level: beginner

1089: .seealso: DMSTAG, DMStagGetStencilType(), DMStagStencilType, DMDASetStencilType()
1090: @*/
1091: PetscErrorCode DMStagSetStencilType(DM dm,DMStagStencilType stencilType)
1092: {
1093:   DM_Stag * const stag = (DM_Stag*)dm->data;

1098:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1099:   stag->stencilType = stencilType;
1100:   return(0);
1101: }

1103: /*@C
1104:   DMStagSetStencilWidth - set elementwise stencil width

1106:   Logically Collective; stencilWidth must contain common value

1108:   Input Parameters:
1109: + dm - the DMStag object
1110: - stencilWidth - stencil/halo/ghost width in elements

1112:   Level: beginner

1114: .seealso: DMSTAG, DMDASetStencilWidth()
1115: @*/
1116: PetscErrorCode DMStagSetStencilWidth(DM dm,PetscInt stencilWidth)
1117: {
1118:   DM_Stag * const stag = (DM_Stag*)dm->data;

1123:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1124:   if (stencilWidth < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Stencil width must be non-negative");
1125:   stag->stencilWidth = stencilWidth;
1126:   return(0);
1127: }

1129: /*@C
1130:   DMStagSetGlobalSizes - set global element counts in each direction

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

1134:   Input Parameters:
1135: + dm - the DMStag object
1136: - N0,N1,N2 - global elementwise sizes

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

1141:   Level: advanced

1143: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMDASetSizes()
1144: @*/
1145: PetscErrorCode DMStagSetGlobalSizes(DM dm,PetscInt N0,PetscInt N1,PetscInt N2)
1146: {
1147:   PetscErrorCode  ierr;
1148:   DM_Stag * const stag = (DM_Stag*)dm->data;
1149:   PetscInt        dim;

1153:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1154:   DMGetDimension(dm,&dim);
1155:   if (N0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in X direction must be positive");
1156:   if (dim > 1 && N1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Y direction must be positive");
1157:   if (dim > 2 && N2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Z direction must be positive");
1158:   if (N0) stag->N[0] = N0;
1159:   if (N1) stag->N[1] = N1;
1160:   if (N2) stag->N[2] = N2;
1161:   return(0);
1162: }

1164: /*@C
1165:   DMStagSetOwnershipRanges - set elements per rank in each direction

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

1169:   Input Parameters:
1170: + dm - the DMStag object
1171: - lx,ly,lz - element counts for each rank in each direction

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

1176:   Level: developer

1178: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagGetOwnershipRanges(), DMDASetOwnershipRanges()
1179: @*/
1180: PetscErrorCode DMStagSetOwnershipRanges(DM dm,PetscInt const *lx,PetscInt const *ly,PetscInt const *lz)
1181: {
1182:   PetscErrorCode  ierr;
1183:   DM_Stag * const stag = (DM_Stag*)dm->data;
1184:   const PetscInt  *lin[3];
1185:   PetscInt        d,dim;

1189:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1190:   lin[0] = lx; lin[1] = ly; lin[2] = lz;
1191:   DMGetDimension(dm,&dim);
1192:   for (d=0; d<dim; ++d) {
1193:     if (lin[d]) {
1194:       if (stag->nRanks[d] < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of ranks");
1195:       if (!stag->l[d]) {
1196:         PetscMalloc1(stag->nRanks[d], &stag->l[d]);
1197:       }
1198:       PetscArraycpy(stag->l[d], lin[d], stag->nRanks[d]);
1199:     }
1200:   }
1201:   return(0);
1202: }

1204: /*@C
1205:   DMStagSetUniformCoordinates - set DMStag coordinates to be a uniform grid

1207:   Collective

1209:   Input Parameters:
1210: + dm - the DMStag object
1211: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

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

1217:   Level: advanced

1219: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType(), DMGetCoordinateDM(), DMGetCoordinates(), DMDASetUniformCoordinates()
1220: @*/
1221: PetscErrorCode DMStagSetUniformCoordinates(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1222: {
1223:   PetscErrorCode  ierr;
1224:   DM_Stag * const stag = (DM_Stag*)dm->data;
1225:   PetscBool       flg_stag,flg_product;

1229:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1230:   if (!stag->coordinateDMType) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"You must first call DMStagSetCoordinateDMType()");
1231:   PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg_stag);
1232:   PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg_product);
1233:   if (flg_stag) {
1234:     DMStagSetUniformCoordinatesExplicit(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1235:   } else if (flg_product) {
1236:     DMStagSetUniformCoordinatesProduct(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1237:   } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported DM Type %s",stag->coordinateDMType);
1238:   return(0);
1239: }

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

1244:   Collective

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

1250:   Notes:
1251:   DMStag supports 2 different types of coordinate DM: either another DMStag, or a DMProduct.
1252:   If the grid is orthogonal, using DMProduct should be more efficient.
1253:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

1255:   Level: beginner

1257: .seealso: DMSTAG, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType()
1258: @*/
1259: PetscErrorCode DMStagSetUniformCoordinatesExplicit(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1260: {
1261:   PetscErrorCode  ierr;
1262:   DM_Stag * const stag = (DM_Stag*)dm->data;
1263:   PetscInt        dim;
1264:   PetscBool       flg;

1268:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1269:   PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg);
1270:   if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1271:   DMStagSetCoordinateDMType(dm,DMSTAG);
1272:   DMGetDimension(dm,&dim);
1273:   switch (dim) {
1274:     case 1: DMStagSetUniformCoordinatesExplicit_1d(dm,xmin,xmax);                     break;
1275:     case 2: DMStagSetUniformCoordinatesExplicit_2d(dm,xmin,xmax,ymin,ymax);           break;
1276:     case 3: DMStagSetUniformCoordinatesExplicit_3d(dm,xmin,xmax,ymin,ymax,zmin,zmax); break;
1277:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1278:   }
1279:   return(0);
1280: }

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

1285:   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
1286:   coordinates.

1288:   Collective

1290:   Input Parameters:
1291: + dm - the DMStag object
1292: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

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

1297:   Level: intermediate

1299: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetCoordinateDMType()
1300: @*/
1301: PetscErrorCode DMStagSetUniformCoordinatesProduct(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1302: {
1303:   PetscErrorCode  ierr;
1304:   DM_Stag * const stag = (DM_Stag*)dm->data;
1305:   DM              dmc;
1306:   PetscInt        dim,d,dof0,dof1;
1307:   PetscBool       flg;

1311:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1312:   PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg);
1313:   if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1314:   DMStagSetCoordinateDMType(dm,DMPRODUCT);
1315:   DMGetCoordinateDM(dm,&dmc);
1316:   DMGetDimension(dm,&dim);

1318:   /* Create 1D sub-DMs, living on subcommunicators */

1320:   dof0 = 0;
1321:   for (d=0; d<dim; ++d) {          /* Include point dof in the sub-DMs if any non-element strata are active in the original DMStag */
1322:     if (stag->dof[d]) {
1323:       dof0 = 1;
1324:       break;
1325:     }
1326:   }
1327:   dof1 = stag->dof[dim] ? 1 : 0; /*  Include element dof in the sub-DMs if the elements are active in the original DMStag */

1329:   for (d=0; d<dim; ++d) {
1330:     DM                subdm;
1331:     MPI_Comm          subcomm;
1332:     PetscMPIInt       color;
1333:     const PetscMPIInt key = 0; /* let existing rank break ties */

1335:     /* Choose colors based on position in the plane orthogonal to this dim, and split */
1336:     switch (d) {
1337:       case 0: color = (dim > 1 ? stag->rank[1] : 0)  + (dim > 2 ? stag->nRanks[1]*stag->rank[2] : 0); break;
1338:       case 1: color =            stag->rank[0]       + (dim > 2 ? stag->nRanks[0]*stag->rank[2] : 0); break;
1339:       case 2: color =            stag->rank[0]       +            stag->nRanks[0]*stag->rank[1]     ; break;
1340:       default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1341:     }
1342:     MPI_Comm_split(PetscObjectComm((PetscObject)dm),color,key,&subcomm);

1344:     /* Create sub-DMs living on these new communicators (which are destroyed by DMProduct) */
1345:     DMStagCreate1d(subcomm,stag->boundaryType[d],stag->N[d],dof0,dof1,stag->stencilType,stag->stencilWidth,stag->l[d],&subdm);
1346:     DMSetUp(subdm);
1347:     switch (d) {
1348:       case 0:
1349:         DMStagSetUniformCoordinatesExplicit(subdm,xmin,xmax,0,0,0,0);
1350:         break;
1351:       case 1:
1352:         DMStagSetUniformCoordinatesExplicit(subdm,ymin,ymax,0,0,0,0);
1353:         break;
1354:       case 2:
1355:         DMStagSetUniformCoordinatesExplicit(subdm,zmin,zmax,0,0,0,0);
1356:         break;
1357:       default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1358:     }
1359:     DMProductSetDM(dmc,d,subdm);
1360:     DMProductSetDimensionIndex(dmc,d,0);
1361:     DMDestroy(&subdm);
1362:     MPI_Comm_free(&subcomm);
1363:   }
1364:   return(0);
1365: }

1367: /*@C
1368:   DMStagVecGetArray - get access to local array

1370:   Logically Collective

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

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

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

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

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

1390:   Input Parameters:
1391: + dm - the DMStag object
1392: - vec - the Vec object

1394:   Output Parameters:
1395: . array - the array

1397:   Notes:
1398:   DMStagVecRestoreArray() must be called, once finished with the array

1400:   Level: beginner

1402: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArray(), DMDAVecGetArrayDOF()
1403: @*/
1404: PetscErrorCode DMStagVecGetArray(DM dm,Vec vec,void *array)
1405: {
1406:   PetscErrorCode  ierr;
1407:   DM_Stag * const stag = (DM_Stag*)dm->data;
1408:   PetscInt        dim;
1409:   PetscInt        nLocal;

1414:   DMGetDimension(dm,&dim);
1415:   VecGetLocalSize(vec,&nLocal);
1416:   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);
1417:   switch (dim) {
1418:     case 1:
1419:       VecGetArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1420:       break;
1421:     case 2:
1422:       VecGetArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1423:       break;
1424:     case 3:
1425:       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);
1426:       break;
1427:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1428:   }
1429:   return(0);
1430: }

1432: /*@C
1433:   DMStagVecGetArrayRead - get read-only access to a local array

1435:   Logically Collective

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

1439:   Input Parameters:
1440: + dm - the DMStag object
1441: - vec - the Vec object

1443:   Output Parameters:
1444: . array - the read-only array

1446:   Notes:
1447:   DMStagVecRestoreArrayRead() must be called, once finished with the array

1449:   Level: beginner

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

1463:   DMGetDimension(dm,&dim);
1464:   VecGetLocalSize(vec,&nLocal);
1465:   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);
1466:   switch (dim) {
1467:     case 1:
1468:       VecGetArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1469:       break;
1470:     case 2:
1471:       VecGetArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1472:       break;
1473:     case 3:
1474:       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);
1475:       break;
1476:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1477:   }
1478:   return(0);
1479: }

1481: /*@C
1482:   DMStagVecRestoreArray - restore access to a raw array

1484:   Logically Collective

1486:   Input Parameters:
1487: + dm - the DMStag object
1488: - vec - the Vec object

1490:   Output Parameters:
1491: . array - the array

1493:   Level: beginner

1495: .seealso: DMSTAG, DMStagVecGetArray(), DMDAVecRestoreArray(), DMDAVecRestoreArrayDOF()
1496: @*/
1497: PetscErrorCode DMStagVecRestoreArray(DM dm,Vec vec,void *array)
1498: {
1499:   PetscErrorCode  ierr;
1500:   DM_Stag * const stag = (DM_Stag*)dm->data;
1501:   PetscInt        dim;
1502:   PetscInt        nLocal;

1507:   DMGetDimension(dm,&dim);
1508:   VecGetLocalSize(vec,&nLocal);
1509:   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);
1510:   switch (dim) {
1511:     case 1:
1512:       VecRestoreArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1513:       break;
1514:     case 2:
1515:       VecRestoreArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1516:       break;
1517:     case 3:
1518:       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);
1519:       break;
1520:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1521:   }
1522:   return(0);
1523: }

1525: /*@C
1526:   DMStagVecRestoreArrayRead - restore read-only access to a raw array

1528:   Logically Collective

1530:   Input Parameters:
1531: + dm - the DMStag object
1532: - vec - the Vec object

1534:   Output Parameters:
1535: . array - the read-only array

1537:   Level: beginner

1539: .seealso: DMSTAG, DMStagVecGetArrayRead(), DMDAVecRestoreArrayRead(), DMDAVecRestoreArrayDOFRead()
1540: @*/
1541: PetscErrorCode DMStagVecRestoreArrayRead(DM dm,Vec vec,void *array)
1542: {
1543:   PetscErrorCode  ierr;
1544:   DM_Stag * const stag = (DM_Stag*)dm->data;
1545:   PetscInt        dim;
1546:   PetscInt        nLocal;

1551:   DMGetDimension(dm,&dim);
1552:   VecGetLocalSize(vec,&nLocal);
1553:   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);
1554:   switch (dim) {
1555:     case 1:
1556:       VecRestoreArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1557:       break;
1558:     case 2:
1559:       VecRestoreArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1560:       break;
1561:     case 3:
1562:       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);
1563:       break;
1564:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1565:   }
1566:   return(0);
1567: }