Mixed displacement–rotation–pressure formulations for linear elasticity

We propose a new locking-free family of mixed finite element and finite volume element methods for the approximation of linear elastostatics, formulated in terms of displacement, rotation vector, and pressure. The unique solvability of the three-field continuous formulation, as well as the well-definiteness and stability of the proposed Galerkin and Petrov–Galerkin methods, is established thanks to the Babuška–Brezzi theory. Optimal a priori error estimates are derived using norms robust with respect to the Lamé constants, turning these numerical methods particularly appealing for nearly incompressible materials. We exemplify the accuracy (in a suitably weighted norm), as well the applicability of the new formulation and the mixed schemes by conducting a number of computational tests in 2D and 3D, also including cases not covered by our theoretical analysis. c ⃝ 2018 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons. org/licenses/by/4.0/). MSC: 65N30; 65N12; 76D07; 65N15


Introduction
The numerical solution of elasticity-based problems encompasses well-documented difficulties. For instance, for pure-displacement formulations, the use of classical finite element discretisations based on piecewise linear and continuous elements, ensures accuracy only for moderate values of the Poisson ratio ν. As ν → 0.5, that is, when the Lamé dilation modulus λ → ∞, and the elastic material becomes nearly incompressible, the numerical scheme might generate spurious solutions (unphysically small deformations related to the well-known locking phenomenon, see for instance [1]). A number of appropriate formulations together with their associated numerical methods are available to overcome this issue. Notably, choosing a mixed scheme would produce accurate solutions even for nearly incompressible materials, and at the same time, one accommodates the direct approximation of auxiliary variables of interest such as pressure, stress, or rotations.
One of the most common mixed approaches for linear elasticity is the Hu-Washizu formulation [2,3]. Some popular methods based on such formulation include the enhanced assumed strain method [4], the assumed stress method [5], the mixed-enhanced strain method [6], the strain gap method [7], and the B-bar scheme [1]. Some of these methods actually coincide under certain conditions (see the discussions in e.g. [8][9][10]). The well-posedness for this class of formulations has been established in [11], where it is also shown that a modified version of the Hu-Washizu formulation is more amenable for obtaining uniform convergence in the incompressibility limit. Alternatively, other mixed approaches (such as the Hellinger-Reissner principle) can be employed to obtain robust methods with respect to the Lamé constants.
Schemes more closely related to the present contribution, state the problem using stresses and rotations. We mention for instance mixed formulations based on stress [12][13][14], the augmented scheme in [15], a family of pseudostress-based methods from [16], displacement-pressure mixed formulations [17]; and the first-order least squares presented in [18]. More recent least squares schemes in connection with the present context include saddle-point least squares methods [19], mixed approaches also considering anisotropy, large strains and quasiincompressibility, or others applied specifically to plates [20,21]. We refer as well to other locking-free methods for plate models [22], and to the membrane elements introduced in [23], also including the rotation tensor as an additional field.
In contrast to the brief literature survey given above, here we advocate a novel formulation of the elasticity equations in terms of displacement, rotation vector, and pressure (similar ideas in the context of vorticity-based formulations for Stokes and Brinkman equations can be found e.g. in [24][25][26][27]). This three-field formulation has a resemblance with the displacement, pressure, and vorticity momentum formulations for acoustic fluid-structure interaction studied in [28]. However in that reference, the system is solved for the fluid displacement and the vorticity momentum arises as the Lagrange multiplier imposing an irrotationality constraint.
In our case, after regarding the pressure together with the rotation vector as a single auxiliary unknown (defined in an appropriate product functional space), we are able to analyse the solvability of the resulting mixed variational formulation simply appealing to the classical Babuška-Brezzi theory for saddle-point problems. Thanks to a rescaling of the rotation vector norm, the well-posedness result and the continuous dependence on the data turn out to be independent of the Lamé constants. This analysis is valid for (possibly non-homogeneous) Dirichlet boundary conditions, and bounded Lipschitz domains.
Concerning numerical approximation, we first introduce a family of finite elements given by piecewise continuous polynomials of degree k ≥ 1 for the displacement, and piecewise polynomials of degree k − 1 for the rotation and pressure. The unique solvability of the finite element scheme is then established using analogous techniques as in the continuous case, that is, exploiting a weighted norm. In addition, we prove optimal a priori error estimates with constants fully independent of the Lamé coefficient λ; guaranteeing robustness of the method also in the nearly incompressible limit. Nevertheless in the case of full incompressibility, both continuous and discrete problems are not necessarily well-defined.
We remark that in the two-dimensional case, the computational cost of the proposed FE method in its lowest order configuration is 2|N h | (where N h denotes the set of vertices in the mesh and |N h | its cardinality), which is lower than, for instance, the MINI element for displacement-pressure formulations (accounting for 7|N h | local degrees of freedom). Furthermore, even if our method involves two additional unknowns (pressure and rotation vector), these can be statically condensed at the implementation stage without incurring in additional computational cost, thus turning the proposed discretisation into a very competitive method.
A further goal in this contribution is to construct a finite volume element (FVE) method specifically tailored for elasticity equations. FVE schemes correspond to Petrov-Galerkin formulations where the trial space is constructed using a primal partition of the domain, whereas the test space is associated with either a dual mesh or a dual basis. Depending on the particular kind of dual grid, the transfer operator between trial and test spaces possesses different interpolation properties which are used in recasting a preliminary pure finite volume formulation into a Petrov-Galerkin one. In general, these methods enjoy some features shared by finite element and finite volume schemes, including local flux conservation properties, liberty to choose different numerical fluxes and dual partitions associated to unstructured primal meshes; and several others (see for example [29]). Discretisation schemes following this principle have been systematically employed in numerous fluid flow problems, including Stokes, Navier-Stokes (see e.g. [30][31][32][33][34]) and also in coupled flow-transport systems arising from diverse applications (see [35][36][37]). However, up to our knowledge, the only contributions addressing FVE-like discretisations for solid mechanics are the hybridstress finite volume method for linear elasticity on quads studied in [38]; and [39], where two alternative stabilisation approaches based on nodal pressure and dual bases and meshes are applied to construct inf-sup stable approximations for nearly incompressible linear elasticity. The class of finite volume element methods we introduce here is based on the lowest-order mixed finite element method discussed above. As in well-established FVE schemes for Stokes equations (cf. [31,32]), it turns out that the two schemes differ only by the assembly of the forcing term, and therefore straightforward derivation of stability properties and energy estimates in natural norms can be done exploiting the results obtained for the family of mixed finite elements. In addition, the FVE scheme features mass conservativity on the dual control volumes, suitability for irregular domains and unstructured partitions, and robust approximations of displacements. We also observe that the proposed schemes perform very well across the scope of regimes considered in our numerical simulations.
Outline. We have structured the contents of the manuscript in the following manner. To simplify the exposition, a few recurrent notations and useful identities are recalled in the remainder of this section. In Section 2 we lay out the precise form of the linear elasticity equations that we will focus on, we derive a suitable mixed weak formulation, and provide its solvability analysis. A Galerkin method is introduced in Section 3, where we also obtain stability properties and a priori error estimates. Section 4 concentrates on the development of a low-order FVE scheme, and its accuracy is studied in connection with the properties of the FE method. We briefly discuss up to which extent the definition and construction of the proposed FE and FVE schemes needs to be modified in order to accommodate the study of mixed displacement-traction boundary conditions. Some indications on how the analysis could be extended are also addressed. Finally, the convergence and robustness of the proposed methods is illustrated via a set of insightful computational tests collected in Section 5, including also some comparisons against other methods.
we recall the following notation for differential operators: We also recall a version of Green's formula given in e.g. [40,Theorem 2.11]: ∫ 1) and the following useful identity 2. The model problem

Derivation of a displacement-rotation-pressure formulation
We assume that an isotropic and linearly elastic solid occupies a polyhedral bounded domain Ω of R d , with Lipschitz boundary ∂Ω . Determining the deformation of a linearly elastic body subject to a volume load and with given boundary conditions, and adopting the hypothesis of small strains, results in the classical linear elasticity problem, formulated as follows. Given an external forcef and a prescribed boundary motion g, we seek the displacements u such that where ε(u) = 1 2 (∇u + ∇u t ) is the infinitesimal strain tensor, I denotes the d × d-identity matrix, and µ, λ are the Lamé coefficients (material properties of the solid, and here assumed constant).
Next, and following the seminal work [12], one notices that using the identity div(ε(u)) = and dividing the momentum equation by λ + µ, we can rewrite (2.1) in the form of the well-known Cauchy-Navier (or Navier-Lamé) equations where the body load has been rescaled as f = 1 λ+µf . We then proceed to define the auxiliary scaling parameter η := µ λ+µ > 0, and recast (2.2) in a displacement-pressure formulation (considering p = −div u as the solid pressure) as follows At this point, and with the aim of deriving formulations whose stability holds independently of the Lamé coefficient λ, we introduce the field of rescaled rotations ω := √ η curl u, as an additional unknown in the problem. Exploiting (1.2) and the definition of pressure in terms of displacements, we observe that (2.3) is fully equivalent to the following set of governing equations, in their pure-Dirichlet case. Find the displacement u, the rotation ω and the pressure p such that (see [18]): (2.7) The theoretical analysis will be restricted to the case of clamped boundaries, g = 0. The case of non-homogeneous Dirichlet boundary conditions can be analysed in an analogous manner after introducing a suitable displacement lifting. On the other hand, the incorporation of mixed (displacement-traction) boundary conditions will be addressed in Sections 3.2 and 4.2.

Weak form of the governing equations
Let us introduce the functional spaces H := H 1 0 (Ω ) d , Z := L 2 (Ω ) d , and Q := L 2 (Ω ), where Z and Q are endowed with their natural norms, and we recall the definition of the norm in the product space Z × Q as ∥(θ , q)∥ 2 Z×Q := ∥θ ∥ 2 0,Ω + ∥q∥ 2 0,Ω . On the other hand, for H we consider the following η-dependent scaled norm (see for instance, [40, Remark 2.7]): ∥v∥ 2 H := η∥curl v∥ 2 0,Ω + ∥div v∥ 2 0,Ω , which is motivated by the natural energy form of the momentum conservation equation in (2.3), and the rescaled rotation. Notice that the stability of continuous and discrete problems will be stated in terms of this norm.
We proceed to test (2.6) against adequate functions, to integrate by parts in two terms, and to take into account the boundary conditions (2.7) in such a way that the resulting mixed variational formulation reads as follows. Find Introducing the bilinear forms a : for all v ∈ H, ω, θ ∈ Z, and p, q ∈ Q; we realise that the variational problem above can be recast as: Find Remark 1. Note that the natural regularity for displacements in (2.8)-(2.9) is actually H 0 (curl, Ω )∩H 0 (div, Ω ), where H 0 (curl, Ω ) := {v ∈ H(curl, Ω ) : (v × n)| ∂Ω = 0} and H 0 (div, Ω ) := {v ∈ H(div, Ω ) : (v · n)| ∂Ω = 0}. According to [40,Lemma 2.5], an algebraic and topological equivalence between this space and H = H 1 0 (Ω ) d holds under quite general assumptions on the domain: Ω only needs to be bounded and ∂Ω Lipschitz-continuous (see also [40,Remark 2.7]). In other instances (for example in the analysis of vector Laplacians, see e.g. [41, Section 2.3.2]) if tangential and normal components of the displacement are to be fixed on different parts of the boundary, then the domain convexity is also required. However that is not the case in the present study.
Remark 2. In the incompressibility limit ν = 0.5, the problem defined in (2.8)-(2.9) reduces to ∫ However, it is not difficult to see that u satisfies ∇(div u) = −f in Ω , which is not well-posed in the space H. Moreover, note that after rescaling, the body force also goes to zero in the incompressibility limit.

Well-posedness
The unique solvability of problem (2.8)-(2.9), together with the continuous dependence on the data will be established using the well-known Babuška-Brezzi theory.
We first observe that the bilinear forms a(·, ·), b(·, ·) and the linear functional F(·) are all bounded by positive constants independent of η (and therefore independent of the Lamé coefficient λ). In fact, since 0 < η < 1, it is easy to check that In addition, the bilinear form a(·, ·) is (Z × Q)-elliptic, uniformly with respect to the scaling parameter η, as stated in the following result.
Moreover, an inf-sup condition holds for the bilinear form b(·, ·). Lemma 2.2. There exists C > 0, independent of η, such that Proof. Let us consider a generic v ∈ H and definẽ We immediately notice that and from the definition of b(·, ·), we readily obtain which finishes the proof. □ We are now in a position to state the solvability of the continuous problem (2.8)-(2.9).
Theorem 2.1. There exists a unique solution ( (ω, p), u ) ∈ (Z × Q) × H to problem (2.8)-(2.9), which satisfies the following continuous dependence on the data Proof. By virtue of the general theory for saddle-point problems (see e.g. [42]), the desired result follows from a direct application of Lemmas 2.1 and 2.2. □ Owing to the well-known regularity for the elasticity equations (see e.g. [43], [44,Theorem 5.2]), the solution u of (2.8)-(2.9) belongs to H 1+s (Ω ) n , for some s > 0 depending on the geometry of Ω and on the Lamé coefficients (and consequently on η). Moreover, there exists C > 0 independent of f such that

Finite element discretisation
In this section, we introduce a Galerkin scheme associated to (2.8)-(2.9), we specify the finite dimensional subspaces to employ, and analyse the well-posedness of the resulting methods using suitable assumptions on the discrete spaces. The section also contains a derivation of error estimates.

Formulation, solvability, and error bounds
Given an integer k ≥ 1 and a set S ⊂ R d , the space of polynomial functions defined in S and having total degree ≤ k will be denoted by P k (S).
Next, we define the following discrete spaces: which are subspaces of H, Z and Q, respectively; and proceed to state a Galerkin scheme for the continuous variational Our next goal is to establish discrete counterparts of Lemmas 2.1 and 2.2, leading to the solvability and stability of the Galerkin method (3.1)-(3.2). Their proofs are obtained using the same arguments exploited in the continuous case. For completeness we provide the essential steps of the latter result.
Proof. For a generic v h ∈ H h , let us definẽ Then we readily notice that , and so, from the definition of the bilinear form b(·, ·), we arrive at the desired bound We can now state the unique solvability, stability, and convergence properties of the discrete problem (3.1)-(3.2), formulated in form of the three following theorems.
Moreover, there exists a constant C > 0, independent of h and η, such that In addition, the following approximation property is satisfied where ((ω, p), u) ∈ (Z × Q) × H is the unique solution of (2.8)-(2.9).
where s > 0 is such that the bound (2.12) is satisfied, and k ≥ 1 denotes the polynomial degree.
Proof. The result follows from Theorem 3.1 and the standard error estimates for the Scott-Zhang interpolant of u and the vectorial and scalar L 2 -orthogonal projections for ω and p, respectively; together with the additional regularity (2.12). □ To close this section, we observe that the convergence of the displacement approximation can be also measured in the L 2 (Ω ) d -norm, thanks to a classical duality strategy.
be the solutions of the continuous and discrete problems (2.8)-(2.9) and (3.1)-(3.2), respectively. Then, there exist constantss ∈ (0, 1], (depending on Ω and on η), and C > 0 (independent of h and η), such that Proof. Resorting to a duality argument, we first consider the following well-posed problem: Note that, as a consequence of (2.12), the unique solution of (3.3)-(3.4) features additional regularity. More precisely, we can assert that there exists ∈ (0, 1] as in (2.12), andC > 0 (independent of η), such that Next, and thanks to (3.4), we observe that for all where we have also employed (2.8) and (3.1). We then proceed to bound the second term in the right-hand side of (3.6). This is carried out by adding and subtracting (ξ , φ), and applying (3.3) where in the last step we have also used (2.9) and (3.2), valid for all z h ∈ H h . Hence, from (3.6)-(3.7), we can deduce that is the L 2 -orthogonal projection, and we can choose z h as the Scott-Zhang interpolation of z onto the piecewise linear and continuous vector fields. Thus, the proof follows from standard error estimates, the additional regularity (3.5), and Theorem 3.2. □ We point out that the values ∈ (0, 1] is associated to the regularity invoked in (3.5) when the data is

A discrete formulation with mixed displacement-traction boundary conditions
Let us now consider the case where a given displacement g is imposed only on a part of the boundary Γ D ⊂ ∂Ω , and set a given tractiont on the remainder of the boundary, say Γ N = ∂Ω \ Γ D . In this case (2.7) is replaced by where n denotes the outward unit normal on Γ N . This traction condition can be conveniently recast in terms of the field variables in the following way where the rescaled traction is t = 1 λ+µt , and where we have used the well-known identity ε(u)n = (∇u)n− 1 2 curl u×n.
The form of (3.9) readily implies that the displacement should now belong to the space M := H Γ D (curl, Ω )∩H Γ D (div, Ω ). It has been proved in [25, Lemma 1] (restricted to the 2D case) that there exists δ ∈ (1/2, 1] such that M is continuously imbedded in H δ (Ω ) 2 . If we set again homogeneous data on the Dirichlet boundary, we could then, as a first attempt, propose the following discrete modification of (3.1)-(3.2) incorporating mixed displacement-traction boundary conditions: [25]), and the diagonal bilinear form c : Since this bilinear form is non-symmetric and not necessarily semi-positive definite, the analysis of (3.10)-(3.11) does not fall in the same framework as (3.1)-(3.2). A possible way-around would be to define a fixed-point iteration scheme that assumes c(·, ·) as part of the linear functional. Then the solvability analysis can be carried out following e.g. [46]. Alternatively, one could introduce suitable Lagrange multipliers in order to deal with the boundary terms. Further investigation is necessary in this regard, and we simply mention that the implementation and numerical verification of test cases involving (3.10)-(3.12) will be addressed in Section 5, where we observe optimal convergence.

Formulation and main properties
In addition to the mesh T h (from now on, the primal mesh), we introduce another partition of Ω , denoted by T ⋆ h and referred to as the dual mesh, where for each element K ∈ T h we create segments joining its barycentre b K with the midpoints (2D barycentres) m F of each face F ⊂ ∂ K (or the midpoints of each edge, in 2D), forming four polyhedra (or three quadrilaterals, in the 2D case) Q z for z in the set of vertices of K , that is, z ∈ N h ∩ K . Then to each vertex s j ∈ N h , we associate a so-called control volume K ⋆ j consisting of the union of the polyhedra (quadrilaterals in 2D) Q s j sharing the vertex s j . A sketch of the resulting control volume associated to s j is depicted in Fig. 4.1(a).
In its lowest-order version, a FVE method for the approximation of (2.8)-(2.9) can be constructed by associating discrete spaces to a dual partition of the domain and notice that no additional space is introduced for the finite volume approximation of ω or p. Furthermore, we define the T ⋆ h -piecewise lumping map H h : H h → H ⋆ h which relates the primal and conforming dual meshes by for all v h ∈ H h , where χ j is the vectorial characteristic function of the control volume K ⋆ j and {ϕ j } j is the canonical FE basis of H h (cf. [32]). For any v ∈ H, this operator satisfies the interpolation bound (see e.g. [29]) ∥v − H h v∥ 0,Ω ≤ Ch|v| 1,Ω . In addition, since for the type of domains we are considering we can write H := H 1 0 (Ω ) d = H 0 (curl; Ω ) ∩ H 0 (div; Ω ), then [40,Remark 2.7] implies that the operator H h (·) also satisfies which plays a role in the convergence proof for the envisioned FVE method. The discrete FVE formulation is obtained by multiplying (2.4) by v ⋆ h ∈ H ⋆ h and integrating by parts over each K ⋆ j ∈ T ⋆ h , multiplying (2.5) by θ h ∈ Z h and integrating by parts over each K ∈ T h , and multiplying (2.6) by (1 + η)q h , for q h ∈ Q h , and integrating by parts over each K ∈ T h . This, along with identity (1.1), results in a Petrov-Galerkin formulation that reads as follows: Find where the bilinear formb : We also introduce the bilinear form B : which will be used to show that the Petrov-Galerkin formulation (4.2)-(4.3) can be regarded as a Galerkin method. We proceed to establish a relationship between the bilinear forms b(·, ·) and B(·, ·), which will be useful to carry out the error analysis in a finite-element-fashion. For the sake of brevity, only the proof for the two-dimensional case is provided. The proof for the three-dimensional case follows in an analogous manner, where we instead consider polyhedral control volumes and boundary surfaces rather than boundary edges.
Proof. First, let g be a function that is continuous on the interior of each quadrilateral Q j (as shown in Fig. 4.1(c)) with ∫ e g = 0 for any boundary edge e. Using Fig. 4.1(c), it is straightforward to show that the following relation holds: where m j+1 b K m j denotes the union of the line segments m j+1 b K and b K m j . We take m j+3 = m j in the case that the index is out of bound. Next, from the definition of the transfer operator H h (·), we find that In order to arrive at (4.4), we use the definition of B ( ·, · ) in combination with integration by parts and the fact that both q h and v h (s j ) are constant in the interior of each quadrilateral Q j , to obtain ] .
Since q h and θ h are constant on the edges of each element K ∈ T h , we can write Then, after one application of integration by parts and identity (1.1), we can assert that which completes the proof. □ is bounded uniformly with respect to η.

A FVE method with displacement-traction boundary conditions
Analogously to Section 3.2, we discuss here how our FVE scheme can be modified to incorporate mixed boundary conditions. We first define the spacẽ And then test (2.4) against v ⋆ h ∈H ⋆ h and integrate by parts, which leads to Next we reason as in the proof of Lemma 4.1 by considering the edges that coincide with the boundary segment Γ N separately. More precisely, by substituting H h v h and by definition of the traction t, we readily obtain for every v h ∈Ĥ h and every edge e of each K ∈ T h (cf. [32]), and the fact that p h ∈ Q h and ω h ∈ Z h are constant on each element K ∈ T h , which implies that

5)
where we have also used that the union of boundary edges of control volumes and the union of boundary edges of elements coincide. Consequently, we can combine (4.5)-(4.6) with (1.1) to finally obtain the following FVE formulation using mixed displacement-traction boundary conditions: where the newly introduced bilinear form C :Ĥ h ×Ĥ h → R is defined as Moreover, the linearity of u h ∈Ĥ h on each element K ∈ T h implies that This relation states that, also for mixed boundary conditions, the lowest-order FE and FVE schemes only differ by assembly of the right-hand side; which is not necessarily true for all nonsymmetric formulations.

Stability and convergence analysis
Back to the homogeneous Dirichlet case, our next goal is to prove a FVE-counterpart of Lemma 3.2, leading to the solvability and stability of (4.2)-(4.3). Recall that Lemma 3.1 establishes that the bilinear form a ( ·, · ) is (Z h × Q h )elliptic, uniformly with respect to η. Lemmas 4.1 and 3.2 readily imply that B(·, ·) satisfies an inf-sup condition, as stated in the following result.
Analogously to the previous section, the following two theorems formulate the unique solvability, stability, best approximation, and convergence properties of the discrete problem (4.2)-(4.3). Moreover, there exists a constant C > 0, independent of h and η, such that In addition, the following best approximation result is satisfied where ((ω, p), u) ∈ (Z × Q) × H is the unique solution to (2.8)-(2.9).
The next lemma establishes linear convergence of the lowest-order FVE method. (4.9) Next, after applying the inf-sup condition from Lemma 4.2 to Eq. (4.8), standard arguments imply that there exists a constant C 0 > 0, independent of h and η, such that Moreover, combining (4.8) with (4.9), and using (4.1) together with Lemma 3.1 implies that there exists a constant C 1 > 0, independent of h and η, such that Applying the triangle inequality to the convergence bound for the FE method established in Theorem 3.2 in combination with the inequalities (4.10) and (4.11) finishes the proof. □ To close this section, we prove an L 2 -estimate for the displacement error. For this purpose we first state a preliminary result (cf. [32]) that involves the transfer operator H h (·).

Lemma 4.3. For any function z h ∈ H h and any element K ∈ T h , one has
such that employing identity (3.7) and identity (4.12) yields which holds for all z h ∈ H h . In particular, we take the Lagrange interpolant of z, denoted by z I ∈ H h . Moreover, we use f K to denote the average of f on a given K ∈ T h . Then, by virtue of Lemma 4.3 and after integrating over K ∈ T h instead of over control volumes K ⋆ j ∈ T ⋆ h , we find that for some constant C 0 > 0, independent of h and η, Applying triangle inequality, using the estimates for the Lagrange interpolants, and exploiting the additional regularity, we get (4.14) and we arrive at the desired result after taking the L 2 -projections for ξ and φ, using interpolation properties, and employing (4.13) in combination with (4.14). □

Numerical tests
We report in this section some numerical examples which confirm our theoretical results, also including some additional cases not covered by our analysis.
Test 1A (accuracy assessment in 2D). For our first computational example we conduct a convergence test using a sequence of successively refined uniform partitions of the elastic domain Ω = (0, 1) 2 . We arbitrarily choose the Lamé parameters µ = 50, λ = 5000, so that η = 0.0099. This example focuses on the pure-Dirichlet problem (2.4)-(2.7), where we propose the following closed-form solutions satisfying the homogeneous Dirichlet datum, and where the forcing term f is constructed using these smooth functions and the linear momentum equation. The convergence study is performed for the FVE method (4.2)-(4.3) (of lowest order), and for the Galerkin schemes (3.1)-(3.2) of order k = 1 and k = 2. For a generic scalar or vectorial field v, on each nested mesh we will denote computed errors and experimental convergence rates as where e,ê stand for errors generated by methods defined on meshes with meshsizes h,ĥ, respectively; and we recall that ∥ · ∥ H denotes the η-dependent norm. These errors are tabulated by number of degrees of freedom in Table 5.1, which correspond here (and in all subsequent tests) to the dimension of the space Z h × Q h × H h . Apart from the displacement error measured in the L 2 -norm (whose error decays with order h k+1 as anticipated by Theorems 3.3 and 4.3), each individual error exhibits an O(h k ) rate of convergence, as predicted by the a priori error estimates stated in Theorems 3.2 and 4.2. Moreover, the errors produced by the first two methods practically coincide, which is explained by the fact that they only differ in the RHS assembly. For reference, in the top row of Fig. 5.1 we depict approximate solutions generated with the lowest-order FVE scheme.
Test 1B (accuracy in a 3D non-convex domain). We now consider again the pure Dirichlet case, now using a non-homogeneous datum set using the following closed-form displacement defined on an L-shaped domain of width 2, height 1.5 and depth 0.5: sin(π x) sin(π y) cos(π z) −x cos(π x) cos(π y) sin(π z) x y sin(π x) cos(π y) cos(π z) ⎞ ⎠ , which is also used to compute the exact rotation, pressure, and body load. The model constants are taken as in Test 1A, and we generate a sequence of refined meshes (non nested and unstructured) and produce a convergence study reported in Table 5 Test 1C (robustness with respect to η). In addition, these methods are robust with respect to the model parameters, which we confirm by a series of tests where we fix a Young modulus E = 10000, we vary the Poisson ratio ν, and measure the errors produced by the first order finite element method on an unstructured mesh of 33282 elements  Table 5.3). Furthermore, we also construct a different smooth forcing term f = 100(cos(x), cos(y)) t , independent of the model parameters, solve the discrete problem for relatively large Lamé constants (we recall that λ = Eν/[(1 + ν)(1 − 2ν)] and µ = E/(2 + 2ν)), and tabulate in the bottom block of Table 5.3 the obtained norms of the approximate solutions. We evidence stable and robust computations even in the nearly incompressibility limit. We have also arbitrarily set the scaling parameter to a very low value η =1e-15 (even if materials with so large differences between the shear and dilation moduli are rarely encountered), and have reproduced Test 1A (experimental convergence against a manufactured solution), observing that the methods still produce optimal convergence rates. We stress that these rates are optimal when measured in the H-norm. The computational cost associated with solving the linear systems arising from the FE and FVE discretisations can be significantly reduced through static condensation of the pressure and rotation blocks. The relevant systems assume  Test 1A (a, b, c). Approximate solutions computed with a FE scheme of order k = 1 for Test 1B (d, e, f).
with σ h := ( p h , ω h ) t . Since A is symmetric and positive definite, σ h can be eliminated from the first equation of (5.1) using σ h = −A −1 B t u h (recall that A is formed by the pressure and rotation mass matrices, so it is block-diagonal and could be easily inverted). Substituting this equation back into the second equation of (5.1) yields the displacement Schur complement system which is smaller, symmetric, and positive definite. Different methods can be employed to solve the Schur complement problem efficiently, also avoiding assembling S := BA −1 B t (see e.g. [47] for an application in elasticity).
Test 2 (2D beam bending). For the next computational example we study the displacement-rotation-pressure patterns of a rectangular beam (with length L = 10 and height l = 2) subjected to a couple (that is, a prescribed traction ( f (1 − y), 0) t , with f = 200) at one end, as shown in Fig. 5.2(a). We assume that the origin O is fully fixed and that the horizontal displacement is zero along the left edge of the domain Ω . Furthermore, on the remainder of the boundary we consider zero normal stresses incorporated through the bilinear form c(·, ·) (see (3.12)) and we set up a zero body force f = 0. The availability of an exact solution (cf. [48]) makes that this problem is frequently used as a benchmark. In Fig. 5.2 we illustrate the components of the displacement, the rotation and the pressure computed on a mesh consisting of 5120 triangular elements using the mixed FE method corresponding to k = 2, where the rectangular beam we consider has the following material properties: Young's modulus E = 1500, Poisson's ratio ν = 0.49, Lamé constants λ = 24664.4 and µ = 503.356, such that the model parameter equals η = 0.02. In addition, we conduct several tests for the lowest-order mixed FE and FVE methods on different mesh resolutions and report on the error with respect to the analytic solution (5.3).  In addition, although this is in general not true, we mention that the second order FE scheme ensures extremely rapid convergence (explained by the regularity of the true solution (5.3)). For ν = 0.4999, optimal convergence is recovered for finer meshes.
We also perform a series of tests for the lowest-order FE method using different Lamé constants and model parameters in order to test the performance of the methods when approaching the incompressibility limit, where we fix a Young's modulus E = 1500, vary the Poisson ratio ν, and use a mesh consisting of 100000 triangular elements and using 301201 D.o.f. Based on the comparisons in Table 5.4, we observe that the performance is barely modified for large values of λ.  Test 3 (3D beam bending). We also consider a three-dimensional beam problem. The beam occupies the domain Ω = (0, ℓ) × (0, w) × (0, w), with ℓ = 2.5, w = 0.5 (see a sketch in Fig. 5.4(a)); and its elastic properties are characterised by a Young modulus of E = 1000 and a Poisson ratio ν = 0.3, giving Lamé constants λ = 576.923, µ = 384.615, and the coefficient η = 0.4. The body force acts in the direction of gravityf = (0, 0, −ρg) t and it is specified by g = 9.8 and ρ = 0.2. Zero displacements are enforced on the face x = 0, whereas on the remainder of the boundary we consider zero normal stresses incorporated through the term ∫ x>0 2η(∇u − div u I)n · v defining the bilinear form c(·, ·) (see (3.12)). In Fig. 5.4 we illustrate (on the deformed configuration) the displacement, rotation vector, and pressure computed on a mesh of 45221 tetrahedral elements, employing a method of order k = 2. In the case of gravity-induced deflection, the Euler-Bernoulli beam theory predicts a maximum vertical deflection of δ = ρg Aℓ/(8E I ), occurring at the free end of the body, A = w 2 is the area of the cross-section, and I = A 4 /12 is the planar inertial moment. Table 5.5 compares the expected deflection with the vertical displacement measured on the midpoint of the face located at x = ℓ, for different discretisation choices. We also tabulate the norms of the approximate solutions generated with the lowest-order FE method on successively refined meshes (see Table 5.6).
Test 4 (Cook's membrane benchmark). We finalise the set of tests by considering a two-dimensional quadrilateral panel with domain Ω defined as the convex hull of the set {(0, 0), (ℓ, w), (ℓ, ℓ + s), (0, w)}, with ℓ = 48, w = 44, s = 16, and proceed to study its elastic response dominated by bending and shear. This benchmark is known as Cook's membrane problem (cf. [49]). The panel is clamped at the left edge (x = 0) and the body is subjected to a shearing distributed loadt = (0, 1/s) t on the opposite end (at x = ℓ and giving a resulting load of magnitude 1,   Fig. 5.5(a)). This effect is incorporated in the formulation through the term − ∫ x=ℓ t · v ds added to the functional F(·) in the modified weak formulation (3.10)-(3.11). A traction-free condition is applied on the non-vertical boundaries (imposed as in the previous test, using (3.12)), and we set up a zero volume force f = 0 (so that the weight of the membrane is not considered). The elastic plate has Young's modulus E = 1, Poisson ratio ν = 1/3, Lamé constants λ = 3/4, µ = 3/8, giving a scaling constant of η = 1/3. Fig. 5.5 portrays the displacement, rotation and pressure fields on the deformed domain (without amplification of the deformation field). We also conduct several  tests for different mesh resolutions and report on the vertical displacement (deflection) measured at the midpoint of the right end of the domain, (x 0 , y 0 ) = (ℓ, ℓ + s/2). The test results are shown in panel (b) of the figure, where the convergence behaviour of the deflection is observed as a function of the number of points discretising the right edge of the membrane. In the absence of a known closed form solution for this problem, we also include a referential value reported in the literature (according to [50][51][52], under plane stress conditions the maximum vertical displacement at this point should be around 23.92). To conclude we perform again the Cook's membrane test, but focusing on the nearly incompressibility limit. We choose the model parameters E = 250, ν = 0.4999, λ = 416611, µ = 83.3389, and η = 0.0002. As reference value for the maximum deflection at the point (x 0 , y 0 ) we consider 7.505 (see [4,53]), and conduct a convergence analysis portrayed in Fig. 5.5(c) (see also Table 5.7, where we display all individual norms for the numerical approximations via FE schemes of different orders). This time the vertical displacement is plotted against the D.o.f. associated to the underlying discretisation, where we also include a comparison against numerical results obtained with other finite element formulations applied to the original equations (2.1) (a classical pure-displacement formulation discretised with piecewise continuous elements of degree k, the Taylor-Hood finite element for a displacement-pressure formulation, the MINI-element [54], and a stabilised interior-penalty DG method [55]).
These schemes have comparable complexity (but we do not include other mixed methods based on stress or pseudostress formulations, as their associated cost would be much higher). A further comparison between these methods is presented focusing now on their computational cost and measured in terms of CPU time. To do so, we consider the simple test case of a square domain with clamped boundaries where the structured primal mesh has 10000 vertices. We set E = 10000, ν = 0.33 and solve the elastostatics problem using different methods whose performance is shown in Table 5.8. The tabulated results display measured wall CPU time comprising matrix assembly, factorisation, and solution. As direct solvers might not be preferable for large systems and 3D problems, we also include the wall time for the matrix inversion using a Krylov solver, the bi-conjugate gradient stabilised method (BiCGStab) preconditioned with an incomplete Cholesky factorisation. The results indicate that the proposed mixed and FVE methods are preferable. As mentioned above, the matrix system associated with the Schur complement formulation (5.2) is substantially smaller than in the other methods. Unfortunately, a drawback of this formulation over the standard implementation of our mixed FE and FVE schemes (3.1)-(3.2) and (4.2)-(4.3) is that the assembly of the blocks and computation of the action of the inverses (that we do not carry out in an optimal manner) consume most of the CPU time.
Test 5 (mixed boundary conditions in a 2D non-convex domain). For our last test, we investigate numerically the accuracy of the formulations proposed in Sections 3.2 and 4.2. Apart from setting displacement-traction boundary conditions, we again define the problem on a non-convex domain (the unit square from Test 1A now has a hole of  Table 5.9 for the lowest-order methods, where we have set the bottom, right, top and left walls of the domain as the displacement boundary and the non-homogeneous traction condition is imposed on the inner circle. Coarse mesh solutions for displacements are exemplified in Fig. 5.6. The FE scheme exhibits a similar behaviour to the one observed in Test 1B: suboptimal displacement convergence in the L 2 -norm, again due to the non-convexity of the domain, while the remaining errors behave as in the pure Dirichlet case. While we only confirm this behaviour numerically, these computations stand as a motivation to investigate further the theoretical properties of the formulations in the case of mixed boundary conditions.