A low-order nonconforming method for linear elasticity on general meshes

TL;DR

This paper introduces a low-order, nonconforming finite element method for linear elasticity that works on general meshes, ensuring stability and convergence without locking, supported by theoretical proofs and numerical validation.

Contribution

It develops a novel low-order nonconforming method for elasticity on arbitrary meshes, with a new stabilization term and proof of discrete Korn inequality for stability.

Findings

Locking-free error estimates in energy and L2 norms.

Convergence rates of h and h^2 for smooth solutions.

Numerical validation confirms theoretical results.

Abstract

In this work we construct a low-order nonconforming approximation method for linear elasticity problems supporting general meshes and valid in two and three space dimensions. The method is obtained by hacking the Hybrid High-Order method, that requires the use of polynomials of degree $k\ge1$ for stability. Specifically, we show that coercivity can be recovered for $k=0$ by introducing a novel term that penalises the jumps of the displacement reconstruction across mesh faces. This term plays a key role in the fulfillment of a discrete Korn inequality on broken polynomial spaces, for which a novel proof valid for general polyhedral meshes is provided. Locking-free error estimates are derived for both the energy- and the $L^2$-norms of the error, that are shown to convergence, for smooth solutions, as $h$ and $h^2$, respectively (here, $h$ denotes the meshsize). A thorough numerical validation on a complete panel of two- and three-dimensional test cases is provided.

Tables

Table 1. Table 1 : Numerical results for the test of Section 6.1 , structured triangular mesh family.

$N_{\mathrm{dofs},h}$	$N_{\mathrm{nz},h}$	${\Vert {\underset{¯}{𝒖}}_h-{\underset{¯}{𝑰}}_h𝒖\Vert }_{a,h}$	EOC	$\Vert {𝒖}_h-{𝝅}_h^0𝒖\Vert$	EOC
$(\mu ,\lambda )=(1,1)$
144	3680	3.82e+00	–	2.08e-01	–
608	17856	1.96e+00	0.97	6.97e-02	1.58
2496	78080	9.64e-01	1.02	1.87e-02	1.90
10112	326016	4.84e-01	1.00	4.74e-03	1.98
40704	1331840	2.43e-01	1.00	1.19e-03	1.99
$(\mu ,\lambda )=(1,1,000)$
144	3680	5.09e+00	–	2.05e-01	–
608	17856	1.95e+00	1.38	7.15e-02	1.52
2496	78080	9.15e-01	1.09	2.00e-02	1.84
10112	326016	4.52e-01	1.02	5.18e-03	1.95
40704	1331840	2.25e-01	1.00	1.31e-03	1.98
$(\mu ,\lambda )=(1,1\cdot {10}^6)$
144	3680	1.10e+02	–	2.05e-01	–
608	17856	1.48e+01	2.90	7.15e-02	1.52
2496	78080	2.07e+00	2.83	2.00e-02	1.84
10112	326016	5.08e-01	2.03	5.19e-03	1.95
40704	1331840	2.27e-01	1.16	1.31e-03	1.98

Table 2. Table 2 : Numerical results for the test of Section 6.1 , unstructured triangular mesh family.

$N_{\mathrm{dofs},h}$	$N_{\mathrm{nz},h}$	${\Vert {\underset{¯}{𝒖}}_h-{\underset{¯}{𝑰}}_h𝒖\Vert }_{a,h}$	EOC	$\Vert {𝒖}_h-{𝝅}_h^0𝒖\Vert$	EOC
$(\mu ,\lambda )=(1,1)$
234	6572	3.00e+00	–	1.38e-01	–
978	30012	1.60e+00	0.90	4.11e-02	1.75
3986	127372	8.15e-01	0.98	9.37e-03	2.13
15542	505828	4.27e-01	0.93	2.61e-03	1.85
63584	2089920	2.12e-01	1.01	6.65e-04	1.97
249238	8228988	1.08e-01	0.97	1.71e-04	1.96
$(\mu ,\lambda )=(1,1,000)$
234	6572	3.57e+00	–	1.45e-01	–
978	30012	1.60e+00	1.15	4.52e-02	1.68
3986	127372	8.00e-01	1.00	1.07e-02	2.07
15542	505828	4.18e-01	0.94	2.99e-03	1.85
63584	2089920	2.08e-01	1.01	7.63e-04	1.97
249238	8228988	1.06e-01	0.97	1.97e-04	1.96
$(\mu ,\lambda )=(1,1\cdot {10}^6)$
234	6572	6.17e+01	–	1.45e-01	–
978	30012	7.55e+00	3.03	4.52e-02	1.68
3986	127372	1.14e+00	2.72	1.07e-02	2.07
15542	505828	4.33e-01	1.40	2.99e-03	1.85
63584	2089920	2.08e-01	1.06	7.63e-04	1.97
249238	8228988	1.06e-01	0.98	1.97e-04	1.96

Table 3. Table 3 : Numerical results for the test of Section 6.1 , Cartesian orthogonal mesh family.

$N_{\mathrm{dofs},h}$	$N_{\mathrm{nz},h}$	${\Vert {\underset{¯}{𝒖}}_h-{\underset{¯}{𝑰}}_h𝒖\Vert }_{a,h}$	EOC	$\Vert {𝒖}_h-{𝝅}_h^0𝒖\Vert$	EOC
$(\mu ,\lambda )=(1,1)$
80	2768	3.13e+00	–	1.55e-01	–
352	15856	1.84e+00	0.77	4.08e-02	1.93
1472	73904	1.09e+00	0.75	1.04e-02	1.98
6016	317488	5.89e-01	0.89	2.89e-03	1.84
24320	1314608	3.02e-01	0.97	7.73e-04	1.90
$(\mu ,\lambda )=(1,1,000)$
80	2768	3.08e+00	–	1.64e-01	–
352	15856	1.81e+00	0.77	4.72e-02	1.80
1472	73904	1.08e+00	0.75	1.37e-02	1.78
6016	317488	5.81e-01	0.89	3.96e-03	1.79
24320	1314608	2.97e-01	0.97	1.06e-03	1.90
$(\mu ,\lambda )=(1,1\cdot {10}^6)$
80	2768	3.08e+00	–	1.64e-01	–
352	15856	1.81e+00	0.77	4.72e-02	1.80
1472	73904	1.08e+00	0.75	1.37e-02	1.78
6016	317488	5.81e-01	0.89	3.96e-03	1.79
24320	1314608	2.97e-01	0.97	1.06e-03	1.90

Table 4. Table 4 : Numerical results for the test of Section 6.1 , distorted quadrangular mesh family.

$N_{\mathrm{dofs},h}$	$N_{\mathrm{nz},h}$	${\Vert {\underset{¯}{𝒖}}_h-{\underset{¯}{𝑰}}_h𝒖\Vert }_{a,h}$	EOC	$\Vert {𝒖}_h-{𝝅}_h^0𝒖\Vert$	EOC
$(\mu ,\lambda )=(1,1)$
80	2768	3.51e+00	–	1.89e-01	–
352	15856	1.91e+00	0.88	5.45e-02	1.79
1472	73904	1.08e+00	0.82	1.34e-02	2.03
6016	317488	5.83e-01	0.89	3.52e-03	1.93
24320	1314608	2.97e-01	0.97	9.18e-04	1.94
97792	5348656	1.49e-01	0.99	2.33e-04	1.98
$(\mu ,\lambda )=(1,1,000)$
80	2768	3.44e+00	–	1.96e-01	–
352	15856	1.87e+00	0.88	5.89e-02	1.73
1472	73904	1.07e+00	0.81	1.63e-02	1.85
6016	317488	5.74e-01	0.89	4.48e-03	1.86
24320	1314608	2.92e-01	0.97	1.18e-03	1.93
97792	5348656	1.47e-01	0.99	3.00e-04	1.97
$(\mu ,\lambda )=(1,1\cdot {10}^6)$
80	2768	9.12e+00	–	1.96e-01	–
352	15856	2.27e+00	2.00	5.89e-02	1.73
1472	73904	1.08e+00	1.07	1.63e-02	1.85
6016	317488	5.74e-01	0.91	4.48e-03	1.86
24320	1314608	2.92e-01	0.97	1.18e-03	1.93
97792	5348656	1.47e-01	0.99	3.00e-04	1.97

Table 5. Table 5 : Numerical results for the test of Section 6.2 and comparison with the high-order method of [ 18 ] with k = 1 𝑘 1 k=1 . For the latter, the energy norm is the one associated to the corresponding bilinear form without jump stabilisation.

$N_{\mathrm{dofs},h}$	$N_{\mathrm{nz},h}$	${\Vert {\underset{¯}{𝒖}}_h-{\underset{¯}{𝑰}}_h𝒖\Vert }_{a,h}$	EOC	$\Vert {𝒖}_h-{𝝅}_h^0𝒖\Vert$	EOC
Present work   -
256	10616	7.65e-01	–	7.51e-02	–
1088	52728	5.63e-01	0.44	3.34e-02	1.17
4480	232568	3.97e-01	0.50	1.40e-02	1.25
18176	974712	2.76e-01	0.53	5.72e-03	1.29
73216	3988856	1.90e-01	0.54	2.31e-03	1.31
293888	16136568	1.31e-01	0.54	9.29e-04	1.31
HHO method of [18], $k=1$
320	7584	1.07e-01	–	9.40e-03	–
1408	36512	7.32e-02	0.55	3.64e-03	1.37
5888	158880	5.01e-02	0.55	1.41e-03	1.36
24064	661664	3.43e-02	0.55	5.52e-04	1.36
97280	2699424	2.35e-02	0.54	2.17e-04	1.35
391168	10903712	1.61e-02	0.54	8.57e-05	1.34

Table 6. Table 6 : Numerical results for the test of Section 6.3 . The number of degrees of freedom and of nonzero matrix entries for the method of [ 18 ] are also included for comparison (except for the last mesh refinement).

$N_{\mathrm{dofs},h}$		$N_{\mathrm{nz},h}$		${\Vert {\underset{¯}{𝒖}}_h-{\underset{¯}{𝑰}}_h𝒖\Vert }_{a,h}$	EOC	$\Vert {𝒖}_h-{𝝅}_h^0𝒖\Vert$	EOC
$k=0$	($k=1$)	$k=0$	($k=1$)
Cartesian orthogonal mesh sequence
60	(108)	2772	(4860)	2.42e+00	–	1.76e-01	–
624	(1296)	70128	(97200)	2.07e+00	0.23	1.01e-01	0.81
5568	(12096)	831024	(1057536)	1.31e+00	0.65	4.09e-02	1.30
46848	(103680)	7879824	(9673344)	7.19e-01	0.87	1.27e-02	1.68
384000	(857088)	68277456	(82425600)	3.71e-01	0.95	3.46e-03	1.88
3108864	–	567808848	–	1.87e-01	0.98	8.95e-04	1.95
Unstructured simplicial mesh sequence
1584	(3024)	107136	(167184)	1.38e+00	–	4.70e-02	–
13248	(25920)	1008288	(1539648)	7.61e-01	0.85	1.64e-02	1.52
108288	(214272)	8676288	(13125888)	3.96e-01	0.94	4.39e-03	1.91
875520	(1741824)	71860608	(108241920)	2.02e-01	0.97	1.14e-03	1.95
7041024	–	584706816	–	1.02e-01	0.99	2.89e-04	1.98

Equations

$$-\nabla \cdot (\mathbf{\sigma }({\nabla }_{\mathrm{s}}\mathbf{u}))$$

$$\mathbf{u}$$

$$\mathbf{\sigma }(\mathbf{\tau })=2\mu \mathbf{\tau }+\lambda \mathrm{tr}(\mathbf{\tau }){\mathbf{I}}_d,$$

$$2\mu -d{\lambda }^{-}\ge \alpha ,$$

$$(\mathbf{\sigma }({\nabla }_{\mathrm{s}}\mathbf{u}),{\nabla }_{\mathrm{s}}\mathbf{v})=(\mathbf{f},\mathbf{v})\forall \mathbf{v}\in H_0^1(\Omega ;{\mathbb{R}}^d),$$

$${\mathbb{P}}^l({\mathcal{T}}_h;\mathbb{R}):-\left\{v_h\in L^2(\Omega ;\mathbb{R}):v_{h|T}\in {\mathbb{P}}^l(T;\mathbb{R})\forall T\in {\mathcal{T}}_h\right\}.$$

$$H^s({\mathcal{T}}_h;\mathbb{R}):-\left\{v\in L^2(\Omega ;\mathbb{R}):v_{|T}\in H^s(T;\mathbb{R})\forall T\in {\mathcal{T}}_h\right\},$$

$$|v{|}_{H^s({\mathcal{T}}_h;\mathbb{R})}:-{\left(\sum _{T\in {\mathcal{T}}_h}|v{|}_{H^s(T;\mathbb{R})}^2\right)}^{\frac{1}{2}}.$$

$${\pi }_X^{0}v:-\frac{1}{|X|}{\int }_{X}v,$$

$$\Vert v-{\pi }_T^{0}v{\Vert }_{L^2(T;\mathbb{R})}+h_T^{\frac{1}{2}}\Vert v-{\pi }_T^{0}v{\Vert }_{L^2(\partial T;\mathbb{R})}\lesssim h_T|v{|}_{H^1(T;\mathbb{R})},$$

$$({\pi }_h^{0}v{)}_{|T}:-{\pi }_T^0{v}_{|T}\forall T\in {\mathcal{T}}_h.$$

$$\nabla {\varpi }_T^{1}v={\mathbf{\pi }}_T^0(\nabla v)\text{mbox}and\frac{1}{|T|}{\int }_T{\varpi }_T^{1}v=\frac{1}{|T|}{\int }_{T}v.$$

$$\Vert v-{\varpi }_T^{1}v{\Vert }_{L^2(T;\mathbb{R})}+h_T^{\frac{1}{2}}\Vert v-{\varpi }_T^{1}v{\Vert }_{L^2(\partial T;\mathbb{R})}\lesssim h_T^2|v{|}_{H^2(T;\mathbb{R})},$$

$$({\varpi }_h^{1}v{)}_{|T}:-{\varpi }_T^1{v}_{|T}\forall T\in {\mathcal{T}}_h.$$

$$[v{]}_F:-(v_{|T_1}{)}_{|F}-(v_{|T_2}{)}_{|F}.$$

$$[v{]}_F:-v_{|F}.$$

$$(I_{\mathrm{av},h}^l{v}_h)(V):-\{\begin{array}{ll}\frac{1}{\mathrm{card}({\mathfrak{T}}_V)}\sum _{\tau \in {\mathfrak{T}}_V}(v_h{)}_{|\tau }(V)& \text{if V∈Ω},\\ 0& \text{if V∈∂Ω}.\end{array}$$

$$\Vert v_h-I_{\mathrm{av},h}^l{v}_h{\Vert }_T^2\lesssim \sum _{F\in {\mathcal{F}}_{\mathcal{V},T}}{h}_F\Vert [v_h{]}_F{\Vert }_F^2,$$

$$|v_h-I_{\mathrm{av},h}^l{v}_h{|}_{H^1({\mathcal{T}}_h;\mathbb{R})}^2$$

$$\lesssim \sum _{T\in {\mathcal{T}}_h}{h}_T^{-2}\sum _{F\in {\mathcal{F}}_{\mathcal{V},T}}{h}_F\Vert [v_h{]}_F{\Vert }_F^2$$

$$\lesssim \sum _{F\in {\mathcal{F}}_h}\sum _{T\in {\mathcal{T}}_{\mathcal{V},F}}{h}_F^{-1}\Vert [v_h{]}_F{\Vert }_F^2,$$

$$|v_h-I_{\mathrm{av},h}^l{v}_h{|}_{H^1({\mathcal{T}}_h;\mathbb{R})}^2\lesssim \sum _{F\in {\mathcal{F}}_h}{h}_F^{-1}\Vert [v_h{]}_F{\Vert }_F^2.$$

$$\Vert {\mathbf{v}}_h{\Vert }_{\mathbf{\epsilon },h}:-{\left(\Vert {\nabla }_{\mathrm{s},h}{\mathbf{v}}_h{\Vert }^2+|{\mathbf{v}}_h{|}_{\mathrm{j},h}^2\right)}^{\frac{1}{2}}\text{mbox}and|{\mathbf{v}}_h{|}_{\mathrm{j},h}:-{\left(\sum _{F\in {\mathcal{F}}_h}{h}_F^{-1}\Vert [{\mathbf{v}}_h{]}_F{\Vert }_F^2\right)}^{\frac{1}{2}},$$

$$|{\mathbf{v}}_h{|}_{H^1({\mathcal{T}}_h;{\mathbb{R}}^d)}\lesssim \Vert {\mathbf{v}}_h{\Vert }_{\mathbf{\epsilon },h}.$$

$$|{\mathbf{v}}_h{|}_{H^1({\mathcal{T}}_h;{\mathbb{R}}^d)}^2$$

$$\lesssim \Vert {\nabla }_{\mathrm{s}}{\mathbf{I}}_{\mathrm{av},h}^l{\mathbf{v}}_h{\Vert }^2+|{\mathbf{v}}_h{|}_{\mathrm{j},h}^2$$

$$\lesssim \Vert {\nabla }_{\mathrm{s},h}{\mathbf{v}}_h{\Vert }^2+\Vert {\nabla }_{\mathrm{s},h}({\mathbf{I}}_{\mathrm{av},h}^l{\mathbf{v}}_h-{\mathbf{v}}_h){\Vert }^2+|{\mathbf{v}}_h{|}_{\mathrm{j},h}^2$$

$$\lesssim \Vert {\nabla }_{\mathrm{s},h}{\mathbf{v}}_h{\Vert }^2+|{\mathbf{v}}_h{|}_{\mathrm{j},h}^2=\Vert {\mathbf{v}}_h{\Vert }_{\mathbf{\epsilon },h}^2,$$

$$\Vert {\mathbf{v}}_h\Vert \lesssim \Vert {\mathbf{v}}_h{\Vert }_{\mathbf{\epsilon },h},$$

$${\underline{\mathbf{U}}}_h:-\left\{{\underline{\mathbf{v}}}_h=(({\mathbf{v}}_T{)}_{T\in {\mathcal{T}}_h},({\mathbf{v}}_F{)}_{F\in {\mathcal{F}}_h}):{\mathbf{v}}_T\in {\mathbb{P}}^0(T;{\mathbb{R}}^d)\forall T\in {\mathcal{T}}_h\text{mbox}and{\mathbf{v}}_F\in {\mathbb{P}}^0(F;{\mathbb{R}}^d)\forall F\in {\mathcal{F}}_h\right\}.$$

$$({\mathbf{v}}_h{)}_{|T}:-{\mathbf{v}}_T\forall T\in {\mathcal{T}}_h.$$

$${\underline{\mathbf{I}}}_h\mathbf{v}:-(({\mathbf{\pi }}_T^0{\mathbf{v}}_{|T}{)}_{T\in {\mathcal{T}}_h},({\mathbf{\pi }}_F^0{\mathbf{v}}_{|F}{)}_{F\in {\mathcal{F}}_h}).$$

$${\underline{\mathbf{I}}}_T\mathbf{v}=({\mathbf{\pi }}_T^0\mathbf{v},({\mathbf{\pi }}_F^0{\mathbf{v}}_{|F}{)}_{F\in {\mathcal{F}}_T}).$$

$${\underline{\mathbf{U}}}_{h,0}:-\left\{{\underline{\mathbf{v}}}_h\in {\underline{\mathbf{U}}}_h:{\mathbf{v}}_F=0\forall F\in {\mathcal{F}}_h^{\mathrm{b}}\right\}.$$

Full text

Key words. Linear elasticity, Korn’s inequality, locking-free methods, Hybrid High-Order methods, polyhedral meshes

AMS subject classification. 65N08, 65N30, 74B05, 74G15

1 Introduction

Discretisation methods supporting meshes with general, possibly non standard, element shapes have experienced a vigorous growth over the last few years. In the context of solid-mechanics, this feature can be useful for several reasons including, e.g., improved robustness to mesh distortion and fracture, local mesh refinement, or the use of hanging nodes for contact and interface problems. A non-exahustive list of contributions in the context of elasticity problems includes [26, 30, 4, 5, 23, 24, 20, 18, 6, 3, 28, 11, 29, 10]; see also references therein.

For large three-dimensional simulations, or whenever one cannot expect the exact solution to be smooth, low-order methods are often privileged in order to reduce the number of unknowns. It is well-known, however, that low-order Finite Element (FE) approximations are in some cases unsatisfactory: affine conforming FE methods are not robust in the quasi-incompressible limit owing to their inability to represent non-trivial divergence-free displacement fields; nonconforming (Crouzeix–Raviart) FE methods, on the other hand, yield unstable discretisations unless appropriate measures are taken; see, e.g., the discussions in [8, 25]. The underlying reason for this lack of stability is the non-fulfillment of a discrete counterpart of Korn’s inequality owing to a poor control of rigid-body motions at mesh faces. For similar reasons, the stability of Hybrid High-Order (HHO) methods for linear elasticity requires the use of polynomials of degree $k\geq 1$ as unknowns; see [18, Lemma 4]. As a matter of fact, as we show in Section 4.4 below, the stability and consistency requirements on the local HHO stabilisation term are incompatible when $k=0$, that is, when piecewise constant polynomials on the mesh and its skeleton are used as discrete unknowns.

In this paper we highlight a modification of the HHO method which recovers stability for $k=0$. The proposed fix consists in adding a novel term which penalises in a least square sense the jumps of the local affine displacement reconstruction. This modification is inspired by the Korn inequality on broken polynomial spaces proved in Lemma 1 below, which appears to be a novel extension of similar results to general polyhedral meshes. The proof combines the techniques of [7, Lemma 2.2] with the recent results of [17] and [12] concerning, respectively, the node-averaging operator and local inverse inequalities on polyhedral meshes. In the context of Crouzeix–Raviart FE approximations of linear elasticity problems on standard meshes, similar jump penalisation terms have been considered in [25].

The resulting method has several appealing features: it is valid in two and three space dimensions, paving the way to unified implementations; it hinges on a reduced number of unknowns ($15$ for a tetrahedron, $21$ for a hexahedron and, for more general polyhedral shapes, $3$ unknowns per face plus $3$ unknowns inside the element); it is robust in the quasi-incompressible limit; it admits a formulation in terms of conservative numerical tractions, which enables its integration in existing Finite Volume simulators (a particularly relevant feature in the context of industrial applications).

We carry out a complete convergence analysis based on the abstract framework of [14] for methods in fully discrete formulation. Specifically, we show that the energy and $L^{2}$-norms of the error converge, respectively, as $h$ and $h^{2}$ (with $h$ denoting, as usual, the meshsize). As for the original HHO method of [18], the error estimates are additionally shown to be robust in the quasi-incompressible limit. Key to this result is the fact that the gradient of the local displacement reconstruction satisfies a suitable commutation property with the $L^{2}$-orthogonal projector. The theoretical results are supported by a thorough numerical investigation, including two- and three-dimensional test cases, as well as a comparison with the original HHO method of [18] on a test case mimicking a mode 1 fracture.

The rest of the paper is organised as follows. In Section 2 we formulate the continuous problem along with the assumptions on the problem data. In Section 3 we establish the discrete setting: after briefly recalling the notion of regular polyhedral mesh, we introduce local and broken polynomial spaces and projectors thereon, and we prove a discrete counterpart of Korn’s first inequality on broken polynomial spaces. In Section 4 we introduce the space of discrete unknowns, define a local affine displacement reconstruction, formulate the discrete bilinear form, discuss the differences with respect to the original HHO bilinear form of [18], and state the discrete problem. Section 5 addresses the convergence analysis of the method in the energy- and $L^{2}$-norms, while Section 6 contains an exhaustive panel of two- and three-dimensional numerical tests. Finally, in Section 7 we show that the method satisfies local balances with equilibrated tractions, for which an explicit expression is provided.

2 Continuous setting

Consider a body which, in its reference configuration, occupies a given region of space $\Omega\subset\mathbb{R}^{d}$, $d\in\{2,3\}$. In what follows, it is assumed that $\Omega$ is a bounded connected open polygonal (if $d=2$) or polyhedral (if $d=3$) set that does not have cracks, i.e., it lies on one side of its boundary $\partial\Omega$. We are interested in finding the displacement field $\boldsymbol{u}:\Omega\to\mathbb{R}^{d}$ of the body when it is subjected to a given force per unit volume $\boldsymbol{f}:\Omega\to\mathbb{R}^{d}$. We work in what follows under the small deformation assumption which implies, in particular, that the strain tensor $\boldsymbol{\varepsilon}$ is given by the symmetric part of the gradient of the displacement field, i.e., $\boldsymbol{\varepsilon}=\boldsymbol{\nabla}_{\rm s}\boldsymbol{u}$ where, for any vector-valued function $\boldsymbol{z}=(z_{i})_{1\leq i\leq d}$ smooth enough, we have set $\boldsymbol{\nabla}\boldsymbol{z}=(\partial_{j}z_{i})_{1\leq i,j\leq d}$ and $\boldsymbol{\nabla}_{\rm s}\boldsymbol{z}\coloneq\frac{1}{2}\left(\boldsymbol{\nabla}\boldsymbol{z}+\boldsymbol{\nabla}\boldsymbol{z}^{\top}\right)$. We further assume, for the sake of simplicity, that the body is clamped along its boundary $\partial\Omega$. Other standard boundary conditions can be considered up to minor modifications. The displacement field is obtained by solving the following linear elasticity problem, which expresses the equilibrium between internal stresses and external loads: Find $\boldsymbol{u}:\Omega\to\mathbb{R}^{d}$ such that

[TABLE]

where, denoting by $\mathbb{R}_{\rm sym}^{d\times d}$ the set of symmetric real-valued $d\times d$ matrices, the mapping $\boldsymbol{\sigma}:\mathbb{R}_{\rm sym}^{d\times d}\to\mathbb{R}_{\rm sym}^{d\times d}$ represents the strain-stress law. For isotropic homogeneous materials, the strain-stress law is such that, for any $\boldsymbol{\tau}\in\mathbb{R}_{\rm sym}^{d\times d}$,

[TABLE]

where $\operatorname{tr}(\boldsymbol{\tau})\coloneq\sum_{i=1}^{d}\tau_{ii}$ is the trace operator and $\boldsymbol{I}_{d}$ the $d\times d$ identity matrix. The real numbers $\mu$ and $\lambda$, which correspond to the Lamé coefficients when $d=3$, are assumed such that, for a real number $\alpha>0$,

[TABLE]

where $\lambda^{-}\coloneq\frac{1}{2}\left(|\lambda|-\lambda\right)$ denotes the negative part of $\lambda$. In what follows, $\mu$, $\lambda$, the related bound (3), and $\boldsymbol{f}$ will be collectively referred to as the problem data.

For any open bounded set $X\subset\Omega$, we denote by $({\cdot},{\cdot})_{X}$ the usual inner product of the space of scalar-valued, square-integrable functions $L^{2}(X;\mathbb{R})$, by $\|{\cdot}\|_{X}$ the corresponding norm, and we adopt the convention that the subscript is omitted whenever $X=\Omega$. The same notation is used for the spaces of vector- and tensor-valued square-integrable functions $L^{2}(X;\mathbb{R}^{d})$ and $L^{2}(X;\mathbb{R}^{d\times d})$, respectively. With this notation, a classical weak formulation of problem (1) reads: Find $\boldsymbol{u}\in H_{0}^{1}(\Omega;\mathbb{R}^{d})$ such that

[TABLE]

where $H_{0}^{1}(\Omega;\mathbb{R}^{d})$ classically denotes the space of vector-valued functions that are square-integrable along with all their partial derivatives, and whose traces on $\partial\Omega$ vanish.

3 Discrete setting

3.1 Mesh

Throughout the rest of the paper, we will use for the sake of simplicity the three-dimensional nomenclature also when $d=2$, i.e., we will speak of polyhedra and faces rather than polygons and edges. We consider here meshes corresponding to couples $\mathcal{M}_{h}\coloneq(\mathcal{T}_{h},\mathcal{F}_{h})$, where $\mathcal{T}_{h}$ is a finite collection of polyhedral elements $T$ such that $h\coloneq\max_{T\in\mathcal{T}_{h}}h_{T}>0$ with $h_{T}$ denoting the diameter of $T$, while $\mathcal{F}_{h}$ is a finite collection of planar faces $F$. It is assumed henceforth that the mesh $\mathcal{M}_{h}$ matches the geometrical requirements detailed in [22, Definition 7.2]; see also [21, Section 2]. This covers, essentially, any reasonable partition of $\Omega$ into polyhedral sets, not necessarily convex or even star-shaped. For every mesh element $T\in\mathcal{T}_{h}$, we denote by $\mathcal{F}_{T}$ the subset of $\mathcal{F}_{h}$ containing the faces that lie on the boundary $\partial T$ of $T$. Symmetrically, for every face $F\in\mathcal{F}_{h}$, we denote by $\mathcal{T}_{F}$ the subset of $\mathcal{T}_{h}$ containing the (one or two) mesh elements that share $F$. For any mesh element $T\in\mathcal{T}_{h}$ and each face $F\in\mathcal{F}_{T}$, $\boldsymbol{n}_{TF}$ is the constant unit normal vector to $F$ pointing out of $T$. Boundary faces lying on $\partial\Omega$ and internal faces contained in $\Omega$ are collected in the sets $\mathcal{F}_{h}^{{\rm b}}$ and $\mathcal{F}_{h}^{{\rm i}}$, respectively. For any $F\in\mathcal{F}_{h}^{{\rm i}}$, we denote by $T_{1}$ and $T_{2}$ the elements of $\mathcal{T}_{h}$ such that $F\subset\partial T_{1}\cap\partial T_{2}$. The numbering of $T_{1}$ and $T_{2}$ is assumed arbitrary but fixed, and we set $\boldsymbol{n}_{F}\coloneq\boldsymbol{n}_{T_{1}F}$. Our focus is on the $h$-convergence analysis, so we consider a sequence of refined meshes that is regular in the sense of [21, Definition 3]. This implies, in particular, that the diameter $h_{T}$ of a mesh element $T\in\mathcal{T}_{h}$ is comparable to the diameter $h_{F}$ of each face $F\in\mathcal{F}_{T}$ uniformly in $h$, and that the number of faces in $\mathcal{F}_{T}$ is bounded above by an integer $N_{\partial}$ independent of $h$.

3.2 Local and broken spaces and projectors

In order to alleviate the exposition, throughout the rest of the paper we use the abridged notation $a\lesssim b$ for the inequality $a\leq Cb$ with real number $C>0$ independent of the meshsize, possibly on the problem data, and, for local inequalities, on the mesh element or face. We also write $a\simeq b$ for $a\lesssim b$ and $b\lesssim a$. The dependencies of the hidden constant are further specified whenever needed.

Let $X$ denote a mesh element or face. For a given integer $l\geq 0$, we denote by $\mathbb{P}^{l}(X;\mathbb{R})$ the space spanned by the restriction to $X$ of $d$-variate, real-valued polynomials of total degree $\leq l$. The corresponding spaces of vector- and tensor-valued functions are respectively denoted by $\mathbb{P}^{l}(X;\mathbb{R}^{d})$ and $\mathbb{P}^{l}(X;\mathbb{R}^{d\times d})$. A similar notation is used also for the vector and tensor versions of the broken spaces introduced in what follows. At the global level, we denote by $\mathbb{P}^{l}(\mathcal{T}_{h};\mathbb{R})$ the space of broken polynomials on $\mathcal{T}_{h}$ whose restriction to every mesh element $T\in\mathcal{T}_{h}$ lies in $\mathbb{P}^{l}(T;\mathbb{R})$, i.e.,

[TABLE]

We also introduce the broken Sobolev spaces

[TABLE]

which will be used in the error estimates to express the regularity requirements on the exact solution. On $H^{s}(\mathcal{T}_{h};\mathbb{R})$, we define the broken seminorm

[TABLE]

Again denoting by $X$ a mesh element or face, the local $L^{2}$-orthogonal projector $\pi_{X}^{0}:L^{2}(X;\mathbb{R})\to\mathbb{P}^{0}(X;\mathbb{R})$ maps every $v\in L^{2}(X;\mathbb{R})$ onto the constant function equal to its mean value inside $T$, that is,

[TABLE]

with $|X|$ denoting the Hausdorff measure of $X$. The vector and tensor versions of the $L^{2}$-projector, both denoted by $\boldsymbol{\pi}_{X}^{0}$, are obtained applying $\pi_{X}^{0}$ component-wise. From [12, Lemmas 3.4 and 3.6], it can be deduced that, for any mesh element $T\in\mathcal{T}_{h}$ and any function $v\in H^{1}(T;\mathbb{R})$, the following approximation properties hold:

[TABLE]

where $\partial T$ denotes the boundary of $T$ and the hidden constant is independent of $h$, $T$, and $v$. The global $L^{2}$-orthogonal projector $\pi_{h}^{0}:L^{2}(\Omega;\mathbb{R})\to\mathbb{P}^{0}(\mathcal{T}_{h};\mathbb{R})$ is such that, for any $v\in L^{2}(\Omega;\mathbb{R}^{d})$,

[TABLE]

The vector and tensor versions, both denoted by $\boldsymbol{\pi}_{h}^{0}$, are obtained applying $\pi_{h}^{0}$ component-wise.

We will also need the elliptic projector $\varpi_{T}^{1}:H^{1}(T;\mathbb{R})\to\mathbb{P}^{1}(T;\mathbb{R})$ such that, for all $v\in H^{1}(T;\mathbb{R})$,

[TABLE]

The first relation makes sense since $\boldsymbol{\nabla}\mathbb{P}^{1}(T;\mathbb{R})=\mathbb{P}^{0}(T;\mathbb{R}^{d})$, and it defines $\varpi_{T}^{1}v$ up to a constant, which is then fixed by the second relation. Also in this case, the vector version $\boldsymbol{\varpi}_{T}^{1}$ of the projector is obtained applying the scalar version component-wise. The following approximation properties for the elliptic projector are a special case of [13, Theorems 1.1 and 1.2]: For all $T\in\mathcal{T}_{h}$ and all $v\in H^{2}(T;\mathbb{R})$,

[TABLE]

where the hidden constant is independent of $h$, $T$, and $v$. For further use, we also define the global elliptic projector $\varpi_{h}^{1}:H^{1}(\mathcal{T}_{h};\mathbb{R})\to\mathbb{P}^{1}(\mathcal{T}_{h};\mathbb{R})$ such that, for any $v\in H^{1}(\mathcal{T}_{h};\mathbb{R})$,

[TABLE]

The vector version $\boldsymbol{\varpi}_{h}^{1}$ of the global elliptic projector is obtained applying $\varpi_{h}^{1}$ component-wise.

3.3 Discrete Korn inequality on broken polynomial spaces

The stability of our method hinges on a discrete counterpart of Korn’s inequality in discrete polynomial spaces stating that the $H^{1}$-seminorm of a vector-valued broken polynomial function is controlled by a suitably defined strain norm. The goal of this section is to prove this inequality.

Let us start with some preliminary results. Recalling that, for any $F\in\mathcal{F}_{h}^{{\rm i}}$, we have denoted by $T_{1}$ and $T_{2}$ the elements sharing $F$ and assumed that the ordering is arbitrary but fixed, we introduce the jump operator such that, for any function $v$ smooth enough to admit a (possibly two-valued) trace on $F$,

[TABLE]

When applied to vector-valued functions, the jump operator acts componentwise.

Let now $\mathfrak{T}_{h}$ denote a matching simplicial submesh of $\mathcal{M}_{h}$ in the sense of [21, Definition 4.2], and let $\mathfrak{F}_{h}$ be the corresponding set of simplicial faces. Given an integer $l\geq 1$, we define the node-averaging operator $I_{{\rm av},h}^{l}:\mathbb{P}^{l}(\mathcal{T}_{h};\mathbb{R})\to\mathbb{P}^{l}(\mathcal{T}_{h};\mathbb{R})\cap H_{0}^{1}(\Omega)$ such that, for any function $v_{h}\in\mathbb{P}^{l}(\mathcal{T}_{h};\mathbb{R})$ and any Lagrange node $V$ of $\mathfrak{T}_{h}$, denoting by $\mathfrak{T}_{V}$ the set of simplices sharing $V$,

[TABLE]

The vector-version, denoted by $\boldsymbol{I}_{{\rm av},h}^{l}$, acts component-wise. Adapting the reasoning of [17, Section 5.5.2] (based in turn on [27]), we infer that it holds, for all $T\in\mathcal{T}_{h}$,

[TABLE]

where $\mathcal{F}_{\mathcal{V},T}$ denotes the set of faces whose closure has nonempty intersection with the closure of $T$ and the hidden constant is independent of $h$, $T$, and $v_{h}$. Combining this result with an inverse inequality (see [12, Remark A.2]) we obtain, with hidden constants as before,

[TABLE]

where we have used (11) to pass to the second line while, to pass to the third line, we have invoked the mesh regularity to write $h_{F}h_{T}^{-2}\lesssim h_{F}^{-1}$ and we have exchanged the order of the sums after introducing the notation $\mathcal{T}_{\mathcal{V},F}$ for the set of mesh elements whose closure has nonzero intersection with the closure of $F$. Using again mesh regularity to infer that $\operatorname{card}(\mathcal{T}_{\mathcal{V},F})$ is bounded uniformly in $h$, we arrive at

[TABLE]

We are now ready to prove the discrete Korn inequality.

Lemma 1** (Discrete Korn inequality).**

Let an integer $l\geq 1$ be fixed and set, for all $\boldsymbol{v}_{h}\in\mathbb{P}^{l}(\mathcal{T}_{h};\mathbb{R}^{d})$,

[TABLE]

where $\boldsymbol{\nabla}_{{\rm s},h}:H^{1}(\mathcal{T}_{h};\mathbb{R}^{d})\to L^{2}(\Omega;\mathbb{R}_{\rm sym}^{d\times d})$ is the broken symmetric gradient such that $(\boldsymbol{\nabla}_{{\rm s},h}\boldsymbol{v})_{|T}=\boldsymbol{\nabla}_{\rm s}\boldsymbol{v}_{|T}$ for any $T\in\mathcal{T}_{h}$. Then, for all $\boldsymbol{v}_{h}\in\mathbb{P}^{l}(\mathcal{T}_{h};\mathbb{R}^{d})$, it holds with hidden constant depending only on $\Omega$, $d$, and the mesh regularity parameter:

[TABLE]

Proof.

The proof adapts the arguments of [7, Lemma 2.2]. We can write

[TABLE]

where we have inserted $\pm\boldsymbol{I}_{{\rm av},h}^{l}\boldsymbol{v}_{h}$ into the seminorm and used a triangle inequality in the first line, we have applied the first Korn inequality in $H^{1}(\Omega;\mathbb{R}^{d})$ to the first term and invoked (12) for the second term after recalling the definition (13) of the jump seminorm in the second line, we have inserted $\pm\boldsymbol{\nabla}_{{\rm s},h}\boldsymbol{I}_{{\rm av},h}^{l}\boldsymbol{v}_{h}$ and used a triangle inequality to pass to the third line, we have invoked again (12) to estimate the second term in the right-hand side to pass to the fourth line, and we have used the definition (13) of the strain norm to conclude. ∎

Remark 2* (Korn–Poincaré inequality).*

Combining the discrete Poincaré inequality resulting from [16, Theorem 6.1] (see also [17, Theorem 5.3 and Corollary 5.4]) with (14), we infer that it holds, for all $\boldsymbol{v}_{h}\in\mathbb{P}^{l}(\mathcal{T}_{h};\mathbb{R}^{d})$,

[TABLE]

with hidden constant independent of $h$ and $\boldsymbol{v}_{h}$.

4 Discretisation

4.1 Discrete space

Given a mesh $\mathcal{M}_{h}=(\mathcal{T}_{h},\mathcal{F}_{h})$, we define the following space of discrete unknowns:

[TABLE]

For all $\underline{\boldsymbol{v}}_{h}\in\underline{\boldsymbol{U}}_{h}$, we denote by $\boldsymbol{v}_{h}\in\mathbb{P}^{0}(\mathcal{T}_{h};\mathbb{R}^{d})$ the piecewise constant function obtained by patching element-based unknowns, that is,

[TABLE]

The restrictions of $\underline{\boldsymbol{U}}_{h}$ and $\underline{\boldsymbol{v}}_{h}\in\underline{\boldsymbol{U}}_{h}$ to a generic mesh element $T\in\mathcal{T}_{h}$ are respectively denoted by $\underline{\boldsymbol{U}}_{T}$ and $\underline{\boldsymbol{v}}_{T}=(\boldsymbol{v}_{T},(\boldsymbol{v}_{F})_{F\in\mathcal{F}_{T}})$. The vector of discrete variables corresponding to a smooth function on $\Omega$ is obtained via the global interpolation operator $\underline{\boldsymbol{I}}_{h}:H^{1}(\Omega;\mathbb{R}^{d})\to\underline{\boldsymbol{U}}_{h}$ such that, for all $\boldsymbol{v}\in H^{1}(\Omega;\mathbb{R}^{d})$,

[TABLE]

Its restriction to a generic mesh element $T\in\mathcal{T}_{h}$ is the local interpolator $\underline{\boldsymbol{I}}_{T}:H^{1}(T;\mathbb{R}^{d})\to\underline{\boldsymbol{U}}_{T}$ such that, for all $\boldsymbol{v}\in H^{1}(T;\mathbb{R}^{d})$,

[TABLE]

The displacement is sought in the following subspace of $\underline{\boldsymbol{U}}_{h}$ that strongly incorporates the homogeneous Dirichlet boundary condition:

[TABLE]

4.2 Displacement reconstruction

Let a mesh element $T\in\mathcal{T}_{h}$ be fixed. We define the local displacement reconstruction operator $\boldsymbol{p}_{T}^{1}:\underline{\boldsymbol{U}}_{T}\to\mathbb{P}^{1}(T;\mathbb{R}^{d})$ such that, for all $\underline{\boldsymbol{v}}_{T}\in\underline{\boldsymbol{U}}_{T}$,

[TABLE]

Remark 3* (Explicit expression for the displacement reconstruction operator).*

From (18), one can infer the following explicit expression for the displacement reconstruction operator: For all $\boldsymbol{x}\in T$,

[TABLE]

where $\overline{\boldsymbol{x}}_{T}\coloneq\frac{1}{|T|}\int_{T}\boldsymbol{x}$ denotes the centroid of $T$.

Proposition 4** (Commutation properties for the displacement reconstruction).**

It holds, for all $\boldsymbol{v}\in H^{1}(T;\mathbb{R}^{d})$,

[TABLE]

Proof.

Let $\boldsymbol{v}\in H^{1}(T;\mathbb{R}^{d})$. Recalling the definition (17) of the local interpolator, we have that

[TABLE]

where we have used the definition (18) of the local displacement reconstruction with $\underline{\boldsymbol{v}}_{T}=\underline{\boldsymbol{I}}_{T}\boldsymbol{v}$ in the first line, the definition (5) of the $L^{2}$-orthogonal projector $\boldsymbol{\pi}_{F}^{0}$ along with the fact that $\boldsymbol{\pi}_{T}^{0}\boldsymbol{v}\otimes\boldsymbol{n}_{TF}$ is constant over $F$ to pass to the second line, the Stokes theorem to pass to the third line and the fact that $\boldsymbol{\pi}_{T}^{0}\boldsymbol{v}$ is constant inside $T$ to cancel the second term therein. This proves the first relation in (20). The second relation in (20) immediately follows accounting for the first and recalling the definition (8) of the elliptic projector after observing that the second relation in (18) gives

[TABLE]

To close this section, we define the global displacement reconstruction operator $\boldsymbol{p}_{h}^{1}:\underline{\boldsymbol{U}}_{h}\to\mathbb{P}^{1}(\mathcal{T}_{h};\mathbb{R}^{d})$ obtained by patching the local reconstructions: For all $\underline{\boldsymbol{v}}_{h}\in\underline{\boldsymbol{U}}_{h}$,

[TABLE]

4.3 Discrete bilinear form

We define the bilinear form $\mathrm{a}_{h}:\underline{\boldsymbol{U}}_{h}\times\underline{\boldsymbol{U}}_{h}\to\mathbb{R}$ such that, for all $\underline{\boldsymbol{w}}_{h},\underline{\boldsymbol{v}}_{h}\in\underline{\boldsymbol{U}}_{h}$,

[TABLE]

In the above expression, $\mathrm{j}_{h}:H^{1}(\mathcal{T}_{h};\mathbb{R}^{d})\times H^{1}(\mathcal{T}_{h};\mathbb{R}^{d})\to\mathbb{R}$ is the jump penalisation bilinear form such that, for all $\boldsymbol{w},\boldsymbol{v}\in H^{1}(\mathcal{T}_{h};\mathbb{R}^{d})$,

[TABLE]

while $\mathrm{s}_{h}:\underline{\boldsymbol{U}}_{h}\times\underline{\boldsymbol{U}}_{h}\to\mathbb{R}$ is a stabilisation bilinear form defined from local contributions as follows:

[TABLE]

In the above expression, for all $T\in\mathcal{T}_{h}$ and all $F\in\mathcal{F}_{T}$, we have introduced the boundary difference operator $\boldsymbol{\delta}_{TF}:\underline{\boldsymbol{U}}_{T}\to\mathbb{R}^{d}$ is such that, for any $\underline{\boldsymbol{v}}_{T}\in\underline{\boldsymbol{U}}_{T}$,

[TABLE]

It can be proved that the stabilisation bilinear form enjoys the following consistency property: For all $\boldsymbol{w}\in H^{1}(\Omega;\mathbb{R}^{d})\cap H^{2}(\mathcal{T}_{h};\mathbb{R}^{d})$,

[TABLE]

with hidden constant independent of both $h$ and $\boldsymbol{w}$.

4.4 Comparison with the original HHO method and role of the jump penalisation term

Compared with the original HHO bilinear form defined by [18, Eqs. (24)–(26) and (38)] and written for $k=0$, the bilinear form (22) includes a novel jump penalisation contribution inspired by the discrete Korn inequality of Lemma 1. This term is needed for stability which, for HHO discretisations of the linear elasticity problem, cannot be achieved through local stabilisation terms for $k=0$. As a matter of fact, following the ideas of [21, Section 4.3.1.4], stability would require the use in (23) of a family of local symmetric, positive semidefinite stabilisation bilinear forms $\left\{\mathrm{s}_{T}\;:\;T\in\mathcal{T}_{h}\right\}$ satisfying the following properties:

1. (i)

Local stability and boundedness. For all $T\in\mathcal{T}_{h}$ and all $\underline{\boldsymbol{v}}_{T}\in\underline{\boldsymbol{U}}_{T}$, with hidden constants independent of $h$, $T$, and $\underline{\boldsymbol{v}}_{T}$,

[TABLE] 2. (ii)

Polynomial consistency. For all $\boldsymbol{w}\in\mathbb{P}^{k+1}(T;\mathbb{R}^{d})$,

[TABLE]

Actually, as noticed in [15, Chapter 7], properties (26) and (27) are incompatible. To see it, assume (27), consider a rigid-body motion $\boldsymbol{v}_{{\rm rbm}}$, that is, a function over $\overline{T}$ for which there exist a vector $\boldsymbol{t}_{\boldsymbol{v}}\in\mathbb{R}^{d}$ and a skew-symmetric matrix $\boldsymbol{R}_{\boldsymbol{v}}\in\mathbb{R}^{d\times d}$ such that, for any $\boldsymbol{x}\in\overline{T}$, $\boldsymbol{v}_{{\rm rbm}}(\boldsymbol{x})=\boldsymbol{t}_{\boldsymbol{v}}+\boldsymbol{R}_{\boldsymbol{v}}\boldsymbol{x}$. Take now $\underline{\boldsymbol{v}}_{T}=\underline{\boldsymbol{I}}_{T}\boldsymbol{v}_{{\rm rbm}}$. Since $\boldsymbol{v}_{{\rm rbm}}\in\mathbb{P}^{1}(T;\mathbb{R}^{d})$, the first relation in (20) shows that $\boldsymbol{\nabla}\boldsymbol{p}_{T}^{1}\underline{\boldsymbol{v}}_{T}=\boldsymbol{\pi}_{T}^{0}(\boldsymbol{\nabla}\boldsymbol{v}_{{\rm rbm}})=\boldsymbol{\nabla}\boldsymbol{v}_{{\rm rbm}}=\boldsymbol{R}_{\boldsymbol{v}}$ so that, in particular, $\boldsymbol{\nabla}\boldsymbol{p}_{T}^{1}\underline{\boldsymbol{v}}_{T}$ is skew-symmetric. Hence, $\boldsymbol{\nabla}_{\rm s}\boldsymbol{p}_{T}^{1}\underline{\boldsymbol{v}}_{T}=\boldsymbol{0}$. Moreover, by (27), $\mathrm{s}_{T}(\underline{\boldsymbol{v}}_{T},\underline{\boldsymbol{v}}_{T})=\mathrm{s}_{T}(\underline{\boldsymbol{I}}_{T}^{0}\boldsymbol{v}_{{\rm rbm}},\underline{\boldsymbol{v}}_{T})=0$, again because $\boldsymbol{v}_{\rm rbm}\in\mathbb{P}^{1}(T;\mathbb{R}^{d})$. Hence, the left-hand side of (26) vanishes for all $\underline{\boldsymbol{v}}_{T}=\underline{\boldsymbol{I}}_{T}\boldsymbol{v}_{{\rm rbm}}$ with $\boldsymbol{v}_{{\rm rbm}}$ rigid-body motion. It is, however, easy to construct a rigid-body motion $\boldsymbol{v}_{{\rm rbm}}$ such that the right-hand side does not vanish, which shows that (26) cannot hold. For this reason, the assumption that the discrete unknowns are at least piecewise affine is required in the original HHO method; see [18, Section 4]. Notice that the choice of $\mathrm{s}_{T}$ in (23) retains the polynomial consistency property (27), which is crucial to prove (25).

We next discuss how the stability property modifies for $k=0$. To this end, recalling the definitions (13) of the double-bar strain norm $\|{\cdot}\|_{\boldsymbol{\varepsilon},h}$ and (23) of the stabilisation bilinear form we introduce the triple-bar strain norm such that, for any $\underline{\boldsymbol{v}}_{h}\in\underline{\boldsymbol{U}}_{h}$,

[TABLE]

Lemma 5** (Global stability and boundedness).**

For all $\underline{\boldsymbol{v}}_{h}\in\underline{\boldsymbol{U}}_{h,0}$ it holds

[TABLE]

with hidden constant independent of both $h$ and $\underline{\boldsymbol{v}}_{h}$.

Proof.

It follows from [19, Lemma 4] that

[TABLE]

where $\boldsymbol{\nabla}_{h}:H^{1}(\mathcal{T}_{h};\mathbb{R}^{d})\to L^{2}(\Omega;\mathbb{R}^{d\times d})$ is the broken gradient such that $(\boldsymbol{\nabla}_{h}\boldsymbol{v})_{|T}=\boldsymbol{\nabla}\boldsymbol{v}_{|T}$ for any $T\in\mathcal{T}_{h}$. On the other hand, using the definition of the symmetric gradient for the first bound and Korn’s inequality (14) for the second, we can write

[TABLE]

Combining (31) with (30) yields the result. ∎

4.5 Discrete problem

The low-order scheme for the approximation of problem (1) reads: Find $\underline{\boldsymbol{u}}_{h}\in\underline{\boldsymbol{U}}_{h,0}$ such that

[TABLE]

Using the coercivity of the bilinear form $\mathrm{a}_{h}$ proved in Lemma 7 below together with the discrete Korn inequality (15), we infer that the discrete problem is well-posed and the a priori estimate $|\kern-1.07639pt|\kern-1.07639pt|\underline{\boldsymbol{u}}_{h}|\kern-1.07639pt|\kern-1.07639pt|_{\boldsymbol{\varepsilon},h}\lesssim\alpha^{-\frac{1}{2}}\|\boldsymbol{f}\|$ holds for the discrete solution, with hidden constant independent of both $h$ and of the problem data, and triple-bar strain seminorm defined by (28) below.

Remark 6* (Static condensation for problem (32)).*

The jump stabilisation introduces a direct link among discrete unknowns attached to neighbouring mesh elements. As a result, static condensation of element-based unknowns no longer appears to be an interesting option.

5 Convergence analysis

In this section, after studying the properties of the discrete bilinear form $\mathrm{a}_{h}$, we prove a priori estimates for the error in the energy- and $L^{2}$-norms.

5.1 Properties of the discrete bilinear form

Lemma 7** (Properties of $\mathrm{a}_{h}$).**

The bilinear form $\mathrm{a}_{h}$ enjoys the following properties:

1. (i)

Stability and boundedness.* Recalling the definition (28) of the triple-bar strain norm and the bound (3) on Lamé’s coefficients, for all $\underline{\boldsymbol{v}}_{h}\in\underline{\boldsymbol{U}}_{h}$ it holds*

[TABLE]

where the hidden constants are independent of both $h$ and the problem data. 2. (ii)

Consistency.* It holds for all $\boldsymbol{w}\in H_{0}^{1}(\Omega;\mathbb{R}^{d})\cap H^{2}(\mathcal{T}_{h};\mathbb{R}^{d})$ such that $\boldsymbol{\nabla}{\cdot}\boldsymbol{\sigma}(\boldsymbol{\nabla}_{\rm s}\boldsymbol{w})\in L^{2}(\Omega;\mathbb{R}^{d})$,*

[TABLE]

where the hidden constant is independent of $\boldsymbol{w}$, $h$, and of the problem data, the linear form $\mathcal{E}_{h}({\boldsymbol{w}};{{\cdot}}):\underline{\boldsymbol{U}}_{h,0}\to\mathbb{R}$ representing the consistency error is such that, for all $\underline{\boldsymbol{v}}_{h}\in\underline{\boldsymbol{U}}_{h,0}$,

[TABLE]

and its dual norm is given by

[TABLE]

Proof.

(i) Stability and boundedness. Let $\underline{\boldsymbol{v}}_{h}\in\underline{\boldsymbol{U}}_{h}$. We recall the Frobenius product such that, for all $\boldsymbol{\tau},\boldsymbol{\eta}\in\mathbb{R}^{d\times d}$, $\boldsymbol{\tau}{:}\boldsymbol{\eta}\coloneq\sum_{i=1}^{d}\sum_{j=1}^{d}\tau_{ij}\eta_{ij}$ with corresponding norm $\|\boldsymbol{\tau}\|_{\rm F}\coloneq(\boldsymbol{\tau}{:}\boldsymbol{\tau})^{\frac{1}{2}}$. Writing (22) for $\underline{\boldsymbol{w}}_{h}=\underline{\boldsymbol{v}}_{h}$, using the assumption (3) on Lamé’s parameters to infer that $\boldsymbol{\sigma}(\boldsymbol{\tau}){:}\boldsymbol{\tau}\geq\alpha\|\boldsymbol{\tau}\|_{\rm F}^{2}$ for any $\boldsymbol{\tau}\in\mathbb{R}_{\rm sym}^{d\times d}$, recalling the definitions (13) and (28) of the double- and triple-bar strain norms, and observing that $2\mu\geq\alpha$, the first inequality in (33) follows. The second inequality can be obtained in a similar way: write (22) for $\underline{\boldsymbol{w}}_{h}=\underline{\boldsymbol{v}}_{h}$, observe that $|\boldsymbol{\sigma}(\boldsymbol{\tau}){:}\boldsymbol{\tau}|\leq(2\mu+d|\lambda|)\|\boldsymbol{\tau}\|_{\rm F}^{2}$ for any $\boldsymbol{\tau}\in\mathbb{R}_{\rm sym}^{d\times d}$, and use again (13) and (28).

(ii) Consistency. Let $\underline{\boldsymbol{v}}_{h}\in\underline{\boldsymbol{U}}_{h,0}$. We reformulate the components of the consistency error. Integrating by parts element by element, we infer that

[TABLE]

where we have used the continuity of normal tractions across interfaces together with the fact that boundary unknowns are set to zero in $\underline{\boldsymbol{U}}_{h,0}$ to insert $\boldsymbol{v}_{F}$ into the right-hand side. To reformulate the second term in (35), in the expression (22) of $\mathrm{a}_{h}$ we use the first property in (20) together with the linearity of the strain-stress law $\boldsymbol{\sigma}$ to write, for all $T\in\mathcal{T}_{h}$, $\boldsymbol{\sigma}(\boldsymbol{\nabla}_{\rm s}\boldsymbol{p}_{T}^{1}\underline{\boldsymbol{I}}_{T}\boldsymbol{w})=\boldsymbol{\sigma}(\boldsymbol{\pi}_{T}^{0}(\boldsymbol{\nabla}_{\rm s}\boldsymbol{w}))=\boldsymbol{\pi}_{T}^{0}(\boldsymbol{\sigma}(\boldsymbol{\nabla}_{\rm s}\boldsymbol{w}))$ and obtain

[TABLE]

After expanding, for all $T\in\mathcal{T}_{h}$, $\boldsymbol{\nabla}_{\rm s}\boldsymbol{p}_{T}^{1}\underline{\boldsymbol{v}}_{T}$ according to its definition (18), we deduce that

[TABLE]

Plugging the above relations into the expression (35) of the consistency error, passing to absolute values, using a generalised Hölder inequality with exponents $(2,\infty,2)$ along with $\|\boldsymbol{n}_{TF}\|_{L^{\infty}(F;\mathbb{R}^{d})}\leq 1$ and $h_{F}\leq h_{T}$ for the first term in the right-hand side, and Cauchy–Schwarz inequalities for the remaining terms, we get

[TABLE]

where we have used the second inequality in (29) together with the definition (28) of the triple-bar strain norm to conclude. Recalling the expression (2) of the strain-stress law, we get, for any $T\in\mathcal{T}_{h}$,

[TABLE]

where we have used the approximation properties (6) of the $L^{2}$-orthogonal projector along with $h_{T}\leq h$ to conclude. Using the above estimate, we infer for the first term

[TABLE]

Moving to the second term, we start by observing that

[TABLE]

where we have used, in this order, the second relation in (20), the definition (13) of the jump seminorm, the fact that the jumps of $\boldsymbol{w}$ vanish across any $F\in\mathcal{F}_{h}$, the definition (10) of the jump operator together with the triangle inequality, and the relation

[TABLE]

to exchange the sums over elements and faces. Hence, using the approximation properties (9) of the elliptic projector, $h_{T}\leq h$, and taking the square root, we arrive at

[TABLE]

For the third term, (25) readily gives

[TABLE]

Plugging (38), (40), and (41) into (36) and passing to the supremum yields (34). ∎

5.2 Energy error estimate

Theorem 8** (Energy error estimate).**

Let $\boldsymbol{u}\in H_{0}^{1}(\Omega;\mathbb{R}^{d})$ denote the unique solution to (4), for which we assume the additional regularity $\boldsymbol{u}\in H^{2}(\mathcal{T}_{h};\mathbb{R}^{d})$. For all $h\in{\cal H}$, let $\underline{\boldsymbol{u}}_{h}\in\underline{\boldsymbol{U}}_{h,0}$ denote the unique solution to (32). Then,

[TABLE]

with hidden constant independent of $h$, $\boldsymbol{u}$, and of the problem data.

Proof.

Applying to the present setting the results of [14, Theorem 10] gives the abstract estimate

[TABLE]

Using the assumed regularity for the exact solution to invoke (34), (42) follows. ∎

Remark 9* (Robustness in the quasi-incompressible limit).*

In the numerical approximation of linear elasticity problems, a key point consists in devising schemes that are robust in the quasi incompressible limit corresponding to $\frac{\lambda}{2\mu}\gg 1$ (which requires, in particular $\lambda\geq 0$). From a mathematical perspective, this property is expressed by the fact that the error estimates are uniform in $\lambda$. For $d=2$ and $\Omega$ convex, it is proved, e.g., in [8, Lemma 2.2] that

[TABLE]

with hidden constant possibly depending on $\Omega$ and $\mu$ but independent of $\lambda$. This result can be extended to $d=3$ reasoning as in the above reference and accounting for the regularity estimates for the Stokes problem derived in [2, Theorem 3]. Plugging (43) into (42) and observing that, when $\lambda\geq 0$, we can take $\alpha=2\mu$ (cf. (3)), we can write, with hidden constant independent of both $h$ and $\lambda$,

[TABLE]

which shows that our error estimate (42) is uniform in $\lambda$. The key point to obtain robustness is the first commutation property in (20), which is used to estimate the term $\mathcal{T}_{1}$ in the proof of Lemma 7.

Remark 10* (Quasi-optimality of the error estimate).*

It follows from the second inequality in (33) that the bilinear form $\mathrm{a}_{h}$ is bounded with boundedness constant independent of $h$. Hence, following [14, Remark 11], the error estimate (42) is quasi-optimal.

Remark 11* (Energy estimate in the $\|{\cdot}\|_{\mathrm{a},h}$-norm for $\lambda\geq 0$).*

When $\lambda\geq 0$, a consistency estimate in $h$ holds for $\|\mathcal{E}_{h}({\boldsymbol{w}};{\cdot})\|_{\mathrm{a},h,*}$, the norm of the consistency error linear form dual to $\|{\cdot}\|_{\mathrm{a},h}$ (see (33)). To see it, observe that, from (36) together with $(2\mu)^{\frac{1}{2}}|\kern-1.07639pt|\kern-1.07639pt|\underline{\boldsymbol{v}}_{h}|\kern-1.07639pt|\kern-1.07639pt|_{\boldsymbol{\varepsilon},h}\leq\|\underline{\boldsymbol{v}}_{h}\|_{\mathrm{a},h}$ (a consequence of the assumption $\lambda\geq 0$), it follows $\left|\mathcal{E}_{h}({\boldsymbol{w}};{\underline{\boldsymbol{v}}_{h}})\right|\lesssim\left(\mathfrak{T}_{1}+\mathfrak{T}_{2}+\mathfrak{T}_{3}\right)(2\mu)^{-\frac{1}{2}}\|\underline{\boldsymbol{v}}_{h}\|_{\mathrm{a},h}$. Hence, passing to the supremum over \big{\{}\underline{\boldsymbol{v}}_{h}\in\underline{\boldsymbol{U}}_{h,0}\;:\;\|\underline{\boldsymbol{v}}_{h}\|_{\mathrm{a},h}=1\big{\}}, we infer

[TABLE]

Invoking again [14, Theorem 10], this time with $\underline{\boldsymbol{U}}_{h,0}$ equipped with the $\|{\cdot}\|_{\mathrm{a},h}$-norm, it is inferred

[TABLE]

with hidden constant having the same dependencies as in (42).

5.3 Improved $L^{2}$-error estimate

Combining the discrete Korn–Poincaré inequality (15) with the error estimate (42), we can infer an estimate in $h$ for the $L^{2}$-norm of the displacement error $\|\boldsymbol{u}_{h}-\boldsymbol{\pi}_{h}^{0}\boldsymbol{u}\|$, where we remind the reader that $\boldsymbol{u}_{h}$ is defined according to (16) as the broken polynomial obtained patching element unknowns, while $\boldsymbol{\pi}_{h}^{0}\boldsymbol{u}$ is the $L^{2}$-orthogonal projection of the exact solution on $\mathbb{P}^{0}(\mathcal{T}_{h};\mathbb{R}^{d})$. It is well-known, however, that improved $L^{2}$-error estimates can be derived in the context of HHO methods when elliptic regularity holds. In this section, we show that the same is true for the low-order method considered in this work. For the sake of simplicity, we assume throughout this section that

[TABLE]

This assumption could be removed, but we keep it here to simplify the discussion and point out the robustness in the quasi-incompressible limit. Recalling the discussion in Remark 9, elliptic regularity for our problem entails that, for all $\boldsymbol{g}\in L^{2}(\Omega;\mathbb{R}^{d})$, the unique solution of the (dual) problem: Find $\boldsymbol{z}_{\boldsymbol{g}}\in H_{0}^{1}(\Omega;\mathbb{R}^{d})$ such that

[TABLE]

satisfies the a priori estimate

[TABLE]

with hidden constant only depending on $\Omega$ and $\mu$.

Theorem 12** (Improved $L^{2}$-error estimate).**

Under the assumptions and notations of Theorem 8, and further assuming $\lambda\geq 0$, elliptic regularity, and $\boldsymbol{f}\in H^{1}(\mathcal{T}_{h};\mathbb{R}^{d})$, it holds that

[TABLE]

where the hidden constant is independent of both $h$ and $\lambda$ (but possibly depends on $\mu$).

Proof.

Inside the proof, hidden constants have the same dependencies as in (47). Applying the results of [14, Theorem 13] to the present setting gives the basic estimate

[TABLE]

We proceed to bound the addends in the right-hand side, denoted for the sake of brevity $\mathfrak{T}_{1}$ and $\mathfrak{T}_{2}$.

(i) Estimate of $\mathfrak{T}_{1}$. Since $\boldsymbol{z}_{\boldsymbol{g}}\in H_{0}^{1}(\Omega;\mathbb{R}^{d})\cap H^{2}(\Omega;\mathbb{R}^{d})$, the consistency estimate (34) followed by the elliptic regularity bound (46) yield, for any $\boldsymbol{g}\in L^{2}(\Omega;\mathbb{R}^{d})$,

[TABLE]

Combined with the energy error estimate (44), this yields

[TABLE]

(ii) Estimate of $\mathfrak{T}_{2}$. Recalling the expression (35) of the consistency error, expanding the bilinear form $\mathrm{a}_{h}$ according to its definition (22) with $\underline{\boldsymbol{w}}_{h}=\underline{\boldsymbol{I}}_{h}\boldsymbol{u}$ and $\underline{\boldsymbol{v}}_{h}=\underline{\boldsymbol{I}}_{h}\boldsymbol{z}_{\boldsymbol{g}}$, and invoking (20) to replace $\boldsymbol{p}_{h}^{1}\underline{\boldsymbol{I}}_{h}$ with $\boldsymbol{\varpi}_{h}^{1}$ and $\boldsymbol{\nabla}_{{\rm s},h}\boldsymbol{p}_{h}^{1}\underline{\boldsymbol{I}}_{h}$ with $\boldsymbol{\pi}_{h}^{0}\boldsymbol{\nabla}_{\rm s}$, we can write

[TABLE]

We have that

[TABLE]

where we have used the fact that (1a) holds almost everywhere in $\Omega$ to replace $-\boldsymbol{\nabla}{\cdot}\boldsymbol{\sigma}(\boldsymbol{\nabla}_{\rm s}\boldsymbol{u})$ with $\boldsymbol{f}$ along with the definitions (7) and (5) of the global and local $L^{2}$-orthogonal projectors in the first line, we have added the quantity $(\boldsymbol{f},\boldsymbol{z}_{\boldsymbol{g}})-(\boldsymbol{\sigma}(\boldsymbol{\nabla}_{\rm s}\boldsymbol{u}),\boldsymbol{\nabla}_{\rm s}\boldsymbol{z}_{\boldsymbol{g}})=0$ (see (4)) in the second line, while, to pass to the third line, we have used the fact that, by definition of $\boldsymbol{\pi}_{h}^{0}$, the function $(\boldsymbol{\pi}_{h}^{0}\boldsymbol{f}-\boldsymbol{f})$ is $L^{2}(\Omega;\mathbb{R}^{d})$-orthogonal to $\mathbb{P}^{0}(\mathcal{T}_{h};\mathbb{R}^{d})$ to insert $\boldsymbol{\pi}_{h}^{0}\boldsymbol{z}_{\boldsymbol{g}}$ into the first term. The Cauchy–Schwarz inequality and the approximation property (6) of the $L^{2}$-orthogonal projector inside each mesh element yield for the first term in the right-hand side

[TABLE]

where we have used a standard estimate on $|\boldsymbol{z}_{\boldsymbol{g}}|_{H^{1}(\Omega;\mathbb{R}^{d})}$ obtained letting $\boldsymbol{v}=\boldsymbol{z}_{\boldsymbol{g}}$ in (45) and using the Cauchy–Schwarz and Korn inequalities to bound the right-hand side. On the other hand, using the definitions (7) and (5) of the global and local $L^{2}$-orthogonal projectors, we can write

[TABLE]

where we have used a Cauchy–Schwarz inequality to pass to the third line, (37) with $\boldsymbol{w}=\boldsymbol{u}$ together with the approximation property (6) of the $L^{2}$-orthogonal projector to pass to the fourth line, and the elliptic regularity bound (46) to conclude. Finally, using Cauchy–Schwarz inequalities, we can write

[TABLE]

where, to pass to the second line, we have used (25) for the terms involving $\mathrm{s}_{h}$ and we have proceeded as in the estimate of $\mathfrak{T}_{2}$ in the proof of point (ii) of Lemma 7 for the terms involving $\mathrm{j}_{h}$ while, to conclude, we have invoked (46). Taking absolute values in (50) and using the estimates (51), (52), (53) yields $\left|\mathcal{E}_{h}({\boldsymbol{u}};{\underline{\boldsymbol{I}}_{h}\boldsymbol{z}_{\boldsymbol{g}}})\right|\lesssim h^{2}\|\boldsymbol{f}\|_{H^{1}(\mathcal{T}_{h};\mathbb{R}^{d})}\|\boldsymbol{g}\|$. Hence, passing to the supremum, we obtain

[TABLE]

Plugging (49) and (54) into (48) concludes the proof. ∎

6 Numerical tests

In what follows we verify, through numerical examples, the results stated in the previous section.

6.1 Two-dimensional quasi-incompressible case

The first test case is inspired by [9]: we solve on the unit square $\Omega=(0,1)^{2}$ the homogeneous Dirichlet problem corresponding to the exact solution such that

[TABLE]

The corresponding forcing term is

[TABLE]

We take $\mu=1$ and, in order to assess the robustness of the method in the quasi-incompressible limit, we let $\lambda$ vary in $\{1,10^{3},10^{6}\}$. For the numerical solution, we consider structured and unstructured triangular, Cartesian orthogonal, and deformed quadrangular mesh families; see Figure 1. The solutions corresponding to $\lambda=1$ and $\lambda=10^{6}$ on the finest Cartesian orthogonal mesh are represented in Figure 2, where we have plotted the components of the global displacement reconstruction obtained from the discrete solution according to (21).

The numerical results are collected in Tables 1–4, where the following quantities are monitored: $N_{{\rm dofs},h}$, the number of degrees of freedom; $N_{{\rm nz},h}$, the number of non-zero entries in the problem matrix; $|\kern-1.07639pt|\kern-1.07639pt|\underline{\boldsymbol{u}}_{h}-\underline{\boldsymbol{I}}_{h}\boldsymbol{u}|\kern-1.07639pt|\kern-1.07639pt|_{\mathrm{a},h}$, the energy-norm of the error; $\|\boldsymbol{u}_{h}-\boldsymbol{\pi}_{h}^{0}\boldsymbol{u}\|$, the $L^{2}$-norm of the error estimated in Theorem 12. Notice that, in view of Remark 11, in this and in the following numerical tests the energy error is measured using the $\|{\cdot}\|_{\mathrm{a},h}$-norm, whose computation can be done using the already assembled problem matrix. We additionally display the Estimated Order of Convergence (EOC) which, denoting by $e_{i}$ an error measure on the $i$th mesh refinement with meshsize $h_{i}$, is computed as

[TABLE]

In all the cases, the asymptotic EOC match the ones predicted by the theory, that is, $1$ for the energy-norm of the error and $2$ for the $L^{2}$-norm. The results additionally highlight the robustness of the method in the quasi-incompressible limit (see Remark 9) and with respect to the mesh, showing errors of comparable magnitude irrespectively of the value of $\lambda$ and of the selected mesh family.

6.2 Two-dimensional singular case

We next consider the solution of [1, Section 5.1] which, in polar coordinates $(r,\theta)$, reads

[TABLE]

where the various parameters take the following numerical values: $\mu=0.65$, $\lambda=0.98$, $G=\frac{5}{13}$, $\kappa=\frac{9}{5}$, $L=0.5444837367825$, $Q=0.5430755788367$. The forcing term in this case is equal to zero, while the Dirichlet boundary condition is inferred from the exact solution.

The domain $\Omega$ is illustrated in Figure 4, while the solution on the finest computational mesh considered here is depicted in Figure 4. This test case is representative of real-life situations corresponding to a mode 1 fracture in a plain strain problem. The solution exhibits a singularity in the origin, which prevents the method from attaining the full orders of convergence predicted for smooth solutions.

For the numerical resolution, we consider a sequence of refined structured quadrangular meshes. The numerical results collected in the top half of Table 5 show an asymptotic EOC in the energy-norm of about $0.54$, while the asymptotic EOC in the $L^{2}$-norm is about $1.31$. For the sake of completeness, we show, in the bottom half of Table 5, a comparison with the original HHO method of [18] with $k=1$. Also in this case, the EOC are limited by the regularity of the solution, and coincide with those observed for the method studied in this work. As expected, the number of unknowns on a given mesh is larger for the method of [18] compared to the method proposed here, despite the fact that static condensation is applied in the former case. It has to be noticed, however, that the reduction in the number of unknowns is balanced by the increased number of nonzero entries in the matrix, due to both the absence of static condensation (see Remark 6) and the presence of the jump penalisation term. This phenomenon is specific to the two-dimensional case: in dimension $d=3$, the matrix corresponding to the method of [18] with $k=1$ is generally more dense; see, e.g., Table 6. The errors in the energy norm appear to be smaller for the HHO method of [18], but this is in part due to the fact that the natural energy norm associated with the corresponding bilinear form does not contain the norm of the jumps.

6.3 Three-dimensional compressible case

To test the performance of the method in three space dimensions, we solve on the unit cube domain $\Omega=(0,1)^{3}$ the homogeneous Dirichlet problem corresponding to the exact solution $\boldsymbol{u}=(u_{i})_{1\leq i\leq d}$ such that

[TABLE]

The corresponding forcing term is

[TABLE]

For the numerical solution, we take $\mu=\lambda=1$. Table 6 collects the numerical results on Cartesian orthogonal and unstructured simplicial mesh families. The monitored quantities are the same as for the other test cases to which we add, for the sake of comparison, the number of unknowns and of nonzero matrix entries for the method of [18] with $k=1$. For both mesh families, the asymptotic EOC for both the energy- and the $L^{2}$-norms of the error agree with the ones predicted. Specifically, on the simplicial mesh family an EOC close to 1 is attained starting from the third mesh refinement in the energy norm, whereas an EOC close to 2 is already observed starting from the second mesh refinement; on the Cartesian orthogonal mesh family, on the other hand, the orders of convergence take longer to settle to the corresponding asymptotic values, likely because the first computational meshes are very coarse.

7 Local balances and continuity of numerical tractions

In this section we show that our method satisfies local force balances with equilibrated face tractions. This property can be exploited, e.g., to derive a posteriori error estimates by flux equilibration, and it makes the proposed method suitable for integration into existing Finite Volume codes.

Lemma 13** (Traction formulation of the discrete bilinear form).**

We have the following reformulation of the discrete bilinear form $\mathrm{a}_{h}$ defined by (22): For all $\underline{\boldsymbol{w}}_{h},\underline{\boldsymbol{v}}_{h}\in\underline{\boldsymbol{U}}_{h,0}$,

[TABLE]

where, for all $T\in\mathcal{T}_{h}$ and all $F\in\mathcal{F}_{T}$, we have introduced the numerical traction $\boldsymbol{\Phi}_{TF}:\underline{\boldsymbol{U}}_{T}\to\mathbb{P}^{0}(F;\mathbb{R}^{d})$ such that

[TABLE]

with jump penalisation and stabilisation contributions respectively defined as

[TABLE]

where, for any $X$ mesh element or face, we have denoted by $\overline{\boldsymbol{x}}_{X}\coloneq\frac{1}{|X|}\int_{X}\boldsymbol{x}$ its centroid and, for any $T\in\mathcal{T}_{h}$ and any $F\in\mathcal{F}_{T}$, $\epsilon_{TF}\coloneq\boldsymbol{n}_{TF}{\cdot}\boldsymbol{n}_{F}$ defines the orientation of $F$ relative to $T$.

Proof.

We proceed to reformulate the three terms in the right-hand side of (22) in order to highlight the corresponding contribution to the numerical traction. For the consistency term, we can write

[TABLE]

where we have used the fact that, for any $T\in\mathcal{T}_{h}$, both $\boldsymbol{\sigma}(\boldsymbol{\nabla}_{\rm s}\boldsymbol{p}_{T}^{1}\underline{\boldsymbol{w}}_{T})$ and $\boldsymbol{\nabla}_{\rm s}\boldsymbol{p}_{T}^{1}\underline{\boldsymbol{v}}_{T}$ are constant inside $T$ along with the fact that $\boldsymbol{\sigma}(\boldsymbol{\nabla}_{\rm s}\boldsymbol{p}_{T}^{1}\underline{\boldsymbol{w}}_{T})$ is symmetric to replace $\boldsymbol{\nabla}_{\rm s}$ with $\boldsymbol{\nabla}$ in the first line, the first relation in (18) to pass to the second line, and we have rearranged the products and sums to conclude.

For the jump penalisation term, we can start by observing that

[TABLE]

where we have used the definition of the jump operator to pass to the second line and exchanged the sums over elements and faces according to (39) to conclude. Using the explicit expression (19) of the local displacement reconstruction, we can go on writing

[TABLE]

where, to insert $\boldsymbol{v}_{F}$ into the first term in the second line, we have used the fact that $[\boldsymbol{p}_{h}^{1}\underline{\boldsymbol{w}}_{h}]_{F}$ is single-valued at interfaces together with $\boldsymbol{v}_{F}=\boldsymbol{0}$ on boundary faces, to pass to the third line we have used the fact that the discrete unknowns in $\underline{\boldsymbol{v}}_{h}$ are constant over mesh elements to take them out of the integrals over faces while, to conclude, we have observed that $(\boldsymbol{x}-\overline{\boldsymbol{x}}_{T}){\cdot}\boldsymbol{n}_{TF}=(\overline{\boldsymbol{x}}_{F}-\overline{\boldsymbol{x}}_{T}){\cdot}\boldsymbol{n}_{TF}$ for all $\boldsymbol{x}\in F$ and we have used the definition of $\boldsymbol{\Phi}_{\mathrm{j},TF}(\underline{\boldsymbol{w}}_{h})$ after switching the names of the mute variables $F$ and $G$ in the second term of the third line.

Moving to the stabilisation term, we can write

[TABLE]

where we have used the definition (24) of the boundary difference operator together with the explicit expression (19) of the local displacement reconstruction to pass to the second line, we have rearranged the terms to pass to the third line, and we have used the definition of $\boldsymbol{\Phi}_{\mathrm{s},TF}(\underline{\boldsymbol{w}}_{T})$ after switching the names of the mute variables $F$ and $G$ in the second term of the third line to conclude. ∎

Corollary 14** (Local balances and equilibrated tractions).**

Under the assumptions and notations of Lemma 13, we have that $\underline{\boldsymbol{u}}_{h}\in\underline{\boldsymbol{U}}_{h,0}$ solves the discrete problem (32) if and only if: For all $T\in\mathcal{T}_{h}$ the following balance holds

[TABLE]

and, for any interface $F\in\mathcal{F}_{h}^{{\rm i}}$ shared by the mesh elements $T_{1}$ and $T_{2}$, it holds that

[TABLE]

Proof.

Plugging the flux reformulation (55) of the bilinear form $\mathrm{a}_{h}$ into the discrete problem (32), and recalling (16), we infer that it is equivalent to: Find $\underline{\boldsymbol{u}}_{h}\in\underline{\boldsymbol{U}}_{h,0}$ such that

[TABLE]

Taking, for a given mesh element $T\in\mathcal{T}_{h}$, $\underline{\boldsymbol{v}}_{h}$ such that $\boldsymbol{v}_{T^{\prime}}=\boldsymbol{0}$ for all $T^{\prime}\in\mathcal{T}_{h}\setminus\{T\}$, $\boldsymbol{v}_{F}=\boldsymbol{0}$ for all $F\in\mathcal{F}_{h}$, and letting $\boldsymbol{v}_{T}$ span $\mathbb{P}^{0}(T;\mathbb{R}^{d})$, (58) reduces to (56). Similarly, given an interface $F\in\mathcal{F}_{h}^{{\rm i}}$ shared by the mesh elements $T_{1}$ and $T_{2}$, taking in (58) $\underline{\boldsymbol{v}}_{h}$ such that $\boldsymbol{v}_{T}=\boldsymbol{0}$ for all $T\in\mathcal{T}_{h}$, $\boldsymbol{v}_{F^{\prime}}=\boldsymbol{0}$ for all $F^{\prime}\in\mathcal{F}_{h}\setminus\{F\}$, and letting $\boldsymbol{v}_{F}$ span $\mathbb{P}^{0}(F;\mathbb{R}^{d})$, (58) reduces to (57) after recalling that the numerical tractions are constant over $F$. ∎

Acknowledgements

The work of the second author was partially supported by Agence Nationale de la Recherche grants HHOMM (ANR-15-CE40-0005) and fast4hho (ANR-17-CE23-0019).