1 static char help[] = "Element closure restrictions in tensor/lexicographic/spectral-element ordering using DMPlex\n\n";
2
3 #include <petscdmplex.h>
4 #include <petscds.h>
5
ViewOffsets(DM dm,Vec X)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;
12
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 }
28
main(int argc,char ** argv)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;
36
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();
45
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_FALSE, PETSC_TRUE));
55 PetscCall(PetscFEDestroy(&fe_coords));
56 }
57 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
58 PetscCall(DMGetDimension(dm, &dim));
59
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;
89
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;
99
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;
120
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, cell_dm;
132 Vec X;
133 PetscInt cdim;
134 PetscBool sparseLocalize;
135
136 PetscCall(DMGetCoordinatesLocalSetUp(dm));
137 PetscCall(DMGetCoordinateDim(dm, &cdim));
138 PetscCall(DMGetCoordinateDM(dm, &cdm));
139 PetscCall(PetscObjectSetName((PetscObject)cdm, "coords"));
140 if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL));
141 PetscCall(DMLocalizeCoordinates(dm));
142 PetscCall(DMGetSparseLocalize(dm, &sparseLocalize));
143 for (c = cStart; c < cEnd; ++c) {
144 const PetscScalar *array;
145 PetscScalar *x = NULL;
146 PetscInt ndof;
147 PetscBool isDG;
148
149 PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &ndof, &array, &x));
150 PetscCheck(ndof % cdim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "ndof not divisible by cdim");
151 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT " coordinates\n", c - cStart));
152 for (PetscInt i = 0; i < ndof; i += cdim) PetscCall(PetscScalarView(cdim, &x[i], PETSC_VIEWER_STDOUT_SELF));
153 PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &ndof, &array, &x));
154 }
155 PetscCall(ViewOffsets(dm, NULL));
156 PetscCall(DMGetCoordinatesLocal(dm, &X));
157 PetscCall(ViewOffsets(cdm, X));
158 PetscCall(DMGetCellCoordinateDM(dm, &cell_dm));
159 PetscCall(PetscObjectSetName((PetscObject)cell_dm, "cell coords"));
160 if (cell_dm && !sparseLocalize) {
161 PetscCall(DMGetCellCoordinatesLocal(dm, &X));
162 PetscCall(ViewOffsets(cell_dm, X));
163 }
164 }
165 PetscCall(DMDestroy(&dm));
166 PetscCall(PetscFinalize());
167 return 0;
168 }
169
170 /*TEST
171
172 test:
173 suffix: 1d_q2
174 args: -dm_plex_dim 1 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2
175 test:
176 suffix: 2d_q1
177 args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2
178 test:
179 suffix: 2d_q1d
180 args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 -petscdualspace_lagrange_continuity 0
181 test:
182 suffix: 2d_q1_q1d
183 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
184 test:
185 suffix: 2d_q2
186 args: -dm_plex_dim 2 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2
187 test:
188 suffix: 2d_q2_q1
189 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
190 test:
191 suffix: 2d_q2_p0
192 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
193 test:
194 suffix: 2d_q2_q1d
195 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -num_fields 2 \
196 -f0_petscspace_degree 2 \
197 -f1_petscspace_degree 1 -f1_petscdualspace_lagrange_continuity 0
198 test:
199 suffix: 2d_q2_p1d
200 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -num_fields 2 \
201 -f0_petscspace_degree 2 \
202 -f1_petscspace_degree 1 -f1_petscspace_poly_tensor 0 -f1_petscdualspace_lagrange_continuity 0
203 test:
204 suffix: 2d_q3
205 args: -dm_plex_dim 2 -petscspace_degree 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1
206
207 test:
208 suffix: 3d_q1
209 args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
210 test:
211 suffix: 1d_q1_periodic
212 requires: !complex
213 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
214 testset:
215 suffix: 2d_q1_periodic
216 requires: !complex
217 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 3,2 -dm_plex_box_bd periodic,none \
218 -petscspace_degree 1 -dm_view -view_coord
219
220 test:
221 test:
222 suffix: sparse
223 args: -dm_sparse_localize false -dm_localize 0
224
225 test:
226 suffix: 3d_q1_periodic
227 requires: !complex
228 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
229 test:
230 suffix: 3d_q1_periodic_project
231 requires: !complex
232 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
233
234 test:
235 suffix: 3d_q2_periodic # not actually periodic because only 2 cells
236 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
237
238 TEST*/
239