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: } AppCtx;
12: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
13: {
14: PetscFunctionBeginUser;
15: options->useFE = PETSC_TRUE;
16: options->check_face = 1;
17: options->closure_tensor = PETSC_FALSE;
18: PetscOptionsBegin(comm, "", "Dof Ordering Options", "DMPLEX");
19: PetscCall(PetscOptionsBool("-use_fe", "Use FE or FV discretization", "ex49.c", options->useFE, &options->useFE, NULL));
20: PetscCall(PetscOptionsInt("-check_face", "Face set to report on", "ex49.c", options->check_face, &options->check_face, NULL));
21: PetscCall(PetscOptionsBool("-closure_tensor", "Use DMPlexSetClosurePermutationTensor()", "ex49.c", options->closure_tensor, &options->closure_tensor, NULL));
22: PetscOptionsEnd();
23: PetscFunctionReturn(PETSC_SUCCESS);
24: }
26: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
27: {
28: PetscFunctionBeginUser;
29: PetscCall(DMCreate(comm, dm));
30: PetscCall(DMSetType(*dm, DMPLEX));
31: PetscCall(DMSetFromOptions(*dm));
32: PetscCall(DMSetApplicationContext(*dm, user));
33: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
38: {
39: DM cdm = dm;
40: PetscInt dim;
42: PetscFunctionBeginUser;
43: PetscCall(DMGetDimension(dm, &dim));
44: if (user->useFE) {
45: PetscFE fe;
46: DMPolytopeType ct;
47: PetscInt cStart;
49: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
50: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
51: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, NULL, -1, &fe));
52: PetscCall(PetscObjectSetName((PetscObject)fe, "scalar"));
53: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
54: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe));
55: PetscCall(PetscFEDestroy(&fe));
56: } else {
57: PetscFV fv;
59: PetscCall(PetscFVCreate(PETSC_COMM_SELF, &fv));
60: PetscCall(PetscFVSetType(fv, PETSCFVLEASTSQUARES));
61: PetscCall(PetscFVSetNumComponents(fv, dim));
62: PetscCall(PetscFVSetSpatialDimension(fv, dim));
63: PetscCall(PetscFVSetFromOptions(fv));
64: PetscCall(PetscFVSetUp(fv));
65: PetscCall(PetscObjectSetName((PetscObject)fv, "vector"));
66: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fv));
67: PetscCall(PetscFVDestroy(&fv));
68: }
69: PetscCall(DMCreateDS(dm));
70: while (cdm) {
71: PetscCall(DMCopyDisc(dm, cdm));
72: PetscCall(DMGetCoarseDM(cdm, &cdm));
73: }
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: static PetscErrorCode CheckOffsets(DM dm, AppCtx *user, const char *domain_name, PetscInt label_value, PetscInt height)
78: {
79: const char *height_name[] = {"cells", "faces"};
80: DMLabel domain_label = NULL;
81: DM cdm;
82: PetscInt Nf, f;
83: ISLocalToGlobalMapping ltog;
85: PetscFunctionBeginUser;
86: if (domain_name) PetscCall(DMGetLabel(dm, domain_name, &domain_label));
87: 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)" : ""));
88: if (domain_name && !domain_label) PetscFunctionReturn(PETSC_SUCCESS);
89: if (user->closure_tensor) PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL));
90: // Offsets for cell closures
91: PetscCall(DMGetNumFields(dm, &Nf));
92: for (f = 0; f < Nf; ++f) {
93: PetscObject obj;
94: PetscClassId id;
95: char name[PETSC_MAX_PATH_LEN];
97: PetscCall(DMGetField(dm, f, NULL, &obj));
98: PetscCall(PetscObjectGetClassId(obj, &id));
99: if (id == PETSCFE_CLASSID) {
100: IS offIS;
101: PetscInt *offsets, Ncell, Ncl, Nc, n;
103: PetscCall(DMPlexGetLocalOffsets(dm, domain_label, label_value, height, f, &Ncell, &Ncl, &Nc, &n, &offsets));
104: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Ncell * Ncl, offsets, PETSC_OWN_POINTER, &offIS));
105: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Field %" PetscInt_FMT " Offsets", f));
106: PetscCall(PetscObjectSetName((PetscObject)offIS, name));
107: PetscCall(ISViewFromOptions(offIS, NULL, "-offsets_view"));
108: PetscCall(ISDestroy(&offIS));
109: } else if (id == PETSCFV_CLASSID) {
110: IS offIS;
111: PetscInt *offsets, *offsetsNeg, *offsetsPos, Nface, Nc, n, i = 0;
113: PetscCall(DMPlexGetLocalOffsetsSupport(dm, domain_label, label_value, &Nface, &Nc, &n, &offsetsNeg, &offsetsPos));
114: PetscCall(PetscMalloc1(Nface * Nc * 2, &offsets));
115: for (PetscInt f = 0; f < Nface; ++f) {
116: for (PetscInt c = 0; c < Nc; ++c) offsets[i++] = offsetsNeg[f] + c;
117: for (PetscInt c = 0; c < Nc; ++c) offsets[i++] = offsetsPos[f] + c;
118: }
119: PetscCheck(i == Nface * Nc * 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total offsets %" PetscInt_FMT " != %" PetscInt_FMT, i, Nface * Nc * 2);
120: PetscCall(PetscFree(offsetsNeg));
121: PetscCall(PetscFree(offsetsPos));
122: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Nface * Nc * 2, offsets, PETSC_OWN_POINTER, &offIS));
123: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Field %" PetscInt_FMT " Offsets", f));
124: PetscCall(PetscObjectSetName((PetscObject)offIS, name));
125: PetscCall(ISViewFromOptions(offIS, NULL, "-offsets_view"));
126: PetscCall(ISDestroy(&offIS));
127: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Unrecognized type for DM field %" PetscInt_FMT, f);
128: }
129: PetscCall(DMGetLocalToGlobalMapping(dm, <og));
130: PetscCall(ISLocalToGlobalMappingViewFromOptions(ltog, NULL, "-ltog_view"));
132: // Offsets for coordinates
133: {
134: Vec X;
135: PetscSection s;
136: const PetscScalar *x;
137: const char *cname;
138: PetscInt cdim, *offsets, Ncell, Ncl, Nc, n;
139: PetscBool isDG = PETSC_FALSE;
141: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
142: if (!cdm) {
143: PetscCall(DMGetCoordinateDM(dm, &cdm));
144: cname = "Coordinates";
145: PetscCall(DMGetCoordinatesLocal(dm, &X));
146: } else {
147: isDG = PETSC_TRUE;
148: cname = "DG Coordinates";
149: PetscCall(DMGetCellCoordinatesLocal(dm, &X));
150: }
151: if (isDG && height) PetscFunctionReturn(PETSC_SUCCESS);
152: if (domain_name) PetscCall(DMGetLabel(cdm, domain_name, &domain_label));
153: if (user->closure_tensor) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL));
154: PetscCall(DMPlexGetLocalOffsets(cdm, domain_label, label_value, height, 0, &Ncell, &Ncl, &Nc, &n, &offsets));
155: PetscCall(DMGetCoordinateDim(dm, &cdim));
156: PetscCheck(Nc == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Geometric dimension %" PetscInt_FMT " should be %" PetscInt_FMT, Nc, cdim);
157: PetscCall(DMGetLocalSection(cdm, &s));
158: PetscCall(VecGetArrayRead(X, &x));
159: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s by element in %s order\n", cname, user->closure_tensor ? "tensor" : "bfs"));
160: for (PetscInt c = 0; c < Ncell; ++c) {
161: for (PetscInt v = 0; v < Ncl; ++v) {
162: PetscInt off = offsets[c * Ncl + v], dgdof;
163: const PetscScalar *vx = &x[off];
165: if (isDG) {
166: PetscCall(PetscSectionGetDof(s, c, &dgdof));
167: PetscCheck(Ncl * Nc == dgdof, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Offset size %" PetscInt_FMT " should be %" PetscInt_FMT, Ncl * Nc, dgdof);
168: }
169: switch (cdim) {
170: case 1:
171: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT "] %" PetscInt_FMT " <-- %2" PetscInt_FMT " (% 4.2f)\n", c, v, off, (double)PetscRealPart(vx[0])));
172: break;
173: case 2:
174: 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])));
175: break;
176: case 3:
177: 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])));
178: }
179: }
180: }
181: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout));
182: PetscCall(VecRestoreArrayRead(X, &x));
183: PetscCall(PetscFree(offsets));
184: PetscCall(DMGetLocalToGlobalMapping(cdm, <og));
185: PetscCall(ISLocalToGlobalMappingViewFromOptions(ltog, NULL, "-coord_ltog_view"));
186: {
187: DM clonedm;
188: Vec cloneX, X;
189: PetscInt clone_num_x, num_x;
190: const PetscScalar *clonex, *x;
192: PetscCall(DMClone(dm, &clonedm));
193: { // Force recreation of local coordinate vector
194: Vec X_global;
196: PetscCall(DMGetCoordinates(dm, &X_global));
197: PetscCall(DMSetCoordinates(clonedm, X_global));
198: }
199: PetscCall(DMGetCoordinatesLocal(dm, &X));
200: PetscCall(DMGetCoordinatesLocal(clonedm, &cloneX));
201: PetscCall(VecGetLocalSize(X, &num_x));
202: PetscCall(VecGetLocalSize(cloneX, &clone_num_x));
203: 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);
205: PetscCall(VecGetArrayRead(X, &x));
206: PetscCall(VecGetArrayRead(cloneX, &clonex));
208: for (PetscInt i = 0; i < num_x; i++) {
209: 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]));
210: }
212: PetscCall(VecRestoreArrayRead(X, &x));
213: PetscCall(VecRestoreArrayRead(cloneX, &clonex));
214: PetscCall(DMDestroy(&clonedm));
215: }
216: }
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: int main(int argc, char **argv)
221: {
222: DM dm;
223: AppCtx user;
224: PetscInt depth;
226: PetscFunctionBeginUser;
227: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
228: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
229: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
230: PetscCall(SetupDiscretization(dm, &user));
231: PetscCall(CheckOffsets(dm, &user, NULL, 0, 0));
232: PetscCall(DMPlexGetDepth(dm, &depth));
233: if (depth > 1) PetscCall(CheckOffsets(dm, &user, "Face Sets", user.check_face, 1));
234: PetscCall(DMDestroy(&dm));
235: PetscCall(PetscFinalize());
236: return 0;
237: }
239: /*TEST
241: test:
242: suffix: 0
243: requires: triangle
244: args: -dm_refine 1 -petscspace_degree 1 -dm_view -offsets_view
246: test:
247: suffix: 1
248: args: -dm_plex_simplex 0 -dm_plex_box_bd periodic,none -dm_plex_box_faces 3,3 -dm_sparse_localize 0 -petscspace_degree 1 \
249: -dm_view -offsets_view
251: test:
252: suffix: cg_2d
253: args: -dm_plex_simplex 0 -dm_plex_box_bd none,none -dm_plex_box_faces 3,3 -petscspace_degree 1 \
254: -dm_view -offsets_view
256: test:
257: suffix: 1d_sfc
258: args: -dm_plex_simplex 0 -dm_plex_dim 1 -dm_plex_shape zbox -dm_plex_box_faces 3 1 -dm_view -coord_ltog_view
260: test:
261: suffix: 2d_sfc
262: nsize: 2
263: 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
265: test:
266: suffix: 2d_sfc_periodic
267: nsize: 2
268: 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
270: testset:
271: 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
272: nsize: 2
273: test:
274: suffix: 2d_sfc_periodic_stranded
275: args: -dm_distribute 0 -dm_plex_box_bd none,periodic
276: test:
277: suffix: 2d_sfc_periodic_stranded_dist
278: args: -dm_distribute 1 -petscpartitioner_type simple -dm_plex_box_bd none,periodic
279: test:
280: suffix: 2d_sfc_biperiodic_stranded
281: args: -dm_distribute 0 -dm_plex_box_bd periodic,periodic
282: test:
283: suffix: 2d_sfc_biperiodic_stranded_dist
284: args: -dm_distribute 1 -petscpartitioner_type simple -dm_plex_box_bd periodic,periodic
286: test:
287: suffix: fv_0
288: requires: triangle
289: args: -dm_refine 1 -use_fe 0 -dm_view -offsets_view
291: TEST*/