Stable Mixed Finite Elements for Linear Elasticity with Thin Inclusions

We consider mechanics of composite materials in which thin inclusions are modeled by lower-dimensional manifolds. By successively applying the dimensional reduction to junctions and intersections within the material, a geometry of hierarchically connected manifolds is formed which we refer to as mixed-dimensional. The governing equations with respect to linear elasticity are then defined on this mixed-dimensional geometry. The resulting system of partial differential equations is also referred to as mixed-dimensional, since functions defined on domains of multiple dimensionalities are considered in a fully coupled manner. With the use of a semi-discrete differential operator, we obtain the variational formulation of this system in terms of both displacements and stresses. The system is then analyzed and shown to be well-posed with respect to appropriately weighted norms. Numerical discretization schemes are proposed using well-known mixed finite elements in all dimensions. The schemes conserve linear momentum locally while relaxing the symmetry condition on the stress tensor. Stability and convergence are shown using a priori error estimates.


Introduction
Thin inclusions in elastic materials arise in a variety of scientific fields, including geo-physics, bio-mechanics, and the study of composite materials. The subsurface, for example, typically includes rock layers with significantly larger horizontal extent compared to their height. Since it is often infeasible to resolve such small heights for large-scale simulations, we consider the setting where the layer, or aquifer, is represented by a lower-dimensional manifold. The governing equations on this manifold can be derived using vertical integration [18], see e.g. [6] for an application with respect to CO 2 storage.
Secondly, membranes occur frequently in the study of bio-mechanical systems. Examples range from cell walls in plants to the heart sac and dermal layer in human physiology. As a modeling assumption, each of these membranes can be represented by lower-dimensional manifolds. Their influence on the coupled mechanical system can then be incorporated by assigning significantly different material properties compared to the surroundings.
A third application concerns the study of composite or reinforced materials. In this context, the lower-dimensional manifolds correspond to the stiffer plates embedded in the material for strengthening purposes. This can be expanded to connected, two-dimensional objects such as H-beams and T-beams. The junctions are then considered one-dimensional manifolds, with inherited or separately defined material properties. We note that this work is limited to manifolds of codimension one and thus does not treat the case of embedded, one-dimensional rods in three dimensions.
The thin features are considered lower-dimensional and have elastic properties, yet a slightly different setting is presented than the conventional theory of thin shells [12]. The main difference is that we focus on a strong coupling of a thin inclusion with a surrounding, elastic medium. The interest of this work is therefore more closely aligned to [11], in which rigid, thin inclusions are considered.
The structure of the derived equations fits well with the mixed-dimensional framework derived in [17,8]. We aim to preserve this structure and retain a local conservation of linear momentum after discretization with the use of conforming, mixed finite elements. The construction of stable finite element pairs representing displacements and symmetric stresses is involved, and typically leads to higher-order elements for the stress space [4]. By relaxing the symmetry condition on the stress tensor as in [2,5], these difficulties can be mitigated.
This article is structured as follows. Section 2 introduces the notational conventions and the decomposition of the geometry according to dimension. On this geometry, Section 3 introduces the governing equations of the model in all dimensions. We introduce the relevant function spaces and present the derivation of the variational formulation in Section 4. The resulting system of equations is proven to be well-posed in Section 5. Finally, we propose conforming discretization schemes in Section 6 for which stability and convergence are shown.

Geometry and Notation
In this section, we introduce the mixed-dimensional geometry and establish notation. Here, we follow the conventions introduced in [9,8].
Let us consider an n-dimensional domain Y that contains thin, embedded structures, represented by lower-dimensional manifolds. In general, we consider n = 3, the two-dimensional case being simpler. Let Ω di i be such a manifold, with i the unique index from a global set I and d i its dimension. The superscript is frequently omitted for brevity. For 1 ≤ d ≤ n − 1, we successively identify intersections between d-manifolds as (d−1)-manifolds. All Ω i are open sets and mutually disjoint. For simplicity, we restrict this work to the case in which all Ω i have zero curvature, i.e. are flat.
As an example, let us consider the two-dimensional set-up in Figure 1 (left). Here, an embedded H-beam is described using two zero-dimensional intersection points and five one-dimensional line segments. The open set corresponding to the surrounding medium is given by . We refer to the union i∈I Ω i as the mixed-dimensional geometry Ω. Let I d be the set of indices corresponding to d-manifolds and let Ω d be the collection of such manifolds. In short, we denote The interface between manifolds of codimension one will play an important role, and we adopt a separate notation for these. Let J be the set of indices such that each j ∈ J corresponds to an interface Γ j between Ω i for some i ∈ I and an adjacent domain of dimension (d i + 1). Γ j physically coincides with Ω i and we assume that a unique ∈ I exists such that Γ j ⊆ ∂Ω.
To distinguish these different interfaces, we define the following index sets for i ∈ I:Ĵ An example is shown in Figure 1 (right), which emphasizes that for j 1 , j 2 ∈Ĵ 4 with j 1 = j 2 , we have 1 = 2 = 8. In other words, we allow for a manifold Ω to border on multiple sides of Ω i and assign a unique index j ∈Ĵ i to each side. Finally, we remark thatĴ i is void for all i ∈ I n , by definition.
Using the same summation convention per dimension as above, we denote Each Γ j is equipped with a unit normal vector n j , from the tangent space of Ω, oriented outward with respect to Ω. The subscript on n is omitted for brevity. In reference to the vector(s) normal to a manifold Ω i , we will use the check notationň i . The boundary of the domain is given by the disjoint union ∂ σ Y ∂ u Y on which different boundary conditions will be imposed. In particular, we assume that the displacement is given on ∂ u Y and the normal stress on ∂ σ Y . We denote for i ∈ I, For analysis purposes, we assume that |∂ u Ω i | > 0 for all i ∈ I n , i.e. each subdomain of dimension n is connected to a part of the boundary on which the displacement is prescribed. By omission of the subscript, we use ∂ u Ω and ∂ σ Ω to refer to the corresponding boundaries of the mixed-dimensional geometry.
Given a function f defined on the mixed-dimensional geometry Ω, let f i denote its restriction to Ω i , i.e. f i = (f )| Ωi . Furthermore, we employ the hat and check notation to distinguish instances of f inherited from different domains onto the interface Γ . This means that on Γ j with j ∈Ĵ i , we denotě Note that the definition off involves a trace of f onto Γ j .

Model Formulation
In this section, we consider the governing equations and introduce the model problem. Starting with the mathematical formulation of linear elasticity in the surrounding medium, we continue with the generalized equations on lowerdimensional manifolds to derive the strong form of the mechanics problem. The variational formulation is considered afterward in Section 4.

Governing Equations in the Surrounding Medium
Let us start by presenting the governing equations for linear elasticity in the surrounding medium Ω n . For i ∈ I n , let σ i denote the elastic stress and u i the displacement. Assuming infinitesimal strain, the stress-strain relationship has the general form: In case of homogeneous and isotropic media, the operator A describes Hooke's law and is given by Tr σ i I , (3.1) in which λ and µ are the Lamé parameters and Tr is the matrix trace operator.
In the variational formulation presented Section 4, the symmetry of the stress tensor σ i will be enforced in a weak sense. In preparation, we introduce the antisymmetric tensor χ i = asym ∇u i such that: With the addition of linear and angular momentum conservation, the following system of equations is formed in each Ω i with i ∈ I n .
With f i the body forces acting on Ω i . Note that the balance of angular momentum is enforced as the symmetry of the stress tensor σ i in (3.2c). The associated boundary conditions are given by with g u a given function. We limit the exposition to homogeneous stress boundary conditions, noting that this can readily be extended to the general case.

Geometrical Scaling and Constraints
The governing equations on the lower-dimensional manifolds will be significantly influenced by the small width of the thin inclusions. We therefore devote this section to defining two key parameters, γ and and the constraints these parameters adhere to. Let γ be a virtual parameter representing the relatively small length from the interface of the higher-dimensional domain to the central plane (2D), line (1D), or point (0D) of the physical inclusion. We assume that on each interface Γ j , γ is constant and positive.
The small value of γ will introduce a scaling in the equations. In the formulation of the problem, it is advantageous if this scaling appears in the coupling terms between variables, rather than on the main diagonal of the system. For this purpose, we introduce a key concept from the context of subsurface flow models [1,9]. As shown in those works, the desired scaling of the system is achieved by employing an appropriately scaled flux variable. In analogy, we employ a scaled stress by introducing the scaling parameter on Ω. This parameter is defined as the square root of the (small) cross-sectional length (2D), area (1D), or volume (0D) of the dimensionally reducible feature. is assumed to be constant and positive for each manifold Ω i . The following relationship is imposed with V i the cross-sectional measure of the feature corresponding to Ω i and i = 1 for i ∈ I n . The relation a b (respectively and ) implies that a constant C > 0 exists, independent of , γ, and the mesh size h, such that a = Cb (respectively ≤ and ≥). Since is assumed to be small, we often use the relationship 1. Next, we relate the different values of between the dimensions. For each Ω i , letˆ max be the maximal value of in the adjacent higher-dimensional manifolds: In Ω n , we defineˆ max as unity. Using this definition, we add a constraint to the geometry by assuming thatˆ max bounds from above, i.e.
With the parameter defined, we continue with the scaling of the stress variable on the lower-dimensional manifolds. In Ω i with d i < n, let σ avg i contain the columns of the Cauchy stress tensor associated with the tangent bundle of Ω i , averaged over the cross-section of the physical feature. The integrated stress tensor σ int i , on the other hand, is obtained after multiplication with the cross-sectional measure V i .
Using the factor i := √ V i , the scaled stress σ i on Ω i is defined as The columns and first d i rows of σ i correspond to the basis vectors from the tangent bundle of Ω i . The final (n − d i ) rows relate to the directions normal to the manifold. Thus, σ i ∈ R n×di by definition, making it undefined in the intersection points Ω 0 . The displacement u i in Ω i remains unscaled and is naturally in R n . Again, we omit the subscript i to refer to the mixeddimensional entities, i.e.
We emphasize that the integrated and average stress quantities can readily be recovered from the scaled stress σ after appropriate post-processing with the known quantity .

Mixed-Dimensional Equations
With the given scaling from the previous subsection, let us consider the governing equations in the lower-dimensional manifolds. In this generalization to the mixed-dimensional geometry, a structure similar to the system (3.2) is uncovered. We start by introducing the linear momentum balance equation, followed by the stress-strain relationships and finish with the conservation of angular momentum. The balance of linear momentum (3.2b) is generalized first. After integrating the conservation law in the direction(s) normal to the inclusion (see e.g. [9,15] for the analogue in fracture flow models), we obtain Here, f i is the body force acting on Ω i , averaged over the cross-section with measure 2 . The in-plane divergence (∇·) on Ω i and the normal trace operator (n·) onto Γ j are applied row-wise. Hence, the divergence in Ω i maps from R n×di to R n and the normal trace on Γ j maps R n×d to R n . For the zero-dimensional manifolds, there are no divergence operator or σ i available and the balance law is completely given by the sum of forces from To shorten notation, we introduce the jump operator · i which maps functions defined on the interface Γ j with j ∈Ĵ i to the central manifold Ω i such that ∀i ∈ I.
Following [9], we introduce the mixed-dimensional divergence operator (D·) as D · σ := ∇ · σ − n ·σ and rewrite the conservation equation to the concise form We now continue by defining the stress-strain relationships in the lowerdimensional manifolds in analogy with (3.2a). For that, we first introduce the gradient operator D as We emphasize that the gradient ∇ relates to the tangential direction(s) and is applied row-wise. Since we have Du defined on both Ω and Γ , we need to provide stress-strain relationships inside and on the boundaries of the domains.
The stress-strain relationship is then described by an operator A acting on the averaged stress. Here, we pay special attention to the scaling with . Thus, recalling that the averaged stress is denoted by σ avg = −1 σ, the stress-strain relationships are given by on Ω × Γ with χ to be defined. To obtain a symmetric system, we scale this equation with . Noting that and A do not necessarily commute, we introduce A := A −1 to obtain the generalized version of the stress-strain relationship: The restrictions of A to the manifolds Ω and interfaces Γ are respectively denoted by The variable χ is the generalization of the asymmetric χ i from section 3.1, given by We interpret (χ)| Γ = 0 and χ i = 0 for i ∈ I 0 ∪ I 1 .
Example 1 We provide an explicit example of A using a fictitious material. In this material, we assume that the stress-strain relationships in tangential and normal directions are independent. This assumption leads to a model which captures in-plane shearing whereas out-of-plane stress components follow a one-dimensional Hooke's law. This behavior is described by the following constitutive laws Here, µ and λ (respectively µ ⊥ and λ ⊥ ) are the Lamé parameters describing the stress-strain relationship tangential (and normal) to the manifold. The inverse relations, mapping stresses to strains, are then given by Here, I d is the identity tensor in R d×d . We remark that in this example, we Finally, we consider the symmetry of the stress tensor. Since the lowerdimensional manifolds model objects with finite width, the limit argument used to prove symmetry of the stress tensor is only valid within manifolds (and not transversely). Consequently, symmetry of the stress tensor is imposed within each manifold, expressed as: We remark that for i ∈ I 0 ∪ I 1 , this equation is trivial since either σ i does not exist or is a vector. For i ∈ I 2 , this equation evaluates the asymmetry of the in-plane components. Gathering (3.5), (3.6), and (3.8), we arrive at the strong form of the generalized system of equations: We emphasize that is defined as the square root of the cross-sectional measure, leading to the appearance of 2 in the second equation. To close the system, the boundary conditions are given by System (3.9) has a structure similar to (3.2) in that it is composed of constitutive law(s) complemented with a differential and algebraic constraint. This structure is common in models concerning linear elasticity with relaxed symmetry [2,5]. We will show in the next section that the system indeed corresponds to a symmetric saddle-point problem.

Variational Formulation
With the goal of obtaining a mixed finite element discretization, this section presents the weak formulation of the continuous problem. In order to do this, we introduce several analytical tools. First, the relevant function spaces are defined as well as the notational conventions concerning inner products. Next, we derive the variational formulation of (3.9) and show that it corresponds to a symmetric saddle point problem.

Function Spaces
The function spaces relevant for this problem are constructed as products of familiar function spaces on the d-dimensional manifolds. In particular, we define where Σ denotes the function space for the stress, U contains the displacement, and R is the function space for the Lagrange multiplier enforcing symmetry of the stress tensor. The exponent k d is given by [2,5].
The mixed-dimensional L 2 -inner products on Ω and Γ are defined as the sum of inner products over all corresponding manifolds: Here, the implicit assumption is made that the contribution is zero for all manifolds on which f is undefined. For example, for σ, τ ∈ Σ, the inner product (σ, τ ) Ω has no contribution on Ω 0 . Likewise for functions in R, the inner product is zero on manifolds Ω i with i ∈ I 0 ∪ I 1 .
For functions σ, τ ∈ Σ, we note that they are defined on both Ω and Γ . For convenience, we introduce the combined inner product which, in the case of the operator A, is understood as The inner products naturally induce the L 2 -type norms · Ω , · Γ , and · Ω×Γ . With these norms, we assume that A is continuous and coercive with respect to the norm · Ω×Γ . Thus, for all σ, τ ∈ Σ, we have: We emphasize that the constants within these bounds are independent of .

Identifying the Symmetric Saddle Point Problem
In this section, we make two key observations which allow us to derive a variational formulation of (3.9) which is symmetric. First let us consider the terms containing u in the stress-strain relationships (3.9a) and (3.9a). We multiply these terms with τ ∈ Σ and n · τ ∈ L 2 (Γ ), respectively, and integrate to obtain the following integration by parts formula: Here D· is the mixed-dimensional divergence operator from (3.5). Secondly, we introduce the operator skw which evaluates the asymmetric part of a matrix. More specifically, for a matrix B ∈ R n×d with components b ij , let This operator is naturally lifted to skw : Σ → R. Next, we turn our attention to the term in (3.9a) containing the asymmetric variable χ. Let us multiply this term with τ ∈ Σ and integrate over Ω 2 ∪ Ω 3 . With the introduction of r = 1 2 skwχ ∈ R, we obtain By employing test functions (τ, v, s) ∈ Σ × U × R, the integration by parts formula (4.3), and the operator skw, we obtain the following variational formulation of the problem (3.9): We identify system (4.4) as a saddle point problem by introducing the bilinear forms a : Σ × Σ → R and b : Σ × (U × R) → R: The problem (4.4) can then be rewritten to the following, equivalent formulation: for all (τ, v, s) ∈ Σ × U × R.

Well-Posedness
In this section, we show well-posedness of the continuous formulation (4.4).
The key is to associate appropriately weighted norms to the function spaces introduced in the previous section. In the mixed-dimensional setting considered here, let us endow Σ, U , and R with the following norms The proof of well-posedness consists of proving sufficient conditions on the bilinear forms a and b from (4.5) to invoke standard saddle-point theory. First, we show continuity of the operators, followed by ellipticity of a and inf-sup on b.
Lemma 1 (Continuity) The bilinear forms a and b from (4.5) are continuous with respect to the norms given by (5.1).
Proof The continuity of a follows from (4.2). For the blinear form b, we derive Next, we focus on the bilinear form a. For the purposes of our analysis, it suffices to show that a is elliptic on a specific subspace of Σ generated by b. This is formally considered in the following lemma.
Theorem 1 (Ellipticity) Given the bilinear forms a and b from (4.5). If then the following ellipticity bound holds a(σ; σ) σ 2 Σ . Proof We set s = 0 in condition (5.2). The assumption holds for all v ∈ U , thus noting that D · σ ∈ i∈I L 2 (Ω i ) = U andˆ max > 0, we obtain 3) The proof is concluded by combining (5.3) with the coercivity of A from (4.2).
With the properties of a proven, we continue by considering an inf-sup condition on the bilinear form b. This is shown in the following theorem, which relies on the constructions from Lemmas 2 and 4, presented afterwards.
Proof The proof consists of constructing a suitable τ ∈ Σ for a given pair (u, r) ∈ U × R. Its construction is based on constructing two auxiliary functions η, ξ ∈ Σ using the techniques from Lemmas 2 and 4. Setting τ as the sum of these two functions then yields the result. First, Lemma 2 allows us to construct η ∈ Σ such that Secondly, we choose ξ ∈ Σ using Lemma 4 with given (r − −1 skwη) ∈ R such that skwξ = r − skwη (5.5a) D · ξ = 0 (5.5b) Furthermore, the bound on τ is derived using (5.4) and (5.5c) The proof is concluded by combining (5.6) and (5.7).
Lemma 2 For each u ∈ U , a function η ∈ Σ exists such that Proof Considering u ∈ U given, the function η is constructed hierarchically. For each dimension d, we first set an interface value φ on Γ d , followed by a suitable extension into Ω d .
0. Given i ∈ I 0 , we construct the adjacent interface functions φ j ∈ L 2 (Γ j ) such φ j = −ˆ max u i for a chosen j ∈Ĵ i with (ˆ )| Γj = (ˆ max )| Ωi and zero for all other j ∈Ĵ i . Repeating this construction for all i ∈ I 0 , it follows that 1. We continue with i ∈ I 1 and perform the following two steps. First, the function η i is constructed as the bounded H(div, Ω i )-extension of the given φ j with j ∈J i . We use the extension operator as described in [19] (Section 4.1.2), giving us the properties Repeating these two steps for all i ∈ I 1 gives us the bound in which the second and third inequalities follow from (3.4) and (5.10c). 2. For n = 3, repeat the previous step for all i ∈ I 2 to obtain η i and φ j with j ∈Ĵ i . 3. The construction of η is finalized with its top-dimensional components η i with i ∈ I n . Let the pair (η i ,ṽ i ) ∈ (H(div, Ω i )) n × (L 2 (Ω i )) n be the weak solution to the Poisson problem: This problem is solved for all i ∈ I n . We then recall that =ˆ max = 1 in Ω n and exploit the elliptic regularity of (5.13) (see e.g. [13]) to obtain η Ω n + ˆ −1 max ∇ · η Ω n = η Ω n + ∇ · η Ω n φ Γ n−1 + u Ω n . (5.14) With η constructed, we consider its two main properties. First, by (5.9a), (5.11), and (5.13b), we have Secondly, we find the following bound from (5.9b), (5.10c), (5.12), and (5.14): thereby concluding the proof.
Before introducing the second ingredient used in the proof of Theorem 2, we require several key analytical tools, organized in the following diagram: The function spaces (Θ and W ) and mappings (Ξ, D·, and D×) are defined next. Let the auxiliary space Θ be given by We emphasize that for an element θ ∈ Θ, this definition implies that θ i is a 2-vector for i ∈ I 2 and a 3 × 3 tensor for i ∈ I 3 . Next, we follow [5] by introducing the mapping Ξ and its right-inverse Ξ −1 as: with w i, the tangential components of w i with respect to Ω i . W is defined as the space of functions that lie in the image of the inverse operator and have a mixed-dimensional curl in Σ, i.e.
We remark that for w ∈ W , w i is a 3-vector for i ∈ I 3 and a 3 × 3 tensor for i ∈ I 3 .
Next, we introduce a divergence-like operator D· : Θ → R given by Here,ň i is the unique unit vector normal to Ω i that forms a positive orientation with the chosen basis of the tangential bundle. By definition, this divergence operator maps from Θ to R, and we emphasize that D · θ is a vector for d = 3 and a scalar for d = 2. Finally, the mixed-dimensional curl (D×) of w ∈ W (see e.g. [14,8]) is given by , e.g. ∇ ⊥ is a rotated gradient operator. We note that all differential operations are performed rowwise. Hence, for n = 3, the mixed-dimensional curl maps to a 3 × 3 tensor in Ω 3 , a 3 × 2 tensor in Ω 2 (in local coordinates) and a 3-vector in Ω 1 (in local coordinates). Thus, an exact correspondence with the function space Σ is obtained, as reflected in the diagram.

Lemma 3
The operators in diagram (5.17) enjoy the following two properties for all w ∈ W : Proof The top row of (5.17) uses the differential operators from the mixeddimensional De Rham complex [8]. The first equality then follows from the fact that exact forms are closed. It remains to show commutativity. By the definition of Ξ, see e.g. [5,7], we have in which the subscript d denotes a restriction of the operator to Ω d . Furthermore, we note that on Ω i with i ∈ I 2 , the skw operator evaluates the asymmetry with respect to the tangent bundle of Ω i . In turn, we have for M ∈ R 3×3 that (skw 2 M )| Ωi =ň i · skw 3 M . This gives us Lemma 4 Given r ∈ R, a function ξ ∈ Σ exists such that Proof We give the proof for n = 3, the case n = 2 being simpler. The strategy is to exploit the properties shown in Lemma 3 and first construct a bounded θ ∈ Θ such that D · θ = 2 r. Then, by setting w = Ξ −1 θ and ξ = −1 D × w, we obtain two of the desired properties The estimate will then follow from the boundedness of θ. The construction of θ proceeds according to the following three steps, consisting of an interface function φ ∈ H 1 (Γ 2 ) which serves as a source function for θ i with i ∈ I 2 and a boundary condition for θ i with i ∈ I 3 .
1. We start by defining a scalar function φ in the trace space H 1 (Γ 2 ). We let φ vanish at all intersections and extremities, i.e. φ is in the function space Φ given by Now, let φ be the solution to the following minimization problem: with Π Ri the projection onto constants on Ω i . Due to the regularity of this problem and the imposed constraint, we have 2. For each i ∈ I 2 , we construct a function θ i using φ as a source function. Specifically, let (θ i , p i ) ∈ (H 1 0 (Ω i )) 2 × L 2 (Ω i ) be the weak solution to the Stokes problem: The following bound is then satisfied from the regularity of (5.28a) (see e.g. [13]) combined with (5.27b) 3. To finalize θ ∈ Θ, we create θ i for i ∈ I 3 using φ from the first step as a boundary condition. Let θ i ∈ (H 1 (Ω i )) 3×3 and an auxiliary pressure variable p i ∈ (L 2 (Ω i )) 3 be the weak solution to the following Stokes problem: Recall thatň i , the unique normal vector of Ω i , and n j , the normal vector defined on Γ j , are equal up to sign. This problem is well-posed since ∂Ω i \Γ 2 has positive measure, for each i ∈ I 3 , by assumption. We have the following bound due to the regularity of the Stokes problems and the fact that = 1 in Ω 3 Combining all θ i from the final two steps gives us θ ∈ Θ. We first note that the properties (5.27a), (5.28b), and (5.30b) result in Hence, we have D · θ = 2 r and it follows from (5.24) that setting w = Ξ −1 θ and ξ = −1 D × w provides the first two properties. The bound follows due to (5.27b), (5.29) and (5.31) With the proven properties of the bilinear forms a and b, the main result of this section is summarized by the following theorem: 4) is well-posed with respect to the norms (5.1). That is, a unique solution exists satisfying the bound Proof It suffices to show continuity of the right-hand side of (4.4) with respect to the norms above. For that purpose, we derive the following bound on the first term using Cauchy-Schwarz and a trace inequality: Moreover, from (3.4), the second term is bounded as follows The result now follows readily from these two estimates, Theorems 1 and 2, and standard saddle point theory [7].

Discretization
In this section, we discretize the continuous problem (4.4) using conforming finite elements. The notion of conformity and the choice of mixed finite element spaces is discussed in Section 6.1 and the resulting discretized problem is presented and analyzed in Section 6.2.

Discrete Spaces
For each i ∈ I, we introduce a shape-regular, simplicial grid Ω i,h which tessellates Ω i . The union of meshes of a given dimension d (with 0 ≤ d ≤ n) is denoted by Ω d h = ∪ i∈I d Ω i,h and we define Ω h = ∪ i∈I Ω i,h . We let the grid respect all lower-dimensional features and be matching across all interfaces. The tesselation of Γ is thus given by Γ h = Γ ∩ ∂Ω h . Moreover, there is an equivalence between Γ j,h and Ω i,h for all j ∈Ĵ i . The typical mesh size is denoted by h and we use h as a subscript to indicate the discretized counterpart of functions and function spaces.
With the aim of obtaining a stable and conforming method, we search for a discrete solution in subspaces of the function spaces defined in Section 4. We choose discrete function spaces (Σ h , U h , R h ) on the grid Ω h according to the following three conditions: (S1) The finite element spaces are conforming, i.e.
(S2) Σ h and U h are such that for each i ∈ I: forms a stable pair for the two-dimensional Stokes problem for each i ∈ I 2 . (S4) Σ i,h × U i,h × R i,h forms a stable triplet for three-dimensional, mixed elasticity for each i ∈ I 3 .
We provide an exemplary family of finite elements satisfying all four conditions. This choice is most concisely described using the notation of finite element exterior calculus [3]. A translation to more conventional nomenclature is provided afterwards, for convenience. Given a polynomial degree k ≥ 0, let In other words, for n = 3: Σ i,h with i ∈ I 3 corresponds to three rows of Nedelec elements of the second kind (N 2 f k+1 [16]) with degrees of freedom on the faces. For i ∈ I 2 , it is three rows of Brezzi-Douglas-Marini elements (BDM k+n−1 [10]). Finally Σ i,h with i ∈ I 1 is given by a triplet of continuous Lagrange elements (P k+2 ). The spaces U i,h and R i,h are defined for i ∈ I as three and k di rows, respectively, of discontinuous Lagrange elements (P −(k+n−di) ).
In this case, the auxiliary space W h of (S3) is explicitly given by For n = 3, W i,h is thus given by three rows of (first kind) edge-based Nedelec element (N 1 e k+2 ) for i ∈ I 3 and by three instances of Lagrange elements (P k+n ) for i ∈ I 2 . The lowest order choice in this family, i.e. with k = 0, is presented in Table 1.
A reduced family of finite elements arises by noting that all stability conditions remain valid after the polynomial order of the trace onto Γ is reduced by one. Table 2 presents the lowest order member of this family.
Properties (S1)-(S3) can be verified with the use of the presented tables. Finally, these families correspond to the stable triplets analyzed in [3] (Sections 11.6-11.7), hence (S4) holds as well.

Discrete Problem
Since the finite elements described above are contained in the continuous spaces from Section 5 by (S1), the discrete formulation of the model prob- Table 1 The finite element spaces chosen for n = 2 and n = 3 of lowest order within the family (6.1). The negative orders in the subscript denote discontinuous Lagrange elements. On the zero-dimensional manifolds, the polynomial order is redundant since any finite element space corresponds to point evaluation there. Table 2 The lowest-order finite element spaces of the reduced family for n = 2 and n = 3.
The superscript minus indicates that the trace of the finite element space onto Γ h is reduced by one order.
lem is a direct translation of (4.4): We note that the saddle-point structure of this problem has not changed, and can readily be uncovered using the bilinear forms from (4.5).

Stability
We continue with the analysis concerning the well-posedness of (6.2). Let us recall the norms from (5.1) for convenience Theorem 4 (Ellipticity) If the discrete spaces satisfy (S2) and then the following ellipticity bound holds Proof By (S2), we have D · Σ h ⊆ U h and the same arguments are used as in Theorem 1.
The next step is to consider the inf-sup condition for the bilinear form b in the discrete case.
Theorem 5 (Inf-Sup) If the discrete spaces satisfy conditions (S1)-(S4), Proof With u h ∈ U h and r h ∈ R h given, we follow a similar strategy as in the proof of Theorem 2. Here, we rely on Lemmas 5 and 6 to provide η h ∈ Σ h to gain control of u h and a divergence-free function and satisfy the bound Following the same steps as in Theorem 2, we define The proof is concluded by combining the above.
Proof We use the same steps as in Lemma 2 to hierarchically construct η h ∈ Σ h . A concise exposition follows, starting with d = 0. For each i ∈ I d , we use (S2) to first construct a discrete φ j,h in the trace space (n · Σ h )| Γj for all j ∈Ĵ i such that in which ∇ · η i,h and φ j,h are understood as zero for i ∈ I 0 , j ∈J i . For i ∈ I d+1 , the function η i,h is then defined as the bounded, discrete H(div, Ω i,h )extension from [19] (Section 4.1.2). These steps are repeated by incrementing d until φ h is completely defined on Γ h . Finally, we solve a discrete Poisson problem for each i ∈ I 3 , in analogy with (5.13), to complete η h ∈ Σ h . It follows by the same arguments as in Theorem 2 that the constructed η h satisfies (6.4).

Lemma 6
Given r h ∈ R h , a function ξ h ∈ Σ h exists such that Proof In this proof, we make extensive use of the discrete space W h from stability requirement (S3). In particular, we will first introduce w h ∈ W h such that D×w h controls r i,h for i ∈ I 2 . Then, a correction is introduced using (S4) in order to control r i,h for i ∈ I 3 as well. For brevity, we omit the subscript h on all variables within this proof. As in Lemma 4, we start with an interface function φ defined on the trace mesh Γ 2 h which serves first as a source function and second as a boundary condition. We proceed according to the following four steps.
1. We consider functions on Γ 2 h in the trace space of (W h )| Ω 3 that vanish at all intersections and extremities. Let us therefore introduce the function space Φ h as It is important to note that the two instances of n in this definition are different. In particular, the former is defined as normal to ∂Ω i , of which Γ j is a subset, whereas the latter is normal with respect to the boundary ∂Γ j . We note that, since D × W h ⊆ Σ h by (S3), we have W i,h ⊆ (H(curl, Ω i )) 3 for i ∈ I 3 . In turn, it follows in the discrete setting that Φ j,h ∈ (H(div, Γ j )) 3 for all j ∈J i , with the divergence tangential to Γ j . Using this observation, we let φ solve the following minimization problem subject to Π Ri (skw 2 ϕ i + 2 r i ) = 0, ∀i ∈ I 2 . (6.8) In other words, a finite-dimensional problem is solved for each i ∈ I 2 to obtain a bounded distribution that has an average asymmetry corresponding to 2 r i . In particular, we obtain the following two properties We generate w i for i ∈ I 2 using the distribution φ from the first step as a source term. Introducing Θ i,h = ΞW i,h , it follows from (S3) that Θ i,h ×R i,h is a stable pair for the discretization of the Stokes problem (5.28a). Hence, we can find θ i ∈ Θ i,h such that 3. For i ∈ I 3 , let w i be given by any bounded extension of φ into H(curl, Ω i ), i.e. w i is chosen for all i ∈ I 3 such that 4. Finally, we gain control of r i for i ∈ I 3 . We recall that by stability condition (S4), the spaces Σ i,h × U i,h × R i,h form a stable triplet for the mixed formulation of elasticity with relaxed symmetry. Considering Γ 2 as a zero traction boundary condition, we use the inf-sup condition associated to this stability to form a function ψ i ∈ Σ i,h such that ∇ · ψ i = 0 (6.14a) To complete the mixed-dimensional function ψ ∈ Σ h , we set ψ i = 0 for i ∈ I 3 .
Using the above ingredients, we set ξ = ψ + −1 D × w and obtain the first two desired properties: The bound now follows by the estimates given in each step: The two previous theorems provide the sufficient ingredients to show stability of the discretization, formally presented in the following theorem.
Theorem 6 (Stability) If the discrete spaces satisfy the conditions (S1)-(S4), then the resulting mixed finite element method is stable, i.e. a unique solution exists satisfying the bound with a constant independent of the grid size h.
Proof Using Theorems 4 and 5, this result follows from standard saddle point theory [7].

Convergence
By consistency of the discretized problem (6.2) with respect to the continuous formulation (4.4) and stability from Theorem 6, we have shown that the proposed mixed finite element discretization is convergent. In turn, this section is devoted to obtaining the rates of convergence through a priori error estimation. Let Π R h and Π U h be the L 2 -projection operators onto the finite element spaces R h and U h . Under the assumption of sufficient regularity, we introduce the canonical projection operator Π Σ h such that the commutativity properties hold for σ ∈ Σ: Π U i,h ∇ · σ i = ∇ · Π Σ i,h σ i , Π U i,h (n · σ)| Γj = (n · Π Σ ,h σ)| Γj , ∀i ∈ I, j ∈Ĵ i .
A direct consequence of these two properties is that for sufficiently regular σ ∈ Σ, we have the commuting property Using · ρ,Ω as short-hand notation for the H ρ (Ω)-norm, the projection operator Σ h has the following approximation properties for 0 ≤ ρ ≤ k i + 1:
Theorem 7 (Convergence) Given (σ, u, r) as the solution to (4.4). Under the assumption of sufficient regularity, then the scheme converges optimally, i.e. the finite element solution (σ h , u h , r h ) to (6.2) satisfies the following estimate h ki h σ ki+1,Ωi + ˆ −1 max D · σ ki,Ωi + j∈Ĵi n · σ ki,Γj + ˆ max u ki,Ωi + r ki,Ωi Proof We restrict our choice of test functions (τ h , v h , s h ) to the finite element spaces and subtract the discrete equations (6.2) from the continuous equations (4.4). We then obtain The projection operators from (6.17) are then used to project the true solution onto the finite element spaces. Introducing σ Π := Π Σ h σ, u Π := Π U h u, and r Π := Π R h r, we rewrite the above equation to Let τ b,h be the discrete stress from the construction in Theorem 5 with the following properties The discrete test functions are then chosen to be with δ > 0 a constant to be determined later. Substituting this choice of functions into the left-hand side of (6.18) gives us Next, due to (6.16) and (S2), we have D · (σ Π − σ h ) = Π U h D · (σ − σ h ) = 0. Hence, the coercivity of A from (4.2) allows us to bound a(σ Π − σ h ; σ Π − σ h ) from below by σ Π − σ h 2 Σ . We then obtain the following bound with respect to the right-hand side of (6.18): Next, we use the continuity of the forms a and b from Lemma 1 to bound the right-hand side further An application of Young's inequality and rearranging terms then gives Using (6.19b) and setting δ sufficiently small leads us to With the triangle inequality, we thus obtain the estimate An application of the approximation properties (6.17) then finishes the proof.