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, &ltogmf);
439:   DMGetLocalToGlobalMapping(dmc, &ltogmc);

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, &ltogmf);
480:   DMGetLocalToGlobalMapping(dmc, &ltogmc);

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, &ltogmap);
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: }