(matrix-free-mixed-fields)=
# Generic weak form (mixed fields)

Similar to {ref}`matrix-free-single-field` for single field, we present a generic weak formulation for a general Dirichlet problem for general mixed problem $\bm{R}(\bm{u}, p) = [\bm{R}_u(\bm{u}, p) ~~ \bm{R}_p(\bm{u}, p)]^T = \bm 0$: Find $(\bm u, p) \in \mathcal{V} \times \mathcal{Q} $, such that

$$
  \begin{aligned}
    \langle\bm{v},  \bm{R}_u(\bm{u}, p) \rangle &= \int_{\Omega_0}{\bm{v} \cdot \bm{f}_0 + \nabla \bm{v} \tcolon \bm{f}_1} \, dV = 0, \quad \forall \bm{v} \in \mathcal{V}, \\
		\langle q,  \bm{R}_p(\bm{u}, p) \rangle &= \int_{\Omega_0}{q \cdot \bm{g}_0 + \nabla q \tcolon \bm{g}_1 } \, dV = 0, \quad \forall q \in \mathcal{Q},
  \end{aligned}
$$ (mixed-residual-matrix-free)

where the operators $\bm{f}_0, \bm{f}_1, \bm{g}_0, \bm{g}_1$ contain all possible sources in the problem i.e., $\bm{u}, \nabla \bm{u}, p, \nabla p$. It should be noted that the gradient in the {eq}`mixed-residual-matrix-free` depends on the configuration system and it could be with respect to initial configuration $\bm{X}$ i.e.,  $\nabla_X \bm{v}$ (Lagrangian approach) or current configuration $\bm{x}$ i.e.,  $\nabla_x \bm{v}$ (Eulerian approach).

In order to solve {eq}`mixed-residual-matrix-free` with Newton-Krylov iterative solvers we need its Jacobian form: Find $(\diff \bm{u}, \diff p) \in \mathcal{V} \times \mathcal{Q} $ such that

$$
  \begin{aligned}
    \langle\bm{v},  \bm{J}_u(\bm{u}, p) \diff \bm{U} \rangle &= \int_{\Omega_0}{\bm{v} \cdot \diff \bm{f}_0 + \nabla \bm{v} \tcolon \diff \bm{f}_1} \, dV, \\
		\langle q,  \bm{J}_p(\bm{u}, p) \diff \bm{U} \rangle &= \int_{\Omega_0}{q \cdot \diff \bm{g}_0 + \nabla q \tcolon \diff \bm{g}_1} \, dV.
  \end{aligned}
$$ (mixed-jacobian-matrix-free)

where $\diff \bm{U} = [\diff \bm{u} ~~ \diff p]^T$ and the linearization of operators $\bm{f}_0, \bm{f}_1, \bm{g}_0, \bm{g}_1$ are

$$
  \begin{aligned}
    \diff \bm{f}_i &= \frac{\partial \bm{f}_i}{\partial \bm{u}} \diff \bm{u} + \frac{\partial \bm{f}_i}{\partial \nabla \bm{u}} \nabla \diff \bm{u} + \frac{\partial \bm{f}_i}{\partial p} \diff p + \frac{\partial \bm{f}_i}{\partial \nabla p} \nabla \diff p, \quad i=0,1 \\
		\diff \bm{g}_i &= \frac{\partial \bm{g}_i}{\partial \bm{u}} \diff \bm{u} + \frac{\partial \bm{g}_i}{\partial \nabla \bm{u}} \nabla \diff \bm{u} + \frac{\partial \bm{g}_i}{\partial p} \diff p + \frac{\partial \bm{g}_i}{\partial \nabla p} \nabla \diff p, \quad i=0,1
  \end{aligned}
$$

Compare with governing equations derived in pervious sections for linear and large deformation, it is easy to verify that

* For mixed linear elasticity described in {eq}`mixed-linear-weak-form` we have

$$
  \begin{aligned}
    \bm f_0 &= -\rho \bm g, \quad \bm f_1 = \bm{\sigma}(\bm u, p), \quad \bm g_0 =-\trace \bm \varepsilon -  \frac{p}{\bulk-\bulk_p}, \quad \bm g_1 = \bm 0, \\
			\diff \bm f_1 &= \bm f_1(\diff \bm{u}, \diff p), ~ ~ \diff \bm g_0 = \bm g_0(\diff \bm{u}, \diff p)
  \end{aligned}
$$

* For mixed hyperelastic in initial configuration described in {eq}`mixed-hyperelastic-weak-form-initial` and {eq}`mixed-jacobian-weak-form-initial` we have

$$
  \begin{aligned}
		\bm f_0 &= -\rho_0 \bm g, \quad \bm f_1 = \bm F \bm S, \quad \bm g_0 = L J, \quad \bm g_1 = \bm 0, \\
		\diff \bm f_1 &= \diff \bm F \bm S + \bm F \diff \bm S, ~ ~ \diff \bm g_0 = \diff L J + L \diff J
  \end{aligned}
$$

* For mixed hyperelastic in current configuration derived in {eq}`mixed-jacobian-weak-form-current`

$$
  \begin{aligned}
		\bm f_0 &= -\rho_0 \bm g, \quad \bm f_1 = \bm \tau, \quad \bm g_0 = L J, \quad \bm g_1 = \bm 0, \\
		\diff \bm f_1 &= \diff \bm{\tau} - \bm{\tau} \left( \nabla_x \diff \bm{u} \right)^T, ~ ~ \diff \bm g_0 = \diff L J + L \diff J
  \end{aligned}
$$

* For mixed hyperelastic with variational approach in initial configuration given in {eq}`mixed-hyperelastic-weak-form-initial-PL` and {eq}`mixed-jacobian-weak-form-initial-PL`

$$
  \begin{aligned}
		\bm f_0 &= -\rho_0 \bm g, \quad \bm f_1 = \bm F \bm S, \quad \bm g_0 = -{U}(J) - \frac{p}{\bulk - \bulk_p}, \quad \bm g_1 = \bm 0, \\
		\diff \bm f_1 &= \diff \bm F \bm S + \bm F \diff \bm S, ~ ~ \diff \bm g_0 = -\diff {U}(J) - \frac{\diff p}{\bulk - \bulk_p}
  \end{aligned}
$$

We solve the above mixed equations using static and quasistatic solver in Ratel as explained in elasticity [static and quasistatic](static-quasistatic) section. For mixed elastodynamic, we only need to add the same terms (see {eq}`elastodynamics-residual` and {eq}`elastodynamics-jacobian`) to the residual and jacobian of the first equation of mixed matrix-free formulation described in {eq}`mixed-residual-matrix-free` and {eq}`mixed-jacobian-matrix-free`, respectively.

:::{note}
For solving incompressible elasticity we need to use different polynomial orders for displacement and pressure fields. These polynomial orders cannot be chosen arbitrarily and must fulfill the Ladyzhenskaya–Babuška–Brezzi (LBB) or inf–sup condition to ensure stability and optimal convergence {cite}`BrezziFortin2012, ChapelleBathe1993`. If we rewrite the linear weak form {eq}`mixed-linear-weak-form` as

$$
  \begin{aligned}
    a(\bm v, \bm{u}) + b(\bm v, p) &= f(\bm v),  \quad \forall \bm v \in \mathcal{V}^0 \\
    b(q, \bm{u}) + c(q, p) &= 0, \quad \forall q \in \mathcal{Q}
  \end{aligned}
$$

then,the inf-sup constant is computable through a set of eigenvalue problems {cite}`FEniCS2012, arnold2009stability`: Find $\lambda \in \mathbb{R}$, $0 \neq (\bm u_h, p_h) \in \mathcal{V}_h \times \mathcal{Q}_h $, such that for all $(\bm v, q) \in \mathcal{V}_h^0 \times \mathcal{Q}_h$

$$
  \langle \bm v, \bm u_h \rangle_{\mathcal{V}_h} + b(\bm v, p_h) + b(q, \bm u_h) = -\lambda \langle q, p_h \rangle_{\mathcal{Q}_h},
$$ (infsup-eigenvalue-problem)

where $b(\bm v, p_h)$ is a bilinear form and $\langle \bm v, \bm u_h \rangle_{\mathcal{V}_h}$ is the inner product of the associated space $\mathcal{V}_h$. The Brezzi stability conditions state that there exists a unique solution for mixed problem if and only if inf-sup constant, $\beta_h = \sqrt{\lambda_{\text{min}}}$, is bounded below for all mesh size $h > 0$.

```{figure} /doc/img/mixed-element.svg
---
align: center
name: mixed-element-fig
---
2D mixed element examples with discontinuous and continuous (last column) pressure.
```

{ref}`mixed-element-fig` shows some mixed elements where the pressure DoF can be continuous or discontinuous. We performed the inf-sup analysis by solving the eigen problem {eq}`infsup-eigenvalue-problem` for various mixed elements on square $[0,1]^2$ and stretched meshes $[0,1] \times [0, 0.1]$ and displayed on {ref}`inf-sup-fig`. All mixed elements with continuous pressure are inf-sup stable on element with aspect ratio 1 and 10. However, for the discontinuous pressure, the inf-sup constant for widely-used $Q_1P_0$ and $Q_nQ_{n-1}$ elements is decreasing under mesh refinement which make them unstable. Moreover, $Q_2P_1$ element which is stable on square mesh, becomes unstable on the stretched meshes. In fact, the inf-sup constant of $Q_2P_1$ is decreasing by factor $1/\sqrt{a}$ where $a$ is element aspect ratio. Thus, for stretched structure one should avoid using this element. Results exhibited below is created using the open source Julia programming language and is reproducible as given in this [demo code](https://github.com/rezgarshakeri/fe-demo/blob/main/MixedElasticity2D.ipynb).

```{figure} /doc/img/inf-sup.png
---
align: center
name: inf-sup-fig
---
Inf-sup constant for various mixed-elements with continuous (top row) and discontinuous (bottom row) pressure on square and stretched meshes.
```

It should be noted element that satisfies the inf-sup condition may exhibit instability in large deformation but one should at least use the element that satisfies LBB condition for linear elasticity. To solve a system with equal-order interpolation, the mixed model must be used in conjunction with stabilization techniques like Streamline Upwind Petrov–Galerkin (SUPG) and Galerkin Least Squares (GLS).

(mixed_element_options)=
:::{list-table} Mixed element Options
:header-rows: 1
:widths: 40 50 10
* - Option
  - Description
  - Default Value

* - `-orders [int list]`
  - Polynomial order of solution basis functions (displacement, pressure)
  -

* - `-[prefix]_petscdualspace_lagrange_continuity [0 | 1]`
  - Enable/disable continuity of basis functions.
  - `1` (continuous element)

* - `-[prefix]_petscspace_poly_tensor [0 | 1]`
  - Enable/disable tensor product of basis functions
  - `1` (tensor basis)
:::

:::