1 static char help[] = "Test for mesh reordering\n\n";
2
3 #include <petscdmplex.h>
4
5 typedef struct {
6 PetscInt dim; /* The topological mesh dimension */
7 PetscReal refinementLimit; /* Maximum volume of a refined cell */
8 PetscInt numFields; /* The number of section fields */
9 PetscInt *numComponents; /* The number of field components */
10 PetscInt *numDof; /* The dof signature for the section */
11 PetscInt numGroups; /* If greater than 1, use grouping in test */
12 } AppCtx;
13
ProcessOptions(AppCtx * options)14 PetscErrorCode ProcessOptions(AppCtx *options)
15 {
16 PetscInt len;
17 PetscBool flg;
18
19 PetscFunctionBegin;
20 options->numFields = 1;
21 options->numComponents = NULL;
22 options->numDof = NULL;
23 options->numGroups = 0;
24
25 PetscOptionsBegin(PETSC_COMM_SELF, "", "Meshing Problem Options", "DMPLEX");
26 PetscCall(PetscOptionsBoundedInt("-num_fields", "The number of section fields", "ex10.c", options->numFields, &options->numFields, NULL, 1));
27 if (options->numFields) {
28 len = options->numFields;
29 PetscCall(PetscCalloc1(len, &options->numComponents));
30 PetscCall(PetscOptionsIntArray("-num_components", "The number of components per field", "ex10.c", options->numComponents, &len, &flg));
31 PetscCheck(!flg || !(len != options->numFields), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Length of components array is %" PetscInt_FMT " should be %" PetscInt_FMT, len, options->numFields);
32 }
33 PetscCall(PetscOptionsBoundedInt("-num_groups", "Group permutation by this many label values", "ex10.c", options->numGroups, &options->numGroups, NULL, 0));
34 PetscOptionsEnd();
35 PetscFunctionReturn(PETSC_SUCCESS);
36 }
37
CleanupContext(AppCtx * user)38 PetscErrorCode CleanupContext(AppCtx *user)
39 {
40 PetscFunctionBegin;
41 PetscCall(PetscFree(user->numComponents));
42 PetscCall(PetscFree(user->numDof));
43 PetscFunctionReturn(PETSC_SUCCESS);
44 }
45
46 /* This mesh comes from~\cite{saad2003}, Fig. 2.10, p. 70. */
CreateTestMesh(MPI_Comm comm,DM * dm,AppCtx * options)47 PetscErrorCode CreateTestMesh(MPI_Comm comm, DM *dm, AppCtx *options)
48 {
49 const PetscInt cells[16 * 3] = {6, 7, 8, 7, 9, 10, 10, 11, 12, 11, 13, 14, 0, 6, 8, 6, 2, 7, 1, 8, 7, 1, 7, 10, 2, 9, 7, 10, 9, 4, 1, 10, 12, 10, 4, 11, 12, 11, 3, 3, 11, 14, 11, 4, 13, 14, 13, 5};
50 const PetscReal coords[15 * 2] = {0, -3, 0, -1, 2, -1, 0, 1, 2, 1, 0, 3, 1, -2, 1, -1, 0, -2, 2, 0, 1, 0, 1, 1, 0, 0, 1, 2, 0, 2};
51
52 PetscFunctionBegin;
53 PetscCall(DMPlexCreateFromCellListPetsc(comm, 2, 16, 15, 3, PETSC_FALSE, cells, 2, coords, dm));
54 PetscFunctionReturn(PETSC_SUCCESS);
55 }
56
TestReordering(DM dm,AppCtx * user)57 PetscErrorCode TestReordering(DM dm, AppCtx *user)
58 {
59 DM pdm;
60 IS perm;
61 Mat A, pA;
62 PetscInt bw, pbw;
63 MatOrderingType order = MATORDERINGRCM;
64
65 PetscFunctionBegin;
66 PetscCall(DMPlexGetOrdering(dm, order, NULL, &perm));
67 PetscCall(DMPlexPermute(dm, perm, &pdm));
68 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pdm, "perm_"));
69 PetscCall(DMSetFromOptions(pdm));
70 PetscCall(ISDestroy(&perm));
71 PetscCall(DMViewFromOptions(dm, NULL, "-orig_dm_view"));
72 PetscCall(DMViewFromOptions(pdm, NULL, "-dm_view"));
73 PetscCall(DMCreateMatrix(dm, &A));
74 PetscCall(DMCreateMatrix(pdm, &pA));
75 PetscCall(MatComputeBandwidth(A, 0.0, &bw));
76 PetscCall(MatComputeBandwidth(pA, 0.0, &pbw));
77 PetscCall(MatViewFromOptions(A, NULL, "-orig_mat_view"));
78 PetscCall(MatViewFromOptions(pA, NULL, "-perm_mat_view"));
79 PetscCall(MatDestroy(&A));
80 PetscCall(MatDestroy(&pA));
81 PetscCall(DMDestroy(&pdm));
82 if (pbw > bw) {
83 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Ordering method %s increased bandwidth from %" PetscInt_FMT " to %" PetscInt_FMT "\n", order, bw, pbw));
84 } else {
85 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Ordering method %s reduced bandwidth from %" PetscInt_FMT " to %" PetscInt_FMT "\n", order, bw, pbw));
86 }
87 PetscFunctionReturn(PETSC_SUCCESS);
88 }
89
CreateGroupLabel(DM dm,PetscInt numGroups,DMLabel * label,AppCtx * options)90 PetscErrorCode CreateGroupLabel(DM dm, PetscInt numGroups, DMLabel *label, AppCtx *options)
91 {
92 const PetscInt groupA[10] = {15, 3, 13, 12, 2, 10, 7, 6, 0, 4};
93 const PetscInt groupB[6] = {14, 11, 9, 1, 8, 5};
94 PetscInt c;
95
96 PetscFunctionBegin;
97 if (numGroups < 2) {
98 *label = NULL;
99 PetscFunctionReturn(PETSC_SUCCESS);
100 }
101 PetscCheck(numGroups == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Test only coded for 2 groups, not %" PetscInt_FMT, numGroups);
102 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "groups", label));
103 for (c = 0; c < 10; ++c) PetscCall(DMLabelSetValue(*label, groupA[c], 101));
104 for (c = 0; c < 6; ++c) PetscCall(DMLabelSetValue(*label, groupB[c], 1001));
105 PetscFunctionReturn(PETSC_SUCCESS);
106 }
107
TestReorderingByGroup(DM dm,AppCtx * user)108 PetscErrorCode TestReorderingByGroup(DM dm, AppCtx *user)
109 {
110 DM pdm;
111 DMLabel label;
112 Mat A, pA;
113 MatOrderingType order = MATORDERINGRCM;
114 IS perm;
115
116 PetscFunctionBegin;
117 PetscCall(CreateGroupLabel(dm, user->numGroups, &label, user));
118 PetscCall(DMPlexGetOrdering(dm, order, label, &perm));
119 PetscCall(DMLabelDestroy(&label));
120 PetscCall(DMPlexPermute(dm, perm, &pdm));
121 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pdm, "perm_"));
122 PetscCall(DMSetFromOptions(pdm));
123 PetscCall(DMViewFromOptions(dm, NULL, "-orig_dm_view"));
124 PetscCall(DMViewFromOptions(pdm, NULL, "-perm_dm_view"));
125 PetscCall(ISDestroy(&perm));
126 PetscCall(DMCreateMatrix(dm, &A));
127 PetscCall(DMCreateMatrix(pdm, &pA));
128 PetscCall(MatViewFromOptions(A, NULL, "-orig_mat_view"));
129 PetscCall(MatViewFromOptions(pA, NULL, "-perm_mat_view"));
130 PetscCall(MatDestroy(&A));
131 PetscCall(MatDestroy(&pA));
132 PetscCall(DMDestroy(&pdm));
133 PetscFunctionReturn(PETSC_SUCCESS);
134 }
135
main(int argc,char ** argv)136 int main(int argc, char **argv)
137 {
138 DM dm;
139 PetscSection s;
140 AppCtx user;
141 PetscInt dim;
142
143 PetscFunctionBeginUser;
144 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
145 PetscCall(ProcessOptions(&user));
146 if (user.numGroups < 1) {
147 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
148 PetscCall(DMSetType(dm, DMPLEX));
149 } else {
150 PetscCall(CreateTestMesh(PETSC_COMM_WORLD, &dm, &user));
151 }
152 PetscCall(DMSetFromOptions(dm));
153 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
154 PetscCall(DMGetDimension(dm, &dim));
155 {
156 PetscInt len = (dim + 1) * PetscMax(1, user.numFields);
157 PetscBool flg;
158
159 PetscCall(PetscCalloc1(len, &user.numDof));
160 PetscOptionsBegin(PETSC_COMM_SELF, "", "Meshing Problem Options", "DMPLEX");
161 PetscCall(PetscOptionsIntArray("-num_dof", "The dof signature for the section", "ex10.c", user.numDof, &len, &flg));
162 if (flg) PetscCheck(len == ((dim + 1) * PetscMax(1, user.numFields)), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Length of dof array is %" PetscInt_FMT " should be %" PetscInt_FMT, len, (dim + 1) * PetscMax(1, user.numFields));
163 PetscOptionsEnd();
164 }
165 if (user.numGroups < 1) {
166 PetscCall(DMSetNumFields(dm, user.numFields));
167 PetscCall(DMCreateDS(dm));
168 PetscCall(DMPlexCreateSection(dm, NULL, user.numComponents, user.numDof, 0, NULL, NULL, NULL, NULL, &s));
169 PetscCall(DMSetLocalSection(dm, s));
170 PetscCall(PetscSectionDestroy(&s));
171 PetscCall(TestReordering(dm, &user));
172 } else {
173 PetscCall(DMSetNumFields(dm, user.numFields));
174 PetscCall(DMCreateDS(dm));
175 PetscCall(DMPlexCreateSection(dm, NULL, user.numComponents, user.numDof, 0, NULL, NULL, NULL, NULL, &s));
176 PetscCall(DMSetLocalSection(dm, s));
177 PetscCall(PetscSectionDestroy(&s));
178 PetscCall(TestReorderingByGroup(dm, &user));
179 }
180 PetscCall(DMDestroy(&dm));
181 PetscCall(CleanupContext(&user));
182 PetscCall(PetscFinalize());
183 return 0;
184 }
185
186 /*TEST
187
188 # Two cell tests 0-3
189 test:
190 suffix: 0
191 requires: triangle
192 args: -dm_plex_simplex 1 -num_dof 1,0,0 -mat_view -dm_coord_space 0
193 test:
194 suffix: 1
195 args: -dm_plex_simplex 0 -num_dof 1,0,0 -mat_view -dm_coord_space 0
196 test:
197 suffix: 2
198 requires: ctetgen
199 args: -dm_plex_dim 3 -dm_plex_simplex 1 -num_dof 1,0,0,0 -mat_view -dm_coord_space 0
200 test:
201 suffix: 3
202 args: -dm_plex_dim 3 -dm_plex_simplex 0 -num_dof 1,0,0,0 -mat_view -dm_coord_space 0
203 # Refined tests 4-7
204 test:
205 suffix: 4
206 requires: triangle
207 args: -dm_plex_simplex 1 -dm_refine_volume_limit_pre 0.00625 -num_dof 1,0,0
208 test:
209 suffix: 5
210 args: -dm_plex_simplex 0 -dm_refine 1 -num_dof 1,0,0
211 test:
212 suffix: 6
213 requires: ctetgen
214 args: -dm_plex_dim 3 -dm_plex_simplex 1 -dm_refine_volume_limit_pre 0.00625 -num_dof 1,0,0,0
215 test:
216 suffix: 7
217 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_refine 1 -num_dof 1,0,0,0
218 # Parallel tests
219 # Grouping tests
220 test:
221 suffix: group_1
222 args: -num_groups 1 -num_dof 1,0,0 -is_view -orig_mat_view -perm_mat_view
223 test:
224 suffix: group_2
225 args: -num_groups 2 -num_dof 1,0,0 -is_view -perm_mat_view
226
227 TEST*/
228