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