Solid mechanics mini-app

This example is located in the subdirectory examples/solids. It solves the steady-state static momentum balance equations using unstructured high-order finite/spectral element spatial discretizations. As for the Compressible Navier-Stokes mini-app case, the solid mechanics elasticity example has been developed using PETSc, so that the pointwise physics (defined at quadrature points) is separated from the parallelization and meshing concerns.

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. We provide the strong and weak forms of static balance of linear momentum in the small strain and finite strain regimes. The stress-strain relationship (constitutive law) for each of the material models is provided. Due to the nonlinearity of material models in Neo-Hookean hyperelasticity, the Newton linearization of the material models is provided.


Linear elasticity and small-strain hyperelasticity can both by obtained from the finite-strain hyperelastic formulation by linearization of geometric and constitutive nonlinearities. The effect of these linearizations is sketched in the diagram below, where \(\bm \sigma\) and \(\bm \epsilon\) are stress and strain, respectively, in the small strain regime, while \(\bm S\) and \(\bm E\) are their finite-strain generalizations (second Piola-Kirchoff tensor and Green-Lagrange strain tensor, respectively) defined in the reference configuration, and \(\mathsf C\) is a linearized constitutive model.

(23)\[\begin{CD} {\overbrace{\bm S(\bm E)}^{\text{Finite Strain Hyperelastic}}} @>{\text{constitutive}}>{\text{linearization}}> {\overbrace{\bm S = \mathsf C \bm E}^{\text{St. Venant-Kirchoff}}} \\ @V{\text{geometric}}V{\begin{smallmatrix}\bm E \to \bm \epsilon \\ \bm S \to \bm \sigma \end{smallmatrix}}V @V{\begin{smallmatrix}\bm E \to \bm \epsilon \\ \bm S \to \bm \sigma \end{smallmatrix}}V{\text{geometric}}V \\ {\underbrace{\bm \sigma(\bm \epsilon)}_\text{Small Strain Hyperelastic}} @>{\text{constitutive}}>\text{linearization}> {\underbrace{\bm \sigma = \mathsf C \bm \epsilon}_\text{Linear Elastic}} \end{CD}\]

Running the mini-app

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

Table 6 Mandatory Runtime Options



-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]


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 hyperFS -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


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

Table 7 Additional Runtime Options



Default value


CEED resource specifier



Number of extra quadrature points



Run in test mode


Problem to solve (linElas, hyperSS or hyperFS)



Forcing term option (none, constant, or mms)



Forcing vector



Multigrid coarsening to use (logarithmic, uniform or none)


-nu_smoother [real]

Poisson’s ratio for multigrid smoothers, \(\nu < 0.5\)


Number of load increments for continuation method

1 if linElas else 10


Output solution at each load increment for viewing


Output solution at final load increment for viewing


View PETSc SNES nonlinear solver configuration


View PETSc performance log


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.


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.

Table 8 (Non)dimensionalization options



Default value


1 meter in scaled length units



1 second in scaled time units



1 kilogram in scaled mass units


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

Table 9 Characteristic units for metals


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.

Table 10 Diagnostic quantities


Linear Elasticity

Hyperelasticity, Small Strain

Hyperelasticity, Finite Strain


\(\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\)

Linear Elasticity

The strong form of the static balance of linear momentum at small strain for the three-dimensional linear elasticity problem is given by [Hug12]:

(24)\[\nabla \cdot \bm{\sigma} + \bm{g} = \bm{0}\]

where \(\bm{\sigma}\) and \(\bm{g}\) are stress and forcing functions, respectively. We multiply (24) by a test function \(\bm v\) and integrate the divergence term by parts to arrive at the weak form: find \(\bm u \in \mathcal V \subset H^1(\Omega)\) such that

(25)\[\int_{\Omega}{ \nabla \bm{v} \colon \bm{\sigma}} \, dV - \int_{\partial \Omega}{\bm{v} \cdot \left(\bm{\sigma} \cdot \hat{\bm{n}}\right)} \, dS - \int_{\Omega}{\bm{v} \cdot \bm{g}} \, dV = 0, \quad \forall \bm v \in \mathcal V,\]

where \(\bm{\sigma} \cdot \hat{\bm{n}}|_{\partial \Omega}\) is replaced by an applied force/traction boundary condition written in terms of the reference configuration.

Constitutive modeling

In their most general form, constitutive models define \(\bm \sigma\) in terms of state variables. In the model taken into consideration in the present mini-app, the state variables are constituted by the vector displacement field \(\bm u\), and its gradient \(\nabla \bm u\). We begin by defining the symmetric (small/infintesimal) strain tensor as

(26)\[\bm{\epsilon} = \dfrac{1}{2}\left(\nabla \bm{u} + \nabla \bm{u}^T \right).\]

This constitutive model \(\bm \sigma(\bm \epsilon)\) is a linear tensor-valued function of a tensor-valued input, but we will consider the more general nonlinear case in other models below. In these cases, an arbitrary choice of such a function will generally not be invariant under orthogonal transformations and thus will not admissible as a physical model must not depend on the coordinate system chosen to express it. In particular, given an orthogonal transformation \(Q\), we desire

(27)\[Q \bm \sigma(\bm \epsilon) Q^T = \bm \sigma(Q \bm \epsilon Q^T),\]

which means that we can change our reference frame before or after computing \(\bm \sigma\), and get the same result either way. Constitutive relations in which \(\bm \sigma\) is uniquely determined by \(\bm \epsilon\) while satisfying the invariance property (27) are known as Cauchy elastic materials. Here, we define a strain energy density functional \(\Phi(\bm \epsilon) \in \mathbb R\) and obtain the strain energy from its gradient,

(28)\[\bm \sigma(\bm \epsilon) = \frac{\partial \Phi}{\partial \bm \epsilon}.\]


The strain energy density functional cannot be an arbitrary function \(\Phi(\bm \epsilon)\); it can only depend on invariants, scalar-valued functions \(\gamma\) satisfying

\[\gamma(\bm \epsilon) = \gamma(Q \bm \epsilon Q^T) \]

for all orthogonal matrices \(Q\).

For the linear elasticity model, the strain energy density is given by

\[\bm{\Phi} = \frac{\lambda}{2} (\operatorname{trace} \bm{\epsilon})^2 + \mu \bm{\epsilon} : \bm{\epsilon} .\]

The constitutive law (stress-strain relationship) is therefore given by its gradient,

\[\bm\sigma = \lambda (\operatorname{trace} \bm\epsilon) \bm I_3 + 2 \mu \bm\epsilon, \]

where \(\bm I_3\) is the \(3 \times 3\) identity matrix, the colon represents a double contraction (over both indices of \(\bm \epsilon\)), and the Lamé parameters are given by

\[\begin{aligned} \lambda &= \frac{E \nu}{(1 + \nu)(1 - 2 \nu)} \\ \mu &= \frac{E}{2(1 + \nu)} \end{aligned}. \]

The constitutive law (stress-strain relationship) can also be written as

(29)\[\bm{\sigma} = \mathsf{C} \!:\! \bm{\epsilon}.\]

For notational convenience, we express the symmetric second order tensors \(\bm \sigma\) and \(\bm \epsilon\) as vectors of length 6 using the Voigt notation. Hence, the fourth order elasticity tensor \(\mathsf C\) (also known as elastic moduli tensor or material stiffness tensor) can be represented as

(30)\[\mathsf C = \begin{pmatrix} \lambda + 2\mu & \lambda & \lambda & & & \\ \lambda & \lambda + 2\mu & \lambda & & & \\ \lambda & \lambda & \lambda + 2\mu & & & \\ & & & \mu & & \\ & & & & \mu & \\ & & & & & \mu \end{pmatrix}.\]

Note that the incompressible limit \(\nu \to \frac 1 2\) causes \(\lambda \to \infty\), and thus \(\mathsf C\) becomes singular.

Hyperelasticity at Small Strain

The strong and weak forms given above, in (24) and (25), are valid for Neo-Hookean hyperelasticity at small strain. However, the strain energy density differs and is given by

\[\bm{\Phi} = \lambda (1 + \operatorname{trace} \bm{\epsilon}) (\log(1 + \operatorname{trace} \bm\epsilon) - 1) + \mu \bm{\epsilon} : \bm{\epsilon} .\]

As above, we have the corresponding constitutive law given by

(31)\[\bm{\sigma} = \lambda \log(1 + \operatorname{trace} \bm\epsilon) \bm{I}_3 + 2\mu \bm{\epsilon}\]

where \(\bm{\epsilon}\) is defined as in (26).

Newton linearization

Due to nonlinearity in the constitutive law, we require a Newton linearization of (31). To derive the Newton linearization, we begin by expressing the derivative,

\[\diff \bm{\sigma} = \dfrac{\partial \bm{\sigma}}{\partial \bm{\epsilon}} \colon \diff \bm{\epsilon}\]


\[\diff \bm{\epsilon} = \dfrac{1}{2}\left( \nabla \diff \bm{u} + \nabla \diff \bm{u}^T \right)\]


\[\diff \nabla \bm{u} = \nabla \diff \bm{u} .\]


(32)\[\diff \bm{\sigma} = \bar{\lambda} \cdot \operatorname{trace} \diff \bm{\epsilon} \cdot \bm{I}_3 + 2\mu \diff \bm{\epsilon}\]

where we have introduced the symbol

\[\bar{\lambda} = \dfrac{\lambda}{1 + \epsilon_v }\]

where volumetric strain is given by \(\epsilon_v = \sum_i \epsilon_{ii}\).

Equation (32) can be written in Voigt matrix notation as follows:

(33)\[\begin{pmatrix} \diff \sigma_{11} \\ \diff \sigma_{22} \\ \diff \sigma_{33} \\ \diff \sigma_{23} \\ \diff \sigma_{13} \\ \diff \sigma_{12} \end{pmatrix} = \begin{pmatrix} 2 \mu +\bar{\lambda} & \bar{\lambda} & \bar{\lambda} & & & \\ \bar{\lambda} & 2 \mu +\bar{\lambda} & \bar{\lambda} & & & \\ \bar{\lambda} & \bar{\lambda} & 2 \mu +\bar{\lambda} & & & \\ & & & \mu & & \\ & & & & \mu & \\ & & & & & \mu \\ \end{pmatrix} \begin{pmatrix} \diff \epsilon_{11} \\ \diff \epsilon_{22} \\ \diff \epsilon_{33} \\ 2 \diff \epsilon_{23} \\ 2 \diff \epsilon_{13} \\ 2 \diff \epsilon_{12} \end{pmatrix}.\]

Hyperelasticity at Finite Strain

In the total Lagrangian approach for the Neo-Hookean hyperelasticity problem, the discrete equations are formulated with respect to the reference configuration. In this formulation, we solve for displacement \(\bm u(\bm X)\) in the reference frame \(\bm X\). The notation for elasticity at finite strain is inspired by [Hol00] to distinguish between the current and reference configurations. As explained in the Common notation section, we denote by capital letters the reference frame and by small letters the current one.

The strong form of the static balance of linear-momentum at finite strain (total Lagrangian) is given by:

(34)\[- \nabla_X \cdot \bm{P} - \rho_0 \bm{g} = \bm{0}\]

where the \(_X\) in \(\nabla_X\) indicates that the gradient is calculated with respect to the reference configuration in the finite strain regime. \(\bm{P}\) and \(\bm{g}\) are the first Piola-Kirchhoff stress tensor and the prescribed forcing function, respectively. \(\rho_0\) is known as the reference mass density. The tensor \(\bm P\) is not symmetric, living in the current configuration on the left and the reference configuration on the right.

\(\bm{P}\) can be decomposed as

(35)\[\bm{P} = \bm{F} \, \bm{S},\]

where \(\bm S\) is the second Piola-Kirchhoff stress tensor, a symmetric tensor defined entirely in the reference configuration, and \(\bm{F} = \bm I_3 + \nabla_X \bm u\) is the deformation gradient. Different constitutive models can define \(\bm S\).

Constitutive modeling

For the constitutive modeling of hyperelasticity at finite strain, we begin by defining two symmetric tensors in the reference configuration, the right Cauchy-Green tensor

\[\bm C = \bm F^T \bm F \]

and the Green-Lagrange strain tensor

(36)\[\bm E = \frac 1 2 (\bm C - \bm I_3) = \frac 1 2 \Big( \nabla_X \bm u + (\nabla_X \bm u)^T + (\nabla_X \bm u)^T \nabla_X \bm u \Big),\]

the latter of which converges to the linear strain tensor \(\bm \epsilon\) in the small-deformation limit. The constitutive models considered, appropriate for large deformations, express \(\bm S\) as a function of \(\bm E\), similar to the linear case, shown in equation (29), which expresses the relationship between \(\bm\sigma\) and \(\bm\epsilon\).

Recall that the strain energy density functional can only depend upon invariants. We will assume without loss of generality that \(\bm E\) is diagonal and take its set of eigenvalues as the invariants. It is clear that there can be only three invariants, and there are many alternate choices, such as \(\operatorname{trace}(\bm E), \operatorname{trace}(\bm E^2), \lvert \bm E \rvert\), and combinations thereof. It is common in the literature for invariants to be taken from \(\bm C = \bm I_3 + 2 \bm E\) instead of \(\bm E\).

For example, if we take the compressible Neo-Hookean model,

(37)\[\begin{aligned} \Phi(\bm E) &= \frac{\lambda}{2}(\log J)^2 + \frac \mu 2 (\operatorname{trace} \bm C - 3) - \mu \log J \\ &= \frac{\lambda}{2}(\log J)^2 + \mu \operatorname{trace} \bm E - \mu \log J, \end{aligned}\]

where \(J = \lvert \bm F \rvert = \sqrt{\lvert \bm C \rvert}\) is the determinant of deformation (i.e., volume change) and \(\lambda\) and \(\mu\) are the Lamé parameters in the infinitesimal strain limit.

To evaluate (28), we make use of

\[\frac{\partial J}{\partial \bm E} = \frac{\partial \sqrt{\lvert \bm C \rvert}}{\partial \bm E} = \lvert \bm C \rvert^{-1/2} \lvert \bm C \rvert \bm C^{-1} = J \bm C^{-1}, \]

where the factor of \(\frac 1 2\) has been absorbed due to \(\bm C = \bm I_3 + 2 \bm E.\) Carrying through the differentiation (28) for the model (37), we arrive at

(38)\[\bm S = \lambda \log J \bm C^{-1} + \mu (\bm I_3 - \bm C^{-1}).\]


An equivalent form of (38) is

\[\bm S = \lambda \log J \bm C^{-1} + 2 \mu \bm C^{-1} \bm E, \]

which is more numerically stable for small \(\bm E\), and thus preferred for computation. Note that the product \(\bm C^{-1} \bm E\) is also symmetric, and that \(\bm E\) should be computed using (36).

Similarly, it is preferable to compute \(\log J\) using log1p, especially in case of nearly incompressible materials. To sketch this idea, suppose we have the \(2\times 2\) symmetric matrix \(C = \left( \begin{smallmatrix} 1 + e_{00} & e_{01} \\ e_{01} & 1 + e_{11} \end{smallmatrix} \right)\). Then we compute

\[\log \sqrt{\lvert C \rvert} = \frac 1 2 \mathtt{log1p}(e_{00} + e_{11} + e_{00} e_{11} - e_{01}^2). \]

which gives accurate results even in the limit when the entries \(e_{ij}\) are very small. For example, if \(e_{ij} \sim 10^{-8}\), then naive computation of \(\bm I_3 - \bm C^{-1}\) and \(\log J\) will have a relative accuracy of order \(10^{-8}\) in double precision and no correct digits in single precision. When using the stable choices above, these quantities retain full \(\varepsilon_{\text{machine}}\) relative accuracy.


One can linearize (38) around \(\bm E = 0\), for which \(\bm C = \bm I_3 + 2 \bm E \to \bm I_3\) and \(J \to 1 + \operatorname{trace} \bm E\), therefore (38) reduces to

(39)\[\bm S = \lambda (\operatorname{trace} \bm E) \bm I_3 + 2 \mu \bm E,\]

which is the St. Venant-Kirchoff model (constitutive linearization without geometric linearization; see (23)).

This model can be used for geometrically nonlinear mechanics (e.g., snap-through of thin structures), but is inappropriate for large strain.

Alternatively, one can drop geometric nonlinearities, \(\bm E \to \bm \epsilon\) and \(\bm C \to \bm I_3\), while retaining the nonlinear dependence on \(J \to 1 + \operatorname{trace} \bm \epsilon\), thereby yielding (31) (see (23)).

Weak form

We multiply (34) by a test function \(\bm v\) and integrate by parts to obtain the weak form for finite-strain hyperelasticity: find \(\bm u \in \mathcal V \subset H^1(\Omega_0)\) such that

(40)\[ \int_{\Omega_0}{\nabla_X \bm{v} \colon \bm{P}} \, dV - \int_{\Omega_0}{\bm{v} \cdot \rho_0 \bm{g}} \, dV - \int_{\partial \Omega_0}{\bm{v} \cdot (\bm{P} \cdot \hat{\bm{N}})} \, dS = 0, \quad \forall \bm v \in \mathcal V,\]

where \(\bm{P} \cdot \hat{\bm{N}}|_{\partial\Omega}\) is replaced by any prescribed force/traction boundary condition written in terms of the reference configuration. This equation contains material/constitutive nonlinearities in defining \(\bm S(\bm E)\), as well as geometric nonlinearities through \(\bm P = \bm F\, \bm S\), \(\bm E(\bm F)\), and the body force \(\bm g\), which must be pulled back from the current configuration to the reference configuration. Discretization of (40) produces a finite-dimensional system of nonlinear algebraic equations, which we solve using Newton-Raphson methods. One attractive feature of Galerkin discretization is that we can arrive at the same linear system by discretizing the Newton linearization of the continuous form; that is, discretization and differentiation (Newton linearization) commute.

Newton linearization

To derive a Newton linearization of (40), we begin by expressing the derivative of (35) in incremental form,

(41)\[\diff \bm P = \frac{\partial \bm P}{\partial \bm F} \!:\! \diff \bm F = \diff \bm F\, \bm S + \bm F \underbrace{\frac{\partial \bm S}{\partial \bm E} \!:\! \diff \bm E}_{\diff \bm S}\]


\[\diff \bm E = \frac{\partial \bm E}{\partial \bm F} \!:\! \diff \bm F = \frac 1 2 \Big( \diff \bm F^T \bm F + \bm F^T \diff \bm F \Big). \]

The quantity \({\partial \bm S} / {\partial \bm E}\) is known as the incremental elasticity tensor, and is analogous to the linear elasticity tensor \(\mathsf C\) of (30). We now evaluate \(\diff \bm S\) for the Neo-Hookean model (38),

(42)\[\diff\bm S = \frac{\partial \bm S}{\partial \bm E} \!:\! \diff \bm E = \lambda (\bm C^{-1} \!:\! \diff\bm E) \bm C^{-1} + 2 (\mu - \lambda \log J) \bm C^{-1} \diff\bm E \, \bm C^{-1},\]

where we have used

\[\diff \bm C^{-1} = \frac{\partial \bm C^{-1}}{\partial \bm E} \!:\! \diff\bm E = -2 \bm C^{-1} \diff \bm E \, \bm C^{-1} . \]


In the small-strain limit, \(\bm C \to \bm I_3\) and \(\log J \to 0\), thereby reducing (42) to the St. Venant-Kirchoff model (39).


Some cancellation is possible (at the expense of symmetry) if we substitute (42) into (41),

(43)\[\begin{aligned} \diff \bm P &= \diff \bm F\, \bm S + \lambda (\bm C^{-1} : \diff \bm E) \bm F^{-T} + 2(\mu - \lambda \log J) \bm F^{-T} \diff\bm E \, \bm C^{-1} \\ &= \diff \bm F\, \bm S + \lambda (\bm F^{-T} : \diff \bm F) \bm F^{-T} + (\mu - \lambda \log J) \bm F^{-T} (\bm F^T \diff \bm F + \diff \bm F^T \bm F) \bm C^{-1} \\ &= \diff \bm F\, \bm S + \lambda (\bm F^{-T} : \diff \bm F) \bm F^{-T} + (\mu - \lambda \log J) \Big( \diff \bm F\, \bm C^{-1} + \bm F^{-T} \diff \bm F^T \bm F^{-T} \Big), \end{aligned}\]

where we have exploited \(\bm F \bm C^{-1} = \bm F^{-T}\) and

\[\begin{aligned} \bm C^{-1} \!:\! \diff \bm E = \bm C_{IJ}^{-1} \diff \bm E_{IJ} &= \frac 1 2 \bm F_{Ik}^{-1} \bm F_{Jk}^{-1} (\bm F_{\ell I} \diff \bm F_{\ell J} + \diff \bm F_{\ell I} \bm F_{\ell J}) \\ &= \frac 1 2 \Big( \delta_{\ell k} \bm F_{Jk}^{-1} \diff \bm F_{\ell J} + \delta_{\ell k} \bm F_{Ik}^{-1} \diff \bm F_{\ell I} \Big) \\ &= \bm F_{Ik}^{-1} \diff \bm F_{kI} = \bm F^{-T} \!:\! \diff \bm F. \end{aligned} \]

We prefer to compute with (42) because (43) is more expensive, requiring access to (non-symmetric) \(\bm F^{-1}\) in addition to (symmetric) \(\bm C^{-1} = \bm F^{-1} \bm F^{-T}\), having fewer symmetries to exploit in contractions, and being less numerically stable.

It is sometimes useful to express (42) in index notation,

(44)\[\begin{aligned} \diff\bm S_{IJ} &= \frac{\partial \bm S_{IJ}}{\partial \bm E_{KL}} \diff \bm E_{KL} \\ &= \lambda (\bm C^{-1}_{KL} \diff\bm E_{KL}) \bm C^{-1}_{IJ} + 2 (\mu - \lambda \log J) \bm C^{-1}_{IK} \diff\bm E_{KL} \bm C^{-1}_{LJ} \\ &= \underbrace{\Big( \lambda \bm C^{-1}_{IJ} \bm C^{-1}_{KL} + 2 (\mu - \lambda \log J) \bm C^{-1}_{IK} \bm C^{-1}_{JL} \Big)}_{\mathsf C_{IJKL}} \diff \bm E_{KL} \,, \end{aligned}\]

where we have identified the effective elasticity tensor \(\mathsf C = \mathsf C_{IJKL}\). It is generally not desirable to store \(\mathsf C\), but rather to use the earlier expressions so that only \(3\times 3\) tensors (most of which are symmetric) must be manipulated. That is, given the linearization point \(\bm F\) and solution increment \(\diff \bm F = \nabla_X (\diff \bm u)\) (which we are solving for in the Newton step), we compute \(\diff \bm P\) via

  1. recover \(\bm C^{-1}\) and \(\log J\) (either stored at quadrature points or recomputed),

  2. proceed with \(3\times 3\) matrix products as in (42) or the second line of (44) to compute \(\diff \bm S\) while avoiding computation or storage of higher order tensors, and

  3. conclude by (41), where \(\bm S\) is either stored or recomputed from its definition exactly as in the nonlinear residual evaluation.


The decision of whether to recompute or store functions of the current state \(\bm F\) depends on a roofline analysis [WWP09,Brown10] of the computation and the cost of the constitutive model. For low-order elements where flops tend to be in surplus relative to memory bandwidth, recomputation is likely to be preferable, where as the opposite may be true for high-order elements. Similarly, analysis with a simple constitutive model may see better performance while storing little or nothing while an expensive model such as Arruda-Boyce [AB93], which contains many special functions, may be faster when using more storage to avoid recomputation. In the case where complete linearization is preferred, note the symmetry \(\mathsf C_{IJKL} = \mathsf C_{KLIJ}\) evident in (44), thus \(\mathsf C\) can be stored as a symmetric \(6\times 6\) matrix, which has 21 unique entries. Along with 6 entries for \(\bm S\), this totals 27 entries of overhead compared to computing everything from \(\bm F\). This compares with 13 entries of overhead for direct storage of \(\{ \bm S, \bm C^{-1}, \log J \}\), which is sufficient for the Neo-Hookean model to avoid all but matrix products.