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