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 DMCreateFieldDecomposition_Stag(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
11: {
12: PetscInt f0, f1, f2, f3, dof0, dof1, dof2, dof3, n_entries, k, d, cnt, n_fields, dim;
13: DMStagStencil *stencil0, *stencil1, *stencil2, *stencil3;
15: DMGetDimension(dm, &dim);
16: DMStagGetDOF(dm, &dof0, &dof1, &dof2, &dof3);
17: DMStagGetEntriesPerElement(dm, &n_entries);
19: f0 = 1;
20: f1 = f2 = f3 = 0;
21: if (dim == 1) {
22: f1 = 1;
23: } else if (dim == 2) {
24: f1 = 2;
25: f2 = 1;
26: } else if (dim == 3) {
27: f1 = 3;
28: f2 = 3;
29: f3 = 1;
30: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
32: PetscCalloc1(f0 * dof0, &stencil0);
33: PetscCalloc1(f1 * dof1, &stencil1);
34: if (dim >= 2) PetscCalloc1(f2 * dof2, &stencil2);
35: if (dim >= 3) PetscCalloc1(f3 * dof3, &stencil3);
36: for (k = 0; k < f0; ++k) {
37: for (d = 0; d < dof0; ++d) {
38: stencil0[dof0 * k + d].i = 0;
39: stencil0[dof0 * k + d].j = 0;
40: stencil0[dof0 * k + d].j = 0;
41: }
42: }
43: for (k = 0; k < f1; ++k) {
44: for (d = 0; d < dof1; ++d) {
45: stencil1[dof1 * k + d].i = 0;
46: stencil1[dof1 * k + d].j = 0;
47: stencil1[dof1 * k + d].j = 0;
48: }
49: }
50: if (dim >= 2) {
51: for (k = 0; k < f2; ++k) {
52: for (d = 0; d < dof2; ++d) {
53: stencil2[dof2 * k + d].i = 0;
54: stencil2[dof2 * k + d].j = 0;
55: stencil2[dof2 * k + d].j = 0;
56: }
57: }
58: }
59: if (dim >= 3) {
60: for (k = 0; k < f3; ++k) {
61: for (d = 0; d < dof3; ++d) {
62: stencil3[dof3 * k + d].i = 0;
63: stencil3[dof3 * k + d].j = 0;
64: stencil3[dof3 * k + d].j = 0;
65: }
66: }
67: }
69: n_fields = 0;
70: if (dof0 != 0) ++n_fields;
71: if (dof1 != 0) ++n_fields;
72: if (dim >= 2 && dof2 != 0) ++n_fields;
73: if (dim >= 3 && dof3 != 0) ++n_fields;
74: if (len) *len = n_fields;
76: if (islist) {
77: PetscMalloc1(n_fields, islist);
79: if (dim == 1) {
80: /* face, element */
81: for (d = 0; d < dof0; ++d) {
82: stencil0[d].loc = DMSTAG_LEFT;
83: stencil0[d].c = d;
84: }
85: for (d = 0; d < dof1; ++d) {
86: stencil1[d].loc = DMSTAG_ELEMENT;
87: stencil1[d].c = d;
88: }
89: } else if (dim == 2) {
90: /* vertex, edge(down,left), element */
91: for (d = 0; d < dof0; ++d) {
92: stencil0[d].loc = DMSTAG_DOWN_LEFT;
93: stencil0[d].c = d;
94: }
95: /* edge */
96: cnt = 0;
97: for (d = 0; d < dof1; ++d) {
98: stencil1[cnt].loc = DMSTAG_DOWN;
99: stencil1[cnt].c = d;
100: ++cnt;
101: }
102: for (d = 0; d < dof1; ++d) {
103: stencil1[cnt].loc = DMSTAG_LEFT;
104: stencil1[cnt].c = d;
105: ++cnt;
106: }
107: /* element */
108: for (d = 0; d < dof2; ++d) {
109: stencil2[d].loc = DMSTAG_ELEMENT;
110: stencil2[d].c = d;
111: }
112: } else if (dim == 3) {
113: /* vertex, edge(down,left), face(down,left,back), element */
114: for (d = 0; d < dof0; ++d) {
115: stencil0[d].loc = DMSTAG_BACK_DOWN_LEFT;
116: stencil0[d].c = d;
117: }
118: /* edges */
119: cnt = 0;
120: for (d = 0; d < dof1; ++d) {
121: stencil1[cnt].loc = DMSTAG_BACK_DOWN;
122: stencil1[cnt].c = d;
123: ++cnt;
124: }
125: for (d = 0; d < dof1; ++d) {
126: stencil1[cnt].loc = DMSTAG_BACK_LEFT;
127: stencil1[cnt].c = d;
128: ++cnt;
129: }
130: for (d = 0; d < dof1; ++d) {
131: stencil1[cnt].loc = DMSTAG_DOWN_LEFT;
132: stencil1[cnt].c = d;
133: ++cnt;
134: }
135: /* faces */
136: cnt = 0;
137: for (d = 0; d < dof2; ++d) {
138: stencil2[cnt].loc = DMSTAG_BACK;
139: stencil2[cnt].c = d;
140: ++cnt;
141: }
142: for (d = 0; d < dof2; ++d) {
143: stencil2[cnt].loc = DMSTAG_DOWN;
144: stencil2[cnt].c = d;
145: ++cnt;
146: }
147: for (d = 0; d < dof2; ++d) {
148: stencil2[cnt].loc = DMSTAG_LEFT;
149: stencil2[cnt].c = d;
150: ++cnt;
151: }
152: /* elements */
153: for (d = 0; d < dof3; ++d) {
154: stencil3[d].loc = DMSTAG_ELEMENT;
155: stencil3[d].c = d;
156: }
157: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
159: cnt = 0;
160: if (dof0 != 0) {
161: DMStagCreateISFromStencils(dm, f0 * dof0, stencil0, &(*islist)[cnt]);
162: ++cnt;
163: }
164: if (dof1 != 0) {
165: DMStagCreateISFromStencils(dm, f1 * dof1, stencil1, &(*islist)[cnt]);
166: ++cnt;
167: }
168: if (dim >= 2 && dof2 != 0) {
169: DMStagCreateISFromStencils(dm, f2 * dof2, stencil2, &(*islist)[cnt]);
170: ++cnt;
171: }
172: if (dim >= 3 && dof3 != 0) {
173: DMStagCreateISFromStencils(dm, f3 * dof3, stencil3, &(*islist)[cnt]);
174: ++cnt;
175: }
176: }
178: if (namelist) {
179: PetscMalloc1(n_fields, namelist);
180: cnt = 0;
181: if (dim == 1) {
182: if (dof0 != 0) {
183: PetscStrallocpy("vertex", &(*namelist)[cnt]);
184: ++cnt;
185: }
186: if (dof1 != 0) {
187: PetscStrallocpy("element", &(*namelist)[cnt]);
188: ++cnt;
189: }
190: } else if (dim == 2) {
191: if (dof0 != 0) {
192: PetscStrallocpy("vertex", &(*namelist)[cnt]);
193: ++cnt;
194: }
195: if (dof1 != 0) {
196: PetscStrallocpy("face", &(*namelist)[cnt]);
197: ++cnt;
198: }
199: if (dof2 != 0) {
200: PetscStrallocpy("element", &(*namelist)[cnt]);
201: ++cnt;
202: }
203: } else if (dim == 3) {
204: if (dof0 != 0) {
205: PetscStrallocpy("vertex", &(*namelist)[cnt]);
206: ++cnt;
207: }
208: if (dof1 != 0) {
209: PetscStrallocpy("edge", &(*namelist)[cnt]);
210: ++cnt;
211: }
212: if (dof2 != 0) {
213: PetscStrallocpy("face", &(*namelist)[cnt]);
214: ++cnt;
215: }
216: if (dof3 != 0) {
217: PetscStrallocpy("element", &(*namelist)[cnt]);
218: ++cnt;
219: }
220: }
221: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
222: if (dmlist) {
223: PetscMalloc1(n_fields, dmlist);
224: cnt = 0;
225: if (dof0 != 0) {
226: DMStagCreateCompatibleDMStag(dm, dof0, 0, 0, 0, &(*dmlist)[cnt]);
227: ++cnt;
228: }
229: if (dof1 != 0) {
230: DMStagCreateCompatibleDMStag(dm, 0, dof1, 0, 0, &(*dmlist)[cnt]);
231: ++cnt;
232: }
233: if (dim >= 2 && dof2 != 0) {
234: DMStagCreateCompatibleDMStag(dm, 0, 0, dof2, 0, &(*dmlist)[cnt]);
235: ++cnt;
236: }
237: if (dim >= 3 && dof3 != 0) {
238: DMStagCreateCompatibleDMStag(dm, 0, 0, 0, dof3, &(*dmlist)[cnt]);
239: ++cnt;
240: }
241: }
242: PetscFree(stencil0);
243: PetscFree(stencil1);
244: if (dim >= 2) PetscFree(stencil2);
245: if (dim >= 3) PetscFree(stencil3);
246: return 0;
247: }
249: static PetscErrorCode DMClone_Stag(DM dm, DM *newdm)
250: {
251: /* Destroy the DM created by generic logic in DMClone() */
252: if (*newdm) DMDestroy(newdm);
253: DMStagDuplicateWithoutSetup(dm, PetscObjectComm((PetscObject)dm), newdm);
254: DMSetUp(*newdm);
255: return 0;
256: }
258: static PetscErrorCode DMCoarsen_Stag(DM dm, MPI_Comm comm, DM *dmc)
259: {
260: const DM_Stag *const stag = (DM_Stag *)dm->data;
261: PetscInt d, dim;
263: DMStagDuplicateWithoutSetup(dm, comm, dmc);
264: DMSetOptionsPrefix(*dmc, ((PetscObject)dm)->prefix);
265: DMGetDimension(dm, &dim);
267: DMStagSetGlobalSizes(*dmc, stag->N[0] / 2, stag->N[1] / 2, stag->N[2] / 2);
268: {
269: PetscInt *l[DMSTAG_MAX_DIM];
270: for (d = 0; d < dim; ++d) {
271: PetscInt i;
272: PetscMalloc1(stag->nRanks[d], &l[d]);
273: for (i = 0; i < stag->nRanks[d]; ++i) {
275: l[d][i] = stag->l[d][i] / 2; /* Just halve everything */
276: }
277: }
278: DMStagSetOwnershipRanges(*dmc, l[0], l[1], l[2]);
279: for (d = 0; d < dim; ++d) PetscFree(l[d]);
280: }
281: DMSetUp(*dmc);
283: if (dm->coordinates[0].dm) { /* Note that with product coordinates, dm->coordinates = NULL, so we check the DM */
284: DM coordinate_dm, coordinate_dmc;
285: PetscBool isstag, isprod;
287: DMGetCoordinateDM(dm, &coordinate_dm);
288: PetscObjectTypeCompare((PetscObject)coordinate_dm, DMSTAG, &isstag);
289: PetscObjectTypeCompare((PetscObject)coordinate_dm, DMPRODUCT, &isprod);
290: if (isstag) {
291: DMStagSetUniformCoordinatesExplicit(*dmc, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0); /* Coordinates will be overwritten */
292: DMGetCoordinateDM(*dmc, &coordinate_dmc);
293: DMStagRestrictSimple(coordinate_dm, dm->coordinates[0].x, coordinate_dmc, (*dmc)->coordinates[0].x);
294: } else if (isprod) {
295: DMStagSetUniformCoordinatesProduct(*dmc, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0); /* Coordinates will be overwritten */
296: DMGetCoordinateDM(*dmc, &coordinate_dmc);
297: for (d = 0; d < dim; ++d) {
298: DM subdm_coarse, subdm_coord_coarse, subdm_fine, subdm_coord_fine;
300: DMProductGetDM(coordinate_dm, d, &subdm_fine);
301: DMGetCoordinateDM(subdm_fine, &subdm_coord_fine);
302: DMProductGetDM(coordinate_dmc, d, &subdm_coarse);
303: DMGetCoordinateDM(subdm_coarse, &subdm_coord_coarse);
304: DMStagRestrictSimple(subdm_coord_fine, subdm_fine->coordinates[0].xl, subdm_coord_coarse, subdm_coarse->coordinates[0].xl);
305: }
306: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown coordinate DM type");
307: }
308: return 0;
309: }
311: static PetscErrorCode DMRefine_Stag(DM dm, MPI_Comm comm, DM *dmc)
312: {
313: const DM_Stag *const stag = (DM_Stag *)dm->data;
315: DMStagDuplicateWithoutSetup(dm, comm, dmc);
316: DMSetOptionsPrefix(*dmc, ((PetscObject)dm)->prefix);
317: DMStagSetGlobalSizes(*dmc, stag->N[0] * 2, stag->N[1] * 2, stag->N[2] * 2);
318: {
319: PetscInt dim, d;
320: PetscInt *l[DMSTAG_MAX_DIM];
321: DMGetDimension(dm, &dim);
322: for (d = 0; d < dim; ++d) {
323: PetscInt i;
324: PetscMalloc1(stag->nRanks[d], &l[d]);
325: for (i = 0; i < stag->nRanks[d]; ++i) { l[d][i] = stag->l[d][i] * 2; /* Just double everything */ }
326: }
327: DMStagSetOwnershipRanges(*dmc, l[0], l[1], l[2]);
328: for (d = 0; d < dim; ++d) PetscFree(l[d]);
329: }
330: DMSetUp(*dmc);
331: /* Note: For now, we do not refine coordinates */
332: return 0;
333: }
335: static PetscErrorCode DMDestroy_Stag(DM dm)
336: {
337: DM_Stag *stag;
338: PetscInt i;
340: stag = (DM_Stag *)dm->data;
341: for (i = 0; i < DMSTAG_MAX_DIM; ++i) PetscFree(stag->l[i]);
342: VecScatterDestroy(&stag->gtol);
343: VecScatterDestroy(&stag->ltog_injective);
344: PetscFree(stag->neighbors);
345: PetscFree(stag->locationOffsets);
346: PetscFree(stag->coordinateDMType);
347: PetscFree(stag);
348: return 0;
349: }
351: static PetscErrorCode DMCreateGlobalVector_Stag(DM dm, Vec *vec)
352: {
353: DM_Stag *const stag = (DM_Stag *)dm->data;
356: VecCreate(PetscObjectComm((PetscObject)dm), vec);
357: VecSetSizes(*vec, stag->entries, PETSC_DETERMINE);
358: VecSetType(*vec, dm->vectype);
359: VecSetDM(*vec, dm);
360: /* Could set some ops, as DMDA does */
361: VecSetLocalToGlobalMapping(*vec, dm->ltogmap);
362: return 0;
363: }
365: static PetscErrorCode DMCreateLocalVector_Stag(DM dm, Vec *vec)
366: {
367: DM_Stag *const stag = (DM_Stag *)dm->data;
370: VecCreate(PETSC_COMM_SELF, vec);
371: VecSetSizes(*vec, stag->entriesGhost, PETSC_DETERMINE);
372: VecSetType(*vec, dm->vectype);
373: VecSetBlockSize(*vec, stag->entriesPerElement);
374: VecSetDM(*vec, dm);
375: return 0;
376: }
378: /* Helper function to check for the limited situations for which interpolation
379: and restriction functions are implemented */
380: static PetscErrorCode CheckTransferOperatorRequirements_Private(DM dmc, DM dmf)
381: {
382: PetscInt dim, stencilWidthc, stencilWidthf, nf[DMSTAG_MAX_DIM], nc[DMSTAG_MAX_DIM], doff[DMSTAG_MAX_STRATA], dofc[DMSTAG_MAX_STRATA];
384: DMGetDimension(dmc, &dim);
385: DMStagGetStencilWidth(dmc, &stencilWidthc);
387: DMStagGetStencilWidth(dmf, &stencilWidthf);
389: DMStagGetLocalSizes(dmf, &nf[0], &nf[1], &nf[2]);
390: DMStagGetLocalSizes(dmc, &nc[0], &nc[1], &nc[2]);
391: for (PetscInt d = 0; d < dim; ++d)
393: DMStagGetDOF(dmc, &dofc[0], &dofc[1], &dofc[2], &dofc[3]);
394: DMStagGetDOF(dmf, &doff[0], &doff[1], &doff[2], &doff[3]);
395: for (PetscInt d = 0; d < dim + 1; ++d)
397: return 0;
398: }
400: /* Since the interpolation uses MATMAIJ for dof > 0 we convert requests for non-MATAIJ baseded matrices to MATAIJ.
401: This is a bit of a hack; the reason for it is partially because -dm_mat_type defines the
402: matrix type for both the operator matrices and the interpolation matrices so that users
403: can select matrix types of base MATAIJ for accelerators
405: Note: The ConvertToAIJ() code below *has been copied from dainterp.c*! ConvertToAIJ() should perhaps be placed somewhere
406: in mat/utils to avoid code duplication, but then the DMStag and DMDA code would need to include the private Mat headers.
407: Since it is only used in two places, I have simply duplicated the code to avoid the need to exposure the private
408: Mat routines in parts of DM. If we find a need for ConvertToAIJ() elsewhere, then we should consolidate it to one
409: place in mat/utils.
410: */
411: static PetscErrorCode ConvertToAIJ(MatType intype, MatType *outtype)
412: {
413: PetscInt i;
414: char const *types[3] = {MATAIJ, MATSEQAIJ, MATMPIAIJ};
415: PetscBool flg;
417: *outtype = MATAIJ;
418: for (i = 0; i < 3; i++) {
419: PetscStrbeginswith(intype, types[i], &flg);
420: if (flg) {
421: *outtype = intype;
422: break;
423: }
424: }
425: return 0;
426: }
428: static PetscErrorCode DMCreateInterpolation_Stag(DM dmc, DM dmf, Mat *A, Vec *vec)
429: {
430: PetscInt dim, entriesf, entriesc, doff[DMSTAG_MAX_STRATA];
431: ISLocalToGlobalMapping ltogmf, ltogmc;
432: MatType mattype;
434: CheckTransferOperatorRequirements_Private(dmc, dmf);
436: DMStagGetEntries(dmf, &entriesf);
437: DMStagGetEntries(dmc, &entriesc);
438: DMGetLocalToGlobalMapping(dmf, <ogmf);
439: DMGetLocalToGlobalMapping(dmc, <ogmc);
441: MatCreate(PetscObjectComm((PetscObject)dmc), A);
442: MatSetSizes(*A, entriesf, entriesc, PETSC_DECIDE, PETSC_DECIDE);
443: ConvertToAIJ(dmc->mattype, &mattype);
444: MatSetType(*A, mattype);
445: MatSetLocalToGlobalMapping(*A, ltogmf, ltogmc);
447: DMGetDimension(dmc, &dim);
448: DMStagGetDOF(dmf, &doff[0], &doff[1], &doff[2], &doff[3]);
449: if (dim == 1) {
450: DMStagPopulateInterpolation1d_a_b_Private(dmc, dmf, *A);
451: } else if (dim == 2) {
452: if (doff[0] == 0) {
453: DMStagPopulateInterpolation2d_0_a_b_Private(dmc, dmf, *A);
454: } else
455: SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "No default interpolation available between 2d DMStag objects with %" PetscInt_FMT " dof/vertex, %" PetscInt_FMT " dof/face and %" PetscInt_FMT " dof/element", doff[0], doff[1], doff[2]);
456: } else if (dim == 3) {
457: if (doff[0] == 0 && doff[1] == 0) {
458: DMStagPopulateInterpolation3d_0_0_a_b_Private(dmc, dmf, *A);
459: } else
460: SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "No default interpolation available between 3d DMStag objects with %" PetscInt_FMT " dof/vertex, %" PetscInt_FMT " dof/edge, %" PetscInt_FMT " dof/face and %" PetscInt_FMT " dof/element", doff[0], doff[1], doff[2], doff[3]);
461: } else SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
462: MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY);
463: MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY);
465: if (vec) *vec = NULL;
466: return 0;
467: }
469: static PetscErrorCode DMCreateRestriction_Stag(DM dmc, DM dmf, Mat *A)
470: {
471: PetscInt dim, entriesf, entriesc, doff[DMSTAG_MAX_STRATA];
472: ISLocalToGlobalMapping ltogmf, ltogmc;
473: MatType mattype;
475: CheckTransferOperatorRequirements_Private(dmc, dmf);
477: DMStagGetEntries(dmf, &entriesf);
478: DMStagGetEntries(dmc, &entriesc);
479: DMGetLocalToGlobalMapping(dmf, <ogmf);
480: DMGetLocalToGlobalMapping(dmc, <ogmc);
482: MatCreate(PetscObjectComm((PetscObject)dmc), A);
483: MatSetSizes(*A, entriesc, entriesf, PETSC_DECIDE, PETSC_DECIDE); /* Note transpose wrt interpolation */
484: ConvertToAIJ(dmc->mattype, &mattype);
485: MatSetType(*A, mattype);
486: MatSetLocalToGlobalMapping(*A, ltogmc, ltogmf); /* Note transpose wrt interpolation */
488: DMGetDimension(dmc, &dim);
489: DMStagGetDOF(dmf, &doff[0], &doff[1], &doff[2], &doff[3]);
490: if (dim == 1) {
491: DMStagPopulateRestriction1d_a_b_Private(dmc, dmf, *A);
492: } else if (dim == 2) {
493: if (doff[0] == 0) {
494: DMStagPopulateRestriction2d_0_a_b_Private(dmc, dmf, *A);
495: } else
496: SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "No default restriction available between 2d DMStag objects with %" PetscInt_FMT " dof/vertex, %" PetscInt_FMT " dof/face and %" PetscInt_FMT " dof/element", doff[0], doff[1], doff[2]);
497: } else if (dim == 3) {
498: if (doff[0] == 0 && doff[0] == 0) {
499: DMStagPopulateRestriction3d_0_0_a_b_Private(dmc, dmf, *A);
500: } else
501: SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "No default restriction available between 3d DMStag objects with %" PetscInt_FMT " dof/vertex, %" PetscInt_FMT " dof/edge, %" PetscInt_FMT " dof/face and %" PetscInt_FMT " dof/element", doff[0], doff[1], doff[2], doff[3]);
502: } else SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
504: MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY);
505: MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY);
506: return 0;
507: }
509: static PetscErrorCode DMCreateMatrix_Stag(DM dm, Mat *mat)
510: {
511: MatType mat_type;
512: PetscBool is_shell, is_aij;
513: PetscInt dim, entries;
514: ISLocalToGlobalMapping ltogmap;
517: DMGetDimension(dm, &dim);
518: DMGetMatType(dm, &mat_type);
519: DMStagGetEntries(dm, &entries);
520: MatCreate(PetscObjectComm((PetscObject)dm), mat);
521: MatSetSizes(*mat, entries, entries, PETSC_DETERMINE, PETSC_DETERMINE);
522: MatSetType(*mat, mat_type);
523: MatSetUp(*mat);
524: DMGetLocalToGlobalMapping(dm, <ogmap);
525: MatSetLocalToGlobalMapping(*mat, ltogmap, ltogmap);
526: MatSetDM(*mat, dm);
528: /* Compare to similar and perhaps superior logic in DMCreateMatrix_DA, which creates
529: the matrix first and then performs this logic by checking for preallocation functions */
530: PetscObjectBaseTypeCompare((PetscObject)*mat, MATAIJ, &is_aij);
531: if (!is_aij) { PetscObjectBaseTypeCompare((PetscObject)*mat, MATSEQAIJ, &is_aij); }
532: if (!is_aij) { PetscObjectBaseTypeCompare((PetscObject)*mat, MATMPIAIJ, &is_aij); }
533: PetscStrcmp(mat_type, MATSHELL, &is_shell);
534: if (is_aij) {
535: Mat preallocator;
536: PetscInt m, n;
537: const PetscBool fill_with_zeros = PETSC_FALSE;
539: MatCreate(PetscObjectComm((PetscObject)dm), &preallocator);
540: MatSetType(preallocator, MATPREALLOCATOR);
541: MatGetLocalSize(*mat, &m, &n);
542: MatSetSizes(preallocator, m, n, PETSC_DECIDE, PETSC_DECIDE);
543: MatSetLocalToGlobalMapping(preallocator, ltogmap, ltogmap);
544: MatSetUp(preallocator);
545: switch (dim) {
546: case 1:
547: DMCreateMatrix_Stag_1D_AIJ_Assemble(dm, preallocator);
548: break;
549: case 2:
550: DMCreateMatrix_Stag_2D_AIJ_Assemble(dm, preallocator);
551: break;
552: case 3:
553: DMCreateMatrix_Stag_3D_AIJ_Assemble(dm, preallocator);
554: break;
555: default:
556: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
557: }
558: MatPreallocatorPreallocate(preallocator, fill_with_zeros, *mat);
559: MatDestroy(&preallocator);
561: if (!dm->prealloc_only) {
562: /* Bind to CPU before assembly, to prevent unnecessary copies of zero entries from CPU to GPU */
563: MatBindToCPU(*mat, PETSC_TRUE);
564: switch (dim) {
565: case 1:
566: DMCreateMatrix_Stag_1D_AIJ_Assemble(dm, *mat);
567: break;
568: case 2:
569: DMCreateMatrix_Stag_2D_AIJ_Assemble(dm, *mat);
570: break;
571: case 3:
572: DMCreateMatrix_Stag_3D_AIJ_Assemble(dm, *mat);
573: break;
574: default:
575: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
576: }
577: MatBindToCPU(*mat, PETSC_FALSE);
578: }
579: } else if (is_shell) {
580: /* nothing more to do */
581: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for Mattype %s", mat_type);
582: return 0;
583: }
585: static PetscErrorCode DMGetCompatibility_Stag(DM dm, DM dm2, PetscBool *compatible, PetscBool *set)
586: {
587: const DM_Stag *const stag = (DM_Stag *)dm->data;
588: const DM_Stag *const stag2 = (DM_Stag *)dm2->data;
589: PetscInt dim, dim2, i;
590: MPI_Comm comm;
591: PetscMPIInt sameComm;
592: DMType type2;
593: PetscBool sameType;
595: DMGetType(dm2, &type2);
596: PetscStrcmp(DMSTAG, type2, &sameType);
597: if (!sameType) {
598: PetscInfo((PetscObject)dm, "DMStag compatibility check not implemented with DM of type %s\n", type2);
599: *set = PETSC_FALSE;
600: return 0;
601: }
603: PetscObjectGetComm((PetscObject)dm, &comm);
604: MPI_Comm_compare(comm, PetscObjectComm((PetscObject)dm2), &sameComm);
605: if (sameComm != MPI_IDENT) {
606: PetscInfo((PetscObject)dm, "DMStag objects have different communicators: %" PETSC_MPI_COMM_FMT " != %" PETSC_MPI_COMM_FMT "\n", comm, PetscObjectComm((PetscObject)dm2));
607: *set = PETSC_FALSE;
608: return 0;
609: }
610: DMGetDimension(dm, &dim);
611: DMGetDimension(dm2, &dim2);
612: if (dim != dim2) {
613: PetscInfo((PetscObject)dm, "DMStag objects have different dimensions");
614: *set = PETSC_TRUE;
615: *compatible = PETSC_FALSE;
616: return 0;
617: }
618: for (i = 0; i < dim; ++i) {
619: if (stag->N[i] != stag2->N[i]) {
620: PetscInfo((PetscObject)dm, "DMStag objects have different global numbers of elements in dimension %" PetscInt_FMT ": %" PetscInt_FMT " != %" PetscInt_FMT "\n", i, stag->n[i], stag2->n[i]);
621: *set = PETSC_TRUE;
622: *compatible = PETSC_FALSE;
623: return 0;
624: }
625: if (stag->n[i] != stag2->n[i]) {
626: PetscInfo((PetscObject)dm, "DMStag objects have different local numbers of elements in dimension %" PetscInt_FMT ": %" PetscInt_FMT " != %" PetscInt_FMT "\n", i, stag->n[i], stag2->n[i]);
627: *set = PETSC_TRUE;
628: *compatible = PETSC_FALSE;
629: return 0;
630: }
631: if (stag->boundaryType[i] != stag2->boundaryType[i]) {
632: PetscInfo((PetscObject)dm, "DMStag objects have different boundary types in dimension %" PetscInt_FMT ": %s != %s\n", i, DMBoundaryTypes[stag->boundaryType[i]], DMBoundaryTypes[stag2->boundaryType[i]]);
633: *set = PETSC_TRUE;
634: *compatible = PETSC_FALSE;
635: return 0;
636: }
637: }
638: /* Note: we include stencil type and width in the notion of compatibility, as this affects
639: the "atlas" (local subdomains). This might be irritating in legitimate cases
640: of wanting to transfer between two other-wise compatible DMs with different
641: stencil characteristics. */
642: if (stag->stencilType != stag2->stencilType) {
643: PetscInfo((PetscObject)dm, "DMStag objects have different ghost stencil types: %s != %s\n", DMStagStencilTypes[stag->stencilType], DMStagStencilTypes[stag2->stencilType]);
644: *set = PETSC_TRUE;
645: *compatible = PETSC_FALSE;
646: return 0;
647: }
648: if (stag->stencilWidth != stag2->stencilWidth) {
649: PetscInfo((PetscObject)dm, "DMStag objects have different ghost stencil widths: %" PetscInt_FMT " != %" PetscInt_FMT "\n", stag->stencilWidth, stag->stencilWidth);
650: *set = PETSC_TRUE;
651: *compatible = PETSC_FALSE;
652: return 0;
653: }
654: *set = PETSC_TRUE;
655: *compatible = PETSC_TRUE;
656: return 0;
657: }
659: static PetscErrorCode DMHasCreateInjection_Stag(DM dm, PetscBool *flg)
660: {
663: *flg = PETSC_FALSE;
664: return 0;
665: }
667: /*
668: Note there are several orderings in play here.
669: 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).
670: Also in all cases, only subdomains which are the last in their dimension have partial elements.
672: 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.
673: 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.
674: 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.
675: */
677: static PetscErrorCode DMLocalToGlobalBegin_Stag(DM dm, Vec l, InsertMode mode, Vec g)
678: {
679: DM_Stag *const stag = (DM_Stag *)dm->data;
681: if (mode == ADD_VALUES) {
682: VecScatterBegin(stag->gtol, l, g, mode, SCATTER_REVERSE);
683: } else if (mode == INSERT_VALUES) {
684: if (stag->ltog_injective) {
685: VecScatterBegin(stag->ltog_injective, l, g, mode, SCATTER_FORWARD);
686: } else {
687: VecScatterBegin(stag->gtol, l, g, mode, SCATTER_REVERSE_LOCAL);
688: }
689: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported InsertMode");
690: return 0;
691: }
693: static PetscErrorCode DMLocalToGlobalEnd_Stag(DM dm, Vec l, InsertMode mode, Vec g)
694: {
695: DM_Stag *const stag = (DM_Stag *)dm->data;
697: if (mode == ADD_VALUES) {
698: VecScatterEnd(stag->gtol, l, g, mode, SCATTER_REVERSE);
699: } else if (mode == INSERT_VALUES) {
700: if (stag->ltog_injective) {
701: VecScatterEnd(stag->ltog_injective, l, g, mode, SCATTER_FORWARD);
702: } else {
703: VecScatterEnd(stag->gtol, l, g, mode, SCATTER_REVERSE_LOCAL);
704: }
705: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported InsertMode");
706: return 0;
707: }
709: static PetscErrorCode DMGlobalToLocalBegin_Stag(DM dm, Vec g, InsertMode mode, Vec l)
710: {
711: DM_Stag *const stag = (DM_Stag *)dm->data;
713: VecScatterBegin(stag->gtol, g, l, mode, SCATTER_FORWARD);
714: return 0;
715: }
717: static PetscErrorCode DMGlobalToLocalEnd_Stag(DM dm, Vec g, InsertMode mode, Vec l)
718: {
719: DM_Stag *const stag = (DM_Stag *)dm->data;
721: VecScatterEnd(stag->gtol, g, l, mode, SCATTER_FORWARD);
722: return 0;
723: }
725: /*
726: If a stratum is active (non-zero dof), make it active in the coordinate DM.
727: */
728: static PetscErrorCode DMCreateCoordinateDM_Stag(DM dm, DM *dmc)
729: {
730: DM_Stag *const stag = (DM_Stag *)dm->data;
731: PetscInt dim;
732: PetscBool isstag, isproduct;
737: DMGetDimension(dm, &dim);
738: PetscStrcmp(stag->coordinateDMType, DMSTAG, &isstag);
739: PetscStrcmp(stag->coordinateDMType, DMPRODUCT, &isproduct);
740: if (isstag) {
741: DMStagCreateCompatibleDMStag(dm, stag->dof[0] > 0 ? dim : 0, stag->dof[1] > 0 ? dim : 0, stag->dof[2] > 0 ? dim : 0, stag->dof[3] > 0 ? dim : 0, dmc);
742: } else if (isproduct) {
743: DMCreate(PETSC_COMM_WORLD, dmc);
744: DMSetType(*dmc, DMPRODUCT);
745: DMSetDimension(*dmc, dim);
746: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate DM type %s", stag->coordinateDMType);
747: return 0;
748: }
750: static PetscErrorCode DMGetNeighbors_Stag(DM dm, PetscInt *nRanks, const PetscMPIInt *ranks[])
751: {
752: DM_Stag *const stag = (DM_Stag *)dm->data;
753: PetscInt dim;
755: DMGetDimension(dm, &dim);
756: switch (dim) {
757: case 1:
758: *nRanks = 3;
759: break;
760: case 2:
761: *nRanks = 9;
762: break;
763: case 3:
764: *nRanks = 27;
765: break;
766: default:
767: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Get neighbors not implemented for dim = %" PetscInt_FMT, dim);
768: }
769: *ranks = stag->neighbors;
770: return 0;
771: }
773: static PetscErrorCode DMView_Stag(DM dm, PetscViewer viewer)
774: {
775: DM_Stag *const stag = (DM_Stag *)dm->data;
776: PetscBool isascii, viewAllRanks;
777: PetscMPIInt rank, size;
778: PetscInt dim, maxRanksToView, i;
780: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
781: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
782: DMGetDimension(dm, &dim);
783: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
784: if (isascii) {
785: PetscViewerASCIIPrintf(viewer, "Dimension: %" PetscInt_FMT "\n", dim);
786: switch (dim) {
787: case 1:
788: PetscViewerASCIIPrintf(viewer, "Global size: %" PetscInt_FMT "\n", stag->N[0]);
789: break;
790: case 2:
791: PetscViewerASCIIPrintf(viewer, "Global sizes: %" PetscInt_FMT " x %" PetscInt_FMT "\n", stag->N[0], stag->N[1]);
792: PetscViewerASCIIPrintf(viewer, "Parallel decomposition: %" PetscInt_FMT " x %" PetscInt_FMT " ranks\n", stag->nRanks[0], stag->nRanks[1]);
793: break;
794: case 3:
795: PetscViewerASCIIPrintf(viewer, "Global sizes: %" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT "\n", stag->N[0], stag->N[1], stag->N[2]);
796: PetscViewerASCIIPrintf(viewer, "Parallel decomposition: %" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT " ranks\n", stag->nRanks[0], stag->nRanks[1], stag->nRanks[2]);
797: break;
798: default:
799: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "not implemented for dim==%" PetscInt_FMT, dim);
800: }
801: PetscViewerASCIIPrintf(viewer, "Boundary ghosting:");
802: for (i = 0; i < dim; ++i) PetscViewerASCIIPrintf(viewer, " %s", DMBoundaryTypes[stag->boundaryType[i]]);
803: PetscViewerASCIIPrintf(viewer, "\n");
804: PetscViewerASCIIPrintf(viewer, "Elementwise ghost stencil: %s", DMStagStencilTypes[stag->stencilType]);
805: if (stag->stencilType != DMSTAG_STENCIL_NONE) {
806: PetscViewerASCIIPrintf(viewer, ", width %" PetscInt_FMT "\n", stag->stencilWidth);
807: } else {
808: PetscViewerASCIIPrintf(viewer, "\n");
809: }
810: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per vertex (0D)\n", stag->dof[0]);
811: if (dim == 3) PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per edge (1D)\n", stag->dof[1]);
812: if (dim > 1) PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per face (%" PetscInt_FMT "D)\n", stag->dof[dim - 1], dim - 1);
813: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per element (%" PetscInt_FMT "D)\n", stag->dof[dim], dim);
814: if (dm->coordinates[0].dm) PetscViewerASCIIPrintf(viewer, "Has coordinate DM\n");
815: maxRanksToView = 16;
816: viewAllRanks = (PetscBool)(size <= maxRanksToView);
817: if (viewAllRanks) {
818: PetscViewerASCIIPushSynchronized(viewer);
819: switch (dim) {
820: case 1:
821: PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local elements : %" PetscInt_FMT " (%" PetscInt_FMT " with ghosts)\n", rank, stag->n[0], stag->nGhost[0]);
822: break;
823: case 2:
824: PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Rank coordinates (%d,%d)\n", rank, stag->rank[0], stag->rank[1]);
825: PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local elements : %" PetscInt_FMT " x %" PetscInt_FMT " (%" PetscInt_FMT " x %" PetscInt_FMT " with ghosts)\n", rank, stag->n[0], stag->n[1], stag->nGhost[0], stag->nGhost[1]);
826: break;
827: case 3:
828: PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Rank coordinates (%d,%d,%d)\n", rank, stag->rank[0], stag->rank[1], stag->rank[2]);
829: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local elements : %" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT " (%" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT " with ghosts)\n", rank, stag->n[0], stag->n[1],
830: stag->n[2], stag->nGhost[0], stag->nGhost[1], stag->nGhost[2]));
831: break;
832: default:
833: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "not implemented for dim==%" PetscInt_FMT, dim);
834: }
835: PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local native entries: %" PetscInt_FMT "\n", rank, stag->entries);
836: PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local entries total : %" PetscInt_FMT "\n", rank, stag->entriesGhost);
837: PetscViewerFlush(viewer);
838: PetscViewerASCIIPopSynchronized(viewer);
839: } else {
840: PetscViewerASCIIPrintf(viewer, "(Per-rank information omitted since >%" PetscInt_FMT " ranks used)\n", maxRanksToView);
841: }
842: }
843: return 0;
844: }
846: static PetscErrorCode DMSetFromOptions_Stag(DM dm, PetscOptionItems *PetscOptionsObject)
847: {
848: DM_Stag *const stag = (DM_Stag *)dm->data;
849: PetscInt dim;
851: DMGetDimension(dm, &dim);
852: PetscOptionsHeadBegin(PetscOptionsObject, "DMStag Options");
853: PetscOptionsInt("-stag_grid_x", "Number of grid points in x direction", "DMStagSetGlobalSizes", stag->N[0], &stag->N[0], NULL);
854: if (dim > 1) PetscOptionsInt("-stag_grid_y", "Number of grid points in y direction", "DMStagSetGlobalSizes", stag->N[1], &stag->N[1], NULL);
855: if (dim > 2) PetscOptionsInt("-stag_grid_z", "Number of grid points in z direction", "DMStagSetGlobalSizes", stag->N[2], &stag->N[2], NULL);
856: PetscOptionsInt("-stag_ranks_x", "Number of ranks in x direction", "DMStagSetNumRanks", stag->nRanks[0], &stag->nRanks[0], NULL);
857: if (dim > 1) PetscOptionsInt("-stag_ranks_y", "Number of ranks in y direction", "DMStagSetNumRanks", stag->nRanks[1], &stag->nRanks[1], NULL);
858: if (dim > 2) PetscOptionsInt("-stag_ranks_z", "Number of ranks in z direction", "DMStagSetNumRanks", stag->nRanks[2], &stag->nRanks[2], NULL);
859: PetscOptionsInt("-stag_stencil_width", "Elementwise stencil width", "DMStagSetStencilWidth", stag->stencilWidth, &stag->stencilWidth, NULL);
860: PetscOptionsEnum("-stag_stencil_type", "Elementwise stencil stype", "DMStagSetStencilType", DMStagStencilTypes, (PetscEnum)stag->stencilType, (PetscEnum *)&stag->stencilType, NULL);
861: PetscOptionsEnum("-stag_boundary_type_x", "Treatment of (physical) boundaries in x direction", "DMStagSetBoundaryTypes", DMBoundaryTypes, (PetscEnum)stag->boundaryType[0], (PetscEnum *)&stag->boundaryType[0], NULL);
862: PetscOptionsEnum("-stag_boundary_type_y", "Treatment of (physical) boundaries in y direction", "DMStagSetBoundaryTypes", DMBoundaryTypes, (PetscEnum)stag->boundaryType[1], (PetscEnum *)&stag->boundaryType[1], NULL);
863: PetscOptionsEnum("-stag_boundary_type_z", "Treatment of (physical) boundaries in z direction", "DMStagSetBoundaryTypes", DMBoundaryTypes, (PetscEnum)stag->boundaryType[2], (PetscEnum *)&stag->boundaryType[2], NULL);
864: PetscOptionsInt("-stag_dof_0", "Number of dof per 0-cell (vertex)", "DMStagSetDOF", stag->dof[0], &stag->dof[0], NULL);
865: PetscOptionsInt("-stag_dof_1", "Number of dof per 1-cell (element in 1D, face in 2D, edge in 3D)", "DMStagSetDOF", stag->dof[1], &stag->dof[1], NULL);
866: PetscOptionsInt("-stag_dof_2", "Number of dof per 2-cell (element in 2D, face in 3D)", "DMStagSetDOF", stag->dof[2], &stag->dof[2], NULL);
867: PetscOptionsInt("-stag_dof_3", "Number of dof per 3-cell (element in 3D)", "DMStagSetDOF", stag->dof[3], &stag->dof[3], NULL);
868: PetscOptionsHeadEnd();
869: return 0;
870: }
872: /*MC
873: DMSTAG - `"stag"` - A DM object representing a "staggered grid" or a structured cell complex.
875: See the [manual chapter on DMStag](/docs/manual/dmstag).
877: This implementation parallels the `DMDA` implementation in many ways, but allows degrees of freedom
878: to be associated with all "strata" in a logically-rectangular grid.
880: Each stratum can be characterized by the dimension of the entities ("points", to borrow the `DMPLEX`
881: terminology), from 0- to 3-dimensional.
883: In some cases this numbering is used directly, for example with `DMStagGetDOF()`.
884: To allow easier reading and to some extent more similar code between different-dimensional implementations
885: of the same problem, we associate canonical names for each type of point, for each dimension of DMStag.
887: * 1-dimensional DMStag objects have vertices (0D) and elements (1D).
888: * 2-dimensional DMStag objects have vertices (0D), faces (1D), and elements (2D).
889: * 3-dimensional DMStag objects have vertices (0D), edges (1D), faces (2D), and elements (3D).
891: This naming is reflected when viewing a DMStag object with `DMView()`, and in forming
892: convenient options prefixes when creating a decomposition with `DMCreateFieldDecomposition()`.
894: Level: beginner
896: .seealso: `DM`, `DMPRODUCT`, `DMDA`, `DMPLEX`, `DMStagCreate1d()`, `DMStagCreate2d()`, `DMStagCreate3d()`, `DMType`, `DMCreate()`, `DMSetType()`
897: M*/
899: PETSC_EXTERN PetscErrorCode DMCreate_Stag(DM dm)
900: {
901: DM_Stag *stag;
902: PetscInt i, dim;
905: PetscNew(&stag);
906: dm->data = stag;
908: stag->gtol = NULL;
909: stag->ltog_injective = NULL;
910: for (i = 0; i < DMSTAG_MAX_STRATA; ++i) stag->dof[i] = 0;
911: for (i = 0; i < DMSTAG_MAX_DIM; ++i) stag->l[i] = NULL;
912: stag->stencilType = DMSTAG_STENCIL_NONE;
913: stag->stencilWidth = 0;
914: for (i = 0; i < DMSTAG_MAX_DIM; ++i) stag->nRanks[i] = -1;
915: stag->coordinateDMType = NULL;
917: DMGetDimension(dm, &dim);
920: PetscMemzero(dm->ops, sizeof(*(dm->ops)));
921: dm->ops->createcoordinatedm = DMCreateCoordinateDM_Stag;
922: dm->ops->createglobalvector = DMCreateGlobalVector_Stag;
923: dm->ops->createlocalvector = DMCreateLocalVector_Stag;
924: dm->ops->creatematrix = DMCreateMatrix_Stag;
925: dm->ops->hascreateinjection = DMHasCreateInjection_Stag;
926: dm->ops->refine = DMRefine_Stag;
927: dm->ops->coarsen = DMCoarsen_Stag;
928: dm->ops->createinterpolation = DMCreateInterpolation_Stag;
929: dm->ops->createrestriction = DMCreateRestriction_Stag;
930: dm->ops->destroy = DMDestroy_Stag;
931: dm->ops->getneighbors = DMGetNeighbors_Stag;
932: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Stag;
933: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Stag;
934: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Stag;
935: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Stag;
936: dm->ops->setfromoptions = DMSetFromOptions_Stag;
937: switch (dim) {
938: case 1:
939: dm->ops->setup = DMSetUp_Stag_1d;
940: break;
941: case 2:
942: dm->ops->setup = DMSetUp_Stag_2d;
943: break;
944: case 3:
945: dm->ops->setup = DMSetUp_Stag_3d;
946: break;
947: default:
948: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
949: }
950: dm->ops->clone = DMClone_Stag;
951: dm->ops->view = DMView_Stag;
952: dm->ops->getcompatibility = DMGetCompatibility_Stag;
953: dm->ops->createfielddecomposition = DMCreateFieldDecomposition_Stag;
954: return 0;
955: }