Actual source code: stagutils.c

petsc-master 2019-11-13
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: /*@C
 35:   DMStagGet1dCoordinateArraysDOFRead - extract 1D coordinate arrays

 37:   Logically Collective

 39:   A high-level helper function to quickly extract raw 1D local coordinate arrays.
 40:   Checks that the coordinate DM is a DMProduct or 1D DMStags, with the same number of dof.
 41:   Check on the number of dof and dimension ensures that the elementwise data
 42:   is the same for each, so the same indexing can be used on the arrays.

 44:   Input Parameter:
 45: . dm - the DMStag object

 47:   Output Parameters:
 48: . arrX,arrY,arrX - local 1D coordinate arrays

 50:   Level: intermediate

 52: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagGet1dCoordinateLocationSlot()
 53: @*/
 54: PetscErrorCode DMStagGet1dCoordinateArraysDOFRead(DM dm,void* arrX,void* arrY,void* arrZ)
 55: {
 57:   PetscInt       dim,d,dofCheck[DMSTAG_MAX_STRATA],s;
 58:   DM             dmCoord;
 59:   void*          arr[DMSTAG_MAX_DIM];

 63:   DMGetDimension(dm,&dim);
 64:   if (dim > DMSTAG_MAX_DIM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for %D dimensions",dim);
 65:   arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
 66:   DMGetCoordinateDM(dm,&dmCoord);
 67:   if (!dmCoord) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM does not have a coordinate DM");
 68:   {
 69:     PetscBool isProduct;
 70:     DMType    dmType;
 71:     DMGetType(dmCoord,&dmType);
 72:     PetscStrcmp(DMPRODUCT,dmType,&isProduct);
 73:     if (!isProduct) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is not of type DMPRODUCT");
 74:   }
 75:   for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
 76:   for (d=0; d<dim; ++d) {
 77:     DM        subDM;
 78:     DMType    dmType;
 79:     PetscBool isStag;
 80:     PetscInt  dof[DMSTAG_MAX_STRATA],subDim;
 81:     Vec       coord1d;
 82:     DMProductGetDM(dmCoord,d,&subDM);
 83:     if (!subDM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is missing sub DM %D",d);
 84:     DMGetDimension(subDM,&subDim);
 85:     if (subDim != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of dimension 1");
 86:     DMGetType(subDM,&dmType);
 87:     PetscStrcmp(DMSTAG,dmType,&isStag);
 88:     if (!isStag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of type DMSTAG");
 89:     DMStagGetDOF(subDM,&dof[0],&dof[1],&dof[2],&dof[3]);
 90:     if (d == 0) {
 91:       for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
 92:     } else {
 93:       for (s=0; s<DMSTAG_MAX_STRATA; ++s) {
 94:         if (dofCheck[s] != dof[s]) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DMs have different dofs");
 95:       }
 96:     }
 97:     DMGetCoordinatesLocal(subDM,&coord1d);
 98:     DMStagVecGetArrayDOFRead(subDM,coord1d,arr[d]);
 99:   }
100:   return(0);
101: }

103: /*@C
104:   DMStagGet1dCoordinateLocationSlot - get slot for use with local 1D coordinate arrays

106:   High-level helper function to get slot ids for 1D coordinate DMs.
107:   For use with DMStagGetIDCoordinateArraysDOFRead() and related functions.

109:   Not Collective

111:   Input Parameters:
112: + dm - the DMStag object
113: - loc - the grid location

115:   Output Parameter:
116: . slot - the index to use in local arrays

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

122:   Level: intermediate

124: .seealso: DMSTAG, DMPRODUCT, DMStagGet1dCoordinateArraysDOFRead(), DMStagSetUniformCoordinates()
125: @*/
126: PETSC_EXTERN PetscErrorCode DMStagGet1dCoordinateLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt *slot)
127: {
129:   DM             dmCoord;
130:   PetscInt       dim,dofCheck[DMSTAG_MAX_STRATA],s,d;

134:   DMGetDimension(dm,&dim);
135:   DMGetCoordinateDM(dm,&dmCoord);
136:   if (!dmCoord) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM does not have a coordinate DM");
137:   {
138:     PetscBool isProduct;
139:     DMType    dmType;
140:     DMGetType(dmCoord,&dmType);
141:     PetscStrcmp(DMPRODUCT,dmType,&isProduct);
142:     if (!isProduct) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is not of type DMPRODUCT");
143:   }
144:   for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
145:   for (d=0; d<dim; ++d) {
146:     DM        subDM;
147:     DMType    dmType;
148:     PetscBool isStag;
149:     PetscInt  dof[DMSTAG_MAX_STRATA],subDim;
150:     DMProductGetDM(dmCoord,d,&subDM);
151:     if (!subDM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is missing sub DM %D",d);
152:     DMGetDimension(subDM,&subDim);
153:     if (subDim != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of dimension 1");
154:     DMGetType(subDM,&dmType);
155:     PetscStrcmp(DMSTAG,dmType,&isStag);
156:     if (!isStag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of type DMSTAG");
157:     DMStagGetDOF(subDM,&dof[0],&dof[1],&dof[2],&dof[3]);
158:     if (d == 0) {
159:       const PetscInt component = 0;
160:       for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
161:       DMStagGetLocationSlot(subDM,loc,component,slot);
162:     } else {
163:       for (s=0; s<DMSTAG_MAX_STRATA; ++s) {
164:         if (dofCheck[s] != dof[s]) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DMs have different dofs");
165:       }
166:     }
167:   }
168:   return(0);
169: }

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

174:   Not Collective

176:   Input Parameter:
177: . dm - the DMStag object

179:   Output Parameters:
180: + x,y,z - starting element indices in each direction
181: . m,n,p - element widths in each direction
182: - nExtrax,nExtray,nExtraz - number of extra partial elements in each direction.

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

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

191:   Level: beginner

193: .seealso: DMSTAG, DMStagGetGhostCorners(), DMDAGetCorners()
194: @*/
195: PetscErrorCode DMStagGetCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *nExtrax,PetscInt *nExtray,PetscInt *nExtraz)
196: {
197:   const DM_Stag * const stag = (DM_Stag*)dm->data;

201:   if (x) *x = stag->start[0];
202:   if (y) *y = stag->start[1];
203:   if (z) *z = stag->start[2];
204:   if (m) *m = stag->n[0];
205:   if (n) *n = stag->n[1];
206:   if (p) *p = stag->n[2];
207:   if (nExtrax) *nExtrax = stag->boundaryType[0] != DM_BOUNDARY_PERIODIC && stag->lastRank[0] ? 1 : 0;
208:   if (nExtray) *nExtray = stag->boundaryType[1] != DM_BOUNDARY_PERIODIC && stag->lastRank[1] ? 1 : 0;
209:   if (nExtraz) *nExtraz = stag->boundaryType[2] != DM_BOUNDARY_PERIODIC && stag->lastRank[2] ? 1 : 0;
210:   return(0);
211: }

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

216:   Not Collective

218:   Input Parameter:
219: . dm - the DMStag object

221:   Output Parameters:
222: + dof0 - the number of points per 0-cell (vertex/node)
223: . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
224: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
225: - dof3 - the number of points per 3-cell (element in 3D)

227:   Level: beginner

229: .seealso: DMSTAG, DMStagGetCorners(), DMStagGetGhostCorners(), DMStagGetGlobalSizes(), DMStagGetStencilWidth(), DMStagGetBoundaryTypes(), DMStagGetLocationDof(), DMDAGetDof()
230: @*/
231: PetscErrorCode DMStagGetDOF(DM dm,PetscInt *dof0,PetscInt *dof1,PetscInt *dof2,PetscInt *dof3)
232: {
233:   const DM_Stag * const stag = (DM_Stag*)dm->data;

237:   if (dof0) *dof0 = stag->dof[0];
238:   if (dof1) *dof1 = stag->dof[1];
239:   if (dof2) *dof2 = stag->dof[2];
240:   if (dof3) *dof3 = stag->dof[3];
241:   return(0);
242: }

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

247:   Not Collective

249:   Input Argument:
250: . dm - the DMStag object

252:   Output Arguments:
253: + x,y,z - starting element indices in each direction
254: - m,n,p - element widths in each direction

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

259:   Level: beginner

261: .seealso: DMSTAG, DMStagGetCorners(), DMDAGetGhostCorners()
262: @*/
263: PetscErrorCode DMStagGetGhostCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
264: {
265:   const DM_Stag * const stag = (DM_Stag*)dm->data;

269:   if (x) *x = stag->startGhost[0];
270:   if (y) *y = stag->startGhost[1];
271:   if (z) *z = stag->startGhost[2];
272:   if (m) *m = stag->nGhost[0];
273:   if (n) *n = stag->nGhost[1];
274:   if (p) *p = stag->nGhost[2];
275:   return(0);
276: }

278: /*@C
279:   DMStagGetGlobalSizes - get global element counts

281:   Not Collective

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

286:   Output Parameters:
287: . M,N,P - global element counts in each direction

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

292:   Level: beginner

294: .seealso: DMSTAG, DMStagGetLocalSizes(), DMDAGetInfo()
295: @*/
296: PetscErrorCode DMStagGetGlobalSizes(DM dm,PetscInt* M,PetscInt* N,PetscInt* P)
297: {
298:   const DM_Stag * const stag = (DM_Stag*)dm->data;

302:   if (M) *M = stag->N[0];
303:   if (N) *N = stag->N[1];
304:   if (P) *P = stag->N[2];
305:   return(0);
306: }

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

311:   Not Collective

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

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

319:   Level: intermediate

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

324: .seealso: DMSTAG, DMStagGetIsLastRank()
325: @*/
326: PetscErrorCode DMStagGetIsFirstRank(DM dm,PetscBool *isFirstRank0,PetscBool *isFirstRank1,PetscBool *isFirstRank2)
327: {
328:   const DM_Stag * const stag = (DM_Stag*)dm->data;

332:   if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
333:   if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
334:   if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
335:   return(0);
336: }

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

341:   Not Collective

343:   Input Parameter:
344: . dm - the DMStag object

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

349:   Level: intermediate

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

355: .seealso: DMSTAG, DMStagGetIsFirstRank()
356: @*/
357: PetscErrorCode DMStagGetIsLastRank(DM dm,PetscBool *isLastRank0,PetscBool *isLastRank1,PetscBool *isLastRank2)
358: {
359:   const DM_Stag * const stag = (DM_Stag*)dm->data;

363:   if (isLastRank0) *isLastRank0 = stag->lastRank[0];
364:   if (isLastRank1) *isLastRank1 = stag->lastRank[1];
365:   if (isLastRank2) *isLastRank2 = stag->lastRank[2];
366:   return(0);
367: }

369: /*@C
370:   DMStagGetLocalSizes - get local elementwise sizes

372:   Not Collective

374:   Input Parameter:
375: . dm - the DMStag object

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

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

384:   Level: beginner

386: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetDOF(), DMStagGetNumRanks(), DMDAGetLocalInfo()
387: @*/
388: PetscErrorCode DMStagGetLocalSizes(DM dm,PetscInt* m,PetscInt* n,PetscInt* p)
389: {
390:   const DM_Stag * const stag = (DM_Stag*)dm->data;

394:   if (m) *m = stag->n[0];
395:   if (n) *n = stag->n[1];
396:   if (p) *p = stag->n[2];
397:   return(0);
398: }

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

403:   Not Collective

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

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

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

415:   Level: beginner

417: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetLocalSize(), DMStagSetNumRank(), DMDAGetInfo()
418: @*/
419: PetscErrorCode DMStagGetNumRanks(DM dm,PetscInt *nRanks0,PetscInt *nRanks1,PetscInt *nRanks2)
420: {
421:   const DM_Stag * const stag = (DM_Stag*)dm->data;

425:   if (nRanks0) *nRanks0 = stag->nRanks[0];
426:   if (nRanks1) *nRanks1 = stag->nRanks[1];
427:   if (nRanks2) *nRanks2 = stag->nRanks[2];
428:   return(0);
429: }

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

434:   Not Collective

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

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

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

446:   Level: developer

448: .seealso: DMSTAG, DMStagGetDOF()
449: @*/
450: PetscErrorCode DMStagGetEntriesPerElement(DM dm,PetscInt *entriesPerElement)
451: {
452:   const DM_Stag * const stag = (DM_Stag*)dm->data;

456:   if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
457:   return(0);
458: }

460: /*@C
461:   DMStagGetStencilType - get elementwise ghost/halo stencil type

463:   Not Collective

465:   Input Parameter:
466: . dm - the DMStag object

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

471:   Level: beginner

473: .seealso: DMSTAG, DMStagSetStencilType(), DMStagStencilType, DMDAGetInfo()
474: @*/
475: PetscErrorCode DMStagGetStencilType(DM dm,DMStagStencilType *stencilType)
476: {
477:   DM_Stag * const stag = (DM_Stag*)dm->data;

481:   *stencilType = stag->stencilType;
482:   return(0);
483: }

485: /*@C
486:   DMStagGetStencilWidth - get elementwise stencil width

488:   Not Collective

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

493:   Output Parameters:
494: . stencilWidth - stencil/halo/ghost width in elements

496:   Level: beginner

498: .seealso: DMSTAG, DMDAGetStencilWidth(), DMDAGetInfo()
499: @*/
500: PetscErrorCode DMStagGetStencilWidth(DM dm,PetscInt *stencilWidth)
501: {
502:   const DM_Stag * const stag = (DM_Stag*)dm->data;

506:   if (stencilWidth) *stencilWidth = stag->stencilWidth;
507:   return(0);
508: }

510: /*@C
511:   DMStagGetOwnershipRanges - get elements per rank in each direction

513:   Not Collective

515:   Input Parameter:
516: .     dm - the DMStag object

518:   Output Parameters:
519: +     lx - ownership along x direction (optional)
520: .     ly - ownership along y direction (optional)
521: -     lz - ownership along z direction (optional)

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

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

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

531:   Level: intermediate

533: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagSetOwnershipRanges(), DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDAGetOwnershipRanges()
534: @*/
535: PetscErrorCode DMStagGetOwnershipRanges(DM dm,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
536: {
537:   const DM_Stag * const stag = (DM_Stag*)dm->data;

541:   if (lx) *lx = stag->l[0];
542:   if (ly) *ly = stag->l[1];
543:   if (lz) *lz = stag->l[2];
544:   return(0);
545: }

547: /*@C
548:   DMStagCreateCompatibleDMStag - create a compatible DMStag with different dof/stratum

550:   Collective

552:   Input Parameters
553: + dm - the DMStag object
554: - dof0,dof1,dof2,dof3 - number of dof on each stratum in the new DMStag

556:   Output Parameters:
557: . newdm - the new, compatible DMStag

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

563:   Level: intermediate

565: .seealso: DMSTAG, DMDACreateCompatibleDMDA(), DMGetCompatibility(), DMStagMigrateVec()
566: @*/
567: PetscErrorCode DMStagCreateCompatibleDMStag(DM dm,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DM *newdm)
568: {
569:   PetscErrorCode  ierr;
570:   const DM_Stag * const stag = (DM_Stag*)dm->data;
571:   PetscInt        dim;

575:   DMGetDimension(dm,&dim);
576:   switch (dim) {
577:     case 1:
578:       DMStagCreate1d(PetscObjectComm((PetscObject)dm),stag->boundaryType[0],stag->N[0],dof0,dof1,stag->stencilType,stag->stencilWidth,NULL,newdm);
579:       break;
580:     case 2:
581:       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);
582:       break;
583:     case 3:
584:       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);
585:       break;
586:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
587:   }
588:   DMSetUp(*newdm);
589:   return(0);
590: }

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

595:   Not Collective

597:   Input Parameters:
598: + dm - the DMStag object
599: . loc - location relative to an element
600: - c - component

602:   Output Parameter:
603: . slot - index to use

605:   Notes:
606:   Provides an appropriate index to use with DMStagVecGetArrayDOF() and friends.
607:   This is required so that the user doesn't need to know about the ordering of
608:   dof associated with each local element.

610:   Level: beginner

612: .seealso: DMSTAG, DMStagVecGetArrayDOF(), DMStagVecGetArrayDOFRead(), DMStagGetDOF(), DMStagGetEntriesPerElement()
613: @*/
614: PetscErrorCode DMStagGetLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt c,PetscInt *slot)
615: {
616:   DM_Stag * const stag = (DM_Stag*)dm->data;

620: #if defined(PETSC_USE_DEBUG)
621:   {
623:     PetscInt       dof;
624:     DMStagGetLocationDOF(dm,loc,&dof);
625:     if (dof < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Location %s has no dof attached",DMStagStencilLocations[loc]);
626:     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);
627:   }
628: #endif
629:   *slot = stag->locationOffsets[loc] + c;
630:   return(0);
631: }

633: /*@C
634:   DMStagMigrateVec - transfer a vector associated with a DMStag to a vector associated with a compatible DMStag

636:   Collective

638:   Input Parameters:
639: + dm - the source DMStag object
640: . vec - the source vector, compatible with dm
641: . dmTo - the compatible destination DMStag object
642: - vecTo - the destination vector, compatible with dmTo

644:   Notes:
645:   Extra dof are ignored, and unfilled dof are zeroed.
646:   Currently only implemented to migrate global vectors to global vectors.

648:   Level: advanced

650: .seealso: DMSTAG, DMStagCreateCompatibleDMStag(), DMGetCompatibility(), DMStagVecSplitToDMDA()
651: @*/
652: PetscErrorCode DMStagMigrateVec(DM dm,Vec vec,DM dmTo,Vec vecTo)
653: {
654:   PetscErrorCode    ierr;
655:   DM_Stag * const   stag = (DM_Stag*)dm->data;
656:   DM_Stag * const   stagTo = (DM_Stag*)dmTo->data;
657:   PetscInt          nLocalTo,nLocal,dim,i,j,k;
658:   PetscInt          start[DMSTAG_MAX_DIM],startGhost[DMSTAG_MAX_DIM],n[DMSTAG_MAX_DIM],nExtra[DMSTAG_MAX_DIM],offset[DMSTAG_MAX_DIM];
659:   Vec               vecToLocal,vecLocal;
660:   PetscBool         compatible,compatibleSet;
661:   const PetscScalar *arr;
662:   PetscScalar       *arrTo;
663:   const PetscInt    epe   = stag->entriesPerElement;
664:   const PetscInt    epeTo = stagTo->entriesPerElement;

671:   DMGetCompatibility(dm,dmTo,&compatible,&compatibleSet);
672:   if (!compatibleSet || !compatible) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"DMStag objects must be shown to be compatible");
673:   DMGetDimension(dm,&dim);
674:   VecGetLocalSize(vec,&nLocal);
675:   VecGetLocalSize(vecTo,&nLocalTo);
676:   if (nLocal != stag->entries|| nLocalTo !=stagTo->entries) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Vector migration only implemented for global vector to global vector.");
677:   DMStagGetCorners(dm,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
678:   DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],NULL,NULL,NULL);
679:   for (i=0; i<DMSTAG_MAX_DIM; ++i) offset[i] = start[i]-startGhost[i];

681:   /* Proceed by transferring to a local vector, copying, and transferring back to a global vector */
682:   DMGetLocalVector(dm,&vecLocal);
683:   DMGetLocalVector(dmTo,&vecToLocal);
684:   DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal);
685:   DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal);
686:   VecGetArrayRead(vecLocal,&arr);
687:   VecGetArray(vecToLocal,&arrTo);
688:   /* Note that some superfluous copying of entries on partial dummy elements is done */
689:   if (dim == 1) {
690:     for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
691:       PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
692:       PetscInt si;
693:       for (si=0; si<2; ++si) {
694:         b   += stag->dof[si];
695:         bTo += stagTo->dof[si];
696:         for (; d < b && dTo < bTo; ++d,++dTo) arrTo[i*epeTo + dTo] = arr[i*epe + d];
697:         for (; dTo < bTo         ;     ++dTo) arrTo[i*epeTo + dTo] = 0.0;
698:         d = b;
699:       }
700:     }
701:   } else if (dim == 2) {
702:     const PetscInt epr   = stag->nGhost[0] * epe;
703:     const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
704:     for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
705:       for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
706:         const PetscInt base   = j*epr   + i*epe;
707:         const PetscInt baseTo = j*eprTo + i*epeTo;
708:         PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
709:         const PetscInt s[4] = {0,1,1,2}; /* Dimensions of points, in order */
710:         PetscInt si;
711:         for (si=0; si<4; ++si) {
712:             b   += stag->dof[s[si]];
713:             bTo += stagTo->dof[s[si]];
714:             for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
715:             for (;          dTo < bTo;     ++dTo) arrTo[baseTo + dTo] = 0.0;
716:             d = b;
717:         }
718:       }
719:     }
720:   } else if (dim == 3) {
721:     const PetscInt epr   = stag->nGhost[0]   * epe;
722:     const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
723:     const PetscInt epl   = stag->nGhost[1]   * epr;
724:     const PetscInt eplTo = stagTo->nGhost[1] * eprTo;
725:     for (k=offset[2]; k<offset[2] + n[2] + nExtra[2]; ++k) {
726:       for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
727:         for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
728:           PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
729:           const PetscInt base   = k*epl   + j*epr   + i*epe;
730:           const PetscInt baseTo = k*eplTo + j*eprTo + i*epeTo;
731:           const PetscInt s[8] = {0,1,1,2,1,2,2,3}; /* dimensions of points, in order */
732:           PetscInt is;
733:           for (is=0; is<8; ++is) {
734:             b   += stag->dof[s[is]];
735:             bTo += stagTo->dof[s[is]];
736:             for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
737:             for (;          dTo < bTo;     ++dTo) arrTo[baseTo + dTo] = 0.0;
738:             d = b;
739:           }
740:         }
741:       }
742:     }
743:   } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
744:   VecRestoreArrayRead(vecLocal,&arr);
745:   VecRestoreArray(vecToLocal,&arrTo);
746:   DMRestoreLocalVector(dm,&vecLocal);
747:   DMLocalToGlobalBegin(dmTo,vecToLocal,INSERT_VALUES,vecTo);
748:   DMLocalToGlobalEnd(dmTo,vecToLocal,INSERT_VALUES,vecTo);
749:   DMRestoreLocalVector(dmTo,&vecToLocal);
750:   return(0);
751: }

753: /*@C
754:   DMStagPopulateLocalToGlobalInjective - populate an internal 1-to-1 local-to-global map

756:   Collective

758:   Creates an internal object which explicitly maps a single local degree of
759:   freedom to each global degree of freedom. This is used, if populated,
760:   instead of SCATTER_REVERSE_LOCAL with the (1-to-many, in general)
761:   global-to-local map, when DMLocalToGlobal() is called with INSERT_VALUES.
762:   This allows usage, for example, even in the periodic, 1-rank case, where
763:   the inverse of the global-to-local map, even when restricted to on-rank
764:   communication, is non-injective. This is at the cost of storing an additional
765:   VecScatter object inside each DMStag object.

767:   Input Parameter:
768: . dm - the DMStag object

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

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

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

781:   Level: developer

783: .seealso: DMSTAG, DMLocalToGlobal(), VecScatter
784: @*/
785: PetscErrorCode DMStagPopulateLocalToGlobalInjective(DM dm)
786: {
787:   PetscErrorCode  ierr;
788:   PetscInt        dim;
789:   DM_Stag * const stag  = (DM_Stag*)dm->data;

793:   if (stag->ltog_injective) return(0); /* Don't re-populate */
794:   DMGetDimension(dm,&dim);
795:   switch (dim) {
796:     case 1: DMStagPopulateLocalToGlobalInjective_1d(dm); break;
797:     case 2: DMStagPopulateLocalToGlobalInjective_2d(dm); break;
798:     case 3: DMStagPopulateLocalToGlobalInjective_3d(dm); break;
799:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
800:   }
801:   return(0);
802: }

804: /*@C
805:   DMStagRestore1dCoordinateArraysDOFRead - restore local array access

807:   Logically Collective

809:   Input Parameter:
810: . dm - the DMStag object

812:   Output Parameters:
813: . arrX,arrY,arrX - local 1D coordinate arrays

815:   Level: intermediate

817: .seealso: DMSTAG, DMStagGet1dCoordinateArraysDOFRead()
818: @*/
819: PetscErrorCode DMStagRestore1dCoordinateArraysDOFRead(DM dm,void *arrX,void *arrY,void *arrZ)
820: {
821:   PetscErrorCode  ierr;
822:   PetscInt        dim,d;
823:   void*           arr[DMSTAG_MAX_DIM];
824:   DM              dmCoord;

828:   DMGetDimension(dm,&dim);
829:   if (dim > DMSTAG_MAX_DIM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for %D dimensions",dim);
830:   arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
831:   DMGetCoordinateDM(dm,&dmCoord);
832:   for (d=0; d<dim; ++d) {
833:     DM        subDM;
834:     Vec       coord1d;
835:     DMProductGetDM(dmCoord,d,&subDM);
836:     DMGetCoordinatesLocal(subDM,&coord1d);
837:     DMStagVecRestoreArrayDOFRead(subDM,coord1d,arr[d]);
838:   }
839:   return(0);
840: }

842: /*@C
843:   DMStagSetBoundaryTypes - set DMStag boundary types

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

847:   Input Parameters:
848: + dm - the DMStag object
849: - boundaryType0,boundaryType1,boundaryType2 - boundary types in each direction

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

854:   Level: advanced

856: .seealso: DMSTAG, DMBoundaryType, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDASetBoundaryType()
857: @*/
858: PetscErrorCode DMStagSetBoundaryTypes(DM dm,DMBoundaryType boundaryType0,DMBoundaryType boundaryType1,DMBoundaryType boundaryType2)
859: {
860:   PetscErrorCode  ierr;
861:   DM_Stag * const stag  = (DM_Stag*)dm->data;
862:   PetscInt        dim;

869:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
870:   DMGetDimension(dm,&dim);
871:   if (boundaryType0           ) stag->boundaryType[0] = boundaryType0;
872:   if (boundaryType1 && dim > 1) stag->boundaryType[1] = boundaryType1;
873:   if (boundaryType2 && dim > 2) stag->boundaryType[2] = boundaryType2;
874:   return(0);
875: }

877: /*@C
878:   DMStagSetCoordinateDMType - set DM type to store coordinates

880:   Logically Collective; dmtype must contain common value

882:   Input Parameters:
883: + dm - the DMStag object
884: - dmtype - DMtype for coordinates, either DMSTAG or DMPRODUCT

886:   Level: advanced

888: .seealso: DMSTAG, DMPRODUCT, DMGetCoordinateDM(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMType
889: @*/
890: PetscErrorCode DMStagSetCoordinateDMType(DM dm,DMType dmtype)
891: {
892:   PetscErrorCode  ierr;
893:   DM_Stag * const stag = (DM_Stag*)dm->data;

897:   PetscFree(stag->coordinateDMType);
898:   PetscStrallocpy(dmtype,(char**)&stag->coordinateDMType);
899:   return(0);
900: }

902: /*@C
903:   DMStagSetDOF - set dof/stratum

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

907:   Input Parameters:
908: + dm - the DMStag object
909: - dof0,dof1,dof2,dof3 - dof per stratum

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

914:   Level: advanced

916: .seealso: DMSTAG, DMDASetDof()
917: @*/
918: PetscErrorCode DMStagSetDOF(DM dm,PetscInt dof0, PetscInt dof1,PetscInt dof2,PetscInt dof3)
919: {
920:   PetscErrorCode  ierr;
921:   DM_Stag * const stag = (DM_Stag*)dm->data;
922:   PetscInt        dim;

930:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
931:   DMGetDimension(dm,&dim);
932:   if (dof0 < 0)            SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof0 cannot be negative");
933:   if (dof1 < 0)            SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof1 cannot be negative");
934:   if (dim > 1 && dof2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof2 cannot be negative");
935:   if (dim > 2 && dof3 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof3 cannot be negative");
936:                stag->dof[0] = dof0;
937:                stag->dof[1] = dof1;
938:   if (dim > 1) stag->dof[2] = dof2;
939:   if (dim > 2) stag->dof[3] = dof3;
940:   return(0);
941: }

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

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

948:   Input Parameters:
949: + dm - the DMStag object
950: - nRanks0,nRanks1,nRanks2 - number of ranks in each direction

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

955:   Level: developer

957: .seealso: DMSTAG, DMDASetNumProcs()
958: @*/
959: PetscErrorCode DMStagSetNumRanks(DM dm,PetscInt nRanks0,PetscInt nRanks1,PetscInt nRanks2)
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 (nRanks0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in X direction cannot be less than 1");
973:   if (dim > 1 && nRanks1 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in Y direction cannot be less than 1");
974:   if (dim > 2 && nRanks2 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in Z direction cannot be less than 1");
975:   if (nRanks0) stag->nRanks[0] = nRanks0;
976:   if (nRanks1) stag->nRanks[1] = nRanks1;
977:   if (nRanks2) stag->nRanks[2] = nRanks2;
978:   return(0);
979: }

981: /*@C
982:   DMStagSetStencilType - set elementwise ghost/halo stencil type

984:   Logically Collective; stencilType must contain common value

986:   Input Parameters:
987: + dm - the DMStag object
988: - stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE

990:   Level: beginner

992: .seealso: DMSTAG, DMStagGetStencilType(), DMStagStencilType, DMDASetStencilType()
993: @*/
994: PetscErrorCode DMStagSetStencilType(DM dm,DMStagStencilType stencilType)
995: {
996:   DM_Stag * const stag = (DM_Stag*)dm->data;

1001:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1002:   stag->stencilType = stencilType;
1003:   return(0);
1004: }

1006: /*@C
1007:   DMStagSetStencilWidth - set elementwise stencil width

1009:   Logically Collective; stencilWidth must contain common value

1011:   Input Parameters:
1012: + dm - the DMStag object
1013: - stencilWidth - stencil/halo/ghost width in elements

1015:   Level: beginner

1017: .seealso: DMSTAG, DMDASetStencilWidth()
1018: @*/
1019: PetscErrorCode DMStagSetStencilWidth(DM dm,PetscInt stencilWidth)
1020: {
1021:   DM_Stag * const stag = (DM_Stag*)dm->data;

1026:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1027:   if (stencilWidth < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Stencil width must be non-negative");
1028:   stag->stencilWidth = stencilWidth;
1029:   return(0);
1030: }

1032: /*@C
1033:   DMStagSetGlobalSizes - set global element counts in each direction

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

1037:   Input Parameters:
1038: + dm - the DMStag object
1039: - N0,N1,N2 - global elementwise sizes

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

1044:   Level: advanced

1046: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMDASetSizes()
1047: @*/
1048: PetscErrorCode DMStagSetGlobalSizes(DM dm,PetscInt N0,PetscInt N1,PetscInt N2)
1049: {
1050:   PetscErrorCode  ierr;
1051:   DM_Stag * const stag = (DM_Stag*)dm->data;
1052:   PetscInt        dim;

1056:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1057:   DMGetDimension(dm,&dim);
1058:   if (N0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in X direction must be positive");
1059:   if (dim > 1 && N1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Y direction must be positive");
1060:   if (dim > 2 && N2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Z direction must be positive");
1061:   if (N0) stag->N[0] = N0;
1062:   if (N1) stag->N[1] = N1;
1063:   if (N2) stag->N[2] = N2;
1064:   return(0);
1065: }

1067: /*@C
1068:   DMStagSetOwnershipRanges - set elements per rank in each direction

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

1072:   Input Parameters:
1073: + dm - the DMStag object
1074: - lx,ly,lz - element counts for each rank in each direction

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

1079:   Level: developer

1081: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagGetOwnershipRanges(), DMDASetOwnershipRanges()
1082: @*/
1083: PetscErrorCode DMStagSetOwnershipRanges(DM dm,PetscInt const *lx,PetscInt const *ly,PetscInt const *lz)
1084: {
1085:   PetscErrorCode  ierr;
1086:   DM_Stag * const stag = (DM_Stag*)dm->data;
1087:   const PetscInt  *lin[3];
1088:   PetscInt        d;

1092:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1093:   lin[0] = lx; lin[1] = ly; lin[2] = lz;
1094:   for (d=0; d<3; ++d) {
1095:     if (lin[d]) {
1096:       if (stag->nRanks[d] < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of ranks");
1097:       if (!stag->l[d]) {
1098:         PetscMalloc1(stag->nRanks[d], &stag->l[d]);
1099:       }
1100:       PetscArraycpy(stag->l[d], lin[d], stag->nRanks[d]);
1101:     }
1102:   }
1103:   return(0);
1104: }

1106: /*@C
1107:   DMStagSetUniformCoordinates - set DMStag coordinates to be a uniform grid

1109:   Collective

1111:   Input Parameters:
1112: + dm - the DMStag object
1113: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

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

1119:   Level: advanced

1121: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType(), DMGetCoordinateDM(), DMGetCoordinates(), DMDASetUniformCoordinates()
1122: @*/
1123: PetscErrorCode DMStagSetUniformCoordinates(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1124: {
1125:   PetscErrorCode  ierr;
1126:   DM_Stag * const stag = (DM_Stag*)dm->data;
1127:   PetscBool       flg_stag,flg_product;

1131:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1132:   if (!stag->coordinateDMType) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"You must first call DMStagSetCoordinateDMType()");
1133:   PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg_stag);
1134:   PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg_product);
1135:   if (flg_stag) {
1136:     DMStagSetUniformCoordinatesExplicit(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1137:   } else if (flg_product) {
1138:     DMStagSetUniformCoordinatesProduct(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1139:   } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported DM Type %s",stag->coordinateDMType);
1140:   return(0);
1141: }

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

1146:   Collective

1148:   Input Parameters:
1149: + dm - the DMStag object
1150: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

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

1157:   Level: beginner

1159: .seealso: DMSTAG, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType()
1160: @*/
1161: PetscErrorCode DMStagSetUniformCoordinatesExplicit(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1162: {
1163:   PetscErrorCode  ierr;
1164:   DM_Stag * const stag = (DM_Stag*)dm->data;
1165:   PetscInt        dim;
1166:   PetscBool       flg;

1170:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1171:   PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg);
1172:   if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1173:   DMStagSetCoordinateDMType(dm,DMSTAG);
1174:   DMGetDimension(dm,&dim);
1175:   switch (dim) {
1176:     case 1: DMStagSetUniformCoordinatesExplicit_1d(dm,xmin,xmax);                     break;
1177:     case 2: DMStagSetUniformCoordinatesExplicit_2d(dm,xmin,xmax,ymin,ymax);           break;
1178:     case 3: DMStagSetUniformCoordinatesExplicit_3d(dm,xmin,xmax,ymin,ymax,zmin,zmax); break;
1179:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1180:   }
1181:   return(0);
1182: }

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

1187:   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
1188:   coordinates.

1190:   Collective

1192:   Input Parameters:
1193: + dm - the DMStag object
1194: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

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

1199:   Level: intermediate

1201: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetCoordinateDMType()
1202: @*/
1203: PetscErrorCode DMStagSetUniformCoordinatesProduct(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1204: {
1205:   PetscErrorCode  ierr;
1206:   DM_Stag * const stag = (DM_Stag*)dm->data;
1207:   DM              dmc;
1208:   PetscInt        dim,d,dof0,dof1;
1209:   PetscBool       flg;

1213:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1214:   PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg);
1215:   if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1216:   DMStagSetCoordinateDMType(dm,DMPRODUCT);
1217:   DMGetCoordinateDM(dm,&dmc);
1218:   DMGetDimension(dm,&dim);

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

1222:   dof0 = 0;
1223:   for (d=0; d<dim; ++d) {          /* Include point dof in the sub-DMs if any non-element strata are active in the original DMStag */
1224:     if (stag->dof[d]) {
1225:       dof0 = 1;
1226:       break;
1227:     }
1228:   }
1229:   dof1 = stag->dof[dim] ? 1 : 0; /*  Include element dof in the sub-DMs if the elements are active in the original DMStag */

1231:   for (d=0; d<dim; ++d) {
1232:     DM                subdm;
1233:     MPI_Comm          subcomm;
1234:     PetscMPIInt       color;
1235:     const PetscMPIInt key = 0; /* let existing rank break ties */

1237:     /* Choose colors based on position in the plane orthogonal to this dim, and split */
1238:     switch (d) {
1239:       case 0: color = (dim > 1 ? stag->rank[1] : 0)  + (dim > 2 ? stag->nRanks[1]*stag->rank[2] : 0); break;
1240:       case 1: color =            stag->rank[0]       + (dim > 2 ? stag->nRanks[0]*stag->rank[2] : 0); break;
1241:       case 2: color =            stag->rank[0]       +            stag->nRanks[0]*stag->rank[1]     ; break;
1242:       default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1243:     }
1244:     MPI_Comm_split(PetscObjectComm((PetscObject)dm),color,key,&subcomm);

1246:     /* Create sub-DMs living on these new communicators (which are destroyed by DMProduct) */
1247:     DMStagCreate1d(subcomm,stag->boundaryType[d],stag->N[d],dof0,dof1,stag->stencilType,stag->stencilWidth,stag->l[d],&subdm);
1248:     DMSetUp(subdm);
1249:     switch (d) {
1250:       case 0:
1251:         DMStagSetUniformCoordinatesExplicit(subdm,xmin,xmax,0,0,0,0);
1252:         break;
1253:       case 1:
1254:         DMStagSetUniformCoordinatesExplicit(subdm,ymin,ymax,0,0,0,0);
1255:         break;
1256:       case 2:
1257:         DMStagSetUniformCoordinatesExplicit(subdm,zmin,zmax,0,0,0,0);
1258:         break;
1259:       default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1260:     }
1261:     DMProductSetDM(dmc,d,subdm);
1262:     DMProductSetDimensionIndex(dmc,d,0);
1263:     DMDestroy(&subdm);
1264:     MPI_Comm_free(&subcomm);
1265:   }
1266:   return(0);
1267: }

1269: /*@C
1270:   DMStagVecGetArrayDOF - get access to raw local array

1272:   Logically Collective

1274:   Input Parameters:
1275: + dm - the DMStag object
1276: - vec - the Vec object

1278:   Output Parameters:
1279: . array - the array

1281:   Notes:
1282:   Indexing is array[k][j][i][idx].
1283:   Obtain idx with DMStagGetLocationSlot().

1285:   Level: beginner

1287: .seealso: DMSTAG, DMStagVecGetArrayDOFRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArrayDOF()
1288: @*/
1289: PetscErrorCode DMStagVecGetArrayDOF(DM dm,Vec vec,void *array)
1290: {
1291:   PetscErrorCode  ierr;
1292:   DM_Stag * const stag = (DM_Stag*)dm->data;
1293:   PetscInt        dim;
1294:   PetscInt        nLocal;

1299:   DMGetDimension(dm,&dim);
1300:   VecGetLocalSize(vec,&nLocal);
1301:   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);
1302:   switch (dim) {
1303:     case 1:
1304:       VecGetArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1305:       break;
1306:     case 2:
1307:       VecGetArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1308:       break;
1309:     case 3:
1310:       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);
1311:       break;
1312:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1313:   }
1314:   return(0);
1315: }

1317: /*@C
1318:   DMStagVecGetArrayDOFRead - get read-only access to raw local array

1320:   Logically Collective

1322:   Input Parameters:
1323: + dm - the DMStag object
1324: - vec - the Vec object

1326:   Output Parameters:
1327: . array - read-only the array

1329:   Notes:
1330:   Indexing is array[k][j][i][idx].
1331:   Obtain idx with DMStagGetLocationSlot()

1333:   Level: beginner

1335: .seealso: DMSTAG, DMStagVecGetArrayDOFRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArrayDOFRead()
1336: @*/
1337: PetscErrorCode DMStagVecGetArrayDOFRead(DM dm,Vec vec,void *array)
1338: {
1339:   PetscErrorCode  ierr;
1340:   DM_Stag * const stag = (DM_Stag*)dm->data;
1341:   PetscInt        dim;
1342:   PetscInt        nLocal;

1347:   DMGetDimension(dm,&dim);
1348:   VecGetLocalSize(vec,&nLocal);
1349:   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);
1350:   switch (dim) {
1351:     case 1:
1352:       VecGetArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1353:       break;
1354:     case 2:
1355:       VecGetArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1356:       break;
1357:     case 3:
1358:       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);
1359:       break;
1360:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1361:   }
1362:   return(0);
1363: }

1365: /*@C
1366:   DMStagVecRestoreArrayDOF - restore read-only access to a raw array

1368:   Logically Collective

1370:   Input Parameters:
1371: + dm - the DMStag object
1372: - vec - the Vec object

1374:   Output Parameters:
1375: . array - the array

1377:   Level: beginner

1379: .seealso: DMSTAG, DMStagVecGetArrayDOF(), DMDAVecRestoreArrayDOFRead()
1380: @*/
1381: PetscErrorCode DMStagVecRestoreArrayDOF(DM dm,Vec vec,void *array)
1382: {
1383:   PetscErrorCode  ierr;
1384:   DM_Stag * const stag = (DM_Stag*)dm->data;
1385:   PetscInt        dim;
1386:   PetscInt        nLocal;

1391:   DMGetDimension(dm,&dim);
1392:   VecGetLocalSize(vec,&nLocal);
1393:   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);
1394:   switch (dim) {
1395:     case 1:
1396:       VecRestoreArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1397:       break;
1398:     case 2:
1399:       VecRestoreArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1400:       break;
1401:     case 3:
1402:       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);
1403:       break;
1404:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1405:   }
1406:   return(0);
1407: }

1409: /*@C
1410:   DMStagVecRestoreArrayDOFRead - restore read-only access to a raw array

1412:   Logically Collective

1414:   Input Parameters:
1415: + dm - the DMStag object
1416: - vec - the Vec object

1418:   Output Parameters:
1419: . array - the read-only array

1421:   Level: beginner

1423: .seealso: DMSTAG, DMStagVecGetArrayDOFRead(), DMDAVecRestoreArrayDOFRead()
1424: @*/
1425: PetscErrorCode DMStagVecRestoreArrayDOFRead(DM dm,Vec vec,void *array)
1426: {
1427:   PetscErrorCode  ierr;
1428:   DM_Stag * const stag = (DM_Stag*)dm->data;
1429:   PetscInt        dim;
1430:   PetscInt        nLocal;

1435:   DMGetDimension(dm,&dim);
1436:   VecGetLocalSize(vec,&nLocal);
1437:   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);
1438:   switch (dim) {
1439:     case 1:
1440:       VecRestoreArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1441:       break;
1442:     case 2:
1443:       VecRestoreArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1444:       break;
1445:     case 3:
1446:       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);
1447:       break;
1448:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1449:   }
1450:   return(0);
1451: }