Actual source code: stagutils.c

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

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

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

108:   Not Collective

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

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

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

121:   Level: intermediate

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

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

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

173:   Not Collective

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

178:   Output Parameters:
179: + x,y,z - starting element indices in each direction
180: . m,n,p - element widths in each direction
181: - nExtrax,nExtray,nExtraz - number of extra partial elements in each direction. The number is 1 on the right, top, and front boundaries of the grid, otherwise 0.

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

186:   Level: beginner

188: .seealso: DMSTAG, DMStagGetGhostCorners()
189: @*/
190: PetscErrorCode DMStagGetCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *nExtrax,PetscInt *nExtray,PetscInt *nExtraz)
191: {
192:   const DM_Stag * const stag = (DM_Stag*)dm->data;

196:   if (x) *x = stag->start[0];
197:   if (y) *y = stag->start[1];
198:   if (z) *z = stag->start[2];
199:   if (m) *m = stag->n[0];
200:   if (n) *n = stag->n[1];
201:   if (p) *p = stag->n[2];
202:   if (nExtrax) *nExtrax = stag->lastRank[0] ? 1 : 0;
203:   if (nExtray) *nExtray = stag->lastRank[1] ? 1 : 0;
204:   if (nExtraz) *nExtraz = stag->lastRank[2] ? 1 : 0;
205:   return(0);
206: }

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

211:   Not Collective

213:   Input Parameter:
214: . dm - the DMStag object

216:   Output Parameters:
217: + dof0 - the number of points per 0-cell (vertex/node)
218: . dof1 - the number of points per 1-cell (elementin 1D, edge in 2D and 3D)
219: . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
220: - dof3 - the number of points per 3-cell (elementin 3D)

222:   Level: beginner

224: .seealso: DMSTAG, DMStagGetCorners(), DMStagGetGhostCorners(), DMStagGetGlobalSizes(), DMStagGetStencilWidth(), DMStagGetBoundaryTypes()
225: @*/
226: PetscErrorCode DMStagGetDOF(DM dm,PetscInt *dof0,PetscInt *dof1,PetscInt *dof2,PetscInt *dof3)
227: {
228:   const DM_Stag * const stag = (DM_Stag*)dm->data;

232:   if (dof0) *dof0 = stag->dof[0];
233:   if (dof1) *dof1 = stag->dof[1];
234:   if (dof2) *dof2 = stag->dof[2];
235:   if (dof3) *dof3 = stag->dof[3];
236:   return(0);
237: }

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

242:   Not Collective

244:   Input Argument:
245: + dm - the DMStag object

247:   Output Arguments:
248: + x,y,z - starting element indices in each direction
249: - m,n,p - element widths in each direction

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

254:   Level: beginner

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

264:   if (x) *x = stag->startGhost[0];
265:   if (y) *y = stag->startGhost[1];
266:   if (z) *z = stag->startGhost[2];
267:   if (m) *m = stag->nGhost[0];
268:   if (n) *n = stag->nGhost[1];
269:   if (p) *p = stag->nGhost[2];
270:   return(0);
271: }

273: /*@C
274:   DMStagGetGlobalSizes - get global element counts

276:   Not Collective

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

281:   Output Parameters:
282: . M,N,P - global element counts in each direction

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

287:   Level: beginner

289: .seealso: DMSTAG, DMStagGetLocalSizes()
290: @*/
291: PetscErrorCode DMStagGetGlobalSizes(DM dm,PetscInt* M,PetscInt* N,PetscInt* P)
292: {
293:   const DM_Stag * const stag = (DM_Stag*)dm->data;

297:   if (M) *M = stag->N[0];
298:   if (N) *N = stag->N[1];
299:   if (P) *P = stag->N[2];
300:   return(0);
301: }

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

306:   Not Collective

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

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

314:   Level: intermediate

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

319: .seealso: DMSTAG, DMStagGetIsLastRank()
320: @*/
321: PetscErrorCode DMStagGetIsFirstRank(DM dm,PetscBool *isFirstRank0,PetscBool *isFirstRank1,PetscBool *isFirstRank2)
322: {
323:   const DM_Stag * const stag = (DM_Stag*)dm->data;

327:   if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
328:   if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
329:   if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
330:   return(0);
331: }

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

336:   Not Collective

338:   Input Parameter:
339: . dm - the DMStag object

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

344:   Level: intermediate

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

350: .seealso: DMSTAG, DMStagGetIsFirstRank()
351: @*/
352: PetscErrorCode DMStagGetIsLastRank(DM dm,PetscBool *isLastRank0,PetscBool *isLastRank1,PetscBool *isLastRank2)
353: {
354:   const DM_Stag * const stag = (DM_Stag*)dm->data;

358:   if (isLastRank0) *isLastRank0 = stag->lastRank[0];
359:   if (isLastRank1) *isLastRank1 = stag->lastRank[1];
360:   if (isLastRank2) *isLastRank2 = stag->lastRank[2];
361:   return(0);
362: }

364: /*@C
365:   DMStagGetLocalSizes - get local elementwise sizes

367:   Not Collective

369:   Input Parameter:
370: . dm - the DMStag object

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

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

379:   Level: beginner

381: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetDOF(), DMStagGetNumRanks()
382: @*/
383: PetscErrorCode DMStagGetLocalSizes(DM dm,PetscInt* m,PetscInt* n,PetscInt* p)
384: {
385:   const DM_Stag * const stag = (DM_Stag*)dm->data;

389:   if (m) *m = stag->n[0];
390:   if (n) *n = stag->n[1];
391:   if (p) *p = stag->n[2];
392:   return(0);
393: }

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

398:   Not Collective

400:   Input Parameter:
401: . dm - the DMStag object

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

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

410:   Level: beginner

412: .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetLocalSize(), DMStagSetNumRank()
413: @*/
414: PetscErrorCode DMStagGetNumRanks(DM dm,PetscInt *nRanks0,PetscInt *nRanks1,PetscInt *nRanks2)
415: {
416:   const DM_Stag * const stag = (DM_Stag*)dm->data;

420:   if (nRanks0) *nRanks0 = stag->nRanks[0];
421:   if (nRanks1) *nRanks1 = stag->nRanks[1];
422:   if (nRanks2) *nRanks2 = stag->nRanks[2];
423:   return(0);
424: }

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

429:   Not Collective

431:   Input Parameter:
432: . dm - the DMStag object

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

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

441:   Level: developer

443: .seealso: DMSTAG, DMStagGetDOF()
444: @*/
445: PetscErrorCode DMStagGetEntriesPerElement(DM dm,PetscInt *entriesPerElement)
446: {
447:   const DM_Stag * const stag = (DM_Stag*)dm->data;

451:   if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
452:   return(0);
453: }

455: /*@C
456:   DMStagGetStencilType - get elementwise ghost/halo stencil type

458:   Not Collective

460:   Input Parameter:
461: . dm - the DMStag object

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

466:   Level: beginner

468: .seealso: DMSTAG, DMStagSetStencilType(), DMStagStencilType
469: @*/
470: PetscErrorCode DMStagGetStencilType(DM dm,DMStagStencilType *stencilType)
471: {
472:   DM_Stag * const stag = (DM_Stag*)dm->data;

476:   *stencilType = stag->stencilType;
477:   return(0);
478: }

480: /*@C
481:   DMStagGetStencilWidth - get elementwise stencil width

483:   Not Collective

485:   Input Parameter:
486: . dm - the DMStag object

488:   Output Parameters:
489: . stencilWidth - stencil/halo/ghost width in elements

491:   Level: beginner

493: .seealso: DMSTAG
494: @*/
495: PetscErrorCode DMStagGetStencilWidth(DM dm,PetscInt *stencilWidth)
496: {
497:   const DM_Stag * const stag = (DM_Stag*)dm->data;

501:   if (stencilWidth) *stencilWidth = stag->stencilWidth;
502:   return(0);
503: }

505: /*@C
506:   DMStagGetOwnershipRanges - get elements per rank in each direction

508:   Not Collective

510:   Input Parameter:
511: .     dm - the DMStag object

513:   Output Parameters:
514: +     lx - ownership along x direction (optional)
515: .     ly - ownership along y direction (optional)
516: -     lz - ownership along z direction (optional)

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

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

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

526:   Level: intermediate

528: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagSetOwnershipRanges(), DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d()
529: @*/
530: PetscErrorCode DMStagGetOwnershipRanges(DM dm,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
531: {
532:   const DM_Stag * const stag = (DM_Stag*)dm->data;

536:   if (lx) *lx = stag->l[0];
537:   if (ly) *ly = stag->l[1];
538:   if (lz) *lz = stag->l[2];
539:   return(0);
540: }

542: /*@C
543:   DMStagCreateCompatibleDMStag - create a compatible DMStag with different dof/stratum

545:   Collective

547:   Input Parameters
548: + dm - the DMStag object
549: - dof0,dof1,dof2,dof3 - number of dof on each stratum in the new DMStag

551:   Output Parameters:
552: . newdm - the new, compatible DMStag

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

558:   Level: intermediate

560: .seealso: DMSTAG, DMDACreateCompatibleDMDA(), DMGetCompatibility(), DMStagMigrateVec()
561: @*/
562: PetscErrorCode DMStagCreateCompatibleDMStag(DM dm,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DM *newdm)
563: {
564:   PetscErrorCode  ierr;
565:   const DM_Stag * const stag = (DM_Stag*)dm->data;
566:   PetscInt        dim;

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

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

590:   Not Collective

592:   Input Parameters:
593: + dm - the DMStag object
594: . loc - location relative to an element
595: - c - component

597:   Output Parameter:
598: . slot - index to use

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

605:   Level: beginner

607: .seealso: DMSTAG, DMStagVecGetArrayDOF(), DMStagVecGetArrayDOFRead(), DMStagGetDOF(), DStagMGetEntriesPerElement()
608: @*/
609: PetscErrorCode DMStagGetLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt c,PetscInt *slot)
610: {
611:   DM_Stag * const stag = (DM_Stag*)dm->data;

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

628: /*@C
629:   DMStagMigrateVec - transfer a vector associated with a DMStag to a vector associated with a compatible DMStag

631:   Collective

633:   Input Parameters:
634: + dm - the source DMStag object
635: . vec - the source vector, compatible with dm
636: . dmTo - the compatible destination DMStag object
637: - vecTo - the destination vector, compatible with dmTo

639:   Notes:
640:   Extra dof are ignored, and unfilled dof are zeroed.
641:   Currently only implemented to migrate global vectors to global vectors.

643:   Level: advanced

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

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

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

748: /*@C
749:   DMStagRestore1dCoordinateArraysDOFRead - restore local array access

751:   Logically Collective

753:   Input Parameter:
754: . dm - the DMStag object

756:   Output Parameters:
757: . arrX,arrY,arrX - local 1D coordinate arrays

759:   Level: intermediate

761: .seealso: DMSTAG, DMStagGet1dCoordinateArraysDOFRead()
762: @*/
763: PetscErrorCode DMStagRestore1dCoordinateArraysDOFRead(DM dm,void *arrX,void *arrY,void *arrZ)
764: {
765:   PetscErrorCode  ierr;
766:   PetscInt        dim,d;
767:   void*           arr[DMSTAG_MAX_DIM];
768:   DM              dmCoord;

772:   DMGetDimension(dm,&dim);
773:   arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
774:   DMGetCoordinateDM(dm,&dmCoord);
775:   for (d=0; d<dim; ++d) {
776:     DM        subDM;
777:     Vec       coord1d;
778:     DMProductGetDM(dmCoord,d,&subDM);
779:     DMGetCoordinatesLocal(subDM,&coord1d);
780:     DMStagVecRestoreArrayDOFRead(subDM,coord1d,arr[d]);
781:   }
782:   return(0);
783: }

785: /*@C
786:   DMStagSetBoundaryTypes - set DMStag boundary types

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

790:   Input Parameters:
791: + dm - the DMStag object
792: - boundaryType0,boundaryType1,boundaryType2 - boundary types in each direction

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

797:   Level: advanced

799: .seealso: DMSTAG, DMBoundaryType, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d()
800: @*/
801: PetscErrorCode DMStagSetBoundaryTypes(DM dm,DMBoundaryType boundaryType0,DMBoundaryType boundaryType1,DMBoundaryType boundaryType2)
802: {
803:   PetscErrorCode  ierr;
804:   DM_Stag * const stag  = (DM_Stag*)dm->data;
805:   PetscInt        dim;

812:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
813:   DMGetDimension(dm,&dim);
814:   if (boundaryType0           ) stag->boundaryType[0] = boundaryType0;
815:   if (boundaryType1 && dim > 1) stag->boundaryType[1] = boundaryType1;
816:   if (boundaryType2 && dim > 2) stag->boundaryType[2] = boundaryType2;
817:   return(0);
818: }

820: /*@C
821:   DMStagSetCoordinateDMType - set DM type to store coordinates

823:   Logically Collective; dmtype must contain common value

825:   Input Parameters:
826: + dm - the DMStag object
827: - dmtype - DMtype for coordinates, either DMSTAG or DMPRODUCT

829:   Level: advanced

831: .seealso: DMSTAG, DMPRODUCT, DMGetCoordinateDM(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMType
832: @*/
833: PetscErrorCode DMStagSetCoordinateDMType(DM dm,DMType dmtype)
834: {
835:   DM_Stag * const stag = (DM_Stag*)dm->data;

839:   stag->coordinateDMType = dmtype;
840:   return(0);
841: }

843: /*@C
844:   DMStagSetDOF - set dof/stratum

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

848:   Input Parameters:
849: + dm - the DMStag object
850: - dof0,dof1,dof2,dof3 - dof per stratum

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

855:   Level: advanced

857: .seealso: DMSTAG
858: @*/
859: PetscErrorCode DMStagSetDOF(DM dm,PetscInt dof0, PetscInt dof1,PetscInt dof2,PetscInt dof3)
860: {
861:   PetscErrorCode  ierr;
862:   DM_Stag * const stag = (DM_Stag*)dm->data;
863:   PetscInt        dim;

871:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
872:   DMGetDimension(dm,&dim);
873:   if (dof0 < 0)            SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof0 cannot be negative");
874:   if (dof1 < 0)            SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof1 cannot be negative");
875:   if (dim > 1 && dof2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof2 cannot be negative");
876:   if (dim > 2 && dof3 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof3 cannot be negative");
877:                stag->dof[0] = dof0;
878:                stag->dof[1] = dof1;
879:   if (dim > 1) stag->dof[2] = dof2;
880:   if (dim > 2) stag->dof[3] = dof3;
881:   return(0);
882: }

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

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

889:   Input Parameters:
890: + dm - the DMStag object
891: - nRanks0,nRanks1,nRanks2 - number of ranks in each direction

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

896:   Level: developer

898: .seealso: DMSTAG
899: @*/
900: PetscErrorCode DMStagSetNumRanks(DM dm,PetscInt nRanks0,PetscInt nRanks1,PetscInt nRanks2)
901: {
902:   PetscErrorCode  ierr;
903:   DM_Stag * const stag = (DM_Stag*)dm->data;
904:   PetscInt        dim;

911:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
912:   DMGetDimension(dm,&dim);
913:   if (nRanks0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in X direction cannot be less than 1");
914:   if (dim > 1 && nRanks1 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in Y direction cannot be less than 1");
915:   if (dim > 2 && nRanks2 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in Z direction cannot be less than 1");
916:   if (nRanks0) stag->nRanks[0] = nRanks0;
917:   if (nRanks1) stag->nRanks[1] = nRanks1;
918:   if (nRanks2) stag->nRanks[2] = nRanks2;
919:   return(0);
920: }

922: /*@C
923:   DMStagSetStencilType - set elementwise ghost/halo stencil type

925:   Logically Collective; stencilType must contain common value

927:   Input Parameters:
928: + dm - the DMStag object
929: - stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE

931:   Level: beginner

933: .seealso: DMSTAG, DMStagGetStencilType(), DMStagStencilType
934: @*/
935: PetscErrorCode DMStagSetStencilType(DM dm,DMStagStencilType stencilType)
936: {
937:   DM_Stag * const stag = (DM_Stag*)dm->data;

942:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
943:   stag->stencilType = stencilType;
944:   return(0);
945: }

947: /*@C
948:   DMStagSetStencilWidth - set elementwise stencil width

950:   Logically Collective; stencilWidth must contain common value

952:   Input Parameters:
953: + dm - the DMStag object
954: - stencilWidth - stencil/halo/ghost width in elements

956:   Level: beginner

958: .seealso: DMSTAG
959: @*/
960: PetscErrorCode DMStagSetStencilWidth(DM dm,PetscInt stencilWidth)
961: {
962:   DM_Stag * const stag = (DM_Stag*)dm->data;

967:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
968:   if (stencilWidth < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Stencil width must be non-negative");
969:   stag->stencilWidth = stencilWidth;
970:   return(0);
971: }

973: /*@C
974:   DMStagSetGlobalSizes - set global element counts in each direction

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

978:   Input Parameters:
979: + dm - the DMStag object
980: - N0,N1,N2 - global elementwise sizes

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

985:   Level: advanced

987: .seealso: DMSTAG, DMStagGetGlobalSizes()
988: @*/
989: PetscErrorCode DMStagSetGlobalSizes(DM dm,PetscInt N0,PetscInt N1,PetscInt N2)
990: {
991:   PetscErrorCode  ierr;
992:   DM_Stag * const stag = (DM_Stag*)dm->data;
993:   PetscInt        dim;

997:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
998:   DMGetDimension(dm,&dim);
999:   if (N0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in X direction must be positive");
1000:   if (dim > 1 && N1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Y direction must be positive");
1001:   if (dim > 2 && N2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Z direction must be positive");
1002:   if (N0) stag->N[0] = N0;
1003:   if (N1) stag->N[1] = N1;
1004:   if (N2) stag->N[2] = N2;
1005:   return(0);
1006: }

1008: /*@C
1009:   DMStagSetOwnershipRanges - set elements per rank in each direction

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

1013:   Input Parameters:
1014: + dm - the DMStag object
1015: - lx,ly,lz - element counts for each rank in each direction

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

1020:   Level: developer

1022: .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagGetOwnershipRanges()
1023: @*/
1024: PetscErrorCode DMStagSetOwnershipRanges(DM dm,PetscInt const *lx,PetscInt const *ly,PetscInt const *lz)
1025: {
1026:   PetscErrorCode  ierr;
1027:   DM_Stag * const stag = (DM_Stag*)dm->data;
1028:   const PetscInt  *lin[3];
1029:   PetscInt        d;

1033:   if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1034:   lin[0] = lx; lin[1] = ly; lin[2] = lz;
1035:   for (d=0; d<3; ++d) {
1036:     if (lin[d]) {
1037:       if (stag->nRanks[d] < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of ranks");
1038:       if (!stag->l[d]) {
1039:         PetscMalloc1(stag->nRanks[d], &stag->l[d]);
1040:       }
1041:       PetscMemcpy(stag->l[d], lin[d], stag->nRanks[d]*sizeof(PetscInt));
1042:     }
1043:   }
1044:   return(0);
1045: }

1047: /*@C
1048:   DMStagSetUniformCoordinates - set DMStag coordinates to be a uniform grid

1050:   Collective

1052:   Input Parameters:
1053: + dm - the DMStag object
1054: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

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

1060:   Level: advanced

1062: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType(), DMGetCoordinateDM(), DMGetCoordinates()
1063: @*/
1064: PetscErrorCode DMStagSetUniformCoordinates(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1065: {
1066:   PetscErrorCode  ierr;
1067:   DM_Stag * const stag = (DM_Stag*)dm->data;
1068:   PetscBool       flg_stag,flg_product;

1072:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1073:   if (!stag->coordinateDMType) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"You must first call DMStagSetCoordinateDMType()");
1074:   PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg_stag);
1075:   PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg_product);
1076:   if (flg_stag) {
1077:     DMStagSetUniformCoordinatesExplicit(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1078:   } else if (flg_product) {
1079:     DMStagSetUniformCoordinatesProduct(dm,xmin,xmax,ymin,ymax,zmin,zmax);
1080:   } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported DM Type %s",stag->coordinateDMType);
1081:   return(0);
1082: }

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

1087:   Collective

1089:   Input Parameters:
1090: + dm - the DMStag object
1091: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

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

1098:   Level: beginner

1100: .seealso: DMSTAG, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType()
1101: @*/
1102: PetscErrorCode DMStagSetUniformCoordinatesExplicit(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1103: {
1104:   PetscErrorCode  ierr;
1105:   DM_Stag * const stag = (DM_Stag*)dm->data;
1106:   PetscInt        dim;
1107:   PetscBool       flg;

1111:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1112:   PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg);
1113:   if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1114:   DMStagSetCoordinateDMType(dm,DMSTAG);
1115:   DMGetDimension(dm,&dim);
1116:   switch (dim) {
1117:     case 1: DMStagSetUniformCoordinatesExplicit_1d(dm,xmin,xmax);                     break;
1118:     case 2: DMStagSetUniformCoordinatesExplicit_2d(dm,xmin,xmax,ymin,ymax);           break;
1119:     case 3: DMStagSetUniformCoordinatesExplicit_3d(dm,xmin,xmax,ymin,ymax,zmin,zmax); break;
1120:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1121:   }
1122:   return(0);
1123: }

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

1128:   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
1129:   coordinates.

1131:   Collective

1133:   Input Parameters:
1134: + dm - the DMStag object
1135: - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values

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

1140:   Level: intermediate

1142: .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetCoordinateDMType()
1143: @*/
1144: PetscErrorCode DMStagSetUniformCoordinatesProduct(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1145: {
1146:   PetscErrorCode  ierr;
1147:   DM_Stag * const stag = (DM_Stag*)dm->data;
1148:   DM              dmc;
1149:   PetscInt        dim,d,dof0,dof1;
1150:   PetscBool       flg;

1154:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1155:   PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg);
1156:   if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1157:   DMStagSetCoordinateDMType(dm,DMPRODUCT);
1158:   DMGetCoordinateDM(dm,&dmc);
1159:   DMGetDimension(dm,&dim);

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

1163:   dof0 = 0;
1164:   for (d=0; d<dim; ++d) {          /* Include point dof in the sub-DMs if any non-element strata are active in the original DMStag */
1165:     if (stag->dof[d]) {
1166:       dof0 = 1;
1167:       break;
1168:     }
1169:   }
1170:   dof1 = stag->dof[dim] ? 1 : 0; /*  Include element dof in the sub-DMs if the elements are active in the original DMStag */

1172:   for (d=0; d<dim; ++d) {
1173:     DM                subdm;
1174:     MPI_Comm          subcomm;
1175:     PetscMPIInt       color;
1176:     const PetscMPIInt key = 0; /* let existing rank break ties */

1178:     /* Choose colors based on position in the plane orthogonal to this dim, and split */
1179:     switch (d) {
1180:       case 0: color = (dim > 1 ? stag->rank[1] : 0)  + (dim > 2 ? stag->nRanks[1]*stag->rank[2] : 0); break;
1181:       case 1: color =            stag->rank[0]       + (dim > 2 ? stag->nRanks[0]*stag->rank[2] : 0); break;
1182:       case 2: color =            stag->rank[0]       +            stag->nRanks[0]*stag->rank[1]     ; break;
1183:       default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1184:     }
1185:     MPI_Comm_split(PetscObjectComm((PetscObject)dm),color,key,&subcomm);

1187:     /* Create sub-DMs living on these new communicators (which are destroyed by DMProduct) */
1188:     DMStagCreate1d(subcomm,stag->boundaryType[d],stag->N[d],dof0,dof1,stag->stencilType,stag->stencilWidth,stag->l[d],&subdm);
1189:     DMSetUp(subdm);
1190:     switch (d) {
1191:       case 0:
1192:         DMStagSetUniformCoordinatesExplicit(subdm,xmin,xmax,0,0,0,0);
1193:         break;
1194:       case 1:
1195:         DMStagSetUniformCoordinatesExplicit(subdm,ymin,ymax,0,0,0,0);
1196:         break;
1197:       case 2:
1198:         DMStagSetUniformCoordinatesExplicit(subdm,zmin,zmax,0,0,0,0);
1199:         break;
1200:       default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1201:     }
1202:     DMProductSetDM(dmc,d,subdm);
1203:     DMProductSetDimensionIndex(dmc,d,0);
1204:     DMDestroy(&subdm);
1205:     MPI_Comm_free(&subcomm);
1206:   }
1207:   return(0);
1208: }

1210: /*@C
1211:   DMStagVecGetArrayDOF - get access to raw local array

1213:   Logically Collective

1215:   Input Parameters:
1216: + dm - the DMStag object
1217: - vec - the Vec object

1219:   Output Parameters:
1220: . array - the array

1222:   Notes:
1223:   Indexing is array[k][j][i][idx].
1224:   Obtain idx with DMStagGetLocationSlot().

1226:   Level: beginner

1228: .seealso: DMSTAG, DMStagVecGetArrayDOFRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector()
1229: @*/
1230: PetscErrorCode DMStagVecGetArrayDOF(DM dm,Vec vec,void *array)
1231: {
1232:   PetscErrorCode  ierr;
1233:   DM_Stag * const stag = (DM_Stag*)dm->data;
1234:   PetscInt        dim;
1235:   PetscInt        nLocal;

1240:   DMGetDimension(dm,&dim);
1241:   VecGetLocalSize(vec,&nLocal);
1242:   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);
1243:   switch (dim) {
1244:     case 1:
1245:       VecGetArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1246:       break;
1247:     case 2:
1248:       VecGetArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1249:       break;
1250:     case 3:
1251:       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);
1252:       break;
1253:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1254:   }
1255:   return(0);
1256: }

1258: /*@C
1259:   DMStagVecGetArrayDOFRead - get read-only access to raw local array

1261:   Logically Collective

1263:   Input Parameters:
1264: + dm - the DMStag object
1265: - vec - the Vec object

1267:   Output Parameters:
1268: . array - read-only the array

1270:   Notes:
1271:   Indexing is array[k][j][i][idx].
1272:   Obtain idx with DMStagGetLocationSlot()

1274:   Level: beginner

1276: .seealso: DMSTAG, DMStagVecGetArrayDOFRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector()
1277: @*/
1278: PetscErrorCode DMStagVecGetArrayDOFRead(DM dm,Vec vec,void *array)
1279: {
1280:   PetscErrorCode  ierr;
1281:   DM_Stag * const stag = (DM_Stag*)dm->data;
1282:   PetscInt        dim;
1283:   PetscInt        nLocal;

1288:   DMGetDimension(dm,&dim);
1289:   VecGetLocalSize(vec,&nLocal);
1290:   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);
1291:   switch (dim) {
1292:     case 1:
1293:       VecGetArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1294:       break;
1295:     case 2:
1296:       VecGetArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1297:       break;
1298:     case 3:
1299:       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);
1300:       break;
1301:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1302:   }
1303:   return(0);
1304: }

1306: /*@C
1307:   DMStagVecRestoreArrayDOF - restore read-only access to a raw array

1309:   Logically Collective

1311:   Input Parameters:
1312: + dm - the DMStag object
1313: - vec - the Vec object

1315:   Output Parameters:
1316: . array - the array

1318:   Level: beginner

1320: .seealso: DMSTAG, DMStagVecGetArrayDOF()
1321: @*/
1322: PetscErrorCode DMStagVecRestoreArrayDOF(DM dm,Vec vec,void *array)
1323: {
1324:   PetscErrorCode  ierr;
1325:   DM_Stag * const stag = (DM_Stag*)dm->data;
1326:   PetscInt        dim;
1327:   PetscInt        nLocal;

1332:   DMGetDimension(dm,&dim);
1333:   VecGetLocalSize(vec,&nLocal);
1334:   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);
1335:   switch (dim) {
1336:     case 1:
1337:       VecRestoreArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1338:       break;
1339:     case 2:
1340:       VecRestoreArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1341:       break;
1342:     case 3:
1343:       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);
1344:       break;
1345:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1346:   }
1347:   return(0);
1348: }

1350: /*@C
1351:   DMStagVecRestoreArrayDOFRead - restore read-only access to a raw array

1353:   Logically Collective

1355:   Input Parameters:
1356: + dm - the DMStag object
1357: - vec - the Vec object

1359:   Output Parameters:
1360: . array - the read-only array

1362:   Level: beginner

1364: .seealso: DMSTAG, DMStagVecGetArrayDOFRead()
1365: @*/
1366: PetscErrorCode DMStagVecRestoreArrayDOFRead(DM dm,Vec vec,void *array)
1367: {
1368:   PetscErrorCode  ierr;
1369:   DM_Stag * const stag = (DM_Stag*)dm->data;
1370:   PetscInt        dim;
1371:   PetscInt        nLocal;

1376:   DMGetDimension(dm,&dim);
1377:   VecGetLocalSize(vec,&nLocal);
1378:   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);
1379:   switch (dim) {
1380:     case 1:
1381:       VecRestoreArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);
1382:       break;
1383:     case 2:
1384:       VecRestoreArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);
1385:       break;
1386:     case 3:
1387:       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);
1388:       break;
1389:     default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1390:   }
1391:   return(0);
1392: }