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 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 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 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 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 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