Actual source code: ex8.c
1: static char help[] = "Element closure restrictions in tensor/lexicographic/spectral-element ordering using DMPlex\n\n";
3: #include <petscdmplex.h>
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;
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: }
28: int main(int argc, char **argv)
29: {
30: DM dm;
31: PetscSection section;
32: PetscFE fe;
33: PetscInt dim, c, cStart, cEnd;
34: PetscBool view_coord = PETSC_FALSE, tensor = PETSC_TRUE, project = PETSC_FALSE;
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 DMProjectCoordinates explicitly", "ex8.c", project, &project, NULL));
41: PetscCall(PetscOptionsBool("-view_coord", "View coordinates of element closures", "ex8.c", view_coord, &view_coord, NULL));
42: PetscOptionsEnd();
44: PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
45: PetscCall(DMSetType(dm, DMPLEX));
46: PetscCall(DMSetFromOptions(dm));
47: if (project) {
48: PetscFE fe_coords;
49: PetscInt cdim;
50: PetscCall(DMGetCoordinateDim(dm, &cdim));
51: PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, cdim, cdim, PETSC_FALSE, 1, 1, &fe_coords));
52: PetscCall(DMProjectCoordinates(dm, fe_coords));
53: PetscCall(PetscFEDestroy(&fe_coords));
54: }
55: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
56: PetscCall(DMGetDimension(dm, &dim));
58: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, PETSC_DETERMINE, &fe));
59: PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
60: PetscCall(DMCreateDS(dm));
61: if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL));
62: PetscCall(DMGetLocalSection(dm, §ion));
63: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
64: for (c = cStart; c < cEnd; c++) {
65: PetscInt numindices, *indices;
66: PetscCall(DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, &numindices, &indices, NULL, NULL));
67: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT "\n", c - cStart));
68: PetscCall(PetscIntView(numindices, indices, PETSC_VIEWER_STDOUT_SELF));
69: PetscCall(DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &numindices, &indices, NULL, NULL));
70: }
71: if (view_coord) {
72: DM cdm;
73: Vec X;
74: PetscInt cdim;
76: PetscCall(DMGetCoordinatesLocalSetUp(dm));
77: PetscCall(DMGetCoordinateDim(dm, &cdim));
78: PetscCall(DMGetCoordinateDM(dm, &cdm));
79: PetscCall(PetscObjectSetName((PetscObject)cdm, "coords"));
80: if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL));
81: for (c = cStart; c < cEnd; ++c) {
82: const PetscScalar *array;
83: PetscScalar *x = NULL;
84: PetscInt ndof;
85: PetscBool isDG;
87: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &ndof, &array, &x));
88: PetscCheck(ndof % cdim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "ndof not divisible by cdim");
89: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT " coordinates\n", c - cStart));
90: for (PetscInt i = 0; i < ndof; i += cdim) PetscCall(PetscScalarView(cdim, &x[i], PETSC_VIEWER_STDOUT_SELF));
91: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &ndof, &array, &x));
92: }
93: PetscCall(ViewOffsets(dm, NULL));
94: PetscCall(DMGetCoordinatesLocal(dm, &X));
95: PetscCall(ViewOffsets(cdm, X));
96: }
97: PetscCall(PetscFEDestroy(&fe));
98: PetscCall(DMDestroy(&dm));
99: PetscCall(PetscFinalize());
100: return 0;
101: }
103: /*TEST
105: test:
106: suffix: 1d_q2
107: args: -dm_plex_dim 1 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2
108: test:
109: suffix: 2d_q1
110: args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2
111: test:
112: suffix: 2d_q2
113: args: -dm_plex_dim 2 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2
114: test:
115: suffix: 2d_q3
116: args: -dm_plex_dim 2 -petscspace_degree 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1
117: test:
118: suffix: 3d_q1
119: args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
120: test:
121: suffix: 1d_q1_periodic
122: requires: !complex
123: 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
124: test:
125: suffix: 2d_q1_periodic
126: requires: !complex
127: 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
128: test:
129: suffix: 3d_q1_periodic
130: requires: !complex
131: 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
132: test:
133: suffix: 3d_q1_periodic_project
134: requires: !complex
135: 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
137: test:
138: suffix: 3d_q2_periodic # not actually periodic because only 2 cells
139: 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
141: TEST*/