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;
47: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, -1, &fe));
48: PetscCall(PetscObjectSetName((PetscObject)fe, "scalar"));
49: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
50: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe));
51: PetscCall(PetscFEDestroy(&fe));
52: } else {
53: PetscFV fv;
55: PetscCall(PetscFVCreate(PETSC_COMM_SELF, &fv));
56: PetscCall(PetscFVSetType(fv, PETSCFVLEASTSQUARES));
57: PetscCall(PetscFVSetNumComponents(fv, dim));
58: PetscCall(PetscFVSetSpatialDimension(fv, dim));
59: PetscCall(PetscFVSetFromOptions(fv));
60: PetscCall(PetscFVSetUp(fv));
61: PetscCall(PetscObjectSetName((PetscObject)fv, "vector"));
62: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fv));
63: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fv));
64: PetscCall(PetscFVDestroy(&fv));
65: }
66: PetscCall(DMCreateDS(dm));
67: while (cdm) {
68: PetscCall(DMCopyDisc(dm, cdm));
69: PetscCall(DMGetCoarseDM(cdm, &cdm));
70: }
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: static PetscErrorCode CheckOffsets(DM dm, AppCtx *user, const char *domain_name, PetscInt label_value, PetscInt height)
75: {
76: const char *height_name[] = {"cells", "faces"};
77: DMLabel domain_label = NULL;
78: DM cdm;
79: PetscInt Nf, f;
80: ISLocalToGlobalMapping ltog;
82: PetscFunctionBeginUser;
83: if (domain_name) PetscCall(DMGetLabel(dm, domain_name, &domain_label));
84: 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)" : ""));
85: if (domain_name && !domain_label) PetscFunctionReturn(PETSC_SUCCESS);
86: if (user->closure_tensor) PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL));
87: // Offsets for cell closures
88: PetscCall(DMGetNumFields(dm, &Nf));
89: for (f = 0; f < Nf; ++f) {
90: PetscObject obj;
91: PetscClassId id;
92: char name[PETSC_MAX_PATH_LEN];
94: PetscCall(DMGetField(dm, f, NULL, &obj));
95: PetscCall(PetscObjectGetClassId(obj, &id));
96: if (id == PETSCFE_CLASSID) {
97: IS offIS;
98: PetscInt *offsets, Ncell, Ncl, Nc, n;
100: PetscCall(DMPlexGetLocalOffsets(dm, domain_label, label_value, height, f, &Ncell, &Ncl, &Nc, &n, &offsets));
101: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Ncell * Ncl, offsets, PETSC_OWN_POINTER, &offIS));
102: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Field %" PetscInt_FMT " Offsets", f));
103: PetscCall(PetscObjectSetName((PetscObject)offIS, name));
104: PetscCall(ISViewFromOptions(offIS, NULL, "-offsets_view"));
105: PetscCall(ISDestroy(&offIS));
106: } else if (id == PETSCFV_CLASSID) {
107: IS offIS;
108: PetscInt *offsets, *offsetsNeg, *offsetsPos, Nface, Nc, n;
110: PetscCall(DMPlexGetLocalOffsetsSupport(dm, domain_label, label_value, &Nface, &Nc, &n, &offsetsNeg, &offsetsPos));
111: PetscCall(PetscMalloc1(Nface * Nc * 2, &offsets));
112: for (PetscInt f = 0, i = 0; f < Nface; ++f) {
113: for (PetscInt c = 0; c < Nc; ++c) offsets[i++] = offsetsNeg[f * Nc + c];
114: for (PetscInt c = 0; c < Nc; ++c) offsets[i++] = offsetsPos[f * Nc + c];
115: }
116: PetscCall(PetscFree(offsetsNeg));
117: PetscCall(PetscFree(offsetsPos));
118: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Nface * Nc * 2, offsets, PETSC_OWN_POINTER, &offIS));
119: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Field %" PetscInt_FMT " Offsets", f));
120: PetscCall(PetscObjectSetName((PetscObject)offIS, name));
121: PetscCall(ISViewFromOptions(offIS, NULL, "-offsets_view"));
122: PetscCall(ISDestroy(&offIS));
123: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Unrecognized type for DM field %" PetscInt_FMT, f);
124: }
125: PetscCall(DMGetLocalToGlobalMapping(dm, <og));
126: PetscCall(ISLocalToGlobalMappingViewFromOptions(ltog, NULL, "-ltog_view"));
128: // Offsets for coordinates
129: {
130: Vec X;
131: PetscSection s;
132: const PetscScalar *x;
133: const char *cname;
134: PetscInt cdim, *offsets, Ncell, Ncl, Nc, n;
135: PetscBool isDG = PETSC_FALSE;
137: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
138: if (!cdm) {
139: PetscCall(DMGetCoordinateDM(dm, &cdm));
140: cname = "Coordinates";
141: PetscCall(DMGetCoordinatesLocal(dm, &X));
142: } else {
143: isDG = PETSC_TRUE;
144: cname = "DG Coordinates";
145: PetscCall(DMGetCellCoordinatesLocal(dm, &X));
146: }
147: if (isDG && height) PetscFunctionReturn(PETSC_SUCCESS);
148: if (domain_name) PetscCall(DMGetLabel(cdm, domain_name, &domain_label));
149: if (user->closure_tensor) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL));
150: PetscCall(DMPlexGetLocalOffsets(cdm, domain_label, label_value, height, 0, &Ncell, &Ncl, &Nc, &n, &offsets));
151: PetscCall(DMGetCoordinateDim(dm, &cdim));
152: PetscCheck(Nc == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Geometric dimension %" PetscInt_FMT " should be %" PetscInt_FMT, Nc, cdim);
153: PetscCall(DMGetLocalSection(cdm, &s));
154: PetscCall(VecGetArrayRead(X, &x));
155: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s by element in %s order\n", cname, user->closure_tensor ? "tensor" : "bfs"));
156: for (PetscInt c = 0; c < Ncell; ++c) {
157: for (PetscInt v = 0; v < Ncl; ++v) {
158: PetscInt off = offsets[c * Ncl + v], dgdof;
159: const PetscScalar *vx = &x[off];
161: if (isDG) {
162: PetscCall(PetscSectionGetDof(s, c, &dgdof));
163: PetscCheck(Ncl * Nc == dgdof, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Offset size %" PetscInt_FMT " should be %" PetscInt_FMT, Ncl * Nc, dgdof);
164: }
165: switch (cdim) {
166: case 1:
167: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT "] %" PetscInt_FMT " <-- %2" PetscInt_FMT " (% 4.2f)\n", c, v, off, (double)PetscRealPart(vx[0])));
168: break;
169: case 2:
170: 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])));
171: break;
172: case 3:
173: 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])));
174: }
175: }
176: }
177: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout));
178: PetscCall(VecRestoreArrayRead(X, &x));
179: PetscCall(PetscFree(offsets));
180: PetscCall(DMGetLocalToGlobalMapping(cdm, <og));
181: PetscCall(ISLocalToGlobalMappingViewFromOptions(ltog, NULL, "-coord_ltog_view"));
182: }
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: int main(int argc, char **argv)
187: {
188: DM dm;
189: AppCtx user;
190: PetscInt depth;
192: PetscFunctionBeginUser;
193: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
194: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
195: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
196: PetscCall(SetupDiscretization(dm, &user));
197: PetscCall(CheckOffsets(dm, &user, NULL, 0, 0));
198: PetscCall(DMPlexGetDepth(dm, &depth));
199: if (depth > 1) PetscCall(CheckOffsets(dm, &user, "Face Sets", user.check_face, 1));
200: PetscCall(DMDestroy(&dm));
201: PetscCall(PetscFinalize());
202: return 0;
203: }
205: /*TEST
207: test:
208: suffix: 0
209: requires: triangle
210: args: -dm_refine 1 -petscspace_degree 1 -dm_view -offsets_view
212: test:
213: suffix: 1
214: args: -dm_plex_simplex 0 -dm_plex_box_bd periodic,none -dm_plex_box_faces 3,3 -dm_sparse_localize 0 -petscspace_degree 1 \
215: -dm_view -offsets_view
217: test:
218: suffix: cg_2d
219: args: -dm_plex_simplex 0 -dm_plex_box_bd none,none -dm_plex_box_faces 3,3 -petscspace_degree 1 \
220: -dm_view -offsets_view
222: test:
223: suffix: 1d_sfc
224: args: -dm_plex_simplex 0 -dm_plex_dim 1 -dm_plex_shape zbox -dm_plex_box_faces 3 1 -dm_view -coord_ltog_view
226: test:
227: suffix: 2d_sfc
228: nsize: 2
229: 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
231: test:
232: suffix: 2d_sfc_periodic
233: nsize: 2
234: 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
236: testset:
237: args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_shape zbox -dm_plex_box_faces 3,2 -petscspace_degree 1 -dm_plex_box_bd none,periodic -dm_view ::ascii_info_detail -closure_tensor
238: nsize: 2
239: test:
240: suffix: 2d_sfc_periodic_stranded
241: args: -dm_distribute 0
242: test:
243: suffix: 2d_sfc_periodic_stranded_dist
244: args: -dm_distribute 1 -petscpartitioner_type simple
246: test:
247: suffix: fv_0
248: requires: triangle
249: args: -dm_refine 1 -use_fe 0 -dm_view -offsets_view
251: TEST*/