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