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