1 // Copyright (c) 2017-2026, 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 // -----------------------------------------------------------------------------
CreateBCLabel(DM dm,const char name[])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
CreateDistributedDM(MPI_Comm comm,AppCtx app_ctx,DM * dm)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
SetupDMByDegree(DM dm,AppCtx app_ctx,PetscInt order,PetscBool boundary,PetscInt num_comp_u)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