A Staggered Discontinuous Galerkin Method for linear elasticity problem on Polytopal Meshes

Abstract: This paper develops a novel staggered discontinuous Galerkin (SDG) method for linear elasticity based on the Hellinger-Reissner variational principle. We construct symmetric stress spaces with normal continuity across element boundaries on arbitrary polytopal meshes, while approximating the displacement field using piecewise polynomial functions defined on the same meshes. The method is locking-free and satisfies a local balance of linear momentum and angular momentum. We present a comprehensive theoretical analysis, including proofs of stability and error estimates. The formulation admits a hybridizable structure, which significantly simplifies the numerical implementation. Numerical experiments validate the theoretical results and demonstrate the effectiveness of the proposed approach.

• The paper develops a staggered discontinuous Galerkin method using the Hellinger–Reissner formulation, with strongly symmetric, normally continuous stresses and locally preserved linear and angular momentum.
• The method is locking-free as the Lamé parameter grows, with stability and error bounds uniform in the incompressible limit and observed displacement convergence of order 2 for k=0 and k+2 generally.
• Hybridization enables efficient implementation without changing the SDG solution, while the construction excludes 2D rectangular and parallelogram meshes and leaves solver conditioning and mixed boundary conditions open.

Motivation and contribution

This paper develops a staggered discontinuous Galerkin (SDG) method for linear elasticity based on the Hellinger–Reissner mixed variational principle, targeting arbitrary polytopal meshes in dimension $d \geq 2$. The Hellinger–Reissner formulation seeks $$\boldsymbol{\sigma} \in H(\mathrm{div},\Omega;\mathbb{S})$$ and $\boldsymbol{u} \in L^2(\Omega;\mathbb{R}^d)$, and its discretization requires (i) an inf-sup condition, (ii) stress symmetry and normal continuity, and (iii) optimal convergence. Classical symmetric mixed elements (e.g., Arnold–Winther-type families) carry vertex degrees of freedom (DoFs), which obstruct hybridization. The authors' central design decision is to reverse the continuity assignment used in earlier SDG elasticity methods [Lee–Kim; Zhao–Park]: the stress is normal-continuous across primal element boundaries while remaining discontinuous elsewhere, and the displacement is continuous across dual-element boundaries. Crucially, the stress is strongly symmetric, whereas prior SDG elasticity schemes enforce only weak symmetry. A further consequence of the design is that both balance laws—linear momentum ($\mathrm{div}\,\boldsymbol{\sigma}=\boldsymbol{f}$) and angular momentum (symmetry of $\boldsymbol{\sigma}$)—are preserved locally, element by element.

The construction follows the primal–dual mesh paradigm: given a shape-regular polytopal primal mesh $\mathcal{K}_h$ of star-shaped elements, each element is refined by connecting its interior point $\boldsymbol{x}_K$ to its vertices, producing a simplicial mesh $\mathcal{T}_h$; the dual mesh $\mathcal{K}_h^*$ consists of the diamond-shaped patches $\omega_F$ around each primal face $$\boldsymbol{\sigma} \in H(\mathrm{div},\Omega;\mathbb{S})$$0. This extends the authors' recent SDG construction for the Poisson equation, but the extension to elasticity is nontrivial because it requires symmetric-tensor DoFs that are normal-continuous across polytopal boundaries.

Finite element spaces and degrees of freedom

For $$\boldsymbol{\sigma} \in H(\mathrm{div},\Omega;\mathbb{S})$$1, the local stress space $$\boldsymbol{\sigma} \in H(\mathrm{div},\Omega;\mathbb{S})$$2 is the discontinuous space $$\boldsymbol{\sigma} \in H(\mathrm{div},\Omega;\mathbb{S})$$3 on each sub-element $$\boldsymbol{\sigma} \in H(\mathrm{div},\Omega;\mathbb{S})$$4 of the refinement. The key structural result is a geometric decomposition separating the bubble space $$\boldsymbol{\sigma} \in H(\mathrm{div},\Omega;\mathbb{S})$$5 (tensors with vanishing normal trace on $$\boldsymbol{\sigma} \in H(\mathrm{div},\Omega;\mathbb{S})$$6) from face-trace contributions $$\boldsymbol{\sigma} \in H(\mathrm{div},\Omega;\mathbb{S})$$7, where $$\boldsymbol{\sigma} \in H(\mathrm{div},\Omega;\mathbb{S})$$8. This yields unisolvence via DoFs consisting of normal-trace moments on each face and interior moments against the bubble space, with a further splitting of the bubble DoFs that is used only to prove the discrete inf-sup condition.

A mesh assumption plays a genuine role here: the tangential symmetric matrices $$\boldsymbol{\sigma} \in H(\mathrm{div},\Omega;\mathbb{S})$$9 over the faces must span the full space $\boldsymbol{u} \in L^2(\Omega;\mathbb{R}^d)$0. The authors verify this for any 2D polygon with at least three pairwise non-collinear edges and for any 3D polyhedron (via a Gauss–Bonnet/angle-defect argument). They concede that the assumption fails for rectangles and parallelograms in 2D, since two independent tangential directions cannot span $\boldsymbol{u} \in L^2(\Omega;\mathbb{R}^d)$1; rectangular meshes therefore require the separate construction of Hu. This is a real restriction on the mesh class, though a mild one for general polytopal discretizations.

The lowest-order case $\boldsymbol{u} \in L^2(\Omega;\mathbb{R}^d)$2 requires special care because the naive trace space is too small. The authors enrich $\boldsymbol{u} \in L^2(\Omega;\mathbb{R}^d)$3 to $\boldsymbol{u} \in L^2(\Omega;\mathbb{R}^d)$4, where $\boldsymbol{u} \in L^2(\Omega;\mathbb{R}^d)$5 is built from the rigid-motion space. Two properties of this enrichment drive the analysis: the trace of $\boldsymbol{u} \in L^2(\Omega;\mathbb{R}^d)$6 is exactly $\boldsymbol{u} \in L^2(\Omega;\mathbb{R}^d)$7, and $\boldsymbol{u} \in L^2(\Omega;\mathbb{R}^d)$8, which makes the $\boldsymbol{u} \in L^2(\Omega;\mathbb{R}^d)$9 scheme's right-hand-side stabilization vanish. A piecewise Korn inequality, with jump terms projected onto $\mathrm{div}\,\boldsymbol{\sigma}=\boldsymbol{f}$0, supplies the norm equivalence $\mathrm{div}\,\boldsymbol{\sigma}=\boldsymbol{f}$1 used throughout the stability analysis.

The SDG scheme, stability, and locking-freeness

The discrete problem finds $\mathrm{div}\,\boldsymbol{\sigma}=\boldsymbol{f}$2, $\mathrm{div}\,\boldsymbol{\sigma}=\boldsymbol{f}$3 (piecewise $\mathrm{div}\,\boldsymbol{\sigma}=\boldsymbol{f}$4 displacements on the primal mesh) satisfying a mixed system with bilinear forms $\mathrm{div}\,\boldsymbol{\sigma}=\boldsymbol{f}$5 and $\mathrm{div}\,\boldsymbol{\sigma}=\boldsymbol{f}$6. The form $\mathrm{div}\,\boldsymbol{\sigma}=\boldsymbol{f}$7 contains the compliance term plus two stabilizers—$\mathrm{div}\,\boldsymbol{\sigma}=\boldsymbol{f}$8 and a jump penalty on interior dual faces. The authors emphasize that this stabilizer exists solely to prevent locking (coercivity of $\mathrm{div}\,\boldsymbol{\sigma}=\boldsymbol{f}$9 degenerates as $\boldsymbol{\sigma}$0), not to compensate for nonconformity as in virtual element or weak Galerkin methods; this distinction is what allows the stabilizer to be independent of the approximation space.

Well-posedness rests on three pillars. First, a discrete inf-sup condition $\boldsymbol{\sigma}$1 is proved by an explicit construction using the modified DoFs. Second, a trace bound shows that $\boldsymbol{\sigma}$2 is controlled by $\boldsymbol{\sigma}$3 plus the stabilization terms, uniformly in the Lamé constants; combined with the decomposition of $\boldsymbol{\sigma}$4 into deviatoric and trace parts, this yields discrete coercivity on the kernel space $\boldsymbol{\sigma}$5 with constants independent of $\boldsymbol{\sigma}$6 and $\boldsymbol{\sigma}$7. The implication is that the method is locking-free in the incompressible limit $\boldsymbol{\sigma}$8. Third, testing the first equation with $\boldsymbol{\sigma}$9 gives $\mathcal{K}_h$0, closing the argument. The stability estimate $\mathcal{K}_h$1 follows.

Hybridization and error estimates

The method admits a hybridizable structure that substantially simplifies implementation. Introducing the space $\mathcal{K}_h$2 with interior components $\mathcal{K}_h$3 and boundary components $\mathcal{K}_h$4, the authors define weak symmetric strain and weak divergence operators that are exact discrete adjoints of one another. Notably, computing $\mathcal{K}_h$5 on the enriched discontinuous space $\mathcal{K}_h$6 rather than $\mathcal{K}_h$7 eliminates the stabilization term that weak Galerkin methods require. The hybridized formulation on $\mathcal{K}_h$8 is shown to be well-posed, and its solution satisfies $\mathcal{K}_h$9 with $\boldsymbol{x}_K$0 coinciding exactly with the SDG solution—so hybridization is an implementation device, not an approximation.

The error analysis gives, for $\boldsymbol{x}_K$1 and $\boldsymbol{x}_K$2,

$\boldsymbol{x}_K$3

and, via a standard duality argument assuming the regularity estimate $\boldsymbol{x}_K$4 (stated without proof), an improved $\boldsymbol{x}_K$5 rate

$\boldsymbol{x}_K$6

Both estimates are uniform in the Lamé constants.

Numerical experiments

Two experiments on the unit square with a manufactured trigonometric solution confirm the theory. On triangular meshes at $\boldsymbol{x}_K$7, the $\boldsymbol{x}_K$8 displacement error converges at order 2.00 and, critically, the errors for $\boldsymbol{x}_K$9 are essentially identical at each mesh size (e.g., $\mathcal{T}_h$0 at the finest level for all four values of $\mathcal{T}_h$1), empirically validating locking-freeness. On genuinely polygonal meshes, the observed rates are $\mathcal{T}_h$2 for $\mathcal{T}_h$3 and $\mathcal{T}_h$4 and $\mathcal{T}_h$5 for $\mathcal{T}_h$6, matching the predicted orders at both $\mathcal{T}_h$7 and $\mathcal{T}_h$8. Experiments were implemented in the MATLAB package iFEM.

Limitations and open questions

Several restrictions should be noted. The mesh assumption excluding rectangles and parallelograms in 2D means the lowest-order element as constructed does not cover Cartesian meshes, for which separate elements are needed. The $\mathcal{T}_h$9 error estimate relies on a convex-domain-type dual regularity assumption that is asserted rather than proved, so the superconvergent rate is conditional on domain smoothness. The stability and error analysis is carried out for homogeneous Dirichlet conditions; extension to mixed boundary conditions is not treated. Finally, the paper does not address preconditioning or the conditioning of the hybridized Schur complement system, which is where the claimed computational savings would need to be quantified.

Conclusion

The paper constructs a locking-free, hybridizable SDG method for linear elasticity on polytopal meshes in which the stress is strongly symmetric and normal-continuous across primal faces, with both linear and angular momentum balanced locally. Optimal-order error estimates uniform in the Lamé parameters are proved, and numerical results on triangular and polygonal meshes confirm the predicted convergence rates and robustness for $\mathcal{K}_h^*$0 up to $\mathcal{K}_h^*$1. The remaining questions concern mesh classes excluded by the spanning assumption, non-homogeneous boundary conditions, and solver-level performance of the hybridized system.