Actual source code: ex49.c
1: static char help[] = "Tests dof numberings for external integrators such as LibCEED.\n\n";
3: #include <petscdmplex.h>
4: #include <petscds.h>
6: typedef struct {
7: PetscBool useFE;
8: PetscInt check_face;
9: PetscBool closure_tensor;
10: PetscBool bc;
11: } AppCtx;
13: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
14: {
15: PetscFunctionBeginUser;
16: options->useFE = PETSC_TRUE;
17: options->check_face = 1;
18: options->closure_tensor = PETSC_FALSE;
19: options->bc = PETSC_FALSE;
20: PetscOptionsBegin(comm, "", "Dof Ordering Options", "DMPLEX");
21: PetscCall(PetscOptionsBool("-use_fe", "Use FE or FV discretization", "ex49.c", options->useFE, &options->useFE, NULL));
22: PetscCall(PetscOptionsInt("-check_face", "Face set to report on", "ex49.c", options->check_face, &options->check_face, NULL));
23: PetscCall(PetscOptionsBool("-closure_tensor", "Use DMPlexSetClosurePermutationTensor()", "ex49.c", options->closure_tensor, &options->closure_tensor, NULL));
24: PetscCall(PetscOptionsBool("-bc", "Add an essential boundary condition on field 0 over marker 1, so the section carries constrained dofs", "ex49.c", options->bc, &options->bc, NULL));
25: PetscOptionsEnd();
26: PetscFunctionReturn(PETSC_SUCCESS);
27: }
29: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
30: {
31: PetscFunctionBeginUser;
32: PetscCall(DMCreate(comm, dm));
33: PetscCall(DMSetType(*dm, DMPLEX));
34: PetscCall(DMSetFromOptions(*dm));
35: PetscCall(DMSetApplicationContext(*dm, user));
36: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
37: PetscFunctionReturn(PETSC_SUCCESS);
38: }
40: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
41: {
42: DM cdm = dm;
43: PetscInt dim;
45: PetscFunctionBeginUser;
46: PetscCall(DMGetDimension(dm, &dim));
47: if (user->useFE) {
48: PetscFE fe;
49: DMPolytopeType ct;
50: PetscInt cStart;
52: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
53: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
54: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, NULL, PETSC_DETERMINE, &fe));
55: PetscCall(PetscObjectSetName((PetscObject)fe, "scalar"));
56: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
57: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe));
58: PetscCall(PetscFEDestroy(&fe));
59: } else {
60: PetscFV fv;
62: PetscCall(PetscFVCreate(PETSC_COMM_SELF, &fv));
63: PetscCall(PetscFVSetType(fv, PETSCFVLEASTSQUARES));
64: PetscCall(PetscFVSetNumComponents(fv, dim));
65: PetscCall(PetscFVSetSpatialDimension(fv, dim));
66: PetscCall(PetscFVSetFromOptions(fv));
67: PetscCall(PetscFVSetUp(fv));
68: PetscCall(PetscObjectSetName((PetscObject)fv, "vector"));
69: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fv));
70: PetscCall(PetscFVDestroy(&fv));
71: }
72: PetscCall(DMCreateDS(dm));
73: if (user->bc) {
74: DMLabel label;
75: PetscInt id = 1;
77: PetscCall(DMGetLabel(dm, "marker", &label));
78: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, NULL, NULL, NULL, NULL));
79: }
80: while (cdm) {
81: PetscCall(DMCopyDisc(dm, cdm));
82: PetscCall(DMGetCoarseDM(cdm, &cdm));
83: }
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: static PetscErrorCode CheckOffsets(DM dm, AppCtx *user, const char *domain_name, PetscInt label_value, PetscInt height)
88: {
89: const char *height_name[] = {"cells", "faces"};
90: DMLabel domain_label = NULL;
91: DM cdm;
92: PetscInt Nf;
93: ISLocalToGlobalMapping ltog;
95: PetscFunctionBeginUser;
96: if (domain_name) PetscCall(DMGetLabel(dm, domain_name, &domain_label));
97: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "## %s: '%s' {%" PetscInt_FMT "}%s\n", height_name[height], domain_name ? domain_name : "default", label_value, domain_name && !domain_label ? " (null label)" : ""));
98: if (domain_name && !domain_label) PetscFunctionReturn(PETSC_SUCCESS);
99: if (user->closure_tensor) PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL));
100: // Offsets for cell closures
101: PetscCall(DMGetNumFields(dm, &Nf));
102: for (PetscInt f = 0; f < Nf; ++f) {
103: PetscObject obj;
104: PetscClassId id;
105: char name[PETSC_MAX_PATH_LEN];
107: PetscCall(DMGetField(dm, f, NULL, &obj));
108: PetscCall(PetscObjectGetClassId(obj, &id));
109: if (id == PETSCFE_CLASSID) {
110: IS offIS;
111: PetscInt *offsets, Ncell, Ncl, Nc, n;
113: PetscCall(DMPlexGetLocalOffsets(dm, domain_label, label_value, height, f, &Ncell, &Ncl, &Nc, &n, &offsets));
114: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Ncell * Ncl, offsets, PETSC_OWN_POINTER, &offIS));
115: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Field %" PetscInt_FMT " Offsets", f));
116: PetscCall(PetscObjectSetName((PetscObject)offIS, name));
117: PetscCall(ISViewFromOptions(offIS, NULL, "-offsets_view"));
118: PetscCall(ISDestroy(&offIS));
119: } else if (id == PETSCFV_CLASSID) {
120: IS offIS;
121: PetscInt *offsets, *offsetsNeg, *offsetsPos, Nface, Nc, n, i = 0;
123: PetscCall(DMPlexGetLocalOffsetsSupport(dm, domain_label, label_value, &Nface, &Nc, &n, &offsetsNeg, &offsetsPos));
124: PetscCall(PetscMalloc1(Nface * Nc * 2, &offsets));
125: for (PetscInt f = 0; f < Nface; ++f) {
126: for (PetscInt c = 0; c < Nc; ++c) offsets[i++] = offsetsNeg[f] + c;
127: for (PetscInt c = 0; c < Nc; ++c) offsets[i++] = offsetsPos[f] + c;
128: }
129: PetscCheck(i == Nface * Nc * 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total offsets %" PetscInt_FMT " != %" PetscInt_FMT, i, Nface * Nc * 2);
130: PetscCall(PetscFree(offsetsNeg));
131: PetscCall(PetscFree(offsetsPos));
132: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Nface * Nc * 2, offsets, PETSC_OWN_POINTER, &offIS));
133: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Field %" PetscInt_FMT " Offsets", f));
134: PetscCall(PetscObjectSetName((PetscObject)offIS, name));
135: PetscCall(ISViewFromOptions(offIS, NULL, "-offsets_view"));
136: PetscCall(ISDestroy(&offIS));
137: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Unrecognized type for DM field %" PetscInt_FMT, f);
138: }
139: PetscCall(DMGetLocalToGlobalMapping(dm, <og));
140: PetscCall(ISLocalToGlobalMappingViewFromOptions(ltog, NULL, "-ltog_view"));
142: // Offsets for coordinates
143: {
144: Vec X;
145: PetscSection s;
146: const PetscScalar *x;
147: const char *cname;
148: PetscInt cdim, *offsets, Ncell, Ncl, Nc, n;
149: PetscBool isDG = PETSC_FALSE;
151: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
152: if (!cdm) {
153: PetscCall(DMGetCoordinateDM(dm, &cdm));
154: cname = "Coordinates";
155: PetscCall(DMGetCoordinatesLocal(dm, &X));
156: } else {
157: isDG = PETSC_TRUE;
158: cname = "DG Coordinates";
159: PetscCall(DMGetCellCoordinatesLocal(dm, &X));
160: }
161: if (isDG && height) PetscFunctionReturn(PETSC_SUCCESS);
162: if (domain_name) PetscCall(DMGetLabel(cdm, domain_name, &domain_label));
163: if (user->closure_tensor) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL));
164: PetscCall(DMPlexGetLocalOffsets(cdm, domain_label, label_value, height, 0, &Ncell, &Ncl, &Nc, &n, &offsets));
165: PetscCall(DMGetCoordinateDim(dm, &cdim));
166: PetscCheck(Nc == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Geometric dimension %" PetscInt_FMT " should be %" PetscInt_FMT, Nc, cdim);
167: PetscCall(DMGetLocalSection(cdm, &s));
168: PetscCall(VecGetArrayRead(X, &x));
169: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s by element in %s order\n", cname, user->closure_tensor ? "tensor" : "bfs"));
170: for (PetscInt c = 0; c < Ncell; ++c) {
171: for (PetscInt v = 0; v < Ncl; ++v) {
172: PetscInt off = offsets[c * Ncl + v], dgdof;
173: const PetscScalar *vx = &x[off];
175: if (isDG) {
176: PetscCall(PetscSectionGetDof(s, c, &dgdof));
177: PetscCheck(Ncl * Nc == dgdof, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Offset size %" PetscInt_FMT " should be %" PetscInt_FMT, Ncl * Nc, dgdof);
178: }
179: switch (cdim) {
180: case 1:
181: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT "] %" PetscInt_FMT " <-- %2" PetscInt_FMT " (% 4.2f)\n", c, v, off, (double)PetscRealPart(vx[0])));
182: break;
183: case 2:
184: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT "] %" PetscInt_FMT " <-- %2" PetscInt_FMT " (% 4.2f, % 4.2f)\n", c, v, off, (double)PetscRealPart(vx[0]), (double)PetscRealPart(vx[1])));
185: break;
186: case 3:
187: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT "] %" PetscInt_FMT " <-- %2" PetscInt_FMT " (% 4.2f, % 4.2f, % 4.2f)\n", c, v, off, (double)PetscRealPart(vx[0]), (double)PetscRealPart(vx[1]), (double)PetscRealPart(vx[2])));
188: }
189: }
190: }
191: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout));
192: PetscCall(VecRestoreArrayRead(X, &x));
193: PetscCall(PetscFree(offsets));
194: PetscCall(DMGetLocalToGlobalMapping(cdm, <og));
195: PetscCall(ISLocalToGlobalMappingViewFromOptions(ltog, NULL, "-coord_ltog_view"));
196: {
197: DM clonedm;
198: Vec cloneX, X;
199: PetscInt clone_num_x, num_x;
200: const PetscScalar *clonex, *x;
202: PetscCall(DMClone(dm, &clonedm));
203: { // Force recreation of local coordinate vector
204: Vec X_global;
206: PetscCall(DMGetCoordinates(dm, &X_global));
207: PetscCall(DMSetCoordinates(clonedm, X_global));
208: }
209: PetscCall(DMGetCoordinatesLocal(dm, &X));
210: PetscCall(DMGetCoordinatesLocal(clonedm, &cloneX));
211: PetscCall(VecGetLocalSize(X, &num_x));
212: PetscCall(VecGetLocalSize(cloneX, &clone_num_x));
213: PetscCheck(num_x == clone_num_x, PETSC_COMM_WORLD, PETSC_ERR_ARG_SIZ, "Cloned DM coordinate size (%" PetscInt_FMT ") different from original DM coordinate size (%" PetscInt_FMT ")", clone_num_x, num_x);
215: PetscCall(VecGetArrayRead(X, &x));
216: PetscCall(VecGetArrayRead(cloneX, &clonex));
218: for (PetscInt i = 0; i < num_x; i++) {
219: PetscCheck(PetscIsCloseAtTolScalar(x[i], clonex[i], 1e-13, 1e-13), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Original coordinate (%4.2f) and cloned coordinate (%4.2f) are different", (double)PetscRealPart(x[i]), (double)PetscRealPart(clonex[i]));
220: }
222: PetscCall(VecRestoreArrayRead(X, &x));
223: PetscCall(VecRestoreArrayRead(cloneX, &clonex));
224: PetscCall(DMDestroy(&clonedm));
225: }
226: }
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: int main(int argc, char **argv)
231: {
232: DM dm;
233: AppCtx user;
234: PetscInt depth;
235: PetscBool flg;
237: PetscFunctionBeginUser;
238: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
239: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
240: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
241: PetscCall(SetupDiscretization(dm, &user));
242: PetscCall(CheckOffsets(dm, &user, NULL, 0, 0));
243: PetscCall(DMPlexGetDepth(dm, &depth));
244: if (depth > 1) PetscCall(CheckOffsets(dm, &user, "Face Sets", user.check_face, 1));
246: PetscCall(PetscOptionsHasName(NULL, NULL, "-view_mat", &flg));
247: if (flg) {
248: Mat A;
249: PetscCall(DMCreateMatrix(dm, &A));
250: PetscCall(MatViewFromOptions(A, NULL, "-view_mat"));
251: PetscCall(MatDestroy(&A));
252: }
253: PetscCall(PetscOptionsHasName(NULL, NULL, "-test_section_reset", &flg));
254: if (flg) {
255: PetscSection sec;
256: ISLocalToGlobalMapping ltog, ltognew;
258: /* Setting a section must invalidate a previously built local-to-global mapping. Hold a
259: reference to the old mapping so a rebuilt one cannot alias its address. */
260: PetscCall(DMGetLocalToGlobalMapping(dm, <og));
261: PetscCall(PetscObjectReference((PetscObject)ltog));
262: PetscCall(DMGetLocalSection(dm, &sec));
263: PetscCall(DMSetLocalSection(dm, sec));
264: PetscCall(DMGetLocalToGlobalMapping(dm, <ognew));
265: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "l2g map %s by DMSetLocalSection()\n", ltognew != ltog ? "invalidated" : "NOT invalidated"));
266: PetscCall(ISLocalToGlobalMappingDestroy(<og));
267: ltog = ltognew;
268: PetscCall(PetscObjectReference((PetscObject)ltog));
269: PetscCall(DMGetGlobalSection(dm, &sec));
270: PetscCall(DMSetGlobalSection(dm, sec));
271: PetscCall(DMGetLocalToGlobalMapping(dm, <ognew));
272: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "l2g map %s by DMSetGlobalSection()\n", ltognew != ltog ? "invalidated" : "NOT invalidated"));
273: PetscCall(ISLocalToGlobalMappingDestroy(<og));
274: }
275: PetscCall(DMDestroy(&dm));
276: PetscCall(PetscFinalize());
277: return 0;
278: }
280: /*TEST
282: test:
283: suffix: 0
284: requires: triangle
285: args: -dm_refine 1 -petscspace_degree 1 -dm_view -offsets_view
287: test:
288: suffix: 1
289: args: -dm_plex_simplex 0 -dm_plex_box_bd periodic,none -dm_plex_box_faces 3,3 -dm_sparse_localize 0 -petscspace_degree 1 \
290: -dm_view -offsets_view
292: test:
293: suffix: cg_2d
294: args: -dm_plex_simplex 0 -dm_plex_box_bd none,none -dm_plex_box_faces 3,3 -petscspace_degree 1 \
295: -dm_view -offsets_view
297: # A chart permutation makes the local storage slot of a point differ from a running count over
298: # the chart, so the local-to-global map must be indexed by the local section's offsets. Serial
299: # and unconstrained, the map must still come out as the identity.
300: test:
301: suffix: reorder_ltog
302: args: -dm_plex_simplex 0 -dm_plex_box_faces 2,2 -petscspace_degree 1 \
303: -dm_reorder_section -dm_reorder_section_type reverse -offsets_view -ltog_view -test_section_reset
305: # The parallel companion pins a non-identity map: ghost points have local slots that differ from
306: # their global indices, so indexing by the wrong offset shows up here. No -offsets_view: those
307: # per-rank views interleave nondeterministically on stdout
308: test:
309: suffix: reorder_ltog_par
310: nsize: 2
311: args: -dm_plex_simplex 0 -dm_plex_box_faces 2,2 -petscspace_degree 1 -petscpartitioner_type simple \
312: -dm_reorder_section -dm_reorder_section_type reverse -ltog_view
314: # Adds constrained dofs (essential bc on field 0 only, so points mix constrained and unconstrained
315: # fields), combined with the chart permutation and parallel ghost points: the map must encode the
316: # constrained dofs as negative entries at the permuted local slots, while unconstrained ghost dofs
317: # carry their positive owner-side global indices
318: test:
319: suffix: reorder_ltog_bc_par
320: nsize: 2
321: args: -dm_plex_simplex 0 -dm_plex_box_faces 2,2 -petscspace_degree 1 -petscpartitioner_type simple \
322: -dm_reorder_section -dm_reorder_section_type reverse -bc -ltog_view
324: test:
325: suffix: 1d_sfc
326: args: -dm_plex_simplex 0 -dm_plex_dim 1 -dm_plex_shape zbox -dm_plex_box_faces 3 1 -dm_view -coord_ltog_view
328: test:
329: suffix: 1d_sfc_periodic
330: args: -dm_plex_simplex 0 -dm_plex_dim 1 -dm_plex_shape zbox -dm_plex_box_faces 3 1 -dm_view -coord_ltog_view -petscspace_degree 1 -view_mat -dm_plex_box_bd periodic
332: test:
333: suffix: 2d_sfc
334: nsize: 2
335: args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_shape zbox -dm_plex_box_faces 4,3 -dm_distribute 0 -petscspace_degree 1 -dm_view
337: test:
338: suffix: 2d_sfc_periodic
339: nsize: 2
340: args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_shape zbox -dm_plex_box_faces 4,3 -dm_distribute 0 -petscspace_degree 1 -dm_plex_box_bd periodic,none -dm_view ::ascii_info_detail -view_mat
342: test:
343: suffix: 2d_sfc_periodic_mat
344: args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_shape zbox -dm_plex_box_faces 1,1 -dm_distribute 0 -petscspace_degree 2 -dm_plex_box_bd periodic,periodic -dm_view ::ascii_info_detail -view_mat
346: testset:
347: args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_shape zbox -dm_plex_box_faces 3,2 -petscspace_degree 1 -dm_view ::ascii_info_detail -closure_tensor
348: nsize: 2
349: test:
350: suffix: 2d_sfc_periodic_stranded
351: args: -dm_distribute 0 -dm_plex_box_bd none,periodic
352: test:
353: suffix: 2d_sfc_periodic_stranded_dist
354: args: -dm_distribute 1 -petscpartitioner_type simple -dm_plex_box_bd none,periodic
355: test:
356: suffix: 2d_sfc_biperiodic_stranded
357: args: -dm_distribute 0 -dm_plex_box_bd periodic,periodic
358: test:
359: suffix: 2d_sfc_biperiodic_stranded_dist
360: args: -dm_distribute 1 -petscpartitioner_type simple -dm_plex_box_bd periodic,periodic
361: test:
362: suffix: 2d_sfc_biperiodic_stranded_dist_box_label
363: args: -dm_distribute 1 -petscpartitioner_type simple -dm_plex_box_label_bd periodic,periodic -dm_plex_box_label
365: test:
366: suffix: fv_0
367: requires: triangle
368: args: -dm_refine 1 -use_fe 0 -dm_view -offsets_view
370: TEST*/