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:   PetscFunctionBegin;
 16:   PetscCall(DMGetDimension(dm, &dim));
 17:   PetscCall(DMStagGetDOF(dm, &dof0, &dof1, &dof2, &dof3));
 18:   PetscCall(DMStagGetEntriesPerElement(dm, &n_entries));

 20:   f0 = 1;
 21:   f1 = f2 = f3 = 0;
 22:   if (dim == 1) {
 23:     f1 = 1;
 24:   } else if (dim == 2) {
 25:     f1 = 2;
 26:     f2 = 1;
 27:   } else if (dim == 3) {
 28:     f1 = 3;
 29:     f2 = 3;
 30:     f3 = 1;
 31:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);

 33:   PetscCall(PetscCalloc1(f0 * dof0, &stencil0));
 34:   PetscCall(PetscCalloc1(f1 * dof1, &stencil1));
 35:   if (dim >= 2) PetscCall(PetscCalloc1(f2 * dof2, &stencil2));
 36:   if (dim >= 3) PetscCall(PetscCalloc1(f3 * dof3, &stencil3));
 37:   for (k = 0; k < f0; ++k) {
 38:     for (d = 0; d < dof0; ++d) {
 39:       stencil0[dof0 * k + d].i = 0;
 40:       stencil0[dof0 * k + d].j = 0;
 41:       stencil0[dof0 * k + d].j = 0;
 42:     }
 43:   }
 44:   for (k = 0; k < f1; ++k) {
 45:     for (d = 0; d < dof1; ++d) {
 46:       stencil1[dof1 * k + d].i = 0;
 47:       stencil1[dof1 * k + d].j = 0;
 48:       stencil1[dof1 * k + d].j = 0;
 49:     }
 50:   }
 51:   if (dim >= 2) {
 52:     for (k = 0; k < f2; ++k) {
 53:       for (d = 0; d < dof2; ++d) {
 54:         stencil2[dof2 * k + d].i = 0;
 55:         stencil2[dof2 * k + d].j = 0;
 56:         stencil2[dof2 * k + d].j = 0;
 57:       }
 58:     }
 59:   }
 60:   if (dim >= 3) {
 61:     for (k = 0; k < f3; ++k) {
 62:       for (d = 0; d < dof3; ++d) {
 63:         stencil3[dof3 * k + d].i = 0;
 64:         stencil3[dof3 * k + d].j = 0;
 65:         stencil3[dof3 * k + d].j = 0;
 66:       }
 67:     }
 68:   }

 70:   n_fields = 0;
 71:   if (dof0 != 0) ++n_fields;
 72:   if (dof1 != 0) ++n_fields;
 73:   if (dim >= 2 && dof2 != 0) ++n_fields;
 74:   if (dim >= 3 && dof3 != 0) ++n_fields;
 75:   if (len) *len = n_fields;

 77:   if (islist) {
 78:     PetscCall(PetscMalloc1(n_fields, islist));

 80:     if (dim == 1) {
 81:       /* face, element */
 82:       for (d = 0; d < dof0; ++d) {
 83:         stencil0[d].loc = DMSTAG_LEFT;
 84:         stencil0[d].c   = d;
 85:       }
 86:       for (d = 0; d < dof1; ++d) {
 87:         stencil1[d].loc = DMSTAG_ELEMENT;
 88:         stencil1[d].c   = d;
 89:       }
 90:     } else if (dim == 2) {
 91:       /* vertex, edge(down,left), element */
 92:       for (d = 0; d < dof0; ++d) {
 93:         stencil0[d].loc = DMSTAG_DOWN_LEFT;
 94:         stencil0[d].c   = d;
 95:       }
 96:       /* edge */
 97:       cnt = 0;
 98:       for (d = 0; d < dof1; ++d) {
 99:         stencil1[cnt].loc = DMSTAG_DOWN;
100:         stencil1[cnt].c   = d;
101:         ++cnt;
102:       }
103:       for (d = 0; d < dof1; ++d) {
104:         stencil1[cnt].loc = DMSTAG_LEFT;
105:         stencil1[cnt].c   = d;
106:         ++cnt;
107:       }
108:       /* element */
109:       for (d = 0; d < dof2; ++d) {
110:         stencil2[d].loc = DMSTAG_ELEMENT;
111:         stencil2[d].c   = d;
112:       }
113:     } else if (dim == 3) {
114:       /* vertex, edge(down,left), face(down,left,back), element */
115:       for (d = 0; d < dof0; ++d) {
116:         stencil0[d].loc = DMSTAG_BACK_DOWN_LEFT;
117:         stencil0[d].c   = d;
118:       }
119:       /* edges */
120:       cnt = 0;
121:       for (d = 0; d < dof1; ++d) {
122:         stencil1[cnt].loc = DMSTAG_BACK_DOWN;
123:         stencil1[cnt].c   = d;
124:         ++cnt;
125:       }
126:       for (d = 0; d < dof1; ++d) {
127:         stencil1[cnt].loc = DMSTAG_BACK_LEFT;
128:         stencil1[cnt].c   = d;
129:         ++cnt;
130:       }
131:       for (d = 0; d < dof1; ++d) {
132:         stencil1[cnt].loc = DMSTAG_DOWN_LEFT;
133:         stencil1[cnt].c   = d;
134:         ++cnt;
135:       }
136:       /* faces */
137:       cnt = 0;
138:       for (d = 0; d < dof2; ++d) {
139:         stencil2[cnt].loc = DMSTAG_BACK;
140:         stencil2[cnt].c   = d;
141:         ++cnt;
142:       }
143:       for (d = 0; d < dof2; ++d) {
144:         stencil2[cnt].loc = DMSTAG_DOWN;
145:         stencil2[cnt].c   = d;
146:         ++cnt;
147:       }
148:       for (d = 0; d < dof2; ++d) {
149:         stencil2[cnt].loc = DMSTAG_LEFT;
150:         stencil2[cnt].c   = d;
151:         ++cnt;
152:       }
153:       /* elements */
154:       for (d = 0; d < dof3; ++d) {
155:         stencil3[d].loc = DMSTAG_ELEMENT;
156:         stencil3[d].c   = d;
157:       }
158:     } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);

160:     cnt = 0;
161:     if (dof0 != 0) {
162:       PetscCall(DMStagCreateISFromStencils(dm, f0 * dof0, stencil0, &(*islist)[cnt]));
163:       ++cnt;
164:     }
165:     if (dof1 != 0) {
166:       PetscCall(DMStagCreateISFromStencils(dm, f1 * dof1, stencil1, &(*islist)[cnt]));
167:       ++cnt;
168:     }
169:     if (dim >= 2 && dof2 != 0) {
170:       PetscCall(DMStagCreateISFromStencils(dm, f2 * dof2, stencil2, &(*islist)[cnt]));
171:       ++cnt;
172:     }
173:     if (dim >= 3 && dof3 != 0) {
174:       PetscCall(DMStagCreateISFromStencils(dm, f3 * dof3, stencil3, &(*islist)[cnt]));
175:       ++cnt;
176:     }
177:   }

179:   if (namelist) {
180:     PetscCall(PetscMalloc1(n_fields, namelist));
181:     cnt = 0;
182:     if (dim == 1) {
183:       if (dof0 != 0) {
184:         PetscCall(PetscStrallocpy("vertex", &(*namelist)[cnt]));
185:         ++cnt;
186:       }
187:       if (dof1 != 0) {
188:         PetscCall(PetscStrallocpy("element", &(*namelist)[cnt]));
189:         ++cnt;
190:       }
191:     } else if (dim == 2) {
192:       if (dof0 != 0) {
193:         PetscCall(PetscStrallocpy("vertex", &(*namelist)[cnt]));
194:         ++cnt;
195:       }
196:       if (dof1 != 0) {
197:         PetscCall(PetscStrallocpy("face", &(*namelist)[cnt]));
198:         ++cnt;
199:       }
200:       if (dof2 != 0) {
201:         PetscCall(PetscStrallocpy("element", &(*namelist)[cnt]));
202:         ++cnt;
203:       }
204:     } else if (dim == 3) {
205:       if (dof0 != 0) {
206:         PetscCall(PetscStrallocpy("vertex", &(*namelist)[cnt]));
207:         ++cnt;
208:       }
209:       if (dof1 != 0) {
210:         PetscCall(PetscStrallocpy("edge", &(*namelist)[cnt]));
211:         ++cnt;
212:       }
213:       if (dof2 != 0) {
214:         PetscCall(PetscStrallocpy("face", &(*namelist)[cnt]));
215:         ++cnt;
216:       }
217:       if (dof3 != 0) {
218:         PetscCall(PetscStrallocpy("element", &(*namelist)[cnt]));
219:         ++cnt;
220:       }
221:     }
222:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
223:   if (dmlist) {
224:     PetscCall(PetscMalloc1(n_fields, dmlist));
225:     cnt = 0;
226:     if (dof0 != 0) {
227:       PetscCall(DMStagCreateCompatibleDMStag(dm, dof0, 0, 0, 0, &(*dmlist)[cnt]));
228:       ++cnt;
229:     }
230:     if (dof1 != 0) {
231:       PetscCall(DMStagCreateCompatibleDMStag(dm, 0, dof1, 0, 0, &(*dmlist)[cnt]));
232:       ++cnt;
233:     }
234:     if (dim >= 2 && dof2 != 0) {
235:       PetscCall(DMStagCreateCompatibleDMStag(dm, 0, 0, dof2, 0, &(*dmlist)[cnt]));
236:       ++cnt;
237:     }
238:     if (dim >= 3 && dof3 != 0) {
239:       PetscCall(DMStagCreateCompatibleDMStag(dm, 0, 0, 0, dof3, &(*dmlist)[cnt]));
240:       ++cnt;
241:     }
242:   }
243:   PetscCall(PetscFree(stencil0));
244:   PetscCall(PetscFree(stencil1));
245:   if (dim >= 2) PetscCall(PetscFree(stencil2));
246:   if (dim >= 3) PetscCall(PetscFree(stencil3));
247:   PetscFunctionReturn(PETSC_SUCCESS);
248: }

250: static PetscErrorCode DMClone_Stag(DM dm, DM *newdm)
251: {
252:   PetscFunctionBegin;
253:   /* Destroy the DM created by generic logic in DMClone() */
254:   if (*newdm) PetscCall(DMDestroy(newdm));
255:   PetscCall(DMStagDuplicateWithoutSetup(dm, PetscObjectComm((PetscObject)dm), newdm));
256:   PetscCall(DMSetUp(*newdm));
257:   PetscFunctionReturn(PETSC_SUCCESS);
258: }

260: static PetscErrorCode DMCoarsen_Stag(DM dm, MPI_Comm comm, DM *dmc)
261: {
262:   const DM_Stag *const stag = (DM_Stag *)dm->data;
263:   PetscInt             d, dim;

265:   PetscFunctionBegin;
266:   PetscCall(DMStagDuplicateWithoutSetup(dm, comm, dmc));
267:   PetscCall(DMSetOptionsPrefix(*dmc, ((PetscObject)dm)->prefix));
268:   PetscCall(DMGetDimension(dm, &dim));
269:   for (d = 0; d < dim; ++d) PetscCheck(stag->N[d] % 2 == 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "coarsening not supported except for even numbers of elements in each dimension ");
270:   PetscCall(DMStagSetGlobalSizes(*dmc, stag->N[0] / 2, stag->N[1] / 2, stag->N[2] / 2));
271:   {
272:     PetscInt *l[DMSTAG_MAX_DIM];
273:     for (d = 0; d < dim; ++d) {
274:       PetscInt i;
275:       PetscCall(PetscMalloc1(stag->nRanks[d], &l[d]));
276:       for (i = 0; i < stag->nRanks[d]; ++i) {
277:         PetscCheck(stag->l[d][i] % 2 == 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "coarsening not supported except for an even number of elements in each direction on each rank");
278:         l[d][i] = stag->l[d][i] / 2; /* Just halve everything */
279:       }
280:     }
281:     PetscCall(DMStagSetOwnershipRanges(*dmc, l[0], l[1], l[2]));
282:     for (d = 0; d < dim; ++d) PetscCall(PetscFree(l[d]));
283:   }
284:   PetscCall(DMSetUp(*dmc));

286:   if (dm->coordinates[0].dm) { /* Note that with product coordinates, dm->coordinates = NULL, so we check the DM */
287:     DM        coordinate_dm, coordinate_dmc;
288:     PetscBool isstag, isprod;

290:     PetscCall(DMGetCoordinateDM(dm, &coordinate_dm));
291:     PetscCall(PetscObjectTypeCompare((PetscObject)coordinate_dm, DMSTAG, &isstag));
292:     PetscCall(PetscObjectTypeCompare((PetscObject)coordinate_dm, DMPRODUCT, &isprod));
293:     if (isstag) {
294:       PetscCall(DMStagSetUniformCoordinatesExplicit(*dmc, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)); /* Coordinates will be overwritten */
295:       PetscCall(DMGetCoordinateDM(*dmc, &coordinate_dmc));
296:       PetscCall(DMStagRestrictSimple(coordinate_dm, dm->coordinates[0].x, coordinate_dmc, (*dmc)->coordinates[0].x));
297:     } else if (isprod) {
298:       PetscCall(DMStagSetUniformCoordinatesProduct(*dmc, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)); /* Coordinates will be overwritten */
299:       PetscCall(DMGetCoordinateDM(*dmc, &coordinate_dmc));
300:       for (d = 0; d < dim; ++d) {
301:         DM subdm_coarse, subdm_coord_coarse, subdm_fine, subdm_coord_fine;

303:         PetscCall(DMProductGetDM(coordinate_dm, d, &subdm_fine));
304:         PetscCall(DMGetCoordinateDM(subdm_fine, &subdm_coord_fine));
305:         PetscCall(DMProductGetDM(coordinate_dmc, d, &subdm_coarse));
306:         PetscCall(DMGetCoordinateDM(subdm_coarse, &subdm_coord_coarse));
307:         PetscCall(DMStagRestrictSimple(subdm_coord_fine, subdm_fine->coordinates[0].xl, subdm_coord_coarse, subdm_coarse->coordinates[0].xl));
308:       }
309:     } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown coordinate DM type");
310:   }
311:   PetscFunctionReturn(PETSC_SUCCESS);
312: }

314: static PetscErrorCode DMRefine_Stag(DM dm, MPI_Comm comm, DM *dmc)
315: {
316:   const DM_Stag *const stag = (DM_Stag *)dm->data;

318:   PetscFunctionBegin;
319:   PetscCall(DMStagDuplicateWithoutSetup(dm, comm, dmc));
320:   PetscCall(DMSetOptionsPrefix(*dmc, ((PetscObject)dm)->prefix));
321:   PetscCall(DMStagSetGlobalSizes(*dmc, stag->N[0] * 2, stag->N[1] * 2, stag->N[2] * 2));
322:   {
323:     PetscInt  dim, d;
324:     PetscInt *l[DMSTAG_MAX_DIM];
325:     PetscCall(DMGetDimension(dm, &dim));
326:     for (d = 0; d < dim; ++d) {
327:       PetscInt i;
328:       PetscCall(PetscMalloc1(stag->nRanks[d], &l[d]));
329:       for (i = 0; i < stag->nRanks[d]; ++i) { l[d][i] = stag->l[d][i] * 2; /* Just double everything */ }
330:     }
331:     PetscCall(DMStagSetOwnershipRanges(*dmc, l[0], l[1], l[2]));
332:     for (d = 0; d < dim; ++d) PetscCall(PetscFree(l[d]));
333:   }
334:   PetscCall(DMSetUp(*dmc));
335:   /* Note: For now, we do not refine coordinates */
336:   PetscFunctionReturn(PETSC_SUCCESS);
337: }

339: static PetscErrorCode DMDestroy_Stag(DM dm)
340: {
341:   DM_Stag *stag;
342:   PetscInt i;

344:   PetscFunctionBegin;
345:   stag = (DM_Stag *)dm->data;
346:   for (i = 0; i < DMSTAG_MAX_DIM; ++i) PetscCall(PetscFree(stag->l[i]));
347:   PetscCall(VecScatterDestroy(&stag->gtol));
348:   PetscCall(VecScatterDestroy(&stag->ltog_injective));
349:   PetscCall(PetscFree(stag->neighbors));
350:   PetscCall(PetscFree(stag->locationOffsets));
351:   PetscCall(PetscFree(stag->coordinateDMType));
352:   PetscCall(PetscFree(stag));
353:   PetscFunctionReturn(PETSC_SUCCESS);
354: }

356: static PetscErrorCode DMCreateGlobalVector_Stag(DM dm, Vec *vec)
357: {
358:   DM_Stag *const stag = (DM_Stag *)dm->data;

360:   PetscFunctionBegin;
361:   PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called after DMSetUp()");
362:   PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), vec));
363:   PetscCall(VecSetSizes(*vec, stag->entries, PETSC_DETERMINE));
364:   PetscCall(VecSetType(*vec, dm->vectype));
365:   PetscCall(VecSetDM(*vec, dm));
366:   /* Could set some ops, as DMDA does */
367:   PetscCall(VecSetLocalToGlobalMapping(*vec, dm->ltogmap));
368:   PetscFunctionReturn(PETSC_SUCCESS);
369: }

371: static PetscErrorCode DMCreateLocalVector_Stag(DM dm, Vec *vec)
372: {
373:   DM_Stag *const stag = (DM_Stag *)dm->data;

375:   PetscFunctionBegin;
376:   PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called after DMSetUp()");
377:   PetscCall(VecCreate(PETSC_COMM_SELF, vec));
378:   PetscCall(VecSetSizes(*vec, stag->entriesGhost, PETSC_DETERMINE));
379:   PetscCall(VecSetType(*vec, dm->vectype));
380:   PetscCall(VecSetBlockSize(*vec, stag->entriesPerElement));
381:   PetscCall(VecSetDM(*vec, dm));
382:   PetscFunctionReturn(PETSC_SUCCESS);
383: }

385: /* Helper function to check for the limited situations for which interpolation
386:    and restriction functions are implemented */
387: static PetscErrorCode CheckTransferOperatorRequirements_Private(DM dmc, DM dmf)
388: {
389:   PetscInt dim, stencilWidthc, stencilWidthf, nf[DMSTAG_MAX_DIM], nc[DMSTAG_MAX_DIM], doff[DMSTAG_MAX_STRATA], dofc[DMSTAG_MAX_STRATA];

391:   PetscFunctionBegin;
392:   PetscCall(DMGetDimension(dmc, &dim));
393:   PetscCall(DMStagGetStencilWidth(dmc, &stencilWidthc));
394:   PetscCheck(stencilWidthc >= 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "DMCreateRestriction not implemented for coarse grid stencil width < 1");
395:   PetscCall(DMStagGetStencilWidth(dmf, &stencilWidthf));
396:   PetscCheck(stencilWidthf >= 1, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "DMCreateRestriction not implemented for fine grid stencil width < 1");
397:   PetscCall(DMStagGetLocalSizes(dmf, &nf[0], &nf[1], &nf[2]));
398:   PetscCall(DMStagGetLocalSizes(dmc, &nc[0], &nc[1], &nc[2]));
399:   for (PetscInt d = 0; d < dim; ++d)
400:     PetscCheck(nf[d] == 2 * nc[d], PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "No support for fine to coarse ratio other than 2 (it is %" PetscInt_FMT " to %" PetscInt_FMT " in dimension %" PetscInt_FMT ")", nf[d], nc[d], d);
401:   PetscCall(DMStagGetDOF(dmc, &dofc[0], &dofc[1], &dofc[2], &dofc[3]));
402:   PetscCall(DMStagGetDOF(dmf, &doff[0], &doff[1], &doff[2], &doff[3]));
403:   for (PetscInt d = 0; d < dim + 1; ++d)
404:     PetscCheck(dofc[d] == doff[d], PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "No support for different numbers of dof per stratum between coarse and fine DMStag objects: dof%" PetscInt_FMT " is %" PetscInt_FMT " (fine) but %" PetscInt_FMT "(coarse))", d, doff[d], dofc[d]);
405:   PetscFunctionReturn(PETSC_SUCCESS);
406: }

408: /* Since the interpolation uses MATMAIJ for dof > 0 we convert requests for non-MATAIJ baseded matrices to MATAIJ.
409:    This is a bit of a hack; the reason for it is partially because -dm_mat_type defines the
410:    matrix type for both the operator matrices and the interpolation matrices so that users
411:    can select matrix types of base MATAIJ for accelerators

413:    Note: The ConvertToAIJ() code below *has been copied from dainterp.c*! ConvertToAIJ() should perhaps be placed somewhere
414:    in mat/utils to avoid code duplication, but then the DMStag and DMDA code would need to include the private Mat headers.
415:    Since it is only used in two places, I have simply duplicated the code to avoid the need to exposure the private
416:    Mat routines in parts of DM. If we find a need for ConvertToAIJ() elsewhere, then we should consolidate it to one
417:    place in mat/utils.
418: */
419: static PetscErrorCode ConvertToAIJ(MatType intype, MatType *outtype)
420: {
421:   PetscInt    i;
422:   char const *types[3] = {MATAIJ, MATSEQAIJ, MATMPIAIJ};
423:   PetscBool   flg;

425:   PetscFunctionBegin;
426:   *outtype = MATAIJ;
427:   for (i = 0; i < 3; i++) {
428:     PetscCall(PetscStrbeginswith(intype, types[i], &flg));
429:     if (flg) {
430:       *outtype = intype;
431:       break;
432:     }
433:   }
434:   PetscFunctionReturn(PETSC_SUCCESS);
435: }

437: static PetscErrorCode DMCreateInterpolation_Stag(DM dmc, DM dmf, Mat *A, Vec *vec)
438: {
439:   PetscInt               dim, entriesf, entriesc, doff[DMSTAG_MAX_STRATA];
440:   ISLocalToGlobalMapping ltogmf, ltogmc;
441:   MatType                mattype;

443:   PetscFunctionBegin;
444:   PetscCall(CheckTransferOperatorRequirements_Private(dmc, dmf));

446:   PetscCall(DMStagGetEntries(dmf, &entriesf));
447:   PetscCall(DMStagGetEntries(dmc, &entriesc));
448:   PetscCall(DMGetLocalToGlobalMapping(dmf, &ltogmf));
449:   PetscCall(DMGetLocalToGlobalMapping(dmc, &ltogmc));

451:   PetscCall(MatCreate(PetscObjectComm((PetscObject)dmc), A));
452:   PetscCall(MatSetSizes(*A, entriesf, entriesc, PETSC_DECIDE, PETSC_DECIDE));
453:   PetscCall(ConvertToAIJ(dmc->mattype, &mattype));
454:   PetscCall(MatSetType(*A, mattype));
455:   PetscCall(MatSetLocalToGlobalMapping(*A, ltogmf, ltogmc));

457:   PetscCall(DMGetDimension(dmc, &dim));
458:   PetscCall(DMStagGetDOF(dmf, &doff[0], &doff[1], &doff[2], &doff[3]));
459:   if (dim == 1) {
460:     PetscCall(DMStagPopulateInterpolation1d_a_b_Private(dmc, dmf, *A));
461:   } else if (dim == 2) {
462:     if (doff[0] == 0) {
463:       PetscCall(DMStagPopulateInterpolation2d_0_a_b_Private(dmc, dmf, *A));
464:     } else
465:       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]);
466:   } else if (dim == 3) {
467:     if (doff[0] == 0 && doff[1] == 0) {
468:       PetscCall(DMStagPopulateInterpolation3d_0_0_a_b_Private(dmc, dmf, *A));
469:     } else
470:       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]);
471:   } else SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
472:   PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
473:   PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));

475:   if (vec) *vec = NULL;
476:   PetscFunctionReturn(PETSC_SUCCESS);
477: }

479: static PetscErrorCode DMCreateRestriction_Stag(DM dmc, DM dmf, Mat *A)
480: {
481:   PetscInt               dim, entriesf, entriesc, doff[DMSTAG_MAX_STRATA];
482:   ISLocalToGlobalMapping ltogmf, ltogmc;
483:   MatType                mattype;

485:   PetscFunctionBegin;
486:   PetscCall(CheckTransferOperatorRequirements_Private(dmc, dmf));

488:   PetscCall(DMStagGetEntries(dmf, &entriesf));
489:   PetscCall(DMStagGetEntries(dmc, &entriesc));
490:   PetscCall(DMGetLocalToGlobalMapping(dmf, &ltogmf));
491:   PetscCall(DMGetLocalToGlobalMapping(dmc, &ltogmc));

493:   PetscCall(MatCreate(PetscObjectComm((PetscObject)dmc), A));
494:   PetscCall(MatSetSizes(*A, entriesc, entriesf, PETSC_DECIDE, PETSC_DECIDE)); /* Note transpose wrt interpolation */
495:   PetscCall(ConvertToAIJ(dmc->mattype, &mattype));
496:   PetscCall(MatSetType(*A, mattype));
497:   PetscCall(MatSetLocalToGlobalMapping(*A, ltogmc, ltogmf)); /* Note transpose wrt interpolation */

499:   PetscCall(DMGetDimension(dmc, &dim));
500:   PetscCall(DMStagGetDOF(dmf, &doff[0], &doff[1], &doff[2], &doff[3]));
501:   if (dim == 1) {
502:     PetscCall(DMStagPopulateRestriction1d_a_b_Private(dmc, dmf, *A));
503:   } else if (dim == 2) {
504:     if (doff[0] == 0) {
505:       PetscCall(DMStagPopulateRestriction2d_0_a_b_Private(dmc, dmf, *A));
506:     } else
507:       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]);
508:   } else if (dim == 3) {
509:     if (doff[0] == 0 && doff[0] == 0) {
510:       PetscCall(DMStagPopulateRestriction3d_0_0_a_b_Private(dmc, dmf, *A));
511:     } else
512:       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]);
513:   } else SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);

515:   PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
516:   PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
517:   PetscFunctionReturn(PETSC_SUCCESS);
518: }

520: static PetscErrorCode DMCreateMatrix_Stag(DM dm, Mat *mat)
521: {
522:   MatType                mat_type;
523:   PetscBool              is_shell, is_aij;
524:   PetscInt               dim, entries;
525:   ISLocalToGlobalMapping ltogmap;

527:   PetscFunctionBegin;
528:   PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called after DMSetUp()");
529:   PetscCall(DMGetDimension(dm, &dim));
530:   PetscCall(DMGetMatType(dm, &mat_type));
531:   PetscCall(DMStagGetEntries(dm, &entries));
532:   PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), mat));
533:   PetscCall(MatSetSizes(*mat, entries, entries, PETSC_DETERMINE, PETSC_DETERMINE));
534:   PetscCall(MatSetType(*mat, mat_type));
535:   PetscCall(MatSetUp(*mat));
536:   PetscCall(DMGetLocalToGlobalMapping(dm, &ltogmap));
537:   PetscCall(MatSetLocalToGlobalMapping(*mat, ltogmap, ltogmap));
538:   PetscCall(MatSetDM(*mat, dm));

540:   /* Compare to similar and perhaps superior logic in DMCreateMatrix_DA, which creates
541:      the matrix first and then performs this logic by checking for preallocation functions */
542:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)*mat, MATAIJ, &is_aij));
543:   if (!is_aij) PetscCall(PetscObjectBaseTypeCompare((PetscObject)*mat, MATSEQAIJ, &is_aij));
544:   if (!is_aij) PetscCall(PetscObjectBaseTypeCompare((PetscObject)*mat, MATMPIAIJ, &is_aij));
545:   PetscCall(PetscStrcmp(mat_type, MATSHELL, &is_shell));
546:   if (is_aij) {
547:     Mat             preallocator;
548:     PetscInt        m, n;
549:     const PetscBool fill_with_zeros = PETSC_FALSE;

551:     PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), &preallocator));
552:     PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
553:     PetscCall(MatGetLocalSize(*mat, &m, &n));
554:     PetscCall(MatSetSizes(preallocator, m, n, PETSC_DECIDE, PETSC_DECIDE));
555:     PetscCall(MatSetLocalToGlobalMapping(preallocator, ltogmap, ltogmap));
556:     PetscCall(MatSetUp(preallocator));
557:     switch (dim) {
558:     case 1:
559:       PetscCall(DMCreateMatrix_Stag_1D_AIJ_Assemble(dm, preallocator));
560:       break;
561:     case 2:
562:       PetscCall(DMCreateMatrix_Stag_2D_AIJ_Assemble(dm, preallocator));
563:       break;
564:     case 3:
565:       PetscCall(DMCreateMatrix_Stag_3D_AIJ_Assemble(dm, preallocator));
566:       break;
567:     default:
568:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
569:     }
570:     PetscCall(MatPreallocatorPreallocate(preallocator, fill_with_zeros, *mat));
571:     PetscCall(MatDestroy(&preallocator));

573:     if (!dm->prealloc_only) {
574:       /* Bind to CPU before assembly, to prevent unnecessary copies of zero entries from CPU to GPU */
575:       PetscCall(MatBindToCPU(*mat, PETSC_TRUE));
576:       switch (dim) {
577:       case 1:
578:         PetscCall(DMCreateMatrix_Stag_1D_AIJ_Assemble(dm, *mat));
579:         break;
580:       case 2:
581:         PetscCall(DMCreateMatrix_Stag_2D_AIJ_Assemble(dm, *mat));
582:         break;
583:       case 3:
584:         PetscCall(DMCreateMatrix_Stag_3D_AIJ_Assemble(dm, *mat));
585:         break;
586:       default:
587:         SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
588:       }
589:       PetscCall(MatBindToCPU(*mat, PETSC_FALSE));
590:     }
591:   } else if (is_shell) {
592:     /* nothing more to do */
593:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for Mattype %s", mat_type);
594:   PetscFunctionReturn(PETSC_SUCCESS);
595: }

597: static PetscErrorCode DMGetCompatibility_Stag(DM dm, DM dm2, PetscBool *compatible, PetscBool *set)
598: {
599:   const DM_Stag *const stag  = (DM_Stag *)dm->data;
600:   const DM_Stag *const stag2 = (DM_Stag *)dm2->data;
601:   PetscInt             dim, dim2, i;
602:   MPI_Comm             comm;
603:   PetscMPIInt          sameComm;
604:   DMType               type2;
605:   PetscBool            sameType;

607:   PetscFunctionBegin;
608:   PetscCall(DMGetType(dm2, &type2));
609:   PetscCall(PetscStrcmp(DMSTAG, type2, &sameType));
610:   if (!sameType) {
611:     PetscCall(PetscInfo((PetscObject)dm, "DMStag compatibility check not implemented with DM of type %s\n", type2));
612:     *set = PETSC_FALSE;
613:     PetscFunctionReturn(PETSC_SUCCESS);
614:   }

616:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
617:   PetscCallMPI(MPI_Comm_compare(comm, PetscObjectComm((PetscObject)dm2), &sameComm));
618:   if (sameComm != MPI_IDENT) {
619:     PetscCall(PetscInfo((PetscObject)dm, "DMStag objects have different communicators: %" PETSC_INTPTR_T_FMT " != %" PETSC_INTPTR_T_FMT "\n", (PETSC_INTPTR_T)comm, (PETSC_INTPTR_T)PetscObjectComm((PetscObject)dm2)));
620:     *set = PETSC_FALSE;
621:     PetscFunctionReturn(PETSC_SUCCESS);
622:   }
623:   PetscCall(DMGetDimension(dm, &dim));
624:   PetscCall(DMGetDimension(dm2, &dim2));
625:   if (dim != dim2) {
626:     PetscCall(PetscInfo((PetscObject)dm, "DMStag objects have different dimensions\n"));
627:     *set        = PETSC_TRUE;
628:     *compatible = PETSC_FALSE;
629:     PetscFunctionReturn(PETSC_SUCCESS);
630:   }
631:   for (i = 0; i < dim; ++i) {
632:     if (stag->N[i] != stag2->N[i]) {
633:       PetscCall(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]));
634:       *set        = PETSC_TRUE;
635:       *compatible = PETSC_FALSE;
636:       PetscFunctionReturn(PETSC_SUCCESS);
637:     }
638:     if (stag->n[i] != stag2->n[i]) {
639:       PetscCall(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]));
640:       *set        = PETSC_TRUE;
641:       *compatible = PETSC_FALSE;
642:       PetscFunctionReturn(PETSC_SUCCESS);
643:     }
644:     if (stag->boundaryType[i] != stag2->boundaryType[i]) {
645:       PetscCall(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]]));
646:       *set        = PETSC_TRUE;
647:       *compatible = PETSC_FALSE;
648:       PetscFunctionReturn(PETSC_SUCCESS);
649:     }
650:   }
651:   /* Note: we include stencil type and width in the notion of compatibility, as this affects
652:      the "atlas" (local subdomains). This might be irritating in legitimate cases
653:      of wanting to transfer between two other-wise compatible DMs with different
654:      stencil characteristics. */
655:   if (stag->stencilType != stag2->stencilType) {
656:     PetscCall(PetscInfo((PetscObject)dm, "DMStag objects have different ghost stencil types: %s != %s\n", DMStagStencilTypes[stag->stencilType], DMStagStencilTypes[stag2->stencilType]));
657:     *set        = PETSC_TRUE;
658:     *compatible = PETSC_FALSE;
659:     PetscFunctionReturn(PETSC_SUCCESS);
660:   }
661:   if (stag->stencilWidth != stag2->stencilWidth) {
662:     PetscCall(PetscInfo((PetscObject)dm, "DMStag objects have different ghost stencil widths: %" PetscInt_FMT " != %" PetscInt_FMT "\n", stag->stencilWidth, stag->stencilWidth));
663:     *set        = PETSC_TRUE;
664:     *compatible = PETSC_FALSE;
665:     PetscFunctionReturn(PETSC_SUCCESS);
666:   }
667:   *set        = PETSC_TRUE;
668:   *compatible = PETSC_TRUE;
669:   PetscFunctionReturn(PETSC_SUCCESS);
670: }

672: static PetscErrorCode DMHasCreateInjection_Stag(DM dm, PetscBool *flg)
673: {
674:   PetscFunctionBegin;
676:   PetscAssertPointer(flg, 2);
677:   *flg = PETSC_FALSE;
678:   PetscFunctionReturn(PETSC_SUCCESS);
679: }

681: /*
682: Note there are several orderings in play here.
683: 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).
684: Also in all cases, only subdomains which are the last in their dimension have partial elements.

686: 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.
687: 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.
688: 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.
689: */

691: static PetscErrorCode DMLocalToGlobalBegin_Stag(DM dm, Vec l, InsertMode mode, Vec g)
692: {
693:   DM_Stag *const stag = (DM_Stag *)dm->data;

695:   PetscFunctionBegin;
696:   if (mode == ADD_VALUES) {
697:     PetscCall(VecScatterBegin(stag->gtol, l, g, mode, SCATTER_REVERSE));
698:   } else if (mode == INSERT_VALUES) {
699:     if (stag->ltog_injective) {
700:       PetscCall(VecScatterBegin(stag->ltog_injective, l, g, mode, SCATTER_FORWARD));
701:     } else {
702:       PetscCall(VecScatterBegin(stag->gtol, l, g, mode, SCATTER_REVERSE_LOCAL));
703:     }
704:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported InsertMode");
705:   PetscFunctionReturn(PETSC_SUCCESS);
706: }

708: static PetscErrorCode DMLocalToGlobalEnd_Stag(DM dm, Vec l, InsertMode mode, Vec g)
709: {
710:   DM_Stag *const stag = (DM_Stag *)dm->data;

712:   PetscFunctionBegin;
713:   if (mode == ADD_VALUES) {
714:     PetscCall(VecScatterEnd(stag->gtol, l, g, mode, SCATTER_REVERSE));
715:   } else if (mode == INSERT_VALUES) {
716:     if (stag->ltog_injective) {
717:       PetscCall(VecScatterEnd(stag->ltog_injective, l, g, mode, SCATTER_FORWARD));
718:     } else {
719:       PetscCall(VecScatterEnd(stag->gtol, l, g, mode, SCATTER_REVERSE_LOCAL));
720:     }
721:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported InsertMode");
722:   PetscFunctionReturn(PETSC_SUCCESS);
723: }

725: static PetscErrorCode DMGlobalToLocalBegin_Stag(DM dm, Vec g, InsertMode mode, Vec l)
726: {
727:   DM_Stag *const stag = (DM_Stag *)dm->data;

729:   PetscFunctionBegin;
730:   PetscCall(VecScatterBegin(stag->gtol, g, l, mode, SCATTER_FORWARD));
731:   PetscFunctionReturn(PETSC_SUCCESS);
732: }

734: static PetscErrorCode DMGlobalToLocalEnd_Stag(DM dm, Vec g, InsertMode mode, Vec l)
735: {
736:   DM_Stag *const stag = (DM_Stag *)dm->data;

738:   PetscFunctionBegin;
739:   PetscCall(VecScatterEnd(stag->gtol, g, l, mode, SCATTER_FORWARD));
740:   PetscFunctionReturn(PETSC_SUCCESS);
741: }

743: /*
744: If a stratum is active (non-zero dof), make it active in the coordinate DM.
745: */
746: static PetscErrorCode DMCreateCoordinateDM_Stag(DM dm, DM *dmc)
747: {
748:   DM_Stag *const stag = (DM_Stag *)dm->data;
749:   PetscInt       dim;
750:   PetscBool      isstag, isproduct;
751:   const char    *prefix;

753:   PetscFunctionBegin;

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

757:   PetscCall(DMGetDimension(dm, &dim));
758:   PetscCall(PetscStrcmp(stag->coordinateDMType, DMSTAG, &isstag));
759:   PetscCall(PetscStrcmp(stag->coordinateDMType, DMPRODUCT, &isproduct));
760:   if (isstag) {
761:     PetscCall(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));
762:   } else if (isproduct) {
763:     PetscCall(DMCreate(PETSC_COMM_WORLD, dmc));
764:     PetscCall(DMSetType(*dmc, DMPRODUCT));
765:     PetscCall(DMSetDimension(*dmc, dim));
766:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate DM type %s", stag->coordinateDMType);
767:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
768:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dmc, prefix));
769:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)*dmc, "cdm_"));
770:   PetscFunctionReturn(PETSC_SUCCESS);
771: }

773: static PetscErrorCode DMGetNeighbors_Stag(DM dm, PetscInt *nRanks, const PetscMPIInt *ranks[])
774: {
775:   DM_Stag *const stag = (DM_Stag *)dm->data;
776:   PetscInt       dim;

778:   PetscFunctionBegin;
779:   PetscCall(DMGetDimension(dm, &dim));
780:   switch (dim) {
781:   case 1:
782:     *nRanks = 3;
783:     break;
784:   case 2:
785:     *nRanks = 9;
786:     break;
787:   case 3:
788:     *nRanks = 27;
789:     break;
790:   default:
791:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Get neighbors not implemented for dim = %" PetscInt_FMT, dim);
792:   }
793:   *ranks = stag->neighbors;
794:   PetscFunctionReturn(PETSC_SUCCESS);
795: }

797: static PetscErrorCode DMView_Stag(DM dm, PetscViewer viewer)
798: {
799:   DM_Stag *const stag = (DM_Stag *)dm->data;
800:   PetscBool      isascii, viewAllRanks;
801:   PetscMPIInt    rank, size;
802:   PetscInt       dim, maxRanksToView, i;

804:   PetscFunctionBegin;
805:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
806:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
807:   PetscCall(DMGetDimension(dm, &dim));
808:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
809:   if (isascii) {
810:     PetscCall(PetscViewerASCIIPrintf(viewer, "Dimension: %" PetscInt_FMT "\n", dim));
811:     switch (dim) {
812:     case 1:
813:       PetscCall(PetscViewerASCIIPrintf(viewer, "Global size: %" PetscInt_FMT "\n", stag->N[0]));
814:       break;
815:     case 2:
816:       PetscCall(PetscViewerASCIIPrintf(viewer, "Global sizes: %" PetscInt_FMT " x %" PetscInt_FMT "\n", stag->N[0], stag->N[1]));
817:       PetscCall(PetscViewerASCIIPrintf(viewer, "Parallel decomposition: %" PetscInt_FMT " x %" PetscInt_FMT " ranks\n", stag->nRanks[0], stag->nRanks[1]));
818:       break;
819:     case 3:
820:       PetscCall(PetscViewerASCIIPrintf(viewer, "Global sizes: %" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT "\n", stag->N[0], stag->N[1], stag->N[2]));
821:       PetscCall(PetscViewerASCIIPrintf(viewer, "Parallel decomposition: %" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT " ranks\n", stag->nRanks[0], stag->nRanks[1], stag->nRanks[2]));
822:       break;
823:     default:
824:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "not implemented for dim==%" PetscInt_FMT, dim);
825:     }
826:     PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary ghosting:"));
827:     for (i = 0; i < dim; ++i) PetscCall(PetscViewerASCIIPrintf(viewer, " %s", DMBoundaryTypes[stag->boundaryType[i]]));
828:     PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
829:     PetscCall(PetscViewerASCIIPrintf(viewer, "Elementwise ghost stencil: %s", DMStagStencilTypes[stag->stencilType]));
830:     if (stag->stencilType != DMSTAG_STENCIL_NONE) {
831:       PetscCall(PetscViewerASCIIPrintf(viewer, ", width %" PetscInt_FMT "\n", stag->stencilWidth));
832:     } else {
833:       PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
834:     }
835:     PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per vertex (0D)\n", stag->dof[0]));
836:     if (dim == 3) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per edge (1D)\n", stag->dof[1]));
837:     if (dim > 1) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per face (%" PetscInt_FMT "D)\n", stag->dof[dim - 1], dim - 1));
838:     PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per element (%" PetscInt_FMT "D)\n", stag->dof[dim], dim));
839:     if (dm->coordinates[0].dm) PetscCall(PetscViewerASCIIPrintf(viewer, "Has coordinate DM\n"));
840:     maxRanksToView = 16;
841:     viewAllRanks   = (PetscBool)(size <= maxRanksToView);
842:     if (viewAllRanks) {
843:       PetscCall(PetscViewerASCIIPushSynchronized(viewer));
844:       switch (dim) {
845:       case 1:
846:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local elements : %" PetscInt_FMT " (%" PetscInt_FMT " with ghosts)\n", rank, stag->n[0], stag->nGhost[0]));
847:         break;
848:       case 2:
849:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Rank coordinates (%d,%d)\n", rank, stag->rank[0], stag->rank[1]));
850:         PetscCall(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]));
851:         break;
852:       case 3:
853:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Rank coordinates (%d,%d,%d)\n", rank, stag->rank[0], stag->rank[1], stag->rank[2]));
854:         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],
855:                                                      stag->n[2], stag->nGhost[0], stag->nGhost[1], stag->nGhost[2]));
856:         break;
857:       default:
858:         SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "not implemented for dim==%" PetscInt_FMT, dim);
859:       }
860:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local native entries: %" PetscInt_FMT "\n", rank, stag->entries));
861:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local entries total : %" PetscInt_FMT "\n", rank, stag->entriesGhost));
862:       PetscCall(PetscViewerFlush(viewer));
863:       PetscCall(PetscViewerASCIIPopSynchronized(viewer));
864:     } else {
865:       PetscCall(PetscViewerASCIIPrintf(viewer, "(Per-rank information omitted since >%" PetscInt_FMT " ranks used)\n", maxRanksToView));
866:     }
867:   }
868:   PetscFunctionReturn(PETSC_SUCCESS);
869: }

871: static PetscErrorCode DMSetFromOptions_Stag(DM dm, PetscOptionItems *PetscOptionsObject)
872: {
873:   DM_Stag *const stag = (DM_Stag *)dm->data;
874:   PetscInt       dim;

876:   PetscFunctionBegin;
877:   PetscCall(DMGetDimension(dm, &dim));
878:   PetscOptionsHeadBegin(PetscOptionsObject, "DMStag Options");
879:   PetscCall(PetscOptionsInt("-stag_grid_x", "Number of grid points in x direction", "DMStagSetGlobalSizes", stag->N[0], &stag->N[0], NULL));
880:   if (dim > 1) PetscCall(PetscOptionsInt("-stag_grid_y", "Number of grid points in y direction", "DMStagSetGlobalSizes", stag->N[1], &stag->N[1], NULL));
881:   if (dim > 2) PetscCall(PetscOptionsInt("-stag_grid_z", "Number of grid points in z direction", "DMStagSetGlobalSizes", stag->N[2], &stag->N[2], NULL));
882:   PetscCall(PetscOptionsInt("-stag_ranks_x", "Number of ranks in x direction", "DMStagSetNumRanks", stag->nRanks[0], &stag->nRanks[0], NULL));
883:   if (dim > 1) PetscCall(PetscOptionsInt("-stag_ranks_y", "Number of ranks in y direction", "DMStagSetNumRanks", stag->nRanks[1], &stag->nRanks[1], NULL));
884:   if (dim > 2) PetscCall(PetscOptionsInt("-stag_ranks_z", "Number of ranks in z direction", "DMStagSetNumRanks", stag->nRanks[2], &stag->nRanks[2], NULL));
885:   PetscCall(PetscOptionsInt("-stag_stencil_width", "Elementwise stencil width", "DMStagSetStencilWidth", stag->stencilWidth, &stag->stencilWidth, NULL));
886:   PetscCall(PetscOptionsEnum("-stag_stencil_type", "Elementwise stencil stype", "DMStagSetStencilType", DMStagStencilTypes, (PetscEnum)stag->stencilType, (PetscEnum *)&stag->stencilType, NULL));
887:   PetscCall(PetscOptionsEnum("-stag_boundary_type_x", "Treatment of (physical) boundaries in x direction", "DMStagSetBoundaryTypes", DMBoundaryTypes, (PetscEnum)stag->boundaryType[0], (PetscEnum *)&stag->boundaryType[0], NULL));
888:   PetscCall(PetscOptionsEnum("-stag_boundary_type_y", "Treatment of (physical) boundaries in y direction", "DMStagSetBoundaryTypes", DMBoundaryTypes, (PetscEnum)stag->boundaryType[1], (PetscEnum *)&stag->boundaryType[1], NULL));
889:   PetscCall(PetscOptionsEnum("-stag_boundary_type_z", "Treatment of (physical) boundaries in z direction", "DMStagSetBoundaryTypes", DMBoundaryTypes, (PetscEnum)stag->boundaryType[2], (PetscEnum *)&stag->boundaryType[2], NULL));
890:   PetscCall(PetscOptionsInt("-stag_dof_0", "Number of dof per 0-cell (vertex)", "DMStagSetDOF", stag->dof[0], &stag->dof[0], NULL));
891:   PetscCall(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));
892:   PetscCall(PetscOptionsInt("-stag_dof_2", "Number of dof per 2-cell (element in 2D, face in 3D)", "DMStagSetDOF", stag->dof[2], &stag->dof[2], NULL));
893:   PetscCall(PetscOptionsInt("-stag_dof_3", "Number of dof per 3-cell (element in 3D)", "DMStagSetDOF", stag->dof[3], &stag->dof[3], NULL));
894:   PetscOptionsHeadEnd();
895:   PetscFunctionReturn(PETSC_SUCCESS);
896: }

898: /*MC
899:   DMSTAG - `"stag"` - A `DM` object representing a "staggered grid" or a structured cell complex.

901:   Level: beginner

903:   Notes:
904:   This implementation parallels the `DMDA` implementation in many ways, but allows degrees of freedom
905:   to be associated with all "strata" in a logically-rectangular grid.

907:   Each stratum can be characterized by the dimension of the entities ("points", to borrow the `DMPLEX`
908:   terminology), from 0- to 3-dimensional.

910:   In some cases this numbering is used directly, for example with `DMStagGetDOF()`.
911:   To allow easier reading and to some extent more similar code between different-dimensional implementations
912:   of the same problem, we associate canonical names for each type of point, for each dimension of DMStag.

914:   * 1-dimensional `DMSTAG` objects have vertices (0D) and elements (1D).
915:   * 2-dimensional `DMSTAG` objects have vertices (0D), faces (1D), and elements (2D).
916:   * 3-dimensional `DMSTAG` objects have vertices (0D), edges (1D), faces (2D), and elements (3D).

918:   This naming is reflected when viewing a `DMSTAG` object with `DMView()`, and in forming
919:   convenient options prefixes when creating a decomposition with `DMCreateFieldDecomposition()`.

921: .seealso: [](ch_stag), `DM`, `DMPRODUCT`, `DMDA`, `DMPLEX`, `DMStagCreate1d()`, `DMStagCreate2d()`, `DMStagCreate3d()`, `DMType`, `DMCreate()`,
922:           `DMSetType()`, `DMStagVecSplitToDMDA()`
923: M*/

925: PETSC_EXTERN PetscErrorCode DMCreate_Stag(DM dm)
926: {
927:   DM_Stag *stag;
928:   PetscInt i, dim;

930:   PetscFunctionBegin;
931:   PetscAssertPointer(dm, 1);
932:   PetscCall(PetscNew(&stag));
933:   dm->data = stag;

935:   stag->gtol           = NULL;
936:   stag->ltog_injective = NULL;
937:   for (i = 0; i < DMSTAG_MAX_STRATA; ++i) stag->dof[i] = 0;
938:   for (i = 0; i < DMSTAG_MAX_DIM; ++i) stag->l[i] = NULL;
939:   stag->stencilType  = DMSTAG_STENCIL_NONE;
940:   stag->stencilWidth = 0;
941:   for (i = 0; i < DMSTAG_MAX_DIM; ++i) stag->nRanks[i] = -1;
942:   stag->coordinateDMType = NULL;

944:   PetscCall(DMGetDimension(dm, &dim));
945:   PetscCheck(dim == 1 || dim == 2 || dim == 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMSetDimension() must be called to set a dimension with value 1, 2, or 3");

947:   PetscCall(PetscMemzero(dm->ops, sizeof(*(dm->ops))));
948:   dm->ops->createcoordinatedm  = DMCreateCoordinateDM_Stag;
949:   dm->ops->createglobalvector  = DMCreateGlobalVector_Stag;
950:   dm->ops->createlocalvector   = DMCreateLocalVector_Stag;
951:   dm->ops->creatematrix        = DMCreateMatrix_Stag;
952:   dm->ops->hascreateinjection  = DMHasCreateInjection_Stag;
953:   dm->ops->refine              = DMRefine_Stag;
954:   dm->ops->coarsen             = DMCoarsen_Stag;
955:   dm->ops->createinterpolation = DMCreateInterpolation_Stag;
956:   dm->ops->createrestriction   = DMCreateRestriction_Stag;
957:   dm->ops->destroy             = DMDestroy_Stag;
958:   dm->ops->getneighbors        = DMGetNeighbors_Stag;
959:   dm->ops->globaltolocalbegin  = DMGlobalToLocalBegin_Stag;
960:   dm->ops->globaltolocalend    = DMGlobalToLocalEnd_Stag;
961:   dm->ops->localtoglobalbegin  = DMLocalToGlobalBegin_Stag;
962:   dm->ops->localtoglobalend    = DMLocalToGlobalEnd_Stag;
963:   dm->ops->setfromoptions      = DMSetFromOptions_Stag;
964:   switch (dim) {
965:   case 1:
966:     dm->ops->setup = DMSetUp_Stag_1d;
967:     break;
968:   case 2:
969:     dm->ops->setup = DMSetUp_Stag_2d;
970:     break;
971:   case 3:
972:     dm->ops->setup = DMSetUp_Stag_3d;
973:     break;
974:   default:
975:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
976:   }
977:   dm->ops->clone                    = DMClone_Stag;
978:   dm->ops->view                     = DMView_Stag;
979:   dm->ops->getcompatibility         = DMGetCompatibility_Stag;
980:   dm->ops->createfielddecomposition = DMCreateFieldDecomposition_Stag;
981:   PetscFunctionReturn(PETSC_SUCCESS);
982: }