Actual source code: stag.c
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: if (PetscUnlikely(!dm->setupcalled)) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
51: VecCreateMPI(PetscObjectComm((PetscObject)dm),stag->entries,PETSC_DECIDE,vec);
52: VecSetDM(*vec,dm);
53: /* Could set some ops, as DMDA does */
54: VecSetLocalToGlobalMapping(*vec,dm->ltogmap);
55: return(0);
56: }
58: static PetscErrorCode DMCreateLocalVector_Stag(DM dm,Vec *vec)
59: {
60: PetscErrorCode ierr;
61: DM_Stag * const stag = (DM_Stag*)dm->data;
64: if (PetscUnlikely(!dm->setupcalled)) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
65: VecCreateSeq(PETSC_COMM_SELF,stag->entriesGhost,vec);
66: VecSetBlockSize(*vec,stag->entriesPerElement);
67: VecSetDM(*vec,dm);
68: return(0);
69: }
71: static PetscErrorCode DMCreateMatrix_Stag(DM dm,Mat *mat)
72: {
73: PetscErrorCode ierr;
74: MatType matType;
75: PetscBool isaij,isshell;
76: PetscInt entries,width,nNeighbors,dim,dof[DMSTAG_MAX_STRATA],stencilWidth;
77: DMStagStencilType stencilType;
78: ISLocalToGlobalMapping ltogmap;
81: if (PetscUnlikely(!dm->setupcalled)) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
82: DMGetDimension(dm,&dim);
83: DMGetMatType(dm,&matType);
84: PetscStrcmp(matType,MATAIJ,&isaij);
85: PetscStrcmp(matType,MATSHELL,&isshell);
86: DMStagGetEntries(dm,&entries);
87: DMStagGetDOF(dm,&dof[0],&dof[1],&dof[2],&dof[3]);
88: DMStagGetStencilWidth(dm,&stencilWidth);
89: DMStagGetStencilType(dm,&stencilType);
91: if (isaij) {
92: /* This implementation gives a very dense stencil, which is likely unsuitable for
93: real applications. */
94: switch (stencilType) {
95: case DMSTAG_STENCIL_NONE:
96: nNeighbors = 1;
97: break;
98: case DMSTAG_STENCIL_STAR:
99: switch (dim) {
100: case 1 :
101: nNeighbors = 2*stencilWidth + 1;
102: break;
103: case 2 :
104: nNeighbors = 4*stencilWidth + 3;
105: break;
106: case 3 :
107: nNeighbors = 6*stencilWidth + 5;
108: break;
109: default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %d",dim);
110: }
111: break;
112: case DMSTAG_STENCIL_BOX:
113: switch (dim) {
114: case 1 :
115: nNeighbors = (2*stencilWidth + 1);
116: break;
117: case 2 :
118: nNeighbors = (2*stencilWidth + 1) * (2*stencilWidth + 1);
119: break;
120: case 3 :
121: nNeighbors = (2*stencilWidth + 1) * (2*stencilWidth + 1) * (2*stencilWidth + 1);
122: break;
123: default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %d",dim);
124: }
125: break;
126: default : SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported stencil");
127: }
128: width = (dof[0] + dof[1] + dof[2] + dof[3]) * nNeighbors;
129: MatCreateAIJ(PetscObjectComm((PetscObject)dm),entries,entries,PETSC_DETERMINE,PETSC_DETERMINE,width,NULL,width,NULL,mat);
130: } else if (isshell) {
131: MatCreate(PetscObjectComm((PetscObject)dm),mat);
132: MatSetSizes(*mat,entries,entries,PETSC_DETERMINE,PETSC_DETERMINE);
133: MatSetType(*mat,MATSHELL);
134: MatSetUp(*mat);
135: } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for Mattype %s",matType);
137: DMGetLocalToGlobalMapping(dm,<ogmap);
138: MatSetLocalToGlobalMapping(*mat,ltogmap,ltogmap);
139: MatSetDM(*mat,dm);
140: return(0);
141: }
143: static PetscErrorCode DMGetCompatibility_Stag(DM dm,DM dm2,PetscBool *compatible,PetscBool *set)
144: {
145: PetscErrorCode ierr;
146: const DM_Stag * const stag = (DM_Stag*)dm->data;
147: const DM_Stag * const stag2 = (DM_Stag*)dm2->data;
148: PetscInt dim,dim2,i;
149: MPI_Comm comm;
150: PetscMPIInt sameComm;
151: DMType type2;
152: PetscBool sameType;
155: DMGetType(dm2,&type2);
156: PetscStrcmp(DMSTAG,type2,&sameType);
157: if (!sameType) {
158: PetscInfo1((PetscObject)dm,"DMStag compatibility check not implemented with DM of type %s\n",type2);
159: *set = PETSC_FALSE;
160: return(0);
161: }
163: PetscObjectGetComm((PetscObject)dm,&comm);
164: MPI_Comm_compare(comm,PetscObjectComm((PetscObject)dm2),&sameComm);
165: if (sameComm != MPI_IDENT) {
166: PetscInfo2((PetscObject)dm,"DMStag objects have different communicators: %d != %d\n",comm,PetscObjectComm((PetscObject)dm2));
167: *set = PETSC_FALSE;
168: return(0);
169: }
170: DMGetDimension(dm ,&dim);
171: DMGetDimension(dm2,&dim2);
172: if (dim != dim2) {
173: PetscInfo((PetscObject)dm,"DMStag objects have different dimensions");
174: *set = PETSC_TRUE;
175: *compatible = PETSC_FALSE;
176: return(0);
177: }
178: for (i=0; i<dim; ++i) {
179: if (stag->N[i] != stag2->N[i]) {
180: PetscInfo3((PetscObject)dm,"DMStag objects have different global numbers of elements in dimension %D: %D != %D\n",i,stag->n[i],stag2->n[i]);
181: *set = PETSC_TRUE;
182: *compatible = PETSC_FALSE;
183: return(0);
184: }
185: if (stag->n[i] != stag2->n[i]) {
186: PetscInfo3((PetscObject)dm,"DMStag objects have different local numbers of elements in dimension %D: %D != %D\n",i,stag->n[i],stag2->n[i]);
187: *set = PETSC_TRUE;
188: *compatible = PETSC_FALSE;
189: return(0);
190: }
191: if (stag->boundaryType[i] != stag2->boundaryType[i]) {
192: PetscInfo3((PetscObject)dm,"DMStag objects have different boundary types in dimension %d: %s != %s\n",i,stag->boundaryType[i],stag2->boundaryType[i]);
193: *set = PETSC_TRUE;
194: *compatible = PETSC_FALSE;
195: return(0);
196: }
197: }
198: /* Note: we include stencil type and width in the notion of compatibility, as this affects
199: the "atlas" (local subdomains). This might be irritating in legitimate cases
200: of wanting to transfer between two other-wise compatible DMs with different
201: stencil characteristics. */
202: if (stag->stencilType != stag2->stencilType) {
203: PetscInfo2((PetscObject)dm,"DMStag objects have different ghost stencil types: %s != %s\n",DMStagStencilTypes[stag->stencilType],DMStagStencilTypes[stag2->stencilType]);
204: *set = PETSC_TRUE;
205: *compatible = PETSC_FALSE;
206: return(0);
207: }
208: if (stag->stencilWidth != stag2->stencilWidth) {
209: PetscInfo2((PetscObject)dm,"DMStag objects have different ghost stencil widths: %D != %D\n",stag->stencilWidth,stag->stencilWidth);
210: *set = PETSC_TRUE;
211: *compatible = PETSC_FALSE;
212: return(0);
213: }
214: *set = PETSC_TRUE;
215: *compatible = PETSC_TRUE;
216: return(0);
217: }
219: /*
220: Note there are several orderings in play here.
221: 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).
222: Also in all cases, only subdomains which are the last in their dimension have partial elements.
224: 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.
225: 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.
226: 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.
227: */
229: static PetscErrorCode DMLocalToGlobalBegin_Stag(DM dm,Vec l,InsertMode mode,Vec g)
230: {
231: PetscErrorCode ierr;
232: DM_Stag * const stag = (DM_Stag*)dm->data;
235: if (mode == ADD_VALUES) {
236: VecScatterBegin(stag->gtol,l,g,mode,SCATTER_REVERSE);
237: } else if (mode == INSERT_VALUES) {
238: if (stag->ltog_injective) {
239: VecScatterBegin(stag->ltog_injective,l,g,mode,SCATTER_FORWARD);
240: } else {
241: VecScatterBegin(stag->gtol,l,g,mode,SCATTER_REVERSE_LOCAL);
242: }
243: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported InsertMode");
244: return(0);
245: }
247: static PetscErrorCode DMLocalToGlobalEnd_Stag(DM dm,Vec l,InsertMode mode,Vec g)
248: {
249: PetscErrorCode ierr;
250: DM_Stag * const stag = (DM_Stag*)dm->data;
253: if (mode == ADD_VALUES) {
254: VecScatterEnd(stag->gtol,l,g,mode,SCATTER_REVERSE);
255: } else if (mode == INSERT_VALUES) {
256: if (stag->ltog_injective) {
257: VecScatterEnd(stag->ltog_injective,l,g,mode,SCATTER_FORWARD);
258: } else {
259: VecScatterEnd(stag->gtol,l,g,mode,SCATTER_REVERSE_LOCAL);
260: }
261: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported InsertMode");
262: return(0);
263: }
265: static PetscErrorCode DMGlobalToLocalBegin_Stag(DM dm,Vec g,InsertMode mode,Vec l)
266: {
267: PetscErrorCode ierr;
268: DM_Stag * const stag = (DM_Stag*)dm->data;
271: VecScatterBegin(stag->gtol,g,l,mode,SCATTER_FORWARD);
272: return(0);
273: }
275: static PetscErrorCode DMGlobalToLocalEnd_Stag(DM dm,Vec g,InsertMode mode,Vec l)
276: {
277: PetscErrorCode ierr;
278: DM_Stag * const stag = (DM_Stag*)dm->data;
281: VecScatterEnd(stag->gtol,g,l,mode,SCATTER_FORWARD);
282: return(0);
283: }
285: /*
286: If a stratum is active (non-zero dof), make it active in the coordinate DM.
287: */
288: static PetscErrorCode DMCreateCoordinateDM_Stag(DM dm,DM *dmc)
289: {
290: PetscErrorCode ierr;
291: DM_Stag * const stag = (DM_Stag*)dm->data;
292: PetscInt dim;
293: PetscBool isstag,isproduct;
297: if (!stag->coordinateDMType) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Before creating a coordinate DM, a type must be specified with DMStagSetCoordinateDMType()");
299: DMGetDimension(dm,&dim);
300: PetscStrcmp(stag->coordinateDMType,DMSTAG,&isstag);
301: PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&isproduct);
302: if (isstag) {
303: DMStagCreateCompatibleDMStag(dm,
304: stag->dof[0] > 0 ? dim : 0,
305: stag->dof[1] > 0 ? dim : 0,
306: stag->dof[2] > 0 ? dim : 0,
307: stag->dof[3] > 0 ? dim : 0,
308: dmc);
309: } else if (isproduct) {
310: DMCreate(PETSC_COMM_WORLD,dmc);
311: DMSetType(*dmc,DMPRODUCT);
312: DMSetDimension(*dmc,dim);
313: } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported coordinate DM type %s",stag->coordinateDMType);
314: return(0);
315: }
317: static PetscErrorCode DMGetNeighbors_Stag(DM dm,PetscInt *nRanks,const PetscMPIInt *ranks[])
318: {
319: PetscErrorCode ierr;
320: DM_Stag * const stag = (DM_Stag*)dm->data;
321: PetscInt dim;
324: DMGetDimension(dm,&dim);
325: switch (dim) {
326: case 1: *nRanks = 3; break;
327: case 2: *nRanks = 9; break;
328: case 3: *nRanks = 27; break;
329: default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Get neighbors not implemented for dim = %D",dim);
330: }
331: *ranks = stag->neighbors;
332: return(0);
333: }
335: static PetscErrorCode DMView_Stag(DM dm,PetscViewer viewer)
336: {
337: PetscErrorCode ierr;
338: DM_Stag * const stag = (DM_Stag*)dm->data;
339: PetscBool isascii,viewAllRanks;
340: PetscMPIInt rank,size;
341: PetscInt dim,maxRanksToView,i;
344: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
345: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
346: DMGetDimension(dm,&dim);
347: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
348: if (isascii) {
349: PetscViewerASCIIPrintf(viewer,"Dimension: %D\n",dim);
350: switch (dim) {
351: case 1:
352: PetscViewerASCIIPrintf(viewer,"Global size: %D\n",stag->N[0]);
353: break;
354: case 2:
355: PetscViewerASCIIPrintf(viewer,"Global sizes: %D x %D\n",stag->N[0],stag->N[1]);
356: PetscViewerASCIIPrintf(viewer,"Parallel decomposition: %D x %D ranks\n",stag->nRanks[0],stag->nRanks[1]);
357: break;
358: case 3:
359: PetscViewerASCIIPrintf(viewer,"Global sizes: %D x %D x %D\n",stag->N[0],stag->N[1],stag->N[2]);
360: PetscViewerASCIIPrintf(viewer,"Parallel decomposition: %D x %D x %D ranks\n",stag->nRanks[0],stag->nRanks[1],stag->nRanks[2]);
361: break;
362: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"not implemented for dim==%D",dim);
363: }
364: PetscViewerASCIIPrintf(viewer,"Boundary ghosting:");
365: for (i=0; i<dim; ++i) {
366: PetscViewerASCIIPrintf(viewer," %s",DMBoundaryTypes[stag->boundaryType[i]]);
367: }
368: PetscViewerASCIIPrintf(viewer,"\n");
369: PetscViewerASCIIPrintf(viewer,"Elementwise ghost stencil: %s, width %D\n",DMStagStencilTypes[stag->stencilType],stag->stencilWidth);
370: PetscViewerASCIIPrintf(viewer,"Stratum dof:");
371: for (i=0; i<dim+1; ++i) {
372: PetscViewerASCIIPrintf(viewer," %D:%D",i,stag->dof[i]);
373: }
374: PetscViewerASCIIPrintf(viewer,"\n");
375: if (dm->coordinateDM) {
376: PetscViewerASCIIPrintf(viewer,"Has coordinate DM\n");
377: }
378: maxRanksToView = 16;
379: viewAllRanks = (PetscBool)(size <= maxRanksToView);
380: if (viewAllRanks) {
381: PetscViewerASCIIPushSynchronized(viewer);
382: switch (dim) {
383: case 1:
384: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local elements : %D (%D with ghosts)\n",rank,stag->n[0],stag->nGhost[0]);
385: break;
386: case 2:
387: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Rank coordinates (%d,%d)\n",rank,stag->rank[0],stag->rank[1]);
388: 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]);
389: break;
390: case 3:
391: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Rank coordinates (%d,%d,%d)\n",rank,stag->rank[0],stag->rank[1],stag->rank[2]);
392: 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]);
393: break;
394: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"not implemented for dim==%D",dim);
395: }
396: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local native entries: %d\n",rank,stag->entries);
397: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local entries total : %d\n",rank,stag->entriesGhost);
398: PetscViewerFlush(viewer);
399: PetscViewerASCIIPopSynchronized(viewer);
400: } else {
401: PetscViewerASCIIPrintf(viewer,"(Per-rank information omitted since >%D ranks used)\n",maxRanksToView);
402: }
403: }
404: return(0);
405: }
407: static PetscErrorCode DMSetFromOptions_Stag(PetscOptionItems *PetscOptionsObject,DM dm)
408: {
409: PetscErrorCode ierr;
410: DM_Stag * const stag = (DM_Stag*)dm->data;
411: PetscInt dim;
414: DMGetDimension(dm,&dim);
415: PetscOptionsHead(PetscOptionsObject,"DMStag Options");
416: PetscOptionsInt("-stag_grid_x","Number of grid points in x direction","DMStagSetGlobalSizes",stag->N[0],&stag->N[0],NULL);
417: if (dim > 1) { PetscOptionsInt("-stag_grid_y","Number of grid points in y direction","DMStagSetGlobalSizes",stag->N[1],&stag->N[1],NULL); }
418: if (dim > 2) { PetscOptionsInt("-stag_grid_z","Number of grid points in z direction","DMStagSetGlobalSizes",stag->N[2],&stag->N[2],NULL); }
419: PetscOptionsInt("-stag_ranks_x","Number of ranks in x direction","DMStagSetNumRanks",stag->nRanks[0],&stag->nRanks[0],NULL);
420: if (dim > 1) { PetscOptionsInt("-stag_ranks_y","Number of ranks in y direction","DMStagSetNumRanks",stag->nRanks[1],&stag->nRanks[1],NULL); }
421: if (dim > 2) { PetscOptionsInt("-stag_ranks_z","Number of ranks in z direction","DMStagSetNumRanks",stag->nRanks[2],&stag->nRanks[2],NULL); }
422: PetscOptionsInt("-stag_stencil_width","Elementwise stencil width","DMStagSetStencilWidth",stag->stencilWidth,&stag->stencilWidth,NULL);
423: PetscOptionsEnum("-stag_stencil_type","Elementwise stencil stype","DMStagSetStencilType",DMStagStencilTypes,(PetscEnum)stag->stencilType,(PetscEnum*)&stag->stencilType,NULL);
424: PetscOptionsEnum("-stag_boundary_type_x","Treatment of (physical) boundaries in x direction","DMStagSetBoundaryTypes",DMBoundaryTypes,(PetscEnum)stag->boundaryType[0],(PetscEnum*)&stag->boundaryType[0],NULL);
425: PetscOptionsEnum("-stag_boundary_type_y","Treatment of (physical) boundaries in y direction","DMStagSetBoundaryTypes",DMBoundaryTypes,(PetscEnum)stag->boundaryType[1],(PetscEnum*)&stag->boundaryType[1],NULL);
426: PetscOptionsEnum("-stag_boundary_type_z","Treatment of (physical) boundaries in z direction","DMStagSetBoundaryTypes",DMBoundaryTypes,(PetscEnum)stag->boundaryType[2],(PetscEnum*)&stag->boundaryType[2],NULL);
427: PetscOptionsInt("-stag_dof_0","Number of dof per 0-cell (vertex/corner)","DMStagSetDOF",stag->dof[0],&stag->dof[0],NULL);
428: PetscOptionsInt("-stag_dof_1","Number of dof per 1-cell (edge)", "DMStagSetDOF",stag->dof[1],&stag->dof[1],NULL);
429: PetscOptionsInt("-stag_dof_2","Number of dof per 2-cell (face)", "DMStagSetDOF",stag->dof[2],&stag->dof[2],NULL);
430: PetscOptionsInt("-stag_dof_3","Number of dof per 3-cell (hexahedron)", "DMStagSetDOF",stag->dof[3],&stag->dof[3],NULL);
431: PetscOptionsTail();
432: return(0);
433: }
435: /*MC
436: DMSTAG = "stag" - A DM object representing a "staggered grid" or a structured cell complex.
438: This implementation parallels the DMDA implementation in many ways, but allows degrees of freedom
439: to be associated with all "strata" in a logically-rectangular grid: vertices, edges, faces, and elements.
441: Level: beginner
443: .seealso: DM, DMPRODUCT, DMDA, DMPLEX, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMType, DMCreate(), DMSetType()
444: M*/
446: PETSC_EXTERN PetscErrorCode DMCreate_Stag(DM dm)
447: {
449: DM_Stag *stag;
450: PetscInt i,dim;
454: PetscNewLog(dm,&stag);
455: dm->data = stag;
457: stag->gtol = NULL;
458: stag->ltog_injective = NULL;
459: for (i=0; i<DMSTAG_MAX_STRATA; ++i) stag->dof[i] = 0;
460: for (i=0; i<DMSTAG_MAX_DIM; ++i) stag->l[i] = NULL;
461: stag->stencilType = DMSTAG_STENCIL_NONE;
462: stag->stencilWidth = 0;
463: for (i=0; i<DMSTAG_MAX_DIM; ++i) stag->nRanks[i] = -1;
464: stag->coordinateDMType = NULL;
466: DMGetDimension(dm,&dim);
467: 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");
469: PetscMemzero(dm->ops,sizeof(*(dm->ops)));
470: dm->ops->createcoordinatedm = DMCreateCoordinateDM_Stag;
471: dm->ops->createglobalvector = DMCreateGlobalVector_Stag;
472: dm->ops->createinterpolation = NULL;
473: dm->ops->createlocalvector = DMCreateLocalVector_Stag;
474: dm->ops->creatematrix = DMCreateMatrix_Stag;
475: dm->ops->destroy = DMDestroy_Stag;
476: dm->ops->getneighbors = DMGetNeighbors_Stag;
477: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Stag;
478: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Stag;
479: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Stag;
480: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Stag;
481: dm->ops->setfromoptions = DMSetFromOptions_Stag;
482: switch (dim) {
483: case 1: dm->ops->setup = DMSetUp_Stag_1d; break;
484: case 2: dm->ops->setup = DMSetUp_Stag_2d; break;
485: case 3: dm->ops->setup = DMSetUp_Stag_3d; break;
486: default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %d",dim);
487: }
488: dm->ops->clone = DMClone_Stag;
489: dm->ops->view = DMView_Stag;
490: dm->ops->getcompatibility = DMGetCompatibility_Stag;
491: return(0);
492: }