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