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