# Nonlinear elasticity complex and a finite element diagram chase

Kaibo Hu  Mathematical Institute, University of Oxford, Andrew Wiles Building,
Radcliffe Observatory Quarter, Oxford, OX2 6GG, UK<EMAIL_ADDRESS>
(Date: . Manuscript prepared for proceedings of the INdAM conference
“Approximation Theory and Numerical Analysis meet Algebra, Geometry,
Topology”, which was held in September 2022 at Cortona, Italy.)

###### Abstract.

In this paper, we present a nonlinear version of the linear elasticity
(Calabi, Kröner, Riemannian deformation) complex which encodes isometric
embedding, metric, curvature and the Bianchi identity. We reformulate the
rigidity theorem and a fundamental theorem of Riemannian geometry as the
exactness of this complex. Then we generalize an algebraic approach for
constructing finite elements for the Bernstein-Gelfand-Gelfand (BGG)
complexes. In particular, we discuss the reduction of degrees of freedom with
injective connecting maps in the BGG diagrams. We derive a strain complex in
two space dimensions with a diagram chase.

##  1\. Introduction

A complex refers to a sequence of spaces $V^{k},k=0,\cdots,n$ and operators
$d^{k}:V^{k}\to V^{k+1},k=0,\cdots,n-1$, such that $d^{k+1}\circ d^{k}=0$ for
any $k$. In this paper, we focus on differential complexes, where
$d^{{\scriptscriptstyle\bullet}}$ are differential operators. Differential
complexes play an important role in geometry, analysis and numerical partial
differential equations [6, 8, 9, 10, 15, 17]. Finite element exterior calculus
(FEEC) [6, 8, 9] was inspired by finite element methods for electromagnetic
and elasticity problems. A key idea of FEEC is to discretize problems with
spaces that fit in complexes. For electromagneitc problems, the de Rham
complex encodes the key structures of the continuous problem; while for
elasticity problems, the elasticity (Calabi, Kröner, Riemannian deformation)
complex [7, 27, 28, 44] plays a similar role.

Most existing literature focuses on complexes with linear differential
operators. With linear operators, cohomology can be defined and carries rich
information (for example, finite dimensional cohomology further implies
various analytic results [10]). In some cases, e.g., the elasticity complex,
the operators are obtained as a linearization of nonlinear versions. The space
and nonlinear operators also form a complex. Nevertheless, these original
complexes with nonlinear operators (referred to as nonlinear complexes below)
have received less attention.

It is of interest to investigate the nonlinear complexes, as they

  1. (1)

incorporate important structures for nonlinear models and PDEs;

  2. (2)

provide a motivation for the linear versions;

  3. (3)

may play a role in clarifying some problems in numerical analysis of nonlinear
problems (e.g., the stability/instability of different choices of finite
elements for nonlinear elasticity [2, 13]).

Moreover, concepts are more clear in the context of nonlinear elasticity, such
as pullbacks and frame indifference. In the first part of this paper, we
present the nonlinear version of the elasticity complex and show that the
exactness of this complex corresponds to rigidity and a fundamental theorem of
Riemannian geometry [23] (thus clarifying (1) and (2) above).

In the second part of the paper, we generalize an approach for constructing
finite element versions of the linear elasticity complex (the strain and
stress complexes in two space dimensions). The elasticity complex is a special
case of the Bernstein-Gelfand-Gelfand (BGG) sequences [10, 16, 17, 28]. The
BGG machinery establishes a link between de Rham complexes and more
complicated ones, e.g., the elasticity complex [10, 16]. Specifically, one
connects several copies of de Rham complexes by algebraic connecting maps and
eliminates some components from the diagram. On the discrete level, one may
mimic the construction and fit finite element spaces in diagrams. We will
refer to this approach of deriving more finite elements from de Rham complexes
with various kinds of regularity as a finite element diagram chase. The main
difficulty of this approach is to construct finite element de Rham complexes
that fulfil the requirement, as the input complexes should have different
kinds of regularity. This can be seen in the simplest case, where one connects
two 1D de Rham complexes and get a BGG complex with a second order operator
(see Figure 1 and Figure 2). In higher dimensions, the diagram chase requires
an input of more delicate finite element de Rham complexes.

Figure 1. BGG diagram for deriving a complex with a second order operator. The
two rows are finite element de Rham complexes with different kinds of
regularity. Figure 2. The derived BGG complex with a second order operator.

A finite element diagram chase was first used in [7] as a re-interpretation of
the Arnold-Winther elasticity element. Recently, this approach has been
extended to other complexes and finite elements [19, 20, 22]. In applications
related to the Hellinger-Reissner principle (e.g., [7]), the connecting maps
are surjective. In this case, one may define a discrete version of the
connecting map by rewriting the degrees of freedom [7]. In other examples, the
connecting maps can be injective. The reduction of spaces and degrees of
freedom is more subtle. A finite element elasticity complex on the Alfeld
split in 3D was constructed in [19]. The first connecting map
$\operatorname{mskw}$ is injective. A discrete BGG diagram plays an important
role in [19]. Nevertheless, the strain space was derived directly, rather than
following a systematic diagram chase.

In this paper, we demonstrate the reduction of spaces and degrees of freedom
by closely following the BGG diagrams. We use the 2D strain complex as an
example. The derived space is not symmetric but is so in a discrete sense (see
Section 3 for more detailed explanations). We also discuss a stress complex on
quadrilateral grids.

Discretizing the entire BGG diagrams is not merely a way to derive discrete
BGG sequences, but it may have the benefit of incorporating more physics. With
each BGG diagram, we can gather spaces from all rows into a twisted complex
where the operators have an off-diagonal part mixing the components from
different rows [10, 16]. The twisted complexes encode more information than
the BGG sequences, as the latter can be obtained by eliminating some spaces
from the former (and the main conclusion of the BGG machinery is that this
elimination is cohomology-preserving). This elimination has a physical
meaning. For example, for the elasticity complex, the Hodge-Laplacian problem
of the twisted complex corresponds to a linear Cosserat model which involves
both displacement and a pointwise rotation. Then to derive the elasticity
complex, one eliminates the rotational degree of freedom and obtains standard
elasticity with only displacement as the main unknown. In [16], it was
observed that such a correspondence exists in 2D for the Reissner-Mindlin
plate and the Kirchhoff-Love plate, and in 3D for the Cosserat and standard
elasticity. Here we further observe that the Hodge-Laplacian problems for the
1D twisted and BGG complexes (continuous versions of Figure 1 and Figure 2),
with the energy functionals
$\|\partial_{x}u-w\|_{A}^{2}+\|\partial_{x}w\|_{B}^{2}$ and
$\|\partial_{x}^{2}w\|^{2}_{C}$, correspond to the linearized Timoshenko beam
and Euler-Bernoulli beam, respectively. Here $u$ is the displacement; $w$ is
the deflection; $\|\cdot\|_{{\scriptscriptstyle\bullet}}$ are weighted norms.
Therefore deriving a discrete version of the entire diagram may also shed a
light on these problems.

Differential complexes also provide a connection between scalar splines and
tensor-valued spaces with relatively lower regularity. In particular,
exactness of a complex
$(V^{{\scriptscriptstyle\bullet}},d^{{\scriptscriptstyle\bullet}})$ gives an
equality of dimensions $\sum_{k}(-1)^{k}\dim V^{k}=0$. This connection also
inspires the construction of tensor-valued finite elements. Specifically, one
may also start with a smooth scalar spline space in a complex, differentiate
it and obtain tensor-valued spaces that complete the sequence. Table 1
summarizes some scalar and tensor-valued spaces and complexes that connect
them.

The rest of the paper is organized as the following. In Section 2, we discuss
a complex with nonlinear differential operators. In Section 3, we derive
finite element strain and stress complexes following a diagram chase. In
Section 4, we provide concluding remarks.

2D Stokes (2013) [29] | Argyris (1968) [4] | Falk-Neilan (2013) [29]
---|---|---
2D Stokes (2018) [21] | Clough-Tocher (1965) [26] | Arnold-Qin (1992) [11]
2D Stokes (2020) [32] | Powell-Sabin (1977) [49] | Zhang (2008) [49]
2D elasticity (2002) [12] | Argyris (1968) [4] | Arnold-Winther (2002) [12]
2D elasticity (2018) [20] | Argyris (1968) [4] | Hu-Zhang (2014) [34]
2D elasticity (2022) [22] | Clough-Tocher (1965) [26] | Johnson-Mercier (1978) [37]
3D Stokes (2015) [42] | Zěníšek (1973) [47] | Neilan (2015) [42]
3D Stokes (2020) [30] | Alfeld (1984) [1] | Zhang (2005) [48]
3D Stokes (2022)[33] | Worsey-Farin (1987) [46] | Zhang (2011) [50]
3D elasticity (2020) [19] | Alfeld (1984) [1] |  | Christiansen-Gopalakrishnan
---
-Guzman-Hu (2020) [19]
Table 1. Complexes connect scalar/vector splines with tensor-valued spaces.
The second column of the table includes scalar/vector spline spaces, and the
third contains the corresponding tensor-valued spaces. The first column
indicates the complexes that contain these spaces. See also [43] for a review
of finite element Stokes complexes and [35] for some BGG complexes.

##  2\. Nonlinear complex

In this section, we present a nonlinear version of the linear elasticity
complex.

Equip $\mathbb{R}^{3}$ with a flat metric $g_{0}$. Let
$\Omega\subset\mathbb{R}^{3}$ be a bounded domain, and
$\varphi_{0}:\Omega\to\varphi_{0}(\Omega)\subset\mathbb{R}^{3}$ be a given
invertible smooth vector-valued function. We refer to the following sequence
as a nonlinear elasticity complex:

(2.1)
${\mathbb{RM}}$${C^{\infty}({\Omega};\mathbb{R}^{3})}$${\Gamma_{+}(S^{2}T^{\ast}\Omega)}$${\Gamma(S^{2}T^{\ast}\Omega)}$${T^{\ast}\Omega}$${0.}$$\scriptstyle{\subset}$$\scriptstyle{\iota}$$\scriptstyle{\mathrm{Ric}}$$\scriptstyle{\mathrm{Bian}}$

Here $\mathbb{RM}$ is the space of rigid body motions composed with
$\varphi_{0}$, i.e.,

$\mathbb{RM}:=\\{c+Q\cdot\varphi_{0}(x),~{}c\in\mathbb{R}^{3},Q\in SO(3)\\}.$

The space $C^{\infty}({\Omega};\mathbb{R}^{3})$ has the geometric meaning of
embedding of $\Omega$ into $\mathbb{R}^{3}$. Each such map $\varphi$ induces a
pullback of metric
$\varphi^{\ast}:\Gamma(S^{2}T^{\ast}\mathbb{R}^{3})\to\Gamma(S^{2}T^{\ast}\Omega)$,
where $\Gamma(S^{2}T^{\ast}M)$ denotes the section of symmetric $(0,2)$ tensor
fields on a manifold $M$. The operator $\iota$ denotes the change of the
metric, i.e.,
$\iota(\varphi):=(\varphi)^{\ast}g_{0}-(\varphi_{0})^{\ast}g_{0}$, and in
coordinate forms, $\iota(\varphi):=D\varphi\cdot
g_{0}\cdot(D\varphi)^{t}-D\varphi_{0}\cdot g_{0}\cdot(D\varphi_{0})^{t}$,
where $(D\varphi)_{ij}:=\partial_{i}\varphi_{j}$. For any $r\in\mathbb{RM}$,
$r$ induces the same map on the metric as $\varphi_{0}$, as $r$ is a
composition of $\varphi_{0}$ with a rigid body motion. Then it is readily seen
that $\iota(r)=0$ for any $r\in\mathbb{RM}$. The space
$\Gamma_{+}(S^{2}T^{\ast}\Omega)$ is defined by

$\Gamma_{+}(S^{2}T^{\ast}\Omega):=\\{g\in\Gamma(S^{2}T^{\ast}\Omega):g+(\varphi_{0})^{\ast}g_{0}>0\\}.$

Here a matrix $g+(\varphi_{0})^{\ast}g_{0}>0$ means that
$g+(\varphi_{0})^{\ast}g_{0}$ is positive definite. In (2.1), $\mathrm{Ric}$
is the operator mapping a metric to the Ricci curvature. In three dimensions,
the Riemannian, Ricci and Einstein tensors encode the same information.
Moreover, $\mathrm{Bian}$ is the operator in the (contracted) differential
Bianchi identity $R_{ij}\mapsto\nabla^{i}(R_{ij}-\frac{1}{2}Rg_{ij})$, where
$R_{ij}$ is the Ricci tensor, $R:=g^{ij}R_{ij}$ is the scalar curvature, and
$g_{ij}$ is the metric tensor. Note that $R_{ij}-\frac{1}{2}Rg_{ij}$ is the
Einstein tensor.

In continuum mechanics, $\varphi$ represents a map from the reference
configuration $\Omega$ to the current configuration $\varphi(\Omega)$. When
$g_{0}=I$ is the Euclidean metric of $\mathbb{R}^{3}$,
$\varphi^{\ast}(g_{0})$, or in the matrix form, $D\varphi\cdot(D\varphi)^{t}$,
is the right Cauchy-Green tensor [41, Definition 3.5].

######  Theorem 1.

The sequence (2.1) is a complex in the sense that composing two successive
operators gives zero.

###### Proof.

We have $\mathrm{Ric}\circ\iota=0$ as $\varphi_{0}$ is a flat metric by
assumption, and $\operatorname{Bian}\circ\operatorname{Ric}=0$ is the Bianchi
identity. ∎

######  Remark 1.

Another complex related to nonlinear elasticity was investigated in [3]. The
complex includes the embedding $\varphi$ and the deformation tensor
$D\varphi$, and therefore corresponds to a de Rham complex. The operators in
the complex in [3] are linear. Thus structures of multiplication are not
explicitly included in the sequences. This is different from (2.1), which
contains $D\varphi\cdot(D\varphi)^{t}$.

In the rest of this section, we discuss exactness of the nonlinear complex and
show how it corresponds to some classical results. For simplicity, hereafter
we assume that $g_{0}$ is the Euclidean metric, or in coordinates, $g_{0}=I$.

The exactness at index zero (at the space
$C^{\infty}({\Omega};\mathbb{R}^{3})$) states that if
$D\varphi\cdot(D\varphi)^{T}-D\varphi_{0}\cdot(D\varphi_{0})^{T}=0$, then
$\varphi$ is in $\mathbb{RM}$, i.e., there exists $c,Q$ such that
$\varphi=c+Q\varphi_{0}$. This is exactly the rigidity theorem [23, Theorem
8.7-1] which is stated for $C^{1}$ maps and connected open subsets of
$\mathbb{R}^{n}$. The connectedness is a natural condition, as in the linear
case exactness relies on trivial topology.

The exactness at index one states that for any
$g\in\Gamma_{+}(S^{2}T^{\ast}\Omega)$, if the Ricci (Riemann, Einstein) tensor
vanishes, then there exists $\varphi_{0}$ such that
$g=D\varphi\cdot(D\varphi)^{t}-D\varphi_{0}\cdot(D\varphi_{0})^{t}$. This is
equivalent to applying the theorem on “the existence of an immersion with a
prescribed metric tensor” [23, Theorem 8.6-1] to
$g+D\varphi_{0}\cdot(D\varphi_{0})^{t}$. In fact,
$g\in\Gamma_{+}(S^{2}T^{\ast}\Omega)$ implies that
$g+D\varphi_{0}\cdot(D\varphi_{0})^{t}$ is positive definite. Moreover, $g$ is
flat implies that $g+D\varphi_{0}\cdot(D\varphi_{0})^{T}$ also has vanishing
Riemnnian tensor. Then we have all the conditions for applying [23, Theorem
8.6-1], which asserts that the exactness holds if $\Omega$ is a simply
connected open set in $\mathbb{R}^{n}$ and $g\in C^{2}$ (therefore the
embedding $\phi$ in $C^{3}$). This again corresponds to the topological
conditions for the exactness of the linear elasticity complex.

As shown above, the exactness of (2.1) at index zero corresponds to the
rigidity and at index one is a special case of the fundamental theorem of
Riemannian geometry. It is tempting to investigate the exactness at index two.
Generally speaking, the question is that assuming that a tensor satisfies the
symmetries of a Ricci (Riemannian) tensor, is it a curvature tensor of some
metric? However, this is not a well-posed question, as to formulate the
Bianchi identity, one needs a metric given beforehand. This issue was
discussed in [45], where the author investigated another possibility of
imposing the compatibility condition as “there exists some metric $g$, such
that $\operatorname{Bian}_{g}\sigma=0$”, where $\operatorname{Bian}_{g}$ is
the Bianchi operator defined by $g$. Nevertheless, with this condition,
$\sigma$ may not be the Riemannian tensor of any metric, i.e., the “exactness”
does not hold.

The linearization of (2.1) around the Euclidean metric gives the linear
elasticity complex

(2.2)
${\mathrm{RM}}$${C^{\infty}\otimes\mathbb{V}}$${C^{\infty}\otimes\mathbb{S}}$${C^{\infty}\otimes\mathbb{S}}$${C^{\infty}\otimes\mathbb{V}}$${0,}$$\scriptstyle{{\operatorname{def}}}$$\scriptstyle{\operatorname{inc}}$$\scriptstyle{\operatorname{div}}$

where we write $\mathbb{V}=\mathbb{R}^{3}$ for the space of vectors, and
$C^{\infty}\otimes\mathbb{V}=C^{\infty}(\Omega)\otimes\mathbb{V}$ is the space
of smooth vector fields. Similarly, $\mathbb{S}$ and
$C^{\infty}\otimes\mathbb{S}$ denote the spaces of $3\times 3$ symmetric
matrices and smooth matrix fields, respectively. Here $\mathrm{RM}$ is the
space of infinitesimal rigid body motions

$\mathrm{RM}:=\\{c+K\cdot x,~{}c\in\mathbb{R}^{3},Q\in{SO}(3)\\}=\\{c+b\times
x,~{}c,b\in\mathbb{R}^{3}\\}.$

The operators in the elasticity complex are the deformation or linearized
strain operator $\operatorname{def}=\operatorname{sym}\operatorname{grad}$,
the incompatibility operator
$\operatorname{inc}=\operatorname{curl}\circ\mathrm{T}\circ\operatorname{curl}$
(where $\mathrm{T}$ denotes the transpose operation and $\operatorname{curl}$
acts on a matrix field by columns), and the column-wise divergence operator
map matrices to vectors.

##  3\. Finite element diagram chase

In this section, we present a diagram chase to construct a finite element
version of elasticity complexes in 2D.

Following [10], we introduce some notations in 2D. We use $\mathbb{R}$,
$\mathbb{V}$, $\mathbb{M}$, $\mathbb{S}$, and $\mathbb{K}$ to denote the
(finite dimensional) spaces of scalars, vectors, matrices, symmetric matrices
and skew-symmetric matrices, respectively. That is,
$\mathbb{V}:=\mathbb{R}^{2}$, $\mathbb{M}:=\mathbb{R}^{2\times 2}$, and
$\mathbb{S}:=\mathbb{R}^{2\times 2}_{\operatorname{sym}}$. The operators
$\operatorname{skw}:\mathbb{M}\to\mathbb{K}$ and
$\operatorname{sym}:\mathbb{M}\to{\mathbb{S}}$ are the skew-symmtrization and
symmetrization, respectively. In $\mathbb{R}^{2}$, a skew symmetric matrix can
be identified with a scalar. Let $\operatorname{mskw}:\mathbb{R}\to\mathbb{K}$
be this identification, i.e.,

$\operatorname{mskw}(u):=\left(\begin{array}[]{cc}0&u\\\
-u&0\end{array}\right)\quad\mbox{in }\mathbb{R}^{2}.$

We also let
$\operatorname{sskw}=\operatorname{mskw}^{-1}\circ\operatorname{skw}:\mathbb{M}\to\mathbb{R}$
be the map taking the skew part of a matrix and identifying it with a scalar.

We use $H^{q}$ to denote the Sobolev space for which all derivatives of order
less than or equal to $q$ are square integrable. For a linear operator $D$,
define $H(D):=\\{u\in L^{2}:Du\in L^{2}\\}$.

###  3.1. Stress and strain complexes: an overview

In 2D, there are two complexes related to (2.2). Related to the Hellinger-
Reissner principle, one has the stress complex

(3.1)
${0}$${H^{2}}$${H(\operatorname{div},{\mathbb{S}})}$${L^{2}\otimes\mathbb{V}}$${0,}$$\scriptstyle{\operatorname{curl}\operatorname{curl}}$$\scriptstyle{\operatorname{div}}$

where $H(\operatorname{div},{\mathbb{S}}):=\\{\sigma\in
L^{2}\otimes\mathbb{S}:\operatorname{div}\sigma\in L^{2}\otimes\mathbb{V}\\}$.
The stress complex may be viewed as a restriction of the last part of the
complex (2.2) to a two dimensional face (thus $\operatorname{inc}$ becomes
$\operatorname{curl}\operatorname{curl}$ and $\operatorname{div}$ remains).
The stress complex can be derived from the following diagram:

(3.2)
${0}$${H^{2}}$${H^{1}\otimes\mathbb{V}}$${L^{2}}$${0}$${0}$${H^{1}\otimes\mathbb{V}}$${H(\operatorname{div})\otimes\mathbb{V}}$${L^{2}\otimes{\mathbb{V}}}$${0.}$$\scriptstyle{\operatorname{curl}}$$\scriptstyle{\operatorname{div}}$$\scriptstyle{\operatorname{curl}}$$\scriptstyle{I}$$\scriptstyle{\operatorname{div}}$$\scriptstyle{\mathrm{-2\operatorname{sskw}}}$

Applications for the Hellinger-Reissner principle of elasticity use the last
part of the diagram. On the continuous level, the idea of deriving (3.1) from
the diagram (3.2) is to eliminate spaces connected by the algebraic diagonal
maps as much as possible from the diagram. For example, we remove the two
spaces connected by $I$ and the components connected by $\operatorname{sskw}$
(thus $L^{2}$ is removed from the top row and the symmetric part of
$H(\operatorname{div})\otimes\mathbb{V}$ is left in the bottom row). We refer
to [10, 16] for a precise discussion of the BGG machinery.

To derive finite elements, we intend to fit finite element spaces in the
diagram (3.2):

(3.3)
${0}$${W^{0}}$${W^{1}}$${W^{2}}$${0}$${0}$${Y^{0}}$${Y^{1}}$${Y^{2}}$${0.}$$\scriptstyle{\operatorname{curl}}$$\scriptstyle{\operatorname{div}}$$\scriptstyle{\operatorname{curl}}$$\scriptstyle{I}$$\scriptstyle{\operatorname{div}}$$\scriptstyle{\mathrm{-2\operatorname{sskw}}}$

One first constructs a finite element Stokes complex
$W^{{\scriptscriptstyle\bullet}}$ (de Rham complexes with enhanced regularity)
for the first row. Then one chooses $Y^{0}$ to be identical to $W^{1}$, so
that they are matched precisely by the connecting map $I$.

Once we have a diagram with discrete spaces (3.3), a similar argument as the
continuous level leads to a discretization of (3.1). The spaces in the derived
complex are kernels and cokernels of the connecting algebraic maps. The recipe
does not immediately implies the degrees of freedom. Nevertheless, one may
rewrite the degrees of freedom for $Y^{1}$, identify some of them with the
degrees of freedom for $W^{2}$ and set them to zero. The remaining ones give a
set of degrees of freedom for symmetric tensors. We refer to [7, 20] for
details of this argument.

Another version of elasticity complexes in 2D, the strain complex

(3.4)
${0}$${H^{2}\otimes\mathbb{V}}$${H^{1}(\operatorname{rot},{\mathbb{S}})}$${L^{2}}$${0,}$$\scriptstyle{\operatorname{def}}$$\scriptstyle{\operatorname{rot}\operatorname{rot}}$

is a restriction of the first several spaces of the 3D complex (2.2) to a 2D
plane (with higher Sobolev regularity, c.f. [22, (64)]). Here
$H^{1}(\operatorname{rot},{\mathbb{S}}):=\\{u\in
H^{1}\otimes\mathbb{S}:\operatorname{rot}u\in H^{1}\otimes\mathbb{V}\\}$. The
strain complex (3.4) contains displacement and linearized Gauss curvature (see
appendix). The complex (3.4) can be derived from the following diagram

(3.5)
${0}$${H^{2}\otimes\mathbb{V}}$${H^{1}(\operatorname{rot},\mathbb{M})}$${H^{1}\otimes{\mathbb{V}}}$${0}$${0}$${H^{2}}$${H^{1}\otimes{\mathbb{V}}}$${L^{2}\otimes{\mathbb{V}}}$${0,}$$\scriptstyle{\operatorname{grad}}$$\scriptstyle{\operatorname{rot}}$$\scriptstyle{\operatorname{grad}}$$\scriptstyle{\operatorname{mskw}}$$\scriptstyle{\operatorname{rot}}$$\scriptstyle{I}$

where $H^{1}(\operatorname{rot},\mathbb{M}):=\\{e\in
H^{1}\otimes\mathbb{M}:\operatorname{rot}e\in H^{1}\otimes\mathbb{V}\\}$. On
the discrete level, to fit finite element spaces in a discrete version of
(3.5),

(3.6)
${0}$${Z^{0}}$${Z^{1}}$${Z^{2}}$${0}$${0}$${V^{0}}$${V^{1}}$${V^{2}}$${0.}$$\scriptstyle{\operatorname{grad}}$$\scriptstyle{\operatorname{rot}}$$\scriptstyle{\operatorname{grad}}$$\scriptstyle{\operatorname{mskw}}$$\scriptstyle{\operatorname{rot}}$$\scriptstyle{I}$

we may follow the same idea as for the stress complex (3.2) by using identical
spaces for the two slots connected by $I$, i.e., $V^{1}=Z^{2}$. However, the
elimination of $V^{0}$ from $Z^{1}$ is more subtle than the stress complex. To
see this, we first observe that in (3.3), a diagram chase shows that as long
as $\operatorname{sskw}$ maps $Y^{1}$ to $W^{2}$ (which is often the case as
$W^{2}$ requires minimal regularity) and both rows are exact, then
$\operatorname{sskw}$ is also surjective on the discrete level, i.e., $W^{2}$
precisely corresponds to the skew-symmetric part of $Y^{1}$. Contrarily, the
subtlety for (3.5) is that there is no canonical way to choose the two spaces
connected by $\operatorname{mskw}$ such that $\operatorname{mskw}V^{0}$ is
precisely the skew-symmetric part of $Z^{1}$. Actually, in typical examples,
$\operatorname{mskw}$ maps $Z^{0}$ to $V^{1}$, but $V^{0}$ is often too small
compared to the skew-symmetric part of $Z^{1}$.

In fact, this difference between (3.3) and (3.6) is already seen on the
continuous level in the framework of [10]. In [10], a BGG complex is derived
from two input complexes. There exists some $J>0$, such that the connecting
maps $S^{j}$ are injective for $j\geq J$ and surjective for $j\geq J$
(therefore bijective for $j=J$) [10, (19)]. The derivation for the two parts
follows slightly different forms. For example, the derived operators [10,
(26)] contain projections to the cokernel for small $j$ (e.g.,
$\operatorname{sym}\operatorname{grad}=P_{\mathcal{R}(\operatorname{mskw})^{\perp}}\operatorname{grad}$),
while for large $j$ the operators directly follow from the input complexes
(e.g., $\operatorname{div}$). This partly explains the difference and
difficulty in dealing with (3.4) and (3.3) compared to existing results for
(3.1) [7, 20].

###  3.2. A discrete stress complex on quadrilateral meshes

We derive a discrete stress complex from a BGG diagram. The idea is similar to
existing results on other types of meshes [7, 20, 22]. Therefore we will only
give a brief outline. Nevertheless, the resulting complex seems to be new in
the literature.

In the construction, we split a quadrilateral by connecting the opposite
vertices. This leads to four sub-triangles. All the finite element spaces are
defined on this subdivision (referred to as the criss-cross grid; see Figure
3).

$\mathbb{R}$$\operatorname{curl}$+10$\operatorname{div}$+11$I$$\mathbb{R}^{2}$+10$\operatorname{curl}$$\mathrm{sskw}$+16$\operatorname{div}$+8
$\operatorname{curl}\operatorname{curl}$+5$\operatorname{div}$+8 Figure 3.
Diagram chase for a quadrilateral macroelement. We denote the top row by
$(W^{{\scriptscriptstyle\bullet}},d^{{\scriptscriptstyle\bullet}})$ and the
bottom row by
$(Y^{{\scriptscriptstyle\bullet}},d^{{\scriptscriptstyle\bullet}})$, following
the notation in (3.3).

We will use the notation in (3.3) to define the spaces. The sequence in the
top row is studied in [36]. The scalar spline space $W^{0}$ can be found in
[38]. The space consists of piecewise cubic polynomials with continuous
derivatives. The last two spaces in the first row were discussed by Arnold and
Qin [11] as a Stokes pair. Specifically, $W^{1}$ consists of continuous
piecewise quadratic polynomials. Therefore there are 10 interior degrees of
freedom (2 at each of the 5 Lagrange points inside the refined quadrilateral).
It is clear that $\operatorname{div}$ maps $W^{1}$ to piecewise linear
polynomials, but this map is not surjective. In fact, for any $u\in W^{1}$,
$\operatorname{div}u$ satisfies a relation at the intersection of the two
diagonal lines. Figure 4 shows a special configuration where the two lines are
the $x$\- and $y$-axises, respectively. Then one verifies that

$(\operatorname{div}u|_{T_{1}}-\operatorname{div}u|_{T_{2}}+\operatorname{div}u|_{T_{3}}-\operatorname{div}u|_{T_{4}})(z)=0.$

The situation for general configurations where the two lines are not
orthogonal is similar.

Figure 4. A compatibility condition is imposed at each singular vertex $z$.

With this constraint, piecewise linear polynomials have dimension 11 on each
quadrilateral. Vertices joined by two straight lines are referred to as
singular vertices [11]. A similar situation appears in the construction of
finite elements on the Powell-Sabin split, where two straight lines intersect
on edges of the triangulation [32].

The second row
$(Y^{{\scriptscriptstyle\bullet}},d^{{\scriptscriptstyle\bullet}})$ in (3.3)
is a vector-valued (two copies of) finite element de Rham complex consisting
of quadratic Lagrange elements, piecewise linear Brezzi-Douglas-Marini (BDM)
elements [14], and piecewise constants on the criss-cross grids.

Similar to the continuous level, to get the BGG complex, we start from
$W^{0}$, and connect the two rows to obtain a second order operator
$\operatorname{curl}\operatorname{curl}$. Then we restrict to
$\mathcal{N}(\operatorname{sskw})\subset Y^{1}$. The degrees of freedom of
$\mathcal{N}(\operatorname{sskw})$ can be obtained by eliminating the degrees
of freedom of $W^{2}$ from $Y^{1}$. More precisely, all the degrees of freedom
on the edges are retained, while we get $5=16-11$ interior degrees of freedom.
The finite element stress complex reads

(3.7)
${0}$${\mathcal{P}_{1}}$${W^{0}}$${\mathcal{N}(\operatorname{sskw})}$${Y^{2}}$${0.}$$\scriptstyle{\subset}$$\scriptstyle{\operatorname{curl}\operatorname{curl}}$$\scriptstyle{\operatorname{div}}$

Standard argument with a dimension count implies the following.

######  Theorem 2.

The sequence (3.7) is a complex. On contractible domains, (3.7) is exact
except for index zero, where
$\mathcal{N}(\operatorname{curl}\operatorname{curl})=\mathcal{P}_{1}$, the
space of linear polynomials.

Johnson and Mercier [37] constructed a stable finite element elasticity pair
on the criss-cross grid on quadrilaterals. In [37], each component of the
displacement (corresponding to $Y^{2}$) consists of linear polynomials on
quadrilaterals (rather than the refined triangles). Therefore the local shape
function space on each quadrilateral has dimension 6. The stress space has the
same edge degrees of freedom as $\mathcal{N}(\operatorname{sskw})$ in (3.7),
but has 3 interior degrees of freedom, instead of 5 for
$\mathcal{N}(\operatorname{sskw})$. The pair in [37] thus does not fit in a
complex as derivatives of piecewise functions on a quadrilateral are not a
polynomial in general.

###  3.3. A new finite element strain complex

In this section, we construct finite elements fitting in (3.6), and derive the
resulting discrete BGG complex for (3.4). In particular, we discuss a discrete
reduction of the skew-symmetric components of the strain tensor. The finite
elements for (3.6) are summarized in Figure 5.

Figure 5. Finite element diagram (3.6) for deriving a discrete strain complex.

All these finite elements will be defined on a Clough-Tocher split of a
triangle, which refers to a split where we choose an interior point in the
triangle and connect it to the three vertices. Therefore a triangle is split
into three smaller ones (c.f. [38]). A Clough-Tocher mesh consists of
triangles with Clough-Tocher splits.

In the discussions below, we will use $\mathcal{P}_{r}(S)$ to denote the space
of polynomials of degree less than or equal to $r$ on a domain $S$. Typically,
$S$ is a triangle or its edges. We use $\mathcal{T}_{h}$ to denote a Clough-
Tocher mesh, and $T\in\mathcal{T}_{h}$ is a macroelement consisting of three
triangles. We use $T_{i}\in T_{\mathrm{CT}},i=1,2,3,$ to denote the three sub-
triangles of $T$. Moreover, $\Delta_{k}(\mathcal{T}_{h})$, or simply
$\Delta_{k}$, is the set of all $k$-dimensional subsimplicies of
$\mathcal{T}_{h}$. If $e\in\Delta_{1}$ is an edge of $\mathcal{T}_{h}$, then
$\tau_{e}$ and $n_{e}$ are the unit tangent and normal vectors, respectively.
We use $\Delta_{k}(T)$ to denote the set of $k$-simplices in $T$.

The complex
$(Z^{{\scriptscriptstyle\bullet}},d^{{\scriptscriptstyle\bullet}})$ on the top
row of (3.6) consists of a Clough-Tocher $H^{2}$ element, an
$H^{1}(\operatorname{rot})$ element, and a Lagrange element on the Clough-
Tocher mesh. The lowest order version of this complex (starting with piecewise
cubic polynomials) can be found in [21, Proposition 3].

– $Z^{0}$. Note that $Z^{0}$ is vector-valued. Each component of $Z^{0}$ is
the Clough-Tocher $C^{1}$ element with piecewise quartic polynomials. The
lowest order of such elements with piecewise cubic polynomials can be found
in, e.g., [38], while high order versions on a triangle can be found in [30].
The local shape functions on a macroelement $T$ are

$\mathcal{P}_{Z_{0}}(T):=\\{u_{h}\in
C^{1}(T):u_{h}|_{T_{i}}\in\mathcal{P}_{4}(T_{i})\otimes\mathbb{V},\,T_{i}\subset
T_{\mathrm{CT}},~{}i=1,2,3\\}.$

The degrees of freedom for $w\in Z^{0}$ can be given by

  * •

$w(x),\partial_{l}w(x),\quad x\in\Delta_{0}(T),~{}l=1,2$,

  * •

$\int_{e}\partial_{n}wq\,ds,\quad
q\in\mathcal{P}_{1}(e),~{}~{}e\in\Delta_{1}(T)$,

  * •

3 interior degrees of freedom.

Here $\partial_{n}w:=n\cdot\operatorname{grad}w$.

– $Z^{1}$. The space $Z^{1}$ consists of two copies of an
$H^{1}(\operatorname{rot})$-conforming element (continuous fields with
continuous $\operatorname{rot}$). The local shape function space for a single
copy (each column of the matrix) is

$\mathcal{P}_{Z_{1}}(T):=\\{u_{h}\in C^{0}(T):\operatorname{rot}u_{h}\in
C^{0}(T),u_{h}|_{T_{i}}\in\mathcal{P}_{3}(T_{i})\otimes\mathbb{M},\,T_{i}\subset
T_{\mathrm{CT}},~{}i=1,2,3\\}.$

The degrees of freedom for $\mathcal{P}_{Z_{1}}(T)$ can be given by the
following:

(3.8a) $\displaystyle u(a),\quad\forall a\in\Delta_{0}(T),$ (3.8b)
$\displaystyle\operatorname{rot}u(a),\quad\forall a\in\Delta_{0}(T),$ (3.8c)
$\displaystyle\int_{e}uq,\quad\forall q\in\mathcal{P}_{1}(e),\quad\forall
e\in\Delta_{1}(T),$ (3.8d)
$\displaystyle\int_{e}\operatorname{rot}u\cdot\tau_{e},\quad\forall
e\in\Delta_{1}(T),$ (3.8e) $\displaystyle\int_{e}\operatorname{rot}u\cdot
n_{e},\quad\forall e\in\Delta_{1}(T),$ (3.8f)
$\displaystyle\int_{T}u:q,\quad\forall
q\in\mathcal{C}(T),\,T\in\Delta_{2}(T),$

where $\mathcal{C}(T)$ denotes the bubble space on $T$.

The concrete form of the interior degrees of freedom is less important in the
sense that they do not affect the interelement continuity. The dimension of
$Z^{i},i=0,1,2$ and the bubble spaces (the dimension of interior degrees of
freedom) can be found in [30, Section 3].

– $Z^{2}$. The space $Z^{2}$ consists of two copies of the scalar quadratic
Lagrange elements on the Clough-Tocher split.

######  Theorem 3.

The degrees of freedom in (3.8) are unisolvent, and the complex
$(Z^{{\scriptscriptstyle\bullet}},d^{{\scriptscriptstyle\bullet}})$ consisting
of global finite element spaces is exact on contractible domains (except for
index 0).

###### Proof.

The unisolvency of $Z^{i}$ is a direct consequence of the dimension formulas
in [30, Section 3]. In fact, we can show that if all the vertex and edge
degrees of freedom vanish, then the field consists of bubbles in the sense of
[30]. Then vanishing integral against bubbles implies that the field actually
vanishes.

To show the exactness, we first observe that for any $u\in Z^{1}$, $u$ is also
in the Nédélec space. Then $\operatorname{rot}u=0$ implies that
$u=\operatorname{grad}\phi$ for some $\phi$ in the Lagrange space (exactness
of standard finite element de Rham sequence). As $u$ is $C^{0}$, we conclude
that $\phi$ is $C^{1}$, and therefore $\phi\in Z^{0}$. The exactness then
follows from a dimension count $1=\dim Z^{0}-\dim Z^{1}+\dim Z^{2}$ with the
Euler number $1=V-E+T$, where $V$, $E$ and $T$ are the numbers of vertices,
edges and triangles, respectively. ∎

The cohomology of the bottom row
$(V^{{\scriptscriptstyle\bullet}},d^{{\scriptscriptstyle\bullet}})$ can be
found in [21].

The bottom complex in (3.6) is the Clough-Tocher finite element Stokes complex
[21, (59)]. The first space is a scalar Clough-Tocher $C^{1}$ element, the
second is a vector-valued quadratic Lagrange element, and the last consists of
piecewise linear polynomials.

In the rest of this section, we start from (3.6) and derive a discrete version
of the strain complex

(3.9)
${0}$${Z^{0}}$${U^{1}}$${V^{2}}$${0.}$$\scriptstyle{\operatorname{def}_{h}}$$\scriptstyle{\operatorname{rot}\operatorname{rot}}$

The first and the last spaces in (3.14) are directly read from (3.6). The
space $U^{1}$ is $Z^{1}$ with reduced symmetry (roughly speaking, obtained by
eliminating the $V^{0}$ component from $Z^{1}$). The key difference compared
to the reduction of skew-symmetric components in existing constructions for
the stress complex [7, 20] is that $U^{1}$ is symmetric only in a discrete
sense, and $\operatorname{def}_{h}$ is different from the $\operatorname{def}$
operator on the continuous level which leads to symmetric fields. See (3.12)
and (3.13) below for precise definitions.

We first observe the following identity.

######  Lemma 1.

We have
$\operatorname{rot}\operatorname{skw}u=-\operatorname{grad}\operatorname{sskw}u$,
and therefore
$\operatorname{rot}u=-\operatorname{grad}\operatorname{sskw}u+\operatorname{rot}\operatorname{sym}u$.

Note that $\operatorname{sskw}(Z^{1})$ is not in $V^{0}$. Inspired by the
identities in Lemma 1, we define a discrete $\operatorname{sskw}$ operator,
$\operatorname{sskw}_{h}:Z^{1}\to V^{0}$, by setting the degrees of freedom of
$V^{0}$:

(3.10a)
$\displaystyle\operatorname{sskw}_{h}u(a):=\operatorname{sskw}u(a),\quad\forall
a\in\Delta_{0},$ (3.10b)
$\displaystyle\operatorname{grad}\operatorname{sskw}_{h}u(a):=-\operatorname{rot}u(a),\quad\forall
a\in\Delta_{0},$ (3.10c)
$\displaystyle\int_{e}\partial_{n}\operatorname{sskw}_{h}u:=-\int_{e}n\cdot\operatorname{rot}u,\quad\forall
e\in\Delta_{1}.$

We have $\operatorname{sskw}_{h}\operatorname{mskw}=I$, since

(3.11a)
$\displaystyle\operatorname{sskw}_{h}\operatorname{mskw}v(a)=\operatorname{sskw}\operatorname{mskw}v(a)=v(a),\quad\forall
a\in\Delta_{0},$ (3.11b)
$\displaystyle\operatorname{grad}\operatorname{sskw}_{h}\operatorname{mskw}v(a)=-\operatorname{rot}\operatorname{mskw}v(a)=\operatorname{grad}v(a),\quad\forall
a\in\Delta_{0},$ (3.11c)
$\displaystyle\int_{e}\partial_{n}\operatorname{sskw}_{h}\operatorname{mskw}v=-\int_{e}n\cdot\operatorname{rot}\operatorname{mskw}v=\int_{e}\partial_{n}v,\quad\forall
e\in\Delta_{1},$

i.e., for any $v\in V^{0}$, the degrees of freedom have the same value on
$\operatorname{sskw}_{h}\operatorname{mskw}v$ and $v$.

Define the shape function space of $U^{1}$ as the kernel of
$\operatorname{sskw}_{h}$ (in $Z^{1}$)

(3.12) $U^{1}:=\mathcal{N}(\operatorname{sskw}_{h},Z^{1}),$

and

(3.13)
$\operatorname{def}_{h}:=(I-\operatorname{mskw}\operatorname{sskw}_{h})\operatorname{grad}.$

Then we get a sequence

(3.14)
${0}$${Z^{0}}$${U^{1}}$${V^{2}}$${0.}$$\scriptstyle{\operatorname{def}_{h}}$$\scriptstyle{\operatorname{rot}\operatorname{rot}}$

We note that $\operatorname{def}_{h}Z^{0}\subseteq U^{1}$ since
$\operatorname{sskw}_{h}\operatorname{def}_{h}=\operatorname{sskw}_{h}(I-\operatorname{mskw}\operatorname{sskw}_{h})\operatorname{grad}=0$.
Moreover,
$\operatorname{rot}\operatorname{rot}\operatorname{def}_{h}=\operatorname{rot}\operatorname{rot}(I-\operatorname{mskw}\operatorname{sskw}_{h})\operatorname{grad}=0$,
since $\operatorname{rot}\operatorname{rot}\operatorname{mskw}=0$ by a diagram
chase on (3.6). So (3.14) is a complex.

Next, we show the exactness of (3.14) (except for index zero). First, we
observe that $\operatorname{rot}\operatorname{rot}:Z^{1}\to V^{2}$ is
surjective, as $\operatorname{rot}:Z^{1}\to Z^{2}$ and
$\operatorname{rot}:V^{1}\to V^{2}$ are both so, and $Z^{2}=V^{1}$. The
following theorem states the exactness at $U^{1}$.

######  Theorem 4.

$\mathcal{N}(\operatorname{rot}\operatorname{rot},U^{1})=\operatorname{def}_{h}Z^{0}$.

###### Proof.

From the exactness of the following complex (see [19, Proposition 2.3])

(3.15)
${Z^{0}}$${Z^{1}}$${V^{2}}$${0}$${V^{0}}$$\scriptstyle{\operatorname{grad}}$$\scriptstyle{\operatorname{rot}\operatorname{rot}}$$\scriptstyle{\operatorname{mskw}}$

we have

$\mathcal{N}(\operatorname{rot}\operatorname{rot},Z^{1})=\operatorname{grad}Z^{0}+\operatorname{mskw}V^{0}.$

This implies that for any $u\in U^{1}$,
$\operatorname{rot}\operatorname{rot}u=0$, there exists $\phi\in Z^{0}$,
$\psi\in V^{0}$, such that

(3.16) $u=\operatorname{grad}\phi+\operatorname{mskw}\psi.$

Therefore

$0=\operatorname{sskw}_{h}u=\operatorname{sskw}_{h}\operatorname{grad}\phi+\operatorname{sskw}_{h}\operatorname{mskw}\psi=\operatorname{sskw}_{h}\operatorname{grad}\phi+\psi.$

Consequently, $\psi=-\operatorname{sskw}_{h}\operatorname{grad}\phi$.
Substituting this identity into (3.16),

$u=\operatorname{grad}\phi-\operatorname{mskw}\operatorname{sskw}_{h}\operatorname{grad}\phi=(I-\operatorname{mskw}\operatorname{sskw}_{h})\operatorname{grad}\phi=\operatorname{def}_{h}\phi.$

∎

Next, we define the degrees of freedom of $U^{1}$ by eliminating degrees of
freedom from $Z^{1}$ that correspond to the skew-symmetric components
(eliminating $V^{0}$ from $Z^{1}$).

The degrees of freedom of $U^{1}$ can be given as follows:

(3.17a) $\displaystyle\operatorname{sym}u(a),\quad\forall a\in\Delta_{0}(T),$
(3.17b) $\displaystyle\int_{e}uq,\quad\forall
q\in\mathcal{P}_{1}(e),\quad\forall e\in\Delta_{1}(T),$ (3.17c)
$\displaystyle\int_{e}\operatorname{rot}u\cdot\tau_{e},\quad\forall
e\in\Delta_{1}(T),$ (3.17d) $\displaystyle\int_{T}u:q,\quad\forall
q\in\mathcal{C}(T),\,T\in\Delta_{2}(T).$

######  Theorem 5.

The degrees of freedom in (3.17) are unisolvent for $U^{1}$.

###### Proof.

First we show that if all the degrees of freedom in (3.17) vanish on $u\in
U^{1}$, then $u=0$. To see this, we will use (3.17) and the fact that
$u\in\mathcal{N}(\operatorname{sskw}_{h})$ to verify that all the degrees of
freedom in (3.8) (i.e., functionals in (3.10)) vanish on $u$. Then from the
unisolvence of (3.8) on $Z^{1}$, we obtain $u=0$ since $U^{1}\subset Z^{1}$.
To this end, we verify the following facts:

  * •

$u(a)=0,\quad\forall a\in\Delta_{0}(T)$, by (3.17a) and (3.10a)

  * •

$\operatorname{rot}u(a)=-\operatorname{grad}\operatorname{sskw}_{h}u(a)=0,\quad\forall
a\in\Delta_{0}(T)$, by (3.10b)

  * •

$\int_{e}uq=0,\quad\forall q\in\mathcal{P}_{1}(e),\quad\forall
e\in\Delta_{1}(T)$, by (3.17b)

  * •

$\int_{e}\operatorname{rot}u\cdot\tau_{e}=0,\quad\forall e\in\Delta_{1}(T)$,
by (3.17c)

  * •

$\int_{e}\operatorname{rot}u\cdot
n_{e}=-\int_{e}\partial_{n}\operatorname{sskw}_{h}u=0,\quad\forall
e\in\Delta_{1}(T)$, by (3.10c)

  * •

interior degrees of freedom vanish by definition. by (3.17d)

By the unisolvence of $Z^{1}$, we get $u=0$.

Then it remains to show that the number of degrees of freedom is equal to the
dimension of the local shape function spaces. To this end, we note that
$\operatorname{mskw}:V^{0}\to Z^{1}$ in (3.6) is injective. Therefore the
dimension of the local shape function space of $U^{1}$ satisfies $\dim
U^{1}=\dim Z^{1}-\dim V^{0}$. On the other hand, we observe that the number of
reduced degrees of freedom of (3.17) from (3.8) is also equal to $\dim V^{0}$.
Therefore $\dim U^{1}$ is equal to the number of the degrees of freedom in
(3.17).

∎

######  Remark 2.

From the above construction, $U^{1}$ has
$H^{1}(\operatorname{rot})$-conformity since it is a subspace of $Z^{1}$.
However, the shape functions of $U^{1}$ may not be symmetric.

######  Remark 3 (strong and weak symmetry).

In the context of the Hellinger-Reissner principle for linear elasticity,
strong symmetry refers to finite elements which contain symmetric matrix-
valued functions as local shape functions. As the construct of such finite
elements is difficult in some circumstances, another approach was proposed
where one uses a Lagrange multiplier to impose the symmetry conditions [5].
This is referred to as the approach of weak symmetry. The BGG diagram plays a
key role in choosing inf-sup stable finite element pairs for weak symmetry
[5]. In this case, the two rows in (3.3) do not match precisely in the sense
that $Y^{0}$ is not identical to $W^{1}$. Instead, one chooses $Y^{0}$ to be
larger than $W^{1}$ and uses an interpolation to replace $I$ in (3.3). The
weak symmetry approach is related to the strong symmetry in several ways.
First, if the finite elements in a discrete BGG diagram match precisely (e.g.,
$\operatorname{sskw}Y^{1}=W^{2}$), then implementing the formulation of
elasticity with weakly imposed symmetry leads to stress tensors with strong
symmetry. This was discovered in [31] referred to as “serendipitous exact
symmetry”. Second, in the case of “serendipitous exact symmetry”, on the
finite element level, one may further eliminate shape functions and degrees of
freedom which correspond to the skew-symmetric components to obtain finite
elements with strong symmetry [12, 20].

The finite element diagram chase in this paper fits in this picture, as we
eliminate shape functions and degrees of freedom on the finite element level.
Nevertheless, the strain complex calls for new techniques as the connecting
map is injective. We cannot repeat the argument in [12, 20] or [5] which
requires surjectivity. Extending formulations with weakly imposed symmetry to
the part of the BGG diagrams where the connecting maps are injective remains
open.

##  4\. Conclusions

In this paper, we presented a nonlinear version of the elasticity complex. The
spaces and operators encode key structures in elasticity theory and Riemannian
geometry. The exactness at index zero and one corresponds to a rigidity
theorem and a special case of equal-dimensional embedding of Riemannian
manifolds. The Korn inequality also has nonlinear versions [24, 25]. Such
results were established in the literature with various regularity
assumptions, c.f., [23, 40]. For the linearized versions, the BGG machinery
provides a systematic approach for establishing the corresponding results.
Therefore, as a future direction, it is of interest to investigate insights
that the BGG machinery and the linearized theory may bring to the nonlinear
case.

In the second part of the paper, we constructed finite element elasticity
complexes in 2D following a diagram chase. We focused on examples with the
lowest order polynomials. The main purpose was to demonstrate the diagram
chase with injective connecting maps. The results may be generalized to higher
order polynomials and higher dimensions (c.f., [19]).

It is also of interest to investigate the role of nonlinear complexes in
computation for nonlinear elasticity and geometric problems, e.g., [2, 13]. It
remains open to fit discrete (finite element) spaces in nonlinear complexes
and preserve the exactness (rigidity and embedding theorems at the discrete
level). A significant challenge is to clarify the algebraic varieties defined
by the operators in the nonlinear complex.

## Appendix: coordinate forms of linearized operators

In this appendix, we present the coordinate (vector/matrix) form of some
linearized operators in geometry. In particular, this clarifies the terms in
the elasticity complex. Similar formulas can be found in, e.g., [18, 39].

Let $g_{ij}$ be a Riemannian metric. The Christoffel symbols have the form

$\Gamma_{jk}^{i}=\frac{1}{2}g^{il}\left(\frac{\partial g_{lj}}{\partial
x^{k}}+\frac{\partial g_{lk}}{\partial x^{j}}-\frac{\partial g_{jk}}{\partial
x^{l}}\right),$

and their contracted versions

$\Gamma^{i}=\frac{1}{2}g^{jk}g^{il}\left(\frac{\partial g_{lj}}{\partial
x^{k}}+\frac{\partial g_{lk}}{\partial x^{j}}-\frac{\partial g_{jk}}{\partial
x^{l}}\right).$

The linearization around the Euclidean metric is thus

$\displaystyle\Gamma^{i}$
$\displaystyle\sim\frac{1}{2}\delta^{jk}\delta^{il}\left(\frac{\partial
g_{lj}}{\partial x^{k}}+\frac{\partial g_{lk}}{\partial x^{j}}-\frac{\partial
g_{jk}}{\partial x^{l}}\right)$ $\displaystyle=\frac{\partial g_{ij}}{\partial
x^{j}}-\frac{1}{2}\frac{\partial g_{j}^{j}}{\partial
x^{l}}=\operatorname{div}g-\frac{1}{2}\nabla\operatorname{tr}(g)=\operatorname{div}S^{-1}g.$

This shows that the linearized contracted Christoffel symbol corresponds to
the operator $\operatorname{div}S^{-1}$.

The linearized Ricci tensor around the Euclidean metric:

$R_{ij}=1/2(-\partial_{k}\partial^{k}g_{ij}-\partial_{i}\partial_{j}g_{k}^{k}+\partial_{i}\partial^{k}g_{kj}+\partial_{j}\partial^{k}g_{ik}),$

and in terms of vector/matrix notation,

$Ric=1/2(-\Delta
g-\nabla\nabla\operatorname{tr}(g)+2\operatorname{def}\operatorname{div}g).$

Taking the trace,

$R=-\Delta\operatorname{tr}(g)+\operatorname{div}\operatorname{div}g=-\operatorname{div}\operatorname{div}Sg.$

This corresponds to the scalar (Hamiltonian) constrain in the linearized
Einstein equations [39, (5.7)]. The linearization of the Einstein tensor
$G_{ij}:=R_{ij}-1/2Rg_{ij}$ becomes

$G_{ij}=1/2\left[-\partial_{k}\partial^{k}g_{ij}-\partial_{i}\partial_{j}g_{k}^{k}+\partial_{i}\partial^{k}g_{kj}+\partial_{j}\partial^{k}g_{ik}+(\partial_{k}\partial^{k}g_{l}^{l}-\partial^{k}\partial^{l}g_{kl})\eta_{ij}\right],$

where $\eta_{ij}$ is the Euclidean metric, and in terms of vector/matrix
notation,

(4.1) $\displaystyle G$ $\displaystyle=1/2[-\Delta
g-\operatorname{hess}\operatorname{tr}(g)+2\operatorname{def}\operatorname{div}g+(\Delta\operatorname{tr}(g)-\operatorname{div}\operatorname{div}g)I].$

In 2D, by a straightforward calculation, we have

(4.2)
$R=-\operatorname{div}\operatorname{div}Sg=\partial_{1}^{2}g_{22}+\partial_{2}^{2}g_{11}-2\partial_{1}\partial_{2}g_{12}=\operatorname{rot}\operatorname{rot}g.$

In 3D, we have

$G=\frac{1}{2}\operatorname{inc}g,$

and

$Ric=S^{-1}G=1/2S^{-1}\operatorname{inc}g.$

From (4.1), we have

$\operatorname{tr}(G)=\frac{1}{2}[\Delta\operatorname{tr}(g)-\operatorname{div}\operatorname{div}g].$

We have used the space dimension $n=3$ in the above calculation. Then we have

(4.3) $Ric=S^{-1}G=1/2[-\Delta
g-\operatorname{hess}\operatorname{tr}(g)+2\operatorname{def}\operatorname{div}g+1/2(\Delta\operatorname{tr}(g)-\operatorname{div}\operatorname{div}g)I]$

######  Lemma 2.

In three space dimensions, $R_{ijkl}=\frac{1}{2}G_{ijkl}$, where
$G_{ijkl}:=\epsilon_{ij}^{\quad s}\epsilon_{kl}^{\quad
t}(\operatorname{inc}g)_{st}$ is the fourth-order tensor version of
$\operatorname{inc}g$ (obtained by identifying each of the two indices of
$\operatorname{inc}g$ with a “skew-symmetric matrix”), and $R_{ijkl}$ is the
linearized Riemannian tensor of $g$. Moreover, the linearized Einstein tensor
$Ein=1/2\operatorname{inc}g$.

###### Proof.

$\displaystyle G_{pqcd}$ $\displaystyle=\epsilon_{pq}^{\quad
i}\epsilon_{cd}^{\quad
j}\epsilon_{iab}\epsilon_{jst}\partial^{a}\partial^{s}g^{bt}$
$\displaystyle=\left(\delta_{pa}\delta_{qb}-\delta_{pb}\delta_{qa}\right)\left(\delta_{cs}\delta_{dt}-\delta_{ct}\delta_{ds}\right)\partial^{a}\partial^{s}g^{bt}$
$\displaystyle=\left(\delta_{cs}\delta_{dt}-\delta_{ct}\delta_{ds}\right)\left(\partial_{p}\partial^{s}g_{q}^{t}-\partial_{q}\partial^{s}g_{p}^{t}\right)$
$\displaystyle=2\left(\delta_{cs}\delta_{dt}-\delta_{ct}\delta_{ds}\right)\left(\partial_{[p}\partial^{s}g_{q]}^{~{}~{}t}\right)$
$\displaystyle=4\partial_{[p}\partial_{[[c}g_{q]d]]},$

where $\partial_{[p}\partial_{[[c}g_{q]d]]}$ denotes the skew-symmetrization
with $p$, $q$ and $c$, $d$, respectively.

On the other hand, we have the linearized Riemamian tensor
$R_{ij~{}~{}~{}l}^{~{}~{}~{}\,\,k}=2\partial_{[i}\Gamma_{j]l}^{k}$, and

$\displaystyle\Gamma_{jl}^{k}$
$\displaystyle=\frac{1}{2}g^{ks}\left(\frac{\partial{g}_{sj}}{\partial{x}^{l}}+\frac{\partial{g}_{sl}}{\partial{x}^{j}}-\frac{\partial{g}_{jl}}{\partial{x}^{s}}\right)\sim\frac{1}{2}\left(\frac{\partial
g^{k}_{~{}~{}j}}{\partial x^{l}}+\frac{\partial g^{k}_{~{}~{}l}}{\partial
x^{j}}-\frac{\partial g_{jl}}{\partial x^{k}}\right)\sim
g_{j}^{~{}~{}[k}\nabla_{l]}+\frac{1}{2}\frac{\partial
g^{k}_{~{}~{}l}}{\partial x^{j}}.$

This implies that
$R_{ij~{}~{}~{}l}^{~{}~{}~{}\,\,k}=2\partial_{[i}\Gamma_{j]l}^{k}\sim
2\partial_{[i}g_{j]}^{~{}~{}[[k}\nabla_{l]]}$, i.e.,
$Riem=2\operatorname{mskw}\operatorname{inc}g$, where $\operatorname{mskw}$
maps each index of $\operatorname{inc}g$ to two indices corresponding to a
skew-symmetric 2-tensor.

Taking trace on the identity
$R_{pq,cd}=1/2\epsilon_{pq}^{~{}~{}~{}i}\epsilon^{~{}~{}~{}j}_{cd}\epsilon_{iab}\epsilon_{jst}\partial^{a}\partial^{s}g^{bt}$,
we get

$\displaystyle R_{qd}$ $\displaystyle=-\frac{1}{2}[\epsilon_{pq}^{\quad
i}\epsilon_{~{}~{}d}^{p~{}~{}j}\epsilon_{iab}\epsilon_{jst}\partial^{a}\partial^{s}g^{bt}]$
$\displaystyle=-\frac{1}{2}[\left(\delta_{qd}\delta^{ij}-\delta_{q}^{j}\delta_{d}^{i}\right)\epsilon_{iab}\epsilon_{jst}\partial^{a}\partial^{s}g^{bt}]$
$\displaystyle=-\frac{1}{2}[\epsilon_{iab}\epsilon^{i}_{~{}~{}st}\partial^{a}\partial^{s}g^{bt}\delta_{qd}-\epsilon_{dab}\epsilon_{qst}\partial^{a}\partial^{s}g^{bt}]$
$\displaystyle=-\frac{1}{2}[(\delta_{as}\delta_{bt}-\delta_{at}\delta_{bs})\partial^{a}\partial^{s}g^{bt}\delta_{qd}-\epsilon_{dab}\epsilon_{qst}\partial^{a}\partial^{s}g^{bt}]$
$\displaystyle=-\frac{1}{2}[(\Delta\operatorname{tr}g-\operatorname{div}\operatorname{div}g)I-\operatorname{inc}g]_{qd}.$

We note that

$\displaystyle\operatorname{tr}\operatorname{inc}g$
$\displaystyle=\epsilon^{d}_{~{}~{}ab}\epsilon_{dst}\partial^{a}\partial^{s}g^{bt}$
$\displaystyle=(\delta_{as}\delta_{bt}-\delta_{at}\delta_{bs})\partial^{a}\partial^{s}g^{bt}$
$\displaystyle=\Delta\operatorname{tr}g-\operatorname{div}\operatorname{div}g.$

This implies that

$2Ric=S\operatorname{inc}g,$

and

$Ein=S^{-1}Ric=\frac{1}{2}\operatorname{inc}g.$

∎

## Acknowledgement

The work was supported by a Hooke Research Fellowship and a Royal Society
University Research Fellowship (URF$\backslash$R1$\backslash$221398). The
author would like to thank Douglas N. Arnold for pointing out the reference
[45].

## References

  * [1] P. Alfeld, A trivariate Clough-Tocher scheme for tetrahedral data, Computer Aided Geometric Design, 1 (1984), pp. 169–181.
  * [2] A. Angoshtari, M. F. Shojaei, and A. Yavari, Compatible-strain mixed finite element methods for 2D compressible nonlinear elasticity, Computer Methods in Applied Mechanics and Engineering, 313 (2017), pp. 596–631.
  * [3] A. Angoshtari and A. Yavari, Hilbert complexes of nonlinear elasticity, Zeitschrift für angewandte Mathematik und Physik, 67 (2016), pp. 1–30.
  * [4] J. H. Argyris, I. Fried, and D. W. Scharpf, The TUBA family of plate elements for the matrix displacement method, The Aeronautical Journal, 72 (1968), pp. 701–709.
  * [5] D. Arnold, R. Falk, and R. Winther, Mixed finite element methods for linear elasticity with weakly imposed symmetry, Mathematics of Computation, 76 (2007), pp. 1699–1723.
  * [6] D. N. Arnold, Finite element exterior calculus, SIAM, 2018.
  * [7] D. N. Arnold, R. S. Falk, and R. Winther, Differential complexes and stability of finite element methods II: The elasticity complex, Compatible spatial discretizations, (2006), pp. 47–67.
  * [8] D. N. Arnold, R. S. Falk, and R. Winther, Finite element exterior calculus, homological techniques, and applications, Acta numerica, 15 (2006), p. 1.
  * [9] , Finite element exterior calculus: from Hodge theory to numerical stability, Bulletin of the American Mathematical Society, 47 (2010), pp. 281–354.
  * [10] D. N. Arnold and K. Hu, Complexes from complexes, Foundations of Computational Mathematics, (2021), pp. 1–36.
  * [11] D. N. Arnold and J. Qin, Quadratic velocity/linear pressure Stokes elements, Advances in computer methods for partial differential equations, 7 (1992), pp. 28–34.
  * [12] D. N. Arnold and R. Winther, Mixed finite elements for elasticity, Numerische Mathematik, 92 (2002), pp. 401–419.
  * [13] F. Auricchio, L. Beirão da Veiga, C. Lovadina, A. Reali, R. L. Taylor, and P. Wriggers, Approximation of incompressible large deformation elastic problems: some unresolved issues, Computational Mechanics, 52 (2013), pp. 1153–1167.
  * [14] F. Brezzi, J. Douglas Jr, and L. D. Marini, Two families of mixed finite elements for second order elliptic problems, Numerische Mathematik, 47 (1985), pp. 217–235.
  * [15] J. Brüning and M. Lesch, Hilbert complexes, Journal of Functional Analysis, 108 (1992), pp. 88–132.
  * [16] A. Čap and K. Hu, BGG sequences with weak regularity and applications, arXiv preprint arXiv:2203.01300, (2022).
  * [17] A. Čap, J. Slovák, and V. Souček, Bernstein-Gelfand-Gelfand sequences, Annals of Mathematics, 154 (2001), pp. 97–113.
  * [18] S. H. Christiansen, On the linearization of Regge calculus, Numerische Mathematik, 119 (2011), pp. 613–640.
  * [19] S. H. Christiansen, J. Gopalakrishnan, J. Guzmán, and K. Hu, A discrete elasticity complex on three-dimensional Alfeld splits, arXiv preprint arXiv:2009.07744, (2020).
  * [20] S. H. Christiansen, J. Hu, and K. Hu, Nodal finite element de Rham complexes, Numerische Mathematik, 139 (2018), pp. 411–446.
  * [21] S. H. Christiansen and K. Hu, Generalized finite element systems for smooth differential forms and Stokes’ problem, Numerische Mathematik, 140 (2018), pp. 327–371.
  * [22] , Finite element systems for vector bundles: elasticity and curvature, Foundations of Computational Mathematics, (2022), pp. 1–52.
  * [23] P. G. Ciarlet, Linear and nonlinear functional analysis with applications, vol. 130, Siam, 2013.
  * [24] P. G. Ciarlet, L. Gratie, and C. Mardare, A nonlinear Korn inequality on a surface, Journal de mathématiques pures et appliquées, 85 (2006), pp. 2–16.
  * [25] P. G. Ciarlet and C. Mardare, Nonlinear Korn inequalities, Journal de Mathématiques Pures et Appliquées, 104 (2015), pp. 1119–1134.
  * [26] R. W. Clough and J. L. Tocher, Finite element stiffness matrices for analysis of plates in bending, Proceedings of the Conference on Matrix Methods in Structural Mechanics, (1965), pp. 515–545.
  * [27] M. Eastwood, A complex from linear elasticity, in Proceedings of the 19th Winter School” Geometry and Physics”, Circolo Matematico di Palermo, 2000, pp. 23–29.
  * [28] M. G. Eastwood, Variations on the de Rham complex, Notices AMS, 46 (1999), pp. 1368–1376.
  * [29] R. S. Falk and M. Neilan, Stokes complexes and the construction of stable finite elements with pointwise mass conservation, SIAM Journal on Numerical Analysis, 51 (2013), pp. 1308–1326.
  * [30] G. Fu, J. Guzmán, and M. Neilan, Exact smooth piecewise polynomial sequences on Alfeld splits, Mathematics of Computation, 89 (2020), pp. 1059–1091.
  * [31] J. Gopalakrishnan and J. Guzmán, A second elasticity element using the matrix bubble, IMA Journal of Numerical Analysis, 32 (2012), pp. 352–372.
  * [32] J. Guzmán, A. Lischke, and M. Neilan, Exact sequences on Powell–Sabin splits, Calcolo, 57 (2020), pp. 1–25.
  * [33] J. Guzmán, A. Lischke, and M. Neilan, Exact sequences on Worsey–Farin splits, Mathematics of Computation, (2022).
  * [34] J. Hu and S. Zhang, A family of conforming mixed finite elements for linear elasticity on triangular grids, arXiv preprint arXiv:1406.7457, (2014).
  * [35] K. Hu, Oberwolfach report: Discretization of Hilbert complexes, arXiv preprint arXiv:2208.03420, (2022).
  * [36] K. Hu, J. Sun, and Q. Zhang, Quadratic and cubic Lagrange finite elements for mixed formulation of Laplace eigenvalue problem on criss-cross meshes, in preparation, (2023).
  * [37] C. Johnson and B. Mercier, Some equilibrium finite element methods for two-dimensional elasticity problems, Numerische Mathematik, 30 (1978), pp. 103–116.
  * [38] M.-J. Lai and L. L. Schumaker, Spline functions on triangulations, no. 110, Cambridge University Press, 2007.
  * [39] L. Li, Regge finite elements with applications in solid mechanics and relativity, PhD thesis, University of Minnesota, 2018.
  * [40] S. Mardare, On isometric immersions of a Riemannian space with little regularity, Analysis and applications, 2 (2004), pp. 193–226.
  * [41] J. E. Marsden and T. J. Hughes, Mathematical foundations of elasticity, Courier Corporation, 1994.
  * [42] M. Neilan, Discrete and conforming smooth de Rham complexes in three dimensions, Mathematics of Computation, (2015).
  * [43] , The Stokes complex: A review of exactly divergence-free finite element pairs for incompressible flows, in 75 Years of Mathematics of Computation: Symposium on Celebrating 75 Years of Mathematics of Computation, November 1-3, 2018, the Institute for Computational and Experimental Research in Mathematics (ICERM), vol. 754, American Mathematical Soc., 2020, p. 141.
  * [44] D. Pauly and W. Zulehner, The elasticity complex: compact embeddings and regular decompositions, arXiv preprint arXiv:2001.11007, (2020).
  * [45] A. D. Rendall, Insufficiency of the Ricci and Bianchi identities for characterising curvature, Journal of Geometry and Physics, 6 (1989), pp. 159–163.
  * [46] A. Worsey and G. Farin, An n-dimensional Clough-Tocher interpolant, Constructive Approximation, 3 (1987), pp. 99–110.
  * [47] A. Ženíšek, Polynomial approximation on tetrahedrons in the finite element method, Journal of Approximation Theory, 7 (1973), pp. 334–351.
  * [48] S. Zhang, A new family of stable mixed finite elements for the 3D Stokes equations, Mathematics of computation, 74 (2005), pp. 543–554.
  * [49] , On the P1 Powell-Sabin divergence-free finite element for the Stokes equations, Journal of Computational Mathematics, (2008), pp. 456–470.
  * [50] , Quadratic divergence-free finite elements on Powell–Sabin tetrahedral grids, Calcolo, 48 (2011), pp. 211–244.