Actual source code: stag.c

petsc-master 2020-02-20
Report Typos and Errors
  1: /*
  2:    Implementation of DMStag, defining dimension-independent functions in the
  3:    DM API. stag1d.c, stag2d.c, and stag3d.c may include dimension-specific
  4:    implementations of DM API functions, and other files here contain additional
  5:    DMStag-specific API functions, as well as internal functions.
  6: */
  7:  #include <petsc/private/dmstagimpl.h>
  8:  #include <petscsf.h>

 10: static PetscErrorCode DMClone_Stag(DM dm,DM *newdm)
 11: {

 15:   /* Destroy the DM created by generic logic in DMClone() */
 16:   if (*newdm) {
 17:     DMDestroy(newdm);
 18:   }
 19:   DMStagDuplicateWithoutSetup(dm,PetscObjectComm((PetscObject)dm),newdm);
 20:   DMSetUp(*newdm);
 21:   return(0);
 22: }

 24: static PetscErrorCode DMDestroy_Stag(DM dm)
 25: {
 27:   DM_Stag        *stag;
 28:   PetscInt       i;

 31:   stag = (DM_Stag*)dm->data;
 32:   for (i=0; i<DMSTAG_MAX_DIM; ++i) {
 33:     PetscFree(stag->l[i]);
 34:   }
 35:   VecScatterDestroy(&stag->gtol);
 36:   VecScatterDestroy(&stag->ltog_injective);
 37:   PetscFree(stag->neighbors);
 38:   PetscFree(stag->locationOffsets);
 39:   PetscFree(stag->coordinateDMType);
 40:   PetscFree(stag);
 41:   return(0);
 42: }

 44: static PetscErrorCode DMCreateGlobalVector_Stag(DM dm,Vec *vec)
 45: {
 46:   PetscErrorCode  ierr;
 47:   DM_Stag * const stag = (DM_Stag*)dm->data;

 50:   VecCreateMPI(PetscObjectComm((PetscObject)dm),stag->entries,PETSC_DECIDE,vec);
 51:   VecSetDM(*vec,dm);
 52:   /* Could set some ops, as DMDA does */
 53:   VecSetLocalToGlobalMapping(*vec,dm->ltogmap);
 54:   return(0);
 55: }

 57: static PetscErrorCode DMCreateLocalVector_Stag(DM dm,Vec *vec)
 58: {
 59:   PetscErrorCode  ierr;
 60:   DM_Stag * const stag = (DM_Stag*)dm->data;

 63:   VecCreateSeq(PETSC_COMM_SELF,stag->entriesGhost,vec);
 64:   VecSetBlockSize(*vec,stag->entriesPerElement);
 65:   VecSetDM(*vec,dm);
 66:   return(0);
 67: }

 69: static PetscErrorCode DMCreateMatrix_Stag(DM dm,Mat *mat)
 70: {
 71:   PetscErrorCode         ierr;
 72:   const DM_Stag * const  stag = (DM_Stag*)dm->data;
 73:   MatType                matType;
 74:   PetscBool              isaij,isshell;
 75:   PetscInt               width,nNeighbors,dim;
 76:   ISLocalToGlobalMapping ltogmap;

 79:   DMGetDimension(dm,&dim);
 80:   DMGetMatType(dm,&matType);
 81:   PetscStrcmp(matType,MATAIJ,&isaij);
 82:   PetscStrcmp(matType,MATSHELL,&isshell);

 84:   if (isaij) {
 85:     /* This implementation gives a very dense stencil, which is likely unsuitable for
 86:        real applications. */
 87:     switch (stag->stencilType) {
 88:       case DMSTAG_STENCIL_NONE:
 89:         nNeighbors = 1;
 90:         break;
 91:       case DMSTAG_STENCIL_STAR:
 92:         switch (dim) {
 93:           case 1 :
 94:             nNeighbors = 2*stag->stencilWidth + 1;
 95:             break;
 96:           case 2 :
 97:             nNeighbors = 4*stag->stencilWidth + 3;
 98:             break;
 99:           case 3 :
100:             nNeighbors = 6*stag->stencilWidth + 5;
101:             break;
102:           default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %d",dim);
103:         }
104:         break;
105:       case DMSTAG_STENCIL_BOX:
106:         switch (dim) {
107:           case 1 :
108:             nNeighbors = (2*stag->stencilWidth + 1);
109:             break;
110:           case 2 :
111:             nNeighbors = (2*stag->stencilWidth + 1) * (2*stag->stencilWidth + 1);
112:             break;
113:           case 3 :
114:             nNeighbors = (2*stag->stencilWidth + 1) * (2*stag->stencilWidth + 1) * (2*stag->stencilWidth + 1);
115:             break;
116:           default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %d",dim);
117:         }
118:         break;
119:       default : SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported stencil");
120:     }
121:     width = (stag->dof[0] + stag->dof[1] + stag->dof[2] + stag->dof[3]) * nNeighbors;
122:     MatCreateAIJ(PETSC_COMM_WORLD,stag->entries,stag->entries,PETSC_DETERMINE,PETSC_DETERMINE,width,NULL,width,NULL,mat);
123:   } else if (isshell) {
124:     MatCreate(PetscObjectComm((PetscObject)dm),mat);
125:     MatSetSizes(*mat,stag->entries,stag->entries,PETSC_DETERMINE,PETSC_DETERMINE);
126:     MatSetType(*mat,MATSHELL);
127:     MatSetUp(*mat);
128:   } else SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not implemented for Mattype %s",matType);

130:   DMGetLocalToGlobalMapping(dm,&ltogmap);
131:   MatSetLocalToGlobalMapping(*mat,ltogmap,ltogmap);
132:   MatSetDM(*mat,dm);
133:   return(0);
134: }

136: static PetscErrorCode DMGetCompatibility_Stag(DM dm,DM dm2,PetscBool *compatible,PetscBool *set)
137: {
138:   PetscErrorCode  ierr;
139:   const DM_Stag * const stag  = (DM_Stag*)dm->data;
140:   const DM_Stag * const stag2 = (DM_Stag*)dm2->data;
141:   PetscInt              dim,dim2,i;
142:   MPI_Comm              comm;
143:   PetscMPIInt           sameComm;
144:   DMType                type2;
145:   PetscBool             sameType;

148:   DMGetType(dm2,&type2);
149:   PetscStrcmp(DMSTAG,type2,&sameType);
150:   if (!sameType) {
151:     PetscInfo1((PetscObject)dm,"DMStag compatibility check not implemented with DM of type %s\n",type2);
152:     *set = PETSC_FALSE;
153:     return(0);
154:   }

156:   PetscObjectGetComm((PetscObject)dm,&comm);
157:   MPI_Comm_compare(comm,PetscObjectComm((PetscObject)dm2),&sameComm);
158:   if (sameComm != MPI_IDENT) {
159:     PetscInfo2((PetscObject)dm,"DMStag objects have different communicators: %d != %d\n",comm,PetscObjectComm((PetscObject)dm2));
160:     *set = PETSC_FALSE;
161:     return(0);
162:   }
163:   DMGetDimension(dm ,&dim );
164:   DMGetDimension(dm2,&dim2);
165:   if (dim != dim2) {
166:     PetscInfo((PetscObject)dm,"DMStag objects have different dimensions");
167:     *set = PETSC_TRUE;
168:     *compatible = PETSC_FALSE;
169:     return(0);
170:   }
171:   for (i=0; i<dim; ++i) {
172:     if (stag->N[i] != stag2->N[i]) {
173:       PetscInfo3((PetscObject)dm,"DMStag objects have different global numbers of elements in dimension %D: %D != %D\n",i,stag->n[i],stag2->n[i]);
174:       *set = PETSC_TRUE;
175:       *compatible = PETSC_FALSE;
176:       return(0);
177:     }
178:     if (stag->n[i] != stag2->n[i]) {
179:       PetscInfo3((PetscObject)dm,"DMStag objects have different local numbers of elements in dimension %D: %D != %D\n",i,stag->n[i],stag2->n[i]);
180:       *set = PETSC_TRUE;
181:       *compatible = PETSC_FALSE;
182:       return(0);
183:     }
184:     if (stag->boundaryType[i] != stag2->boundaryType[i]) {
185:       PetscInfo3((PetscObject)dm,"DMStag objects have different boundary types in dimension %d: %s != %s\n",i,stag->boundaryType[i],stag2->boundaryType[i]);
186:       *set = PETSC_TRUE;
187:       *compatible = PETSC_FALSE;
188:       return(0);
189:     }
190:   }
191:   /* Note: we include stencil type and width in the notion of compatibility, as this affects
192:      the "atlas" (local subdomains). This might be irritating in legitimate cases
193:      of wanting to transfer between two other-wise compatible DMs with different
194:      stencil characteristics. */
195:   if (stag->stencilType != stag2->stencilType) {
196:     PetscInfo2((PetscObject)dm,"DMStag objects have different ghost stencil types: %s != %s\n",DMStagStencilTypes[stag->stencilType],DMStagStencilTypes[stag2->stencilType]);
197:     *set = PETSC_TRUE;
198:     *compatible = PETSC_FALSE;
199:     return(0);
200:   }
201:   if (stag->stencilWidth != stag2->stencilWidth) {
202:     PetscInfo2((PetscObject)dm,"DMStag objects have different ghost stencil widths: %D != %D\n",stag->stencilWidth,stag->stencilWidth);
203:     *set = PETSC_TRUE;
204:     *compatible = PETSC_FALSE;
205:     return(0);
206:   }
207:   *set = PETSC_TRUE;
208:   *compatible = PETSC_TRUE;
209:   return(0);
210: }


213: /*
214: Note there are several orderings in play here.
215: In all cases, non-element dof are associated with the element that they are below/left/behind, and the order in 2D proceeds vertex/bottom edge/left edge/element (with all dof on each together).
216: Also in all cases, only subdomains which are the last in their dimension have partial elements.

218: 1) "Natural" Ordering (not used). Number adding each full or partial (on the right or top) element, starting at the bottom left (i=0,j=0) and proceeding across the entire domain, row by row to get a global numbering.
219: 2) Global ("PETSc") ordering. The same as natural, but restricted to each domain. So, traverse all elements (again starting at the bottom left and going row-by-row) on rank 0, then continue numbering with rank 1, and so on.
220: 3) Local ordering. Including ghost elements (both interior and on the right/top/front to complete partial elements), use the same convention to create a local numbering.
221: */

223: static PetscErrorCode DMLocalToGlobalBegin_Stag(DM dm,Vec l,InsertMode mode,Vec g)
224: {
225:   PetscErrorCode  ierr;
226:   DM_Stag * const stag = (DM_Stag*)dm->data;

229:   if (mode == ADD_VALUES) {
230:     VecScatterBegin(stag->gtol,l,g,mode,SCATTER_REVERSE);
231:   } else if (mode == INSERT_VALUES) {
232:     if (stag->ltog_injective) {
233:       VecScatterBegin(stag->ltog_injective,l,g,mode,SCATTER_FORWARD);
234:     } else {
235:       VecScatterBegin(stag->gtol,l,g,mode,SCATTER_REVERSE_LOCAL);
236:     }
237:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported InsertMode");
238:   return(0);
239: }

241: static PetscErrorCode DMLocalToGlobalEnd_Stag(DM dm,Vec l,InsertMode mode,Vec g)
242: {
243:   PetscErrorCode  ierr;
244:   DM_Stag * const stag = (DM_Stag*)dm->data;

247:   if (mode == ADD_VALUES) {
248:     VecScatterEnd(stag->gtol,l,g,mode,SCATTER_REVERSE);
249:   } else if (mode == INSERT_VALUES) {
250:     if (stag->ltog_injective) {
251:       VecScatterEnd(stag->ltog_injective,l,g,mode,SCATTER_FORWARD);
252:     } else {
253:       VecScatterEnd(stag->gtol,l,g,mode,SCATTER_REVERSE_LOCAL);
254:     }
255:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported InsertMode");
256:   return(0);
257: }

259: static PetscErrorCode DMGlobalToLocalBegin_Stag(DM dm,Vec g,InsertMode mode,Vec l)
260: {
261:   PetscErrorCode  ierr;
262:   DM_Stag * const stag = (DM_Stag*)dm->data;

265:   VecScatterBegin(stag->gtol,g,l,mode,SCATTER_FORWARD);
266:   return(0);
267: }

269: static PetscErrorCode DMGlobalToLocalEnd_Stag(DM dm,Vec g,InsertMode mode,Vec l)
270: {
271:   PetscErrorCode  ierr;
272:   DM_Stag * const stag = (DM_Stag*)dm->data;

275:   VecScatterEnd(stag->gtol,g,l,mode,SCATTER_FORWARD);
276:   return(0);
277: }

279: /*
280: If a stratum is active (non-zero dof), make it active in the coordinate DM.
281: */
282: static PetscErrorCode DMCreateCoordinateDM_Stag(DM dm,DM *dmc)
283: {
284:   PetscErrorCode  ierr;
285:   DM_Stag * const stag = (DM_Stag*)dm->data;
286:   PetscInt        dim;
287:   PetscBool       isstag,isproduct;


291:   if (!stag->coordinateDMType) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Before creating a coordinate DM, a type must be specified with DMStagSetCoordinateDMType()");

293:   DMGetDimension(dm,&dim);
294:   PetscStrcmp(stag->coordinateDMType,DMSTAG,&isstag);
295:   PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&isproduct);
296:   if (isstag) {
297:     DMStagCreateCompatibleDMStag(dm,
298:         stag->dof[0] > 0 ? dim : 0,
299:         stag->dof[1] > 0 ? dim : 0,
300:         stag->dof[2] > 0 ? dim : 0,
301:         stag->dof[3] > 0 ? dim : 0,
302:         dmc);
303:   } else if (isproduct) {
304:     DMCreate(PETSC_COMM_WORLD,dmc);
305:     DMSetType(*dmc,DMPRODUCT);
306:     DMSetDimension(*dmc,dim);
307:   } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported coordinate DM type %s",stag->coordinateDMType);
308:   return(0);
309: }

311: static PetscErrorCode DMGetNeighbors_Stag(DM dm,PetscInt *nRanks,const PetscMPIInt *ranks[])
312: {
313:   PetscErrorCode  ierr;
314:   DM_Stag * const stag = (DM_Stag*)dm->data;
315:   PetscInt        dim;

318:   DMGetDimension(dm,&dim);
319:   switch (dim) {
320:     case 1: *nRanks = 3; break;
321:     case 2: *nRanks = 9; break;
322:     case 3: *nRanks = 27; break;
323:     default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Get neighbors not implemented for dim = %D",dim);
324:   }
325:   *ranks = stag->neighbors;
326:   return(0);
327: }

329: static PetscErrorCode DMView_Stag(DM dm,PetscViewer viewer)
330: {
331:   PetscErrorCode  ierr;
332:   DM_Stag * const stag = (DM_Stag*)dm->data;
333:   PetscBool       isascii,viewAllRanks;
334:   PetscMPIInt     rank,size;
335:   PetscInt        dim,maxRanksToView,i;

338:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
339:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
340:   DMGetDimension(dm,&dim);
341:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
342:   if (isascii) {
343:     PetscViewerASCIIPrintf(viewer,"Dimension: %D\n",dim);
344:     switch (dim) {
345:       case 1:
346:         PetscViewerASCIIPrintf(viewer,"Global size: %D\n",stag->N[0]);
347:         break;
348:       case 2:
349:         PetscViewerASCIIPrintf(viewer,"Global sizes: %D x %D\n",stag->N[0],stag->N[1]);
350:         PetscViewerASCIIPrintf(viewer,"Parallel decomposition: %D x %D ranks\n",stag->nRanks[0],stag->nRanks[1]);
351:         break;
352:       case 3:
353:         PetscViewerASCIIPrintf(viewer,"Global sizes: %D x %D x %D\n",stag->N[0],stag->N[1],stag->N[2]);
354:         PetscViewerASCIIPrintf(viewer,"Parallel decomposition: %D x %D x %D ranks\n",stag->nRanks[0],stag->nRanks[1],stag->nRanks[2]);
355:         break;
356:       default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"not implemented for dim==%D",dim);
357:     }
358:     PetscViewerASCIIPrintf(viewer,"Boundary ghosting:");
359:     for (i=0; i<dim; ++i) {
360:       PetscViewerASCIIPrintf(viewer," %s",DMBoundaryTypes[stag->boundaryType[i]]);
361:     }
362:     PetscViewerASCIIPrintf(viewer,"\n");
363:     PetscViewerASCIIPrintf(viewer,"Elementwise ghost stencil: %s, width %D\n",DMStagStencilTypes[stag->stencilType],stag->stencilWidth);
364:     PetscViewerASCIIPrintf(viewer,"Stratum dof:");
365:     for (i=0; i<dim+1; ++i) {
366:       PetscViewerASCIIPrintf(viewer," %D:%D",i,stag->dof[i]);
367:     }
368:     PetscViewerASCIIPrintf(viewer,"\n");
369:     if(dm->coordinateDM) {
370:       PetscViewerASCIIPrintf(viewer,"Has coordinate DM\n");
371:     }
372:     maxRanksToView = 16;
373:     viewAllRanks = (PetscBool)(size <= maxRanksToView);
374:     if (viewAllRanks) {
375:       PetscViewerASCIIPushSynchronized(viewer);
376:       switch (dim) {
377:         case 1:
378:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local elements : %D (%D with ghosts)\n",rank,stag->n[0],stag->nGhost[0]);
379:           break;
380:         case 2:
381:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Rank coordinates (%d,%d)\n",rank,stag->rank[0],stag->rank[1]);
382:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local elements : %D x %D (%D x %D with ghosts)\n",rank,stag->n[0],stag->n[1],stag->nGhost[0],stag->nGhost[1]);
383:           break;
384:         case 3:
385:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Rank coordinates (%d,%d,%d)\n",rank,stag->rank[0],stag->rank[1],stag->rank[2]);
386:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local elements : %D x %D x %D (%D x %D x %D with ghosts)\n",rank,stag->n[0],stag->n[1],stag->n[2],stag->nGhost[0],stag->nGhost[1],stag->nGhost[2]);
387:           break;
388:         default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"not implemented for dim==%D",dim);
389:       }
390:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local native entries: %d\n",rank,stag->entries     );
391:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local entries total : %d\n",rank,stag->entriesGhost);
392:       PetscViewerFlush(viewer);
393:       PetscViewerASCIIPopSynchronized(viewer);
394:     } else {
395:       PetscViewerASCIIPrintf(viewer,"(Per-rank information omitted since >%D ranks used)\n",maxRanksToView);
396:     }
397:   }
398:   return(0);
399: }

401: static PetscErrorCode DMSetFromOptions_Stag(PetscOptionItems *PetscOptionsObject,DM dm)
402: {
403:   PetscErrorCode  ierr;
404:   DM_Stag * const stag = (DM_Stag*)dm->data;
405:   PetscInt        dim;

408:   DMGetDimension(dm,&dim);
409:   PetscOptionsHead(PetscOptionsObject,"DMStag Options");
410:   PetscOptionsInt("-stag_grid_x","Number of grid points in x direction","DMStagSetGlobalSizes",stag->N[0],&stag->N[0],NULL);
411:   if (dim > 1) { PetscOptionsInt("-stag_grid_y","Number of grid points in y direction","DMStagSetGlobalSizes",stag->N[1],&stag->N[1],NULL); }
412:   if (dim > 2) { PetscOptionsInt("-stag_grid_z","Number of grid points in z direction","DMStagSetGlobalSizes",stag->N[2],&stag->N[2],NULL); }
413:   PetscOptionsInt("-stag_ranks_x","Number of ranks in x direction","DMStagSetNumRanks",stag->nRanks[0],&stag->nRanks[0],NULL);
414:   if (dim > 1) { PetscOptionsInt("-stag_ranks_y","Number of ranks in y direction","DMStagSetNumRanks",stag->nRanks[1],&stag->nRanks[1],NULL); }
415:   if (dim > 2) { PetscOptionsInt("-stag_ranks_z","Number of ranks in z direction","DMStagSetNumRanks",stag->nRanks[2],&stag->nRanks[2],NULL); }
416:   PetscOptionsInt("-stag_stencil_width","Elementwise stencil width","DMStagSetStencilWidth",stag->stencilWidth,&stag->stencilWidth,NULL);
417:   PetscOptionsEnum("-stag_stencil_type","Elementwise stencil stype","DMStagSetStencilType",DMStagStencilTypes,(PetscEnum)stag->stencilType,(PetscEnum*)&stag->stencilType,NULL);
418:   PetscOptionsEnum("-stag_boundary_type_x","Treatment of (physical) boundaries in x direction","DMStagSetBoundaryTypes",DMBoundaryTypes,(PetscEnum)stag->boundaryType[0],(PetscEnum*)&stag->boundaryType[0],NULL);
419:   PetscOptionsEnum("-stag_boundary_type_y","Treatment of (physical) boundaries in y direction","DMStagSetBoundaryTypes",DMBoundaryTypes,(PetscEnum)stag->boundaryType[1],(PetscEnum*)&stag->boundaryType[1],NULL);
420:   PetscOptionsEnum("-stag_boundary_type_z","Treatment of (physical) boundaries in z direction","DMStagSetBoundaryTypes",DMBoundaryTypes,(PetscEnum)stag->boundaryType[2],(PetscEnum*)&stag->boundaryType[2],NULL);
421:   PetscOptionsInt("-stag_dof_0","Number of dof per 0-cell (vertex/corner)","DMStagSetDOF",stag->dof[0],&stag->dof[0],NULL);
422:   PetscOptionsInt("-stag_dof_1","Number of dof per 1-cell (edge)",         "DMStagSetDOF",stag->dof[1],&stag->dof[1],NULL);
423:   PetscOptionsInt("-stag_dof_2","Number of dof per 2-cell (face)",         "DMStagSetDOF",stag->dof[2],&stag->dof[2],NULL);
424:   PetscOptionsInt("-stag_dof_3","Number of dof per 3-cell (hexahedron)",   "DMStagSetDOF",stag->dof[3],&stag->dof[3],NULL);
425:   PetscOptionsTail();
426:   return(0);
427: }

429: /*MC
430:   DMSTAG = "stag" - A DM object representing a "staggered grid" or a structured cell complex.

432:   This implementation parallels the DMDA implementation in many ways, but allows degrees of freedom
433:   to be associated with all "strata" in a logically-rectangular grid: vertices, edges, faces, and elements.

435:   Level: beginner

437: .seealso: DM, DMPRODUCT, DMDA, DMPLEX, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMType, DMCreate(), DMSetType()
438: M*/

440: PETSC_EXTERN PetscErrorCode DMCreate_Stag(DM dm)
441: {
443:   DM_Stag        *stag;
444:   PetscInt       i,dim;

448:   PetscNewLog(dm,&stag);
449:   dm->data = stag;

451:   stag->gtol                                          = NULL;
452:   stag->ltog_injective                                = NULL;
453:   for (i=0; i<DMSTAG_MAX_STRATA; ++i) stag->dof[i]    = 0;
454:   for (i=0; i<DMSTAG_MAX_DIM;    ++i) stag->l[i]      = NULL;
455:   stag->stencilType                                   = DMSTAG_STENCIL_NONE;
456:   stag->stencilWidth                                  = 0;
457:   for (i=0; i<DMSTAG_MAX_DIM;    ++i) stag->nRanks[i] = -1;
458:   stag->coordinateDMType                              = NULL;

460:   DMGetDimension(dm,&dim);
461: #if defined(PETSC_USE_DEBUG)
462:   if (dim != 1 && dim != 2 && dim != 3) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DMSetDimension() must be called to set a dimension with value 1, 2, or 3");
463: #endif

465:   PetscMemzero(dm->ops,sizeof(*(dm->ops)));
466:   dm->ops->createcoordinatedm  = DMCreateCoordinateDM_Stag;
467:   dm->ops->createglobalvector  = DMCreateGlobalVector_Stag;
468:   dm->ops->createinterpolation = NULL;
469:   dm->ops->createlocalvector   = DMCreateLocalVector_Stag;
470:   dm->ops->creatematrix        = DMCreateMatrix_Stag;
471:   dm->ops->destroy             = DMDestroy_Stag;
472:   dm->ops->getneighbors        = DMGetNeighbors_Stag;
473:   dm->ops->globaltolocalbegin  = DMGlobalToLocalBegin_Stag;
474:   dm->ops->globaltolocalend    = DMGlobalToLocalEnd_Stag;
475:   dm->ops->localtoglobalbegin  = DMLocalToGlobalBegin_Stag;
476:   dm->ops->localtoglobalend    = DMLocalToGlobalEnd_Stag;
477:   dm->ops->setfromoptions      = DMSetFromOptions_Stag;
478:   switch (dim) {
479:     case 1: dm->ops->setup     = DMSetUp_Stag_1d; break;
480:     case 2: dm->ops->setup     = DMSetUp_Stag_2d; break;
481:     case 3: dm->ops->setup     = DMSetUp_Stag_3d; break;
482:     default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %d",dim);
483:   }
484:   dm->ops->clone               = DMClone_Stag;
485:   dm->ops->view                = DMView_Stag;
486:   dm->ops->getcompatibility    = DMGetCompatibility_Stag;
487:   return(0);
488: }