Actual source code: ex8.c
1: static char help[] = "Element closure restrictions in tensor/lexicographic/spectral-element ordering using DMPlex\n\n";
3: #include <petscdmplex.h>
4: #include <petscds.h>
6: static PetscErrorCode ViewOffsets(DM dm, Vec X)
7: {
8: PetscInt num_elem, elem_size, num_comp, num_dof;
9: PetscInt *elem_restr_offsets;
10: const PetscScalar *x = NULL;
11: const char *name;
13: PetscFunctionBegin;
14: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
15: PetscCall(DMPlexGetLocalOffsets(dm, NULL, 0, 0, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets));
16: PetscCall(PetscPrintf(PETSC_COMM_SELF, "DM %s offsets: num_elem %" PetscInt_FMT ", size %" PetscInt_FMT ", comp %" PetscInt_FMT ", dof %" PetscInt_FMT "\n", name, num_elem, elem_size, num_comp, num_dof));
17: if (X) PetscCall(VecGetArrayRead(X, &x));
18: for (PetscInt c = 0; c < num_elem; c++) {
19: PetscCall(PetscIntView(elem_size, &elem_restr_offsets[c * elem_size], PETSC_VIEWER_STDOUT_SELF));
20: if (x) {
21: for (PetscInt i = 0; i < elem_size; i++) PetscCall(PetscScalarView(num_comp, &x[elem_restr_offsets[c * elem_size + i]], PETSC_VIEWER_STDERR_SELF));
22: }
23: }
24: if (X) PetscCall(VecRestoreArrayRead(X, &x));
25: PetscCall(PetscFree(elem_restr_offsets));
26: PetscFunctionReturn(PETSC_SUCCESS);
27: }
29: int main(int argc, char **argv)
30: {
31: DM dm;
32: PetscSection section;
33: PetscFE fe;
34: PetscInt dim, Nf = 1, c, cStart, cEnd;
35: PetscBool view_coord = PETSC_FALSE, tensor = PETSC_TRUE, project = PETSC_FALSE;
37: PetscFunctionBeginUser;
38: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
39: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Tensor closure restrictions", "DMPLEX");
40: PetscCall(PetscOptionsBool("-closure_tensor", "Apply DMPlexSetClosurePermutationTensor", "ex8.c", tensor, &tensor, NULL));
41: PetscCall(PetscOptionsBool("-project_coordinates", "Call DMSetCoordinateDisc() explicitly", "ex8.c", project, &project, NULL));
42: PetscCall(PetscOptionsBool("-view_coord", "View coordinates of element closures", "ex8.c", view_coord, &view_coord, NULL));
43: PetscCall(PetscOptionsBoundedInt("-num_fields", "The number of fields to use", "ex8.c", Nf, &Nf, NULL, 1));
44: PetscOptionsEnd();
46: PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
47: PetscCall(DMSetType(dm, DMPLEX));
48: PetscCall(DMSetFromOptions(dm));
49: if (project) {
50: PetscFE fe_coords;
51: PetscInt cdim;
52: PetscCall(DMGetCoordinateDim(dm, &cdim));
53: PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, cdim, cdim, PETSC_FALSE, 1, 1, &fe_coords));
54: PetscCall(DMSetCoordinateDisc(dm, fe_coords, PETSC_TRUE));
55: PetscCall(PetscFEDestroy(&fe_coords));
56: }
57: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
58: PetscCall(DMGetDimension(dm, &dim));
60: if (Nf == 1) {
61: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, PETSC_DETERMINE, &fe));
62: PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
63: PetscCall(PetscFEDestroy(&fe));
64: } else {
65: for (PetscInt f = 0; f < Nf; ++f) {
66: char prefix[16];
67: PetscCall(PetscSNPrintf(prefix, 16, "f%" PetscInt_FMT "_", f));
68: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, prefix, PETSC_DETERMINE, &fe));
69: PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
70: PetscCall(PetscFEDestroy(&fe));
71: }
72: }
73: PetscCall(DMCreateDS(dm));
74: if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL));
75: PetscCall(DMGetLocalSection(dm, §ion));
76: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
77: for (c = cStart; c < cEnd; c++) {
78: PetscInt numindices, *indices;
79: PetscCall(DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, &numindices, &indices, NULL, NULL));
80: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT "\n", c - cStart));
81: PetscCall(PetscIntView(numindices, indices, PETSC_VIEWER_STDOUT_SELF));
82: PetscCall(DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &numindices, &indices, NULL, NULL));
83: }
84: {
85: DMLabel domain_label;
86: IS valueIS;
87: const PetscInt *values;
88: PetscInt Nv;
90: PetscCall(DMGetLabel(dm, "Face Sets", &domain_label));
91: PetscCall(DMLabelGetNumValues(domain_label, &Nv));
92: // Check that discretization bas any values on faces
93: {
94: PetscDS ds;
95: PetscFE fe;
96: IS fieldIS;
97: const PetscInt *fields;
98: PetscInt Nf, dmf = 0, dsf = -1;
100: PetscCall(DMGetRegionDS(dm, domain_label, &fieldIS, &ds, NULL));
101: // Translate DM field number to DS field number
102: PetscCall(ISGetIndices(fieldIS, &fields));
103: PetscCall(ISGetSize(fieldIS, &Nf));
104: for (PetscInt f = 0; f < Nf; ++f) {
105: if (dmf == fields[f]) {
106: dsf = f;
107: break;
108: }
109: }
110: PetscCall(ISRestoreIndices(fieldIS, &fields));
111: PetscCall(PetscDSGetDiscretization(ds, dsf, (PetscObject *)&fe));
112: PetscCall(PetscFEGetHeightSubspace(fe, 1, &fe));
113: if (!fe) Nv = 0;
114: }
115: PetscCall(DMLabelGetValueIS(domain_label, &valueIS));
116: PetscCall(ISGetIndices(valueIS, &values));
117: for (PetscInt v = 0; v < Nv; ++v) {
118: PetscInt *elem_restr_offsets_face;
119: PetscInt num_elem, elem_size, num_comp, num_dof;
121: PetscCall(DMPlexGetLocalOffsets(dm, domain_label, values[v], 1, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets_face));
122: PetscCall(PetscPrintf(PETSC_COMM_SELF, "========= Face restriction; marker: %" PetscInt_FMT " ========\n", values[v]));
123: PetscCall(PetscPrintf(PETSC_COMM_SELF, "num_elem: %" PetscInt_FMT ", elem_size: %" PetscInt_FMT ", num_dof: %" PetscInt_FMT "\n", num_elem, elem_size, num_dof));
124: for (PetscInt c = 0; c < num_elem; c++) PetscCall(PetscIntView(elem_size, &elem_restr_offsets_face[c * elem_size], PETSC_VIEWER_STDOUT_SELF));
125: PetscCall(PetscFree(elem_restr_offsets_face));
126: }
127: PetscCall(ISRestoreIndices(valueIS, &values));
128: PetscCall(ISDestroy(&valueIS));
129: }
130: if (view_coord) {
131: DM cdm;
132: Vec X;
133: PetscInt cdim;
135: PetscCall(DMGetCoordinatesLocalSetUp(dm));
136: PetscCall(DMGetCoordinateDim(dm, &cdim));
137: PetscCall(DMGetCoordinateDM(dm, &cdm));
138: PetscCall(PetscObjectSetName((PetscObject)cdm, "coords"));
139: if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL));
140: for (c = cStart; c < cEnd; ++c) {
141: const PetscScalar *array;
142: PetscScalar *x = NULL;
143: PetscInt ndof;
144: PetscBool isDG;
146: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &ndof, &array, &x));
147: PetscCheck(ndof % cdim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "ndof not divisible by cdim");
148: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT " coordinates\n", c - cStart));
149: for (PetscInt i = 0; i < ndof; i += cdim) PetscCall(PetscScalarView(cdim, &x[i], PETSC_VIEWER_STDOUT_SELF));
150: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &ndof, &array, &x));
151: }
152: PetscCall(ViewOffsets(dm, NULL));
153: PetscCall(DMGetCoordinatesLocal(dm, &X));
154: PetscCall(ViewOffsets(cdm, X));
155: }
156: PetscCall(DMDestroy(&dm));
157: PetscCall(PetscFinalize());
158: return 0;
159: }
161: /*TEST
163: test:
164: suffix: 1d_q2
165: args: -dm_plex_dim 1 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2
166: test:
167: suffix: 2d_q1
168: args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2
169: test:
170: suffix: 2d_q1d
171: args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 -petscdualspace_lagrange_continuity 0
172: test:
173: suffix: 2d_q1_q1d
174: args: -dm_plex_dim 2 -num_fields 2 -f0_petscspace_degree 1 -f1_petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -f1_petscdualspace_lagrange_continuity 0
175: test:
176: suffix: 2d_q2
177: args: -dm_plex_dim 2 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2
178: test:
179: suffix: 2d_q2_q1
180: args: -dm_plex_dim 2 -num_fields 2 -f0_petscspace_degree 2 -f1_petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,1
181: test:
182: suffix: 2d_q2_p0
183: args: -dm_plex_dim 2 -num_fields 2 -f0_petscspace_degree 2 -f1_petscspace_degree 0 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -f1_petscdualspace_lagrange_continuity 0
184: test:
185: suffix: 2d_q2_q1d
186: args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -num_fields 2 \
187: -f0_petscspace_degree 2 \
188: -f1_petscspace_degree 1 -f1_petscdualspace_lagrange_continuity 0
189: test:
190: suffix: 2d_q2_p1d
191: args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -num_fields 2 \
192: -f0_petscspace_degree 2 \
193: -f1_petscspace_degree 1 -f1_petscspace_poly_tensor 0 -f1_petscdualspace_lagrange_continuity 0
194: test:
195: suffix: 2d_q3
196: args: -dm_plex_dim 2 -petscspace_degree 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1
198: test:
199: suffix: 3d_q1
200: args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
201: test:
202: suffix: 1d_q1_periodic
203: requires: !complex
204: args: -dm_plex_dim 1 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 3 -dm_plex_box_bd periodic -dm_view -view_coord
205: test:
206: suffix: 2d_q1_periodic
207: requires: !complex
208: args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 3,2 -dm_plex_box_bd periodic,none -dm_view -view_coord
209: test:
210: suffix: 3d_q1_periodic
211: requires: !complex
212: args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 3,2,1 -dm_plex_box_bd periodic,none,none -dm_view -view_coord
213: test:
214: suffix: 3d_q1_periodic_project
215: requires: !complex
216: args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,3 -dm_plex_box_bd none,none,periodic -dm_view -view_coord -project_coordinates
218: test:
219: suffix: 3d_q2_periodic # not actually periodic because only 2 cells
220: args: -dm_plex_dim 3 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2,2 -dm_plex_box_bd periodic,none,periodic -dm_view
222: TEST*/