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