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:   PetscFunctionBegin;
 20:   options->numFields     = 1;
 21:   options->numComponents = NULL;
 22:   options->numDof        = NULL;
 23:   options->numGroups     = 0;

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

 38: PetscErrorCode CleanupContext(AppCtx *user)
 39: {
 40:   PetscFunctionBegin;
 41:   PetscCall(PetscFree(user->numComponents));
 42:   PetscCall(PetscFree(user->numDof));
 43:   PetscFunctionReturn(PETSC_SUCCESS);
 44: }

 46: /* This mesh comes from~\cite{saad2003}, Fig. 2.10, p. 70. */
 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};

 52:   PetscFunctionBegin;
 53:   PetscCall(DMPlexCreateFromCellListPetsc(comm, 2, 16, 15, 3, PETSC_FALSE, cells, 2, coords, dm));
 54:   PetscFunctionReturn(PETSC_SUCCESS);
 55: }

 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;

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

 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;

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

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;

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

136: int main(int argc, char **argv)
137: {
138:   DM           dm;
139:   PetscSection s;
140:   AppCtx       user;
141:   PetscInt     dim;

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;

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

186: /*TEST

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

227: TEST*/