1 static char help[] = "Test degenerate near null space";
2
3 #include <petscdmplex.h>
4 #include <petscsnes.h>
5 #include <petscds.h>
6 #include <petscbag.h>
7 #include <petscconvest.h>
8
CreateMesh(MPI_Comm comm,DM * dm)9 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
10 {
11 PetscFunctionBeginUser;
12 PetscCall(DMCreate(comm, dm));
13 PetscCall(DMSetType(*dm, DMPLEX));
14 PetscCall(DMSetFromOptions(*dm));
15 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
16 PetscFunctionReturn(PETSC_SUCCESS);
17 }
18
SetupBoundaries(DM dm)19 static PetscErrorCode SetupBoundaries(DM dm)
20 {
21 DMLabel label;
22 PetscInt id;
23 PetscInt dim;
24
25 PetscFunctionBeginUser;
26 PetscCall(DMGetCoordinateDim(dm, &dim));
27 PetscCall(DMGetLabel(dm, "Face Sets", &label));
28 PetscCheck(label, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_INCOMP, "Must have face sets label");
29
30 if (dim == 2) {
31 PetscInt cmp, cmps_y[] = {0, 1};
32
33 cmp = 0;
34 id = 4;
35 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left", label, 1, &id, 0, 1, &cmp, NULL, NULL, NULL, NULL));
36 cmp = 0;
37 id = 2;
38 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right", label, 1, &id, 0, 1, &cmp, NULL, NULL, NULL, NULL));
39 id = 1;
40 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom", label, 1, &id, 0, 2, cmps_y, NULL, NULL, NULL, NULL));
41 } else if (dim == 3) {
42 PetscInt cmps_xy[] = {0, 1};
43 PetscInt cmps_z[] = {0, 1, 2};
44
45 id = 6;
46 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left", label, 1, &id, 0, 2, cmps_xy, NULL, NULL, NULL, NULL));
47 id = 5;
48 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right", label, 1, &id, 0, 2, cmps_xy, NULL, NULL, NULL, NULL));
49 id = 2;
50 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom", label, 1, &id, 0, 3, cmps_z, NULL, NULL, NULL, NULL));
51 id = 3;
52 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "front", label, 1, &id, 0, 2, cmps_xy, NULL, NULL, NULL, NULL));
53 id = 4;
54 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "back", label, 1, &id, 0, 2, cmps_xy, NULL, NULL, NULL, NULL));
55 }
56 PetscFunctionReturn(PETSC_SUCCESS);
57 }
58
SetupFE(DM dm,const char name[])59 PetscErrorCode SetupFE(DM dm, const char name[])
60 {
61 PetscFE fe;
62 char prefix[PETSC_MAX_PATH_LEN];
63 DMPolytopeType ct;
64 PetscInt dim, cStart;
65
66 PetscFunctionBegin;
67 /* Create finite element */
68 PetscCall(DMGetDimension(dm, &dim));
69 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
70 PetscCall(DMPlexGetCellType(dm, cStart, &ct));
71 PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name));
72 PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, dim, ct, name ? prefix : NULL, -1, &fe));
73 PetscCall(PetscObjectSetName((PetscObject)fe, name));
74 /* Set discretization and boundary conditions for each mesh */
75 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
76 PetscCall(DMCreateDS(dm));
77 PetscCall(SetupBoundaries(dm));
78 PetscCall(PetscFEDestroy(&fe));
79 PetscFunctionReturn(PETSC_SUCCESS);
80 }
81
MakeNullSpaceRigidBody(DM dm)82 PetscErrorCode MakeNullSpaceRigidBody(DM dm)
83 {
84 Mat A;
85 MatNullSpace null_space;
86
87 PetscFunctionBegin;
88 /* Create null space and set onto matrix */
89 PetscCall(DMCreateMatrix(dm, &A));
90 PetscCall(DMPlexCreateRigidBody(dm, 0, &null_space));
91 PetscCall(MatSetNearNullSpace(A, null_space));
92 PetscCall(MatNullSpaceView(null_space, PETSC_VIEWER_STDOUT_WORLD));
93 PetscCall(MatNullSpaceDestroy(&null_space));
94 PetscCall(MatDestroy(&A));
95 PetscFunctionReturn(PETSC_SUCCESS);
96 }
97
main(int argc,char ** argv)98 int main(int argc, char **argv)
99 {
100 DM dm;
101
102 PetscFunctionBeginUser;
103 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
104 PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm));
105 PetscCall(SetupFE(dm, "displacement"));
106 PetscCall(MakeNullSpaceRigidBody(dm));
107 PetscCall(DMDestroy(&dm));
108 PetscCall(PetscFinalize());
109 return 0;
110 }
111
112 /*TEST
113
114 test:
115 suffix: 3d_q1
116 args: -dm_plex_box_faces 1,1,2 -dm_plex_simplex 0 -dm_plex_dim 3 -displacement_petscspace_degree 1
117
118 test:
119 suffix: 3d_q2
120 args: -dm_plex_box_faces 1,1,2 -dm_plex_simplex 0 -dm_plex_dim 3 -displacement_petscspace_degree 2
121
122 test:
123 suffix: 2d_q1
124 args: -dm_plex_box_faces 1,2 -dm_plex_simplex 0 -dm_plex_dim 2 -displacement_petscspace_degree 1
125
126 test:
127 suffix: 2d_q2
128 args: -dm_plex_box_faces 1,2 -dm_plex_simplex 0 -dm_plex_dim 2 -displacement_petscspace_degree 2
129
130 TEST*/
131