libCEED: Solid Mechanics Example

This page provides a description of the solid mechanics example for the libCEED library, based on PETSc. PETSc v3.17 or a development version of PETSc at commit 0e95d842 or later is required.

This code solves the steady-state static momentum balance equations using unstructured high-order finite/spectral element spatial discretizations. In 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. All three of these formulations are for compressible materials.

Build by using:

make

and run with:

./elasticity -mesh [.exo file] -degree [degree] -nu [nu] -E [E] [boundary options] -problem [problem type] -forcing [forcing] -ceed [ceed]

Runtime options

The elasticity mini-app is controlled via command-line options, the following of which are mandatory.

Mandatory Runtime Options

Option	Description
-mesh [filename]	Path to mesh file in any format supported by PETSc.
-degree [int]	Polynomial degree of the finite element basis
-E [real]	Young’s modulus, \(E > 0\)
-nu [real]	Poisson’s ratio, \(\nu < 0.5\)
-bc_clamp [int list]	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]. Note: The default for a clamped face is zero displacement. All displacement is with respect to the initial configuration.
-bc_traction [int list]	List of face sets on which to set traction boundary conditions with the traction vector -bc_traction_[facenumber] [tx,ty,tz]

Note

This solver can use any mesh format that PETSc’s DMPlex can read (Exodus, Gmsh, Med, etc.). Our 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. Note that many mesh formats require PETSc to be configured appropriately; e.g., --download-exodusii for Exodus support.

Consider the specific example of the mesh seen below:

With the sidesets defined in the figure, we provide here an example of a minimal set of command line options:

./elasticity -mesh [.exo file] -degree 4 -E 1e6 -nu 0.3 -bc_clamp 998,999 -bc_clamp_998_translate 0,-0.5,1

In 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\).

As an alternative to specifying a mesh with -mesh, the user may use a DMPlex box mesh by specifying -dm_plex_box_faces [int list], -dm_plex_box_upper [real list], and -dm_plex_box_lower [real list].

As an alternative example exploiting -dm_plex_box_faces, we consider a 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:

./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

Note

If 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.

On each boundary node, the rotation magnitude is computed: theta = (c_0 + c_1 * cx) * loadIncrement where cx = kx * x + ky * y + kz * z, with kx, ky, kz are normalized values.

The command line options just shown are the minimum requirements to run the mini-app, but additional options may also be set as follows

Additional Runtime Options

Option	Description	Default value
-ceed	CEED resource specifier	/cpu/self
-q_extra	Number of extra quadrature points	0
-test	Run in test mode
-problem	Problem to solve (Linear, SS-NH, FSInitial-NH1, etc.)	Linear
-forcing	Forcing term option (none, constant, or mms)	none
-forcing_vec	Forcing vector	0,-1,0
-multigrid	Multigrid coarsening to use (logarithmic, uniform or none)	logarithmic
-nu_smoother [real]	Poisson’s ratio for multigrid smoothers, \(\nu < 0.5\)
-num_steps	Number of load increments for continuation method	1 if Linear else 10
-view_soln	Output solution at each load increment for viewing
-view_final_soln	Output solution at final load increment for viewing
-snes_view	View PETSc SNES nonlinear solver configuration
-log_view	View PETSc performance log
-output_dir	Output directory	.
-help	View comprehensive information about run-time options

To verify the convergence of the linear elasticity formulation on a given mesh with the method of manufactured solutions, run:

./elasticity -mesh [mesh] -degree [degree] -nu [nu] -E [E] -forcing mms

This option attempts to recover a known solution from an analytically computed forcing term.

On algebraic solvers

This mini-app is configured to use the following Newton-Krylov-Multigrid method by default.

• Newton-type methods for the nonlinear solve, with the hyperelasticity models globalized using load increments.

• Preconditioned conjugate gradients to solve the symmetric positive definite linear systems arising at each Newton step.

• Preconditioning via \(p\)-version multigrid coarsening to linear elements, with algebraic multigrid (PETSc’s GAMG) for the coarse solve. The default smoother uses degree 3 Chebyshev with Jacobi preconditioning. (Lower degree is often faster, albeit less robust; try -outer_mg_levels_ksp_max_it 2, for example.) 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).

Many related solvers can be implemented by composing PETSc command-line options.

Nondimensionalization

Quantities such as the Young’s modulus vary over many orders of magnitude, and thus can lead to poorly scaled equations. One can nondimensionalize the model by choosing an alternate system of units, such that displacements and residuals are of reasonable scales.

(Non)dimensionalization options

Option	Description	Default value
-units_meter	1 meter in scaled length units	1
-units_second	1 second in scaled time units	1
-units_kilogram	1 kilogram in scaled mass units	1

For example, consider a problem involving metals subject to gravity.

Characteristic units for metals

Quantity	Typical value in SI units
Displacement, \(\bm u\)	\(1 \,\mathrm{cm} = 10^{-2} \,\mathrm m\)
Young’s modulus, \(E\)	\(10^{11} \,\mathrm{Pa} = 10^{11} \,\mathrm{kg}\, \mathrm{m}^{-1}\, \mathrm s^{-2}\)
Body force (gravity) on volume, \(\int \rho \bm g\)	\(5 \cdot 10^4 \,\mathrm{kg}\, \mathrm m^{-2} \, \mathrm s^{-2} \cdot (\text{volume} \, \mathrm m^3)\)

One can choose units of displacement independently (e.g., -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. This 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.

Diagnostic Quantities

Diagnostic quantities for viewing are provided when the command line options for visualization output, -view_soln or -view_final_soln are used. The 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. The table below summarizes the formulations of each of these quantities for each problem type.

Diagnostic quantities

Quantity	Linear Elasticity	Hyperelasticity, Small Strain	Hyperelasticity, Finite Strain
Pressure	\(\lambda \operatorname{trace} \bm{\epsilon}\)	\(\lambda \log \operatorname{trace} \bm{\epsilon}\)	\(\lambda \log J\)
Volumetric Strain	\(\operatorname{trace} \bm{\epsilon}\)	\(\operatorname{trace} \bm{\epsilon}\)	\(\operatorname{trace} \bm{E}\)
\(\operatorname{trace} \bm{E}^2\)	\(\operatorname{trace} \bm{\epsilon}^2\)	\(\operatorname{trace} \bm{\epsilon}^2\)	\(\operatorname{trace} \bm{E}^2\)
\(\lvert J \rvert\)	\(1 + \operatorname{trace} \bm{\epsilon}\)	\(1 + \operatorname{trace} \bm{\epsilon}\)	\(\lvert J \rvert\)
Strain Energy Density	\(\frac{\lambda}{2} (\operatorname{trace} \bm{\epsilon})^2 + \mu \bm{\epsilon} : \bm{\epsilon}\)	\(\lambda (1 + \operatorname{trace} \bm{\epsilon}) (\log(1 + \operatorname{trace} \bm{\epsilon} ) - 1) + \mu \bm{\epsilon} : \bm{\epsilon}\)	\(\frac{\lambda}{2}(\log J)^2 + \mu \operatorname{trace} \bm{E} - \mu \log J\)