(perturbed-lagrange-multiplier-method)=
# Perturbed Lagrange-multiplier method
## Functional Form
In the previous sections we introduced the strong form of general mixed formulation for small and finite strain and by considering test functions $(\bm v, q)$ as derived in {eq}`mixed-linear-weak-form` and {eq}`mixed-hyperelastic-weak-form-initial`. Alternatively, we can derive the mixed $\bm u- p$ weak formulation based on minimization of two fields functional $\Pi(\bm u, p)$, which it is known as the perturbed Lagrangian approach. For the mixed linear case, we can write

$$
  \begin{aligned}
    \Pi (\bm u, p) &= \int_{\Omega} \left[ \mu \, \bm \varepsilon_{\text{dev}} \tcolon \bm \varepsilon_{\text{dev}} - p\,  \trace \bm \varepsilon + \frac{\bulk_p}{2}  \left(\trace \bm \varepsilon \right)^2 - \frac{p^2}{2 (\bulk-\bulk_p)}  \right] \, dv - \Pi_{\text{ext}} (\bm u) \\
    \Pi_{\text{ext}} (\bm u) &= \int_{\Omega}{\bm{u} \cdot \rho \bm{g}} \, dv + \int_{\partial \Omega}{\bm{u} \cdot \bar{\bm{t}}} \, da
  \end{aligned}
$$ (mixed-linear-functional)

and by invoking the stationarity of $\Pi$ with respect to $\bm u$ and $p$, we obtain

$$
  \begin{aligned}
      \int_{\Omega} \nabla \delta \bm u \tcolon \left[ 2 \mu \, \bm \varepsilon_{\text{dev}} + \left(\bulk_p \trace \bm \varepsilon - p \right) \bm{I} \right] \, dv - \bm L_{\text{ext}}(\delta \bm u) &=0, \\
      \int_{\Omega} \delta p \left( -\trace \bm \varepsilon -  \frac{p}{\bulk-\bulk_p} \right) \, dv &= 0,
  \end{aligned}
$$ (mixed-linear-weak-form-PL)

where $\bm L_{\text{ext}}(\delta \bm u) = \int_{\Omega}{\delta \bm{u} \cdot \rho \bm{g}} \, dv + \int_{\partial \Omega}{\delta \bm{u} \cdot \bar{\bm{t}}} \, da$, and $\delta \bm u$ and $\delta p$ are virtual displacement and pressure and can be seen as the test functions $\bm v, q$, i.e., $\delta \bm u = \bm v, \, \delta p = q$. It is clear that the weak form {eq}`mixed-linear-weak-form-PL` agrees with what we obtained in {eq}`mixed-linear-weak-form`. However, the hyperelastic weak form {eq}`mixed-hyperelastic-weak-form-initial`, can not be derived by minimizing any functional and its linearization is not symmetric.

To write the mixed functional $\Pi(\bm u, p)$ for hyperelastic, we need to consider the strain energy of the form

$$
  \psi \left(\bm{C} \right) = \psi_{\text{vol}}(J) + \psi_{\text{iso}}(\bar{\bm{C}}) = \frac{\bulk}{2} \left(U(J)\right)^2 + \psi_{\text{iso}}(\bar{\bm{C}})
$$ (mixed-strain-energy-PL)

we can write a two fields energy functional as

$$
  \begin{aligned}
  \Pi (\bm u, p) &= \int_{\Omega_0} \left[ \psi_{\text{iso}}(\bar{\bm{C}}) - p \, U(J) - \frac{1}{2}\frac{p^2}{\bulk - \bulk_p} + \frac{\bulk_p}{2} U^2\right] \, dV - \Pi_{\text{ext}} (\bm u) \\
  \Pi_{\text{ext}} (\bm u) &= \int_{\Omega_0}{\bm{u} \cdot \rho_0 \bm{g}} \, dV + \int_{\partial \Omega_0}{\bm{u} \cdot \bar{\bm{t}}} \, dA
  \end{aligned}
$$ (mixed-energy-functional)

Finding the stationary conditions with respect to $\bm{u}$ and $p$ by taking Gateaux derivative gives the weak form

$$
  \begin{aligned}
    \int_{\Omega_0} \bm F \underbrace{\left( \bm{S}_{\text{iso}} + (\bulk_p U - p) \, J \, U' \, \bm{C}^{-1} \right)}_{\bm S} \tcolon \nabla_X \bm v \, dV &= L_{\text{ext}} (\bm v) \\
    \int_{\Omega_0} \left(- U(J) - \frac{p}{\bulk - \bulk_p} \right) q \, dV &= 0
  \end{aligned}
$$ (mixed-hyperelastic-weak-form-initial-PL)

where, $L_{\text{ext}} (\bm v) = \int_{\Omega_0}{\bm{v} \cdot \rho_0 \bm{g}} \, dV + \int_{\partial \Omega_0}{\bm{v} \cdot \bar{\bm{t}}} \, dA $ and we have used

$$
  \begin{aligned}
    \frac{\partial \psi_{\text{iso}}}{\partial \bm u} \cdot {\delta \bm u} &= \frac{\partial \psi_{\text{iso}}}{\partial \bm E} \tcolon \frac{\partial \bm E}{ \partial \bm u}{\delta \bm u} = \bm{S}_{\text{iso}} \tcolon {\delta \bm E} = \bm{S}_{\text{iso}} \tcolon \text{sym} \left(\bm F^T \delta \bm F \right) \\
    \frac{\partial U}{\partial \bm u} \cdot {\delta \bm u} &= \frac{\partial U}{\partial J} \frac{\partial J}{\partial \bm u} {\delta \bm u} = U' {\delta J} = J \, U' \, \bm{C}^{-1} \tcolon \tcolon {\delta \bm E} = J \, U' \, \bm{C}^{-1} \tcolon \text{sym} \left(\bm F^T \delta \bm F \right)
  \end{aligned}
$$

where $\delta \bm{F} = \nabla_X \delta \bm u = \nabla_X \bm v$. The Jacobian for problem {eq}`mixed-hyperelastic-weak-form-initial-PL` can be written as

$$
  \begin{aligned}
    \int_{\Omega_0} \nabla_X \bm{v} \tcolon \left(\bm F \diff \bm{S} + \diff \bm{F} \bm{S} \right) dV &= -R_u^{PL}, \\
    \int_{\Omega_0} q \left( -\diff U - \frac{\diff p}{\bulk - \bulk_p} \right) dV &= -R_p^{PL},
  \end{aligned}
$$ (mixed-jacobian-weak-form-initial-PL)

where

$$
  \begin{aligned}
    \diff \bm{S} &= \diff \bm{S}_{\text{iso}} + \diff \bm{S}_{\text{vol}}^u  + \diff \bm{S}_{\text{vol}}^p, \\
    \diff \bm{S}_{\text{vol}}^u &= \left(\bulk_p U' \diff J \right) J U' \bm{C}^{-1} + \left(\bulk_p U - p \right) \left(\diff J \, U' \bm{C}^{-1} + J\, U^{''} \diff J \bm{C}^{-1} + J U' \diff \bm{C}^{-1} \right), \\
    &= \left[\bulk_p \left( J U' \right)^2 +  \left(\bulk_p U - p \right) \left( J U' + J^2  U^{''} \right) \right]  \left( \bm{C}^{-1} \tcolon \diff \bm E \right) \bm{C}^{-1} \\
    &+ \left(\bulk_p U - p \right) J U' \, \diff \bm{C}^{-1}\\
    \diff \bm{S}_{\text{vol}}^p &= - dp J U' \, \bm{C}^{-1} \\
    \diff U &= U' \diff J = J \, U' \, \bm{C}^{-1} \tcolon \diff \bm E,
  \end{aligned}
$$ (mixed-stationary-point-linearization)

where $\diff \bm{S}_{\text{iso}}$ is derived for neo-Hookean, Mooney-Rivlin and Ogden in {eq}`dS-iso-neo-hookean`, {eq}`dS-iso-mooney-rivlin` and {eq}`dS-iso-ogden`. To complete the derivation we only need $U(J)$ function with condition

$$
  U(J) = 0 \quad \text{if and only if} \quad J = 1.
$$

For the volumetric strain energy function of the form given in {eq}`mixed-neo-hookean-energy`, if we choose $V(J) = \frac{1}{4} \left(J^2 - 1 - 2 \log J \right)$, from $\bulk V(J) = \frac{\bulk}{2} U^2(J)$, we will have

$$
  U(J) = \pm \sqrt{2 V} = \frac{\sign(J-1)}{\sqrt{2}} ( \underbrace{J^2 - 1 - 2 \log J }_{A})^{1/2}
$$

where the derivatives are

$$
  \begin{aligned}
    U' &= \frac{\partial U}{\partial J} = \sign(J-1) \frac{J^2 - 1}{J \sqrt{2}} A^{-1/2}, \\
    U{''} &= \frac{\partial^2 U}{\partial J^2} = \sign(J-1) \frac{1}{J^2 \sqrt{2}} \left((J^2 + 1) A^{-1/2} - (J^2 - 1)^2 A^{-3/2} \right).
  \end{aligned}
$$

## Command-line interface
To enable the perturbed Lagrange model, use the model option `-model elasticity-mixed-neo-hookean-PL-current` or `-model elasticity-mixed-neo-hookean-PL-initial` and set the material parameters listed in {ref}`perturbed-lagrange-options`. Any parameter without a default option is required.

(perturbed-lagrange-options)=

:::{list-table} Perturbed Lagrange-multiplier model options
:header-rows: 1
:widths: 30 55 15
* - Option
  - Description
  - Default value

* - `-model elasticity-mixed-neo-hookean-PL-current | elasticity-mixed-neo-hookean-PL-initial`
  - Required to enable the perturbed Lagrange model, `|` denotes an option for the user.
  -

* - `-E [real]`
  - [Young's modulus](https://en.wikipedia.org/wiki/Young%27s_modulus), $E > 0$
  -

* - `-nu [real]`
  - [Poisson's ratio](https://en.wikipedia.org/wiki/Poisson%27s_ratio), $\nu \leq 0.5$.
  -

* - `-nu_primal [real]`
  - Primal part of the Poisson's ratio for Jacobian operator, $-1 \leq \nu_p < \nu \leq 0.5$.
  This parameter may improve solver performance and should not affect the solution.
  - `-1.0`

* - `-nu_primal_pc [real]`
  - Primal part of the Poisson's ratio for preconditioner operator, $-1 \leq \nu_{pc} < \nu \leq 0.5$.
  This parameter may improve solver performance and should not affect the solution.
  - `-nu_primal` value

:::

An example using the pertrubed Lagrange model can be run via
```console
$ ./bin/ratel-quasistatic -options_file examples/ymls/ex01-static-elasticity-mixed-neo-hookean-PL-current-pcjacobi.yml
```