Mixed Finite Element Method

• Mixed Finite Element Method is a numerical technique that simultaneously discretizes multiple coupled fields in suitable function spaces to form a saddle-point system.
• It enables strong conservation, handles variational inequalities, and enforces weak constraints in applications such as elastoplasticity, porous media, and electromagnetics.
• Recent advancements, including higher-order and hp-adaptivity alongside robust error analysis, have significantly enhanced its accuracy and practical applicability.

The mixed finite element method (MFEM) is a broad class of finite element techniques in which two or more physically or mathematically coupled fields are discretized simultaneously, each in appropriately chosen function spaces, leading to a saddle-point system. MFEMs are particularly powerful for applications in mechanics, porous media, electromagnetics, and interface problems, providing strong conservation properties, the ability to handle variational inequalities or nonsmooth phenomena, and the flexibility to enforce constraints such as incompressibility, symmetry, or interface conditions in a weak (Lagrange multiplier) sense. Recent advances, including higher-order and $hp$-adaptivity, the use of nonconforming or immersed discretizations, and robust stability/error analysis, have greatly expanded the applicability and rigor of these methods.

1. Variational Formulations and Three-Field Mixed Methods

The foundation of mixed finite element methods lies in the simultaneous discretization of coupled field variables, often derived by recasting PDE systems or variational inequalities as saddle-point problems. A canonical example from elastoplasticity is the three-field formulation arising in small-strain models with kinematic hardening and von Mises-type dissipation:

Let $\Omega\subset\mathbb R^d$ be a bounded domain, $V=H^1(\Omega;\mathbb R^d)$ the displacement space with Dirichlet conditions, and $Q=L^2(\Omega; S_{d,0})$ the space of deviatoric plastic strains. The classical, nonsmooth plasticity functional can be dualized by a Lagrange multiplier $\lambda$ taking values in a convex set $$\Lambda=\{\mu\in Q: |\mu(x)|_F\leq \sigma_y \ a.e.\}$$, leading to the following mixed problem (Bammer et al., 17 Jan 2024):

• Find $(u,p,\lambda)\in V\times Q\times\Lambda$ such that
  • $$a((u,p), (v,q)) + \int_\Omega \lambda : q\, dx = \ell(v)$$ for all $(v,q)\in V\times Q$,
  • $\int_\Omega (\mu-\lambda) : p\, dx \leq 0$ for all $\mu\in \Lambda$,
  • where $a$ encodes elastic/plastic moduli, and possible hardening. This construction untangles the non-differentiability of the plastic dissipation and renders the problem amenable to robust discretization and analysis.

2. Conforming and Nonconforming $hp$-Finite Element Spaces

The MFEM requires careful selection of function spaces to discretize each field, ensuring both conformity in the relevant Sobolev norms and compatibility for the saddle-point structure. Higher-order and $hp$-adaptivity allow spectral convergence rates on smooth domains and robust accuracy near singularities.

Given a quadrilateral/hexahedral mesh $\mathcal{T}_h$ with diameters $h_T$, and elementwise polynomial degrees $p_T$, primary choices are (Bammer et al., 17 Jan 2024):

• Displacement: $$V_{hp} = \{v_h \in V: v_h\circ F_T \in \mathbb{P}_{p_T}(\hat T)^d\}$$ (conforming),
• Plastic strain: $$Q_{hp} = \{q_h \in Q: q_h\circ F_T \in \mathbb{P}_{p_T-1}(\hat T)^{d\times d}\}$$ (conforming),
• Multiplier: $$\Lambda_{hp}= \{\mu_h\in Q_{hp}: |\mu_h(x_{k,T})|_F\leq \sigma_y$$ for all Gauss points $x_{k,T}$ on $T\}$ (enforced at quadrature points).

Higher polynomial degree $p_T$ or local mesh refinement $h_T$ can be combined to achieve $hp$-adaptivity. Notably, enforcing set-valued constraints (such as yield or Tresca conditions) at quadrature nodes balances computational efficiency with rigorous approximation properties.

3. Discrete Saddle-Point Stability and Well-Posedness

Stability in MFEM discretizations is fundamentally dictated by an inf-sup (Babuška-Brezzi) condition, ensuring the absence of spurious solutions and uniform control over the discrete solution as $h, p$ vary. For the three-field elastoplasticity MFEM, the functional $b(\mu_h;q_h) = \int_{\Omega} \mu_h : q_h\,dx$ satisfies a discrete inf-sup with constant $1$: $$\|\mu_h\|_{L^2(\Omega)} = \sup_{(v_h, q_h)\ne 0} \frac{b(\mu_h;q_h)}{\|(v_h, q_h)\|},\ \forall \mu_h \in Q_{hp},$$ where $$\|(v, q)\|^2 = \|v\|_{H^1(\Omega)}^2 + \|q\|_{L^2(\Omega)}^2$$. This property, verified by choosing test fields directly, leads to well-posedness and explicit stability estimates: $$\|(u_h,p_h)\| \leq \frac{1}{\alpha}( \|f\|_{V^*} + c_{tr}\|g\|_{H^{-1/2}}),\ \|\lambda_h\|_{L^2} \leq \frac{c_a}{\alpha}( \|f\|_{V^*} + c_{tr}\|g\|_{H^{-1/2}}),$$ with mesh- and $p$-independent constants (Bammer et al., 17 Jan 2024).

4. A Priori Error Analysis and Optimal Convergence Rates

Error analysis in the MFEM context establishes strong convergence and rates depending on both mesh size $h$ and polynomial degree $p$:

• Without extra solution regularity, one guarantees strong convergence of all fields as $h/p\rightarrow 0$.
• For minimal-order elements ($p_T\equiv1$), and solution regularity $(u,p,\lambda)\in H^s\times H^t\times H^l$ with $s\geq1$, $t,l\geq0$, the error admits the sharp estimate:

$$\|u-u_h\|^2_{H^1} + \|p-p_h\|^2_{L^2} + \|\lambda-\lambda_h\|^2_{L^2} \lesssim h^{2\min(1,s-1,t,l)}( |u|_{s}^2 + |p|_{t}^2 + |\lambda|_{l}^2 )$$

yielding $O(h^2)$ in the fully smooth case ($s,t,l\geq 2$).

• For general $hp$-elements and solution smoothness $t,l>d/2$, $s\geq1$:

$$\|u-u_h\|^2_{H^1} + \|p-p_h\|^2_{L^2} + \|\lambda-\lambda_h\|^2_{L^2} \lesssim (h^{\min(p,2s-2,t,l)}\,p^{-\min(2s-2,t,l)})^2$$

This demonstrates that MFEMs achieve spectral ($hp$-) convergence on smooth domains, with the only minor loss coming from the constraint's quadrature enforcement (Bammer et al., 17 Jan 2024).

5. Numerical Experiments and Adaptive Strategies

Empirical studies confirm theoretical convergence rates and highlight the effectiveness of MFEMs under both uniform and adaptive refinement (Bammer et al., 17 Jan 2024):

• Uniform $h$-refinement with $p=1$ achieves EOC $\approx 0.45$ (corresponding to nearly $h^{1/2}$ in 1D, the minimal regularity scenario).
• Higher $p$ with uniform $h$ yields reduced EOCs in the presence of geometric singularities.
• Adaptive $h$-refinement recovers full optimal rates as the local mesh resolves boundary or internal layers.
• Full $hp$-adaptive refinement (elementwise refinement and $p$-enrichment) achieves EOCs $\approx 1.5$ or higher, exemplifying the power of MFEM in resolving singularities and boundary layers efficiently.

All numerical findings validate the saddle-point theory, stability, and a priori estimates developed for the higher-order MFEM framework, and substantiate its practical utility in demanding elastoplastic simulations.

6. Extensions and Significance in Broader Context

The mixed finite element methodology is extensively generalized for:

• Other variational inequalities and saddle-point systems (Contact, friction, Stokes, Maxwell),
• Nonlinear and degenerately constrained problems,
• Multi-physics applications (poroelasticity, magneto-mechanics, interface problems),
• Coupled multipoint or hybridized formulations for cell-centered and domain-decomposition solvers.

The explicit incorporation of additional fields (multiplier, plastic strain, rotations, etc.) in the discrete formulation, combined with inf-sup based stability/data, establishes MFEM as a mathematically rigorous and versatile computational framework, with proven optimality, adaptability, and robustness (Bammer et al., 17 Jan 2024).