1# libCEED: Solid Mechanics Example 2 3This page provides a description of the solid mechanics example for the libCEED library, based on PETSc. 4 5Ratel, a more fully featured solid mechanics library, can be found on [GitLab](https://gitlab.com/micromorph/ratel). 6 7This code solves the steady-state static momentum balance equations using unstructured high-order finite/spectral element spatial discretizations. 8In this mini-app, we consider three formulations used in solid mechanics applications: linear elasticity, Neo-Hookean hyperelasticity at small strain, and Neo-Hookean hyperelasticity at finite strain. 9All three of these formulations are for compressible materials. 10 11Build by using: 12 13``` 14make 15``` 16 17and run with: 18 19``` 20./elasticity -mesh [.exo file] -degree [degree] -nu [nu] -E [E] [boundary options] -problem [problem type] -forcing [forcing] -ceed [ceed] 21``` 22 23## Runtime options 24 25<!-- % inclusion-solids-marker--> 26 27The elasticity mini-app is controlled via command-line options, the following of which are mandatory. 28 29:::{list-table} Mandatory Runtime Options 30:header-rows: 1 31:widths: 3 7 32 33* - Option 34 - Description 35* - `-mesh [filename]` 36 - Path to mesh file in any format supported by PETSc. 37* - `-degree [int]` 38 - Polynomial degree of the finite element basis 39* - `-E [real]` 40 - [Young's modulus](https://en.wikipedia.org/wiki/Young%27s_modulus), $E > 0$ 41* - `-nu [real]` 42 - [Poisson's ratio](https://en.wikipedia.org/wiki/Poisson%27s_ratio), $\nu < 0.5$ 43* - `-bc_clamp [int list]` 44 - List of face sets on which to displace by `-bc_clamp_[facenumber]_translate [x,y,z]` and/or `bc_clamp_[facenumber]_rotate [rx,ry,rz,c_0,c_1]`. 45 Note: The default for a clamped face is zero displacement. 46 All displacement is with respect to the initial configuration. 47* - `-bc_traction [int list]` 48 - List of face sets on which to set traction boundary conditions with the traction vector `-bc_traction_[facenumber] [tx,ty,tz]` 49::: 50 51:::{note} 52This solver can use any mesh format that PETSc's `DMPlex` can read (Exodus, Gmsh, Med, etc.). 53Our tests have primarily been using Exodus meshes created using [CUBIT]; sample meshes used for the example runs suggested here can be found in [this repository]. 54Note that many mesh formats require PETSc to be configured appropriately; e.g., `--download-exodusii` for Exodus support. 55::: 56 57Consider the specific example of the mesh seen below: 58 59```{image} https://github.com/jeremylt/ceedSampleMeshes/raw/master/cylinderDiagram.png 60``` 61 62With the sidesets defined in the figure, we provide here an example of a minimal set of command line options: 63 64``` 65./elasticity -mesh [.exo file] -degree 4 -E 1e6 -nu 0.3 -bc_clamp 998,999 -bc_clamp_998_translate 0,-0.5,1 66``` 67 68In this example, we set the left boundary, face set $999$, to zero displacement and the right boundary, face set $998$, to displace $0$ in the $x$ direction, $-0.5$ in the $y$, and $1$ in the $z$. 69 70As an alternative to specifying a mesh with {code}`-mesh`, the user may use a DMPlex box mesh by specifying {code}`-dm_plex_box_faces [int list]`, {code}`-dm_plex_box_upper [real list]`, and {code}`-dm_plex_box_lower [real list]`. 71 72As an alternative example exploiting {code}`-dm_plex_box_faces`, we consider a {code}`4 x 4 x 4` mesh where essential (Drichlet) boundary condition is placed on all sides. 73Sides 1 through 6 are rotated around $x$-axis: 74 75``` 76./elasticity -problem FS-NH -E 1 -nu 0.3 -num_steps 40 -snes_linesearch_type cp -dm_plex_box_faces 4,4,4 -bc_clamp 1,2,3,4,5,6 -bc_clamp_1_rotate 0,0,1,0,.3 -bc_clamp_2_rotate 0,0,1,0,.3 -bc_clamp_3_rotate 0,0,1,0,.3 -bc_clamp_4_rotate 0,0,1,0,.3 -bc_clamp_5_rotate 0,0,1,0,.3 -bc_clamp_6_rotate 0,0,1,0,.3 77``` 78 79:::{note} 80If the coordinates for a particular side of a mesh are zero along the axis of rotation, it may appear that particular side is clamped zero. 81::: 82 83On each boundary node, the rotation magnitude is computed: {code}`theta = (c_0 + c_1 * cx) * loadIncrement` where {code}`cx = kx * x + ky * y + kz * z`, with {code}`kx`, {code}`ky`, {code}`kz` are normalized values. 84 85The command line options just shown are the minimum requirements to run the mini-app, but additional options may also be set as follows 86 87:::{list-table} Additional Runtime Options 88:header-rows: 1 89 90* - Option 91 - Description 92 - Default value 93 94* - `-ceed` 95 - CEED resource specifier 96 - `/cpu/self` 97 98* - `-q_extra` 99 - Number of extra quadrature points 100 - `0` 101 102* - `-test` 103 - Run in test mode 104 - 105 106* - `-problem` 107 - Problem to solve (`Linear`, `FS-NH`, `FS-MR`, etc.) 108 - `Linear` 109 110* - `-forcing` 111 - Forcing term option (`none`, `constant`, or `mms`) 112 - `none` 113 114* - `-forcing_vec` 115 - Forcing vector 116 - `0,-1,0` 117 118* - `-multigrid` 119 - Multigrid coarsening to use (`logarithmic`, `uniform` or `none`) 120 - `logarithmic` 121 122* - `-nu_smoother [real]` 123 - Poisson's ratio for multigrid smoothers, $\nu < 0.5$ 124 - 125 126* - `-num_steps` 127 - Number of load increments for continuation method 128 - `1` if `Linear` else `10` 129 130* - `-view_soln` 131 - Output solution at each load increment for viewing 132 - 133 134* - `-view_final_soln` 135 - Output solution at final load increment for viewing 136 - 137 138* - `-snes_view` 139 - View PETSc `SNES` nonlinear solver configuration 140 - 141 142* - `-log_view` 143 - View PETSc performance log 144 - 145 146* - `-output_dir` 147 - Output directory 148 - `.` 149 150* - `-help` 151 - View comprehensive information about run-time options 152 - 153::: 154 155To verify the convergence of the linear elasticity formulation on a given mesh with the method of manufactured solutions, run: 156 157``` 158./elasticity -mesh [mesh] -degree [degree] -nu [nu] -E [E] -forcing mms 159``` 160 161This option attempts to recover a known solution from an analytically computed forcing term. 162 163### On algebraic solvers 164 165This mini-app is configured to use the following Newton-Krylov-Multigrid method by default. 166 167- Newton-type methods for the nonlinear solve, with the hyperelasticity models globalized using load increments. 168- Preconditioned conjugate gradients to solve the symmetric positive definite linear systems arising at each Newton step. 169- Preconditioning via $p$-version multigrid coarsening to linear elements, with algebraic multigrid (PETSc's `GAMG`) for the coarse solve. 170 The default smoother uses degree 3 Chebyshev with Jacobi preconditioning. 171 (Lower degree is often faster, albeit less robust; try {code}`-outer_mg_levels_ksp_max_it 2`, for example.) 172 Application of the linear operators for all levels with degree $p > 1$ is performed matrix-free using analytic Newton linearization, while the lowest order $p = 1$ operators are assembled explicitly (using coloring at present). 173 174Many related solvers can be implemented by composing PETSc command-line options. 175 176### Nondimensionalization 177 178Quantities such as the Young's modulus vary over many orders of magnitude, and thus can lead to poorly scaled equations. 179One can nondimensionalize the model by choosing an alternate system of units, such that displacements and residuals are of reasonable scales. 180 181:::{list-table} (Non)dimensionalization options 182:header-rows: 1 183 184* - Option 185 - Description 186 - Default value 187 188* - `-units_meter` 189 - 1 meter in scaled length units 190 - `1` 191 192* - `-units_second` 193 - 1 second in scaled time units 194 - `1` 195 196* - `-units_kilogram` 197 - 1 kilogram in scaled mass units 198 - `1` 199::: 200 201For example, consider a problem involving metals subject to gravity. 202 203:::{list-table} Characteristic units for metals 204:header-rows: 1 205 206* - Quantity 207 - Typical value in SI units 208 209* - Displacement, $\bm u$ 210 - $1 \,\mathrm{cm} = 10^{-2} \,\mathrm m$ 211 212* - Young's modulus, $E$ 213 - $10^{11} \,\mathrm{Pa} = 10^{11} \,\mathrm{kg}\, \mathrm{m}^{-1}\, \mathrm s^{-2}$ 214 215* - Body force (gravity) on volume, $\int \rho \bm g$ 216 - $5 \cdot 10^4 \,\mathrm{kg}\, \mathrm m^{-2} \, \mathrm s^{-2} \cdot (\text{volume} \, \mathrm m^3)$ 217::: 218 219One can choose units of displacement independently (e.g., {code}`-units_meter 100` to measure displacement in centimeters), but $E$ and $\int \rho \bm g$ have the same dependence on mass and time, so cannot both be made of order 1. 220This reflects the fact that both quantities are not equally significant for a given displacement size; the relative significance of gravity increases as the domain size grows. 221 222### Diagnostic Quantities 223 224Diagnostic quantities for viewing are provided when the command line options for visualization output, {code}`-view_soln` or {code}`-view_final_soln` are used. 225The diagnostic quantities include displacement in the $x$ direction, displacement in the $y$ direction, displacement in the $z$ direction, pressure, $\operatorname{trace} \bm{E}$, $\operatorname{trace} \bm{E}^2$, $\lvert J \rvert$, and strain energy density. 226The table below summarizes the formulations of each of these quantities for each problem type. 227 228:::{list-table} Diagnostic quantities 229 :header-rows: 1 230 231 * - Quantity 232 - Linear Elasticity 233 - Hyperelasticity, Small Strain 234 - Hyperelasticity, Finite Strain 235 236 * - Pressure 237 - $\lambda \operatorname{trace} \bm{\epsilon}$ 238 - $\lambda \log \operatorname{trace} \bm{\epsilon}$ 239 - $\lambda \log J$ 240 241 * - Volumetric Strain 242 - $\operatorname{trace} \bm{\epsilon}$ 243 - $\operatorname{trace} \bm{\epsilon}$ 244 - $\operatorname{trace} \bm{E}$ 245 246 * - $\operatorname{trace} \bm{E}^2$ 247 - $\operatorname{trace} \bm{\epsilon}^2$ 248 - $\operatorname{trace} \bm{\epsilon}^2$ 249 - $\operatorname{trace} \bm{E}^2$ 250 251 * - $\lvert J \rvert$ 252 - $1 + \operatorname{trace} \bm{\epsilon}$ 253 - $1 + \operatorname{trace} \bm{\epsilon}$ 254 - $\lvert J \rvert$ 255 256 * - Strain Energy Density 257 - $\frac{\lambda}{2} (\operatorname{trace} \bm{\epsilon})^2 + \mu \bm{\epsilon} : \bm{\epsilon}$ 258 - $\lambda (1 + \operatorname{trace} \bm{\epsilon}) (\log(1 + \operatorname{trace} \bm{\epsilon} ) - 1) + \mu \bm{\epsilon} : \bm{\epsilon}$ 259 - $\frac{\lambda}{2}(\log J)^2 + \mu \operatorname{trace} \bm{E} - \mu \log J$ 260::: 261 262[cubit]: https://cubit.sandia.gov/ 263[this repository]: https://github.com/jeremylt/ceedSampleMeshes 264
