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