1 static char help[] = "Element closure restrictions in tensor/lexicographic/spectral-element ordering using DMPlex\n\n"; 2 3 #include <petscdmplex.h> 4 5 static PetscErrorCode ViewOffsets(DM dm, Vec X) 6 { 7 PetscInt num_elem, elem_size, num_comp, num_dof; 8 PetscInt *elem_restr_offsets; 9 const PetscScalar *x = NULL; 10 const char *name; 11 12 PetscFunctionBegin; 13 PetscCall(PetscObjectGetName((PetscObject)dm, &name)); 14 PetscCall(DMPlexGetLocalOffsets(dm, NULL, 0, 0, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets)); 15 PetscCall(PetscPrintf(PETSC_COMM_SELF, "DM %s offsets: num_elem %" PetscInt_FMT ", size %" PetscInt_FMT ", comp %" PetscInt_FMT ", dof %" PetscInt_FMT "\n", name, num_elem, elem_size, num_comp, num_dof)); 16 if (X) PetscCall(VecGetArrayRead(X, &x)); 17 for (PetscInt c = 0; c < num_elem; c++) { 18 PetscCall(PetscIntView(elem_size, &elem_restr_offsets[c * elem_size], PETSC_VIEWER_STDOUT_SELF)); 19 if (x) { 20 for (PetscInt i = 0; i < elem_size; i++) PetscCall(PetscScalarView(num_comp, &x[elem_restr_offsets[c * elem_size + i]], PETSC_VIEWER_STDERR_SELF)); 21 } 22 } 23 if (X) PetscCall(VecRestoreArrayRead(X, &x)); 24 PetscCall(PetscFree(elem_restr_offsets)); 25 PetscFunctionReturn(PETSC_SUCCESS); 26 } 27 28 int main(int argc, char **argv) 29 { 30 DM dm; 31 PetscSection section; 32 PetscFE fe; 33 PetscInt dim, Nf = 1, c, cStart, cEnd; 34 PetscBool view_coord = PETSC_FALSE, tensor = PETSC_TRUE, project = PETSC_FALSE; 35 36 PetscFunctionBeginUser; 37 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 38 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Tensor closure restrictions", "DMPLEX"); 39 PetscCall(PetscOptionsBool("-closure_tensor", "Apply DMPlexSetClosurePermutationTensor", "ex8.c", tensor, &tensor, NULL)); 40 PetscCall(PetscOptionsBool("-project_coordinates", "Call DMSetCoordinateDisc() explicitly", "ex8.c", project, &project, NULL)); 41 PetscCall(PetscOptionsBool("-view_coord", "View coordinates of element closures", "ex8.c", view_coord, &view_coord, NULL)); 42 PetscCall(PetscOptionsBoundedInt("-num_fields", "The number of fields to use", "ex8.c", Nf, &Nf, NULL, 1)); 43 PetscOptionsEnd(); 44 45 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 46 PetscCall(DMSetType(dm, DMPLEX)); 47 PetscCall(DMSetFromOptions(dm)); 48 if (project) { 49 PetscFE fe_coords; 50 PetscInt cdim; 51 PetscCall(DMGetCoordinateDim(dm, &cdim)); 52 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, cdim, cdim, PETSC_FALSE, 1, 1, &fe_coords)); 53 PetscCall(DMSetCoordinateDisc(dm, fe_coords, PETSC_TRUE)); 54 PetscCall(PetscFEDestroy(&fe_coords)); 55 } 56 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 57 PetscCall(DMGetDimension(dm, &dim)); 58 59 if (Nf == 1) { 60 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, PETSC_DETERMINE, &fe)); 61 PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); 62 PetscCall(PetscFEDestroy(&fe)); 63 } else { 64 for (PetscInt f = 0; f < Nf; ++f) { 65 char prefix[16]; 66 PetscCall(PetscSNPrintf(prefix, 16, "f%" PetscInt_FMT "_", f)); 67 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, prefix, PETSC_DETERMINE, &fe)); 68 PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); 69 PetscCall(PetscFEDestroy(&fe)); 70 } 71 } 72 PetscCall(DMCreateDS(dm)); 73 if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL)); 74 PetscCall(DMGetLocalSection(dm, §ion)); 75 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 76 for (c = cStart; c < cEnd; c++) { 77 PetscInt numindices, *indices; 78 PetscCall(DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, &numindices, &indices, NULL, NULL)); 79 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT "\n", c - cStart)); 80 PetscCall(PetscIntView(numindices, indices, PETSC_VIEWER_STDOUT_SELF)); 81 PetscCall(DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &numindices, &indices, NULL, NULL)); 82 } 83 if (view_coord) { 84 DM cdm; 85 Vec X; 86 PetscInt cdim; 87 88 PetscCall(DMGetCoordinatesLocalSetUp(dm)); 89 PetscCall(DMGetCoordinateDim(dm, &cdim)); 90 PetscCall(DMGetCoordinateDM(dm, &cdm)); 91 PetscCall(PetscObjectSetName((PetscObject)cdm, "coords")); 92 if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL)); 93 for (c = cStart; c < cEnd; ++c) { 94 const PetscScalar *array; 95 PetscScalar *x = NULL; 96 PetscInt ndof; 97 PetscBool isDG; 98 99 PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &ndof, &array, &x)); 100 PetscCheck(ndof % cdim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "ndof not divisible by cdim"); 101 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT " coordinates\n", c - cStart)); 102 for (PetscInt i = 0; i < ndof; i += cdim) PetscCall(PetscScalarView(cdim, &x[i], PETSC_VIEWER_STDOUT_SELF)); 103 PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &ndof, &array, &x)); 104 } 105 PetscCall(ViewOffsets(dm, NULL)); 106 PetscCall(DMGetCoordinatesLocal(dm, &X)); 107 PetscCall(ViewOffsets(cdm, X)); 108 } 109 PetscCall(DMDestroy(&dm)); 110 PetscCall(PetscFinalize()); 111 return 0; 112 } 113 114 /*TEST 115 116 test: 117 suffix: 1d_q2 118 args: -dm_plex_dim 1 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2 119 test: 120 suffix: 2d_q1 121 args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 122 test: 123 suffix: 2d_q2_q1 124 args: -dm_plex_dim 2 -num_fields 2 -f0_petscspace_degree 2 -f1_petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 125 test: 126 suffix: 2d_q1_discontinuous 127 args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 -petscdualspace_lagrange_continuity 0 128 test: 129 suffix: 2d_q1_q1_discontinuous 130 args: -dm_plex_dim 2 -num_fields 2 -f0_petscspace_degree 1 -f1_petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -f1_petscdualspace_lagrange_continuity 0 131 test: 132 suffix: 2d_q2_p0_discontinuous 133 args: -dm_plex_dim 2 -num_fields 2 -f0_petscspace_degree 2 -f1_petscspace_degree 0 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -f1_petscdualspace_lagrange_continuity 0 134 test: 135 suffix: 2d_q2 136 args: -dm_plex_dim 2 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 137 test: 138 suffix: 2d_q3 139 args: -dm_plex_dim 2 -petscspace_degree 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 140 test: 141 suffix: 3d_q1 142 args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 143 test: 144 suffix: 1d_q1_periodic 145 requires: !complex 146 args: -dm_plex_dim 1 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 3 -dm_plex_box_bd periodic -dm_view -view_coord 147 test: 148 suffix: 2d_q1_periodic 149 requires: !complex 150 args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 3,2 -dm_plex_box_bd periodic,none -dm_view -view_coord 151 test: 152 suffix: 3d_q1_periodic 153 requires: !complex 154 args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 3,2,1 -dm_plex_box_bd periodic,none,none -dm_view -view_coord 155 test: 156 suffix: 3d_q1_periodic_project 157 requires: !complex 158 args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,3 -dm_plex_box_bd none,none,periodic -dm_view -view_coord -project_coordinates 159 160 test: 161 suffix: 3d_q2_periodic # not actually periodic because only 2 cells 162 args: -dm_plex_dim 3 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2,2 -dm_plex_box_bd periodic,none,periodic -dm_view 163 164 TEST*/ 165