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 PetscInt num_elem, elem_size, num_comp, num_dof; 7 PetscInt *elem_restr_offsets; 8 const PetscScalar *x = NULL; 9 const char *name; 10 11 PetscFunctionBegin; 12 PetscCall(PetscObjectGetName((PetscObject)dm, &name)); 13 PetscCall(DMPlexGetLocalOffsets(dm, NULL, 0, 0, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets)); 14 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)); 15 if (X) PetscCall(VecGetArrayRead(X, &x)); 16 for (PetscInt c = 0; c < num_elem; c++) { 17 PetscCall(PetscIntView(elem_size, &elem_restr_offsets[c * elem_size], PETSC_VIEWER_STDOUT_SELF)); 18 if (x) { 19 for (PetscInt i = 0; i < elem_size; i++) PetscCall(PetscScalarView(num_comp, &x[elem_restr_offsets[c * elem_size + i]], PETSC_VIEWER_STDERR_SELF)); 20 } 21 } 22 if (X) PetscCall(VecRestoreArrayRead(X, &x)); 23 PetscCall(PetscFree(elem_restr_offsets)); 24 PetscFunctionReturn(0); 25 } 26 27 int main(int argc, char **argv) { 28 DM dm; 29 PetscSection section; 30 PetscFE fe; 31 PetscInt dim, c, cStart, cEnd; 32 PetscBool view_coord = PETSC_FALSE, tensor = PETSC_TRUE, project = PETSC_FALSE; 33 34 PetscFunctionBeginUser; 35 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 36 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Tensor closure restrictions", "DMPLEX"); 37 PetscCall(PetscOptionsBool("-closure_tensor", "Apply DMPlexSetClosurePermutationTensor", "ex8.c", tensor, &tensor, NULL)); 38 PetscCall(PetscOptionsBool("-project_coordinates", "Call DMProjectCoordinates explicitly", "ex8.c", project, &project, NULL)); 39 PetscCall(PetscOptionsBool("-view_coord", "View coordinates of element closures", "ex8.c", view_coord, &view_coord, NULL)); 40 PetscOptionsEnd(); 41 42 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 43 PetscCall(DMSetType(dm, DMPLEX)); 44 PetscCall(DMSetFromOptions(dm)); 45 if (project) { 46 PetscFE fe_coords; 47 PetscInt cdim; 48 PetscCall(DMGetCoordinateDim(dm, &cdim)); 49 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, cdim, cdim, PETSC_FALSE, 1, 1, &fe_coords)); 50 PetscCall(DMProjectCoordinates(dm, fe_coords)); 51 PetscCall(PetscFEDestroy(&fe_coords)); 52 } 53 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 54 PetscCall(DMGetDimension(dm, &dim)); 55 56 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, PETSC_DETERMINE, &fe)); 57 PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); 58 PetscCall(DMCreateDS(dm)); 59 if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL)); 60 PetscCall(DMGetLocalSection(dm, §ion)); 61 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 62 for (c = cStart; c < cEnd; c++) { 63 PetscInt numindices, *indices; 64 PetscCall(DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, &numindices, &indices, NULL, NULL)); 65 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT "\n", c - cStart)); 66 PetscCall(PetscIntView(numindices, indices, PETSC_VIEWER_STDOUT_SELF)); 67 PetscCall(DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &numindices, &indices, NULL, NULL)); 68 } 69 if (view_coord) { 70 DM cdm; 71 Vec X; 72 PetscInt cdim; 73 74 PetscCall(DMGetCoordinateDim(dm, &cdim)); 75 PetscCall(DMGetCoordinateDM(dm, &cdm)); 76 PetscCall(PetscObjectSetName((PetscObject)cdm, "coords")); 77 if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL)); 78 for (c = cStart; c < cEnd; ++c) { 79 const PetscScalar *array; 80 PetscScalar *x = NULL; 81 PetscInt ndof; 82 PetscBool isDG; 83 84 PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &ndof, &array, &x)); 85 PetscCheck(ndof % cdim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "ndof not divisible by cdim"); 86 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT " coordinates\n", c - cStart)); 87 for (PetscInt i = 0; i < ndof; i += cdim) PetscCall(PetscScalarView(cdim, &x[i], PETSC_VIEWER_STDOUT_SELF)); 88 PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &ndof, &array, &x)); 89 } 90 PetscCall(ViewOffsets(dm, NULL)); 91 PetscCall(DMGetCoordinatesLocal(dm, &X)); 92 PetscCall(ViewOffsets(cdm, X)); 93 } 94 PetscCall(PetscFEDestroy(&fe)); 95 PetscCall(DMDestroy(&dm)); 96 PetscCall(PetscFinalize()); 97 return 0; 98 } 99 100 /*TEST 101 102 test: 103 suffix: 1d_q2 104 args: -dm_plex_dim 1 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2 105 test: 106 suffix: 2d_q1 107 args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 108 test: 109 suffix: 2d_q2 110 args: -dm_plex_dim 2 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 111 test: 112 suffix: 2d_q3 113 args: -dm_plex_dim 2 -petscspace_degree 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 114 test: 115 suffix: 3d_q1 116 args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 117 test: 118 suffix: 1d_q1_periodic 119 requires: !complex 120 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 121 test: 122 suffix: 2d_q1_periodic 123 requires: !complex 124 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 125 test: 126 suffix: 3d_q1_periodic 127 requires: !complex 128 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 129 test: 130 suffix: 3d_q1_periodic_project 131 requires: !complex 132 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 133 134 test: 135 suffix: 3d_q2_periodic # not actually periodic because only 2 cells 136 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 137 138 TEST*/ 139