Mixed finite element methods for linear elasticity with weakly imposed symmetry

Abstract

In this paper, we construct new finite element methods for the approximation of the equations of linear elasticity in three space dimensions that produce direct approximations to both stresses and displacements. The methods are based on a modified form of the Hellinger–Reissner variational principle that only weakly imposes the symmetry condition on the stresses. Although this approach has been previously used by a number of authors, a key new ingredient here is a constructive derivation of the elasticity complex starting from the de Rham complex. By mimicking this construction in the discrete case, we derive new mixed finite elements for elasticity in a systematic manner from known discretizations of the de Rham complex. These elements appear to be simpler than the ones previously derived. For example, we construct stable discretizations which use only piecewise linear elements to approximate the stress field and piecewise constant functions to approximate the displacement field.

1. Introduction

The equations of linear elasticity can be written as a system of equations of the form

$$\begin{equation} A\sigma = \operatorname \epsilon u, \qquad \operatorname {div}\sigma = f \quad \text{in }\Omega . {\tag{1.1}} \end{equation}$$

Here the unknowns $\sigma$ and $u$ denote the stress and displacement fields engendered by a body force $f$ acting on a linearly elastic body which occupies a region $\Omega \subset \mathbb{R}^3$. Then $\sigma$ takes values in the space ${\mathbb{S}}:=\mathbb{R}^{3\times 3}_{\text{sym}}$ of symmetric matrices and $u$ takes values in ${\mathbb{V}}:=\mathbb{R}^3$. The differential operator $\operatorname \epsilon$ is the symmetric part of the gradient, the div operator is applied row-wise to a matrix, and the compliance tensor $A=A(x):{\mathbb{S}}\to {\mathbb{S}}$ is a bounded and symmetric, uniformly positive definite operator reflecting the properties of the body. If the body is clamped on the boundary $\partial \Omega$ of $\Omega$, then the proper boundary condition for the system 1.1 is $u = 0$ on $\partial \Omega$. For simplicity, this boundary condition will be assumed here. The issues that arise when other boundary conditions are assumed (e.g., the case of pure traction boundary conditions $\sigma n = g$) are discussed in 9.

The pair $(\sigma ,u)$ can alternatively be characterized as the unique critical point of the Hellinger–Reissner functional

$$\begin{equation} {\mathcal{J}}(\tau ,v)=\int _\Omega \bigl (\frac{1}{2} A\tau :\tau +\operatorname {div}\tau \cdot v-f\cdot v\bigr )\,dx. {\tag{1.2}} \end{equation}$$

The critical point is sought among all $\tau \in H(\operatorname {div},\Omega ;{\mathbb{S}})$, the space of square-integrable symmetric matrix fields with square-integrable divergence, and all $v\in L^2(\Omega ;{\mathbb{V}})$, the space of square-integrable vector fields. Equivalently, $(\sigma ,u)\in H(\operatorname {div},\Omega ;{\mathbb{S}})\times L^2(\Omega ;{\mathbb{V}})$ is the unique solution to the following weak formulation of the system 1.1:

$$\begin{equation} \begin{array}{lrll} \int _\Omega (A\sigma :\tau + \operatorname {div}\tau \cdot u)\,dx &=& 0, \quad &\tau \in H(\operatorname {div},\Omega ;{\mathbb{S}}),\\[5.69046pt] \int _\Omega \operatorname {div}\sigma \cdot v \,dx &=& \int _\Omega f \cdot v\,dx, \quad &v \in L^2(\Omega ;{\mathbb{V}}). \end{array} {\tag{1.3}} \end{equation}$$

A mixed finite element method determines an approximate stress field $\sigma _h$ and an approximate displacement field $u_h$ as the critical point of ${\mathcal{J}}$ over $\Sigma _h\times V_h$ where $\Sigma _h\subset H(\operatorname {div},\Omega ;{\mathbb{S}})$ and $V_h\subset L^2(\Omega ;{\mathbb{V}})$ are suitable piecewise polynomial subspaces. Equivalently, the pair $(\sigma _h,u_h)\in \Sigma _h\times V_h$ is determined by the weak formulation 1.3, with the test space restricted to $\Sigma _h\times V_h$. As is well known, the subspaces $\Sigma _h$ and $V_h$ cannot be chosen arbitrarily. To ensure that a unique critical point exists and that it provides a good approximation of the true solution, they must satisfy the stability conditions from Brezzi’s theory of mixed methods 1213.

Despite four decades of effort, no stable simple mixed finite element spaces for elasticity have been constructed. For the corresponding problem in two space dimensions, stable finite elements were presented in 10. For the lowest order element, the space $\Sigma _h$ is composed of piecewise cubic functions, with 24 degrees of freedom per triangle, while the space $V_h$ consists of piecewise linear functions. Another approach which has been discussed in the two-dimensional case is the use of composite elements, in which $V_h$ consists of piecewise polynomials with respect to one triangulation of the domain, while $\Sigma _h$ consists of piecewise polynomials with respect to a different, more refined, triangulation 5212331. In three dimensions, a partial analogue of the element in 10 has been proposed and shown to be stable in 1. This element uses piecewise quartic stresses with 162 degrees of freedom per tetrahedron, and piecewise linear displacements.

Because of the lack of suitable mixed elasticity elements, several authors have resorted to the use of Lagrangian functionals which are modifications of the Hellinger–Reissner functional given above 24627282930, in which the symmetry of the stress tensor is enforced only weakly or abandoned altogether. In order to discuss these methods, we consider the compliance tensor $A(x)$ as a symmetric and positive definite operator mapping ${\mathbb{M}}$ into ${\mathbb{M}}$, where ${\mathbb{M}}$ is the space of $3 \times 3$ matrices. In the isotropic case, for example, the mapping $\sigma \mapsto A \sigma$ has the form

$$\begin{equation*} A \sigma = \frac{1}{2\mu }\bigl (\sigma - \frac{\lambda }{2\mu + 3 \lambda } \operatorname {tr}(\sigma )I\bigr ), \end{equation*}$$

where $\lambda (x), \mu (x)$ are positive scalar coefficients, the Lamé coefficients. A modification of the variational principle discussed above is obtained if we consider the extended Hellinger–Reissner functional

$$\begin{equation} {\mathcal{J}}_e(\tau ,v,q)= {\mathcal{J}}(\tau ,v) + \int _\Omega \tau : q\,dx {\tag{1.4}} \end{equation}$$

over the space $H(\operatorname {div},\Omega ;{\mathbb{M}})\times L^2(\Omega ;{\mathbb{V}})\times L^2(\Omega ;\mathbb{K})$, where $\mathbb{K}$ denotes the space of skew symmetric matrices. We note that the symmetry condition for the space of matrix fields is now enforced through the introduction of the Lagrange multiplier, $q$. A critical point $(\sigma ,u,p)$ of the functional ${\mathcal{J}}_e$ is characterized as the unique solution of the system

$$\begin{equation} \begin{array}{lrll} \int _\Omega (A\sigma :\tau + \operatorname {div}\tau \cdot u + \tau : p)\,dx &=& 0, \quad &\tau \in H(\operatorname {div},\Omega ;{\mathbb{M}}),\\[5.69046pt] \int _\Omega \operatorname {div}\sigma \cdot v \,dx &=& \int _\Omega f \cdot v\,dx, \quad &v \in L^2(\Omega ;{\mathbb{V}}),\\[5.69046pt] \int _\Omega \sigma :q \,dx &=& 0, \quad &q \in L^2(\Omega ;\mathbb{K}). \end{array} {\tag{1.5}} \end{equation}$$

Clearly, if $(\sigma ,u,p)$ is a solution of this system, then $\sigma$ is symmetric, i.e., $\sigma \in H(\operatorname {div},\Omega ;{\mathbb{S}})$, and therefore the pair $(\sigma ,u) \in H(\operatorname {div},\Omega ;{\mathbb{S}})\times L^2(\Omega ;{\mathbb{V}})$ solves the corresponding system 1.3. On the other hand, if $(u,p)$ solves 1.3, then $u\in H^1(\Omega ;{\mathbb{V}})$ and, if we set $p$ to the skew-symmetric part of $\operatorname {grad}u$, then $(\sigma ,u,p)$ solves 1.5. In this respect, the two systems 1.3 and 1.5 are equivalent. However, the extended system 1.5 leads to new possibilities for discretization. Assume that we choose finite element spaces $\Sigma _h\times V_h \times Q_h \subset H(\operatorname {div},\Omega ;{\mathbb{M}})\times L^2(\Omega ;{\mathbb{V}})\times L^2(\Omega ;\mathbb{K})$ and consider a discrete system corresponding to 1.5. If $(\sigma _h, u_h,p_h) \in \Sigma _h\times V_h \times Q_h$ is a discrete solution, then $\sigma _h$ will not necessarily inherit the symmetry property of $\sigma$. Instead, $\sigma _h$ will satisfy the weak symmetry condition

$$\begin{equation*} \int _\Omega \sigma _h :q \,dx = 0, \qquad \text{for all } q \in Q_h. \end{equation*}$$

Therefore, these solutions in general will not correspond to solutions of the discrete system obtained from 1.3.

Discretizations based on the system 1.5 will be referred to as mixed finite element methods with weakly imposed symmetry. Such discretizations were already introduced by Fraejis de Veubeke in 21 and further developed in 4. In particular, the so-called PEERS element proposed in 4 for the corresponding problem in two space dimensions used a combination of piecewise linear functions and cubic bubble functions, with respect to a triangulation of the domain, to approximate the stress $\sigma$, piecewise constants to approximate the displacements, and continuous piecewise linear functions to approximate the Lagrange multiplier $p$. Prior to the PEERS paper, Amara and Thomas 2 developed methods with weakly imposed symmetry using a dual hybrid approach. The lowest order method they discussed approximates the stresses with quadratic polynomials plus bubble functions and the multiplier by discontinuous constant or linear polynomials. The displacements are approximated on boundary edges by linear functions. Generalizations of the idea of weakly imposed symmetry to other triangular elements, rectangular elements, and three space dimensions were developed in 28, 29, 30 and 24. In 29, a family of elements is developed in both two and three dimensions. The lowest order element in the family uses quadratics plus the curls of quartic bubble functions in two dimensions or quintic bubble functions in three dimensions to approximate the stresses, discontinuous linears to approximate the displacements, and discontinuous quadratics to approximate the multiplier. In addition, a lower order method is introduced that approximates the stress by piecewise linear functions augmented by the curls of cubic bubble functions plus a cubic bubble times the gradient of local rigid motions. The multiplier is approximated by discontinuous piecewise linear functions and the displacement by local rigid motions. Morley 24 extends PEERS to a family of triangular elements, to rectangular elements, and to three dimensions. In addition, the multiplier is approximated by nonconforming rather than continuous piecewise polynomials.

There is a close connection between mixed finite elements for linear elasticity and discretization of an associated differential complex, the elasticity complex, which will be introduced in §3 below. In fact, the importance of this complex was already recognized in 10, where mixed methods for elasticity in two space dimensions were discussed. The new ingredient here is that we utilize a constructive derivation of the elasticity complex starting from the de Rham complex. This construction is described in Eastwood 18 and is based on the the Bernstein–Gelfand–Gelfand resolution; cf. 11 and also 14. By mimicking the construction in the discrete case, we are able to derive new mixed finite elements for elasticity in a systematic manner from known discretizations of the de Rham complex. As a result, we can construct new elements in both two and three space dimensions which are significantly simpler than those derived previously. For example, we will construct stable discretizations of the system 1.5 which only use piecewise linear and piecewise constant functions, as illustrated in the figure below. For simplicity, the entire discussion of the present paper will be given in the three-dimensional case. A detailed discussion in two space dimensions can be found in 8. Besides the methods discussed here, we note that by slightly generalizing the approach of this paper, one can also analyze some of the previously known methods mentioned above that are also based on the weak symmetry formulation (see 19 for details).

Figure 1.

Elements for the stress, displacement, and multiplier in the lowest order case in two dimensions and three dimensions.

An alternative approach to construct finite element methods for linear elasticity is to consider a pure displacement formulation. Since the coefficient $A$ in 1.1 is invertible, the stress $\sigma$ can be eliminated using the first equation in 1.1, the stress-strain relation. This leads to the second order equation

$$\begin{equation} \operatorname {div}A^{-1} \operatorname \epsilon u = f \quad \text{in } \Omega {\tag{1.6}} \end{equation}$$

for the displacement $u$. A weak solution of this equation can be characterized as the global minimizer of the energy functional

$$\begin{equation*} {\mathcal{E}}(u)= \int _\Omega \bigl (\frac{1}{2} A^{-1}\operatorname \epsilon u:\operatorname \epsilon u + f\cdot u\bigr )\,dx \end{equation*}$$

over the Sobolev space $H^1_0(\Omega ;{\mathbb{V}})$. Here $H^1_0(\Omega ;{\mathbb{V}})$ denotes the space of all square integrable vector fields on $\Omega$, with square integrable derivatives, and which vanish on the boundary $\partial \Omega$. A finite element approach based on this formulation, where we seek a minimum over a finite element subspace of $H^1_0(\Omega ;{\mathbb{V}})$ is standard and discussed in textbooks, (e.g., 16). However, for more general models, arising, for example, in viscoelasticity and plasticity (cf. 15), the stress–strain relation is not local and an elimination of the stress $\sigma$ is impossible. For such models, a pure displacement model is excluded, and a mixed approach seems to be an obvious alternative. The construction of stable mixed elements for linear elasticity is an important step in the construction of mixed methods for these more complicated models. Another advantage of the mixed approach is that we automatically obtain schemes which are uniformly stable in the incompressible limit, i.e., as the Lamé parameter $\lambda$ tends to infinity. Since this behavior of mixed methods is well known, we will not focus further on this property here. A more detailed discussion in this direction can, for example, be found in 5.

An outline of the paper is as follows. In §2, we describe the notation to be used, state our main result, and provide some preliminary discussion on the relation between stability of mixed finite element methods and discrete exact complexes. In §3, we present two complexes related to the two mixed formulations of elasticity given by 1.3 and 1.5. In §4, we introduce the framework of differential forms and show how the elasticity complex can be derived from the de Rham complex. In §5, we derive discrete analogues of the elasticity complex beginning from discrete analogues of the de Rham complex and identify the required properties of the discrete spaces necessary for this construction. This procedure is our basic design principle. In §6, we apply the construction of the preceding section to specific discrete analogues of the de Rham complex to obtain a family of discrete elasticity complexes. In §7 we use this family to construct stable finite element schemes for the approximation of the mixed formulation of the equations of elasticity with weakly imposed symmetry. Finally, in §8, we show how a slightly more complicated procedure leads to a simplified elasticity element.

2. Notation, statement of main results, and preliminaries

We begin with some basic notation and hypotheses. We continue to denote by ${\mathbb{V}}=\mathbb{R}^3$ the space of $3$-vectors, by ${\mathbb{M}}$ the space of $3\times 3$ real matrices, and by ${\mathbb{S}}$ and $\mathbb{K}$ the subspaces of symmetric and skew symmetric matrices, respectively. The operators $\operatorname {sym}:{\mathbb{M}}\to {\mathbb{S}}$ and $\operatorname {skw}:{\mathbb{M}}\to \mathbb{K}$ denote the symmetric and skew symmetric parts, respectively. Note that an element of the space $\mathbb{K}$ can be identified with its axial vector in ${\mathbb{V}}$ given by the map $\operatorname {\operatorname {vec}}: \mathbb{K}\to {\mathbb{V}}$:

$$\begin{equation*} \operatorname {\operatorname {vec}}\begin{pmatrix} 0 & -v_3 & v_2 \\ v_3 & 0 & -v_1 \\ -v_2 & v_1 & 0 \end{pmatrix} = \begin{pmatrix} v_1 \\ v_2 \\ v_3 \end{pmatrix}, \end{equation*}$$

i.e., $\operatorname {\operatorname {vec}}^{-1}(v)w= v\times w$ for any vectors $v$ and $w$.

We assume that $\Omega$ is a domain in $\mathbb{R}^3$ with boundary $\partial \Omega$. We shall use the standard function spaces, like the Lebesgue space $L^2(\Omega )$ and the Sobolev space $H^s(\Omega )$. For vector-valued functions, we include the range space in the notation following a semicolon, so $L^2(\Omega ;{\mathbb{X}})$ denotes the space of square integrable functions mapping $\Omega$ into a normed vector space ${\mathbb{X}}$. The space $H(\operatorname {div},\Omega ;{\mathbb{V}})$ denotes the subspace of (vector-valued) functions in $L^2(\Omega ;{\mathbb{V}})$ whose divergence belongs to $L^2(\Omega )$. Similarly, $H(\operatorname {div},\Omega ;{\mathbb{M}})$ denotes the subspace of (matrix-valued) functions in $L^2(\Omega ;{\mathbb{M}})$ whose divergence (by rows) belongs to $L^2(\Omega ;{\mathbb{V}})$.

Assuming that ${\mathbb{X}}$ is an inner product space, then $L^2(\Omega ;{\mathbb{X}})$ has a natural norm and inner product, which will be denoted by $\|\,\cdot \,\|$ and $(\,\cdot \,,\,\cdot \,)$, respectively. For a Sobolev space $H^s(\Omega ;{\mathbb{X}})$, we denote the norm by $\|\,\cdot \,\|_s$ and for $H(\operatorname {div}, \Omega ;{\mathbb{X}})$, the norm is denoted by $\|v \|_{\operatorname {div}} := (\|v\|^2 + \|\operatorname {div}v\|^2)^{1/2}$. The space ${\mathcal{P}}_k(\Omega )$ denotes the space of polynomial functions on $\Omega$ of total degree $\le k$. Usually we abbreviate this to just ${\mathcal{P}}_k$.

In this paper we shall consider mixed finite element approximations derived from 1.5. These schemes take the form:

Find $(\sigma _h, u_h, p_h) \in \Sigma _h \times V_h \times Q_h$ such that

$$\begin{equation} \begin{array}{lrll} \int _\Omega (A\sigma _h:\tau + \operatorname {div}\tau \cdot u_h + \tau : p_h)\,dx &=& 0, \quad &\tau \in \Sigma _h,\\[5.69046pt] \int _\Omega \operatorname {div}\sigma _h \cdot v \,dx &=& \int _\Omega f \cdot v\,dx, \quad &v \in V_h,\\[5.69046pt] \int _\Omega \sigma _h :q \,dx &=& 0, \quad &q \in Q_h, \end{array} {\tag{2.1}} \end{equation}$$

where now $\Sigma _h\subset H(\operatorname {div},\Omega ;{\mathbb{M}})$, $V_h\subset L^2(\Omega ;{\mathbb{V}})$ and $Q_h \in L^2(\Omega ; \mathbb{K})$.

Following the general theory of mixed finite element methods (cf. 1213) the stability of the saddle–point system 2.1 is ensured by the following conditions:

(A$1$)

$\lVert \tau \rVert ^2_{\operatorname {div}} \le c_1(A\tau , \tau )$ whenever $\tau \in \Sigma _h$ satisfies $(\operatorname {div}\tau , v) =0\quad \forall v\in V_h$, and $(\tau , q) =0 \quad \forall q\in Q_h$,

 (A$2$)

for all nonzero $(v,q) \in V_h \times Q_h$, there exists nonzero $\tau \in \Sigma _h$ with $(\operatorname {div}\tau , v) + (\tau , q) \ge c_2 \lVert \tau \rVert _{\operatorname {div}}(\lVert v\rVert + \lVert q\rVert )$,

where $c_1$ and $c_2$ are positive constants independent of $h$.

The main result of this paper, given in Theorem 7.1, is to construct a new family of stable finite element spaces $\Sigma _h$, $V_h$, $Q_h$ that satisfy the stability conditions (A$1$) and (A$2$). We shall show that for $r \ge 0$, the choices of the Nédélec second family of H(div) elements of degree $r+1$ for $\Sigma _h$ (cf. 26) and of discontinuous piecewise polynomials of degree $r$ for $V_h$ and $Q_h$ provide a stable finite element approximation. In contrast to the previous work described in the introduction, no stabilizing bubble functions are needed; nor is interelement continuity imposed on the multiplier. In §8 we also discuss a somewhat simpler lowest order element ($r=0$) in which the local stress space is a strict subspace of the full space of linear matrix fields.

Our approach to the construction of stable mixed elements for elasticity is motivated by the success in developing stable mixed elements for steady heat conduction (i.e., the Poisson problem) based on discretizations of the de Rham complex. We recall (see, e.g., 7) that there is a close connection between the construction of such elements and discretizations of the de Rham complex

$$\begin{equation} \begin{CD} \mathbb{R}@.\,\hookrightarrow \,@. C^\infty (\Omega ) @>\operatorname {grad}>> C^\infty (\Omega ;{\mathbb{V}})@>\operatorname {curl}>> C^\infty (\Omega ; {\mathbb{V}}) @>\operatorname {div}>> C^\infty (\Omega ) @>>> 0. \end{CD} {\tag{2.2}} \end{equation}$$

More specifically, a key to the construction and analysis of stable mixed elements is a commuting diagram of the form

$$\begin{equation} \begin{CD} \mathbb{R}@.\,\hookrightarrow \,@. C^\infty (\Omega ) @>\operatorname {grad}>> C^\infty (\Omega ;{\mathbb{V}})@>\operatorname {curl}>> C^\infty (\Omega ; {\mathbb{V}}) @>\operatorname {div}>> C^\infty (\Omega ) @>>> 0\\@. @. @VV\Pi ^1_hV @VV\Pi _h^cV @VV\Pi ^d_hV @VV\Pi ^0_hV\\\mathbb{R}@.\,\hookrightarrow \,@. W_h @>\operatorname {grad}>> U_h @>\operatorname {curl}>> V_h @>\operatorname {div}>> Q_h @>>> 0. \end{CD} {\tag{2.3}} \end{equation}$$

Here, the spaces $V_h\subset H(\operatorname {div})$ and $Q_h\subset L^2$ are the finite element spaces used to discretize the flux and temperature fields, respectively. The spaces $U_h\subset H(\operatorname {curl})$ and $W_h\subset H^1$ are additional finite element spaces, which can be found for all well-known stable element choices. The bottom row of the diagram is a discrete de Rham complex, which is exact when the de Rham complex is (i.e., when the domain is contractible). The vertical operators are projections determined by the natural degrees of freedom of the finite element spaces. As pointed out in 7, there are many such discretizations of the de Rham complex.

A diagram analogous to 2.3, but with the de Rham complex replaced by the elasticity complex defined just below, will be crucial to our construction of stable mixed elements for elasticity. Discretization of the elasticity complex also gives insight into the difficulties of constructing finite element approximations of the mixed formulation of elasticity with strongly imposed symmetry; cf. 8.

3. The elasticity complex

We now proceed to a description of two elasticity complexes, corresponding to strongly or weakly imposed symmetry of the stress tensor. In the case of strongly imposed symmetry, relevant to the mixed elasticity system 1.3, the characterization of the divergence-free symmetric matrix fields will be needed. In order to give such a characterization, define $\operatorname {curl}: C^\infty (\Omega ; {\mathbb{M}}) \to C^\infty (\Omega ; {\mathbb{M}})$ to be the differential operator defined by taking $\operatorname {curl}$ of each row of the matrix. Then define a second order differential operator $J : C^\infty (\Omega ;{\mathbb{S}}) \to C^\infty (\Omega ;{\mathbb{S}})$ by

$$\begin{equation} J \tau = \operatorname {curl}(\operatorname {curl}\tau )^T, \quad \tau \in C^\infty (\Omega ;{\mathbb{S}}). {\tag{3.1}} \end{equation}$$

It is easy to check that $\operatorname {div}\circ J =0$ and that $J \circ \operatorname \epsilon = 0$. In other words,

$$\begin{equation} \begin{CD} {\mathbb{T}}@.\,\hookrightarrow \,@. C^\infty ({\mathbb{V}}) @>\operatorname \epsilon >> C^\infty ({\mathbb{S}}) @>J>> C^\infty ({\mathbb{S}}) @>\operatorname {div}>> C^\infty ({\mathbb{V}}) @>>> 0 \end{CD} {\tag{3.2}} \end{equation}$$

is a complex. Here the dependence of the domain $\Omega$ is suppressed, i.e., $C^\infty ({\mathbb{S}})= C^\infty (\Omega ; {\mathbb{S}})$, and ${\mathbb{T}}= {\mathbb{T}}(\Omega )$ denotes the six-dimensional space of infinitesimal rigid motions on $\Omega$, i.e., functions of the form $x\mapsto a+Bx$ with $a\in {\mathbb{V}}$ and $B\in \mathbb{K}$. In fact, when $\Omega$ is contractible, then 3.2 is an exact sequence, a fact which will follow from the discussion below. The complex 3.2 will be referred to as the elasticity complex.

A natural approach to the construction of stable mixed finite elements for elasticity would be to extend the complex 3.2 to a complete commuting diagram of the form 2.3, where 3.2 is the top row and the bottom row is a discrete analogue. However, due to the pointwise symmetry requirement on the discrete stresses, this construction requires piecewise polynomials of high order. For the corresponding problem in two space dimensions, such a complex was proposed in 10 with a piecewise cubic stress space; cf. also 8. An analogous complex was derived in the three-dimensional case in 3. It uses a piecewise quartic space, with 162 degrees of freedom on each tetrahedron for the stresses.

We consider the formulation based on weakly imposed symmetry of the stress tensor, i.e., the mixed system 1.5. Then the relevant complex is, instead of 3.2,

$$\begin{equation} \begin{CD} {\mathbb{T}}' @.\,\hookrightarrow \,@. C^\infty ({\mathbb{V}}\times \mathbb{K}) @>(\operatorname {grad}, -I)>> C^\infty ({\mathbb{M}}) @>J>> C^\infty ({\mathbb{M}}) @>(\operatorname {div}, \operatorname {skw})^T >> C^\infty ({\mathbb{V}}\times \mathbb{K}) @>>> 0. \end{CD} {\tag{3.3}} \end{equation}$$

Here,

$$\begin{equation*} {\mathbb{T}}' = \{\,(v,\operatorname {grad}v)\,|\,v\in {\mathbb{T}}\,\}, \end{equation*}$$

and $J : C^\infty (\Omega ;{\mathbb{M}}) \to C^\infty (\Omega ;{\mathbb{M}})$ denotes the extension of the operator defined on $C^\infty (\Omega ;{\mathbb{S}})$ by 3.1 such that $J \equiv 0$ on $C^\infty (\Omega ;\mathbb{K})$. We remark that $J$ may be written

$$\begin{equation} J\tau = \operatorname {curl}\Xi ^{-1}\operatorname {curl}\tau , {\tag{3.4}} \end{equation}$$

where $\Xi :{\mathbb{M}}\to {\mathbb{M}}$ is the algebraic operator

$$\begin{equation} \Xi \mu = \mu ^T - \operatorname {tr}(\mu )\delta , \qquad \Xi ^{-1}\mu = \mu ^T - \frac{1}{2}\operatorname {tr}(\mu )\delta , {\tag{3.5}} \end{equation}$$

with $\delta$ the identity matrix. Indeed, if $\tau$ is symmetric, then $\operatorname {curl}\tau$ is trace free, and therefore the definition 3.4 reduces to 3.1 on $C^\infty (\Omega ;{\mathbb{S}})$. On the other hand, if $\tau$ is skew with axial vector $u$, then $\operatorname {curl}\tau = -\Xi \operatorname {grad}u$, and so $\operatorname {curl}\Xi ^{-1}\operatorname {curl}\tau =0$.

Observe that there is a close connection between 3.2 and 3.3. In fact, 3.2 can be derived from 3.3 by performing a projection step. To see this, consider the diagram

$$\begin{equation} \begin{CD} {\mathbb{T}}'@.\,\hookrightarrow \,@. C^\infty ({\mathbb{V}}\times \mathbb{K}) @>(\operatorname {grad}, -I)>> C^\infty ({\mathbb{M}}) @>J>> C^\infty ({\mathbb{M}}) @>(\operatorname {div}, \operatorname {skw})^T >> C^\infty ({\mathbb{V}}\times \mathbb{K}) @>>> 0 \\@. @. @VV\pi _0V @VV\pi _1V @VV\pi _2V @VV\pi _3V \\{\mathbb{T}}@.\,\hookrightarrow \,@. C^\infty ({\mathbb{V}}) @>\operatorname \epsilon >> C^\infty ({\mathbb{S}}) @>J>> C^\infty ({\mathbb{S}}) @>\operatorname {div}>> C^\infty ({\mathbb{V}}) @>>> 0, \end{CD} {\tag{3.6}} \end{equation}$$

where the projection operators $\pi _k$ are defined by

$$\begin{equation*} \pi _0(u,q) = u, \quad \pi _1(\sigma ) = \pi _2(\sigma ) = \operatorname {sym}(\sigma ), \quad \pi _3(u,q) = u - \operatorname {div}q. \end{equation*}$$

We may identify $C^\infty ({\mathbb{V}})$ with a subspace of $C^\infty ({\mathbb{V}}\times \mathbb{K})$, namely,

$$\begin{equation*} \{(u,q): \, u \in C^\infty ({\mathbb{V}}), q = \operatorname {skw}(\operatorname {grad}u)\}. \end{equation*}$$

Under this identification, ${\mathbb{T}}\subset C^\infty ({\mathbb{V}})$ corresponds to ${\mathbb{T}}'\subset C^\infty ({\mathbb{V}}\times \mathbb{K})$. We identify the $C^\infty ({\mathbb{V}})$ on the right with a different subspace of $C^\infty ({\mathbb{V}}\times \mathbb{K})$, namely,

$$\begin{equation*} \{(u,q): \, u \in C^\infty ({\mathbb{V}}), q = 0\}. \end{equation*}$$

With these identifications, the bottom row is a subcomplex of the top row, and the operators $\pi _k$ are all projections. Furthermore, the diagram commutes. It follows easily that the exactness of the upper row implies exactness of the bottom row.

In the next section, we shall discuss these complexes further. In particular, we show the elasticity complex with weakly imposed symmetry, i.e., 3.3 follows from the de Rham complex 2.2. Hence, as a consequence of the discussion above, both 3.2 and 3.3 will follow from 2.2.

4. From the de Rham to the elasticity complex

The main purpose of this section is to demonstrate the connection between the de Rham complex 2.2 and the elasticity complex 3.3. In particular, we show that whenever 2.2 is exact, 3.3 is exact. This section serves as an introduction to a corresponding construction of a discrete elasticity complex, to be given in the next section. In the following section, the discrete complex will be used to construct stable finite elements for the system 1.5.

The de Rham complex 2.2 is most clearly stated in terms of differential forms. Here we briefly recall the definitions and properties we will need. We use a completely coordinate-free approach. For a slightly more expanded discussion and the expressions in coordinates see, e.g., 7, §4. We let $\Lambda ^k$ denote the space of smooth differential $k$-forms on $\Omega$, i.e. $\Lambda ^k = \Lambda ^k(\Omega ) = C^\infty (\Omega ; \text{Alt}^k{\mathbb{V}})$, where $\text{Alt}^k{\mathbb{V}}$ denotes the vector space of alternating $k$-linear maps on ${\mathbb{V}}$. If $\omega \in \Lambda ^k$ we let $\omega _x \in \text{Alt}^k{\mathbb{V}}$ denote $\omega$ evaluated at $x$, i.e., we use subscripts to indicate the spatial dependence.

Using the inner product on $\text{Alt}^k{\mathbb{V}}$ inherited from the inner product on ${\mathbb{V}}$ (see equation (4.1) of 7, §4), we may also define the Hilbert space $L^2\Lambda ^k(\Omega )=L^2(\Omega ;\text{Alt}^k{\mathbb{V}})$ of square integrable differential forms with norm denoted by $\|\,\cdot \,\|$, and also the $m$th order Sobolev space $H^m\Lambda ^k(\Omega )=H^m(\Omega ;\text{Alt}^k{\mathbb{V}})$, consisting of square integrable $k$-forms for which the norm

$$\begin{equation*} \|\omega \|_m := \bigl ( \sum _{|\alpha |\le m} \|\partial ^\alpha \omega \|^2\bigr )^{1/2} \end{equation*}$$

is finite (where the sum is over multi-indices of degree at most $m$).

Thus, $0$-forms are scalar functions and $1$-forms are covector fields. We will not emphasize the distinction between vectors and covectors, since, given the inner product in ${\mathbb{V}}$, we may identify a $1$-form $\omega$ with the vector field $v$ for which $\omega (p)=v\cdot p$, $p\in {\mathbb{V}}$. In the three-dimensional case, we can identify a $2$-form $\omega$ with a vector field $v$ and a $3$-form $\mu$ with a scalar field $c$ by

$$\begin{equation*} \omega (p,q)=v\cdot p\times q, \quad \mu (p,q,r) = c(p \times q \cdot r), \quad p,q,r\in {\mathbb{V}}. \end{equation*}$$

The exterior derivative $d=d_k: \Lambda ^k \to \Lambda ^{k+1}$ is defined by

$$\begin{equation} \begin{aligned} &d\omega _x(v_1,\ldots ,v_{k+1})\\ &\qquad =\sum _{j=1}^{k+1} (-1)^{j+1} \partial _{v_j}\omega _x(v_1,\ldots ,\hat{v}_j,\ldots ,v_{k+1}),\quad \omega \in \Lambda ^k, v_1,\ldots ,v_{k+1}\in {\mathbb{V}}, \end{aligned} {\tag{4.1}} \end{equation}$$

where the hat is used to indicate a suppressed argument and $\partial _{v}$ denotes the directional derivative in the direction of the vector $v$. It is useful to define

$$\begin{equation*} H\Lambda ^k = \{\, \omega \in L^2(\Omega ;\text{Alt}^k{\mathbb{V}})\,|\, d\omega \in L^2(\Omega ;\text{Alt}^{k+1}{\mathbb{V}})\,\}, \end{equation*}$$

with norm given by $\|\omega \|_{H\Lambda }^2=\|\omega \|^2+\|d\omega \|^2$. Using the identifications given above, the $d_k$ correspond to $\operatorname {grad}$, $\operatorname {curl}$, and $\operatorname {div}$ for $k=0,1,2$, respectively, and the $H\Lambda ^k$ correspond to $H^1$, $H(\operatorname {curl})$, $H(\operatorname {div})$, and, for $k=3$, $L^2$.

The de Rham complex 2.2 can then be written

$$\begin{equation} \begin{CD} \mathbb{R}@.\,\hookrightarrow \,@. \Lambda ^0 @>d>> \Lambda ^1 @>d>> \Lambda ^2@>d>> \Lambda ^3@>>> 0. \end{CD} {\tag{4.2}} \end{equation}$$

It is a complex since $d\circ d=0$.

A differential $k$-form $\omega$ on $\Omega$ may be restricted to a differential $k$-form on any submanifold $M\subset \bar{\Omega }$; at each point of $M$ the restriction of $\omega$ is an alternating linear form on tangent vectors. Moreover, if $\dim M=k$, the integral $\int _M \omega$ is defined.

If ${\mathbb{X}}$ is a vector space, then $\Lambda ^k({\mathbb{X}}) = \Lambda ^k(\Omega ;{\mathbb{X}})$ refers to the $k$-forms with values in ${\mathbb{X}}$, i.e., $\Lambda ^k({\mathbb{X}}) = C^\infty (\Omega ; \text{Alt}^k({\mathbb{V}};{\mathbb{X}}))$, where $\text{Alt}^k({\mathbb{V}};{\mathbb{X}})$ are alternating $k$-linear forms on ${\mathbb{V}}$ with values in ${\mathbb{X}}$. Given an inner product on ${\mathbb{X}}$, we obtain an inner product on $\Lambda ^k({\mathbb{X}})$. Obviously the corresponding complex

$$\begin{equation} \begin{CD} {\mathbb{X}}@.\,\hookrightarrow \,@. \Lambda ^0({\mathbb{X}}) @>d>> \Lambda ^1 ({\mathbb{X}})@>d>> \Lambda ^2({\mathbb{X}}) @>d>> \Lambda ^3({\mathbb{X}})@>>> 0, \end{CD} {\tag{4.3}} \end{equation}$$

is exact whenever the de Rham complex is.

We now construct the elasticity complex as a subcomplex of a complex isomorphic to the de Rham complex with values in the six-dimensional vector space ${\mathbb{W}}:=\mathbb{K}\times {\mathbb{V}}$. First, for any $x\in \mathbb{R}^3$ we define $K_x:{\mathbb{V}}\to \mathbb{K}$ by $K_x v =2\operatorname {skw}(xv^T)$. We then define an operator $K:\Lambda ^k(\Omega ;{\mathbb{V}})\to \Lambda ^k(\Omega ;\mathbb{K})$ by

$$\begin{equation} (K\omega )_x(v_1,\ldots ,v_k) = K_x[\omega _x(v_1,\ldots ,v_k)]. {\tag{4.4}} \end{equation}$$

Next, we define an isomorphism $\Phi :\Lambda ^k({\mathbb{W}})\to \Lambda ^k({\mathbb{W}})$ by

$$\begin{equation*} \Phi (\omega ,\mu ) = (\omega + K\mu , \mu ), \end{equation*}$$

with inverse given by

$$\begin{equation*} \Phi ^{-1}(\omega ,\mu ) = (\omega - K\mu , \mu ). \end{equation*}$$

Next, define the operator ${\mathcal{A}}: \Lambda ^k({\mathbb{W}}) \to \Lambda ^{k+1}({\mathbb{W}})$ by ${\mathcal{A}}= \Phi d \Phi ^{-1}$. Inserting the isomorphisms $\Phi$ in the ${\mathbb{W}}$-valued de Rham sequence, we obtain a complex

$$\begin{equation} \begin{CD} \Phi ({\mathbb{W}}) @.\,\hookrightarrow \,@. \Lambda ^0({\mathbb{W}}) @>{\mathcal{A}}>> \Lambda ^1 ({\mathbb{W}})@>{\mathcal{A}}>> \Lambda ^2({\mathbb{W}}) @>{\mathcal{A}}>> \Lambda ^3({\mathbb{W}})@>>> 0, \end{CD} {\tag{4.5}} \end{equation}$$

which is exact whenever the de Rham complex is.

The operator ${\mathcal{A}}$ has a simple form. Using the definition of $\Phi$, we obtain for $(\omega ,\mu ) \in \Lambda ^k({\mathbb{W}})$,

$$\begin{equation*} {\mathcal{A}}(\omega ,\mu ) = \Phi \circ d (\omega - K\mu ,\mu ) = \Phi (d\omega - d K\mu ,d\mu ) = (d\omega - S \mu , d\mu ), \end{equation*}$$

where $S= S_k:\Lambda ^k({\mathbb{V}}) \to \Lambda ^{k+1}(\mathbb{K})$, $k=0,1,2$ is given by $S=dK - Kd$. Using the definition 4.1 of the exterior derivative, the definition 4.4 of $K$, and the Leibniz rule

$$\begin{equation} d(\omega \wedge \mu ) = d\omega \wedge \mu + (-1)^k \omega \wedge d\mu , \quad \omega \in \Lambda ^k,\nobreakspace \mu \in \Lambda ^\ell , {\tag{4.6}} \end{equation}$$

we obtain

$$\begin{equation} (S\omega )(v_1,\ldots ,v_{k+1}) = \sum _{j=1}^{k+1}(-1)^{j+1} K_{v_j}[\omega (v_1,\ldots ,\hat{v}_j,\ldots v_{k+1})],\quad \omega \in \Lambda ^k(\Omega ;{\mathbb{V}}). {\tag{4.7}} \end{equation}$$

Note that the operator $S$ is purely algebraic, and independent of $x$.

Since $d^2 = 0$, we have

$$\begin{equation*} dS = d^2 K - dKd = - (dK - Kd)d \end{equation*}$$

or

$$\begin{equation} dS = - Sd. {\tag{4.8}} \end{equation}$$

Noting that

$$\begin{multline*} (S_1\mu )(v_1,v_2)=K_{v_1}[\mu (v_2)]-K_{v_2}[\mu (v_1)]=2\operatorname {skw}[v_1\mu (v_2)^T-v_2\mu (v_1)^T], \\ \mu \in \Lambda ^1(\Omega ;{\mathbb{V}}),\nobreakspace v_1,v_2\in {\mathbb{V}}, \end{multline*}$$

we find, using the identity

$$\begin{equation} a\times b=-2\operatorname {\operatorname {vec}}\operatorname {skw}ab^T, {\tag{4.9}} \end{equation}$$

that $S_1$ is invertible with

$$\begin{multline*} (S_1^{-1}\omega )(v_1)\times v_2\cdot v_3=\frac{1}{2}[ \operatorname {\operatorname {vec}}\bigl (\omega (v_2,v_3)\bigr )\cdot v_1 - \operatorname {\operatorname {vec}}\bigl (\omega (v_1,v_2)\bigr )\cdot v_3 + \operatorname {\operatorname {vec}}\bigl (\omega (v_1,v_3)\bigr )\cdot v_2 ], \\ \omega \in \Lambda ^2(\Omega ;\mathbb{K}),v_1, v_2, v_3\in {\mathbb{V}}. \end{multline*}$$

We now define the desired subcomplex. Define

$$\begin{equation*} \Gamma ^1 = \{\,(\omega , \mu ) \in \Lambda ^1(\Omega ;{\mathbb{W}})\,|\, d\omega = S_1\mu \,\}, \quad \Gamma ^2 = \{\,(\omega , \mu ) \in \Lambda ^2(\Omega ;{\mathbb{W}})\,|\, \omega = 0\,\}, \end{equation*}$$

with projections $\pi ^1:\Lambda ^1(\Omega ;{\mathbb{W}})\to \Gamma ^1$ and $\pi ^2:\Lambda ^2(\Omega ;{\mathbb{W}})\to \Gamma ^2$ given by

$$\begin{equation*} \pi ^1(\omega , \mu ) = (\omega , S_1^{-1}d \omega ), \quad \pi ^2(\omega , \mu ) = (0 , \mu + d S_1^{-1} \omega ). \end{equation*}$$

Using 4.8, it is straightforward to check that ${\mathcal{A}}$ maps $\Lambda ^0({\mathbb{W}})$ into $\Gamma ^1$ and $\Gamma ^1$ into $\Gamma ^2$, and that the diagram

$$\begin{equation} \begin{CD} \Phi ({\mathbb{W}}) @.\,\hookrightarrow \,@. \Lambda ^0({\mathbb{W}}) @>{\mathcal{A}}>> \Lambda ^1 ({\mathbb{W}})@>{\mathcal{A}}>> \Lambda ^2({\mathbb{W}}) @>{\mathcal{A}}>> \Lambda ^3({\mathbb{W}})@>>> 0\\@. @. @VVidV @VV\pi ^1V @VV\pi ^2V @VVidV\\\Phi ({\mathbb{W}}) @.\,\hookrightarrow \,@. \Lambda ^0({\mathbb{W}}) @>{\mathcal{A}}>> \Gamma ^1@>{\mathcal{A}}>> \Gamma ^2@>{\mathcal{A}}>> \Lambda ^3({\mathbb{W}})@>>> 0 \end{CD} {\tag{4.10}} \end{equation}$$

commutes, and therefore the subcomplex in the bottom row is exact when the de Rham complex is. This subcomplex is, essentially, the elasticity complex. Indeed, by identifying elements $(\omega ,\mu ) \in \Gamma ^1$ with $\omega \in \Lambda ^1(\mathbb{K})$, and elements $(0,\mu ) \in \Gamma ^2$ with $\mu \in \Lambda ^2({\mathbb{V}})$, the subcomplex becomes

$$\begin{equation} \begin{aligned} \Phi ({\mathbb{W}}) \hookrightarrow \Lambda ^0(\mathbb{K}\times {\mathbb{V}}) &\xrightarrow {(d_0,-S_0)} \Lambda ^1 (\mathbb{K})\xrightarrow {d_1\circ S_1^{-1}\circ d_1} \Lambda ^2({\mathbb{V}}) \\ &\xrightarrow {(-S_2,d_2)^T} \Lambda ^3(\mathbb{K}\times {\mathbb{V}})\to 0. \end{aligned} {\tag{4.11}} \end{equation}$$

This complex may be identified with 3.3. As an initial step of this identification we observe that the algebraic operator $\Xi : C^\infty ({\mathbb{M}})\to C^\infty ({\mathbb{M}})$ appearing in 3.3 via 3.4 and the operator $S_1 : \Lambda ^1({\mathbb{V}}) \to \Lambda ^2(\mathbb{K})$ are connected by the identity

$$\begin{equation} \Xi = \Upsilon _2^{-1}S_1\Upsilon _1, {\tag{4.12}} \end{equation}$$

where $\Upsilon _1: C^\infty ({\mathbb{M}}) \to \Lambda ^1({\mathbb{V}})$ and $\Upsilon _2: C^\infty ({\mathbb{M}}) \to \Lambda ^2(\mathbb{K})$ are given by $\Upsilon _1F(v) =Fv$ and $\Upsilon _2 F(v_1,v_2) = \operatorname {\operatorname {vec}}^{-1}F(v_1 \times v_2)$ for $F \in C^\infty ({\mathbb{M}})$. In fact, using 4.9, we have for any $v_1, v_2 \in {\mathbb{V}}$,

$$\begin{align*} S_1 \Upsilon _1F(v_1,v_2) &= 2\operatorname {skw}[v_1(Fv_2)^T -v_2(Fv_1)^T]\\ &= \operatorname {\operatorname {vec}}^{-1}(v_2 \times Fv_1 - v_1 \times Fv_2). \end{align*}$$

On the other hand,

$$\begin{equation*} \Upsilon _2 \Xi F(v_1,v_2) = \operatorname {\operatorname {vec}}^{-1}[\Xi F(v_1 \times v_2)], \end{equation*}$$

and hence 4.12 follows from the algebraic identity

$$\begin{equation*} \Xi F (v_1 \times v_2) = v_2 \times Fv_1 - v_1 \times Fv_2, \end{equation*}$$

which holds for any $F \in {\mathbb{M}}$.

We may further identify the four spaces of fields in 3.3 with the corresponding spaces of forms in 4.11 in a natural way:

•

$(u,p)\in C^\infty ({\mathbb{V}}\times \mathbb{K})\sim (\operatorname {\operatorname {vec}}^{-1}u,\operatorname {\operatorname {vec}}p)\in \Lambda ^0(\mathbb{K}\times {\mathbb{V}})$.

•

$F\in C^\infty ({\mathbb{M}}) \sim \omega \in \Lambda ^1(\mathbb{K})$ given by $\omega (v)=\operatorname {\operatorname {vec}}^{-1}(Fv)$.

•

$F\in C^\infty ({\mathbb{M}}) \sim \mu \in \Lambda ^2({\mathbb{V}})$ given by $\mu (v_1,v_2)=F(v_1\times v_2)$.

•

$(u,p)\in C^\infty ({\mathbb{V}}\times \mathbb{K})\sim (\omega ,\mu )\in \Lambda ^3(\mathbb{K}\times {\mathbb{V}})$ given by $\omega (v_1,v_2,v_3)=p(v_1\times v_2\cdot v_3)$, $\mu (v_1,v_2,v_3)=u(v_1\times v_2\cdot v_3)$.

Under these identifications, we find that

•

$d_0:\Lambda ^0(\mathbb{K})\to \Lambda ^1(\mathbb{K})$ corresponds to the row-wise gradient $C^\infty ({\mathbb{V}})\to C^\infty ({\mathbb{M}})$.

•

$S_0:\Lambda ^0({\mathbb{V}})\to \Lambda ^1(\mathbb{K})$ corresponds to the inclusion of $C^\infty (\mathbb{K}) \to C^\infty ({\mathbb{M}})$.

•

$d_1\circ S_1^{-1}\circ d_1:\Lambda ^1(\mathbb{K})\to \Lambda ^2({\mathbb{V}})$ corresponds to $J=\operatorname {curl}\Xi ^{-1}\operatorname {curl}:C^\infty ({\mathbb{M}})\to C^\infty ({\mathbb{M}})$.

•

$d_2:\Lambda ^2({\mathbb{V}})\to \Lambda ^3({\mathbb{V}})$ corresponds to the row-wise divergence $C^\infty ({\mathbb{M}})\to C^\infty ({\mathbb{V}})$.

•

$S_2:\Lambda ^2({\mathbb{V}})\to \Lambda ^3(\mathbb{K})$ corresponds to the operator $- 2 \operatorname {skw}:C^\infty ({\mathbb{M}})\to C^\infty (\mathbb{K})$.

Thus, modulo these identifications and the (unimportant) constant factor in the last identification, 3.3 and 4.11 are identical. Hence we have established the following result.

Theorem 4.1.

When the de Rham complex 2.2 is exact, (i.e., the domain is contractible), then so is the elasticity complex 3.3.

To end this section, we return to the operator $S:\Lambda ^k({\mathbb{V}})\to \Lambda ^{k+1}(\mathbb{K})$ defined by $S=dK-Kd$. Let $K':\Lambda ^k(\mathbb{K})\to \Lambda ^k({\mathbb{V}})$ be the adjoint of $K$ (with respect to the Euclidean inner product on ${\mathbb{V}}$ and the Frobenius inner product on $\mathbb{K}$), which is given by $(K'\omega )_x(v_1,\ldots ,v_k)=-2\omega _x(v_1,\ldots ,v_k)x$. Define $S':\Lambda ^k(\mathbb{K})\to \Lambda ^{k+1}({\mathbb{V}})$ by $S'=dK'-K'd$. Recall that the wedge product $\wedge :\Lambda ^k\times \Lambda ^l\to \Lambda ^{k+l}$ is given by

$$\begin{align*} &(\omega \wedge \mu )(v_1,\ldots ,v_{k+l})\\ &\qquad =\sum (\operatorname {sign}\sigma )\omega (v_{\sigma _1},\ldots ,v_{\sigma _k}) \mu (v_{\sigma _{k+1}},\ldots ,v_{\sigma _{k+l}}),\quad \omega \in \Lambda ^k, \mu \in \Lambda ^l, v_i\in V, \end{align*}$$

where the sum is over the set of all permutations of $\{1,\ldots ,k+l\}$, for which $\sigma _1<\sigma _2<\cdots <\sigma _k$ and $\sigma _{k+1}<\sigma _{j+2}<\cdots <\sigma _{k+l}$. This extends as well to differential forms with values in an inner product space, using the inner product to multiply the terms inside the summation. Using the Leibniz rule 4.6, we have

$$\begin{equation} (S\omega )\wedge \mu = (-1)^{k}\omega \wedge S'\mu , \quad \omega \in \Lambda ^k({\mathbb{V}}),\nobreakspace \mu \in \Lambda ^l(\mathbb{K}). {\tag{4.13}} \end{equation}$$

We thus have

$$\begin{equation*} dK\omega \wedge \mu = (-1)^{k+1}K\omega \wedge d\mu +d(K\omega \wedge \mu ) = (-1)^{k+1}\omega \wedge K'd\mu +d(\omega \wedge K'\mu ), \end{equation*}$$

and

$$\begin{equation*} Kd\omega \wedge \mu =d\omega \wedge K'\mu = (-1)^{k+1}\omega \wedge dK'\mu +d(\omega \wedge K'\mu ). \end{equation*}$$

Subtracting these two expressions gives 4.13.

For later reference, we note that, analogously to 4.7, we have

$$\begin{equation} (S'\omega )(v_1,\ldots ,v_{k+1}) = -2\sum _{j=1}^{k+1}(-1)^{j+1} \omega (v_1,\ldots ,\hat{v}_j,\ldots v_{k+1}) v_j,\quad \omega \in \Lambda ^k(\Omega ;\mathbb{K}). {\tag{4.14}} \end{equation}$$

5. The discrete construction

In this section we derive a discrete version of the elasticity sequence by adapting the construction of the previous section. To carry out the construction, we will use two discretizations of the de Rham sequence. For $k=0,1,2,3$, let $\Lambda ^k_h$ denote a finite-dimensional space of $H\Lambda ^k$ for which $d\Lambda ^k_h\subset \Lambda ^{k+1}_h$, and for which there exist projections $\Pi _h=\Pi ^k_h:\Lambda ^k\to \Lambda ^k_h$ which make the following diagram commute:

$$\begin{equation} \begin{CD} \mathbb{R}@.\,\hookrightarrow \,@. \Lambda ^0 @>d>> \Lambda ^1 @>d>> \Lambda ^2 @>d>> \Lambda ^3 @>>> 0\\@. @. @VV\Pi _hV @VV\Pi _hV @VV\Pi _hV @VV\Pi _hV \\\mathbb{R}@.\,\hookrightarrow \,@. \Lambda _h^0 @>d>> \Lambda _h^1 @>d>> \Lambda _h^2 @>d>> \Lambda _h^3 @>>> 0 \end{CD} {\tag{5.1}} \end{equation}$$

This is simply the diagram 2.3 written in the language of differential forms. We do not make a specific choice of the discretization yet, but, as recalled in §2, there exist many such discrete de Rham complexes based on piecewise polynomials. In fact, as explained in 7, for each polynomial degree $r\ge 0$ we may choose $\Lambda ^3_h$ to be the space of all piecewise polynomial $3$-forms with respect to some simplicial decomposition of $\Omega$, and construct four such diagrams. We make the assumption that ${\mathcal{P}}_1(\Omega )\subset \Lambda ^0_h$, which is true in all the cases mentioned.

Let $\tilde{\Lambda }^k_h$ be a second set of finite dimensional spaces with corresponding projection operators $\tilde{\Pi }_h$ enjoying the same properties, giving us a second discretization of the de Rham sequence. Supposing a compatibility condition between these two discretizations, which we describe below, we shall construct a discrete elasticity complex.

We start with the complex

$$\begin{equation} \begin{CD} \mathbb{K}\times {\mathbb{V}}@.\,\hookrightarrow \,@. \Lambda ^0_h(\mathbb{K})\times \tilde{\Lambda }^0_h({\mathbb{V}}) @>d>> \cdots @>d>> \Lambda ^3_h(\mathbb{K})\times \tilde{\Lambda }^3_h({\mathbb{V}}) @>>> 0 \end{CD} {\tag{5.2}} \end{equation}$$

where $\Lambda ^k_h(\mathbb{K})$ denotes the $\mathbb{K}$-valued analogue of $\Lambda ^k_h$ and similarly for $\tilde{\Lambda }^k_h({\mathbb{V}})$. For brevity, we henceforth write $\Lambda ^k_h({\mathbb{W}})$ for $\Lambda ^k_h(\mathbb{K})\times \tilde{\Lambda }^k_h({\mathbb{V}})$. As a discrete analogue of the operator $K$, we define $K_h : \tilde{\Lambda }_h^k({\mathbb{V}}) \to \Lambda _h^k(\mathbb{K})$ by $K_h = \Pi _h K$ where $\Pi _h$ is the interpolation operator onto $\Lambda _h^k(\mathbb{K})$.

Next define $S_h =S_{k,h}: \tilde{\Lambda }_h^k({\mathbb{V}}) \to \Lambda _h^{k+1}(\mathbb{K})$ by $S_h = d K_h - K_hd$, for $k=0,1,2$. Observe that the discrete version of 4.8,

$$\begin{equation} dS_h =- S_hd, {\tag{5.3}} \end{equation}$$

follows exactly as in the continuous case. From the commutative diagram 5.1, we see that

$$\begin{equation*} S_h = d\Pi _h K - \Pi _hKd = \Pi _h(dK - Kd) = \Pi _h S. \end{equation*}$$

Continuing to mimic the continuous case, we define the automorphism $\Phi _h$ on $\Lambda ^k_h({\mathbb{W}})$ by

$$\begin{equation*} \Phi _h(\omega ,\mu ) = (\omega + K_h\mu , \mu ), \end{equation*}$$

and the operator ${\mathcal{A}}_h : \Lambda _h^k({\mathbb{W}}) \to \Lambda _h^{k+1}({\mathbb{W}})$ by ${\mathcal{A}}_h = \Phi _h d \Phi _h^{-1}$, which leads to

$$\begin{equation*} {\mathcal{A}}_h(\omega ,\mu ) = (d\omega - S_h \mu , d\mu ). \end{equation*}$$

Thanks to the assumption that ${\mathcal{P}}_1\subset \Lambda ^0_h$, we have $\Phi _h({\mathbb{W}})=\Phi ({\mathbb{W}})$. Hence, inserting the isomorphisms $\Phi _h$ into 5.2, we obtain

$$\begin{equation} \begin{CD} \Phi ({\mathbb{W}}) @.\,\hookrightarrow \,@. \Lambda _h^0({\mathbb{W}}) @>{\mathcal{A}}_h>> \Lambda _h^1 ({\mathbb{W}})@>{\mathcal{A}}_h>> \Lambda _h^2({\mathbb{W}}) @>{\mathcal{A}}_h>> \Lambda _h^3({\mathbb{W}})@>>> 0. \end{CD} {\tag{5.4}} \end{equation}$$

In analogy to the continuous case, we define

$$\begin{equation*} \Gamma ^1_h = \{\,(\omega , \mu ) \in \Lambda ^1_h({\mathbb{W}})\,|\, d\omega = S_{1,h}\mu \,\} , \quad \Gamma ^2_h = \{\,(\omega , \mu ) \in \Lambda ^2_h({\mathbb{W}})\,|\, \omega = 0\,\}. \end{equation*}$$

As in the continuous case, we can identify $\Gamma _h^2$ with $\tilde{\Lambda }^2_h({\mathbb{V}})$, but, unlike in the continuous case, we cannot identify $\Gamma ^1_h$ with $\Lambda ^1_h(\mathbb{K})$, since we do not require that $S_{1,h}$ be invertible (it is not in the applications). Hence, in order to derive the analogue of the diagram 4.10 we require a surjectivity assumption:

$$\begin{equation} \text{\textit{The operator $S_{1,h}$ maps $\tilde{\Lambda }_h^1({\mathbb{V}})$ onto $\Lambda _h^2(\mathbb{K})$.}} {\tag{5.5}} \end{equation}$$

Under this assumption, the operator $S_h=S_{1,h}$ has a right inverse $S_h^\dagger$ mapping $\Lambda _h^2(\mathbb{K})$ into $\Lambda _h^1({\mathbb{V}})$. This allows us to define discrete counterparts of the projection operators $\pi ^1$ and $\pi ^2$ by

$$\begin{equation*} \pi _h^1(\omega , \mu ) = (\omega , \mu - S_h^{\dagger }S_h \mu + S_h^{\dagger }d \omega ), \quad \pi _h^2(\omega , \mu ) = (0 , \mu + d S_h^{\dagger } \omega ), \end{equation*}$$

and obtain the discrete analogue of 4.10:

$$\begin{equation} \begin{CD} \Phi ({\mathbb{W}}) @.\,\hookrightarrow \,@. \Lambda _h^0({\mathbb{W}}) @>{\mathcal{A}}_h>> \Lambda _h^1 ({\mathbb{W}})@>{\mathcal{A}}_h>> \Lambda _h^2({\mathbb{W}}) @>{\mathcal{A}}_h>> \Lambda _h^3({\mathbb{W}})@>>> 0\\@. @. @VVidV @VV\pi _h^1V @VV\pi _h^2V @VVidV\\\Phi ({\mathbb{W}}) @.\,\hookrightarrow \,@. \Lambda _h^0({\mathbb{W}}) @>{\mathcal{A}}_h>> \Gamma _h^1@>{\mathcal{A}}_h>> \Gamma _h^2@>{\mathcal{A}}_h>> \Lambda _h^3({\mathbb{W}})@>>> 0 \end{CD} {\tag{5.6}} \end{equation}$$

It is straightforward to check that this diagram commutes. For example, if $(\omega ,\mu ) \in \Lambda _h^0({\mathbb{W}})$, then

$$\begin{align*} \pi _h^1 {\mathcal{A}}_h(\omega ,\mu ) &= \pi _h^1 (d\omega - S_h \mu ,d\mu ) = (d\omega - S_h \mu , d\mu - S_h^{\dagger }S_h d\mu + S_h^{\dagger }d [d\omega -S_h \mu ]) \\ &= (d\omega - S_h \mu , d\mu - S_h^{\dagger }[S_h d\mu + d S_h \mu ]) = {\mathcal{A}}_h(\omega ,\mu ), \end{align*}$$

where the last equality follows from 5.3. Thus the bottom row of 5.6 is a subcomplex of the top row, and the vertical maps are commuting projections. In particular, when the top row is exact, so is the bottom. Thus we have established the following result.

Theorem 5.1.

For $k=0, \ldots , 3$, let $\Lambda ^k_h$ be a finite dimensional subspace of $H \Lambda ^k$ for which $d \Lambda ^k_h \subset \Lambda ^{k+1}_h$ and for which there exist projections $\Pi _h = \Pi ^k_h: \Lambda ^k \to \Lambda ^k_h$ that make the diagram 5.1 commute. Let $\tilde{\Lambda }^k_h$ be a second set of finite dimensional spaces with corresponding projection operators $\tilde{\Pi }^k_h$ with the same properties. If $S_{1,h} := d \Pi ^1_h K - \Pi ^2_h K d$ maps $\tilde{\Lambda }^1_h({\mathbb{V}})$ onto $\Lambda ^2_h(\mathbb{K})$, and the bottom row of 5.1 is exact for both sequences $\Lambda ^k_h$ and $\tilde{\Lambda }^k_h$, then the discrete elasticity sequence given by the bottom row of 5.6 is also exact.

The exactness of the bottom row of 5.6 suggests that the following choice of finite element spaces will lead to a stable discretization of 2.1:

$$\begin{equation*} \Sigma _h \sim \tilde{\Lambda }^2_h({\mathbb{V}}), \quad V_h \sim \tilde{\Lambda }^3_h({\mathbb{V}}), \quad Q_h \sim \Lambda ^3_h(\mathbb{K}). \end{equation*}$$

In the next section we will make specific choices for the discrete de Rham complexes, and then verify the stability in the following section.

For use in the next section, we state the following result, giving a sufficient condition for the key requirement that $S_{1,h}$ be surjective.

Proposition 5.2.

If the diagram

$$\begin{equation} \begin{CD} \Lambda ^1({\mathbb{V}}) @>S_1>> \Lambda ^2(\mathbb{K}) \\@V\tilde{\Pi }_h^1VV @V\Pi _h^2VV \\\tilde{\Lambda }_h^1({\mathbb{V}}) @>S_{1,h}>> \Lambda _h^2(\mathbb{K}) \end{CD} {\tag{5.7}} \end{equation}$$

commutes, then $S_{1,h}$ is surjective.

6. A family of discrete elasticity complexes

In this section, we present a family of examples of the general discrete construction presented in the previous section by choosing specific discrete de Rham complexes. These furnish a family of discrete elasticity complexes, indexed by an integer degree $r\ge 0$. In the next section we use these complexes to derive finite elements for elasticity. In the lowest order case, the method will require only piecewise linear functions to approximate stresses and piecewise constants to approximate the displacements and multipliers.

We begin by recalling the two principal families of piecewise polynomial spaces of differential forms, following the presentation in 7. We henceforth assume that the domain $\Omega$ is a contractible polyhedron. Let ${\mathcal{T}}_h$ be a triangulation of $\Omega$. Let ${\mathcal{T}}_h$ be a triangulation of $\Omega$ by tetrahedra, and set

$$\begin{align*} {\mathcal{P}}_r\Lambda ^k({\mathcal{T}}_h)&=\{\,\omega \in H\Lambda ^k(\Omega )\,\,|\,\, \omega _{|_T}\in {\mathcal{P}}_r\Lambda ^k(T)\quad \forall T\in {\mathcal{T}}_h\,\}, \\ {\mathcal{P}}_r^+\Lambda ^k({\mathcal{T}}_h)&=\{\,\omega \in H\Lambda ^k(\Omega )\,\,|\,\, \omega _{|_T}\in {\mathcal{P}}_r^+\Lambda ^k(T)\quad \forall T\in {\mathcal{T}}_h\,\}. \end{align*}$$

Here ${\mathcal{P}}_r^+\Lambda ^k(T):={\mathcal{P}}_r\Lambda ^k(T)+\kappa {\mathcal{P}}_r\Lambda ^{k+1}(T)$ where $\kappa :\Lambda ^{k+1}(T)\to \Lambda ^k(T)$ is the Koszul differential defined by

$$\begin{equation*} (\kappa \omega )_x(v^1, \cdots , v^k) = \omega _x(x,v^1, \cdots , v^k). \end{equation*}$$

The spaces ${\mathcal{P}}_r^+\Lambda ^0({\mathcal{T}}_h)={\mathcal{P}}_{r+1}\Lambda ^0({\mathcal{T}}_h)$ correspond to the usual degree $r+1$ Lagrange piecewise polynomial subspaces of $H^1$, and the spaces ${\mathcal{P}}_r^+\Lambda ^3({\mathcal{T}}_h)={\mathcal{P}}_r\Lambda ^3({\mathcal{T}}_h)$ correspond to the usual degree $r$ subspace of discontinuous piecewise polynomials in $L^2(\Omega )$. For $k=1$ and $2$, the spaces ${\mathcal{P}}_r^+\Lambda ^k({\mathcal{T}}_h)$ correspond to the discretizations of $H(\operatorname {curl})$ and $H(\operatorname {div})$, respectively, presented by Nédélec in 25, and the spaces ${\mathcal{P}}_r\Lambda ^k({\mathcal{T}}_h)$ are the ones presented by Nédélec in 26. An element $\omega \in {\mathcal{P}}_r\Lambda ^k({\mathcal{T}}_h)$ is uniquely determined by the following quantities:

$$\begin{equation} \int _f \omega \wedge \zeta , \quad \zeta \in {\mathcal{P}}_{r-d-1+k}^+\Lambda ^{d-k}(f), \quad f\in \Delta _d({\mathcal{T}}_h), \quad k\le d \le 3. {\tag{6.1}} \end{equation}$$

Here $\Delta _d({\mathcal{T}}_h)$ is the set of vertices, edges, faces, or tetrahedra in the mesh ${\mathcal{T}}_h$, for $d=0,1,2,3$, respectively, and for $r<0$, we interpret ${\mathcal{P}}_r^+\Lambda ^k(T)={\mathcal{P}}_r\Lambda ^k(T)=0$. Note that for $\omega \in \Lambda ^k$, $\omega$ naturally restricts on the face $f$ to an element of $\Lambda ^k(f)$. Therefore, for $\zeta \in \Lambda ^{d-k}(f)$, the wedge product $\omega \wedge \zeta$ belongs to $\Lambda ^d(f)$ and hence the integral of $\omega \wedge \zeta$ on the $d$-dimensional face $f$ of $T$ is a well-defined and natural quantity. Using the quantities in 6.1, we obtain a projection operator from $\Lambda ^k$ to ${\mathcal{P}}_r\Lambda ^k({\mathcal{T}}_h)$.

Similarly, an element $\omega \in {\mathcal{P}}_r^+\Lambda ^k({\mathcal{T}}_h)$ is uniquely determined by

$$\begin{equation} \int _f \omega \wedge \zeta , \quad \zeta \in {\mathcal{P}}_{r-d+k}\Lambda ^{d-k}(f), \quad f\in \Delta _d({\mathcal{T}}_h), \quad k\le d \le 3, {\tag{6.2}} \end{equation}$$

and so these quantities determine a projection.

If ${\mathbb{X}}$ is a vector space, we use the notation ${\mathcal{P}}_r\Lambda ^k({\mathcal{T}}_h;{\mathbb{X}})$ and ${\mathcal{P}}^+_r\Lambda ^k({\mathcal{T}}_h;{\mathbb{X}})$ to denote the corresponding spaces of piecewise polynomial differential forms with values in ${\mathbb{X}}$. Furthermore, if ${\mathbb{X}}$ is an inner product space, the corresponding degrees of freedom are given by 6.1 and 6.2, but where the test spaces are replaced by the corresponding ${\mathbb{X}}$ valued spaces.

To carry out the construction described in the previous section we need to choose the two sets of spaces $\Lambda ^k_h$ and $\tilde{\Lambda }^k_h$ for $k=0,1,2,3$. We fix $r\ge 0$ and set $\Lambda ^k_h={\mathcal{P}}_r^+\Lambda ^k({\mathcal{T}}_h)$, $k=0,1,2,3$, and $\tilde{\Lambda }^0_h={\mathcal{P}}_{r+2}\Lambda ^0({\mathcal{T}}_h)$, $\tilde{\Lambda }^1_h={\mathcal{P}}_{r+1}^+\Lambda ^1({\mathcal{T}}_h)$, $\tilde{\Lambda }^2_h= {\mathcal{P}}_{r+1}\Lambda ^2({\mathcal{T}}_h)$, and $\tilde{\Lambda }^3_h=P_r\Lambda ^3_h$. As explained in 7, both these choices give a discrete de Rham sequence with commuting projections, i.e., a diagram like 5.1 makes sense and is commutative.

We establish the key surjectivity assumption for our choice of spaces by verifying the commutativity of 5.7.

Lemma 6.1.

Let $\tilde{\Lambda }^1_h({\mathbb{V}})={\mathcal{P}}_{r+1}^+\Lambda ^1({\mathcal{T}}_h;{\mathbb{V}})$ and $\Lambda ^2_h(\mathbb{K}) ={\mathcal{P}}_r^+\Lambda ^2({\mathcal{T}}_h;\mathbb{K})$ with projections $\tilde{\Pi }^1_h$, $\Pi ^2_h$ defined via the corresponding vector-valued moments of the form 6.1 and 6.2. If $S_{1,h}=\Pi ^2_h S_1$ then

$$\begin{equation} S_{1,h}\tilde{\Pi }_h^1=\Pi _h^2 S_1, {\tag{6.3}} \end{equation}$$

and so $S_{1,h}$ is surjective.

Proof.

We must show that $(\Pi ^2_h S_1-S_{1,h}\tilde{\Pi }^1_h)\sigma =0$ for $\sigma \in \Lambda ^1({\mathbb{V}})$. Defining $\omega =(I-\tilde{\Pi }^1_h)\sigma$, the required condition becomes $\Pi ^2_hS_1\omega =0$. Since $\tilde{\Pi }^1_h\omega =0$, we have

$$\begin{equation} \int _f \omega \wedge \zeta = 0 , \quad \zeta \in {\mathcal{P}}_{r+2-d}\Lambda ^{d-1}(f;{\mathbb{V}}), \quad f\in \Delta _d({\mathcal{T}}_h), \quad 2\le d\le 3, {\tag{6.4}} \end{equation}$$

(in fact 6.4 holds for $d=1$ as well, but this is not used here). We must show that 6.4 implies

$$\begin{equation} \int _f S_1\omega \wedge \mu = 0 , \quad \mu \in {\mathcal{P}}_{r+2-d}\Lambda ^{d-2}(f;\mathbb{K}),\quad f\in \Delta _d({\mathcal{T}}_h), \quad 2\le d\le 3. {\tag{6.5}} \end{equation}$$

From 4.13, we have $S_1\omega \wedge \mu = - \omega \wedge \zeta$, where $\zeta =S'_{d-2}\mu \in {\mathcal{P}}_{r+2-d}\Lambda ^{d-1}(f;{\mathbb{V}})$ for $\mu \in {\mathcal{P}}_{r+2-d}\Lambda ^{d-2}(f;\mathbb{K})$, as is evident from 4.14. Hence 6.5 follows from 6.4.

■

7. Stable mixed finite elements for elasticity

Based on the discrete elasticity complexes just constructed, we obtain mixed finite element spaces for the formulation 2.1 of the elasticity equations by choosing $\Sigma _h$, $V_h$, and $Q_h$ as the spaces of matrix and vector fields corresponding to appropriate spaces of forms in the $\mathbb{K}$- and ${\mathbb{V}}$-valued de Rham sequences used in the construction. Specifically, these are

$$\begin{equation} \Sigma _h \sim {\mathcal{P}}_{r+1}\Lambda ^2({\mathcal{T}}_h;{\mathbb{V}}), \quad V_h \sim {\mathcal{P}}_r\Lambda ^3({\mathcal{T}}_h;{\mathbb{V}}), \quad Q_h \sim {\mathcal{P}}_r\Lambda ^3({\mathcal{T}}_h;\mathbb{K}). {\tag{7.1}} \end{equation}$$

In other terminology, $\Sigma _h$ may be thought of as the product of three copies of the Nédélec $H(\operatorname {div})$ space of the second kind of degree $r+1$, and $V_h$ and $Q_h$ are spaces of all piecewise polynomials of degree at most $r$ with values in $\mathbb{K}$ and ${\mathbb{V}}$, respectively, with no imposed interelement continuity. In this section, we establish stability and convergence for this finite element method.

The stability of the method requires the two conditions (A$1$) and (A$2$) stated in §2. The first condition is obvious since, by construction, $\operatorname {div}\Sigma _h\subset V_h$, i.e., $d{\mathcal{P}}_{r+1}\Lambda ^2({\mathcal{T}}_h;{\mathbb{V}})\subset {\mathcal{P}}_r\Lambda ^3({\mathcal{T}}_h;{\mathbb{V}})$. The condition (A$2$) is more subtle. We will prove a stronger version, namely,

 (A$2'$)

for all nonzero $(v,q) \in V_h \times Q_h$, there exists nonzero $\tau \in \Sigma _h$ with $\operatorname {div}\tau =v$, $2\Pi _{Q_h}\operatorname {skw}\tau = q$ and$$\begin{equation*} \|\tau \|_{\operatorname {div}} \le c (\|v\|+\|q\|), \end{equation*}$$

where $\Pi _{Q_h}$ is the $L^2$ projection into $Q_h$ and $c$ is a constant.

Recalling that $\Gamma ^2_h=0\times {\mathcal{P}}_{r+1}\Lambda ^2({\mathcal{T}}_h;{\mathbb{V}})$ and ${\mathcal{A}}_h(0, \sigma ) = (-S_{2,h} \sigma , d \sigma )$, and that the operator $S_2$ corresponds on the matrix level to $-2\operatorname {skw}$, we restate condition (A$2'$) in the language of differential forms.

Theorem 7.1.

Given that $(\omega , \mu ) \in {\mathcal{P}}_r\Lambda ^3({\mathcal{T}}_h;\mathbb{K}) \times {\mathcal{P}}_r\Lambda ^3({\mathcal{T}}_h;{\mathbb{V}})$, there exists $\sigma \in {\mathcal{P}}_{r+1}\Lambda ^2({\mathcal{T}}_h;{\mathbb{V}})$ such that ${\mathcal{A}}_h(0, \sigma ) = (\omega , \mu )$ and

$$\begin{equation} \|\sigma \|_{H\Lambda } \le c (\|\omega \| + \|\mu \|), {\tag{7.2}} \end{equation}$$

where the constant $c$ is independent of $\omega ,\mu$ and $h$.

Before proceeding to the proof, we need to establish some bounds on projection operators. We do this for the corresponding scalar-valued spaces. The extensions to vector-valued spaces are straightforward. First we claim that

$$\begin{equation} \|\tilde{\Pi }^2_h\eta \|\le c\|\eta \|_1 \quad \forall \eta \in H^1\Lambda ^2,\quad \|\Pi ^3_h\omega \|\le c\|\omega \|_0\quad \forall \omega \in H^1\Lambda ^3. {\tag{7.3}} \end{equation}$$

Here the constant may depend on the shape regularity of the mesh, but not on the meshsize. The second bound is obvious (with $c=1$), since $\Pi ^3_h$ is just the $L^2$ projection. The first bound follows by a standard scaling argument. Namely, let $\hat{T}$ denote the reference simplex. For any $\hat{\beta }\in {\mathcal{P}}_{r+1}\Lambda ^2(\hat{T})$, we have

$$\begin{equation} \|\hat{\beta }\|_{0,\hat{T}} \le c (\sum _{\hat{f}}\sum _{\hat{\mu }}|\int _{\hat{f}} \hat{\beta }\wedge \hat{\mu }| +\sum _{\hat{\zeta }}|\int _{\hat{T}}\hat{\beta }\wedge \hat{\zeta }|), {\tag{7.4}} \end{equation}$$

where $\hat{f}$ ranges over the faces of $\hat{T}$, $\hat{\mu }$ over a basis for ${\mathcal{P}}^+_r\Lambda ^0(\hat{f})$, and $\hat{\zeta }$ over a basis for ${\mathcal{P}}^+_{r-1}\Lambda ^1(\hat{T})$. This is true because the integrals on the right hand side of 7.4 form a set of degrees of freedom for $\hat{\beta }\in {\mathcal{P}}_{r+1}\Lambda ^2(\hat{T})$ (see 6.1), and so we may use the equivalence of all norms on this finite dimensional space. We apply this result with $\hat{\beta }=\hat{\Pi }^2_h\hat{\eta }$, where $\hat{\Pi }^2_h$ is the projection defined to preserve the integrals on the right hand side of 7.4. It follows that

$$\begin{equation*} \|\hat{\Pi }^2_h\hat{\eta }\|_{0,\hat{T}} \le c(\sum _{\hat{f}}\sum _{\hat{\mu }}|\int _{\hat{f}} \hat{\eta }\wedge \hat{\mu }| +\sum _{\hat{\zeta }}|\int _{\hat{T}}\hat{\eta }\wedge \hat{\zeta }|)\le c\|\hat{\eta }\|_{1,\hat{T}}, \end{equation*}$$

where we have used a standard trace inequality in the last step. Next, if $T$ is an arbitrary simplex and $\eta \in H^1\Lambda ^2(T)$, we map the reference simplex $\hat{T}$ onto $T$ by an affine map $\hat{x}\mapsto B\hat{x} +b$, and define $\hat{\eta }\in H^1\Lambda ^2(\hat{T})$ by

$$\begin{equation*} \hat{\eta }_{\hat{x}}(\hat{v}_1, \hat{v}_2) = \eta _x(B\hat{v}_1, B \hat{v}_2), \end{equation*}$$

for any $x=B\hat{x}+b\in T$ and any vectors $\hat{v}_1, \hat{v}_2$. It is easy to check that $\widehat{\tilde{\Pi }^2_h\eta }=\hat{\Pi }^2_h\hat{\eta }$, and that

$$\begin{equation*} \|\tilde{\Pi }^2_h\eta \|_{0,T}\le c \|\hat{\Pi }^2_h\hat{\eta }\|_{0,\hat{T}} \le c\|\hat{\eta }\|_{1,\hat{T}}\le c (\|\eta \|_{0,T} + h|\eta |_{1,T}) \le c\|\eta \|_{1,T}. \end{equation*}$$

Squaring and adding this over all the simplices in the mesh ${\mathcal{T}}_h$ gives the first bound in 7.3.

We also need a bound on the projection of a form in $H^1\Lambda ^1$ into $\tilde{\Lambda }^1_h={\mathcal{P}}^+_{r+1}\Lambda ^1({\mathcal{T}}_h)$. However, the projection operator $\tilde{\Pi }^1_h$ is not bounded on $H^1$, because its definition involves integrals over edges. A similar problem has arisen before (see, e.g., 10), and we use the same remedy. Namely we start by defining an operator $\tilde{\Pi }^1_{0h}:H^1\Lambda ^1\to {\mathcal{P}}^+_{r+1}\Lambda ^1({\mathcal{T}}_h)$ by the conditions

$$\begin{align} \int _T\tilde{\Pi }^1_{0h}\omega \wedge \zeta &=\int _T\omega \wedge \zeta , \quad \zeta \in {\mathcal{P}}_{r-1}\Lambda ^2(T),\quad T\in {\mathcal{T}}_h, {\tag{7.5}}\\ \int _f\tilde{\Pi }^1_{0h}\omega \wedge \zeta &=\int _f\omega \wedge \zeta , \quad \zeta \in {\mathcal{P}}_r\Lambda ^1(f),\quad f\in \Delta _2({\mathcal{T}}_h), {\tag{7.6}}\\ \int _e\tilde{\Pi }^1_{0h}\omega \wedge \zeta &=0, \quad \zeta \in {\mathcal{P}}_{r+1}\Lambda ^0(e), \quad e\in \Delta _1({\mathcal{T}}_h). {\tag{7.7}} \end{align}$$

Note that, in contrast to $\tilde{\Pi }^1_h$, in the definition of $\tilde{\Pi }^1_{0h}$, we have set the troublesome edge degrees of freedom to zero. Let $\hat{\Pi }^1_0:H^1\Lambda ^1(\hat{T})\to {\mathcal{P}}^+_{r+1}\Lambda ^1(\hat{T})$ be defined analogously on the reference element.

Now for $\hat{\rho }\in H^1\Lambda ^1(\hat{T})$, $\hat{d}\hat{\Pi }^1_0\hat{\rho }\in {\mathcal{P}}_{r+1}\Lambda ^2(\hat{T})$, so

$$\begin{equation*} \|\hat{d}\hat{\Pi }^1_0\hat{\rho }\|_{0,\hat{T}} \le c (\sum _{\hat{f}}\sum _{\hat{\mu }} |\int _{\hat{f}} \hat{d}\hat{\Pi }^1_0\hat{\rho }\wedge \hat{\mu }| +\sum _{\hat{\zeta }}|\int _{\hat{T}}\hat{d}\hat{\Pi }^1_0\hat{\rho }\wedge \hat{\zeta }|), \end{equation*}$$

where again $\hat{f}$ ranges over the faces of $\hat{T}$, $\mu$ over a basis of ${\mathcal{P}}^+_r\Lambda ^0(\hat{f})$, and $\zeta$ over a basis of ${\mathcal{P}}^+_{r-1}\Lambda ^1(\hat{T})$. But

$$\begin{equation*} \int _{\hat{f}} \hat{d}\hat{\Pi }^1_0\hat{\rho }\wedge \hat{\mu }=\int _{\hat{f}}\hat{\Pi }^1_0\hat{\rho }\wedge \hat{d}\hat{\mu }=\int _{\hat{f}}\hat{\rho }\wedge \hat{d}\hat{\mu }, \end{equation*}$$

where we have used Stokes’ theorem and the fact that the vanishing of the edge quantities in the definition of $\hat{\Pi }^1_0$ to obtain the first equality, and the face degrees of freedom entering the definition of $\hat{\Pi }^1_0$ to obtain the second. Similarly,

$$\begin{equation*} \int _{\hat{T}}\hat{d}\hat{\Pi }^1_0\hat{\rho }\wedge \hat{\zeta }=\int _{\hat{T}}\hat{\Pi }^1_0\hat{\rho }\wedge \hat{d}\hat{\zeta }+\int _{\partial \hat{T}}\hat{\Pi }^1_0\hat{\rho }\wedge \hat{\zeta }=\int _{\hat{T}}\hat{\rho }\wedge \hat{d}\hat{\zeta }+\int _{\partial \hat{T}}\hat{\rho }\wedge \hat{\zeta }=\int _{\hat{T}}\hat{d}\hat{\rho }\wedge \hat{\zeta }. \end{equation*}$$

It follows that

$$\begin{equation*} \|\hat{d}\Pi ^1_0\hat{\rho }\|_{0,\hat{T}} \le c \|\hat{\rho }\|_{1,\hat{T}}, \quad \rho \in H^1\Lambda ^1(\hat{T}). \end{equation*}$$

When we scale this result to an arbitrary simplex and add over the mesh, we obtain

$$\begin{equation*} \|d\Pi ^1_{0h}\rho \| \le c (h^{-1}\|\rho \| + \|\rho \|_1),\quad \rho \in H^1\Lambda ^1(\Omega ). \end{equation*}$$

To remove the problematic $h^{-1}$ in the last estimate, we introduce the Clement interpolant $R_h$ mapping $H^1\Lambda ^1$ into continuous piecewise linear $1$-forms (still following 10). Then

$$\begin{equation*} \|\rho -R_h\rho \|+h\|\rho -R_h\rho \|_1 \le c h \|\rho \|_1, \quad \rho \in H^1\Lambda ^1. \end{equation*}$$

Defining $\bar{\Pi }^1_h:H^1\Lambda ^1\to {\mathcal{P}}^+_{r+1}\Lambda ^1_h$ by

$$\begin{equation} \bar{\Pi }^1_h = \tilde{\Pi }^1_{0h}(I-R_h)+ R_h, {\tag{7.8}} \end{equation}$$

we obtain

$$\begin{align*} \|d\bar{\Pi }^1_h\rho \|&\le \|d\tilde{\Pi }^1_{0h}(I-R_h)\rho \| + \|dR_h\rho \|\\ &\le c (h^{-1}\|(I-R_h)\rho \| + \|(I-R_h)\rho \|_1+ \|dR_h\rho \|) \le c\|\rho \|_1. \end{align*}$$

Thus we have shown that

$$\begin{equation} \|d\bar{\Pi }^1_h\rho \|\le c\|\rho \|_1, \quad \rho \in H^1\Lambda ^1. {\tag{7.9}} \end{equation}$$

Having modified $\tilde{\Pi }^1_h$ to obtain the bounded operator $\bar{\Pi }^1_h$, we now verify that the key property 6.3 in Lemma 6.1 carries over to

$$\begin{equation} S_{1,h}\bar{\Pi }_h^1=\Pi _h^2 S_1, {\tag{7.10}} \end{equation}$$

where we now use the vector-valued forms of the projection operators. It follows easily from 7.8, 7.5, and 7.6 that 6.4 holds with $\omega =(I-\bar{\Pi }^1_h)\sigma$, so that the proof of 7.10 is the same as for 6.3.

We can now give the proof of Theorem 7.1.

Proof of Theorem 7.1.

Given $\mu \in {\mathcal{P}}_r\Lambda ^3({\mathcal{T}}_h;{\mathbb{V}})$ there exists $\eta \in H^1\Lambda ^2({\mathbb{V}})$ such that $d\eta =\mu$ with the bound $\|\eta \|_1\le c\|\mu \|$ (since $d$ maps $H^1\Lambda ^2$ onto $L^2\Lambda ^3$). Similarly, given $\omega \in {\mathcal{P}}_r\Lambda ^3({\mathcal{T}}_h;\mathbb{K})$ there exists $\tau \in H^1\Lambda ^2(\mathbb{K})$ with $d\tau =\omega +S_{2,h}\tilde{\Pi }^2_h\eta$ with the bound $\|\tau \|_1\le c\|\omega +S_{2,h}\tilde{\Pi }^2_h\eta \|$. Let $\rho =S_1^{-1}\tau$ (recall that $S_1$ is an isomorphism) and set

$$\begin{equation*} \sigma = d\bar{\Pi }^1_h\rho +\tilde{\Pi }^2_h\eta . \end{equation*}$$

We will now show that ${\mathcal{A}}_h(0, \sigma ) = (\omega , \mu )$. From the definition of $\sigma$, we have

$$\begin{equation*} -S_{2,h}\sigma =-S_{2,h} d\bar{\Pi }^1_h\rho -S_{2,h}\tilde{\Pi }^2_h\eta . \end{equation*}$$

Then, using 5.3, 7.10, and the commutativity $d\Pi _h=\Pi _hd$, we see

$$\begin{multline*} S_{2,h} d\bar{\Pi }^1_h\rho =-d S_{1,h}\bar{\Pi }^1_h\rho =-d \Pi ^2_h S_1\rho \\ =-d\Pi ^2_h\tau =-\Pi ^3_h d\tau = -\Pi ^3_h(\omega +S_{2,h}\tilde{\Pi }^2_h\eta )= -\omega -S_{2,h}\tilde{\Pi }^2_h\eta . \end{multline*}$$

Combining, we get $-S_{2,h}\sigma =\omega$ as desired. Furthermore, from the commutativity $d\tilde{\Pi }_h=\tilde{\Pi }_hd$ and the definition of $\eta$, we get

$$\begin{equation*} d\sigma =d\tilde{\Pi }^2_h\eta =\tilde{\Pi }^3_hd\eta =\tilde{\Pi }^3_h\mu =\mu , \end{equation*}$$

and so we have established that ${\mathcal{A}}_h(0,\sigma )=(\mu ,\omega )$.

It remains to prove the bound 7.2. Using 7.3, we have

$$\begin{equation*} \|S_{2,h}\tilde{\Pi }^2_h\eta \|=\|\Pi ^3_h S_2\tilde{\Pi }^2_h\eta \| \le c \|S_2\tilde{\Pi }^2_h\eta \| \le c\|\tilde{\Pi }^2_h\eta \| \le c\|\eta \|_1\le c\|\mu \|. \end{equation*}$$

Thus $\|\rho \|_1\le c\|\tau \|_1\le c(\|\omega \|+\|\mu \|)$. Using 7.9, we then get $\|d\bar{\Pi }^1_h\rho \|\le c\|\rho \|_1\le c(\|\omega \|+\|\mu \|)$, and, using 7.3, that $\|\tilde{\Pi }^2_h\eta \|\le c\|\eta \|_1\le c\|\mu \|$. Therefore $\|\sigma \|\le c(\|\omega \|+\|\mu \|)$, while $\|d\sigma \|=\|\mu \|$, and thus we have the desired bound 7.2.

■

We have thus verified the stability conditions (A$1$) and (A$2$), and so may apply the standard theory of mixed methods (cf. 12, 13, 17, 20) and standard results about approximation by finite element spaces to obtain convergence and error estimates.

Theorem 7.2.

Suppose $(\sigma , u, p)$ is the solution of the elasticity system 1.5 and $(\sigma _h, u_h, p_h)$ is the solution of discrete system 2.1, where the finite element spaces $\Sigma _h$, $V_h$, and $Q_h$ are given by 7.1 for some integer $r\ge 0$. Then there is a constant $C$, independent of $h$, such that

$$\begin{gather*} \begin{split} &\|\sigma - \sigma _h\|_{\operatorname {div}} + \|u-u_h\| + \|p-p_h\| \\ &\qquad \le C \inf _{\tau _h \in \Sigma _h, v_h \in V_h, q_h \in Q_h} (\|\sigma - \tau _h\|_{\operatorname {div}} + \|u-v_h\| + \|p-q_h\|), \end{split} \\ \|\sigma - \sigma _h\| + \|p-p_h\| + \|u_h - \tilde{\Pi }^n_h u\| \le C (\|\sigma - \tilde{\Pi }^{n-1}_h \sigma \| + \|p - \Pi ^n_h p\|), \\ \|u-u_h\| \le C(\|\sigma - \tilde{\Pi }^{n-1}_h \sigma \| + \|p - \Pi ^n_h p\| + \|u- \tilde{\Pi }^n_h u\|), \\ \|d (\sigma - \sigma _h)\| = \|d \sigma - \tilde{\Pi }^n d \sigma \|. \end{gather*}$$

If $u$ and $\sigma$ are sufficiently smooth, then

$$\begin{equation*} \|\sigma - \sigma _h\| + \|u-u_h\| + \|p-p_h\| \!\le \! Ch^{r+1} \|u\|_{r+2}, \qquad \|\operatorname {div}(\sigma - \sigma _h)\| \!\le \!Ch^{r+1} \|\operatorname {div}\sigma \|_{r+1}. \end{equation*}$$

Remark.

Note that the errors $\|\sigma - \sigma _h\|$ and $\|u_h - \tilde{\Pi }^n_h u\|$ depend on the approximation of both $\sigma$ and $p$. For the choices made here, the approximation of $p$ is one order less than the approximation of $\sigma$, and thus we do not obtain improved estimates, as one does in the case of the approximation of Poisson’s equation, where the extra variable $p$ does not enter.

8. A simplified element

Recall that the lowest order element in the stable family described above, for a discretization based on 1.5, is of the form

$$\begin{equation*} \Sigma _h \sim {\mathcal{P}}_1 \Lambda ^2({\mathcal{T}}_h;{\mathbb{V}}), \qquad V_h \sim {\mathcal{P}}_0 \Lambda ^3({\mathcal{T}}_h;{\mathbb{V}}), \qquad Q_h \sim {\mathcal{P}}_0 \Lambda ^3({\mathcal{T}}_h;\mathbb{K}). \end{equation*}$$

The purpose of this section is to present a stable element which is slightly simpler than this one. More precisely, the spaces $V_h$ and $Q_h$ are unchanged, but $\Sigma _h$ will be simplified from full linears to matrix fields whose tangential–normal components on each two-dimensional face of a tetrahedron are only a reduced space of linears.

We will still adopt the notation of differential forms. By examining the proof of Theorem 7.1, we realize that we do not use the complete sequence 5.2 for the given spaces. We only use the sequences

$$\begin{equation} \begin{array}{rl} \begin{CD} &{\mathcal{P}}_0^+ \Lambda ^2({\mathcal{T}}_h;\mathbb{K}) @>d>> {\mathcal{P}}_0 \Lambda ^3({\mathcal{T}}_h;\mathbb{K})@>>> 0, \end{CD}\\\begin{CD} {\mathcal{P}}_{1}^+\Lambda ^1 ({\mathcal{T}}_h;{\mathbb{V}})@>d>> &{\mathcal{P}}_{1} \Lambda ^2({\mathcal{T}}_h;{\mathbb{V}}) @>d>> {\mathcal{P}}_0 \Lambda ^3({\mathcal{T}}_h;{\mathbb{V}})@>>> 0. \end{CD} \end{array} {\tag{8.1}} \end{equation}$$

The purpose here is to show that it is possible to choose subspaces of some of the spaces in 8.1 such that the desired properties still hold. More precisely, compared to 8.1, the spaces ${\mathcal{P}}_{1}^+\Lambda ^1 ({\mathcal{T}}_h;{\mathbb{V}})$ and ${\mathcal{P}}_{1} \Lambda ^2({\mathcal{T}}_h;{\mathbb{V}})$ are simplified, while the three others remain unchanged. If we denote by ${\mathcal{P}}_{1,-}^+\Lambda ^1 ({\mathcal{T}}_h;{\mathbb{V}})$ and ${\mathcal{P}}_{1,-} \Lambda ^2({\mathcal{T}}_h;{\mathbb{V}})$ the simplifications of the spaces ${\mathcal{P}}_{1}^+\Lambda ^1 ({\mathcal{T}}_h;{\mathbb{V}})$ and ${\mathcal{P}}_{1} \Lambda ^2({\mathcal{T}}_h;{\mathbb{V}})$, respectively, then the properties we need are that

$$\begin{equation} \begin{CD} {\mathcal{P}}_{1,-}^+\Lambda ^1 ({\mathcal{T}}_h;{\mathbb{V}})@>d>> &{\mathcal{P}}_{1,-} \Lambda ^2({\mathcal{T}}_h;{\mathbb{V}}) @>d>> {\mathcal{P}}_0 \Lambda ^3({\mathcal{T}}_h;{\mathbb{V}})@>>> 0 \end{CD} {\tag{8.2}} \end{equation}$$

is a complex and that the surjectivity assumption 5.5 holds, i.e., $S_h = S_{1,h}$ maps the space ${\mathcal{P}}_{1,-}^+\Lambda ^1 ({\mathcal{T}}_h;{\mathbb{V}})$ onto ${\mathcal{P}}_0^+ \Lambda ^2({\mathcal{T}}_h;\mathbb{K})$. Note that if ${\mathcal{P}}_{0}^{+} \Lambda ^2({\mathcal{T}}_h;{\mathbb{V}}) \subset {\mathcal{P}}_{1,-} \Lambda ^2({\mathcal{T}}_h;{\mathbb{V}})$, then $d$ maps ${\mathcal{P}}_{1,-} \Lambda ^2({\mathcal{T}}_h;{\mathbb{V}})$ onto ${\mathcal{P}}_0 \Lambda ^3({\mathcal{T}}_h;{\mathbb{V}})$.

We first show how to construct ${\mathcal{P}}_{1,-}^+\Lambda ^1 ({\mathcal{T}}_h;{\mathbb{V}})$ as a subspace of ${\mathcal{P}}_{1}^+\Lambda ^1 ({\mathcal{T}}_h;{\mathbb{V}})$. Since the construction is done locally on each tetrahedron, we will show how to construct a space ${\mathcal{P}}_{1,-}^+\Lambda ^1 (T;{\mathbb{V}})$ as a subspace of ${\mathcal{P}}_{1}^+\Lambda ^1 (T;{\mathbb{V}})$. We begin by recalling that the face degrees of freedom of ${\mathcal{P}}_{1}^+\Lambda ^1 (T;{\mathbb{V}})$ have the form

$$\begin{equation*} \int _f \omega \wedge \mu , \qquad \mu \in {\mathcal{P}}_0 \Lambda ^1(f, {\mathbb{V}}). \end{equation*}$$

We then observe that this six-dimensional space can be decomposed into

$$\begin{equation*} {\mathcal{P}}_0 \Lambda ^1(f; {\mathbb{V}}) = {\mathcal{P}}_0 \Lambda ^1(f; T_f) + {\mathcal{P}}_0 \Lambda ^1(f; N_f), \end{equation*}$$

i.e., into forms with values in the tangent space to $f$, $T_f$ or the normal space $N_f$. This is a $4+2$-dimensional decomposition. Furthermore,

$$\begin{equation*} {\mathcal{P}}_0 \Lambda ^1(f; T_f) = {\mathcal{P}}_0 \Lambda _{\operatorname {sym}}^1(f; T_f) + {\mathcal{P}}_0 \Lambda _{\operatorname {skw}}^1(f; T_f), \end{equation*}$$

where $\mu \in {\mathcal{P}}_0 \Lambda ^1(f; T_f)$ is in ${\mathcal{P}}_0 \Lambda _{\operatorname {sym}}^1(f; T_f)$ if and only if $\mu (s)\cdot t = \mu (t)\cdot s$ for orthonormal tangent vectors $s$ and $t$. Finally, we obtain a $3+3$-dimensional decomposition

$$\begin{equation*} {\mathcal{P}}_0 \Lambda ^1(f; {\mathbb{V}}) = {\mathcal{P}}_0 \Lambda _{\operatorname {sym}}^1(f;T_f) + {\mathcal{P}}_0 \Lambda _{\operatorname {skw}}^1(f; {\mathbb{V}}), \end{equation*}$$

where

$$\begin{equation*} {\mathcal{P}}_0 \Lambda _{\operatorname {skw}}^1(f; {\mathbb{V}})= {\mathcal{P}}_0 \Lambda _{\operatorname {skw}}^1(f; T_f) + {\mathcal{P}}_0 \Lambda ^1(f; N_f). \end{equation*}$$

In more explicit terms, if $\mu \in {\mathcal{P}}_0 \Lambda ^1(F; {\mathbb{V}})$ has the form

$$\begin{equation*} \mu (q) = (a_1 t + a_2 s + a_3 n) q \cdot t + (a_4 t + a_5 s + a_6 n) q \cdot s, \end{equation*}$$

where $t$ and $s$ are orthonormal tangent vectors on $f$, $n$ is the unit normal to $f$, and $q$ is a tangent vector on $f$, then we can write $\mu = \mu _1 + \mu _2$, with $\mu _1 \in {\mathcal{P}}_0 \Lambda ^1_{\operatorname {sym}}(f;{\mathbb{V}})$ and $\mu _2 \in {\mathcal{P}}_0 \Lambda ^1_{\operatorname {skw}}(f;{\mathbb{V}})$, where

$$\begin{align*} \mu _1(q) &= \left(a_1t + \frac{a_2+a_4}{2} s\right) q \cdot t + \left(\frac{a_2+a_4}{2} t + a_5 s\right) q \cdot s, \\ \mu _2(q) &= \left(\frac{a_2-a_4}{2} s + a_3 n\right) q \cdot t + \left(\frac{a_4-a_2}{2} t + a_6 n\right) q \cdot s. \end{align*}$$

The reason for this particular decomposition of the degrees of freedom is that if we examine the proof of Lemma 6.1, where equation 6.3 is established, we see that the only degrees of freedom that are used for an element $\omega \in {\mathcal{P}}_{1}^+\Lambda ^1 (T;{\mathbb{V}})$ are the subset of the face degrees of freedom given by

$$\begin{equation*} \int _f \omega \wedge (S_0^{\prime } \nu ), \qquad \nu \in {\mathcal{P}}_0 \Lambda ^0(f; \mathbb{K}). \end{equation*}$$

However, for $\nu \in {\mathcal{P}}_0 \Lambda ^0(f;\mathbb{K})$, $\mu = S_0^{\prime } \nu$ is given by $\mu (q) = \nu q$. Since the general element $\nu \in {\mathcal{P}}_0 \Lambda ^0(\mathbb{K})$ can be written in the form $b_1(t s^T - s t^T) + b_2 (n t^T - t n^T) + b_3 (n s^T - s n^T)$, $\nu q = (- b_1 s + b_2 n) q \cdot t + (b_1 t + b_3 n) q \cdot s$ for $q$ a tangent vector, and thus $\mu \in {\mathcal{P}}_0 \Lambda ^1_{\operatorname {skw}}(f; {\mathbb{V}})$. Hence, we have split the degrees of freedom into three on each face that we need to retain for the proof of Lemma 6.1 and three on each face that are not needed. The reduced space ${\mathcal{P}}_{1,-}^+ \Lambda ^1 (T;{\mathbb{V}})$ that we now construct has two properties. The first is that it still contains the space ${\mathcal{P}}_{1} \Lambda ^1 (T;{\mathbb{V}})$ and the second is that the unused face degrees of freedom are eliminated (by setting them equal to zero). We can achieve these conditions by first writing an element $\omega \in {\mathcal{P}}_{1}^+ \Lambda ^1 (T;{\mathbb{V}})$ as $\omega = \Pi _h \omega + (I- \Pi _h) \omega$, where $\Pi _h$ denotes the usual projection operator into ${\mathcal{P}}_{1} \Lambda ^1 (T;{\mathbb{V}})$ defined by the edge degrees of freedom. Then the elements in $(I- \Pi _h) {\mathcal{P}}_{1}^+ \Lambda ^1 (T;{\mathbb{V}})$ will satisfy

$$\begin{equation*} \int _e \omega \wedge \mu =0, \quad \mu \in {\mathcal{P}}_1 \Lambda ^0(e;{\mathbb{V}}), \quad e \in \Delta _1(T), \end{equation*}$$

i.e., their traces on the edges will be zero. Thus, they are completely defined by the face degrees of freedom

$$\begin{equation*} \int _f \omega \wedge \mu , \qquad \mu \in {\mathcal{P}}_0 \Lambda ^1(f;{\mathbb{V}}), \quad f \in \Delta _2(T). \end{equation*}$$

Since this is the case, we henceforth denote $(I- \Pi _h) {\mathcal{P}}_{1}^+ \Lambda ^1 (T;{\mathbb{V}})$ by ${\mathcal{P}}_{1,f}^+ \Lambda ^1 (T;{\mathbb{V}})$.

We then define our reduced space

$$\begin{equation*} {\mathcal{P}}_{1,-}^+ \Lambda ^1 (T;{\mathbb{V}}) = {\mathcal{P}}_{1} \Lambda ^1 (T;{\mathbb{V}}) + {\mathcal{P}}_{1,f,-}^+ \Lambda ^1 (T;{\mathbb{V}}), \end{equation*}$$

where ${\mathcal{P}}_{1,f,-}^+ \Lambda ^1 (T;{\mathbb{V}})$ denotes the set of forms $\omega \in {\mathcal{P}}_{1,f}^+ \Lambda ^1 (T;{\mathbb{V}})$ satisfying

$$\begin{equation*} \int _f \omega \wedge \mu = 0, \qquad \mu \in {\mathcal{P}}_0 \Lambda ^1_{\operatorname {sym}}(f;{\mathbb{V}}), \end{equation*}$$

i.e., we have set the unused degrees of freedom to be zero.

Then

$$\begin{equation*} {\mathcal{P}}_{1,-}^+\Lambda _h^1 ({\mathcal{T}}_h;{\mathbb{V}})= \{ \omega \in {\mathcal{P}}_{1}^+\Lambda ^1({\mathcal{T}}_h;{\mathbb{V}}) : \, \omega |_T \in {\mathcal{P}}_{1,-}^+\Lambda ^1(T;{\mathbb{V}}), \, \forall T \in {\mathcal{T}}_h\}. \end{equation*}$$

The degrees of freedom for this space are then given by

$$\begin{equation} \int _e \omega \wedge \mu , \quad \mu \in {\mathcal{P}}_1 \Lambda ^0(e;{\mathbb{V}}), e \in \Delta _1(T), \quad \int _f \omega \wedge \mu , \quad \mu \in {\mathcal{P}}_0 \Lambda _{skw}^1(f;{\mathbb{V}}), f \in \Delta _2(T). {\tag{8.3}} \end{equation}$$

It is clear from this definition that the space ${\mathcal{P}}_{1,-}^+\Lambda ^1 (T;{\mathbb{V}})$ will have 48 degrees of freedom (36 edge degrees of freedom and 12 face degrees of freedom). The unisolvency of this space follows immediately from the unisolvency of the spaces ${\mathcal{P}}_{1} \Lambda ^1(T;{\mathbb{V}})$ and ${\mathcal{P}}_{1,f,-}^{+} \Lambda ^1(T;{\mathbb{V}})$.

The motivation for this choice of the space ${\mathcal{P}}_{1,-}^+\Lambda _h^1 ({\mathcal{T}}_h;{\mathbb{V}})$ is that it easily leads to a definition of the space ${\mathcal{P}}_{1,-} \Lambda ^2({\mathcal{T}}_h;{\mathbb{V}})$ that satisfies the property that 8.2 is a complex. We begin by defining

$$\begin{equation*} {\mathcal{P}}_{1,-} \Lambda ^2 (T;{\mathbb{V}}) = {\mathcal{P}}_{0}^+ \Lambda ^2 (T;{\mathbb{V}}) + d {\mathcal{P}}_{1,f,-}^+ \Lambda ^1 (T;{\mathbb{V}}). \end{equation*}$$

It is easy to see that this space will have 24 face degrees of freedom. Note this is a reduction of the space ${\mathcal{P}}_{1} \Lambda ^2 (T;{\mathbb{V}})$, since

$$\begin{equation*} {\mathcal{P}}_{1} \Lambda ^2 (T;{\mathbb{V}}) = {\mathcal{P}}_{0}^+ \Lambda ^2 (T;{\mathbb{V}}) + d {\mathcal{P}}_{1,f}^+ \Lambda ^1 (T;{\mathbb{V}}). \end{equation*}$$

We then define

$$\begin{equation*} {\mathcal{P}}_{1,-} \Lambda ^2 ({\mathcal{T}}_h;{\mathbb{V}})= \{ \omega \in {\mathcal{P}}_{1} \Lambda ^2({\mathcal{T}}_h;{\mathbb{V}}) : \, \omega |_T \in {\mathcal{P}}_{1,-} \Lambda ^2(T;{\mathbb{V}}), \, \forall T \in {\mathcal{T}}_h\}. \end{equation*}$$

It is clear that ${\mathcal{P}}_{0}^+ \Lambda ^2({\mathcal{T}}_h;{\mathbb{V}}) \subset {\mathcal{P}}_{1,-} \Lambda ^2({\mathcal{T}}_h;{\mathbb{V}})$. The fact that the complex 8.2 is exact now follows directly from the fact that the complex

$$\begin{equation} \begin{CD} {\mathcal{P}}_{1} \Lambda ^1 (T;{\mathbb{V}})@>d>> &{\mathcal{P}}_{0}^+ \Lambda ^2(T;{\mathbb{V}}) @>d>> {\mathcal{P}}_0 \Lambda ^3(T;{\mathbb{V}})@>>> 0 \end{CD} \tag{8.4} \end{equation}$$

is exact and the definition

$$\begin{equation*} P_{1,-} \Lambda ^1(T;{\mathbb{V}}) = {\mathcal{P}}_{0}^+ \Lambda ^1(T;{\mathbb{V}}) + {\mathcal{P}}_{1,f,-}^+ \Lambda ^1(T;{\mathbb{V}}). \end{equation*}$$

We will define appropriate degrees of freedom for the space ${\mathcal{P}}_{1,-} \Lambda ^2 (T;{\mathbb{V}})$ by using a subset of the 36 degrees of freedom for ${\mathcal{P}}_{1} \Lambda ^2 (T;{\mathbb{V}})$, i.e., of $\int _f \omega \wedge \mu$, $\mu \in {\mathcal{P}}_1 \Lambda ^0(f;{\mathbb{V}})$. In particular, we take as degrees of freedom for ${\mathcal{P}}_{1,-} \Lambda ^2 (T;{\mathbb{V}})$,

$$\begin{equation*} \int _f \omega \wedge \mu ,\quad \mu \in {\mathcal{P}}_{1,skw} \Lambda ^0(f;{\mathbb{V}}), \qquad \forall f \in \Delta _2(T), \end{equation*}$$

where ${\mathcal{P}}_{1,skw} \Lambda ^0(f;{\mathbb{V}})$ denotes the set of $\mu \in {\mathcal{P}}_1 \Lambda ^0(f;{\mathbb{V}})$ that satisfy $d \mu \in {\mathcal{P}}_0 \Lambda ^1_{skw}(f;{\mathbb{V}})$. It is easy to check that such $\mu$ will have the form

$$\begin{equation} \mu = \mu _0 + \alpha _1 (x \cdot t)n + \alpha _2 (x \cdot s)n + \alpha _3 [(x \cdot t)s - (x \cdot s)t], {\tag{8.5}} \end{equation}$$

where $\mu _0 \in {\mathcal{P}}_{0} \Lambda ^0(f;{\mathbb{V}})$.

Since ${\mathcal{P}}_{1,skw} \Lambda ^0(f;{\mathbb{V}})$ is a six-dimensional space on each face, the above quantities specify 24 degrees of freedom for the space ${\mathcal{P}}_{1,-} \Lambda ^2 (T;{\mathbb{V}})$. To see that these are a unisolvent set of degrees of freedom for ${\mathcal{P}}_{1,-} \Lambda ^2 (T;{\mathbb{V}})$, we let $\omega = \omega _0 + d \omega _1$, where $\omega _0 \in {\mathcal{P}}_{0}^+ \Lambda ^2 (T;{\mathbb{V}})$ and $\omega _1 \in {\mathcal{P}}_{1,f,-} \Lambda ^1(T;{\mathbb{V}})$ and set all degrees of freedom equal to zero. Then for $\mu \in {\mathcal{P}}_{0} \Lambda ^0(f;{\mathbb{V}})$, since

$$\begin{equation*} \int _f (\omega _0 + d \omega _1) \wedge \mu =\int _f \omega _0 \wedge \mu , \end{equation*}$$

we see that $\omega _0 =0$ by the unisolvency of the standard degrees of freedom for ${\mathcal{P}}_{0}^+ \Lambda ^2 (T;{\mathbb{V}})$. In addition, for all $\mu \in {\mathcal{P}}_{1,skw} \Lambda ^0(f;{\mathbb{V}})$ and $\omega _0=0$,

$$\begin{equation*} \int _f \omega \wedge \mu = \int _f d \omega _1 \wedge \mu = \int _f \omega _1 \wedge d \mu . \end{equation*}$$

Since $d \mu \in {\mathcal{P}}_0 \Lambda ^1_{skw}(f;{\mathbb{V}})$, $\omega _1 =0$ by the unisolvency of the degrees of freedom of the space ${\mathcal{P}}_{1,f,-} \Lambda ^1(T;{\mathbb{V}})$.

Using an argument completely parallel to that used previously, it is straightforward to show that the simplified spaces also satisfy assumption 5.5, i.e., that $S_h$ is onto.

To translate the degrees of freedom of the space ${\mathcal{P}}_{1,-}\Lambda ^2(T; {\mathbb{V}})$ to more standard finite element degrees of freedom, we use the identification of an element $\omega \in \Lambda ^2(T;{\mathbb{V}})$ with a matrix $F$ given by $\omega (v_1,v_2) = F(v_1 \times v_2)$. Then $\omega (t,s) = F n$ and $\int _f \omega \wedge \mu = \int _f \mu ^T F n \, df$. Since $\mu \in {\mathcal{P}}_{1,skw} \Lambda ^0(f;{\mathbb{V}})$ and hence is of the form 8.5, we get on each face the six degrees of freedom

$$\begin{equation*} \int _f F n \,df, \quad \int _f (x \cdot t)n^T F n \, df, \quad \int _f (x \cdot s)n^T F n \, df, \quad \int _f [(x \cdot t)s^T - (x \cdot s)t^T] F n \, df. \end{equation*}$$

Finally, we note that the analogue of Theorem 7.2 holds with $r=0$ for the simplified spaces.