Actual source code: ex10.c
1: static char help[] = "Test for mesh reordering\n\n";
3: #include <petscdmplex.h>
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: char orderType[256]; /* The ordering passed to DMPlexGetOrdering() */
13: } AppCtx;
15: PetscErrorCode ProcessOptions(AppCtx *options)
16: {
17: PetscInt len;
18: PetscBool flg;
20: PetscFunctionBegin;
21: options->numFields = 1;
22: options->numComponents = NULL;
23: options->numDof = NULL;
24: options->numGroups = 0;
25: PetscCall(PetscStrncpy(options->orderType, MATORDERINGRCM, sizeof(options->orderType)));
27: PetscOptionsBegin(PETSC_COMM_SELF, "", "Meshing Problem Options", "DMPLEX");
28: PetscCall(PetscOptionsBoundedInt("-num_fields", "The number of section fields", "ex10.c", options->numFields, &options->numFields, NULL, 1));
29: if (options->numFields) {
30: len = options->numFields;
31: PetscCall(PetscCalloc1(len, &options->numComponents));
32: PetscCall(PetscOptionsIntArray("-num_components", "The number of components per field", "ex10.c", options->numComponents, &len, &flg));
33: 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);
34: }
35: PetscCall(PetscOptionsBoundedInt("-num_groups", "Group permutation by this many label values", "ex10.c", options->numGroups, &options->numGroups, NULL, 0));
36: PetscCall(PetscOptionsString("-order_type", "The ordering type, for example rcm or morton", "ex10.c", options->orderType, options->orderType, sizeof(options->orderType), NULL));
37: PetscOptionsEnd();
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: PetscErrorCode CleanupContext(AppCtx *user)
42: {
43: PetscFunctionBegin;
44: PetscCall(PetscFree(user->numComponents));
45: PetscCall(PetscFree(user->numDof));
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: /* This mesh comes from~\cite{saad2003}, Fig. 2.10, p. 70. */
50: PetscErrorCode CreateTestMesh(MPI_Comm comm, DM *dm, AppCtx *options)
51: {
52: 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};
53: 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};
55: PetscFunctionBegin;
56: PetscCall(DMPlexCreateFromCellListPetsc(comm, 2, 16, 15, 3, PETSC_FALSE, cells, 2, coords, dm));
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
60: PetscErrorCode TestReordering(DM dm, AppCtx *user)
61: {
62: DM pdm;
63: IS perm;
64: Mat A, pA;
65: PetscInt bw, pbw;
66: MatOrderingType order = user->orderType;
68: PetscFunctionBegin;
69: PetscCall(DMPlexGetOrdering(dm, order, NULL, &perm));
70: PetscCall(DMPlexPermute(dm, perm, &pdm));
71: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pdm, "perm_"));
72: PetscCall(DMSetFromOptions(pdm));
73: PetscCall(ISDestroy(&perm));
74: PetscCall(DMViewFromOptions(dm, NULL, "-orig_dm_view"));
75: PetscCall(DMViewFromOptions(pdm, NULL, "-dm_view"));
76: PetscCall(DMCreateMatrix(dm, &A));
77: PetscCall(DMCreateMatrix(pdm, &pA));
78: PetscCall(MatComputeBandwidth(A, 0.0, &bw));
79: PetscCall(MatComputeBandwidth(pA, 0.0, &pbw));
80: PetscCall(MatViewFromOptions(A, NULL, "-orig_mat_view"));
81: PetscCall(MatViewFromOptions(pA, NULL, "-perm_mat_view"));
82: PetscCall(MatDestroy(&A));
83: PetscCall(MatDestroy(&pA));
84: PetscCall(DMDestroy(&pdm));
85: if (pbw > bw) {
86: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Ordering method %s increased bandwidth from %" PetscInt_FMT " to %" PetscInt_FMT "\n", order, bw, pbw));
87: } else {
88: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Ordering method %s reduced bandwidth from %" PetscInt_FMT " to %" PetscInt_FMT "\n", order, bw, pbw));
89: }
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: PetscErrorCode CreateGroupLabel(DM dm, PetscInt numGroups, DMLabel *label, AppCtx *options)
94: {
95: const PetscInt groupA[10] = {15, 3, 13, 12, 2, 10, 7, 6, 0, 4};
96: const PetscInt groupB[6] = {14, 11, 9, 1, 8, 5};
98: PetscFunctionBegin;
99: if (numGroups < 2) {
100: *label = NULL;
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
103: PetscCheck(numGroups == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Test only coded for 2 groups, not %" PetscInt_FMT, numGroups);
104: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "groups", label));
105: for (PetscInt c = 0; c < 10; ++c) PetscCall(DMLabelSetValue(*label, groupA[c], 101));
106: for (PetscInt c = 0; c < 6; ++c) PetscCall(DMLabelSetValue(*label, groupB[c], 1001));
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: PetscErrorCode TestReorderingByGroup(DM dm, AppCtx *user)
111: {
112: DM pdm;
113: DMLabel label;
114: Mat A, pA;
115: MatOrderingType order = user->orderType;
116: IS perm;
118: PetscFunctionBegin;
119: PetscCall(CreateGroupLabel(dm, user->numGroups, &label, user));
120: PetscCall(DMPlexGetOrdering(dm, order, label, &perm));
121: PetscCall(DMLabelDestroy(&label));
122: PetscCall(DMPlexPermute(dm, perm, &pdm));
123: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pdm, "perm_"));
124: PetscCall(DMSetFromOptions(pdm));
125: PetscCall(DMViewFromOptions(dm, NULL, "-orig_dm_view"));
126: PetscCall(DMViewFromOptions(pdm, NULL, "-perm_dm_view"));
127: PetscCall(ISDestroy(&perm));
128: PetscCall(DMCreateMatrix(dm, &A));
129: PetscCall(DMCreateMatrix(pdm, &pA));
130: PetscCall(MatViewFromOptions(A, NULL, "-orig_mat_view"));
131: PetscCall(MatViewFromOptions(pA, NULL, "-perm_mat_view"));
132: PetscCall(MatDestroy(&A));
133: PetscCall(MatDestroy(&pA));
134: PetscCall(DMDestroy(&pdm));
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: int main(int argc, char **argv)
139: {
140: DM dm;
141: PetscSection s;
142: AppCtx user;
143: PetscInt dim;
145: PetscFunctionBeginUser;
146: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
147: PetscCall(ProcessOptions(&user));
148: if (user.numGroups < 1) {
149: PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
150: PetscCall(DMSetType(dm, DMPLEX));
151: } else {
152: PetscCall(CreateTestMesh(PETSC_COMM_WORLD, &dm, &user));
153: }
154: PetscCall(DMSetFromOptions(dm));
155: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
156: PetscCall(DMGetDimension(dm, &dim));
157: {
158: PetscInt len = (dim + 1) * PetscMax(1, user.numFields);
159: PetscBool flg;
161: PetscCall(PetscCalloc1(len, &user.numDof));
162: PetscOptionsBegin(PETSC_COMM_SELF, "", "Meshing Problem Options", "DMPLEX");
163: PetscCall(PetscOptionsIntArray("-num_dof", "The dof signature for the section", "ex10.c", user.numDof, &len, &flg));
164: 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));
165: PetscOptionsEnd();
166: }
167: if (user.numGroups < 1) {
168: PetscCall(DMSetNumFields(dm, user.numFields));
169: PetscCall(DMCreateDS(dm));
170: PetscCall(DMPlexCreateSection(dm, NULL, user.numComponents, user.numDof, 0, NULL, NULL, NULL, NULL, &s));
171: PetscCall(DMSetLocalSection(dm, s));
172: PetscCall(PetscSectionDestroy(&s));
173: PetscCall(TestReordering(dm, &user));
174: } else {
175: PetscCall(DMSetNumFields(dm, user.numFields));
176: PetscCall(DMCreateDS(dm));
177: PetscCall(DMPlexCreateSection(dm, NULL, user.numComponents, user.numDof, 0, NULL, NULL, NULL, NULL, &s));
178: PetscCall(DMSetLocalSection(dm, s));
179: PetscCall(PetscSectionDestroy(&s));
180: PetscCall(TestReorderingByGroup(dm, &user));
181: }
182: PetscCall(DMDestroy(&dm));
183: PetscCall(CleanupContext(&user));
184: PetscCall(PetscFinalize());
185: return 0;
186: }
188: /*TEST
190: # Space-filling-curve ordering. Uses the same meshes as tests 1 and 3, so it exercises the
191: # DMPLEXCURVEMORTON branch of DMPlexGetOrdering() in 2D and 3D.
192: test:
193: suffix: morton_2d
194: args: -dm_plex_simplex 0 -num_dof 1,0,0 -mat_view -dm_coord_space 0 -order_type morton
195: test:
196: suffix: morton_3d
197: args: -dm_plex_dim 3 -dm_plex_simplex 0 -num_dof 1,0,0,0 -mat_view -dm_coord_space 0 -order_type morton
198: test:
199: suffix: morton_refined
200: args: -dm_plex_simplex 0 -dm_refine 2 -num_dof 1,0,0 -order_type morton
201: # A periodic mesh with localized coordinates keeps a second, per-cell coordinate field. The
202: # centroid of a cell that crosses the periodic boundary must come from that field, so this covers
203: # the other branch of DMPlexGetCellCoordinates(). Four of the 16 cells cross the boundary.
204: test:
205: suffix: morton_periodic
206: args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 -dm_plex_box_bd periodic,none -num_dof 1,0,0 -order_type morton
208: # Two cell tests 0-3
209: test:
210: suffix: 0
211: requires: triangle
212: args: -dm_plex_simplex 1 -num_dof 1,0,0 -mat_view -dm_coord_space 0
213: test:
214: suffix: 1
215: args: -dm_plex_simplex 0 -num_dof 1,0,0 -mat_view -dm_coord_space 0
216: test:
217: suffix: 2
218: requires: ctetgen
219: args: -dm_plex_dim 3 -dm_plex_simplex 1 -num_dof 1,0,0,0 -mat_view -dm_coord_space 0
220: test:
221: suffix: 3
222: args: -dm_plex_dim 3 -dm_plex_simplex 0 -num_dof 1,0,0,0 -mat_view -dm_coord_space 0
223: # Refined tests 4-7
224: test:
225: suffix: 4
226: requires: triangle
227: args: -dm_plex_simplex 1 -dm_refine_volume_limit_pre 0.00625 -num_dof 1,0,0
228: test:
229: suffix: 5
230: args: -dm_plex_simplex 0 -dm_refine 1 -num_dof 1,0,0
231: test:
232: suffix: 6
233: requires: ctetgen
234: args: -dm_plex_dim 3 -dm_plex_simplex 1 -dm_refine_volume_limit_pre 0.00625 -num_dof 1,0,0,0
235: test:
236: suffix: 7
237: args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_refine 1 -num_dof 1,0,0,0
238: # Parallel tests
239: # Grouping tests
240: test:
241: suffix: group_1
242: args: -num_groups 1 -num_dof 1,0,0 -is_view -orig_mat_view -perm_mat_view
243: test:
244: suffix: group_2
245: args: -num_groups 2 -num_dof 1,0,0 -is_view -perm_mat_view
247: TEST*/