Actual source code: stagutils.c

petsc-master 2020-01-28
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];

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

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

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

 91:   Logically Collective

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

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

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

105:   Input Parameter:
106: . dm - the DMStag object

108:   Output Parameters:
109: . arrX,arrY,arrZ - local 1D coordinate arrays

111:   Level: intermediate

113: .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArraysRead(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagGetProductCoordinateLocationSlot()
114: @*/
115: PetscErrorCode DMStagGetProductCoordinateArrays(DM dm,void* arrX,void* arrY,void* arrZ)
116: {

120:   DMStagGetProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_FALSE);
121:   return(0);
122: }

124: /*@C
125:   DMStagGetProductCoordinateArraysRead - extract product coordinate arrays, read-only

127:   Logically Collective

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

131:   Input Parameter:
132: . dm - the DMStag object

134:   Output Parameters:
135: . arrX,arrY,arrZ - local 1D coordinate arrays

137:   Level: intermediate

139: .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArrays(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagGetProductCoordinateLocationSlot()
140: @*/
141: PetscErrorCode DMStagGetProductCoordinateArraysRead(DM dm,void* arrX,void* arrY,void* arrZ)
142: {

146:   DMStagGetProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_TRUE);
147:   return(0);
148: }

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

153:   Not Collective

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

158:   Input Parameters:
159: + dm - the DMStag object
160: - loc - the grid location

162:   Output Parameter:
163: . slot - the index to use in local arrays

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

169:   Level: intermediate

171: .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead(), DMStagSetUniformCoordinates()
172: @*/
173: PETSC_EXTERN PetscErrorCode DMStagGetProductCoordinateLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt *slot)
174: {
176:   DM             dmCoord;
177:   PetscInt       dim,dofCheck[DMSTAG_MAX_STRATA],s,d;

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

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

221:   Not Collective

223:   Input Parameter:
224: . dm - the DMStag object

226:   Output Parameters:
227: + x,y,z - starting element indices in each direction
228: . m,n,p - element widths in each direction
229: - nExtrax,nExtray,nExtraz - number of extra partial elements in each direction.

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

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

238:   Level: beginner

240: .seealso: DMSTAG, DMStagGetGhostCorners(), DMDAGetCorners()
241: @*/
242: PetscErrorCode DMStagGetCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *nExtrax,PetscInt *nExtray,PetscInt *nExtraz)
243: {
244:   const DM_Stag * const stag = (DM_Stag*)dm->data;

248:   if (x) *x = stag->start[0];
249:   if (y) *y = stag->start[1];
250:   if (z) *z = stag->start[2];
251:   if (m) *m = stag->n[0];
252:   if (n) *n = stag->n[1];
253:   if (p) *p = stag->n[2];
254:   if (nExtrax) *nExtrax = stag->boundaryType[0] != DM_BOUNDARY_PERIODIC && stag->lastRank[0] ? 1 : 0;
255:   if (nExtray) *nExtray = stag->boundaryType[1] != DM_BOUNDARY_PERIODIC && stag->lastRank[1] ? 1 : 0;
256:   if (nExtraz) *nExtraz = stag->boundaryType[2] != DM_BOUNDARY_PERIODIC && stag->lastRank[2] ? 1 : 0;
257:   return(0);
258: }

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

263:   Not Collective

265:   Input Parameter:
266: . dm - the DMStag object

268:   Output Parameters:
269: + dof0 - the number of points per 0-cell (vertex/node)
270: . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
271: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
272: - dof3 - the number of points per 3-cell (element in 3D)

274:   Level: beginner

276: .seealso: DMSTAG, DMStagGetCorners(), DMStagGetGhostCorners(), DMStagGetGlobalSizes(), DMStagGetStencilWidth(), DMStagGetBoundaryTypes(), DMStagGetLocationDof(), DMDAGetDof()
277: @*/
278: PetscErrorCode DMStagGetDOF(DM dm,PetscInt *dof0,PetscInt *dof1,PetscInt *dof2,PetscInt *dof3)
279: {
280:   const DM_Stag * const stag = (DM_Stag*)dm->data;

284:   if (dof0) *dof0 = stag->dof[0];
285:   if (dof1) *dof1 = stag->dof[1];
286:   if (dof2) *dof2 = stag->dof[2];
287:   if (dof3) *dof3 = stag->dof[3];
288:   return(0);
289: }

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

294:   Not Collective

296:   Input Argument:
297: . dm - the DMStag object

299:   Output Arguments:
300: + x,y,z - starting element indices in each direction
301: - m,n,p - element widths in each direction

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

306:   Level: beginner

308: .seealso: DMSTAG, DMStagGetCorners(), DMDAGetGhostCorners()
309: @*/
310: PetscErrorCode DMStagGetGhostCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
311: {
312:   const DM_Stag * const stag = (DM_Stag*)dm->data;

316:   if (x) *x = stag->startGhost[0];
317:   if (y) *y = stag->startGhost[1];
318:   if (z) *z = stag->startGhost[2];
319:   if (m) *m = stag->nGhost[0];
320:   if (n) *n = stag->nGhost[1];
321:   if (p) *p = stag->nGhost[2];
322:   return(0);
323: }

325: /*@C
326:   DMStagGetGlobalSizes - get global element counts

328:   Not Collective

330:   Input Parameter:
331: . dm - the DMStag object

333:   Output Parameters:
334: . M,N,P - global element counts in each direction

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

339:   Level: beginner

341: .seealso: DMSTAG, DMStagGetLocalSizes(), DMDAGetInfo()
342: @*/
343: PetscErrorCode DMStagGetGlobalSizes(DM dm,PetscInt* M,PetscInt* N,PetscInt* P)
344: {
345:   const DM_Stag * const stag = (DM_Stag*)dm->data;

349:   if (M) *M = stag->N[0];
350:   if (N) *N = stag->N[1];
351:   if (P) *P = stag->N[2];
352:   return(0);
353: }

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

358:   Not Collective

360:   Input Parameter:
361: . dm - the DMStag object

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

366:   Level: intermediate

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

371: .seealso: DMSTAG, DMStagGetIsLastRank()
372: @*/
373: PetscErrorCode DMStagGetIsFirstRank(DM dm,PetscBool *isFirstRank0,PetscBool *isFirstRank1,PetscBool *isFirstRank2)
374: {
375:   const DM_Stag * const stag = (DM_Stag*)dm->data;

379:   if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
380:   if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
381:   if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
382:   return(0);
383: }

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

388:   Not Collective

390:   Input Parameter:
391: . dm - the DMStag object

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

396:   Level: intermediate

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

402: .seealso: DMSTAG, DMStagGetIsFirstRank()
403: @*/
404: PetscErrorCode DMStagGetIsLastRank(DM dm,PetscBool *isLastRank0,PetscBool *isLastRank1,PetscBool *isLastRank2)
405: {
406:   const DM_Stag * const stag = (DM_Stag*)dm->data;

410:   if (isLastRank0) *isLastRank0 = stag->lastRank[0];
411:   if (isLastRank1) *isLastRank1 = stag->lastRank[1];
412:   if (isLastRank2) *isLastRank2 = stag->lastRank[2];
413:   return(0);
414: }

416: /*@C
417:   DMStagGetLocalSizes - get local elementwise sizes

419:   Not Collective

421:   Input Parameter:
422: . dm - the DMStag object

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

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

431:   Level: beginner

433: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetDOF(), DMStagGetNumRanks(), DMDAGetLocalInfo()
434: @*/
435: PetscErrorCode DMStagGetLocalSizes(DM dm,PetscInt* m,PetscInt* n,PetscInt* p)
436: {
437:   const DM_Stag * const stag = (DM_Stag*)dm->data;

441:   if (m) *m = stag->n[0];
442:   if (n) *n = stag->n[1];
443:   if (p) *p = stag->n[2];
444:   return(0);
445: }

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

450:   Not Collective

452:   Input Parameter:
453: . dm - the DMStag object

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

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

462:   Level: beginner

464: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetLocalSize(), DMStagSetNumRank(), DMDAGetInfo()
465: @*/
466: PetscErrorCode DMStagGetNumRanks(DM dm,PetscInt *nRanks0,PetscInt *nRanks1,PetscInt *nRanks2)
467: {
468:   const DM_Stag * const stag = (DM_Stag*)dm->data;

472:   if (nRanks0) *nRanks0 = stag->nRanks[0];
473:   if (nRanks1) *nRanks1 = stag->nRanks[1];
474:   if (nRanks2) *nRanks2 = stag->nRanks[2];
475:   return(0);
476: }

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

481:   Not Collective

483:   Input Parameter:
484: . dm - the DMStag object

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

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

493:   Level: developer

495: .seealso: DMSTAG, DMStagGetDOF()
496: @*/
497: PetscErrorCode DMStagGetEntriesPerElement(DM dm,PetscInt *entriesPerElement)
498: {
499:   const DM_Stag * const stag = (DM_Stag*)dm->data;

503:   if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
504:   return(0);
505: }

507: /*@C
508:   DMStagGetStencilType - get elementwise ghost/halo stencil type

510:   Not Collective

512:   Input Parameter:
513: . dm - the DMStag object

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

518:   Level: beginner

520: .seealso: DMSTAG, DMStagSetStencilType(), DMStagStencilType, DMDAGetInfo()
521: @*/
522: PetscErrorCode DMStagGetStencilType(DM dm,DMStagStencilType *stencilType)
523: {
524:   DM_Stag * const stag = (DM_Stag*)dm->data;

528:   *stencilType = stag->stencilType;
529:   return(0);
530: }

532: /*@C
533:   DMStagGetStencilWidth - get elementwise stencil width

535:   Not Collective

537:   Input Parameter:
538: . dm - the DMStag object

540:   Output Parameters:
541: . stencilWidth - stencil/halo/ghost width in elements

543:   Level: beginner

545: .seealso: DMSTAG, DMDAGetStencilWidth(), DMDAGetInfo()
546: @*/
547: PetscErrorCode DMStagGetStencilWidth(DM dm,PetscInt *stencilWidth)
548: {
549:   const DM_Stag * const stag = (DM_Stag*)dm->data;

553:   if (stencilWidth) *stencilWidth = stag->stencilWidth;
554:   return(0);
555: }

557: /*@C
558:   DMStagGetOwnershipRanges - get elements per rank in each direction

560:   Not Collective

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

565:   Output Parameters:
566: +     lx - ownership along x direction (optional)
567: .     ly - ownership along y direction (optional)
568: -     lz - ownership along z direction (optional)

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

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

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

578:   Level: intermediate

580: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagSetOwnershipRanges(), DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDAGetOwnershipRanges()
581: @*/
582: PetscErrorCode DMStagGetOwnershipRanges(DM dm,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
583: {
584:   const DM_Stag * const stag = (DM_Stag*)dm->data;

588:   if (lx) *lx = stag->l[0];
589:   if (ly) *ly = stag->l[1];
590:   if (lz) *lz = stag->l[2];
591:   return(0);
592: }

594: /*@C
595:   DMStagCreateCompatibleDMStag - create a compatible DMStag with different dof/stratum

597:   Collective

599:   Input Parameters
600: + dm - the DMStag object
601: - dof0,dof1,dof2,dof3 - number of dof on each stratum in the new DMStag

603:   Output Parameters:
604: . newdm - the new, compatible DMStag

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

610:   Level: intermediate

612: .seealso: DMSTAG, DMDACreateCompatibleDMDA(), DMGetCompatibility(), DMStagMigrateVec()
613: @*/
614: PetscErrorCode DMStagCreateCompatibleDMStag(DM dm,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DM *newdm)
615: {
616:   PetscErrorCode  ierr;
617:   const DM_Stag * const stag = (DM_Stag*)dm->data;
618:   PetscInt        dim;

622:   DMGetDimension(dm,&dim);
623:   switch (dim) {
624:     case 1:
625:       DMStagCreate1d(PetscObjectComm((PetscObject)dm),stag->boundaryType[0],stag->N[0],dof0,dof1,stag->stencilType,stag->stencilWidth,NULL,newdm);
626:       break;
627:     case 2:
628:       DMStagCreate2d(PetscObjectComm((PetscObject)dm),stag->boundaryType[0],stag->boundaryType[1],stag->N[0],stag->N[1],stag->nRanks[0],stag->nRanks[1],dof0,dof1,dof2,stag->stencilType,stag->stencilWidth,NULL,NULL,newdm);
629:       break;
630:     case 3:
631:       DMStagCreate3d(PetscObjectComm((PetscObject)dm),stag->boundaryType[0],stag->boundaryType[1],stag->boundaryType[2],stag->N[0],stag->N[1],stag->N[2],stag->nRanks[0],stag->nRanks[1],stag->nRanks[2],dof0,dof1,dof2,dof3,stag->stencilType,stag->stencilWidth,NULL,NULL,NULL,newdm);
632:       break;
633:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
634:   }
635:   DMSetUp(*newdm);
636:   return(0);
637: }

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

642:   Not Collective

644:   Input Parameters:
645: + dm - the DMStag object
646: . loc - location relative to an element
647: - c - component

649:   Output Parameter:
650: . slot - index to use

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

657:   Level: beginner

659: .seealso: DMSTAG, DMStagVecGetArray(), DMStagVecGetArrayRead(), DMStagGetDOF(), DMStagGetEntriesPerElement()
660: @*/
661: PetscErrorCode DMStagGetLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt c,PetscInt *slot)
662: {
663:   DM_Stag * const stag = (DM_Stag*)dm->data;

667: #if defined(PETSC_USE_DEBUG)
668:   {
670:     PetscInt       dof;
671:     DMStagGetLocationDOF(dm,loc,&dof);
672:     if (dof < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Location %s has no dof attached",DMStagStencilLocations[loc]);
673:     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);
674:   }
675: #endif
676:   *slot = stag->locationOffsets[loc] + c;
677:   return(0);
678: }

680: /*@C
681:   DMStagMigrateVec - transfer a vector associated with a DMStag to a vector associated with a compatible DMStag

683:   Collective

685:   Input Parameters:
686: + dm - the source DMStag object
687: . vec - the source vector, compatible with dm
688: . dmTo - the compatible destination DMStag object
689: - vecTo - the destination vector, compatible with dmTo

691:   Notes:
692:   Extra dof are ignored, and unfilled dof are zeroed.
693:   Currently only implemented to migrate global vectors to global vectors.

695:   Level: advanced

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

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

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

800: /*@C
801:   DMStagPopulateLocalToGlobalInjective - populate an internal 1-to-1 local-to-global map

803:   Collective

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

814:   Input Parameter:
815: . dm - the DMStag object

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

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

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

828:   Level: developer

830: .seealso: DMSTAG, DMLocalToGlobal(), VecScatter
831: @*/
832: PetscErrorCode DMStagPopulateLocalToGlobalInjective(DM dm)
833: {
834:   PetscErrorCode  ierr;
835:   PetscInt        dim;
836:   DM_Stag * const stag  = (DM_Stag*)dm->data;

840:   if (stag->ltog_injective) return(0); /* Don't re-populate */
841:   DMGetDimension(dm,&dim);
842:   switch (dim) {
843:     case 1: DMStagPopulateLocalToGlobalInjective_1d(dm); break;
844:     case 2: DMStagPopulateLocalToGlobalInjective_2d(dm); break;
845:     case 3: DMStagPopulateLocalToGlobalInjective_3d(dm); break;
846:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
847:   }
848:   return(0);
849: }

851: static PetscErrorCode DMStagRestoreProductCoordinateArrays_Private(DM dm,void *arrX,void *arrY,void *arrZ,PetscBool read)
852: {
853:   PetscErrorCode  ierr;
854:   PetscInt        dim,d;
855:   void*           arr[DMSTAG_MAX_DIM];
856:   DM              dmCoord;

860:   DMGetDimension(dm,&dim);
861:   if (dim > DMSTAG_MAX_DIM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for %D dimensions",dim);
862:   arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
863:   DMGetCoordinateDM(dm,&dmCoord);
864:   for (d=0; d<dim; ++d) {
865:     DM  subDM;
866:     Vec coord1d_local;

868:     DMProductGetDM(dmCoord,d,&subDM);
869:     DMGetCoordinatesLocal(subDM,&coord1d_local);
870:     if (read) {
871:       DMStagVecRestoreArrayRead(subDM,coord1d_local,arr[d]);
872:     } else {
873:       DMStagVecRestoreArray(subDM,coord1d_local,arr[d]);
874:     }
875:   }
876:   return(0);
877: }

879: /*@C
880:   DMStagRestoreProductCoordinateArrays - restore local array access

882:   Logically Collective

884:   Input Parameter:
885: . dm - the DMStag object

887:   Output Parameters:
888: . arrX,arrY,arrZ - local 1D coordinate arrays

890:   Level: intermediate

892:   Notes:
893:   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:

895: $   DMGetCoordinateDM(dm,&cdm);
896: $   for (d=0; d<3; ++d) {
897: $     DM  subdm;
898: $     Vec coor,coor_local;

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

908: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
909: @*/
910: PetscErrorCode DMStagRestoreProductCoordinateArrays(DM dm,void *arrX,void *arrY,void *arrZ)
911: {
912:   PetscErrorCode  ierr;

915:   DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_FALSE);
916:   return(0);
917: }

919: /*@C
920:   DMStagRestoreProductCoordinateArraysRead - restore local product array access, read-only

922:   Logically Collective

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

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

930:   Level: intermediate

932: .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
933: @*/
934: PetscErrorCode DMStagRestoreProductCoordinateArraysRead(DM dm,void *arrX,void *arrY,void *arrZ)
935: {
936:   PetscErrorCode  ierr;

939:   DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_TRUE);
940:   return(0);
941: }

943: /*@C
944:   DMStagSetBoundaryTypes - set DMStag boundary types

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

948:   Input Parameters:
949: + dm - the DMStag object
950: - boundaryType0,boundaryType1,boundaryType2 - boundary types in each direction

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

955:   Level: advanced

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

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

978: /*@C
979:   DMStagSetCoordinateDMType - set DM type to store coordinates

981:   Logically Collective; dmtype must contain common value

983:   Input Parameters:
984: + dm - the DMStag object
985: - dmtype - DMtype for coordinates, either DMSTAG or DMPRODUCT

987:   Level: advanced

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

998:   PetscFree(stag->coordinateDMType);
999:   PetscStrallocpy(dmtype,(char**)&stag->coordinateDMType);
1000:   return(0);
1001: }

1003: /*@C
1004:   DMStagSetDOF - set dof/stratum

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

1008:   Input Parameters:
1009: + dm - the DMStag object
1010: - dof0,dof1,dof2,dof3 - dof per stratum

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

1015:   Level: advanced

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

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

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

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

1049:   Input Parameters:
1050: + dm - the DMStag object
1051: - nRanks0,nRanks1,nRanks2 - number of ranks in each direction

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

1056:   Level: developer

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

1071:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1072:   DMGetDimension(dm,&dim);
1073:   if (nRanks0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in X direction cannot be less than 1");
1074:   if (dim > 1 && nRanks1 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in Y direction cannot be less than 1");
1075:   if (dim > 2 && nRanks2 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in Z direction cannot be less than 1");
1076:   if (nRanks0) stag->nRanks[0] = nRanks0;
1077:   if (nRanks1) stag->nRanks[1] = nRanks1;
1078:   if (nRanks2) stag->nRanks[2] = nRanks2;
1079:   return(0);
1080: }

1082: /*@C
1083:   DMStagSetStencilType - set elementwise ghost/halo stencil type

1085:   Logically Collective; stencilType must contain common value

1087:   Input Parameters:
1088: + dm - the DMStag object
1089: - stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE

1091:   Level: beginner

1093: .seealso: DMSTAG, DMStagGetStencilType(), DMStagStencilType, DMDASetStencilType()
1094: @*/
1095: PetscErrorCode DMStagSetStencilType(DM dm,DMStagStencilType stencilType)
1096: {
1097:   DM_Stag * const stag = (DM_Stag*)dm->data;

1102:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1103:   stag->stencilType = stencilType;
1104:   return(0);
1105: }

1107: /*@C
1108:   DMStagSetStencilWidth - set elementwise stencil width

1110:   Logically Collective; stencilWidth must contain common value

1112:   Input Parameters:
1113: + dm - the DMStag object
1114: - stencilWidth - stencil/halo/ghost width in elements

1116:   Level: beginner

1118: .seealso: DMSTAG, DMDASetStencilWidth()
1119: @*/
1120: PetscErrorCode DMStagSetStencilWidth(DM dm,PetscInt stencilWidth)
1121: {
1122:   DM_Stag * const stag = (DM_Stag*)dm->data;

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

1133: /*@C
1134:   DMStagSetGlobalSizes - set global element counts in each direction

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

1138:   Input Parameters:
1139: + dm - the DMStag object
1140: - N0,N1,N2 - global elementwise sizes

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

1145:   Level: advanced

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

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

1168: /*@C
1169:   DMStagSetOwnershipRanges - set elements per rank in each direction

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

1173:   Input Parameters:
1174: + dm - the DMStag object
1175: - lx,ly,lz - element counts for each rank in each direction

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

1180:   Level: developer

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

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

1207: /*@C
1208:   DMStagSetUniformCoordinates - set DMStag coordinates to be a uniform grid

1210:   Collective

1212:   Input Parameters:
1213: + dm - the DMStag object
1214: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

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

1220:   Level: advanced

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

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

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

1247:   Collective

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

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

1258:   Level: beginner

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

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

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

1288:   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
1289:   coordinates.

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:   Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.

1300:   Level: intermediate

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

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

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

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

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

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

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

1370: /*@C
1371:   DMStagVecGetArray - get access to local array

1373:   Logically Collective

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

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

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

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

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

1393:   Input Parameters:
1394: + dm - the DMStag object
1395: - vec - the Vec object

1397:   Output Parameters:
1398: . array - the array

1400:   Notes:
1401:   DMStagVecRestoreArray() must be called, once finished with the array

1403:   Level: beginner

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

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

1435: /*@C
1436:   DMStagVecGetArrayRead - get read-only access to a local array

1438:   Logically Collective

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

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

1446:   Output Parameters:
1447: . array - the read-only array

1449:   Notes:
1450:   DMStagVecRestoreArrayRead() must be called, once finished with the array

1452:   Level: beginner

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

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

1484: /*@C
1485:   DMStagVecRestoreArray - restore access to a raw array

1487:   Logically Collective

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

1493:   Output Parameters:
1494: . array - the array

1496:   Level: beginner

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

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

1528: /*@C
1529:   DMStagVecRestoreArrayRead - restore read-only access to a raw array

1531:   Logically Collective

1533:   Input Parameters:
1534: + dm - the DMStag object
1535: - vec - the Vec object

1537:   Output Parameters:
1538: . array - the read-only array

1540:   Level: beginner

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

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