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