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