1 /// @file 2 /// DM setup for solid mechanics example using PETSc 3 4 #include "../include/boundary.h" 5 #include "../include/setup-dm.h" 6 #include "../include/petsc-macros.h" 7 8 // ----------------------------------------------------------------------------- 9 // Setup DM 10 // ----------------------------------------------------------------------------- 11 PetscErrorCode CreateBCLabel(DM dm, const char name[]) { 12 PetscErrorCode ierr; 13 DMLabel label; 14 15 PetscFunctionBeginUser; 16 17 ierr = DMCreateLabel(dm, name); CHKERRQ(ierr); 18 ierr = DMGetLabel(dm, name, &label); CHKERRQ(ierr); 19 ierr = DMPlexMarkBoundaryFaces(dm, 1, label); CHKERRQ(ierr); 20 ierr = DMPlexLabelComplete(dm, label); CHKERRQ(ierr); 21 22 PetscFunctionReturn(0); 23 }; 24 25 // Read mesh and distribute DM in parallel 26 PetscErrorCode CreateDistributedDM(MPI_Comm comm, AppCtx app_ctx, DM *dm) { 27 PetscErrorCode ierr; 28 const char *filename = app_ctx->mesh_file; 29 PetscBool interpolate = PETSC_TRUE; 30 DM distributed_mesh = NULL; 31 PetscPartitioner part; 32 33 PetscFunctionBeginUser; 34 35 if (!*filename) { 36 PetscInt dim = 3, faces[3] = {3, 3, 3}; 37 ierr = PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", 38 faces, &dim, NULL); CHKERRQ(ierr); 39 if (!dim) dim = 3; 40 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, faces, NULL, 41 NULL, NULL, interpolate, dm); CHKERRQ(ierr); 42 } else { 43 ierr = DMPlexCreateFromFile(comm, filename, NULL, interpolate, dm); 44 CHKERRQ(ierr); 45 } 46 47 // Distribute DM in parallel 48 ierr = DMPlexGetPartitioner(*dm, &part); CHKERRQ(ierr); 49 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 50 ierr = DMPlexDistribute(*dm, 0, NULL, &distributed_mesh); CHKERRQ(ierr); 51 if (distributed_mesh) { 52 ierr = DMDestroy(dm); CHKERRQ(ierr); 53 *dm = distributed_mesh; 54 } 55 ierr = DMViewFromOptions(*dm, NULL, "-dm_view"); CHKERRQ(ierr); 56 57 PetscFunctionReturn(0); 58 }; 59 60 // Setup DM with FE space of appropriate degree 61 PetscErrorCode SetupDMByDegree(DM dm, AppCtx app_ctx, PetscInt order, 62 PetscBool boundary, PetscInt num_comp_u) { 63 PetscErrorCode ierr; 64 MPI_Comm comm; 65 PetscInt dim; 66 PetscFE fe; 67 IS face_set_is; // Index Set for Face Sets 68 const char *name = "Face Sets"; // PETSc internal requirement 69 PetscInt num_face_sets; // Number of FaceSets in face_set_is 70 const PetscInt *face_set_ids; // id of each FaceSet 71 72 PetscFunctionBeginUser; 73 74 // Setup DM 75 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 76 ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); 77 ierr = PetscFECreateLagrange(comm, dim, num_comp_u, PETSC_FALSE, order, order, 78 &fe); CHKERRQ(ierr); 79 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 80 ierr = DMAddField(dm, NULL, (PetscObject)fe); CHKERRQ(ierr); 81 ierr = DMCreateDS(dm); CHKERRQ(ierr); 82 { 83 /* create FE field for coordinates */ 84 PetscFE fe_coords; 85 PetscInt num_comp_coord; 86 ierr = DMGetCoordinateDim(dm, &num_comp_coord); CHKERRQ(ierr); 87 ierr = PetscFECreateLagrange(comm, dim, num_comp_coord, PETSC_FALSE, 1, 1, 88 &fe_coords); CHKERRQ(ierr); 89 ierr = DMProjectCoordinates(dm, fe_coords); CHKERRQ(ierr); 90 ierr = PetscFEDestroy(&fe_coords); CHKERRQ(ierr); 91 } 92 93 // Add Dirichlet (Essential) boundary 94 if (boundary) { 95 if (app_ctx->forcing_choice == FORCE_MMS) { 96 if (app_ctx->test_mode) { 97 // -- Test mode - box mesh 98 PetscBool has_label; 99 PetscInt marker_ids[1] = {1}; 100 ierr = DMHasLabel(dm, "marker", &has_label); CHKERRQ(ierr); 101 if (!has_label) { 102 ierr = CreateBCLabel(dm, "marker"); CHKERRQ(ierr); 103 } 104 DMLabel label; 105 ierr = DMGetLabel(dm, "marker", &label); CHKERRQ(ierr); 106 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "mms", label, "marker", 1, marker_ids, 107 0, 0, NULL, (void(*)(void))BCMMS, NULL, NULL, NULL); 108 CHKERRQ(ierr); 109 } else { 110 // -- ExodusII mesh with MMS 111 ierr = DMGetLabelIdIS(dm, name, &face_set_is); CHKERRQ(ierr); 112 ierr = ISGetSize(face_set_is,&num_face_sets); CHKERRQ(ierr); 113 ierr = ISGetIndices(face_set_is, &face_set_ids); CHKERRQ(ierr); 114 DMLabel label; 115 ierr = DMGetLabel(dm, "Face Sets", &label); CHKERRQ(ierr); 116 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "mms", label, "Face Sets", 117 num_face_sets, face_set_ids, 0, 0, NULL, 118 (void(*)(void))BCMMS, NULL, NULL, NULL); 119 CHKERRQ(ierr); 120 ierr = ISRestoreIndices(face_set_is, &face_set_ids); CHKERRQ(ierr); 121 ierr = ISDestroy(&face_set_is); CHKERRQ(ierr); 122 } 123 } else { 124 // -- Mesh with user specified BCs 125 DMLabel label; 126 ierr = DMGetLabel(dm, "Face Sets", &label); CHKERRQ(ierr); 127 // -- Clamp BCs 128 for (PetscInt i = 0; i < app_ctx->bc_clamp_count; i++) { 129 char bcName[25]; 130 snprintf(bcName, sizeof bcName, "clamp_%d", app_ctx->bc_clamp_faces[i]); 131 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, bcName, label, "Face Sets", 1, 132 &app_ctx->bc_clamp_faces[i], 0, 0, 133 NULL, (void(*)(void))BCClamp, NULL, 134 (void *)&app_ctx->bc_clamp_max[i], NULL); 135 CHKERRQ(ierr); 136 } 137 } 138 } 139 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 140 CHKERRQ(ierr); 141 142 // Cleanup 143 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 144 145 PetscFunctionReturn(0); 146 }; 147