An analysis of the TDNNS method using natural norms

The tangential-displacement normal-normal-stress (TDNNS) method is a finite element method for mixed elasticity. As the name suggests, the tangential component of the displacement vector as well as the normal-normal component of the stress are the degrees of freedom of the finite elements. The TDNNS method was shown to converge of optimal order, and to be robust with respect to shear and volume locking. However, the method is slightly nonconforming, and an analysis with respect to the natural norms of the arising spaces was still missing. We present a sound mathematical theory of the infinite dimensional problem using the space \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf {H}}}(\mathbf {curl})$$\end{document}H(curl) for the displacement. We define the space \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\underline{{\mathbf {H}}}}({\text {div}}\,\mathbf {{div}})$$\end{document}H_(divdiv) for the stresses and provide trace operators for the normal-normal stress. Moreover, the finite element problem is shown to be stable with respect to the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf {H}}}(\mathbf {curl})$$\end{document}H(curl) and a discrete \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\underline{{\mathbf {H}}}}({\text {div}}\,\mathbf {{div}})$$\end{document}H_(divdiv) norm. A-priori error estimates of optimal order with respect to these norms are obtained. Beside providing a new analysis for the elasticity equation, the numerical techniques developed in this paper are a foundation for more complex models from structural mechanics such as Reissner Mindlin plate equations, see Pechstein and Schöberl (Numerische Mathematik 137(3):713–740, 2017).


Introduction
In [19], we introduced the TDNNS method for treating the problem of linear elasticity. The TDNNS method is a finite element method that uses tangential-continuous elements for the displacements and symmetric normal-normal continuous finite elements for the stresses. We showed that the TDNNS method is capable of overcoming shear locking [20] and volume locking [22].
However, the TDNNS method is slightly nonconforming, as the stress finite elements are not in the infinite-dimensional distributional space H(div div), which was introduced in [19]. The analysis of TDNNS finite elements provided in our former work [19,20] is based on discrete, broken norms rather than the natural norms of the infinite-dimensional spaces H(curl) and H(div div). In the present paper, we want to provide an analysis based on the natural norms of the Sobolev spaces. This analysis takes the fact that the stress space is nonconforming into account, and leads to optimal order a-priori error estimates.
The necessity of this new framework becomes evident in the analysis of problems in structural mechanics. In [21], the authors use the new technique to analyze the Reissner-Mindlin plate problem. There, we present finite elements based on the TDNNS formulation, which can be seen as an extension of the Hellan-Herrmann-Johnson formulation [9,10,13] to moderately thick plates. In the Reissner-Mindlin plate problem, the unknowns are the deflection and the rotation of the cross-section. As the plate approaches the limit of zero thickness, the rotation tends to the gradient of the deflection. One speaks of the Kirchhoff constraint, which leads to locking if the gradient of the deflection is not contained in the rotation space. However, this inclusion is satisfied (also for the discrete spaces) if one searches for the deflection in H 1 and the cross-section rotation in H (curl). Compared to the TDNNS formulation, the rotation takes the role of the displacement, while the symmetric stress tensor in H (div div) is replaced by the tensor of bending and twisting moments.

Notation
We shortly present the notation used throughout the paper: Vectors shall be denoted as boldface (e.g. u), while tensors are boldface and underlined (e.g. σ ). On the boundary of some domain A, we use the outer normal vector n. For a vector field u, u n = u · n is the normal component, and u t = u−u n n is the tangential component. For a tensor field σ , let σ n = σ n be the normal component, which is further split into its normal-normal component σ nn = (σ n) · n and its normal-tangential component σ nt = σ n − σ nn n.
Gradient, curl and divergence operators ∇, curl and div operators are defined in the usual way. The gradient ∇ of a vector field is a tensor containing in each row the gradient of the corresponding vector component. The divergence div of a tensor field is a vector, where each component is the divergence of the corresponding row of the tensor.
On some domain A, we use the Lebesgue space L 2 (A) and the standard Sobolev space H 1 (A) of weakly differentiable L 2 functions with gradient in L 2 (A). To indicate vector or tensor valued spaces, we use L 2 (A), H 1 (A) and L 2 (A), H 1 (A), respectively.
The space of tensor-valued symmetric functions with components in L 2 (A) is denoted as L 2 sym (A). The space of smooth functions on the closureĀ is denoted as C ∞ (Ā), and C ∞ 0,Γ (Ā) is the subspace where all derivatives vanish on the boundary part Γ of A. If the domain of interest Ω is concerned, it may be omitted, writing e.g. H 1 for On the boundary of a domain we use differential operators and spaces as introduced in the work of Buffa and Ciarlet [5,6]. For the exact definitions, we refer to their work. We will mostly need the rather well-known trace space H 1/2 (∂ A) and the spaces H 1/2 (Γ ), H 1/2 00 (Γ ) on a part Γ of the boundary, where the latter can be extended by zero to the whole boundary space H 1/2 (∂ A).

Problem geometry
Throughout the paper, we assume the domain of interest Ω ⊂ R 3 to be a bounded, connected, polyhedral domain with Lipschitz boundary ∂Ω. Note that all results can directly be transferred to the two-dimensional case.
The (closed) polygonal faces of the polyhedral domain Ω shall be denoted by Γ i with i ∈ I and I a suitable index set. Different boundary conditions will be prescribed different parts of the boundary ∂Ω. To this end, the boundary is divided into two closed parts Γ D and Γ N = ∂Ω\int(Γ D ). The boundary part Γ D , where the displacement will be prescribed (Dirichlet boundary condition), shall be non-trivial, whereas the boundary part Γ N , where surface tractions are given (Neumann boundary condition), may vanish.
We assume that both Γ D and Γ N are aligned with the boundary faces Γ i , such that they each are a union of boundary faces, In accordance with [12] we assume that for each connected component of the Dirichlet boundary Γ D,i we can find an open Lipschitz domain Ω D,i ⊂ R 3 such that Moreover, Ω D,i and Ω D, j have positive distance for i = j, and the interior of Ω ∪ Ω D,i is Lipschitz.

Linear elasticity
Let u : Ω → R 3 be the displacement vector. In linear elasticity, we use the linearized strain tensor We are interested in finding displacement vector u and symmetric stress tensor σ satisfying Hooke's law (4) connects strain and stress tensor by the compliance tensor C −1 , which is the inverse of the standard elasticity tensor C depending on Young's modulus E and the Poisson ratio ν. We assume that Young's modulus E is bounded, and the Poisson ratio ν is bounded away from 1/2, such that both C and C −1 exist and lie in L ∞ (Ω). Equation (5) is the equilibrium condition. We assume that all boundary conditions are prescribed on the boundary parts Γ D and Γ N introduced above. The displacement shall be given on Γ D , while tractions are given on Γ N ,

Motivation of the TDNNS method
Two different variational formulations are widely known for the partial differential equations (4), (5). Most standard finite element methods rely on a primal formulation, where the stress tensor σ is eliminated. In this formulation, the displacement boundary condition on Γ D is essential, and usually treated by a homogenization approach. To this end, it is necessary to have the existence of an extensionũ D ∈ H 1 (Ω) of the boundary displacement u D to the whole domain Ω. Then one searches for u ∈ũ To be conforming, the displacement finite element space has to be continuous. On the other hand, a dual Hellinger-Reissner formulation can be obtained from system (4), (5). Integration by parts puts all continuity assumptions to the stress tensor. It has to allow for a weak divergence, while only L 2 regularity is required for u. In this case, the traction boundary conditions are essential. One needs an extension of the surface tractions t N to the domain, a tensor fieldσ N ∈ H sym (div) withσ N ,n = t N on Γ N . Inhomogeneous displacement boundary conditions can be included in weak form into the right hand side of equation (9).
One searches for σ ∈σ N + H sym,0,Γ N (div) with the space H sym 0,Γ N (div) = {σ ∈ L 2 sym (Ω) : divσ ∈ L 2 (Ω), σ n = 0 on Γ N } satisfying the homogeneous traction boundary condition and u ∈ L 2 (Ω) such that To define a conforming finite element method, one has to provide stress elements which are symmetric and normal continuous. Such elements have been found [1][2][3], but come only at high computational costs, as they involve at least 24 degrees of freedom per element in two dimensions or 162 in three dimensions. The TDNNS formulation is in between the primal and the dual concept. We want to design a formulation, where the tangential component u t of the displacement and the normal component σ nn of the normal stress vector are essential boundary conditions on the respective boundary parts Γ D and Γ N . In other words, the displacement space has to allow for the definition of a tangential trace, while the stress space allows a normalnormal trace. It will turn out that the displacement space is the space H 0,Γ D (curl) satisfying zero tangential boundary conditions on Γ D .
Below, we formally write the variational formulation. It is of the standard mixed form treated in [4]. We use the stress space and displacement space V, which will be rigorously defined in Sect. 2. Currently, we only state that v ∈ V implies v t = 0 on Γ D and τ ∈ implies τ nn = 0 on Γ N . Accordingly, we need two extensions, one for the tangential component of the displacement and one for the normal-normal component of the stress: on the Dirichlet boundary Γ D someũ D withũ D,t = u D,t and on the Neumann boundary Γ N someσ N withσ N ,nn = t N ,n . We want to find u iñ u D + V and σ ∈σ N + such that For smooth functions, the bilinear forms a(·, ·) and b(·, ·) are given by = In Sect. 2 we will define the function spaces and give a precise meaning to the arising integrals and bilinear forms in a distributional setting. We will determine in which way the boundary terms have to be understood. We will see that the bilinear form b(·, ·) corresponds to a distributional divergence operator.

The variational formulation of the TDNNS method
We shall specify the spaces, in which the variational formulation (11)-(12) is posed. While the displacement space is well-known, the stress space was introduced in [19] and shall be analysed in detail in this work.

The displacement space H(curl)
We use the space This is a Hilbert space equipped with inner product and induced norm It is well known, that the space H(curl) allows for the definition of a tangential trace. According to [6], we may define the subspace of H(curl) satisfying homogeneous tangential boundary conditions on Γ D , which is our TDNNS displacement space The following theorem is taken, with notation adapted to our work, from [6, Theorem 6.6, Remark 6.7]:

Theorem 1
The tangential trace operator v → v t is bounded and surjective as a mapping and The first statement of Theorem 1 ensures the existence of an extensionũ D of a given tangential-displacement boundary value. Additionally, Theorem 1 tells that the surface integral in (12) can be understood as a duality product. We will elaborate on this matter in Sect. 3.
A conforming finite element space for H(curl) has to be tangential continuous, such as the Nédélec spaces introduced in [16,17].
An essential tool in the analysis of the TDNNS formulation is the regular decomposition. Decompositions satisfying homogeneous Dirichlet or Neumann boundary conditions have been shown by [18] and [11], respectively. The case of mixed boundary conditions can be found in [12].
The respective parts can be bounded by with a generic constant c.

The stress space H(div div)
We still need to specify the stress space. Roughly, it is a subspace of L 2 where the (scalar-valued) divergence of the (vector-valued) divergence of the stress tensor lies in the dual space of H 1 0,Γ D . In (23), the norm of the desired space is stated for smooth functions. We will proceed as follows: first, we formally define the space H(div div) as the closure of smooth functions, and give an interpretation of the norm which rectifies the name H(div div). Then, we show that the normal-normal trace can be bounded in this norm in the appropriate setting. Thus we can define the subspace H 0,Γ N (div div) satisfying zero normal-normal boundary conditions on Γ N as the closure of smooth functions vanishing on Γ N . Last, we provide an inverse trace theorem, which allows to extend normal-normal stress distributions from the boundary to the whole space H(div div).
The norm · H(div div) shall be defined for smooth τ ∈ C ∞ (Ω) by Note that, due to the symmetry of the Hessian ∇ 2 ϕ, the symmetric expression ε(∇ϕ) is the same as the conventional notation ∇ 2 ϕ. We use ε(∇ϕ), as it shows the relation to linear elasticity.
We define the space H(div div) as The second term in the definition of the norm (23) is a seminorm and can be interpreted as the norm of div div τ in the dual space of H 1 0,Γ D . Integration by parts of the denominator gives for smooth τ , using that ϕ and ∂ϕ/∂t vanish on Γ D , The supremum in (23) can be interpreted as a dual norm: in the interior, div divτ is in the dual space of In the last term, the tangential derivative ∂ϕ/∂t appears. Since the gradient of H 1 lies in H(curl), this tangential derivative is in H The normal-tangential stress τ nt is thus in the dual of this space, which means [6] We will comment on this restriction in Sect. 4, as it our finite element space is not conforming in this term.
We will now define a space for the normal-normal trace, and show that the normalnormal trace is bounded in the H(div div) norm. To this end, we need the space of traces of the normal derivative of For a polyhedral domain and Γ D = ∅, H 1/2 n (∂Ω) consists of piecewise H 1/2 (Γ i ) without continuity assumptions on the polyhedron edges or vertices, see e.g. [7]. For Γ D = ∂Ω, H 1/2 n (∂Ω) is the subspace of the piecewise H 1/2 00 (Γ i ) spaces without continuity assumptions on the polyhedron edges or vertices, see e.g. [8]. To the best knowledge of the authors, this space has not been analyzed so far for general, nontrivial Γ D . In this work, we only use that the space can be defined piecewise on each polyhedral face.
The normal-normal trace space of H(div div) is then given by An appropriate norm on H −1/2 n (∂Ω) is given by Note that due to the piecewise nature of H 1/2 n (∂Ω), the trace space can be restricted to each polyhedral face Γ i , and extended from each face to the whole boundary by zero.

Theorem 3 The normal-normal trace operator is bounded from H(div div) to H
Thus, it is well defined on H(div div) as the extension from C ∞ sym . For τ ∈ H(div div) there holds the bound with the constant c independent of τ .
Proof Let τ ∈ C ∞ sym be fixed. We first show that the normal-normal trace on a boundary face Γ i can be bounded by the normal-normal trace on the whole boundary. Since H 1/2 n (∂Ω) is a piecewise defined space without continuity assumptions between polyhedron faces, any ϕ ∈ H 1/2 n (Γ i ) can be extended by zero to φ ∈ H 1/2 n (∂Ω). By definition of the dual norm we see We proceed showing the actual trace inequality, where we use that H We see that the supremum from Eq. (38) is already contained in the H(div div)-norm. For the supremum from Eq. (39), we use that Therefore, we arrive at the desired result The trace theorem above allows to define the space . Finally, we provide an inverse trace theorem for the space H(div div), before we proceed to the analysis of the TDNNS elasticity problem. The inverse trace theorem allows to find an extension of a given (scalar) normal-normal stress on the boundary to a (tensor-valued) stress field on the domain.

Theorem 4 For any boundary face
with constant c independent of g.
Proof For boundary face Γ i , let g ∈ H −1/2 n (Γ i ) be given. The extension of g by zero lies in H −1/2 n (∂Ω), which is the dual of the trace space of H 2 ∩ H 1 0,Γ D . This allows to pose the following problem in Solvability of (46) is clear, as we note that ε(∇w) = ∇ 2 w due to the symmetry of the Hessian. By the standard theory of Lax and Milgram we obtain the stability estimate We choose τ = ε(∇w), which is clearly symmetric and bounded in L 2 : Additionally, it satisfies the natural boundary condition τ nn = g in H −1/2 n (Γ i ). It remains to show that our choice of τ lies actually in H(div div), and satisfies the estimate (45). To this end, we still need to bound the supremum term in the definition of the norm. Any ϕ ∈ C ∞ ∩ H 1 0,Γ D is a valid test function for the variational equation (46). This implies Adding up (48) and (49)-(52) leads to the desired bound (45).
With tools concerning H(div div) now at hand, we can proceed to the analysis of the variational problem (11)-(12).

Analysis of the TDNNS problem
In the current section, we show existence and uniqueness of a solution to the TDNNS elasticity problem (11)- (12). We specify the variational spaces and V foreshadowed in Sect. 2, We shall use the theory on mixed problem treated in detail in [4]. First, we concentrate on boundary conditions, then we show stability estimates for the bilinear forms a(·, ·) and b(·, ·). These results allow us to derive existence, uniqueness and stability of a solution to the TDNNS elasticity problem.

Boundary conditions
We assumed boundary conditions u D on Γ D and t N on Γ N to be given. We shall comment on the regularity necessary for these boundary conditions, such that the variational problem is well-defined. We treat the essential boundary conditions on tangential displacement and normal-normal stress first, and proceed to the natural boundary conditions on the normal displacement and normal-tangential stress afterwards.
In the variational formulation (11)-(12), we used extensionsũ D andσ N of the given boundary data. The trace theorems for H(curl) and H(div div) ensure that, given u D,t ∈ H −1/2 ⊥,00 (curl, Γ D ) and t N ,n ∈ H −1/2 n (Γ N ), extensions can be found satisfying Natural boundary conditions on normal displacement and tangential component of normal stress are included into the right hand side of the variational problem (11)- (12). For smooth functions, they are included as surface integrals Γ D u D,n τ nn ds and We will see that both boundary integrals can be understood in the sense of duality products in the respective trace spaces, which makes them well-defined on the whole variational spaces. The trace theorem on H(div div), Theorem 3, ensures that Thus, it is necessary to have the normal displacement u D,n ∈ H 1/2 n (Γ D ). The second statement from Theorem 1, (20) ensures that for v ∈ V the tangential trace v t allows for a surface curl, v t ∈ H −1/2 ⊥ (curl 0 , Γ N ). Therefore, the normaltangential trace of the given surface tractions has to lie in its dual space, which is by [6] Let us shortly comment on the condition τ nt = t N ,t ∈ H −1/2 ,00 (div, Γ N ). This means that the given normal-tangential (shear) stress has to allow for a distributional surface divergence of some kind, which includes a continuity assumption on the in-plane normal across boundary edges. As the normal-tangential (shear) component of the proposed stress finite elements does not satisfy this condition, the finite element method is nonconforming, see Sect. 4. Similarly, also the given shear stress t N ,t does not need to satisfy any continuity assumptions in the finite element setting.

Stability of bilinear forms
To apply the theory on mixed systems by [4], we have to show -boundedness of a(·, ·) and b(·, ·), -coercivity of a(·, ·) on the kernel space Ker(B), and -inf-sup stability of b(·, ·).

There, the kernel space Ker(B) is defined by
The conditions on a(·, ·) will follow rather quickly, assuming the elasticity matrix C −1 to be regular, which is true for compressible materials with ν < ν 0 < 1/2. In the case of nearly incompressible materials with ν → 1/2, a refined analysis comparable to [22,Chapter 5] has to be carried out, which is not done in the scope of the present work. The estimates on b(·, ·) are more involved.

Lemma 1
The bilinear form a(·, ·) is bounded on = H 0,Γ N (div div). Moreover, it is coercive on the kernel space Ker(B) from (60), there exists a constant c a independent of τ such that Proof Boundedness of a(·, ·) in H(div div) is clear, since a(·, ·) is a (C −1 -scaled) inner product on L 2 and H(div div) is a subspace of L 2 . Obviously, a(·, ·) is also coercive with respect to the L 2 norm, To get coercivity with respect to the H(div div) norm, we need to show that the additional supremum term in (23) vanishes for τ ∈ Ker(B). Since for ϕ ∈ H 1 0,Γ D we have ∇ϕ ∈ V = H 0,Γ D (curl), we directly see from the definition of the kernel space (60), This concludes the proof.
Next, we treat boundedness of the bilinear form b(·, ·).

Lemma 2 The bilinear form
Proof Let τ ∈ = H 0,Γ N (div div) and v ∈ V = H 0,Γ D (curl) be smooth, such that the integrals in (14) are well defined. We use the regular decomposition (22) Cauchy's inequality in L 2 and density ensure We use the bound z H 1 + ∇ p L 2 ≤ c v H(curl) , and arrive at This continuity result allows us to extend the bilinear form from smooth functions to the whole of H 0,Γ N (div div) × H 0,Γ D (curl) in the sense of a distributional divergence operator.

Lemma 3 The bilinear form
The constant c b > 0 is independent of v.
Proof Let v ∈ H 0,Γ D (curl) be fixed. Find u ∈ H 1 ∩ H 0,Γ D (curl) as a solution to the primal elasticity problem that for allṽ ∈ This solution satisfies the following "classical" boundary conditions: On Γ N , we have a free boundary, with (ε(u)) n = 0, on Γ D , the tangential displacement u t = 0 is fixed, while the normal displacement satisfies u n = −(ε(u)) nn . In the variational setting, we have the combined boundary condition where the surface integral is to be understood as a duality product. We choose σ := Cε(u). We have to show that σ lies in H 0,Γ N (div div). To this end, we first prove that it satisfies the H(div div)-essential boundary condition σ nn = 0 in We now want to bound σ H(div div) . Standard theory for the primal problem (72) ensures We bound the supremum term in the H(div div) norm to show that σ ∈ H(div div): Again, ∇ϕ is a valid test function for the variational problem (72), Together with (75) this leads to the bound

Finite element spaces and their norms
Let T = {T } be a simplicial, (shape-)regular triangulation of Ω as defined in [15,Def. 5.11]. We denote the set of element faces F = {F}. Any piecewise smooth vector field v ∈ H(curl) has to be tangential continuous on element interfaces, i.e. v t is uniquely defined on each element interface. In [19,22], it was shown that a piecewise smooth tensor τ ∈ H(div div) is normal-normal continuous across element faces, i.e. the normal-normal component τ nn is continuous. So far, the bilinear form b(·, ·) is defined only for smooth vector and tensor fields in (14), (15). This definition can be extended to piecewise smooth fields: Let v ∈ H 0,Γ D (curl) and τ ∈ H 0,Γ N (div div) be piecewise smooth and tangential and normalnormal continuous on the triangularization T , respectively, then In [19,22], normal-normal continuous symmetric stress finite elements were constructed, which were used together with Nédélec elements for the displacement. A stable finite element method was obtained. However, the method is slightly nonconforming: to be in H(div div), a piecewise continuous function has to be normal-normal continuous, and the normal-tangential component τ nt has to lie in the dual space of the trace space of H(curl, T ), which means τ nt ∈ H −1/2 || (div ∂ T ). Then the surface integral in (84) can be understood as duality product and evaluated for arbitrary v ∈ H 0,Γ D (curl). However, this is a continuity restriction on τ nt at element edges, which does not hold for general τ piecewise smooth. In general, the surface vector field τ nt is discontinuous across element edges.
We choose the following finite element spaces for integer k ≥ 1 Additionally, we need an H 1 -conforming scalar finite element space W h of order k +1 satisfying zero boundary conditions, For this choice we have that ∇W h ⊂ V h . Utilizing the finite element spaces above, we may pose the finite element problem (restricting ourselves to the case of trivial essential boundary conditions u t = 0 on Γ D and σ nn = 0 on Γ N ): find σ h ∈ h and u h ∈ V h such that

Discrete stress norm
While the Nédélec space V h and the continuous space W h are endowed with H(curl) and H 1 norms, respectively, we provide a discrete norm for the stress space h , Note that, for finite element functions τ h , the face terms in (91) can be bounded by the L 2 term (see [19]). Thus, the face terms may be omitted for finite element functions, which is often done in this work. However, this is not possible for general piecewise smooth normal-normal continuous tensor fields τ .
In [19], we showed stability of the problem not using the H(curl) and discrete H(div div) norm (91), but using the L 2 norm for the stresses and a broken H 1 norm for the displacements,

The reference element and transformations to the mesh element
We introduce the reference tetrahedronT = {x = (x 1 ,x 2 ,x 3 ) :x i > 0,x 1 +x 2 +x 3 < 1}. Barycentric coordinates onT are given bŷ For any element T ∈ T , let be a smooth one-to-one mapping of the reference tetrahedron to tetrahedron T . The Jacobian of this transformation shall be denoted by F T = ∇Φ T , the Jacobi determinant by J T = det (F T ). The local mesh size is defined as the spectral norm of F T , h T = |F T | s . For a face F and an edge E, let J F , J E denote the transformation of measures of the mappingsF → F,Ê → E. For the normal n F to face F and the tangential vector t E to some edge E we have The finite element basis functions are defined on the reference tetrahedron, and mapped to an element T by a conforming transformation. A conforming transformation has to preserve the degrees of freedom of the finite element. While H 1 conforming elements can be transformed directly, we need the tangential-trace preserving covariant transformation for H(curl) conforming elements, and a transformation which preserves the normal-normal trace for the stress elements, By application of basic calculus one can see that gradient and strain operators transform as

Interpolation operators
The a-priori error analysis of the proposed finite element method relies on interpolation operators for the finite element spaces. Subsequently, we recall the nodal interpolation operators for the spaces W h and V h , as well as the Clément quasi-interpolation operator for the piecewise linear, continuous finite element space. Their definitions and according estimates can be found in [15,Sect. 5.5 and 5.6]. Additionally, we present an error estimate for the H(curl)-interpolant in the broken H 1 norm. In Sect. 5.2 we define an interpolation operator for the stress space and give error estimates in the discrete stress norm · h .

Commuting interpolation operators for H 1 and H(curl) and the Clément quasi-interpolation operator
Let I W , I V be the nodal interpolation operators defined on the finite element spaces W h , V h . These interpolation operators are based on the degrees of freedom of the finite element spaces. They are defined in such a way that they commute with the gradient operator, Note that the interpolation operators are not well-defined for general functions in The tensor fields are given bŷ The interpolation operator I mapping any sufficiently smooth, normal-normalcontinuous tensor field τ to I τ ∈ h , is uniquely defined by the following conditions, (113)

Lemma 4 The interpolation operator
Since τ h,nn is polynomial of order k on each face, this implies that τ h,nn = 0 on each face.
Since the normal-normal component of τ h vanishes on all element interfaces, on each element T ∈ T , τ h is a linear combination of element-local interior shape functions. According to [19], there are two types of interior shape functions on the reference element, p iŜT ,n ,p i basis for P k (T ), n = 1, 2.
Transforming the integrals (112), (113) to the reference element using the H(div div) conforming transformation (98) leads to Proof In [22], it was shown that a very similar interpolation operator approximates in the L 2 /L 2 (F) norm. The same estimates holds for I , which is expected, since the local space is the full polynomial space of order k and I preserves piecewise polynomial finite element functions. The proof relies on a scaling argument and the Bramble-Hilbert lemma applied on the reference element, and is not provided in detail here, To estimate the full norm τ − I τ h , we show that the supremum term vanishes, We observe, due to the definition of the interpolation operator I , The stress interpolation operator I is defined in such a way that b(τ − I τ , ∇w h ) vanishes for any w h ∈ W h . Thus, the interpolation error estimate in the natural norm coincides with the estimate in L 2 norm.

Analysis of the finite element problem
A crucial tool for the analysis of the finite element problem is a discrete version of the regular decomposition from Theorem 2. The following discrete decomposition can be deduced directly from the regular decomposition, see [14] for the case of Γ = Γ D .

Lemma 5 For a finite element vector field
with z ∈ H 1 0,Γ D , curlz = curlv h and p h ∈ W h . The respective parts can be bounded by with a generic constant c.

Continuity of the finite element problem
We are concerned with continuity of the bilinear forms with respect to the discrete norm · h and the H(curl) norm · H(curl) . Obviously, a(·, ·) is continuous, as it is continuous in L 2 . For b(·, ·), showing continuity is more challenging.
Proof Let τ , τ | T ∈ H 1 sym (T ) for all T ∈ T , τ nn | F ∈ L 2 (F) be a normal-normal continuous piecewise smooth tensor field. Note that this includes all finite element tensor fields τ h . For v h ∈ V h , let v h = I V z + ∇ p h be the decomposition from Lemma 5. We have We estimate the two parts separately. For the estimate concerning z, we first use that b(·, ·) is continuous in the L 2 /broken H 1 setting.
Next, we utilize the Clément interpolation operator C, which is continuous in By an inverse inequality for the finite element function I V z − Cz we see Using interpolation error estimates for I V and C, we arrive at The estimate concerning ∇ p h follows directly, Together, (137) and (139) lead to the desired continuity result.

Stability of the finite element problem
According to [4], we need to provide stability of the bilinear forms with respect to the discrete norms, i.e. we need to show discrete kernel-coercivity of a(·, ·) and an inf-sup condition for b(·, ·).

Lemma 7
The bilinear form a(·, ·) is coercive on the discrete kernel with generic constant c independent of the mesh size.
Proof Let τ h ∈ Ker(B h ) be fixed. Since τ h is a finite element function, the facet term in the discrete norm (91) can be omitted, Since ∇W h ⊂ V h , the supremum term in (141) vanishes, Hence, coercivity is implied by coercivity of the compliance tensor,

Lemma 8 The bilinear form b(·, ·) is inf-sup stable on h
Proof The proof is very similar to the one of Lemma 3 in the infinite dimensional setting. Nevertheless, we provide it in detail here.
Let v h ∈ V h be fixed. According to the finite element theory using L 2 and broken H 1 norms, there exists a unique pair (τ h , u h ) ∈ h × V h satisfying Moreover, we have the stability estimate Since ṽ h H 1 ,h ≥ c ṽ h H(curl) , we deduce Testing the second Eq. (147) with a gradient function ∇w h , we see using that Hence, we deduce Adding up squared (149) and (151)

Error estimates
Since the finite element method is slightly nonconforming, h ⊂ , the error cannot be bounded directly using the theory from [4]. Instead, we rely on techniques from Strang's second lemma, where consistency and interpolation error bound the approximation error.
Proof We divide the approximation error into two parts, the interpolation error and a consistency term. To this end, we add and subtract the interpolants I σ , I V u and use the triangle inequality, ≤ ( σ − I σ h + u − I V u H(curl) ) (157) We refer to the terms in (157) as interpolation error, while the terms in (158) are referred to as consistency error. We first elaborate on the consistency error, which can itself be bounded by the interpolation error: Due to the discrete stability, we have We proceed for the first term, adding and subtracting Ω C −1 σ : τ h dx and using that the solution u is sufficiently smooth to have ε(u) = C −1 σ in L 2 ,