# The Johnson-Mercier Elasticity Element in Any Dimensions

Jay Gopalakrishnan  Portland State University, PO Box 751 (MTH), Portland, OR
97207, USA<EMAIL_ADDRESS>,  Johnny Guzmán  Division of Applied Mathematics,
Brown University, Box F, 182 George Street, Providence, RI 02912, USA
<EMAIL_ADDRESS>and  Jeonghun J. Lee  Department of Mathematics,
Baylor University, Sid Richardson Science Building, 1410 S.4th Street, Waco,
TX 76706, USA<EMAIL_ADDRESS>
###### Abstract.

Mixed methods for linear elasticity with strongly symmetric stresses of lowest
order are studied in this paper. On each simplex, the stress space has
piecewise linear components with respect to its Alfeld split (which connects
the vertices to barycenter), generalizing the Johnson-Mercier two-dimensional
element to higher dimensions. Further reductions in the stress space in the
three-dimensional case (to 24 degrees of freedom per tetrahedron) are possible
when the displacement space is reduced to local rigid displacements. Proofs of
optimal error estimates of numerical solutions and improved error estimates
via postprocessing and the duality argument are presented.

This work was supported in part by the National Science Foundation (NSF). The
authors gratefully acknowledge the interaction opportunities provided by ICERM
during their 2024 program.

##  1\. Introduction

The classical Hellinger–Reissner variational principle of linear elasticity
yields saddle point problems simultaneously seeking the stress tensor field
and the displacement vector field. A standard finite element discretization of
the Hellinger–Reissner formulation then needs two finite element spaces,
namely a stress space $\Sigma_{h}$ of symmetric matrix-valued piecewise
polynomials whose row-wise divergence is square integrable, and a displacement
space $V_{h}$ of vector fields with square-integrable components. Decades of
research have illuminated the nontrivial difficulties in finding such a pair
of spaces, $\Sigma_{h}\times V_{h}$, which also satisfy the Babuška–Brezzi
stability conditions of saddle point problems.

To the best of our knowledge, the first such finite element space
$\Sigma_{h}\times V_{h}$ with rigorous mathematical analysis was presented in
[29] for the two-dimensional elastostatic problem (see also [37]). This finite
element on triangular meshes in [29], often called the Johnson–Mercier
element, uses a Clough–Tocher refinement [16] of a given mesh to define the
shape functions of the discrete stress space $\Sigma_{h},$ whereas the
discrete displacement space $V_{h}$ is composed of piecewise linear functions
defined on the original mesh. More precisely, $\Sigma_{h}$ is the space of
symmetric tensors such that divergence of each of its rows is square-
integrable, and each of its components is piecewise linear on the
Clough–Tocher refinement. A family of higher order elasticity elements in two
dimensions with the Clough-Tocher refinement of triangular meshes was
developed in [6], where the displacement spaces $V_{h}$ also have piecewise
polynomials on the Clough–Tocher refinement.

The goal of this work is an investigation of higher dimensional versions of
these works on an Alfeld splits, the natural generalization of Clough-Tocher
split to tetrahedra. At the outset, we were unaware that an extension of the
Johnson–Mercier element to three dimensions (3D) was already done in 1982 by
Křížek in [30]. This is even before Alfeld’s work [1] after which the
tetrahedral split was known as the Alfeld split. A cursory literature search
gives the impression that Křížek’s work is largely unknown at present: it is
not mentioned in influential previous studies like [2, 26, 27, 23, 31, 28].
Nonetheless, a proof of unisolvency of the 3D version of Johnson-Mercier
stress element and error estimates for stress are presented in [30]. In this
paper, we present a different proof of unisolvency of the same 3D stress
element. We show that our new technique of proof generalizes to give a stress
element in higher dimensions. Even in 3D, the new understanding helps us
develop a reduced stress element with only 24 degrees of freedom per element,
which to our knowledge, is the simplest conforming symmetric stress element
known currently. We also prove error estimates for the displacement, prove
superconvergence estimates by duality argument, show robustness in the
incompressible limit, and develop a postprocessing to enhance convergence
rates. We also study the method obtained by pairing the same stress space with
piecewise constant displacements with respect to the Alfeld refinement
(instead of the piecewise linear displacement spaces that Johnson and Mercier
used).

Our approach is based on the analysis of [29, 6], rather than that of [30].
The unisolvency proofs in [29, 6] of stress finite elements $\Sigma_{h}$
proceed by characterizing the divergence-free subspace of $\Sigma_{h}$ that
have vanishing normal components on the edges of a triangular element. They
show that it must be the Airy stress tensor of $C^{1}$ spline on the
Clough–Tocher split. In the more modern language of finite element theory,
this fits in a local discrete complex on Clough–Tocher splits such that the
complex is exact with the corresponding spaces with homogeneous boundary
conditions. Note that ${\mathop{\mathrm{div}\,}}\Sigma_{h}$ is not a subspace
of the displacement space $V_{h}$ used in the Johnson–Mercier element, so one
would have to replace $V_{h}$ with piecewise constant displacements on the
Clough–Tocher triangulation to relate it to the complex. In contrast, the
three-dimensional Johnson–Mercier element in [30] utilizes special geometric
properties of tetrahedra and inter-element continuity of polynomials for
unisolvency proof but this argument seems to be difficult to generalize to
higher dimensions. For Johnson–Mercier element in any dimensions, the key is
to exploit an identity used in the BGG resolution. We give an elementary proof
of this identity is arbitrary dimension. Although we have chosen to present
this generalization without using the language of exterior calculus, it is
indeed motivated by the considerations to make the local stress space fit into
an exact discrete elasticity complex. Indeed, similar ideas were critical in
[14] where the first 3D discrete elasticity complex was developed on Alfeld
splits. Since we only need to focus on the last part of such an
$N$-dimensional complex to obtain the generalization of the Johnson-Mercier
element, a complete study of the full $N$-dimensional complex is outside the
scope of this paper.

Composite stress elements are interesting because they promise liberation from
the necessity of vertex degrees of freedom when using polynomial elements [2].
Due to complications with supersmoothness, this promise was not realized on
Alfeld splits in [14], but was realized on the more involved Worsey-Farin
splits (which splits a tetrahedron to 12 subtetrahedra) later in [23]. Through
our improved understanding of the lowest order case, in this paper we are
finally able to provide a 3D stress element on Alfeld splits without vertex
and edge degrees of freedom. Even in higher dimensions, our composite stress
element does not have degrees of freedom on subsimplices of codimension
greater than one. In particular, this allows for a natural facet-based
hybridization. Another nice property of our elasticity element is that only
low degree piecewise polynomials are used for the stress finite element. Most
symmetric conforming stress elements for 3D elasticity contain piecewise
quadratic or higher degree polynomials for stress (although [23] is an
exception). Only some reduced elements have the stress shape functions which
do not contain all quadratic polynomials [9], but contains other high
polynomial degree shape functions. In this paper, even in $N$-dimensions, only
piecewise linear polynomials are used for the stress.

The paper is organized as follows. We begin by defining the elasticity problem
in Section 2. There we present our unisolvency proof of the 3D Johnson–Mercier
element. We then analyze the stability and the error estimates of numerical
solutions with the Johnson–Mercier element. Optimal error estimates of
numerical solutions and superconvergent error estimates of post-processed
solutions are proved. In Section 3, another numerical method with the 3D
Johnson–Mercier stress element paired with piecewise constant displacement, is
proposed and analyzed. We then present a reduced finite element pair in
Section 4 and prove its error estimates. The connection of the Johnson–Mercier
element and linear elasticity finite element methods with weakly symmetric
stress is discussed in Section 5. Generalization of the Johnson–Mercier
element for any $N$-dimensions ($N\geq 2$) is given in Section 6. Finally, a
further reduced finite element pair and its error estimates are discussed in
Appendix A.

##  2\. The Johnson-Mercier Element for 3D Elasticity

###  2.1. Notation

For a measurable set $D\subset\mathbb{R}^{N}$ $(N\geq 2)$, let $L^{2}(D)$
denote the Lebesgue space with $N$-dimensional Lebesgue measure. For a domain
$D\subset\mathbb{R}^{N}$ with Lipschitz boundary, $H^{k}(D)$, $k\geq 1$,
denotes the standard Sobolev space [19] of functions in $L^{2}(D)$ all of
whose derivatives of up to $k$th order are also in $L^{2}(D)$. Let
$\mathbb{M}$ denote the space of $N\times N$ matrices. For $M\in\mathbb{M}$,
let $\text{sym }M=\frac{1}{2}(M+M^{\prime}),$ $\text{skw
}M=\frac{1}{2}(M-M^{\prime}),$ and let
$\mathbb{S}=\operatorname{sym}\mathbb{M}$ and
$\mathbb{K}=\operatorname{skw}(\mathbb{M})$. For vector spaces $\mathbb{X}$
such as $\mathbb{R},\mathbb{R}^{N},\mathbb{M},\mathbb{S},\mathbb{K},$ etc., we
let $L^{2}(\Omega,\mathbb{X})$ denote the space of $\mathbb{X}$-valued
functions whose components are in $L^{2}(\Omega)$. The inner product of
$L^{2}(D,\mathbb{X})$ and the corresponding norm are denoted
$(\cdot,\cdot)_{D}$ and $\|\cdot\|_{L^{2}(D,\mathbb{X})}$, respectively.
Similarly, the $\mathbb{X}$-valued Sobolev spaces and norms are
$H^{k}(D,\mathbb{X})$ and $\|\cdot\|_{H^{k}(D,\mathbb{X})}$. The $\mathbb{X}$
in the subscripts of (semi)norm notation will be omitted if there is no
concern of confusion. Furthermore, if $D$ is the domain $\Omega$ of our
boundary value problem, then we abbreviate $(\cdot,\cdot)_{\Omega}$ to
$(\cdot,\cdot)$. On matrix-valued functions, the divergence, denoted by
$\operatorname{div}$, is calculated row wise. Doing so in the sense of
distributions, we define

$\displaystyle H(\operatorname{div},D,\mathbb{M}):=\\{\omega\in
L^{2}(D,\mathbb{M}):\operatorname{div}\omega\in L^{2}(D,\mathbb{R}^{N})\\}$

For a set $D$ in $d$-dimensional hyperspace ($d\leq N$) and a nonnegative
integer $k$, $\mathcal{P}_{k}(D)$ is the space of polynomials of degree at
most $k$, restricted to the domain $D$; similarly,
$\mathcal{P}_{k}(D,\mathbb{X})$ denotes the space of $\mathbb{X}$-valued
polynomials of degree $\leq k$.

###  2.2. Elasticity equations

Let $\Omega\subset\mathbb{R}^{3}$ be a domain with a Lipschitz polyhedral
boundary. We denote the spaces of three-dimensional column vectors and the
space of $3\times 3$ symmetric matrices by $\mathbb{R}^{3}$ and $\mathbb{S}$.
We consider the linear elasticity equations with homogeneous Dirichlet
boundary conditions:

(2.1a) $\displaystyle\mathcal{A}\sigma=$ $\displaystyle\epsilon(u)\qquad$
$\displaystyle\text{ in }\Omega,$ (2.1b)
$\displaystyle{\mathop{\mathrm{div}\,}}\sigma=$ $\displaystyle f\qquad$
$\displaystyle\text{ in }\Omega,$ (2.1c) $\displaystyle u=$ $\displaystyle
0\qquad$ $\displaystyle\text{ on }\partial\Omega$

where $u:\Omega\rightarrow\mathbb{R}^{3}$ is the displacement vector field,
$\sigma:\Omega\rightarrow\mathbb{S}$ is the stress tensor field,
$f:\Omega\rightarrow\mathbb{R}^{3}$ is the load vector, and
$\epsilon(u)=\operatorname{sym}(\nabla u)$ is the symmetrized gradient. Here
$\mathcal{A}$ is the compliance tensor, a rank-4 tensor field with entries
$\mathcal{A}_{ijkl}:\Omega\rightarrow\mathbb{R}$, $1\leq i,j,k,l\leq 3$ which
satisfy symmetry conditions

(2.2)
$\displaystyle\mathcal{A}_{ijkl}=\mathcal{A}_{jikl}=\mathcal{A}_{ijlk}=\mathcal{A}_{klij}.$

For $\tau\in\mathbb{S}$ with entries $\tau_{ij}$, $1\leq i,j\leq 3$,
$\mathcal{A}(x)\tau$ for $x\in\Omega$ is the matrix whose $(i,j)$-entry is
defined by $\sum_{1\leq k,l\leq 3}\mathcal{A}_{ijkl}(x)\tau_{kl}$. From (2.2),
one can check that $\mathcal{A}(x)\tau$ also belongs to $\mathbb{S}$, and
$\mathcal{A}(x)$ gives a bilinear form on $\mathbb{S}$ by
$\mathcal{A}(x)\tau:\omega$ for $\tau,\omega\in\mathbb{S}$ where the colon
stands for the Frobenius inner product of matrices. We assume that
$\mathcal{A}$ is positive definite over symmetric matrices and is bounded,
namely there are positive constants $\theta,\gamma>0$ such that

(2.3) $\displaystyle\mathcal{A}(x)\omega:\omega$
$\displaystyle\geq\theta\omega:\omega,\quad\omega\in\mathbb{S},x\in\Omega,$
(2.4) $\displaystyle\|\mathcal{A}\|_{L^{\infty}(\Omega)}$
$\displaystyle\leq\gamma.$

In our results for nearly incompressible elastic materials, we make an
alternate assumption since (2.3) does not generally hold in such cases.
Instead of (2.3), then assume that there exists a $\theta>0$ such that

(2.5) $\displaystyle\mathcal{A}(x)\omega:\omega$
$\displaystyle\geq\theta\,\text{dev}(\omega):\text{dev}(\omega),\quad\omega\in\mathbb{S},x\in\Omega$

where $\text{dev}(\omega):=\omega-\operatorname{tr}(\omega)\mathbb{I}$ with
the $3\times 3$ identity matrix $\mathbb{I}$. For isotropic materials, as the
Lame material parameter $\lambda\to\infty$, we can expect (2.5) to hold, but
not (2.3).

A well-known variational formulation of (2.1) is as follows: Find
$\sigma\in\Sigma$, $u\in V$ such that

(2.6a)
$\displaystyle(\mathcal{A}\sigma,\tau)+(u,{\mathop{\mathrm{div}\,}}\tau)=$
$\displaystyle 0\qquad$ $\displaystyle{\text{ for all }}\tau\in\Sigma,$ (2.6b)
$\displaystyle({\mathop{\mathrm{div}\,}}\sigma,v)=$ $\displaystyle(f,v)\qquad$
$\displaystyle{\text{ for all }}v\in V$

where $\Sigma=H(\operatorname{div},\Omega,\mathbb{S})$,
$V:=[L^{2}(\Omega)]^{3}$.

###  2.3. The finite element

A canonical finite element discretization of (2.6) is to find
$\sigma_{h}\in\Sigma_{h}$, $u_{h}\in V_{h}$ such that

(2.7a)
$\displaystyle(\mathcal{A}\sigma_{h},\tau)+(u_{h},{\mathop{\mathrm{div}\,}}\tau)=$
$\displaystyle 0\qquad$ $\displaystyle{\text{ for all }}\tau\in\Sigma_{h},$
(2.7b) $\displaystyle({\mathop{\mathrm{div}\,}}\sigma_{h},v)=$
$\displaystyle(f,v)\qquad$ $\displaystyle{\text{ for all }}v\in V_{h}$

with suitable finite element spaces $\Sigma_{h}\subset\Sigma$, $V_{h}\subset
V$. It is well-known that the finite element pair $(\Sigma_{h},V_{h})$ needs
to satisfy the Babuška–Brezzi stability conditions for existence and
uniqueness of numerical solutions and for accurate approximation of exact
solutions (see, e.g., [12]).

To construct a finite element subspace $\Sigma_{h}$, we work on Alfeld
simplicial complexes and start by establishing notation associated to an
Alfeld split. Starting with a tetrahedron $T=[x_{0},x_{1},x_{2},x_{3}]$, let
$T^{\text{A}}$ be an Alfeld triangulation of $T$, i.e., we choose an interior
point $z$ of $T$ and we let $T_{0}=[z,x_{1},x_{2},x_{3}]$,
$T_{1}=[z,x_{0},x_{2},x_{3}],$ $T_{2}=[z,x_{0},x_{1},x_{3}],$
$T_{3}=[z,x_{0},x_{1},x_{2}]$ and set
$T^{\text{A}}=\\{T_{0},T_{1},T_{2},T_{3}\\}$.

We let $\mathcal{T}_{h}$ be a conforming triangulation $\Omega$ with
tetrahedra and $\mathcal{T}_{h}^{\text{A}}$ be the resulting mesh after
performing an Alfeld split on each $T\in\mathcal{T}_{h}$. For
$T\in\mathcal{T}_{h}$ we define the local spaces as

$\displaystyle\Sigma_{h}(T):=$ $\displaystyle\\{\omega\in
H({\mathop{\mathrm{div}\,}},T,\mathbb{S}):\omega|_{K}\in\EuScript{P}_{1}(K,\mathbb{S}),\text{
for all }K\in T^{\text{A}}\\}$ (2.8) $\displaystyle V_{h}(T):=$
$\displaystyle[\EuScript{P}_{1}(T)]^{3}.$

Although $\operatorname{div}(\Sigma_{h}(T))$ is not contained in $V_{h}(T)$,
it is contained in

$\displaystyle W_{h}(T):=$
$\displaystyle\\{v\in[L^{2}(T)]^{3}:v\in[\EuScript{P}_{0}(K)]^{3},\text{ for
all }K\in T^{\text{A}}\\}.$

Note that $W_{h}(T)$ and $V_{h}(T)$ have the same dimension. We let $P_{T}$
denote the $L^{2}$-orthogonal projection onto $W_{h}(T)$:

(2.9) $\displaystyle(P_{T}v,w)_{T}=(v,w)_{T}\qquad{\text{ for all }}w\in
W_{h}(T),\;{\text{ for all }}T\in\mathcal{T}_{h}.$

For every $K\in T^{\text{A}}$ let $x_{K}$ be the barycenter of $K$. We let
$I_{T}:W_{h}(T)\mapsto V_{h}(T)$ be defined by

(2.10) $I_{T}w(x_{K})=w(x_{K})\qquad\text{ for }K\in T^{\text{A}}.$

A simple quadrature argument gives that

(2.11) $(I_{T}w,m)_{T}=(w,m)_{T}\qquad{\text{ for all }}w,m\in W_{h}(T),$

which, in particular, shows that $I_{T}$ is an injection. Since $\dim
W_{h}(T)=\dim V_{h}(T),$ injectivity of $I_{T}$ implies bijectivity of
$I_{T}$. Moreover, it is easy to see from (2.9) and (2.11) that

(2.12) $P_{T}I_{T}w=w,\qquad\text{ for all }w\in W_{h}(T),$

so the restriction of $P_{T}$ to $V_{h}(T)$, namely
$P_{T}|_{V_{h}(T)}:V_{h}(T)\to W_{h}(T),$ and $I_{T}:W_{h}(T)\to V_{h}(T)$ are
inverses of each other.

The next lemma is crucial for identifying degrees of freedom of the stress
element $\Sigma_{h}(T)$. First we recall some preliminaries that are standard
(see e.g. [2]). Let $\mathbb{M}$ denote the space of $3\times 3$ matrices and
let $\Xi:\mathbb{M}\to\mathbb{M}$ be defined by

(2.13) $\Xi M=M^{\prime}-\text{tr}(M)\mathbb{I},$

where $(\cdot)^{\prime}$ denotes the transpose. The inverse of this operator
is given by $\Xi^{-1}M=M^{\prime}-\frac{1}{2}\text{tr}(M)\mathbb{I}.$ We
define $\operatorname{mskw}:\mathbb{R}^{3}\to\mathbb{K}$ and vskw by

$\displaystyle\text{mskw}\begin{pmatrix}v_{1}\\\ v_{2}\\\
v_{3}\end{pmatrix}:=\begin{pmatrix}0&-v_{3}&v_{2}\\\ v_{3}&0&-v_{1}\\\
-v_{2}&v_{1}&0\end{pmatrix},\qquad\text{vskw}:=\text{mskw}^{-1}\circ\text{skw
}.$

It is easy to see that

(2.14) $\displaystyle{\mathop{\mathrm{div}\,}}\Xi=2\text{vskw
}{\mathop{\mathrm{curl}\,}}$

holds. Let $\EuScript{P}_{p}(T^{\text{A}})=\\{u:u|_{K}\in\EuScript{P}_{p}(K)$
for all $K\in T^{\text{A}}\\}$ and let
$\mathring{\mathcal{L}}_{p}(T^{\text{A}})=\EuScript{P}_{p}(T^{\text{A}})\cap\mathring{H}^{1}(T)$.
Using this notation, we prove the next lemma showing that the space of
divergence-free elements with vanishing traces of $\Sigma_{h}(T)$ is trivial.
Our proof of this lemma marks the main point of departure from the analysis of
[30]. An alternative to the approach of [30], is to prove this lemma using
local discrete elasticity sequences on Alfeld splits obtained by small
modifications of the recent three-dimensional exact sequences in [14].
Nonetheless, we have chosen to present a new proof that is easier and
generalizable to higher dimensions (as we shall demonstrate in Section 6).
Throughout, $n$ denotes the unit outward normal vector field on the boundary
of a domain under consideration.

######  Lemma 2.1.

The equality
$\\{\omega\in\Sigma_{h}(T):{\mathop{\mathrm{div}\,}}\omega=0,\omega
n|_{\partial T}=0\\}=\\{0\\}$ holds.

###### Proof.

Consider the matrix analogue of the well-known lowest order Brezzi-Douglas-
Marini (BDM) space, namely ${B\\!D\\!M}_{1}(T^{\text{A}})=\\{\omega\in
H(\operatorname{div},T,\mathbb{M}):\omega|_{K}\in\EuScript{P}_{1}(K,\mathbb{M})\text{
for all }K\in T^{\text{A}}\\}$ and let
$\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}})=\\{\omega\in{B\\!D\\!M}_{1}(T^{\text{A}}):\omega
n=0\text{ on }\partial T\\}$. Note that

$\displaystyle\\{\omega\in\Sigma_{h}(T):{\mathop{\mathrm{div}\,}}\omega=0,\omega
n|_{\partial
T}=0\\}=\\{\omega\in\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}}):\;\operatorname{vskw}\omega=0,\;{\mathop{\mathrm{div}\,}}\omega=0\\}.$

Let
$\tilde{\EuScript}{P}_{0}(T^{\text{A}})=\\{u\in\EuScript{P}_{0}(T^{\text{A}}):\int_{T}u=0\\}$
and let

$X=\left\\{(u,v):u\in[\EuScript{P}_{1}(T^{\text{A}})]^{3},\;v\in[\tilde{\EuScript}{P}_{0}(T^{\text{A}})]^{3},\;\int_{T}(u\cdot
c+v\cdot(c\times x))=0\text{ for all }c\in\mathbb{R}^{3}\right\\}.$

The proof is based on the operator
$A:\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}})\to X$ given by

(2.15)
$A\eta=(\operatorname{vskw}\eta,{\mathop{\mathrm{div}\,}}\eta),\qquad\eta\in\mathring{{B\\!D\\!M}}_{1}.$

Observe that for any $c\in\mathbb{R}^{3}$ and any
$\eta\in\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}})$,

(2.16)
$\int_{T}c\cdot\operatorname{vskw}\eta=\int_{T}\eta:\operatorname{mskw}c=\int_{T}\eta:{\mathop{\mathrm{grad}\,}}(c\times
x)=-\int_{T}\operatorname{div}\eta\cdot(c\times x).$

Thus $A$ indeed maps into $X$. We proceed to show that $A$ is surjective.

Let $(u,v)\in X$. Since $v$ has components of zero mean, by a standard exact
sequence property (e.g. [7]), there exists a
$\sigma\in\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}})$ such that

(2.17) ${\mathop{\mathrm{div}\,}}\sigma=v.$

Since $(u,v)\in X$, this implies

$0=\int_{T}(u\cdot c+v\cdot(c\times x))=\int_{T}(u\cdot
c+\operatorname{div}\sigma\cdot(c\times
x))=\int_{T}(u-\operatorname{vskw}\sigma)\cdot c$

for all $c\in\mathbb{R}^{3}$, where we have again used (2.16) to obtain the
last equality. Hence by a well-known exact sequence property ([38, 25, 22]) on
Alfeld splits, there exists an
$\eta\in\mathring{\mathcal{L}}_{2}(T^{\text{A}},\mathbb{M})$ such that

(2.18) ${\mathop{\mathrm{div}\,}}\eta=2(u-\operatorname{vskw}\sigma).$

Now $\tau=\sigma+{\mathop{\mathrm{curl}\,}}\Xi^{-1}\eta$ is obviously in
$\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}})$ and satisfies

$\displaystyle\operatorname{vskw}\tau$
$\displaystyle=\operatorname{vskw}\sigma+\frac{1}{2}{\mathop{\mathrm{div}\,}}\eta=u,$
(by (2.14) and (2.18)) $\displaystyle{\mathop{\mathrm{div}\,}}\tau$
$\displaystyle=v.$ (by (2.17))

Thus we have proved that $A$ is surjective.

To conclude, the surjectivity of $A$ implies that $\operatorname{rank}(A)=\dim
X\geq
3\dim\EuScript{P}_{1}(T^{\text{A}})+3\dim\tilde{\EuScript}{P}_{0}(T^{\text{A}})-3=54$.
Since $\dim\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}})=3\cdot 6\cdot 3=54$ on an
Alfeld split (see e.g., [14, eq. (2.3)]), the rank-nullity theorem implies
that the null space of $A$ is trivial. It means that
$\\{\omega\in\mathring{{B\\!D\\!M}}_{1}:\;\operatorname{vskw}\omega=0,\;{\mathop{\mathrm{div}\,}}\omega=0\\}=\\{0\\}$,
so the conclusion follows. ∎

We can now state unisolvent degrees of freedom (“dofs”) for the space
$\Sigma_{h}(T)$. Let $\triangle_{k}(T)$ denote the set of $k$-subsimplices of
a simplex $T$, and let $\triangle_{k}(\mathcal{T}_{h})$ be the union of
$\triangle_{k}(T)$ for all $T\in\mathcal{T}_{h}$.

######  Theorem 2.2.

The dimension of $\Sigma_{h}(T)$ is 42. Moreover, an element
$\omega\in\Sigma_{h}(T)$ is uniquely determined by the following dofs:

(2.19a) $\displaystyle\int_{F}\omega n\cdot\kappa,\qquad$
$\displaystyle\kappa\in[\EuScript{P}_{1}(F)]^{3},\;F\in\triangle_{2}(T),\qquad$
(36 dofs) (2.19b) $\displaystyle\int_{T}\omega.$ ($6$ dofs)

###### Proof.

Since
$\Sigma_{h}(T)=\\{\omega\in{B\\!D\\!M}_{1}(T^{\text{A}}):\int_{T}\omega:\eta=0\;\text{
for all }\eta\in\EuScript{P}_{1}(T^{\text{A}},\mathbb{K})\\}$,

$\displaystyle\dim\Sigma_{h}(T)$
$\displaystyle\geq\dim{B\\!D\\!M}_{1}(T^{\text{A}})-\dim\EuScript{P}_{1}(T^{\text{A}},\mathbb{K})$
$\displaystyle=90-48=42,$

which is exactly the number of dofs given in (2.19). We now show that
$\dim\Sigma_{h}(T)\leq 42$.

Let $\omega\in\Sigma_{h}(T)$ for which all dofs of (2.19) vanish. We will show
that $\omega=0$. The dofs (2.19a) give that $\omega n=0$ on $\partial T$. Let
$v\in\EuScript{P}_{1}(T)$, then integration by parts gives

(2.20) $\int_{T}\operatorname{div}\omega\cdot v=-\int_{T}\omega:\epsilon(v)=0$

where we used (2.19b) in the last equality. Note that
$\operatorname{div}\omega\in W_{h}(T)$. Hence, substituting
$v=I_{T}\operatorname{div}\omega$ in (2.20) and using (2.11) we get
$0=\int_{T}v\,\operatorname{div}\omega=\int_{T}|\operatorname{div}\omega|^{2}$
which shows that $\operatorname{div}\omega=0$. By Lemma 2.1, $\omega$
vanishes. This proves that $\dim\Sigma_{h}(T)\leq 42$, so
$\dim\Sigma_{h}(T)=42.$ The argument also shows that the dofs (2.19) uniquely
determine an element $\omega\in\Sigma_{h}(T)$. ∎

Unisolvent dofs of a finite element generate an associated canonical
interpolant. Let $\Pi_{T}$ denote this canonical projection into
$\Sigma_{h}(T)$. It satisfies

(2.21a) $\displaystyle\int_{F}\Pi_{T}\omega n\,\cdot\kappa=\int_{F}\omega
n\cdot\kappa,\qquad$
$\displaystyle\kappa\in[\EuScript{P}_{1}(F)]^{3},\;F\in\triangle_{2}(T),$
(2.21b) $\displaystyle\int_{T}\Pi_{T}\omega=\int_{T}\omega$

for any $\omega$ in $D_{\Pi_{T}}:=\\{\omega\in
H(\operatorname{div},T,\mathbb{S}):\;\omega n|_{F}\in[L^{2}(F)]^{3},\text{ for
all }F\in\triangle_{2}(T)\\}$. It follows that

(2.22) $\int_{T}\operatorname{div}(\Pi_{T}\omega-\omega)\cdot v=0\quad{\text{
for all }}v\in V_{h}(T),\;\omega\in D_{\Pi_{T}}$

by integrating by parts and (2.21).

The global discrete spaces are given by

(2.23) $\displaystyle\Sigma_{h}:=$ $\displaystyle\\{\omega\in
H(\operatorname{div},\Omega,\mathbb{S}):\omega|_{T}\in\Sigma_{h}(T),\text{ for
all }T\in\mathcal{T}_{h}\\},$ (2.24) $\displaystyle V_{h}:=$
$\displaystyle\\{v\in[L^{2}(\Omega)]^{3}:v|_{T}\in V_{h}(T),\text{ for all
}T\in\mathcal{T}_{h}\\},$ (2.25) $\displaystyle W_{h}:=$
$\displaystyle\\{w\in[L^{2}(\Omega)]^{3}:w|_{T}\in W_{h}(T),\text{ for all
}T\in\mathcal{T}_{h}\\}.$

We also define the global interpolant $\Pi:D_{\Pi}\mapsto\Sigma_{h}$ by
$\Pi\omega|_{T}=\Pi_{T}\omega$ for all $T\in\mathcal{T}_{h}$ where

(2.26) $\displaystyle D_{\Pi}:=\\{\omega\in
H(\operatorname{div},\Omega,\mathbb{S}):\omega n|_{F}\in[L^{2}(F)]^{3},{\text{
for all }}F\in\triangle_{2}(\mathcal{T}_{h})\\}.$

Note that $H^{1}(\Omega,\mathbb{S})\subset D_{\Pi}$. Finally, we need the
global versions of $P_{T}$ and $I_{T}$: the projection $P$ onto $W_{h}$ and
$I:W_{h}\mapsto V_{h}$ are defined by $Pv|_{T}=P_{T}v$ and $Iw|_{T}=I_{T}w$
for all $T\in\mathcal{T}_{h}$. The next result is an easy consequence of
(2.22).

######  Proposition 2.3.

For all $v\in V_{h}$ and $\omega\in D_{\Pi}$,

(2.27) $(\operatorname{div}(\Pi\omega-\omega),v)=0.$

Moreover, for $s=1,2$, for all $\omega\in H^{s}({\Omega},\mathbb{S})$,

(2.28)
$\displaystyle\|\omega-\Pi\omega\|_{L^{2}(\Omega)}+h\|\operatorname{div}(\omega-\Pi\omega)\|_{L^{2}(\Omega)}$
$\displaystyle\leq Ch^{s}\|\omega\|_{H^{s}(\Omega)}.$

The next result follows easily by combining a well-known technique of proving
wellposedness [20] with a more recent result on regular potentials [17,
Corollary 4.7] that show that given any $z\in H^{\ell}({\Omega})$, for any
real $\ell$, there exists a $u\in[H^{\ell+1}({\Omega})]^{3}$ such that

(2.29) $\operatorname{div}u=z,\qquad\|u\|_{H^{\ell+1}({\Omega})}\leq
C\|z\|_{H^{\ell}({\Omega})}.$

Hereon, $C$ will denote a generic mesh-independent constant whose value at
different occurrences may vary.

######  Lemma 2.4.

Given any $v\in[L^{2}({\Omega})]^{3}$, there exists a $\tau\in
H^{1}({\Omega},\mathbb{S})$ such that

(2.30) $\operatorname{div}\tau=v,\qquad\|\tau\|_{H^{1}(\Omega,\mathbb{S})}\leq
C\|v\|_{L^{2}(\Omega)}.$

###### Proof.

Using (2.29) with $\ell=0$ for each component of $v$, we find a $\tau_{1}\in
H^{1}({\Omega},\mathbb{M})$ such that

(2.31)
$\operatorname{div}\tau_{1}=v,\qquad\|\tau_{1}\|_{H^{1}({\Omega},\mathbb{M})}\leq
C\|v\|_{L^{2}({\Omega})}.$

Next, we apply (2.29) with $\ell=1$ and $z=-2\operatorname{vskw}\tau_{1}$ to
find a $\rho$ such that
$\operatorname{div}(\Xi\rho)=-2\operatorname{vskw}\tau_{1}$. Since
$\operatorname{vskw}\tau_{1}\in[H^{1}({\Omega})]^{3}$, and $\Xi$ is a
pointwise algebraic bijection, the resulting $\rho$ is in
$H^{2}({\Omega},\mathbb{M})$ and

(2.32) $\|\rho\|_{H^{2}({\Omega})}\leq C\|\tau_{1}\|_{H^{1}({\Omega})}.$

Then putting $\tau=\tau_{1}+{\mathop{\mathrm{curl}\,}}\rho$ we find, using
(2.14), that $\operatorname{div}\tau=v$ and

$2\operatorname{vskw}\tau=2\operatorname{vskw}\tau_{1}+2\operatorname{vskw}({\mathop{\mathrm{curl}\,}}\rho)=2\operatorname{vskw}\tau_{1}+\operatorname{div}(\Xi\rho)=0.$

Combining (2.31) and (2.32), we conclude that
$\tau=\tau_{1}+{\mathop{\mathrm{curl}\,}}\rho$ is in
$H^{1}(\Omega,\mathbb{S})$, satisfies $\operatorname{div}\tau=v$, as well as
the estimate in (2.30). ∎

######  Theorem 2.5 (Discrete stability).

There exists a constant $\beta>0$ such that

(2.33)
$\beta\|v\|_{L^{2}({\Omega})}\leq\sup_{0\neq\omega\in\Sigma_{h}}\frac{(\operatorname{div}\omega,v)}{\|\omega\|_{H(\operatorname{div},\Omega)}}\qquad{\text{
for all }}v\in V_{h}.$

###### Proof.

Given a $v\in V_{h}$, we apply Lemma 2.4 to get a $\tau$ in
$H^{1}({\Omega},\mathbb{S})$ satisfying (2.30). Since
$H^{1}({\Omega},\mathbb{S})\subset D_{\Pi}$, the function $\omega=\Pi\tau$
exists in $\Sigma_{h}.$ Then we have

(2.34)
$\|v\|_{L^{2}(\Omega)}^{2}=(\operatorname{div}\tau,v)=(\operatorname{div}\omega,v)$

by (2.27). The result now follows after using (2.28) with $s=1$ and (2.30)
which gives $\|\omega\|_{H(\operatorname{div},\Omega)}\leq
C\|\tau\|_{H^{1}(\Omega)}\leq C\|v\|_{L^{2}(\Omega)}$. ∎

The mixed method with Johnson–Mercier element finds $\sigma_{h}\in\Sigma_{h}$
and $u_{h}\in V_{h}$ satisfying

(2.35a)
$\displaystyle(\mathcal{A}\sigma_{h},\tau)+(u_{h},\operatorname{div}\tau)=$
$\displaystyle 0\qquad$ $\displaystyle{\text{ for all }}\tau\in\Sigma_{h},$
(2.35b) $\displaystyle(\operatorname{div}\sigma_{h},v)=$
$\displaystyle(f,v)\qquad$ $\displaystyle{\text{ for all }}v\in V_{h}.$

###  2.4. Error analysis

We prove error estimates for (2.35). Let
$\|\tau\|_{\mathcal{A}}=(\mathcal{A}\tau,\tau)^{1/2}$.

######  Theorem 2.6.

Let $\sigma,u$ solve (2.6) and $\sigma_{h},u_{h}$ solve (2.35). Then, we have

(2.36)
$\|\Pi\sigma-\sigma_{h}\|_{\mathcal{A}}\leq\|\Pi\sigma-\sigma\|_{\mathcal{A}}$

and

(2.37) $\|u-u_{h}\|_{L^{2}(\Omega)}\leq
C\big{(}\|\sigma-\sigma_{h}\|_{\mathcal{A}}+\|Qu-u\|_{L^{2}(\Omega)}\big{)}$

where $Q$ is the $L^{2}$-orthogonal projection onto $V_{h}$.

###### Proof.

Using (2.27), (2.35b) and (2.6b) we get that

$(\operatorname{div}(\Pi\sigma-\sigma_{h}),v)=0\quad{\text{ for all }}v\in
V_{h}.$

However,
$(\operatorname{div}(\Pi\sigma-\sigma_{h}),Pv)=(\operatorname{div}(\Pi\sigma-\sigma_{h}),v)$.
Since $P:V_{h}\rightarrow W_{h}$ is an isomorphism, this implies that

(2.38) $\operatorname{div}(\Pi\sigma-\sigma_{h})=0.$

Then, (2.36) follows from (2.35a) and (2.6a).

Next, using Lemma 2.4, we obtain a $\tau\in H^{1}(\Omega,\mathbb{S})$
satisfying

(2.39) $\displaystyle\operatorname{div}\tau$ $\displaystyle=Qu-u_{h},\qquad$
(2.40) $\displaystyle\|\tau\|_{H^{1}(\Omega)}$ $\displaystyle\leq C\|Qu-
u_{h}\|_{L^{2}(\Omega)}.$

Combining (2.39) with (2.27),

$\displaystyle\|Qu-u_{h}\|_{L^{2}(\Omega)}^{2}$ $\displaystyle=(Qu-u_{h},Qu-
u_{h})=(u-u_{h},Qu-u_{h})$
$\displaystyle=(u-u_{h},\operatorname{div}\Pi\tau)+(u-u_{h},\operatorname{div}(\tau-\Pi\tau))$
(by (2.39))
$\displaystyle=-(\mathcal{A}(\sigma-\sigma_{h}),\Pi\tau)+(u-u_{h},\operatorname{div}(\tau-\Pi\tau))$
(by (2.6a), (2.35a))
$\displaystyle=-(\mathcal{A}(\sigma-\sigma_{h}),\Pi\tau)+(u-Qu,\operatorname{div}(\tau-\Pi\tau))$
(by (2.27)).

By the Cauchy–Schwarz inequality, (2.28), and (2.40),

$\|Qu-u_{h}\|_{L^{2}(\Omega)}\leq
C\big{(}\|\sigma-\sigma_{h}\|_{\mathcal{A}}+\|Qu-u\|_{L^{2}(\Omega)}\big{)},$

so (2.37) now follows by the triangle inequality. ∎

Using the above result along with (2.28) and the approximation properties of
$Q$ we obtain the following corollary.

######  Corollary 2.7.

Under the hypothesis of Theorem 2.6, assuming also that $\sigma\in
H^{2}(\Omega,\mathbb{S})$ and $u\in[H^{2}(\Omega)]^{3}$, we have

(2.41) $\|\sigma-\sigma_{h}\|_{\mathcal{A}}+\|u-u_{h}\|_{L^{2}(\Omega)}\leq
Ch^{2}(\|\sigma\|_{H^{2}(\Omega)}+\|u\|_{H^{2}(\Omega)}).$

###  2.5. Superconvergence and postprocessing

In this subsection, we prove that a projection of the error $u-u_{h}$
superconverges at a rate $h^{3}$ if we assume full elliptic regularity, i.e.,
if we assume that ${\Omega}$ and $\mathcal{A}$ are such that the inequality

(2.42) $\|u\|_{H^{2}({\Omega})}+\|\sigma\|_{H^{1}({\Omega})}\leq
C\|f\|_{L^{2}({\Omega})}$

holds for any solution $\sigma,u$ of (2.1) obtained with an $f\in
L^{2}({\Omega})$. We will need the following commuting property:

(2.43) $\operatorname{div}(\Pi\omega-\omega)=0,\qquad\omega\in
H^{1}(\Omega,\mathbb{S}),\;\operatorname{div}\omega\in W_{h}.$

Indeed, for such $\omega$, since $\operatorname{div}(\Pi\omega-\omega)\in
W_{h}$, for all $w\in W_{h},$ we have
$(\operatorname{div}(\Pi\omega-\omega),w)=(\operatorname{div}(\Pi\omega-\omega),Iw)$
by (2.11), and the latter vanishes by (2.27).

######  Theorem 2.8.

Let $\sigma,u$ solve (2.6) and $\sigma_{h},u_{h}$ solve (2.35). Assume that
the full regularity estimate (2.42) holds. Then, the following estimate holds

(2.44) $\|P(u-u_{h})\|_{L^{2}(\Omega)}\leq
Ch\|\sigma-\sigma_{h}\|_{\mathcal{A}}+Ch^{2}\|\operatorname{div}(\sigma-\Pi\sigma)\|_{L^{2}(\Omega)}.$

If $\sigma\in H^{2}({\Omega})$, this implies

(2.45) $\|P(u-u_{h})\|_{L^{2}(\Omega)}\leq Ch^{3}|\sigma|_{H^{2}({\Omega})}.$

###### Proof.

Let $\tau$ and $\psi$ solve

(2.46a) $\displaystyle\mathcal{A}\tau=$ $\displaystyle\epsilon(\psi)\qquad$
$\displaystyle\text{ in }\Omega,$ (2.46b)
$\displaystyle\operatorname{div}\tau=$ $\displaystyle P(u-u_{h})\qquad$
$\displaystyle\text{ in }\Omega,$ (2.46c) $\displaystyle\psi=$ $\displaystyle
0\qquad$ $\displaystyle\text{ on }\partial\Omega.$

By the full elliptic regularity assumption (2.42), we have

(2.47) $\|\tau\|_{H^{1}(\Omega)}+\|\psi\|_{H^{2}(\Omega)}\leq
C\|P(u-u_{h})\|_{L^{2}(\Omega)}.$

Then,

$\displaystyle\|P(u-u_{h})\|_{L^{2}(\Omega)}^{2}$
$\displaystyle=(P(u-u_{h}),\operatorname{div}\tau)\qquad$ (by (2.46b))
$\displaystyle=(P(u-u_{h}),{\mathop{\mathrm{div}\,}}\Pi\tau)\qquad$ (by
(2.43)) $\displaystyle=(u-u_{h},{\mathop{\mathrm{div}\,}}\Pi\tau)\qquad$
$\displaystyle=-(\mathcal{A}(\sigma-\sigma_{h}),\Pi\tau)\qquad$ (by (2.35a),
(2.6a))
$\displaystyle=-(\mathcal{A}(\sigma-\sigma_{h}),\Pi\tau-\tau)-(\mathcal{A}(\sigma-\sigma_{h}),\tau).$

We can now simplify the last term

$\displaystyle(\mathcal{A}(\sigma-\sigma_{h}),\tau)$
$\displaystyle=(\sigma-\sigma_{h},\mathcal{A}\tau)\qquad$ (by symmetry of
$\mathcal{A}$) $\displaystyle=(\sigma-\sigma_{h},\epsilon(\psi))\qquad$ (by
(2.46a))
$\displaystyle=-({\mathop{\mathrm{div}\,}}(\sigma-\sigma_{h}),\psi)\qquad$
(integration by parts)
$\displaystyle=-({\mathop{\mathrm{div}\,}}(\sigma-\Pi\sigma),\psi-Q\psi)\qquad$
(by (2.38), (2.27)).

Hence, we obtain

$\|P(u-u_{h})\|_{L^{2}(\Omega)}^{2}=-(\mathcal{A}(\sigma-\sigma_{h}),\Pi\tau-\tau)+({\mathop{\mathrm{div}\,}}(\sigma-\Pi\sigma),\psi-Q\psi).$

The first estimate (2.44) of the theorem now follows after using (2.28),
approximation properties of $Q$, and (2.47). The second follows from the first
and (2.28). ∎

We can use the superconvergence result to construct a better approximation to
$u$ using a local post-processing scheme. Let $\mathcal{R}(T)$ denote the
space of rigid displacements on $T$. Then, we can decompose $W_{h}(T)$ as

(2.48) $W_{h}(T)=P_{T}\mathcal{R}(T)+(P_{T}\mathcal{R}(T))^{\perp}.$

where

$\displaystyle P_{T}\mathcal{R}(T)$
$\displaystyle:=\\{P_{T}v:v\in\mathcal{R}(T)\\}=\\{w\in
W_{h}(T):I_{T}w\in\mathcal{R}(T)\\},$
$\displaystyle(P_{T}\mathcal{R}(T))^{\perp}$ $\displaystyle:=\\{w\in
W_{h}(T):(w,P_{T}v)_{T}=0,\text{ for all }v\in\mathcal{R}(T)\\}.$

Define the projection $P_{T}^{\mathcal{R}}:L^{2}(T)\to P_{T}\mathcal{R}(T)$ by

(2.49) $(P_{T}^{\mathcal{R}}q,w)_{T}=(q,w)_{T},\quad{\text{ for all }}w\in
P_{T}\mathcal{R}(T).$

Recalling that $P_{T}$ is injective on $V_{h}(T)$ due to (2.12), we have

(2.50a) $\displaystyle\|v\|_{L^{2}(T)}$ $\displaystyle\leq
C\|P_{T}v\|_{L^{2}(T)}\qquad$ $\displaystyle{\text{ for all }}v\in V_{h}(T),$
by equivalence of norms on finite dimensional spaces and scaling. Similarly,
(2.50b) $\displaystyle\|q\|_{L^{2}(T)}$ $\displaystyle\leq
Ch_{T}\|\epsilon(q)\|_{L^{2}(T)}+\|P_{T}^{\mathcal{R}}q\|_{L^{2}(T)}\qquad$
$\displaystyle{\text{ for all }}q\in[\EuScript{P}_{2}(T)]^{3}.$

Indeed, to show that the right hand side of (2.50b) is a norm on
$[\EuScript{P}_{2}(T)]^{3},$ suppose $\epsilon(q)=0$ on $T$. Then $q$ must be
a rigid displacement. If furthermore $P_{T}^{\mathcal{R}}q=0,$ then we see
that $P_{T}q=0$ and (2.50a) shows that $q=0$. Together with standard scaling
arguments, this establishes (2.50b). Define the quadratic space

(2.51) $\displaystyle
S_{h}(T):=\\{q\in[\EuScript{P}_{2}(T)]^{3}:(q,v)_{T}=0,{\text{ for all
}}v\in\mathcal{R}(T)\\}.$

Using it, and following Stenberg [36], we define our post-processed
approximation next.

The post-processing system finds $u^{\star}_{h}\in[\EuScript{P}_{2}(T)]^{3}$
satisfying

(2.52a) $\displaystyle(\epsilon(u^{\star}_{h}),\epsilon(v))_{T}=$
$\displaystyle(\mathcal{A}\sigma_{h},\epsilon(v))_{T}\qquad$
$\displaystyle{\text{ for all }}v\in S_{h}(T),$ (2.52b)
$\displaystyle(u^{\star}_{h},w)_{T}=$ $\displaystyle(Pu_{h},w)_{T}\qquad$
$\displaystyle{\text{ for all }}w\in P_{T}\mathcal{R}(T),$

for every $T\in\mathcal{T}_{h}$, given $\sigma_{h},u_{h}$ solving (2.35). It
is obvious from the definition (2.51) of $S_{h}(T)$ that $\dim
S_{h}(T)+\dim\mathcal{R}(T)=\dim[\EuScript{P}_{2}(T)]^{3}$. Since

(2.53) $\dim\mathcal{R}(T)=\dim P_{T}\mathcal{R}(T)$

by the injectivity of $P_{T}$ on $V_{h}(T)$, the system (2.52) is square.
Moreover, (2.52) is uniquely solvable: indeed if its right hand side vanishes,
then, (2.52a) shows that $u^{\star}_{h}$ is a rigid displacement, which then
together with (2.52b) with (2.50b) implies that $u^{\star}_{h}$ must vanish.

Next, we show that $u-u_{h}^{\star}$ converges at a rate higher than that
expected for $u-u_{h}$ from Corollary 2.7. Let $Q_{2}$ denote the element-wise
$L^{2}$-orthogonal projection onto $[\EuScript{P}_{2}(T)]^{3}$.

######  Theorem 2.9.

Let $\sigma,u$ solve (2.6) and $\sigma_{h},u_{h}$ solve (2.35). If
$u^{\star}_{h}$ is given by (2.52), then

$\displaystyle\|u-u^{\star}_{h}\|_{L^{2}(T)}$ $\displaystyle\leq
C(\|u-Q_{2}u\|_{L^{2}(T)}+\|P_{T}(u_{h}-u)\|_{L^{2}(T)})$
$\displaystyle+Ch_{T}\big{(}\|\epsilon(u-Q_{2}u)\|_{L^{2}(T)}+\|\mathcal{A}(\sigma-\sigma_{h})\|_{L^{2}(T)}\big{)}.$

###### Proof.

Using (2.52a) and (2.6a) we obtain

(2.54) $\displaystyle(\epsilon(u^{\star}_{h}-u),\epsilon(v))_{T}=$
$\displaystyle(\mathcal{A}(\sigma_{h}-\sigma),\epsilon(v))_{T}\qquad\text{ for
all }v\in S_{h}(T).$

Since the same equality trivially holds for all $v\in\mathcal{R}(T)$, it
actually holds for all $v\in[\EuScript{P}_{2}(T)]^{3}$. Hence we may choose
$v=u_{h}^{\star}-Q_{2}u$ in $[\EuScript{P}_{2}(T)]^{3}$ and manipulate to get

(2.55)
$\|\epsilon(u^{\star}_{h}-Q_{2}u)\|_{L^{2}(T)}\leq\|\epsilon(u-Q_{2}u)\|_{L^{2}(T)}+\|\mathcal{A}(\sigma-\sigma_{h})\|_{L^{2}(T)}.$

Using (2.50b) we have

(2.56) $\|u^{\star}_{h}-Q_{2}u\|_{L^{2}(T)}\leq
Ch_{T}\|\epsilon(u^{\star}_{h}-Q_{2}u)\|_{L^{2}(T)}+C\|P_{T}^{\mathcal{R}}(u^{\star}_{h}-Q_{2}u)\|_{L^{2}(T)}.$

It is easy to see that
$P_{T}^{\mathcal{R}}Q_{2}u=P_{T}^{\mathcal{R}}P_{T}Q_{2}u$. Also, (2.52b)
implies $P_{T}^{\mathcal{R}}u_{h}^{\star}=P_{T}^{\mathcal{R}}P_{T}u_{h}$.
Therefore,

$\displaystyle\|P_{T}^{\mathcal{R}}(u^{\star}_{h}-Q_{2}u)\|_{L^{2}(T)}$
$\displaystyle\leq\|P_{T}(u_{h}-Q_{2}u)\|_{L^{2}(T)}$ (2.57)
$\displaystyle\leq\|P_{T}(u_{h}-u)\|_{L^{2}(T)}+\|u-Q_{2}u\|_{L^{2}(T)}.$

The result now follows by using the estimates of (2.55) and (2.57) within
(2.56). ∎

######  Corollary 2.10.

Under the hypotheses of Theorems 2.8 and 2.9, assuming also that $\sigma\in
H^{2}(\Omega,\mathbb{S}),u\in[H^{3}(\Omega)]^{3}$, we have

$\displaystyle\|u-u^{\star}_{h}\|_{L^{2}(\Omega)}\leq
Ch^{3}(\|\sigma\|_{H^{2}(\Omega)}+\|u\|_{H^{3}(\Omega)}).$

###### Proof.

Apply (2.44), (2.28), and a standard estimate for projection error of $Q_{2}$,
to further bound the estimate for $\|u-u^{\star}_{h}\|_{L^{2}(\Omega)}$ given
by Theorem 2.9. ∎

###  2.6. Robustness in the incompressible limit

Finally, we investigate the convergence of this method for nearly
incompressible isotropic materials. Recall that $\mathcal{A}$ remains bounded
in this case but may become close to singular, so methods using
$\mathcal{A}^{-1}$ can become problematic. One of the advantages of the mixed
method is that it only needs $\mathcal{A}$. Nonetheless, as
$\|\cdot\|_{\mathcal{A}}$-norm becomes weaker in the incompressible limit, it
is natural to ask if $\sigma_{h}$ converges in a material-independent norm, a
question not answered by Corollary 2.7. To answer this, we now use techniques
similar to those in [6] to prove that the $L^{2}$ error
$\|\sigma-\sigma_{h}\|_{L^{2}(\Omega)}$ has an error bound that is robust for
nearly incompressible materials satisfying the weaker ellipticity inequality
(2.5). In particular, for isotropic materials with $\lambda\to\infty$, the $C$
below does not depend on $\lambda$.

######  Theorem 2.11.

Assume that only (2.5) holds (instead of (2.3)). Let $\sigma,u$ solve (2.6)
and $\sigma_{h},u_{h}$ solve (2.35). Then,

(2.58) $\displaystyle\|\sigma-\sigma_{h}\|_{L^{2}(\Omega)}\leq
C(\|\sigma-\sigma_{h}\|_{\mathcal{A}}+h\|\operatorname{div}(\sigma-\Pi\sigma)\|_{L^{2}(\Omega)}).$

###### Proof.

It is evident by setting $\tau_{h}$ to the $3\times 3$ identity matrix
$\mathbb{I}$ in (2.6) and (2.35) that
$\int_{\Omega}\operatorname{tr}(\sigma-{\sigma}_{h})=0.$ Hence there exists
$\phi\in[H_{0}^{1}(\Omega)]^{3}$ (see, e.g., [17]) such that

(2.59)
${\mathop{\mathrm{div}\,}}\phi=\operatorname{tr}(\sigma-\sigma_{h}),\qquad\|\phi\|_{H^{1}(\Omega)}\leq
C\|\operatorname{tr}(\sigma-\sigma_{h})\|_{L^{2}(\Omega)}.$

Using this,

(2.60)
$\displaystyle\|\operatorname{tr}(\sigma-{\sigma}_{h})\|_{L^{2}(\Omega)}^{2}$
$\displaystyle=(\operatorname{tr}(\sigma-{\sigma}_{h}),\operatorname{div}\phi)=(\operatorname{tr}(\sigma-{\sigma}_{h})\mathbb{I},\epsilon(\phi))$
$\displaystyle=3\big{(}\sigma-{\sigma}_{h},\epsilon(\phi)\big{)}-3(\operatorname{dev}(\sigma-{\sigma}_{h}),\epsilon(\phi)\big{)}$
$\displaystyle=3({\mathop{\mathrm{div}\,}}(\sigma-{\sigma}_{h}),\phi)-3(\operatorname{dev}(\sigma-{\sigma}_{h}),\epsilon(\phi)).$
$\displaystyle=3({\mathop{\mathrm{div}\,}}(\sigma-\Pi\sigma),\phi)-3(\operatorname{dev}(\sigma-{\sigma}_{h}),\epsilon(\phi)),$

where we used (2.38) in the last step. To estimate the first term on the final
right hand side above, we again use (2.38) (and also (2.27), which is implied
by (2.38)) to get

(2.61) $\displaystyle({\mathop{\mathrm{div}\,}}(\sigma-{\sigma}_{h}),\phi)$
$\displaystyle=({\mathop{\mathrm{div}\,}}(\sigma-\Pi\sigma),\phi)=({\mathop{\mathrm{div}\,}}(\sigma-\Pi\sigma),\phi-Q\phi)$
$\displaystyle\leq
Ch\|{\mathop{\mathrm{div}\,}}(\sigma-\Pi\sigma)\|_{L^{2}(\Omega)}\|\phi\|_{H^{1}(\Omega)}$

The last term in (2.60) is easily estimated by (2.5):

$\displaystyle(\operatorname{dev}(\sigma-{\sigma}_{h}),\epsilon(\phi))$
$\displaystyle\leq
C\|\sigma-{\sigma}_{h}\|_{\mathcal{A}}\|\phi\|_{H^{1}(\Omega)}.$

Together with (2.59), these estimates imply that
$\|\operatorname{tr}(\sigma-{\sigma}_{h})\|_{L^{2}(\Omega)}\leq
Ch\|\sigma\|_{H^{1}(\Omega)}$, so the conclusion follows from
$C\|\sigma-\sigma_{h}\|_{L^{2}(\Omega)}\leq\|\operatorname{dev}(\sigma-\sigma_{h})\|_{L^{2}(\Omega)}+\|\operatorname{tr}(\sigma-\sigma_{h})\|_{L^{2}(\Omega)}$
and (2.5). ∎

######  Remark 2.12.

The above proof extends to other cases (as we shall see in later sections) as
long as an analogue of the critical ingredient (2.38) is available. Also note
that to obtain convergence rates, it is enough to simply use (2.28) in (2.58).

##  3\. Piecewise constant displacements on the refinement

In this section we show that we can replace the displacement space $V_{h}$
with $W_{h}$ and obtain very similar results. This replacement is attractive
since $\operatorname{div}\Sigma_{h}=W_{h}$ and the equilibrium equation can
then be exactly satisfied for numerical solutions. Also, the stress and
displacement elements are implemented with respect to the same refined mesh.
However, a price to pay, as we shall see now, is that we are not able to prove
an $O(h^{3})$ superconvergence result (like in (2.45)) unless the right hand
side $f$ belongs to $W_{h}$.

We start by proving unisolvency of slightly different dofs (cf. (2.19)) for
the same stress space $\Sigma_{h}$. It is useful to note at the outset that by
(2.48) and (2.53),

$\dim(P_{T}\mathcal{R}(T))^{\perp}=12-6,$

so the number of dofs in (2.19b) and (3.1b) both equal six.

######  Theorem 3.1.

An element $\omega\in\Sigma_{h}(T)$ is uniquely determined by the following
dofs:

(3.1a) $\displaystyle\int_{F}\omega n\cdot\kappa,\qquad$
$\displaystyle\kappa\in[\EuScript{P}_{1}(F)]^{3},F\in\triangle_{2}(T),\qquad$
(36 dofs) (3.1b) $\displaystyle\int_{T}{\mathop{\mathrm{div}\,}}\omega\cdot
v,\qquad$ $\displaystyle v\in(P_{T}\mathcal{R}(T))^{\perp}.\qquad$ ($6$ dofs)

###### Proof.

Let $\omega\in\Sigma_{h}(T)$ and assume that the dofs in (3.1) vanish. Then we
must show that $\omega=0$. The dofs in (3.1a) show that $\omega n=0$ on
$\partial T$. Then, by integration by parts, we get
${\mathop{\mathrm{div}\,}}\omega\in(P_{T}\mathcal{R}(T))^{\perp}$, and (3.1b)
shows that ${\mathop{\mathrm{div}\,}}\omega=0$. Then, by Lemma 2.1,
$\omega=0$. ∎

The mixed method with piecewise constant displacement finds
$\hat{\sigma}_{h}\in\Sigma_{h}$ and $\hat{u}_{h}\in W_{h}$ satisfying

(3.2a)
$\displaystyle(\mathcal{A}\hat{\sigma}_{h},\tau)+(\hat{u}_{h},{\mathop{\mathrm{div}\,}}\tau)=$
$\displaystyle 0\qquad$ $\displaystyle{\text{ for all }}\tau\in\Sigma_{h},$
(3.2b) $\displaystyle({\mathop{\mathrm{div}\,}}\hat{\sigma}_{h},v)=$
$\displaystyle(f,v)\qquad$ $\displaystyle{\text{ for all }}v\in W_{h}.$

The next result is analogous to Theorem 2.5, so we omit its proof.

######  Theorem 3.2.

There exists a constant $\beta>0$ such that

(3.3) $\beta\leq\inf_{0\not=v\in
W_{h}}\sup_{0\neq\omega\in\Sigma_{h}}\frac{({\mathop{\mathrm{div}\,}}\omega,v)}{\|\omega\|_{H({\mathop{\mathrm{div}\,}},\Omega)}\|v\|_{L^{2}(\Omega)}}.$

Our error analysis of the method (3.2) does not proceed as in the previous
section, but rather through the following auxiliary problem to find
$\tilde{\sigma}\in\Sigma$, $\tilde{u}\in V$ such that

(3.4a)
$\displaystyle(\mathcal{A}\tilde{\sigma},\tau)+(\tilde{u},{\mathop{\mathrm{div}\,}}\tau)=$
$\displaystyle 0\qquad$ $\displaystyle{\text{ for all }}\tau\in\Sigma,$ (3.4b)
$\displaystyle({\mathop{\mathrm{div}\,}}\tilde{\sigma},v)=$
$\displaystyle(Pf,v)\qquad$ $\displaystyle{\text{ for all }}v\in V.$

Now note that ${\mathop{\mathrm{div}\,}}\tilde{\sigma}=Pf\in W_{h}$.

######  Proposition 3.3.

Let $\sigma,u$ solve (2.6) and $\tilde{\sigma},\tilde{u}$ solve (3.4). Then

(3.5)
$\|\sigma-\tilde{\sigma}\|_{\mathcal{A}}+\|u-\tilde{u}\|_{L^{2}(\Omega)}\leq
Ch\|f-Pf\|_{L^{2}(\Omega)}.$

###### Proof.

Equations (2.6) and (3.4) imply

(3.6)
$\displaystyle\|\sigma-\tilde{\sigma}\|_{\mathcal{A}}^{2}=(f-Pf,u-\tilde{u})=(f-Pf,(u-\tilde{u})-P(u-\tilde{u})).$

However, using element-wise Poincare inequality (see, e.g., [10]) and then
using Korn’s first inequality ([15, Theorem 6.3-4]),

$\|(u-\tilde{u})-P(u-\tilde{u})\|_{L^{2}(\Omega)}\leq
Ch\|\epsilon(u-\tilde{u})\|_{L^{2}(\Omega)}=Ch\|\mathcal{A}(\sigma-\tilde{\sigma})\|_{L^{2}(\Omega)}.$

Using this in (3.6) after an application of Cauchy-Schwarz inequality, the
stated bound for $\|\sigma-\tilde{\sigma}\|_{\mathcal{A}}$ follows. The same
bound for $\|u-\tilde{u}\|_{L^{2}({\Omega})}$ now follows from (3.4a), which
implies that
$(u-\tilde{u},\operatorname{div}\tau)=(\mathcal{A}(\tilde{\sigma}-\sigma),\tau),$
after choosing a $\tau$ such that $\operatorname{div}\tau=u-\tilde{u}$ by
Lemma 2.4. ∎

######  Remark 3.4.

At this point, we may choose to proceed using the canonical interpolant $J$
corresponding to the dofs in (3.1). One can prove an analogue of (2.28) for
$J$ and also that for any $\omega\in D_{\Pi}$ such that
${\mathop{\mathrm{div}\,}}\omega\in W_{h}$, we have
$({\mathop{\mathrm{div}\,}}(J\omega-\omega),w)=0$ for all $w\in W_{h},$ which
is enough for an error analysis using $J$. But we pursue the quicker
alternative of using the same $\Pi$ from Section 2.

######  Theorem 3.5.

Let $\sigma,u$ solve (2.6) and $\hat{\sigma}_{h},\hat{u}_{h}$ solve (3.2).
Then,

(3.7)
$\|\Pi\sigma-\hat{\sigma}_{h}\|_{\mathcal{A}}\leq\|\Pi\tilde{\sigma}-\tilde{\sigma}\|_{\mathcal{A}}+\|\Pi(\tilde{\sigma}-\sigma)\|_{\mathcal{A}}$

and

(3.8)
$C\|u-\hat{u}_{h}\|_{L^{2}(\Omega)}\leq\|\sigma-\hat{\sigma}_{h}\|_{\mathcal{A}}+\|Pu-u\|_{L^{2}(\Omega)}.$

###### Proof.

Subtracting (3.2) from (3.4),

(3.9)
$\displaystyle(\mathcal{A}(\tilde{\sigma}-\hat{\sigma}_{h}),\tau)+(\tilde{u}-\hat{u}_{h},{\mathop{\mathrm{div}\,}}\tau)=$
$\displaystyle 0\qquad$ $\displaystyle{\text{ for all }}\tau\in\Sigma_{h},$
(3.10)
$\displaystyle({\mathop{\mathrm{div}\,}}(\tilde{\sigma}-\hat{\sigma}_{h}),v)=$
$\displaystyle 0\qquad$ $\displaystyle{\text{ for all }}v\in W_{h},$

Since ${\mathop{\mathrm{div}\,}}(\tilde{\sigma}-\hat{\sigma}_{h})$ is in
$W_{h}$, equation (3.10) implies

(3.11) ${\mathop{\mathrm{div}\,}}(\tilde{\sigma}-\hat{\sigma}_{h})=0.$

Also, using (2.43),

(3.12) ${\mathop{\mathrm{div}\,}}(\Pi\tilde{\sigma}-\tilde{\sigma})=0.$

Adding (3.11) and (3.12), we have

(3.13) ${\mathop{\mathrm{div}\,}}(\Pi\tilde{\sigma}-\hat{\sigma}_{h})=0.$

Hence, setting $\tau=\Pi\tilde{\sigma}-\hat{\sigma}_{h}$ in (3.9), we
immediately obtain

$\|\Pi\tilde{\sigma}-\hat{\sigma}_{h}\|_{\mathcal{A}}\leq\|\Pi\tilde{\sigma}-\tilde{\sigma}\|_{\mathcal{A}}.$

The estimate (3.8) follows the same argument as the proof of (2.37). ∎

######  Corollary 3.6.

Let $\sigma,u$ solve (2.6), $\tilde{\sigma},\tilde{u}$ solve (3.4) and
$\hat{\sigma}_{h},\hat{u}_{h}$ solve (3.2). If $\sigma,\tilde{\sigma}\in
H^{2}(\Omega,\mathbb{S}),u\in[H^{1}(\Omega)]^{3}$, then

$\|\sigma-\hat{\sigma}_{h}\|_{\mathcal{A}}+h\|u-\hat{u}_{h}\|_{L^{2}(\Omega)}\leq
Ch^{2}(\|\sigma\|_{H^{2}(\Omega)}+\|\tilde{\sigma}\|_{H^{2}(\Omega)}+\|f\|_{H^{1}(\Omega)}+\|u\|_{H^{1}(\Omega)}).$

The next intermediate result gives an $O(h^{3})$ superconvergence result for
$P(\tilde{u}-\hat{u}_{h})$. It will help us prove superconvergence of the
projection of the error $P(u-\hat{u}_{h})$.

######  Lemma 3.7.

Let $\tilde{\sigma},\tilde{u}$ solve (3.4) and $\hat{\sigma}_{h},\hat{u}_{h}$
solve (3.2). Assume that the full regularity estimate (2.42) holds. Then, the
following estimate holds

(3.14) $\|P(\tilde{u}-\hat{u}_{h})\|_{L^{2}(\Omega)}\leq
Ch\|\tilde{\sigma}-\hat{\sigma}_{h}\|_{\mathcal{A}}.$

###### Proof.

Let $\tau$ and $\psi$ solve

(3.15a) $\displaystyle\mathcal{A}\tau$ $\displaystyle=\epsilon(\psi)\qquad$
$\displaystyle\text{ in }\Omega,$ (3.15b)
$\displaystyle{\mathop{\mathrm{div}\,}}\tau$
$\displaystyle=P(\tilde{u}-\hat{u}_{h})\qquad$ $\displaystyle\text{ in
}\Omega,$ (3.15c) $\displaystyle\psi$ $\displaystyle=0\qquad$
$\displaystyle\text{ on }\partial\Omega.$

Since we are assuming full $H^{2}$-elliptic regularity we have:

(3.16) $\|\tau\|_{H^{1}(\Omega)}+\|\psi\|_{H^{2}(\Omega)}\leq
C\|P(\tilde{u}-\hat{u}_{h})\|_{L^{2}(\Omega)}.$

Then,

$\displaystyle\|P(\tilde{u}-\hat{u}_{h})\|_{L^{2}(\Omega)}^{2}=$
$\displaystyle(P(\tilde{u}-\hat{u}_{h}),{\mathop{\mathrm{div}\,}}\tau)\qquad$
(by (3.15b))
$\displaystyle=(P(\tilde{u}-\hat{u}_{h}),{\mathop{\mathrm{div}\,}}\Pi\tau)\qquad$
(by (2.43))
$\displaystyle=(\tilde{u}-\hat{u}_{h},{\mathop{\mathrm{div}\,}}\Pi\tau)$
$\displaystyle=-(\mathcal{A}(\tilde{\sigma}-\hat{\sigma}_{h}),\Pi\tau)\qquad$
(by (3.9))
$\displaystyle=-(\mathcal{A}(\tilde{\sigma}-\hat{\sigma}_{h}),\Pi\tau-\tau)-(\mathcal{A}(\tilde{\sigma}-\hat{\sigma}_{h}),\tau).$

We can now simplify the last term

$\displaystyle(\mathcal{A}(\tilde{\sigma}-\hat{\sigma}_{h}),\tau)$
$\displaystyle=(\tilde{\sigma}-\hat{\sigma}_{h},\mathcal{A}\tau)\qquad$ (by
symmetry of $\mathcal{A}$)
$\displaystyle=-(\tilde{\sigma}-\hat{\sigma}_{h},\epsilon(w))\qquad$ (by
(3.15a))
$\displaystyle=({\mathop{\mathrm{div}\,}}(\tilde{\sigma}-\hat{\sigma}_{h}),w),\qquad$
(by integration by parts)

which vanishes by (3.11). Hence, we obtain

$\|P(\tilde{u}-\hat{u}_{h})\|_{L^{2}(\Omega)}^{2}=-(\mathcal{A}(\tilde{\sigma}-\hat{\sigma}_{h}),\Pi\tau-\tau).$

The result now follows after using (2.28) and (3.16). ∎

######  Corollary 3.8.

Under the same hypothesis of the above theorem we have

(3.17) $\|P(u-\hat{u}_{h})\|_{L^{2}(\Omega)}\leq
Ch^{2}(\|\sigma\|_{H^{2}(\Omega)}+\|\tilde{\sigma}\|_{H^{2}(\Omega)}+\|f\|_{H^{1}(\Omega)}+\|u\|_{H^{1}(\Omega)}).$

If $f\in W_{h}$ then

(3.18) $\|P(u-\hat{u}_{h})\|_{L^{2}(\Omega)}\leq
Ch^{3}(\|\sigma\|_{H^{2}(\Omega)}+\|\tilde{\sigma}\|_{H^{2}(\Omega)}+\|f\|_{H^{1}(\Omega)}+\|u\|_{H^{1}(\Omega)}).$

###### Proof.

The estimate (3.17) follows from (3.14) and Theorem 3.5 and Proposition 3.3.
If $f\in W_{h}$ then $u=\tilde{u}$ and $\sigma=\tilde{\sigma}$ hence the
result follows also from (3.14) and Theorem 3.5 and Proposition 3.3. ∎

######  Theorem 3.9.

Assume that (2.5). Let $\sigma,u$ solve (2.6), and
$\hat{\sigma}_{h},\hat{u}_{h}$ solve (3.2). Then,

$\displaystyle\|\sigma-\hat{\sigma}_{h}\|_{L^{2}(\Omega)}\leq
C(\|\sigma-\hat{\sigma}_{h}\|_{\mathcal{A}}+h\|f-Pf\|_{L^{2}(\Omega)}).$

###### Proof.

The proof of Theorem 2.11 can be repeated with $\sigma$ and $\sigma_{h}$
replaced by $\tilde{\sigma}$ and $\hat{\sigma}_{h}$, respectively, since the
analogue of the critical ingredient there (see Remark 2.12) is available and
given by (3.13). Then we obtain an analogue of (2.58),

(3.19) $\displaystyle C\|\tilde{\sigma}-\hat{\sigma}_{h}\|_{L^{2}(\Omega)}$
$\displaystyle\leq\|\tilde{\sigma}-\hat{\sigma}_{h}\|_{\mathcal{A}}+h\|\operatorname{div}(\tilde{\sigma}-\Pi\tilde{\sigma})\|_{L^{2}(\Omega)}.$

The last term vanishes by (2.43).

To complete the proof, it suffices to estimate
$\|\sigma-\tilde{\sigma}\|_{L^{2}({\Omega})}$. Let
$\phi\in[H_{0}^{1}(\Omega)]^{3}$ be as in (2.59) with $\sigma_{h}$ replaced by
$\tilde{\sigma}$. Then, just as in (2.60), we obtain

$\displaystyle\|\operatorname{tr}(\sigma-\tilde{\sigma})\|_{L^{2}({\Omega})}$
$\displaystyle=3({\mathop{\mathrm{div}\,}}(\sigma-\tilde{\sigma}),\phi)-3(\operatorname{dev}(\sigma-\tilde{\sigma}),\epsilon(\phi)).$

Now,
$({\mathop{\mathrm{div}\,}}(\sigma-\tilde{\sigma}),\phi)=(f-Pf,\phi-P\phi)$,
and by (2.5), $\|\operatorname{dev}(\sigma-\tilde{\sigma})\|\leq
C\|\sigma-\tilde{\sigma}\|_{\mathcal{A}}$. Hence we obtain

$C\|\sigma-\tilde{\sigma}\|_{L^{2}({\Omega})}\leq\|\sigma-\tilde{\sigma}\|_{\mathcal{A}}+h\|f-Pf\|_{L^{2}({\Omega})}.$

The proof is finished by combining this with (3.19). ∎

To summarize the results in this section, we have shown (in Theorem 3.5)
optimal error estimates for both $\hat{\sigma}_{h}$ and $\hat{u}_{h}$ when
$\hat{u}_{h}$ is piecewise constant, superconvergence of order $O(h^{2})$ for
$P(u-\hat{u}_{h})$ for general data $f$, and superconvergence of order
$O(h^{3})$ for $P(u-\hat{u}_{h})$ when $f\in W_{h}$ (in Corollary 3.8). We
have also shown (in Theorem 3.9) the robustness of the method in the
incompressible limit.

##  4\. A reduced space pair

In this section we present a finite element method obtained by reducing
certain judiciously chosen dimensions of the 3D Johnson–Mercier pair. We
replace the space $V_{h}$ by the space of element-wise rigid displacements and
reduce $\Sigma_{h}$ accordingly. The resulting method has optimal error
estimates. It is also robust for nearly incompressible materials (a property
that we cannot prove for an even further reduced pair of spaces given in
Appendix A).

On a facet with normal $n$ we let $v_{t}=n\times(v\times n)$ for vector fields
$v$, while for matrix fields $\omega$, we let $\omega_{nn}=\omega n\cdot n$
and $\omega_{nt}=n\times(\omega n\times n)\equiv(\omega n)_{t}$. Define

$\displaystyle\Sigma_{h}^{R}(T)$ $\displaystyle:=\\{\omega\in\Sigma_{h}(T)$
$\displaystyle:\;$ $\displaystyle{\mathop{\mathrm{div}\,}}\omega\in
P_{T}\mathcal{R}(T),\;\omega_{nt}\in\mathcal{R}(F),{\text{ for all
}}F\in\triangle_{2}(T)\\},$ $\displaystyle V_{h}^{R}(T)$
$\displaystyle:=\mathcal{R}(T).$

Here $\mathcal{R}(F)$ is the space of rigid body motions defined on $F$. We
see that there are $18=6+4\times 3$ constraints imposed on $\Sigma_{h}(T)$ to
obtain $\Sigma_{h}^{R}(T)$. The corresponding global spaces are denoted by
$\Sigma_{h}^{R}$ and $V_{h}^{R}$.

######  Theorem 4.1.

An element $\omega\in\Sigma_{h}^{R}(T)$ is uniquely determined by the
following dofs:

(4.1a) $\displaystyle\int_{F}\omega_{nn}\,\kappa\,{},\qquad$
$\displaystyle\kappa\in\EuScript{P}_{1}(F),F\in\triangle_{2}(T),\qquad$ (12
dofs) (4.1b) $\displaystyle\int_{F}\omega_{nt}\cdot\kappa\,{},\qquad$
$\displaystyle\kappa\in\mathcal{R}(F),F\in\triangle_{2}(T),\qquad$ (12 dofs).

###### Proof.

We see that $\dim\Sigma_{h}^{R}(T)\geq\dim\Sigma_{h}(T)-18=24$ which is
exactly the number of dofs in (4.1). Now suppose that
$\omega\in\Sigma_{h}^{R}(T)$ and the dofs (4.1a), (4.1b) vanish. Then,
obviously, $\omega n=0$ on $\partial T$.

By the definition of the reduced stress space, there is a $v\in\mathcal{R}(T)$
such that $P_{T}v={\mathop{\mathrm{div}\,}}\omega$. Combined with the
vanishing trace $\omega n=0$ on $\partial T$,

$\displaystyle\|{\mathop{\mathrm{div}\,}}\omega\|_{L^{2}(T)}^{2}$
$\displaystyle=\int_{T}{\mathop{\mathrm{div}\,}}\omega\cdot
P_{T}v=\int_{T}{\mathop{\mathrm{div}\,}}\omega\cdot v$
$\displaystyle=\int_{\partial T}\omega n\cdot v=0.$

This shows that ${\mathop{\mathrm{div}\,}}\omega=0$. Then by Lemma 2.1,
$\omega$ must vanish. ∎

The canonical interpolant of the dofs in (4.1), denoted by
$\Pi^{R}_{T}\omega\in\Sigma_{h}^{R}(T)$, satisfies

(4.2a)
$\displaystyle\int_{F}(\Pi^{R}_{T}\omega)_{nn}\,\kappa\,{}=\int_{F}\omega_{nn}\,\kappa\,{},\qquad$
$\displaystyle\kappa\in\EuScript{P}_{1}(F),F\in\triangle_{2}(T),$ (4.2b)
$\displaystyle\int_{F}(\Pi^{R}_{T}\omega)_{nt}\cdot\kappa\,{}=\int_{F}\omega_{nt}\cdot\kappa\,{}\qquad$
$\displaystyle\kappa\in\mathcal{R}(F),F\in\triangle_{2}(T),$

and we denote the corresponding global projection by $\Pi^{R}$.

######  Lemma 4.2.

For any $\omega\in D_{\Pi}$,

(4.3) $({\mathop{\mathrm{div}\,}}(\Pi^{R}\omega-\omega),v)=0\qquad{\text{ for
all }}v\in V_{h}^{R}.$

###### Proof.

Let $v\in\mathcal{R}(T)$. Then,

$\displaystyle({\mathop{\mathrm{div}\,}}(\Pi^{R}\omega-\omega),v)_{T}$
$\displaystyle=\int_{\partial T}(\Pi^{R}\omega-\omega)n\cdot v$
$\displaystyle=\int_{\partial T}(\Pi^{R}\omega-\omega)_{nn}\,v\cdot
n+\int_{\partial T}(\Pi^{R}\omega-\omega)_{nt}\cdot v_{t}=0,$

since $v_{t}\in\mathcal{R}(F)$. ∎

The following two results are similar to their previous counterparts and their
proofs are omitted.

######  Proposition 4.3.

For all $\omega\in H^{1}({\Omega},\mathbb{S}),$

(4.4)
$\displaystyle\|\omega-\Pi^{R}\omega\|_{L^{2}(\Omega)}+h\|{\mathop{\mathrm{div}\,}}(\omega-\Pi^{R}\omega)\|_{L^{2}(\Omega)}$
$\displaystyle\leq Ch\|\omega\|_{H^{1}(\Omega)}.$

######  Theorem 4.4.

There exists a constant $\beta>0$ such that

$\beta\leq\inf_{0\not=v\in
V^{R}_{h}}\sup_{0\neq\omega\in\Sigma_{h}^{R}}\frac{({\mathop{\mathrm{div}\,}}\omega,v)}{\|\omega\|_{H({\mathop{\mathrm{div}\,}},\Omega)}\|v\|_{L^{2}(\Omega)}}.$

The mixed method with the first reduced space pair finds
$\sigma^{R}_{h}\in\Sigma^{R}_{h}$ and $u^{R}_{h}\in V^{R}_{h}$ satisfying

(4.5a)
$\displaystyle(\mathcal{A}\sigma^{R}_{h},\tau)+(u^{R}_{h},{\mathop{\mathrm{div}\,}}\tau)$
$\displaystyle=0\qquad$ $\displaystyle{\text{ for all
}}\tau\in\Sigma^{R}_{h},$ (4.5b)
$\displaystyle({\mathop{\mathrm{div}\,}}\sigma^{R}_{h},v)$
$\displaystyle=(f,v)\qquad$ $\displaystyle{\text{ for all }}v\in V^{R}_{h}.$

We present an a priori error estimate for this method.

######  Theorem 4.5.

Let $\sigma,u$ solve (2.6) and $\sigma^{R}_{h},u^{R}_{h}$ solve (4.5) then the
following holds

(4.6)
$\|\Pi^{R}\sigma-\sigma^{R}_{h}\|_{\mathcal{A}}\leq\|\Pi^{R}\sigma-\sigma\|_{\mathcal{A}},$

and

(4.7) $\|u-u^{R}_{h}\|_{L^{2}(\Omega)}\leq
C\big{(}\|\sigma-\sigma^{R}_{h}\|_{\mathcal{A}}+\|Q^{R}u-u\|_{L^{2}(\Omega)}\big{)}$

where $Q^{R}:[L^{2}(\Omega)]^{3}\mapsto V^{R}_{h}$ is the $L^{2}$-orthogonal
projection onto $V^{R}_{h}$.

###### Proof.

Using (4.3), (4.5b) and (2.6b) we get that

$({\mathop{\mathrm{div}\,}}(\Pi^{R}\sigma-\sigma^{R}_{h}),v)=0\quad{\text{ for
all }}v\in V^{R}_{h}.$

However,
$({\mathop{\mathrm{div}\,}}(\Pi^{R}\sigma-\sigma^{R}_{h}),Pv)=({\mathop{\mathrm{div}\,}}(\Pi^{R}\sigma-\sigma^{R}_{h}),v)$.
Recalling (2.12), since $P_{T}:\mathcal{R}(T)\rightarrow P_{T}\mathcal{R}(T)$
is bijection, we conclude that

(4.8) ${\mathop{\mathrm{div}\,}}(\Pi^{R}\sigma-\sigma^{R}_{h})=0.$

Then, (4.6) follows from (4.5a) and (2.6a). To prove (4.7) one follows the
same lines as the proof of (2.37). We leave the details to the reader. ∎

######  Corollary 4.6.

Let $\sigma,u$ solve (2.6) and $\sigma^{R}_{h},u^{R}_{h}$ solve (4.5). If
$\sigma\in H^{1}(\Omega,\mathbb{S})$ and $u\in[H^{1}(\Omega)]^{3}$, then

(4.9) $\displaystyle\|\sigma-\sigma^{R}_{h}\|_{\mathcal{A}}$
$\displaystyle\leq Ch\|\sigma\|_{H^{1}(\Omega)},$ (4.10)
$\displaystyle\|u-u^{R}_{h}\|_{L^{2}(\Omega)}$ $\displaystyle\leq
Ch(\|\sigma\|_{H^{1}(\Omega)}+\|u\|_{H^{1}(\Omega)}).$

######  Theorem 4.7.

Assume that only (2.5) holds. Let $\sigma,u$ solve (2.6) and
$\sigma_{h}^{R},u_{h}^{R}$ solve (2.35). Then,

(4.11) $\displaystyle\|\sigma-\sigma_{h}^{R}\|_{L^{2}(\Omega)}\leq
C(\|\sigma-\sigma_{h}^{R}\|_{\mathcal{A}}+h\|\operatorname{div}(\sigma-\Pi^{R}\sigma)\|_{L^{2}(\Omega)}).$

###### Proof.

The proof proceeds along the same lines as that of Theorem 2.11 since the
analogue of the critical ingredient mentioned in Remark 2.12 is supplied by
(4.8). ∎

##  5\. Connection to stress elements with weakly imposed symmetry

In this section we discuss connection of the methods proposed in Sections 2
and 3 to mixed finite element methods for elasticity using stresses with
weakly imposed stress symmetry, i.e., the symmetry of $\sigma$ is imposed
weakly by the additional variational equation

$\displaystyle\int_{\Omega}\sigma:\eta=0,\quad\eta\in\Gamma=L^{2}(\Omega,\mathbb{K})$

where $\mathbb{K}$ is the space of $N\times N$ skew-symmetric matrices, with
$N=2$ or $3$. Introducing $\rho={\mathop{\mathrm{grad}\,}}u-\epsilon(u)$, the
skew-symmetric part of ${\mathop{\mathrm{grad}\,}}u$, the alternative mixed
formulation finds $\sigma\in\Sigma^{w}$, $u\in V$, $\rho\in\Gamma$ such that

(5.1a)
$\displaystyle(\bar{\mathcal{A}}\sigma,\tau)+(u,{\mathop{\mathrm{div}\,}}\tau)+(\rho,\tau)$
$\displaystyle=0\qquad$ $\displaystyle{\text{ for all }}\tau\in\Sigma^{w},$
(5.1b) $\displaystyle({\mathop{\mathrm{div}\,}}\sigma,v)$
$\displaystyle=(f,v)\qquad$ $\displaystyle{\text{ for all }}v\in V,$ (5.1c)
$\displaystyle(\sigma,\eta)$ $\displaystyle=0\qquad$ $\displaystyle{\text{ for
all }}\eta\in\Gamma,$

where $\Sigma^{w}:=H({\mathop{\mathrm{div}\,}},\Omega,\mathbb{M})$ and
$\bar{\mathcal{A}}$ is an operator which identical to $\mathcal{A}$ for
$L^{2}(\Omega,\mathbb{S})$ and is bounded and coercive on
$L^{2}(\Omega,\mathbb{K})$. Since $\Sigma^{w}_{h}\subset\Sigma^{w}$ can be a
matrix-valued finite element space such that each row is a standard
$H(\text{div})$ finite element, it is easier to construct stable mixed methods
for (5.1).

Finite element subspaces $\Sigma^{w}_{h}\times
V^{w}_{h}\times\Gamma_{h}\subset\Sigma^{w}\times V\times\Gamma$ giving stable
mixed methods for (5.1) of the form

(5.2a)
$\displaystyle(\bar{\mathcal{A}}\hat{\sigma}_{h}^{w},\tau)+(\hat{u}_{h}^{w},{\mathop{\mathrm{div}\,}}\tau)+(\rho_{h},\tau)$
$\displaystyle=0\qquad$ $\displaystyle{\text{ for all
}}\tau\in\Sigma^{w}_{h},$ (5.2b)
$\displaystyle({\mathop{\mathrm{div}\,}}\hat{\sigma}_{h}^{w},v)$
$\displaystyle=(f,v)\qquad$ $\displaystyle{\text{ for all }}v\in V^{w}_{h},$
(5.2c) $\displaystyle(\hat{\sigma}_{h}^{w},\eta)$ $\displaystyle=0\qquad$
$\displaystyle{\text{ for all }}\eta\in\Gamma_{h},$

has long been studied [5, 35, 34, 21, 3, 11, 24]. In [24] it was proved that
on meshes where each element $T$ has a Clough–Tocher or Alfeld split, the
spaces

(5.3a) $\displaystyle\Sigma^{w}_{h}$ $\displaystyle=\\{\omega\in
H({\mathop{\mathrm{div}\,}},{\Omega};\mathbb{M}):\omega|_{K}\in\EuScript{P}_{k}(K,\mathbb{M}),{\text{
for all }}K\in T^{\text{A}},T\in\mathcal{T}_{h}\\},$ (5.3b) $\displaystyle
V^{w}_{h}$
$\displaystyle=\\{v\in[L^{2}(T)]^{3}:v|_{K}\in[\EuScript{P}_{k-1}(K)]^{d},{\text{
for all }}K\in T^{\text{A}},T\in\mathcal{T}_{h}\\},$ (5.3c)
$\displaystyle\Gamma_{h}$ $\displaystyle=\\{\eta\in
L^{2}(T,\mathbb{K}):\eta|_{K}\in\EuScript{P}_{k}(K,\mathbb{K}),{\text{ for all
}}K\in T^{\text{A}},T\in\mathcal{T}_{h}\\}$

give stable methods for $k\geq N-1$ in the $N$-dimensional elasticity
formulation for $N=2,3$. Note that the polynomial degrees of $\Sigma^{w}_{h}$
and $\Gamma_{h}$ are the same in (5.3). Hence, the numerical solution
$\hat{\sigma}_{h}^{w}$ of (5.3) is exactly symmetric pointwise. In [32,
Subsection 4.2.3], this result was extended to $k\geq 1$ for $N=3$. Note that
when $k=1$, the space $V^{w}_{h}$ in (5.3b) equals $W_{h}$. In particular, it
was proved in [32] that

(5.4) given any $w\in W_{h}$ and $\eta\in\Gamma_{h},$ there exists
$\tau\in\Sigma^{w}_{h}$ satisfying ${\mathop{\mathrm{div}\,}}\tau=w$,
$(\tau,\eta)=\|\eta\|_{L^{2}(\Omega)}^{2}$, and
$C\|\tau\|_{H(\operatorname{div},\Omega)}\leq\|w\|_{L^{2}(\Omega)}+\|\eta\|_{L^{2}(\Omega)}$.

The inf-sup stability of (5.3) follows from (5.4).

######  Proposition 5.1.

The first two components of the unique solution
$(\hat{\sigma}_{h}^{w},\hat{u}_{h}^{w},\rho_{h})$ of (5.3), when $k=1$,
coincide with the solution $(\hat{\sigma}_{h},\hat{u}_{h})$ of (3.2).

###### Proof.

By (5.2c), $\hat{\sigma}_{h}^{w}$ is symmetric, so it is in $\Sigma_{h}$, the
three-dimensional Johnson–Mercier element. By restricting $\tau_{h}$ to
$\Sigma_{h}$, (5.2a) and (5.2b) are exactly same as (3.2) because
$\bar{\mathcal{A}}=\mathcal{A}$ on
$H({\mathop{\mathrm{div}\,}},\Omega,\mathbb{S})$. Therefore,
$\hat{\sigma}_{h}^{w}=\hat{\sigma}_{h}$ and $\hat{u}_{h}^{w}=\hat{u}_{h}$ by
uniqueness of solutions. ∎

Next, consider the case where the displacement is discretized in the space of
piecewise linear functions on the unsplit mesh $\mathcal{T}_{h}$, namely the
space $V_{h}$ from Section 2, i.e.,

(5.5a)
$\displaystyle(\bar{\mathcal{A}}\sigma_{h}^{w},\tau)+(u_{h}^{w},{\mathop{\mathrm{div}\,}}\tau)+(\rho_{h},\tau)$
$\displaystyle=0\qquad$ $\displaystyle{\text{ for all
}}\tau\in\Sigma^{w}_{h},$ (5.5b)
$\displaystyle({\mathop{\mathrm{div}\,}}\sigma_{h}^{w},v)$
$\displaystyle=(f,v)\qquad$ $\displaystyle{\text{ for all }}v\in V_{h},$
(5.5c) $\displaystyle(\sigma_{h}^{w},\eta)$ $\displaystyle=0\qquad$
$\displaystyle{\text{ for all }}\eta\in\Gamma_{h}.$

######  Proposition 5.2.

The method (5.5) is inf-sup stable and the first two components of its unique
solution $(\sigma_{h}^{w},u_{h}^{w},\rho_{h})\in\Sigma^{w}_{h}\times
V_{h}\times\Gamma_{h}$ coincide with $(\sigma_{h}^{w},u_{h}^{w})$ solving
(2.35).

###### Proof.

Given any $(v,\eta)\in V_{h}\times\Gamma_{h}$, applying (5.4) to $(Pv,\eta)\in
W_{h}\times\Gamma_{h}$, we find that there exists $\tau\in\Sigma^{w}_{h}$
satisfying

(5.6)
${\mathop{\mathrm{div}\,}}\tau=Pv,\quad(\tau,\eta)=\|\eta\|_{L^{2}(\Omega)}^{2},\quad
C\|\tau\|_{H(\operatorname{div},\Omega)}\leq\|Pv\|_{L^{2}(\Omega)}+\|\eta\|_{L^{2}(\Omega)}.$

Using this $\tau$,

$\displaystyle\sup_{\omega\in\Sigma^{w}_{h}}\frac{(v,{\mathop{\mathrm{div}\,}}\omega)+(\eta,\omega)}{\|\omega\|_{H(\operatorname{div},{\Omega})}}$
$\displaystyle\geq\frac{(v,{\mathop{\mathrm{div}\,}}\tau)+(\eta,\tau)}{\|\tau\|_{H(\operatorname{div},{\Omega})}}\geq
C(\|Pv\|^{2}_{L^{2}({\Omega})}+\|\eta\|_{L^{2}({\Omega})}^{2})^{1/2}$
$\displaystyle\geq
C(\|v\|^{2}_{L^{2}({\Omega})}+\|\eta\|^{2}_{L^{2}({\Omega})})^{1/2},$

where we have used (2.50a). Thus, (5.5) is inf-sup stable and uniquely
solvable. Now we can use the same uniqueness argument as in the proof of
Proposition 5.1 to show that $(\sigma_{h}^{w},u_{h}^{w})$ and the solution of
(2.35) are equal. ∎

To conclude this section, we have shown that it is possible to compute the
solutions of methods in Sections 2 and 3 using completely standard finite
element spaces via formulations (5.2) and (5.5). However, we do not know how
to arrive at the reduced space $\Sigma_{h}^{R}\times V_{h}^{R}$ of Section 4
through methods with weakly imposed stress symmetry. The space
$\Sigma_{h}^{R}(T)$ does not have full $\EuScript{P}_{1}(K,\mathbb{M})$
polynomials for $K\in T^{\text{A}}$, so finding an appropriate
$\Gamma_{h}^{R}\subset\Gamma$ which gives exact symmetry of stress tensors and
satisfies the inf-sup stability, appears to be nontrivial.

##  6\. The $N$-dimensional case

This section generalizes the Johnson–Mercier element to higher dimensions. To
generalize our arguments in Section 2, we need an $N$-dimensional version of
the identity (2.14). This can be found, described in the language of exterior
calculus, in [7, 8, 20], using the so-called “BGG resolution” [3, 13, 18].
Nonetheless, we start by giving an elementary and self-contained proof of such
an identity using standard vector notation in Subsection 6.1. We then use it
to develop the Johnson–Mercier element in $\mathbb{R}^{N}$ for any $N\geq 2$.

###  6.1. A commuting diagram of tensor fields for $N$-dimensional elasticity

Let $\mathbb{V}=\mathbb{R}^{N}$, let $\mathbb{M}$ denote the space of $N\times
N$ matrices, and let $\mathbb{K}=\operatorname{skw}(\mathbb{M})$ denote the
space of skew symmetric matrices. We use $\otimes$ to denote the tensor
product of vector spaces. Then for (column) vectors $w,v\in\mathbb{V}$,
$w\otimes v=wv^{t}$ equals the outer product. Let $e_{i}$ denote the standard
unit basis of $\mathbb{V}$, $\mathsf{E}_{ij}=2\operatorname{skw}(e_{i}\otimes
e_{j})=e_{i}e_{j}^{t}-e_{j}e_{i}^{t}\in\mathbb{K}$. Let $\Lambda(\mathbb{X})$
be $C^{\infty}$-smooth $\mathbb{X}$-valued fields where $\mathbb{X}$ is one of
$\mathbb{V},\mathbb{M},$ $\mathbb{K}$ or tensor products of these. Henceforth,
we use Einstein’s summation convention.

We continue to use the standard divergence operator on vector fields in
$\Lambda(\mathbb{V})$ and the (previously used) row-wise divergence operator
on matrix fields in $\Lambda(\mathbb{M})$. This is a natural to extend
$\operatorname{div}$ to $\mathbb{K}\otimes\Lambda(\mathbb{V})$ by

${\mathop{\mathrm{div}\,}}(\eta\otimes
w)=\eta\otimes({\mathop{\mathrm{div}\,}}w),\,\qquad\eta\in\mathbb{K},\;w\in\Lambda(\mathbb{V}).$

This motivates the following extension of divergence to the isomorphic space
$\Lambda(\mathbb{K}\otimes\mathbb{V})$: writing an element of
$\Lambda(\mathbb{K}\otimes\mathbb{V})$ as $b=b_{ijk}\mathsf{E}_{ij}\otimes
e_{k}$ for some scalar component functions $b_{ijk}(x)$, we define
$\operatorname{div}:\Lambda(\mathbb{K}\otimes\mathbb{V})\to\Lambda(\mathbb{K})$
by

(6.1) $\operatorname{div}(b_{ijk}\mathsf{E}_{ij}\otimes
e_{k})=(\partial_{k}b_{ijk})\mathsf{E}_{ij}.$

Next, define a differential operator
$\mathsf{d}:\Lambda(\mathbb{K})\rightarrow\Lambda(\mathbb{V})$ by

(6.2)
$\mathsf{d}(\eta_{ij}\mathsf{E}_{ij}):=\partial_{j}(\eta_{ij}-\eta_{ji})e_{i},\qquad\eta=\eta_{ij}\mathsf{E}_{ij}\in\Lambda(\mathbb{K}).$

It is natural to extend $\mathsf{d}$ to $\mathbb{V}\otimes\Lambda(\mathbb{K})$
by

$\mathsf{d}(v\otimes\omega)=v\otimes\mathsf{d}\omega,\qquad
v\in\mathbb{V},\omega\in\Lambda(\mathbb{K}).$

Since $\Lambda(\mathbb{V}\otimes\mathbb{K})$ is isomorphic to
$\mathbb{V}\otimes\Lambda(\mathbb{K})$, this motivates the definition of the
extension of $\mathsf{d}$ in (6.2) to
$\mathsf{d}:\Lambda(\mathbb{V}\otimes\mathbb{K})\to\Lambda(\mathbb{M})$ by

(6.3)
$\mathsf{d}(a_{ijk}e_{i}\otimes\mathsf{E}_{jk})=\partial_{k}(a_{ijk}-a_{ikj})\,e_{i}\otimes
e_{j},\qquad
a_{ijk}e_{i}\otimes\mathsf{E}_{jk}\in\Lambda(\mathbb{V}\otimes\mathbb{K}).$

######  Remark 6.1.

Alternately, we can arrive at these definitions through a generalized
divergence operator [33], $\operatorname{Div}$, that acts on
$\Lambda(\mathbb{V}\otimes\mathbb{V})$ and
$\Lambda(\mathbb{V}\otimes\mathbb{V}\otimes\mathbb{V})$ by

(6.4) $\displaystyle\operatorname{Div}(a_{ij}e_{i}\otimes
e_{j})=(\partial_{j}a_{ij})e_{i},\qquad\operatorname{Div}(a_{ijk}e_{i}\otimes
e_{j}\otimes e_{k})=(\partial_{k}a_{ijk})e_{i}\otimes e_{j}.$

Since $\mathbb{K}\otimes\mathbb{V}$ and $\mathbb{V}\otimes\mathbb{K}$ are both
subspaces of $\mathbb{V}\otimes\mathbb{V}\otimes\mathbb{V}$, the latter
defines $\operatorname{Div}$ on $\Lambda(\mathbb{K}\otimes\mathbb{V})$ and
$\Lambda(\mathbb{V}\otimes\mathbb{K})$. Since

$\displaystyle\eta$
$\displaystyle=\eta_{ij}\mathsf{E}_{ij}=(\eta_{ij}-\eta_{ji})e_{i}\otimes
e_{j}$ $\displaystyle\in\Lambda(\mathbb{K}),$ $\displaystyle a$
$\displaystyle=a_{ijk}e_{i}\otimes\mathsf{E}_{jk}=(a_{ijk}-a_{ikj})e_{i}\otimes
e_{j}\otimes e_{k}$ $\displaystyle\in\Lambda(\mathbb{V}\otimes\mathbb{K}),$
$\displaystyle b$ $\displaystyle=b_{ijk}\mathsf{E}_{ij}\otimes
e_{k}=(b_{ijk}-b_{jik})e_{i}\otimes e_{j}\otimes e_{k}$
$\displaystyle\in\Lambda(\mathbb{K}\otimes\mathbb{V}),$

applying (6.4) we easily see that equations (6.1), (6.2), and (6.3) are the
same as

$\operatorname{div}b=\operatorname{Div}b,\quad\mathsf{d}\eta=\operatorname{Div}\eta,\quad\mathsf{d}a=\operatorname{Div}a,$

respectively.

With these definitions, we have

(6.5) $\displaystyle{\mathop{\mathrm{div}\,}}\mathsf{d}\,\eta$
$\displaystyle=0,$ $\displaystyle{\text{ for all
}}\eta\in\Lambda(\mathbb{K}),$ (6.6)
$\displaystyle{\mathop{\mathrm{div}\,}}\mathsf{d}\,\omega$ $\displaystyle=0,$
$\displaystyle{\text{ for all
}}\omega\in\Lambda(\mathbb{V}\otimes\mathbb{K}).$

Indeed, to see that (6.5) holds, it suffices to observe that
${\mathop{\mathrm{div}\,}}\mathsf{d}\eta={\mathop{\mathrm{div}\,}}(\partial_{j}(\eta_{ij}-\eta_{ji})e_{i})=\partial_{i}\partial_{j}(\eta_{ij}-\eta_{ji})=0.$
Essentially the same argument also shows (6.6). In the language of exterior
calculus, (6.5) and (6.6) imply that the sequences

(6.8)
${\Lambda(\mathbb{K})}$${\Lambda(\mathbb{V})}$${\Lambda(\mathbb{R}),}$$\scriptstyle{\mathsf{d}}$$\scriptstyle{{\mathop{\mathrm{div}\,}}}$
(6.10)
${\Lambda(\mathbb{V}\otimes\mathbb{K})}$${\Lambda(\mathbb{M})}$${\Lambda(\mathbb{V}),}$$\scriptstyle{\mathsf{d}}$$\scriptstyle{{\mathop{\mathrm{div}\,}}}$

form complexes, elucidating the connection to relevant $N$-dimensional de Rham
complexes. Note that the tensor products of $\mathbb{V}$ with function spaces
in (6.8) yield their analogue in (6.10) up to an isomorphism.

Next, we define a linear operator
$\Theta:\Lambda(\mathbb{V}\otimes\mathbb{K})\rightarrow\Lambda(\mathbb{K}\otimes\mathbb{V})$
by defining $\Theta a$ for $a=a_{ijk}e_{i}\otimes e_{j}\otimes e_{k}$ in the
larger set $\Lambda(\mathbb{V}\otimes\mathbb{V}\otimes\mathbb{V})$, namely

(6.11) $\Theta(a_{ijk}e_{i}\otimes e_{j}\otimes
e_{k}):=(a_{ijk}-a_{jik})\;e_{i}\otimes e_{j}\otimes e_{k}$

When $a$ is in $\Lambda(\mathbb{V}\otimes\mathbb{K})$, we have
$a=a_{ijk}e_{i}\otimes e_{j}\otimes
e_{k}=\frac{1}{2}a_{ijk}e_{i}\otimes\mathsf{E}_{jk}$. Hence, (6.11) can
equivalently be expressed as

(6.12)
$\Theta(a_{ijk}e_{i}\otimes\mathsf{E}_{jk})=2a_{ijk}\mathsf{E}_{ij}\otimes
e_{k}.$

######  Lemma 6.2.

The operator
$\Theta:\Lambda(\mathbb{V}\otimes\mathbb{K})\rightarrow\Lambda(\mathbb{K}\otimes\mathbb{V})$
is invertible and its inverse is given by

$\Theta^{-1}\,b=\frac{1}{2}(b_{ijk}-b_{ikj}-b_{jki})\,e_{i}\otimes
e_{j}\otimes e_{k}$

for any $b=b_{ijk}\,e_{i}\otimes e_{j}\otimes
e_{k}\in\Lambda(\mathbb{K}\otimes\mathbb{V})$. Moreover, for all
$\omega\in\Lambda(\mathbb{V}\otimes\mathbb{K})$,

(6.13) $\displaystyle{\mathop{\mathrm{div}\,}}\Theta\,\omega$
$\displaystyle=2\operatorname{skw}\mathsf{d}\,\omega.$

###### Proof.

By direct calculation, it is easy to verify that the given expression for
$\Theta^{-1}$ satisfies $\Theta(\Theta^{-1}b)=b$ for all
$b\in\Lambda(\mathbb{K}\otimes\mathbb{V})$ and $\Theta^{-1}(\Theta a)=a$ for
all $a\in\Lambda(\mathbb{V}\otimes\mathbb{K})$. To detail the latter, letting
$\theta=\Theta a=(a_{ijk}-a_{jik})e_{i}\otimes e_{j}\otimes e_{k},$ observe
that

$\displaystyle 2\Theta^{-1}\Theta a$
$\displaystyle=(\theta_{ijk}-\theta_{ikj}-\theta_{jki})e_{i}\otimes
e_{j}\otimes e_{k}$
$\displaystyle=((a_{ijk}-a_{jik})-(a_{ikj}-a_{kij})-(a_{jki}-a_{kji}))e_{i}\otimes
e_{j}\otimes e_{k}$

When $a\in\Lambda(\mathbb{V}\otimes\mathbb{K})$, its components have skew
symmetry in the last two indices (i.e., $a_{ijk}=-a_{ikj}$), so the last
expression above simplifies to $2a_{ijk}e_{i}\otimes e_{j}\otimes e_{k}=2a$.

Next, to prove (6.13), let
$\omega=\omega_{ijk}e_{i}\otimes\mathsf{E}_{jk}\in\Lambda(\mathbb{V}\otimes\mathbb{K})$.
Then, by (6.1) and (6.12),

$\displaystyle\operatorname{div}(\Theta\omega)$
$\displaystyle=\operatorname{div}(2\omega_{ijk}\mathsf{E}_{ij}\otimes
e_{k})=2(\partial_{k}\omega_{ijk})\,\mathsf{E}_{ij}.$

Also, by (6.3),

$\displaystyle 2\operatorname{skw}\mathsf{d}\omega$
$\displaystyle=\partial_{k}(\omega_{ijk}-\omega_{ikj})\,2\operatorname{skw}(e_{i}\otimes
e_{j})$
$\displaystyle=\partial_{k}(\omega_{ijk}-\omega_{jik})\;\mathsf{E}_{ij}=2(\partial_{k}\omega_{ijk})\,\mathsf{E}_{ij}.$

Thus the left and right hand sides of (6.13) are equal. ∎

In other words, we have shown that the following diagram commutes and that its
top and bottom rows are complexes.

${\Lambda(\mathbb{K}\otimes\mathbb{V})}$${\Lambda(\mathbb{K})}$${\Lambda(\mathbb{V}\otimes\mathbb{K})}$${\Lambda(\mathbb{M})}$${\Lambda(\mathbb{V}).}$$\scriptstyle{{\mathop{\mathrm{div}\,}}}$$\scriptstyle{\mathsf{d}}$$\scriptstyle{\Theta}$$\scriptstyle{{\mathop{\mathrm{div}\,}}}$$\scriptstyle{2\mathrm{skw}}$

For later use, we state the following result about the trace of $\mathsf{d}$
on a hyperplane.

######  Lemma 6.3.

Let $\eta=\eta_{ij}\mathsf{E}_{ij}\in\Lambda(\mathbb{K})$ and suppose that
$\Gamma$ is a hyperplane in which the components of $\eta$ vanish then
$\mathsf{d}\eta\cdot n$ vanishes on $\Gamma$ where $n$ is a unit normal to
$\Gamma$.

###### Proof.

After possibly rotating coordinates we can assume without loss of generality
that $n=e_{N}$. Then,

$\mathsf{d}\eta\cdot n=\partial_{j}(\eta_{ij}-\eta_{ji})e_{i}\cdot
e_{N}=\partial_{j}(\eta_{Nj}-\eta_{jN})=\sum_{j=1}^{N-1}\partial_{j}(\eta_{Nj}-\eta_{jN}).$

The latter vanish on $\Gamma$ since $\eta_{Nj}-\eta_{jN}$ vanish on $\Gamma$
and $\partial_{j}$ (for $1\leq j\leq N-1$) are tangential derivatives on
$\Gamma$. ∎

Of course, we can define the weak version of $\mathsf{d}$ and extend the
domain of $\mathsf{d}$. To this end, let $\eta=\eta_{ij}\mathsf{E}_{ij}$ we
say $\mathsf{d}\eta$ exists if there exists $\omega\in
L^{1}(\Omega,\mathbb{V})$ such that

(6.14)
$(\omega,\phi_{k}e_{k})=-(\eta_{ij}-\eta_{ji},\partial_{j}\phi_{k})\delta_{ik},\qquad{\text{
for all }}\phi=\phi_{k}e_{k}\in C_{0}^{\infty}(\Omega,\mathbb{V}),$

and we set $\mathsf{d}\eta=\omega$. Here $\delta_{ik}$ is the Kronecker delta.
We then define $H(\mathsf{d},\Omega):=\\{\eta\in
L^{2}(\Omega,\mathbb{K}):\mathsf{d}\eta\in L^{2}(\Omega,\mathbb{V})\\}$. The
identities (6.6) and (6.13) hold for $\omega\in\mathbb{V}\otimes
H(\mathsf{d},\Omega)$.

###  6.2. The stress element on the Alfeld split in higher dimensions

For an $N$-simplex $T=[x_{0},x_{1},x_{2},\cdots,x_{N}]$, let $T^{\text{A}}$ be
an Alfeld split of $T$, i.e., choosing an interior point $z$ of $T$, define
$T^{\text{A}}=\\{T_{0},T_{1},T_{2},\cdots,T_{N}\\}$ with $N$-simplices
$T_{i}=[z,x_{0},\cdots,\widehat{x_{i}},\cdots,x_{N}]$ where $\widehat{x_{i}}$
means that $x_{i}$ is not in the vertices of $T_{i}$. For a given
triangulation $\mathcal{T}_{h}$ of $\Omega$, we let
$\mathcal{T}_{h}^{\text{A}}$ be the resulting triangulation after performing
an Alfeld split to each $T\in\mathcal{T}_{h}$. On each macro element
$T\in\mathcal{T}_{h}$, define the local spaces

$\displaystyle\EuScript{P}_{k}(T^{\text{A}},\mathbb{X}):=$
$\displaystyle\\{\omega\in
L^{2}(T,\mathbb{X}):\omega|_{K}\in\EuScript{P}_{k}(K,\mathbb{X}),{\text{ for
all }}K\in T^{\text{A}}\\},$ $\displaystyle\Sigma_{h}(T):=$ $\displaystyle
H(\operatorname{div},T,\mathbb{S})\cap\EuScript{P}_{1}(T^{\text{A}},\mathbb{S}),$
$\displaystyle V_{h}(T):=$ $\displaystyle\EuScript{P}_{1}(T,\mathbb{V}).$

The $N$-dimensional matrix analogue of the well-known BDM space and its
subspace with vanishing normal components on $\partial T$ are namely

$\displaystyle{B\\!D\\!M}_{1}(T^{\text{A}})$
$\displaystyle:=H(\operatorname{div},T,\mathbb{M})\cap\EuScript{P}_{1}(T^{\text{A}},\mathbb{M}),$
$\displaystyle\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}})$
$\displaystyle:=\\{\tau\in{B\\!D\\!M}_{1}(T^{\text{A}}):\tau n|_{\partial
T}=0\\}.$

Notice that the Alfeld split of $T$ has $N+1$ distinct sub-simplices in
$\triangle_{N-1}(T)$ and $\frac{1}{2}N(N+1)$ distinct internal subsimplices in
$\triangle_{N-1}(T^{\text{A}})$, therefore
$\dim{B\\!D\\!M}_{1}(T^{\text{A}})=(N+1)\cdot N^{2}+\frac{1}{2}N(N+1)\cdot
N^{2}=(\frac{N}{2}+1)(N+1)N^{2}$ and
$\dim\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}})=\frac{1}{2}(N+1)N^{3}$. Notice
also that
$\dim\EuScript{P}_{1}(T^{\text{A}},\mathbb{K})=\frac{1}{2}N(N-1)(N+1)^{2}$.

We will also need the space
$\mathring{\mathcal{L}}_{2}(T^{\text{A}}):=\EuScript{P}_{2}(T^{\text{A}})\cap\mathring{H}^{1}(T)$
and $\mathbb{X}$-valued versions
$\mathring{\mathcal{L}}_{2}(T^{\text{A}},\mathbb{X})=\mathring{\mathcal{L}}_{2}(T^{\text{A}})\otimes\mathbb{X}$.
We begin by recalling the next result from [22, Theorem 3.1] (see also [25]).

######  Proposition 6.4.

Let $v\in\EuScript{P}_{1}(T^{A})$ with $\int_{T}v=0$ there exists
$\rho\in\mathring{\mathcal{L}}_{2}(T^{\text{A}},\mathbb{V})$ such that

${\mathop{\mathrm{div}\,}}\rho=v.$

######  Proposition 6.5.

If $\rho\in\mathcal{L}_{2}(T^{\text{A}},\mathbb{V}\otimes\mathbb{K})$ then
$\mathsf{d}\rho\in{B\\!D\\!M}_{1}(T^{\text{A}})$. Moreover, if
$\rho\in\mathring{\mathcal{L}}_{2}(T^{\text{A}},\mathbb{V}\otimes\mathbb{K})$
then $\mathsf{d}\rho\in\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}})$.

###### Proof.

Let $\rho\in\mathcal{L}_{2}(T^{\text{A}},\mathbb{V}\otimes\mathbb{K})$ then
for each $K\in T^{\text{A}}$ $\mathsf{d}\rho|_{K}$ has linear components.
Hence, it is enough to show that $\rho\in H(\mathsf{d},T)\otimes\mathbb{V}$,
but this follows immediately from (6.14). If
$\rho\in\mathring{\mathcal{L}}_{2}(T^{\text{A}},\mathbb{V}\otimes\mathbb{K})$
then we must also show that $(\mathsf{d}\rho)\,n$ vanishes on $\partial T$,
but this follows from Lemma 6.3. ∎

######  Lemma 6.6.

The equality
$\\{\omega\in\Sigma_{h}(T):{\mathop{\mathrm{div}\,}}\omega=0,\omega
n|_{\partial T}=0\\}=\\{0\\}$ continues to hold in the $N$-dimensional case.

###### Proof.

Let
$\tilde{\EuScript}{P}_{0}(T^{\text{A}})=\\{u\in\EuScript{P}_{0}(T^{\text{A}}):\int_{T}u=0\\}$
and set

$X=\left\\{(\eta,v):\eta\in\EuScript{P}_{1}(T^{\text{A}},\mathbb{K}),\;v\in\tilde{\EuScript}{P}_{0}(T^{\text{A}},\mathbb{V}),\;\int_{T}(\eta:\kappa+v\cdot\kappa
x)=0\text{ for all }\kappa\in\mathbb{K}\right\\}$

where $\kappa x$ is the matrix-vector product with the coordinate vector
polynomial $x$. The proof is based on the operator
$A:\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}})\to X$ given by

(6.15)
$A\tau=(-\operatorname{skw}\tau,{\mathop{\mathrm{div}\,}}\tau),\qquad\tau\in\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}}).$

Observe that for any $\kappa\in\mathbb{K}$ and
$\tau\in\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}})$,

(6.16)
$\int_{T}\kappa:\operatorname{skw}\tau=\int_{T}{\mathop{\mathrm{grad}\,}}(\kappa
x):\tau=-\int_{T}(\kappa x)\cdot\operatorname{div}\tau.$

Thus $A$ indeed maps into $X$. We proceed to show that $A$ is surjective.

Let $(\eta,v)\in X$. Since $v$ has components of zero mean, by a standard
exact sequence property [4], there exists a
$\sigma\in\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}})$ such that

(6.17) ${\mathop{\mathrm{div}\,}}\sigma=v.$

Since $(\eta,v)\in X$, this implies

$\int_{T}(\eta-\operatorname{skw}\sigma):\kappa=\int_{T}(\eta:\kappa+(\operatorname{div}\sigma)\cdot\kappa
x)=\int_{T}(\eta:\kappa+v\cdot\kappa x)=0,$

where we have again used (6.16). Letting
$\zeta=\eta-\operatorname{skw}\sigma\in\EuScript{P}_{1}(T^{\text{A}},\mathbb{K}),$
the above shows all the components have zero mean. Hence, by Proposition 6.4
there exists
$\omega\in\mathring{\mathcal{L}}_{2}(T^{\text{A}},\mathbb{V}\otimes\mathbb{K})$
such that

(6.18) ${\mathop{\mathrm{div}\,}}\omega=\eta-\operatorname{skw}\sigma.$

We set $a=\Theta^{-1}\omega$ and note that
$a\in\mathring{\mathcal{L}}_{2}(T^{\text{A}},\mathbb{V}\otimes\mathbb{K})$.
Setting $\tau=\sigma+\mathsf{d}\,a$, then by Proposition 6.5 we have
$\tau\in\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}})$.

By the commuting diagram property (6.13),
$\operatorname{skw}\mathsf{d}\,a=\operatorname{div}\Theta
a=\operatorname{div}\omega$, so we have

$\displaystyle\operatorname{skw}\tau$
$\displaystyle=\operatorname{skw}\sigma+\operatorname{div}\omega=\eta,$ (by
(6.18)) $\displaystyle{\mathop{\mathrm{div}\,}}\tau$ $\displaystyle=v.$ (by
(6.17) and (6.6))

Thus we have proved that $A$ is surjective.

To conclude, the surjectivity of $A$ implies that

$\displaystyle\operatorname{rank}(A)$
$\displaystyle\geq\dim\EuScript{P}_{1}(T^{\text{A}},\mathbb{K})+\dim\tilde{\EuScript}{P}_{0}(T^{\text{A}},\mathbb{V})-\dim\mathbb{K}$
$\displaystyle=\frac{1}{2}N(N-1)(N+1)^{2}+N^{2}-\frac{1}{2}N(N-1)$
$\displaystyle=\frac{1}{2}N^{3}(N+1).$

Since $\dim\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}})=\frac{1}{2}N^{3}(N+1)$,
the rank-nullity theorem implies that the null space of $A$ is trivial, i.e.,
$\\{\tau\in\mathring{{B\\!D\\!M}}_{1}(T^{\text{A}}):\;\operatorname{skw}\tau=0,\;{\mathop{\mathrm{div}\,}}\tau=0\\}=\\{0\\}$.
∎

We now prove the main result of this section. The proof is the same as the
proof of Theorem 2.2.

######  Theorem 6.7.

The dimension of $\Sigma_{h}(T)$ is $(N+\frac{1}{2})N(N+1)$. Moreover, an
element $\omega\in\Sigma_{h}(T)$ is uniquely determined by the following dofs:

(6.19a) $\displaystyle\int_{F}\omega n_{T}\cdot\kappa\,{},\qquad$
$\displaystyle\kappa\in\EuScript{P}_{1}(F,\mathbb{V}),\quad
F\in\triangle_{N-1}(T),\qquad$ ($(N+1)N^{2}$ dofs) (6.19b)
$\displaystyle\int_{T}\omega.$ ($\frac{1}{2}N(N+1)$ dofs)

###### Proof.

Since
$\Sigma_{h}(T)=\\{\omega\in{B\\!D\\!M}_{1}(T^{\text{A}}):\int_{T}\omega:\eta=0\;\text{
for all }\eta\in\EuScript{P}_{1}(T^{\text{A}},\mathbb{K})\\}$,

$\displaystyle\dim\Sigma_{h}(T)$
$\displaystyle\geq\dim{B\\!D\\!M}_{1}(T^{\text{A}})-\dim\EuScript{P}_{1}(T^{\text{A}},\mathbb{K})$
$\displaystyle=(\frac{N}{2}+1)(N+1)N^{2}-\frac{N}{2}(N-1)(N+1)^{2}$
$\displaystyle=\frac{N}{2}(N+1)(2N+1).$

To show $\dim\Sigma_{h}(T)=\frac{N}{2}(N+1)(2N+1)$, it is sufficient to prove
that $\omega=0$ if the dofs (6.19) of $\omega\in\Sigma_{h}(T)$ vanish because
the number of dofs in (6.19) is $\frac{N}{2}(N+1)(2N+1)$. This also proves
that $\Sigma_{h}(T)$ is unisolvent by (6.19).

To this end, let $\omega\in\Sigma_{h}(T)$ have vanishing dofs (6.19). Then,
$\omega n=0$ on $\partial T$ by (6.19a). Moreover, using (6.19b) and
integration by parts gives

(6.20) $\displaystyle\int_{T}\operatorname{div}\omega\cdot
v=-\int_{T}\omega:\epsilon(v)=0\qquad{\text{ for all
}}v\in\EuScript{P}_{1}(T,\mathbb{V}).$

By a similar simple argument using quadrature rules as in the 3D case, we
conclude that $\operatorname{div}\omega=0$. Hence, by Lemma 6.6 we conclude
that $\omega=0$. ∎

##  Appendix A Another reduced space pair

In this section we reduce the spaces even further from the reduced pair in
Section 4. In this further reduced finite element method, $V_{h}$ is replaced
by only piecewise constant elements on $T\in\mathcal{T}_{h}$ and $\Sigma_{h}$
is also further reduced correspondingly. We prove optimal
$\mathcal{A}$-weighted $L^{2}$ error estimate for $\sigma$ and $L^{2}$ error
estimate for $u$. However, for this further reduced method, we are not able to
prove (an analogue of Theorem 2.11 giving) a robust error estimate for the
nearly incompressible case. The reason for this is that we now only have very
weak control of the divergence of the stresses. Hence we are unable to
recommend this method. Nonetheless, we present it as a curiosity, since it has
only three dofs per facet and yet remains stable.

The discrete spaces for this method are

$\displaystyle\tilde{\Sigma}_{h}^{R}(T):$
$\displaystyle=\\{\omega\in\Sigma_{h}(T):{\mathop{\mathrm{div}\,}}\omega\in
P_{T}\mathcal{R}(T),\omega n\in[\EuScript{P}_{0}(F)]^{3},{\text{ for all
}}F\in\triangle_{2}(T)\\},$ $\displaystyle\tilde{V}_{h}^{R}(T):$
$\displaystyle=[\EuScript{P}_{0}(T)]^{3}.$

The corresponding global spaces are denoted by
$\tilde{\Sigma}_{h}^{R},\tilde{V}_{h}^{R}$.

######  Theorem A.1.

An element $\omega\in\tilde{\Sigma}_{h}^{R}(T)$ is uniquely determined by the
following dofs:

(A.1) $\displaystyle\int_{F}\omega n\cdot\kappa\,{},\qquad$
$\displaystyle\kappa\in[\EuScript{P}_{0}(F)]^{3},F\in\triangle_{2}(T),\qquad$
(12 dofs).

###### Proof.

We see that $\dim\tilde{\Sigma}_{h}^{R}(T)\geq\dim\Sigma_{h}(T)-(6+6\cdot
4)=12$ which is exactly the number of dofs in (A.1). Now suppose that
$\omega\in\tilde{\Sigma}_{h}^{R}(T)$ and the dofs (A.1) vanish. Then, we have
$\omega n=0$ on $\partial T$.

Let $v\in\mathcal{R}(T)$ be such that $P_{T}v={\mathop{\mathrm{div}\,}}\omega$
then

$\displaystyle\|{\mathop{\mathrm{div}\,}}\omega\|_{L^{2}(T)}^{2}=\int_{T}{\mathop{\mathrm{div}\,}}\omega\cdot
P_{T}v=\int_{T}{\mathop{\mathrm{div}\,}}\omega\cdot v=\int_{\partial T}\omega
n\cdot v=0.$

Hence, this shows that ${\mathop{\mathrm{div}\,}}\omega=0$. By Lemma 2.1 we
have that $\omega$ vanishes. ∎

The corresponding projection with the degrees of freedom is defined as
follows.

(A.2a) $\displaystyle\int_{F}(\tilde{\Pi}^{R}_{T}\omega
n)\cdot\kappa\,{}=\int_{F}\omega n\cdot\kappa\,{},\qquad$
$\displaystyle\kappa\in[\EuScript{P}_{0}(F)]^{3},F\in\triangle_{2}(T).$

From this we define the global projection $\tilde{\Pi}^{R}$.

######  Lemma A.2.

It holds,

(A.3)
$({\mathop{\mathrm{div}\,}}(\tilde{\Pi}^{R}\omega-\omega),v)=0\qquad{\text{
for all }}v\in\tilde{V}_{h}^{R},\omega\in D_{\Pi}.$

###### Proof.

Let $v\in[\EuScript{P}_{0}(T)]^{3}$. Then,

$\displaystyle({\mathop{\mathrm{div}\,}}(\tilde{\Pi}^{R}\omega-\omega),v)_{T}$
$\displaystyle=\int_{\partial T}(\tilde{\Pi}^{R}\omega-\omega)n\cdot v=0.$

∎

We omit the proofs of the following results analogous to previous sections.

######  Proposition A.3.

For all $\omega\in H^{1}(\Omega,\mathbb{S})$,

(A.4)
$\displaystyle\|\omega-\tilde{\Pi}^{R}\omega\|_{L^{2}(\Omega)}+h\|{\mathop{\mathrm{div}\,}}(\omega-\tilde{\Pi}^{R}\omega)\|_{L^{2}(\Omega)}\leq$
$\displaystyle Ch\|\omega\|_{H^{1}(\Omega)}.$

######  Theorem A.4.

There exists a constant $\beta>0$ such that

(A.5)
$\beta\leq\inf_{0\not=v\in\tilde{V}^{R}_{h}}\sup_{0\neq\omega\in\tilde{\Sigma}_{h}^{R}}\frac{({\mathop{\mathrm{div}\,}}\omega,v)}{\|\omega\|_{H({\mathop{\mathrm{div}\,}},\Omega)}\|v\|_{L^{2}(\Omega)}}.$

The mixed method with the second reduced element finds
$\tilde{\sigma}^{R}_{h}\in\tilde{\Sigma}^{R}_{h}$ and
$\tilde{u}^{R}_{h}\in\tilde{V}^{R}_{h}$ satisfying

(A.6a)
$\displaystyle(\mathcal{A}\tilde{\sigma}^{R}_{h},\tau)+(\tilde{u}^{R}_{h},{\mathop{\mathrm{div}\,}}\tau)$
$\displaystyle=0\qquad$ $\displaystyle{\text{ for all
}}\tau\in\tilde{\Sigma}^{R}_{h},$ (A.6b)
$\displaystyle({\mathop{\mathrm{div}\,}}\tilde{\sigma}^{R}_{h},v)$
$\displaystyle=(f,v)\qquad$ $\displaystyle{\text{ for all
}}v\in\tilde{V}^{R}_{h}.$

The error analysis of this method is more involved than the previous cases
since
${\mathop{\mathrm{div}\,}}(\tilde{\Pi}^{R}\sigma-\tilde{\sigma}^{R}_{h})$ is
not necessarily zero. We proceed by proving a stability result for the method,
as in [7], using the bilinear form

$B(\omega_{h},v_{h};\tau_{h},w_{h}):=(\mathcal{A}\omega_{h},\tau_{h})+(v_{h},{\mathop{\mathrm{div}\,}}\tau_{h})-({\mathop{\mathrm{div}\,}}\omega_{h},w_{h}).$

######  Theorem A.5.

There exists constants $c_{1}$ and $c_{2}$ such that for every
$\omega_{h}\in\tilde{\Sigma}^{R}_{h}$ and $v_{h}\in\tilde{V}_{h}^{R}$, there
exists $\tau_{h}\in\tilde{\Sigma}^{R}_{h}$ and $w_{h}\in\tilde{V}_{h}^{R}$,
satisfying

(A.7)
$\displaystyle\|\omega_{h}\|_{\mathcal{A}}^{2}+\|\tilde{P}^{R}{\mathop{\mathrm{div}\,}}\omega_{h}\|_{L^{2}(\Omega)}^{2}+\|v_{h}\|_{L^{2}(\Omega)}^{2}\leq
c_{1}B(\omega_{h},v_{h};\tau_{h},w_{h}),$

and

$\displaystyle\|\tau_{h}\|_{\mathcal{A}}^{2}+\|\tilde{P}^{R}{\mathop{\mathrm{div}\,}}\tau_{h}\|_{L^{2}(\Omega)}^{2}+\|w_{h}\|_{L^{2}(\Omega)}^{2}$
(A.8) $\displaystyle\leq
c_{2}(\|\omega_{h}\|_{\mathcal{A}}^{2}+\|\tilde{P}^{R}{\mathop{\mathrm{div}\,}}\omega_{h}\|_{L^{2}(\Omega)}^{2}+\|v_{h}\|_{L^{2}(\Omega)}^{2}).$

###### Proof.

By the proof of Theorem A.4 there exists $\rho_{h}\in\tilde{\Sigma}^{R}_{h}$
such that

(A.9) $\displaystyle\|v_{h}\|_{L^{2}(\Omega)}^{2}=$
$\displaystyle({\mathop{\mathrm{div}\,}}\rho_{h},v_{h}),$ (A.10)
$\displaystyle\|\rho_{h}\|_{H({\mathop{\mathrm{div}\,}},\Omega)}\leq$
$\displaystyle\kappa\|v_{h}\|_{L^{2}(\Omega)}.$

Choose $\tau_{h}=\omega_{h}+\frac{1}{\gamma\kappa^{2}}\rho_{h}$ and
$w_{h}=v_{h}-\tilde{P}^{R}{\mathop{\mathrm{div}\,}}\omega_{h}$. Then, we have

$\displaystyle B(\omega_{h},v_{h};\tau_{h},w_{h})=$
$\displaystyle\|\omega_{h}\|_{\mathcal{A}}^{2}+\frac{1}{\gamma\kappa^{2}}(\mathcal{A}\omega_{h},\rho_{h})+(v_{h},{\mathop{\mathrm{div}\,}}\omega_{h})+\frac{1}{\gamma\kappa^{2}}(v_{h},{\mathop{\mathrm{div}\,}}\rho_{h})$
$\displaystyle-({\mathop{\mathrm{div}\,}}\omega_{h},v_{h})+({\mathop{\mathrm{div}\,}}\omega_{h},\tilde{P}^{R}{\mathop{\mathrm{div}\,}}\omega_{h})$
$\displaystyle=\|\omega_{h}\|_{\mathcal{A}}^{2}+\|\tilde{P}^{R}{\mathop{\mathrm{div}\,}}\omega_{h}\|_{L^{2}(\Omega)}^{2}+\frac{1}{\gamma\kappa^{2}}\|v_{h}\|_{L^{2}(\Omega)}^{2}+\frac{1}{\gamma\kappa^{2}}(\mathcal{A}\omega_{h},\rho_{h}),$

where we used (A.9). On the other hand, by (A.10)

$\displaystyle\frac{1}{\gamma\kappa^{2}}(\mathcal{A}\omega_{h},\rho_{h})$
$\displaystyle\geq-\frac{1}{2}\|\omega_{h}\|_{\mathcal{A}}^{2}-\frac{1}{2\gamma^{2}\kappa^{4}}\|\rho_{h}\|_{\mathcal{A}}^{2}$
$\displaystyle\geq-\frac{1}{2}\|\omega_{h}\|_{\mathcal{A}}^{2}-\frac{1}{2\gamma\kappa^{2}}\|v_{h}\|_{L^{2}(\Omega)}^{2}.$

Hence,

$\displaystyle\frac{1}{2}\|\omega_{h}\|_{\mathcal{A}}^{2}+\|\tilde{P}^{R}{\mathop{\mathrm{div}\,}}\omega_{h}\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\gamma\kappa^{2}}\|v_{h}\|_{L^{2}(\Omega)}^{2}\leq
B(\omega_{h},v_{h};\tau_{h},w_{h}).$

This shows (A.7). We clearly have (A.8) if we use (A.10). ∎

Using the above stability result we can prove an a-priori error estimate.

######  Theorem A.6.

Let $\sigma,u$ solve (2.6) and $\tilde{\sigma}^{R}_{h},\tilde{u}^{R}_{h}$
solve (A.6), then

$\displaystyle\|\tilde{\Pi}^{R}\sigma-\tilde{\sigma}^{R}_{h}\|_{\mathcal{A}}+\|\tilde{P}^{R}{\mathop{\mathrm{div}\,}}(\tilde{\Pi}^{R}\sigma-\tilde{\sigma}^{R}_{h})\|_{L^{2}(\Omega)}^{2}+\|\tilde{P}^{R}u-\tilde{u}^{R}_{h}\|_{L^{2}(\Omega)}$
$\displaystyle\leq
C(\|\tilde{\Pi}^{R}\sigma-\sigma\|_{\mathcal{A}}+\|\tilde{P}^{R}u-u\|_{L^{2}(\Omega)}).$

###### Proof.

By Theorem A.5, there exist $\tau_{h}\in\tilde{\Sigma}^{R}_{h}$ and
$w_{h}\in\tilde{V}_{h}^{R}$, satisfying

$\displaystyle\|\tilde{\Pi}^{R}\sigma-\tilde{\sigma}^{R}_{h}\|_{\mathcal{A}}^{2}+\|\tilde{P}^{R}{\mathop{\mathrm{div}\,}}(\tilde{\Pi}^{R}\sigma-\tilde{\sigma}^{R}_{h})\|_{L^{2}(\Omega)}^{2}+\|\tilde{P}^{R}u-\tilde{u}^{R}_{h}\|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq
c_{1}B(\tilde{\Pi}^{R}\sigma-\tilde{\sigma}^{R}_{h},\tilde{P}^{R}u-\tilde{u}^{R}_{h};\tau_{h},w_{h})$

and

$\displaystyle\|\tau_{h}\|_{\mathcal{A}}^{2}+\|\tilde{P}^{R}{\mathop{\mathrm{div}\,}}\tau_{h}\|_{L^{2}(\Omega)}^{2}+\|w_{h}\|_{L^{2}(\Omega)}^{2}$
(A.11) $\displaystyle\leq
c_{2}(\|\tilde{\Pi}^{R}\sigma-\tilde{\sigma}^{R}_{h}\|_{\mathcal{A}}^{2}+\|\tilde{P}^{R}{\mathop{\mathrm{div}\,}}(\tilde{\Pi}^{R}\sigma-\tilde{\sigma}^{R}_{h})\|_{L^{2}(\Omega)}^{2}+\|\tilde{P}^{R}u-\tilde{u}^{R}_{h}\|_{L^{2}(\Omega)}^{2}).$

Then,

$\displaystyle\|\tilde{\Pi}^{R}\sigma-\tilde{\sigma}^{R}_{h}\|_{\mathcal{A}}^{2}+\|\tilde{P}^{R}{\mathop{\mathrm{div}\,}}(\tilde{\Pi}^{R}\sigma-\tilde{\sigma}^{R}_{h})\|_{L^{2}(\Omega)}^{2}+\|\tilde{P}^{R}u-\tilde{u}^{R}_{h}\|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq
c_{1}\Big{(}(\mathcal{A}(\tilde{\Pi}^{R}\sigma-\tilde{\sigma}^{R}_{h},\tau_{h})+(\tilde{P}^{R}u-\tilde{u}^{R}_{h},{\mathop{\mathrm{div}\,}}\tau_{h})\Big{)}$

where we used (2.6b), (A.6b) and (A.3) to say the term
$({\mathop{\mathrm{div}\,}}(\tilde{\Pi}^{R}\sigma-\tilde{\sigma}^{R}_{h}),w_{h})$
vanishes.

On the other hand, using (2.6a), (A.6a) we get

$\displaystyle(\mathcal{A}(\tilde{\Pi}^{R}\sigma-\tilde{\sigma}^{R}_{h}),\tau_{h})+(\tilde{P}^{R}u-\tilde{u}^{R}_{h},{\mathop{\mathrm{div}\,}}\tau_{h})=(\mathcal{A}(\tilde{\Pi}^{R}\sigma-\sigma),\tau_{h})+(\tilde{P}^{R}u-u,{\mathop{\mathrm{div}\,}}\tau_{h}).$

The result now follows by using (A). ∎

As a corollary, we give an error estimate with convergence rate.

######  Corollary A.7.

Let $\sigma,u$ solve (2.6) and $\tilde{\sigma}^{R}_{h},\tilde{u}^{R}_{h}$
solve (A.6). If $\sigma\in H^{1}(\Omega,\mathbb{S}),u\in[H^{1}(\Omega)]^{3}$,
then

(A.12)
$\|\sigma-\tilde{\sigma}^{R}_{h}\|_{\mathcal{A}}+\|u-\tilde{u}^{R}_{h}\|_{L^{2}(\Omega)}\leq
Ch(\|\sigma\|_{H^{1}(\Omega)}+\|u\|_{H^{1}(\Omega)}).$

######  Remark A.8.

The argument in the proof of Theorem 2.11 (on robustness near
incompressibility) does not extend to this reduced space pair because an
analogue of (4.8) is not available for $\tilde{\sigma}_{h}^{R}$ and
$\tilde{\Pi}^{R}$.

## References

  * [1] P. Alfeld.  A trivariate Clough–Tocher scheme for tetrahedral data.  Computer Aided Geometric Design, 1(2):169–181, 1984.
  * [2] D. Arnold, G. Awanou, and R. Winther.  Finite elements for symmetric tensors in three dimensions.  Mathematics of Computation, 77(263):1229–1251, 2008.
  * [3] D. Arnold, R. Falk, and R. Winther.  Mixed finite element methods for linear elasticity with weakly imposed symmetry.  Mathematics of Computation, 76(260):1699–1723, 2007.
  * [4] D. N. Arnold.  Finite element exterior calculus.  SIAM, 2018.
  * [5] D. N. Arnold, F. Brezzi, and J. Douglas, Jr.  PEERS: a new mixed finite element for plane elasticity.  Japan J. Appl. Math., 1(2):347–367, 1984.
  * [6] D. N. Arnold, J. Douglas, and C. P. Gupta.  A family of higher order mixed finite element methods for plane elasticity.  Numerische Mathematik, 45:1–22, 1984.
  * [7] D. N. Arnold, R. S. Falk, and R. Winther.  Finite element exterior calculus, homological techniques, and applications.  Acta numerica, 15:1–155, 2006.
  * [8] D. N. Arnold and K. Hu.  Complexes from complexes.  Found. Comput. Math., 21(6):1739–1774, 2021.
  * [9] D. N. Arnold and R. Winther.  Mixed finite elements for elasticity.  Numerische Mathematik, 92(3):401–419, 2002.
  * [10] M. Bebendorf.  A note on the Poincaré inequality for convex domains.  Z. Anal. Anwendungen, 22(4):751–756, 2003.
  * [11] D. Boffi, F. Brezzi, and M. Fortin.  Reduced symmetry elements in linear elasticity.  Commun. Pure Appl. Anal., 8(1):95–121, 2009.
  * [12] D. Boffi, F. Brezzi, and M. Fortin.  Mixed finite element methods and applications, volume 44 of Springer Series in Computational Mathematics.  Springer, Heidelberg, 2013.
  * [13] E. Calabi.  On compact, Riemannian manifolds with constant curvature I.  Proc. Sympos. Pure Math., Amer. Math. Soc., Providence, RI, 1961, pages 155–180, 1961.
  * [14] S. Christiansen, J. Gopalakrishnan, J. Guzmán, and K. Hu.  A discrete elasticity complex on three-dimensional Alfeld splits.  Numer. Math., 2023.
  * [15] P. G. Ciarlet.  Mathematical elasticity. Vol. I, volume 20 of Studies in Mathematics and its Applications.  North-Holland Publishing Co., Amsterdam, 1988.  Three-dimensional elasticity.
  * [16] P. G. Ciarlet, L. Gratie, and C. Mardare.  Intrinsic methods in elasticity: a mathematical survey.  Discrete & Continuous Dynamical Systems-A, 23(1&2):133, 2009.
  * [17] M. Costabel and A. McIntosh.  On Bogovskiĭ and regularized Poincaré integral operators for de Rham complexes on Lipschitz domains.  Mathematische Zeitschrift, 265(2):297–320, 2010.
  * [18] M. Eastwood.  A complex from linear elasticity.  In Proceedings of the 19th Winter School” Geometry and Physics”, pages 23–29. Circolo Matematico di Palermo, 2000.
  * [19] L. C. Evans.  Partial differential equations, volume 19 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, 1998.
  * [20] R. S. Falk.  Finite element methods for linear elasticity.  In Mixed Finite Elements, Compatibility Conditions, and Applications, volume 1939 of Lecture Notes in Mathematics, pages 160–194. Springer Verlag, 2008.
  * [21] M. Farhloul and M. Fortin.  Dual hybrid methods for the elasticity and the Stokes problems: a unified approach.  Numer. Math., 76(4):419–440, 1997.
  * [22] G. Fu, J. Guzman, and M. Neilan.  Exact smooth piecewise polynomial sequences on Alfeld splits.  Mathematics of Computation, 2020.
  * [23] S. Gong, J. Gopalakrishnan, J. Guzmán, and M. Neilan.  Discrete elasticity exact sequences on Worsey-Farin splits.  ESAIM Math. Model. Numer. Anal., 57(6):3373–3402, 2023.
  * [24] J. Gopalakrishnan and J. Guzmán.  A second elasticity element using the matrix bubble.  IMA Journal of Numerical Analysis, 32(1):352–372, 2012.
  * [25] J. Guzmán and M. Neilan.  Inf-sup stable finite elements on barycentric refinements producing divergence–free approximations in arbitrary dimensions.  SIAM Journal on Numerical Analysis, 56(5):2826–2844, 2018.
  * [26] J. Hu and S. Zhang.  A family of symmetric mixed finite elements for linear elasticity on tetrahedral grids.  Sci. China Math., 58(2):297–307, 2015.
  * [27] J. Hu and S. Zhang.  Finite element approximations of symmetric tensors on simplicial grids in $\mathbb{R}^{n}$: the lower order case.  Math. Models Methods Appl. Sci., 26(9):1649–1669, 2016.
  * [28] X. Huang, C. Zhang, Y. Zhou, and Y. Zhu.  New low-order mixed finite element methods for linear elasticity.  Adv. Comput. Math., 50(2):Paper No. 17, 2024.
  * [29] C. Johnson and B. Mercier.  Some equilibrium finite element methods for two-dimensional elasticity problems.  Numerische Mathematik, 30(1):103–116, 1978.
  * [30] M. Křížek.  An equilibrium finite element method in three-dimensional elasticity.  Apl. Mat., 27(1):46–75, 1982.  With a loose Russian summary.
  * [31] P. L. Lederer and R. Stenberg.  Energy norm analysis of exactly symmetric mixed finite elements for linear elasticity.  Math. Comp., 92(340):583–605, 2023.
  * [32] J. J. Lee.  Towards a unified analysis of mixed methods for elasticity with weakly symmetric stress.  Adv. Comput. Math., 42(2):361–376, 2016.
  * [33] J. E. Marsden and T. J. R. Hughes.  Mathematical foundations of elasticity.  Dover Publications Inc., New York, 1994.
  * [34] M. E. Morley.  A family of mixed finite elements for linear elasticity.  Numer. Math., 55(6):633–666, 1989.
  * [35] R. Stenberg.  A family of mixed finite elements for the elasticity problem.  Numer. Math., 53(5):513–538, 1988.
  * [36] R. Stenberg.  Postprocessing schemes for some mixed finite elements.  RAIRO Modél. Math. Anal. Numér., 25(1):151–167, 1991.
  * [37] V. B. Watwood and B. J. Hartz.  An equilibrium stress field model for finite element solutions of two-dimensional elastostatic problems.  International Journal of Solids and Structures, 4:857–873, 1968\.
  * [38] S. Zhang.  A new family of stable mixed finite elements for the 3D Stokes equations.  Math. Comp., 74(250):543–554, 2005.