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: } AppCtx;

 14: PetscErrorCode ProcessOptions(AppCtx *options)
 15: {
 16:   PetscInt  len;
 17:   PetscBool flg;

 19:   options->numFields     = 1;
 20:   options->numComponents = NULL;
 21:   options->numDof        = NULL;
 22:   options->numGroups     = 0;

 24:   PetscOptionsBegin(PETSC_COMM_SELF, "", "Meshing Problem Options", "DMPLEX");
 25:   PetscOptionsBoundedInt("-num_fields", "The number of section fields", "ex10.c", options->numFields, &options->numFields, NULL, 1);
 26:   if (options->numFields) {
 27:     len = options->numFields;
 28:     PetscCalloc1(len, &options->numComponents);
 29:     PetscOptionsIntArray("-num_components", "The number of components per field", "ex10.c", options->numComponents, &len, &flg);
 31:   }
 32:   PetscOptionsBoundedInt("-num_groups", "Group permutation by this many label values", "ex10.c", options->numGroups, &options->numGroups, NULL, 0);
 33:   PetscOptionsEnd();
 34:   return 0;
 35: }

 37: PetscErrorCode CleanupContext(AppCtx *user)
 38: {
 39:   PetscFree(user->numComponents);
 40:   PetscFree(user->numDof);
 41:   return 0;
 42: }

 44: /* This mesh comes from~\cite{saad2003}, Fig. 2.10, p. 70. */
 45: PetscErrorCode CreateTestMesh(MPI_Comm comm, DM *dm, AppCtx *options)
 46: {
 47:   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};
 48:   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};

 50:   DMPlexCreateFromCellListPetsc(comm, 2, 16, 15, 3, PETSC_FALSE, cells, 2, coords, dm);
 51:   return 0;
 52: }

 54: PetscErrorCode TestReordering(DM dm, AppCtx *user)
 55: {
 56:   DM              pdm;
 57:   IS              perm;
 58:   Mat             A, pA;
 59:   PetscInt        bw, pbw;
 60:   MatOrderingType order = MATORDERINGRCM;

 62:   DMPlexGetOrdering(dm, order, NULL, &perm);
 63:   DMPlexPermute(dm, perm, &pdm);
 64:   PetscObjectSetOptionsPrefix((PetscObject)pdm, "perm_");
 65:   DMSetFromOptions(pdm);
 66:   ISDestroy(&perm);
 67:   DMViewFromOptions(dm, NULL, "-orig_dm_view");
 68:   DMViewFromOptions(pdm, NULL, "-dm_view");
 69:   DMCreateMatrix(dm, &A);
 70:   DMCreateMatrix(pdm, &pA);
 71:   MatComputeBandwidth(A, 0.0, &bw);
 72:   MatComputeBandwidth(pA, 0.0, &pbw);
 73:   MatViewFromOptions(A, NULL, "-orig_mat_view");
 74:   MatViewFromOptions(pA, NULL, "-perm_mat_view");
 75:   MatDestroy(&A);
 76:   MatDestroy(&pA);
 77:   DMDestroy(&pdm);
 78:   if (pbw > bw) {
 79:     PetscPrintf(PetscObjectComm((PetscObject)dm), "Ordering method %s increased bandwidth from %" PetscInt_FMT " to %" PetscInt_FMT "\n", order, bw, pbw);
 80:   } else {
 81:     PetscPrintf(PetscObjectComm((PetscObject)dm), "Ordering method %s reduced bandwidth from %" PetscInt_FMT " to %" PetscInt_FMT "\n", order, bw, pbw);
 82:   }
 83:   return 0;
 84: }

 86: PetscErrorCode CreateGroupLabel(DM dm, PetscInt numGroups, DMLabel *label, AppCtx *options)
 87: {
 88:   const PetscInt groupA[10] = {15, 3, 13, 12, 2, 10, 7, 6, 0, 4};
 89:   const PetscInt groupB[6]  = {14, 11, 9, 1, 8, 5};
 90:   PetscInt       c;

 92:   if (numGroups < 2) {
 93:     *label = NULL;
 94:     return 0;
 95:   }
 97:   DMLabelCreate(PETSC_COMM_SELF, "groups", label);
 98:   for (c = 0; c < 10; ++c) DMLabelSetValue(*label, groupA[c], 101);
 99:   for (c = 0; c < 6; ++c) DMLabelSetValue(*label, groupB[c], 1001);
100:   return 0;
101: }

103: PetscErrorCode TestReorderingByGroup(DM dm, AppCtx *user)
104: {
105:   DM              pdm;
106:   DMLabel         label;
107:   Mat             A, pA;
108:   MatOrderingType order = MATORDERINGRCM;
109:   IS              perm;

111:   CreateGroupLabel(dm, user->numGroups, &label, user);
112:   DMPlexGetOrdering(dm, order, label, &perm);
113:   DMLabelDestroy(&label);
114:   DMPlexPermute(dm, perm, &pdm);
115:   PetscObjectSetOptionsPrefix((PetscObject)pdm, "perm_");
116:   DMSetFromOptions(pdm);
117:   DMViewFromOptions(dm, NULL, "-orig_dm_view");
118:   DMViewFromOptions(pdm, NULL, "-perm_dm_view");
119:   ISDestroy(&perm);
120:   DMCreateMatrix(dm, &A);
121:   DMCreateMatrix(pdm, &pA);
122:   MatViewFromOptions(A, NULL, "-orig_mat_view");
123:   MatViewFromOptions(pA, NULL, "-perm_mat_view");
124:   MatDestroy(&A);
125:   MatDestroy(&pA);
126:   DMDestroy(&pdm);
127:   return 0;
128: }

130: int main(int argc, char **argv)
131: {
132:   DM           dm;
133:   PetscSection s;
134:   AppCtx       user;
135:   PetscInt     dim;

138:   PetscInitialize(&argc, &argv, NULL, help);
139:   ProcessOptions(&user);
140:   if (user.numGroups < 1) {
141:     DMCreate(PETSC_COMM_WORLD, &dm);
142:     DMSetType(dm, DMPLEX);
143:   } else {
144:     CreateTestMesh(PETSC_COMM_WORLD, &dm, &user);
145:   }
146:   DMSetFromOptions(dm);
147:   DMViewFromOptions(dm, NULL, "-dm_view");
148:   DMGetDimension(dm, &dim);
149:   {
150:     PetscInt  len = (dim + 1) * PetscMax(1, user.numFields);
151:     PetscBool flg;

153:     PetscCalloc1(len, &user.numDof);
154:     PetscOptionsBegin(PETSC_COMM_SELF, "", "Meshing Problem Options", "DMPLEX");
155:     PetscOptionsIntArray("-num_dof", "The dof signature for the section", "ex10.c", user.numDof, &len, &flg);
157:     PetscOptionsEnd();
158:   }
159:   if (user.numGroups < 1) {
160:     DMSetNumFields(dm, user.numFields);
161:     DMCreateDS(dm);
162:     DMPlexCreateSection(dm, NULL, user.numComponents, user.numDof, 0, NULL, NULL, NULL, NULL, &s);
163:     DMSetLocalSection(dm, s);
164:     PetscSectionDestroy(&s);
165:     TestReordering(dm, &user);
166:   } else {
167:     DMSetNumFields(dm, user.numFields);
168:     DMCreateDS(dm);
169:     DMPlexCreateSection(dm, NULL, user.numComponents, user.numDof, 0, NULL, NULL, NULL, NULL, &s);
170:     DMSetLocalSection(dm, s);
171:     PetscSectionDestroy(&s);
172:     TestReorderingByGroup(dm, &user);
173:   }
174:   DMDestroy(&dm);
175:   CleanupContext(&user);
176:   PetscFinalize();
177:   return 0;
178: }

180: /*TEST

182:   # Two cell tests 0-3
183:   test:
184:     suffix: 0
185:     requires: triangle
186:     args: -dm_plex_simplex 1 -num_dof 1,0,0 -mat_view -dm_coord_space 0
187:   test:
188:     suffix: 1
189:     args: -dm_plex_simplex 0 -num_dof 1,0,0 -mat_view -dm_coord_space 0
190:   test:
191:     suffix: 2
192:     requires: ctetgen
193:     args: -dm_plex_dim 3 -dm_plex_simplex 1 -num_dof 1,0,0,0 -mat_view -dm_coord_space 0
194:   test:
195:     suffix: 3
196:     args: -dm_plex_dim 3 -dm_plex_simplex 0 -num_dof 1,0,0,0 -mat_view -dm_coord_space 0
197:   # Refined tests 4-7
198:   test:
199:     suffix: 4
200:     requires: triangle
201:     args: -dm_plex_simplex 1 -dm_refine_volume_limit_pre 0.00625 -num_dof 1,0,0
202:   test:
203:     suffix: 5
204:     args: -dm_plex_simplex 0 -dm_refine 1 -num_dof 1,0,0
205:   test:
206:     suffix: 6
207:     requires: ctetgen
208:     args: -dm_plex_dim 3 -dm_plex_simplex 1 -dm_refine_volume_limit_pre 0.00625 -num_dof 1,0,0,0
209:   test:
210:     suffix: 7
211:     args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_refine 1 -num_dof 1,0,0,0
212:   # Parallel tests
213:   # Grouping tests
214:   test:
215:     suffix: group_1
216:     args: -num_groups 1 -num_dof 1,0,0 -is_view -orig_mat_view -perm_mat_view
217:   test:
218:     suffix: group_2
219:     args: -num_groups 2 -num_dof 1,0,0 -is_view -perm_mat_view

221: TEST*/