Convergence of a cell-centered finite volume discretization for linear elasticity

We show convergence of a cell-centered finite volume discretization for linear elasticity. The discretization, termed the MPSA method, was recently proposed in the context of geological applications, where cell-centered variables are often preferred. Our analysis utilizes a hybrid variational formulation, which has previously been used to analyze finite volume discretizations for the scalar diffusion equation. The current analysis deviates significantly from previous in three respects. First, additional stabilization leads to a more complex saddle-point problem. Secondly, a discrete Korn's inequality has to be established for the global discretization. Finally, robustness with respect to the Poisson ratio is analyzed. The stability and convergence results presented herein provide the first rigorous justification of the applicability of cell-centered finite volume methods to problems in linear elasticity.

1. Introduction. We consider the following problem of isotropic (but heterogeneous) linear elasticity [16] ∇ · σ + f = 0 in Ω σ = 2µ∇u + λ(tr ∇u)I in Ω Here the domain Ω is a bounded connected polygonal subset of R d , with boundary ∂Ω = Γ D ∪ Γ N . We have introduced the symmetric gradient operator by the notion ∇u ≡ (∇u + ∇u T )/2. Furthermore, let the parameter functions f ∈ (L 2 ) d and the Lamé parameters 0 < µ ≤ µ(x) ≤μ and λ be bounded and positive, defined almost everywhere. If Γ D has positive measure, equations (1.1) have a unique weak solution in (H 1 (Ω)) d . Otherwise, if Γ N = ∂Ω, equations (1.1) have a unique weak solution in (H 1 (Ω)) d /R(Ω), where R is the space of rigid body motions: R(Ω) = {a + ω ∧ x; a, ω ∈ R d }. In the latter case Ω f dx = ∂Ω g N dS is a necessary compatibility condition on the data. Without loss of generality, we will assume, by subtracting any smooth function satisfying the boundary conditions and correspondingly modifying the right-hand side, that both g D = 0 and g N = 0. We note in particular that we do not consider transformation which is available for the (simpler) case of homogeneous Dirichlet boundary conditions, when equations (1.1) can be recast in a locally coercive form [16].
In the continuation, it will be convenient to refer directly to the weak form of equations (1.1). We will here, and in the following use the convention H 1 (Ω) ≡ (H 1 (Ω)) d , and tacitly assume that the space is restricted to the homogeneous Dirichlet boundary condition. The weak form of equations (1.1) then takes the form (see e.g. [16]): Find u ∈ H 1 (Ω) such that Recently, we proposed to extend cell-centered finite volume methods for the scalar diffusion equation to analyze finite volume methods for elasticity. Cell-centered finite volume methods may be advantageous for problems associated with poro-elastic materials in geological applications. In particular, it is advantageous to A) Exploit a similar grid and data structure between the fluid flow and mechanical discretizations [19,24], B) Share the same restrictions on non-matching grids and hanging nodes for both the flow and mechanical discretizations, and C) Have explicit force balance and traction at surfaces and grid faces (often associated with fractures and faults). We are therefore in particular interested in the generalization of the so-called Multi-Point Flux Approximation (MPFA)-O method for scalar equations [1]. This variant of finite volume method has proved popular in applications and is also amenable to theoretical analysis. Most notably, 2 nd -order convergence in potential and 1 st -order convergence in flux was established using a link to mixed finite element methods already in [18,27], while analysis of the method following the discrete functional framework was presented in [4]. The generalization of the MPFA-O method to linear elasticity was formulated for general anisotropic and heterogeneous problems in [23]. It is there termed the Multi-Point Stress Approximation (MPSA) method, and was supported by extensive numerical experiments indicating robust convergence results for a wide range of grids and Poisson ratios. It is the objective of this paper to provide a theoretical convergence analysis of this method by generalizing the discrete functional framework for finite volume methods.
The discrete functional framework for finite volume methods is detailed in [12]. This approach was utilized by [13] to develop finite volume discretizations for the scalar diffusion equation for which convergence could be proved under quite weak assumptions on the grid and coefficients. Furthermore, the framework was adapted to non-symmetrical discretizations in [4] and [3] to generalize and prove convergence of the MPFA methods.
The main obstacle in order to extend the analysis of discretizations for scalar diffusion equations to discretizations of equations (1.1) is to ensure coercivity of the discretization. In particular, the discrete functional spaces previously used for the scalar problem are conceptually similar to the Crouzeix-Raviart finite element space. This space does not satisfy Korn's inequality and thus lacks stability. Our work therefore extends the spaces to allow for a natural stabilization analogous to discontinuous Galerkin methods [14]. Furthermore, to account for the additional challenges associated with the lack of local coercivity for equation (1.1), we will additionally need to lean on ideas from variational multiscale [15] and discontinuous Galerkin [7] methods in our analysis. We will also address the issue of stability with respect to the Poisson ratio (so-called numerical locking) by reverting to ideas from mixed methods [9,19].
We note previous work on finite volume methods for elasticity. Most work where node-centered [6] and cell-centered [26] finite volume methods are introduced for elasticity contain only numerical validation. This includes also recent work on cellcentered methods [23,10]. When additional variables are introduced, convergence of finite volume methods has been established [20]. Similarly, convergence has recently been established for face-valued finite volume methods [19]. This latter citation is particularly important, as the method is furthermore shown to be locking-free. To the knowledge of this author, this contribution represents the first rigorous convergence proof for cell-centered finite volume methods on general grids and heterogeneities. Furthermore, we establish that the method is locking free for a large class of grids.
The manuscript is structured as follows. In section 2, we establish notation and recall the formulation and main results from the hybrid-variational framework. In section 3, we define our discrete mixed variational problem and establish its connection to the MPSA-O discretization. In section 4, we establish a local coercivity condition which guarantees the global coercivity of the discretization. This represents the key technical obstacle in order to extend previous work and prove stability of the discretization. In section 5, we largely exploit previous work on discrete functional discretizations to obtain the main convergence result. Section 6 provides detailed comments on the method, including application to homogeneous Dirichlet problems, reduced integration on simplex grids, and the corresponding scalar diffusion discretization. Section 7 details how the local coercivity condition simplifies to an explicit condition on the mesh for homogeneous problems, and addresses the issue of numerical locking. Section 8 concludes the paper.
2. Discrete functional framework. In this section we give the definition of our finite volume mesh and discrete variables.
2.1. Finite volume mesh. Following [3], we modify the construction of [12], and denote a finite volume mesh by the triplet D = (T , F , V), representing the mesh Tessellation, Faces, and Vertexes, such that: • T is a non-overlapping partition of the domain Ω. Furthermore, let m K denote the d-dimensional measure of K ∈ T . • F is a set of faces of the partitioning T . We consider only cases where elements σ ∈ F are subsets of d − 1 dimensional hyper-planes of R d , and all elements σ ∈ F we associate the d − 1 dimensional measure m σ . Naturally, the faces must be compatible with the mesh, such that for all K ∈ T there exists a subset F K ⊂ F such that ∂K = σ∈FK σ. • V is a set of vertexes of the partitioning T . Thus for any d faces σ i ∈ F , either their intersection is empty or i σ i = s ∈ V. Note that in the above (and throughout the manuscript), we abuse notation by referring to the object and the index by the same notation. E.g., we will by K ∈ T allow K to denote the index, as in F K , but also the actual subdomain of Ω, such that the expression ∂K is meaningful.
Additionally, we state the following useful subsets of the mesh triplet, which allows us to efficiently sum over neighboring cells, faces or vertexes: • For each cell K ∈ T , we denote the faces that comprise its boundary by F K and the vertexes of K by V K . We will associate with each corner s ∈ V K a subcell of K, identified by (K, s), with a volume m s K such that s∈VK m s K = m K . • For each face σ ∈ F , we denote the neighboring cells T σ and its corners for V σ . Note that for all internal faces T σ will contain exactly two elements, while it contains a single element when σ ⊂ ∂Ω. We will associate with each corner s ∈ V σ a subface of σ, identified by (s, σ), with an area m s σ such that s∈Vσ m s σ = m σ . • For each vertex s ∈ V, we denote the adjacent cells by T s and the adjacent faces by F s . We associate for each element K ∈ T a unique point (cell center) x K ∈ K such that K is star-shaped with respect to x K , and we denote the diameter of K by d K . Furthermore, we denote the distance between cell centers x K and x L as We associate with each face σ its outward normal vector with respect to the cell K ∈ T σ as n K,σ , and the Euclidian distance to the cell center d K,σ . For each subface (s, σ) we denote the subface center as x s σ and the smallest set of Gauss quadrature points sufficient for exact integration of second-order polynomials on (s, σ) as G s σ . For each quadrature point β ∈ G s σ we associate the position x β and weight ω β . We associate for each vertex s ∈ V its coordinate x s ∈ Ω. The above definition covers all 2D grids of interest. However, we place two restrictions on grids in 3D: Firstly, the above definition of a mesh requires all cell faces to be planar. The analysis that follows can be extended, at the cost of extra notation, to the case of non-planar faces can be allowed for as long as the faces allow for a piece-wise planar subdivision associated with the face corners V σ . Secondly, we will require that no more than three faces meet at a vertex. This permits quadrilaterals, simplexes, and all so-called 2.5-D grids (e.g. 2D horizontal grids extended vertically, as common in the petroleum industry). However, as an example we do not permit certain 3D grids such as pyramids. The formulation of the method readily generalizes to this case, however the application of Korn's inequality (section 4.3) becomes more technical.
Regularity assumptions on the discretization D are detailed elsewhere (see e.g. [12]), we will in the interest of simplicity of exposition henceforth assume that the classical grid regularity parameters (grid skewness, internal cell angles, and coordination number of vertexes) do not deteriorate.

Discrete variables and norms.
In contrast to MPFA methods for the scalar diffusion equation [1], it is for the MPSA O-method not sufficient to use socalled continuity points to provide sufficient constraints to yield a unique discretization [23]. Thus, we need to extend the discrete function spaces utilized previously [3,13] to allow for multiple unknowns per subface. We detail the three discrete spaces used in our analysis below.
The following discrete space is classical [12]: As with the dual interpretation of the elements K ∈ T , the space H T (Ω) is isomorphic to the space of discrete variables associated with the cell-center points x K . There should also be no cause for confusion in the following when we work with the vector-valued spaces, still denoted H T .
For the space H T we introduce the inner product and its induced semi-norm Here the operator γ σ u interpolates the piecewise constant values of H T onto the faces of the mesh, weighted by the distances d K,σ .
For Dirichlet boundary edges, σ ∈ Γ D , we take γ σ u = 0. Equivalently, the operator γ σ u can be defined as the value which minimizes the definition of the seminorm |u| T . We note that this semi-norm, and those that follow, are equivalent to full norms when Γ D has positive measure. In the case where Γ N = ∂Ω, rigid body motions must be excluded, as we will note later.
The following space is the discontinuous discrete space which we use to construct the consistent gradient functions of our scheme. It is to our knowledge novel in the context of analysis of finite volume methods, however it is natural in the sense that it is a discrete version of the first-order discontinuous Galerkin space (as we will emphasize later).
Definition 2.2. For the mesh triplet D, let H D be the set of real scalars (u K , u σ,β K,s ), for all K ∈ T , for all (s, σ) ∈ V K × F K and for all β ∈ G s σ . The space H D thus contains one unknown per cell, in addition to multiple unknowns on each interior sub-face. This will be essential to control the space R(Ω). As above, we will immediately take u σ,β K,s = 0 for all σ ∈ Γ D . We denote for all internal subfaces [[u]] σ,β s = u σ,β R,s − u σ,β L,s for u ∈ H D and T σ = {R, L} as the jump in the discrete function u across that edge. We will also need a notion of an average face value, and we denote similarly for all internal sub- As with |u| T , it is straight-forward to also define a proper norm for H D . The above "discontinuous" discrete space generalizes the "continuous" discrete space, which we recall as [3]: Definition 2.3. For the mesh triplet D, let H C be the set of real scalars (u K , u σ s ), for all K ∈ T and for all (s, σ) ∈ V K × F K .
By introducing the natural projection operator Π D : K,s = u σ s for all K ∈ T and for all (s, σ) ∈ V K × F K , we can immediately define the inner product In addition to the projection operator defined above, we shall need a few more operators to move between function spaces.
• Let the operator Π T : H D → H T be defined as (Π T u)(x) = u K for all x ∈ K and K ∈ T . Furthermore, as there should be no reason for confusion we also define Π T : for all K ∈ T and for all (s, σ) ∈ V K × F K . The spaces defined above satisfy the following inequalities.
• Discrete Sobolev inequality [12]: For all u ∈ H T and for all q ∈ [1, 2d/(d−2+ǫ)) Similarly, H T ,s and H C,s are defined through the projection operators defined above. The local spaces have the natural norms, which to be precise are given for all u ∈ H D as And for all u ∈ H T as Such that both 3. The MPSA Finite Volume discretization. In this section, we will utilize the spaces defined in Section 2 to establish a cell-centered finite volume method for elasticity. The method presented herein is a slight generalization of the MPSA Omethod as it was defined in [23]. The construction is inspired by, but generalizes in necessary aspects, the hybrid finite volume methods in [4,3].
3.1. Discrete mixed variational problem. Since we are dealing with the vector equation (1.1), we will seek solutions u in the discrete vector-valued spaces These spaces inherit all the definitions of their scalar counterparts, with the understanding that all inner products are extended with vector inner products in their definitions. The norms are in consequence the root of the square of the component-wise norms.
Following [4], we will use two notions of discrete gradients. However, since we will work with both spaces H C and H D , where the latter is multivalued on internal edges, the precise construction of the gradients differs from those works. The first gradient is the proper negative transpose of the divergence operator, and is constructed to yield a conservative finite volume formulation.
Definition 3.1. For each K ∈ T and each s ∈ V K we define the finite volume gradient for all u ∈ H C : Here, we denote the vector outer product with ⊗, such that each row of the product matrix contains the approximation to the gradient of the corresponding component of u. We construct a second gradient with the property that it is exact for linear variation in u with respect to the underlying physical space.
Definition 3.2. For each K ∈ T and each s ∈ V K we define the consistent gradient for all u ∈ H D : Here we extend the averaging notation in the natural way such that u σ K,s = 1 m σ s β∈G σ s ω β u σ,β K,s . In order to satisfy the desired consistency property, we require that (∇u) s K is exact for linear displacements, therefore g s K,σ are defined by the system of equations: Here I is the d-dimensional second-order identity tensor. For all 2D grids and for all 3D grids where no more than three faces meet at any vertex, equation (3.3) uniquely determines g s K,σ [1,3]. This encompasses all grids considered herein.
In the case of both gradients, the symmetric gradient is obtained in the same way as in the continuous setting by taking the average of the gradient and its transpose. Furthermore, the equivalent discrete divergence is the trace of the gradient, e.g.
We now define our finite volume scheme for linear elasticity, equations (1.1), through the specification of the following three bilinear forms. The first bilinear form embodies a discrete form of Hooke's law together with the finite volume structure of the method, and is analogous to the weak form stated in equation (1.2). Thus we The discrete coefficients are given as subcell averages of their continuous counterparts, e.g. µ K = m −1 K K µ dx. The second bilinear form controls jumps across subcell faces, and is defined for The family of weights α σ s are assumed to be uniformly bounded, 0 < α − ≤ α σ s ≤ α + < ∞. We retain the freedom to specify α σ s to improve the stability of the scheme. Numerically experiments indicate that the weights α σ s should be related to the harmonic mean of the material coefficients µ K and λ K [23].
Finally, we introduce a bilinear form to constrain the discrete unknowns on subfaces u σ,β K,s to the hyperplane given by u K and (∇u) s K ; thus for all (u, The above bilinear forms allow us to define the discrete mixed variational problem: The solution u D ∈ H D contains the solution satisfying the equations of elasticity on finite volume form (3.7), constrained to be piece-wise linear on sub-cells (3.8), while the remaining degrees of freedom are selected to minimize jumps in the solution (3.9). The component y C and y D are Lagrange multipliers for the constrained minimization problem, and will not be of further interest.
We note that equation (3.7) can be seen as a direct finite volume formulation of equations (1.1), wherein these equations hold in an integral sense for each cell K ∈ H T . Conversely, equation (3.7) can be identified as a Petrov-Galerkin discretization of equations (1.2), wherin the test functions are chosen as piece-wise constants on the cells K. In this interpretation, the shape functions are defined implicitly by equations (3.8) and (3.9). We return to the consistency of equations (3.7)-(3.9) in Section 5.

Finite volume formulation.
In this section we identify that the discrete mixed variational problem (3.7)-(3.9) is equivalent to a finite volume scheme. The forces acting on a subface T σ K,s are naturally defined from the bilinear form b D and Hooke's law (equation (1.1)b), thus we define for all for all u ∈ H D the tractions We verify by comparison to equation By considering v from the canonical basis of H C , we can now identify that the discrete variational mixed formulation (3.7)-(3.9) as equivalent to the hybrid finite volume We will see in Section 4 that due to the particular structure of the discrete differential operators, and the local definition of the forces in equation (3.10), the minimization problem in equation (3.16) has a unique solution in terms of the variables u T ∈ H T ⊂ H D , the set of variables associated with the cell centers.
Assume for the moment (we will return to this point in the next section) that for each s ∈ V we can define the local projection operator Π F V,s : H T ,s → H D,s , as the solution of Π F V,s u T ,s = arg min u∈U s a D (u, u) (3.17) Where U s ⊂ H D,s are the spaces that satisfy the constraints of (3.13)-(3.15) with u T ,s given. It follows that we can construct explicit, local expressions for the forces T σ K,s using the expression given in equation (3.10): The local coefficient tensors t K,K ′ ,s,σ are referred to as subface stress weight tensors, and generalize the notion of transmissibilities from the scalar diffusion equation [23]. We infer from equation (3.14) that t K,K ′ ,s,σ = −t K ′ ,K,s,σ . Furthermore, we have from equation (3.13), that also the face stress weight tensors can be calculated, with Combining equations (3.12) and (3.19), we arrive at the cell-centered finite volume scheme expressed in terms of cell-centered variables only. This scheme is identical to the scheme presented as MPSA O-method (general) in [23].
4. Local problems, Korn's inequality, and coercivity. Our goal is to show that the discrete mixed variational problem (3.7)-(3.9) is well-posed. This requires four steps. First, we formalize the discussion in section 3.2 using a variational multiscale framework to state variational problem only in terms of variables in the space H T , exploiting local operators which are defined through local problems. Secondly, we show the stability of these local problems. Thereafter, we arrive at a discrete Korn's inequality through a projection onto piece-wise linear space on each subcell. Finally, we establish that coercivity, and thus wellposedness, of the full problem can be verified based on local coercivity criteria.

4.1.
A non-mixed discrete variational formulation. We use an approach similar to the variational multiscale methods [15] as applied to mixed problems [5,22], in that we split the mixed problem (3.7)-(3.9) into two coupled problems. We introduce H F D ⊂ H D denoting the variables associated with cell faces, such that Identifying now H T as the space of coarse variables, and H F D as the space of fine variables, we thus consider the problem: Note that there is no integral term on the right hand sides of equations (4.2)-(4.4) since Π T {0 P , v} = 0. Furthermore only the fine-scale problem is on mixed form. As observed in section 3, the aim is to resolve the mixed fine problem locally, which in the present context is realized by interchanging sums in the definition of the operators (3.4)-(3.6) to observe that for χ ∈ {a, c} and for (u, v, w) ∈ H D ×H C ×H D where the local bilinear forms are defined as and a D,s (u, w) = Existence and uniqueness of solutions to (4.9)-(4.11) is treated in Lemma 4.1 below. We note that the local projections are linear, and construct the global finite volume projection Π F V : Using the finite-volume projection, we establish the decoupled coarse problem: Equations (4.9)-(4.13) are algebraically equivalent to the original formulation given in equations (3.7)-(3.8). However, the present formulation has the advantage that wellposedness can be established in steps. In particular, the mixed problems (4.9)-(4.11) are all local, so that we avoid considering a global inf-sup type condition. Furthermore, the global problem is similar in structure to a standard Galerkin approximation to the original system (1.1), and for this problem we will establish (local) conditions to ensure coercivity.

4.2.
Well-posedness of local mixed problems. The local mixed problems are given by equations (4.9)-(4.11), and define the finite volume projection.
Lemma 4.1. The solutions of the local mixed problems (4.9)-(4.11) exist and are unique.
Proof. Since a D,s is a symmetric, positive semi-definite bilinear form, the system is equivalent to a constrained minimization problem, to which existence of solutions is guaranteed. In particular, for u T = 0, it is clear that Π F V,s u T = 0 satisfies the constraints (4.9)-(4.10) and has the minimum energy a D,s (Π F V,s 0 T , Π F V,s 0 T ) = 0. Thus for u T = 0, the space of minimizers is isomorphic to the space of continuous piecewise linear functions on T s . However, all cells K ∈ T are star-shaped with respect to x K , thus ( x σ K,s −x K )·n σ,s > 0, from which it follows that the constraints (4.9) are linearly independent to the null-space of a D,s ({0 T , u}, {0 T , u}). Uniqueness follows since it is straight-forward to verify that the null-space of a D,s ({0 T , u}, {0 T , u}) has at most d 2 degrees of freedom, while at least d 2 constraints (4.9) are independent.
Since the local mixed problem (4.9)-(4.11) is finite-dimensional, the norm of the projection operator is finite and we define: Definition 4.2. For every s ∈ V, the local mixed problem (4.9)-(4.11) has a unique solution by lemma 4.1, and we define the norm of the of the solution operator Π F V,s as θ 1,s , such that |Π F V,s u| D,s ≤ θ 1,s |u| T ,s for all u ∈ H T (4.14) The coefficient θ 1,s will in general be dependent on the geometry of the mesh D and parameter functions µ and λ, but independent of the mesh size h due to the scaling invariance of the norm. We define the maximum over the local norms as θ 1 = max s∈V θ 1,s .

4.3.
Korn's inequality. We will need a discrete Korn's inequality to show coercivity of the method. Since the constraints (4.10) guarantee that the projections Π F V,s are consistent with a piecewise linear approximation on the subcells, we are inspired to consider Korn's inequality in this setting: Definition 4.3. Let the projection operator Π L 2 : H D → (L 2 (Ω)) d be defined such that for all (K, s) ∈ T × V and x ∈ K, we have We note that this projection is particularly appropriate for functions in u ∈ Π F V,s H T , since equations (4.7) and (4.10) assure that in this space the face variables satisfy u σ,β K,s = Π L 2 u(x σ,β K,s ) for all cells K ∈ T , faces σ ∈ F K , corners s ∈ V K ∩ V σ and quadrature points β ∈ G s σ . For the space of piecewise linear functions we may thus recall the appropriate Korn's inequality [7] Lemma 4.4.
Proof. The left-hand side is identical to the broken H 1 semi norm of Π L 2 u, while the jump term bounds also the surface L 2 norm of jumps internal to cells due to mesh regularity and continuity at the points x K for all cells K ∈ T . Lemma 4.4 therefore follows from equation (1.18) in [7]. 4.4. Local coercivity conditions. We will derive local coercivity conditions which guarantee the coercivity of the bilinear form on the reduced subspace as given by equation (4.13). The reduced subspace is critical, as b D,s is not coercive on the full space H D (even without the introduction of a full norm). As a consequence, coercivity will depend on the local finite volume projections Π F V,s , and we state the following assumption on the solution of the local mixed problems. Where the local energy semi-norm is associated with the symmetrized bilinear form This assumption can be verified locally while assembling the discretization, and moreover it can be verified a priori for certain classes of meshes, see Section 7. We recall two useful stability estimates. Lemma 4.6. a) There exists a constant c D , only dependent on the regularity of D, such that for all u ∈ H D Proof. The first inequality follows readily from the definition of the semi-norm and the consistent gradient (∇u) s K . The second inequality is shown by contradiction. Let the right-hand side be zero. Then (∇u) s K = 0, and due to (3.8) the function u thus lies in the space of rigid body motions on each subcell (K, s). However, since the second-order Gauss quadrature is exact for linear functions, u is in the space R(Ω). But then also the left-hand side of the inequality is zero. The scaling of the norms then guarantees that the constant c L 2 is independent of h.
We are now prepared to state the global coercivity result for the method Theorem 4.7. For given parameter fields µ, λ, and mesh D, let assumption 4.5 hold. Then the coarse variational problem (4.13) is coercive in the sense that it satisfies if a) Γ D is measurable then for all u ∈ H T and b) if Γ N = ∂Ω then for all Where the constant Θ is dependent on the mesh triplet D but does not scale with h, and is bounded below by: Proof. By the definition of the local bilinear forms, Assumption 4.5 and the lower bounds on the coefficient functions µ and λ we have for all We denote c µ = min(2µ, 1 − ξ), for any 0 < ξ < 1, and now due to Lemma 4.4 and 4.6a) we obtain Utilizing Lemma 4.6b) together with the basic norm inequalities stated in Section 2.1 completes the proof when ξ is chosen to maximize Θ.

Convergence of the method.
Convergence of the scheme as given in Section 3 is now established for all grid sequences satisfying a uniform coercivity bound. In particular, weak convergence to some function u ∈ (H 1 ) d of the discrete solution Π T u D and its gradient ∇ D u D (defined below), follows immediately from previous work [3]. The strong convergence of the gradient ∇ D u D , and establishing that u is a weak solution of Equations (1.1), requires invoking the finite-volume projection Π F V due to the discontinuous nature of the discrete space H D and the mixed formulation of the saddle-point problems used in the current work. We structure this section accordingly.
Definitions 5.1. We consider the following continuous extensions of the cellaverage finite volume gradients for discrete functions u ∈ H D : (5.1)

Furthermore, we consider the continuous extension of the consistent gradient
Note that from its definition, the discrete gradient satisfies the stability estimate Where g 0 = max K∈T ,σ∈FK ,s∈Vσ |g s K,σ |/d K,σ , which depends on the regularity of the grid but not on h.
We now recall the following result: Lemma 5.2. Let D n be a family of regular discretization triplets (in the sense that mesh parameters remain bounded) such that h n → 0, as n → ∞. Furthermore, let θ 1 and Θ be bounded independently of n. Then for all n, the solution u n of equations (3.7)-(3.8) exist and are unique, there exists u ∈ (H 1 (Ω)) d , and up to a subsequence (still denoted by n) Π T u n → u converges in (L q (Ω)) d , for q ∈ [1, 2d/(d − 2 + ǫ)) as h n → 0. Finally, the cell-average finite volume gradient ∇ D u n converges weakly to ∇ u in (L 2 (Ω)) d 2 .
Proof. The proof follows immediately from the coercivity of the scheme (Section 4) and the compactness arguments detailed in [13] and [3].
We also state the strong convergence of the consistent gradient ∇ D u(x). This result was achieved for the scalar diffusion equation in [3], where the discrete solution lies in H C . We can extend their calculation to the present case as summarized below. Lemma 5.3. Consider the same case as in Lemma 5.2. Then the consistent gradient ∇ D u n converges strongly to ∇ u in (L 2 (Ω)) d 2 .
Proof. We need to show that Introducing a function ϕ ∈ (C ∞ (Ω)) d which approximates u, and using the identity we can bound the integral in (5.3) by The second term on the right-hand side converges since the projection operators are exact for linear functions, while the last term vanishes for n → ∞ [3], so it suffices to deal with the first right-hand side term in equation (5.4). However, since u n = {Π T u D,n , Π F V Π T u D,n }, the first term is bounded by the bilinear form b D due to coercivity. A straight-forward calculation exploiting that ϕ approximates u then verifies that all terms on the right-hand side of equation (5.4) converge to zero.
The preceding definitions and lemmas allow us to state the main convergence result for the MPSA method.
Theorem 5.4. Consider the same case as in Lemma 5.2. Then the limit u ∈ (H 1 (Ω)) d of the discrete mixed variational problem (3.7)-(3.9), and consequently the MPSA O-method, is the unique weak solution of the problem (1.1).
Proof. Lemmas 5.2 and 5.3 establish that the limit u ∈ (H 1 (Ω)) d exists, and the appropriate notion of convergence of the discrete gradients. It remains to show that u is a weak solution of problem (1.1). Uniqueness then follows from the uniqueness of weak solutions to (1.1). Recall again that u n = Π F V u T ,n . Then due to the stability of the projections and the strong and weak convergence of ∇ D u n and ∇ D u n , respectively, it follows that for all u, v ∈ (C ∞ (Ω)) d However, since (5.5) and (5.6) are the terms weak form of equation (1.1), as stated in equations (1.2), and since C ∞ is dense in H 1 , it follows that the limit u is the weak solution to (1.1).
6. Comments on the method. In this section we present various comments on the developments of sections 2-6. In particular, we comment on 1) the application to pure Dirichlet problems with homogeneous coefficients, 2) Reduced quadrature on simplex grids, 3) the corresponding finite volume for the scalar diffusion equation, and 4) a related finite difference method.
6.1. Homogeneous Dirichlet problems. For the special case where Γ D = ∂Ω, and the Lamé coefficients µ and λ are constant on Ω, it is well known that equations (1.1) can be re-written by integration-by-parts as [16]: This form is locally coercive, since the symmetrized gradient does not appear. We may proceed to discretize equation (6.1) as in the preceding sections, but with the bilinear form defined in equation (3.4) The resulting numerical method is also locally coercive and the Korn's inequality (section 4.3) is not needed to show coercivity of the method. While the local coercivity assumption 4.5 is still needed, this locally coercive formulation will due to the absence of Korn's inequality have a more favorable global coercivity constant.
6.2. Reduced integration on simplex grids. It is possible to consider using lower-order quadrature by choosing a smaller set of Gauss points G s σ . In particular, choosing a single Gauss quadrature point leads to the MPSA vector analog of the classical MPFA O(η) methods for the scalar equation, where the η is a parameterization of the location of the Gauss point [1,2]. In this setting the quadrature point will act as a point of strong continuity over the faces σ ∈ F between the (linear) solution in the two adjacent subcells T σ . On general grids, this does not lead to a well-posed discretization, since the local mixed system (4.7)-(4.8) is no longer wellposed: There exists for this case a non-trivial element of the kernel of a D,s which satisfies the constraints given by b D,s and c D,s .
However, for simplex grids with strictly acute corners, numerical experiments indicate that the local mixed system remains well-posed [23]. In this case, it is particularly attractive to consider the MPSA O(1/3) method. This parameterization implies that for all (K, s) ∈ T , V K the Gauss quadrature points G s σ are chosen such that the points x K , x β (for all β = G s σ with σ ∈ F K ∩F s ) and the location of the vertex s for a parallelogram in 2D and a parallelepiped in 3D. In this case, a straight-forward calculation shows that the two discrete gradients coincide, since Consequently, the bilinear form b D is symmetric [17,11,3].
6.3. Corresponding discretization for the scalar diffusion equation. The discretization analyzed herein can be directly applied to the scalar diffusion equation, and represents a new method in that context.
For the scalar diffusion equation, the analysis presented herein simplifies somewhat. In particular, the use of Korn's inequality can be omitted (as was also the case in Section 6.1), since the scalar bilinear form is locally coercive. Thus the local coercivity conditions can be stated directly, and the global coercivity of the method can be obtained by simple summation. The convergence of the method then follows by identical arguments to those used in Section 5.
6.4. A related finite difference method. The local coercivity Assumption 4.5 can be avoided by considering the related finite difference method. Indeed, consider the symmetric bilinear form Then we can formulate a finite difference method by considering the discrete mixed variational problem: together with equations (3.8) and (3.9) hold. This finite difference method clearly benefits from all the results given in Sections 4 and 5, however without the requirement of a local coercivity condition. Furthermore, the method provides a symmetric discretization, however these benefits come at the cost of a loss of exact force balance on each grid cell K ∈ T .
7. Local coercivity assumption for special grids. The key assumption in the proof is that there exists a class of grids satisfying Assumption 4.5. In this section, we will verify that such grids exist. To be precise, we give a sufficient condition on the sub-cell geometry to guarantee coercivity. This will allow us to establish a priori coercivity on regular triangulations of hexagons, squares, and equilateral parallelograms. Furthermore, we verify unconditional coercivity of the reduced integration proposed for simplex grids in section 6.2. Finally, we discuss the role of the Poisson ratio and locking in Section 7.3.

Vertex-symmetric meshes.
We introduce the following notion, valid for 2D grids.
Definition 7.1. We refer to a subcell (s, K) ∈ V ×T s as vertex-symmetric if it is symmetric with respect to the line through the points (x K , x s ). Similarly, we refer to a mesh as vertex-symmetric if all subcells in the mesh triplet are vertex-symmetric.
Vertex-symmetric meshes include the regular triangulations of hexagons, squares, and equilateral parallelograms.
We first state a preliminary lemma which resembles the coercivity assumption: Proof. Since, due to (4.10), the space Π F V,s H T is locally linear on subcells, and the consistent gradient is by Definition 3.2 exact on linear functions, Lemma 7.2 therefore reduces to showing that for all matrices P , such that P = 1 and for all (K, s) ∈ T , V K , it holds that µ K (P + P T ) : ( ∇(P x)) s K + λ K tr P ( ∇ · (P x)) s K ≥ θ 2,s µ K (P + P T ) : (P + P T ) + λ K (tr P ) 2 (7.1) We consider first the case λ K = 0, and consider the contradiction (P + P T ) : ( ∇(P x)) s K ≤ 0 (7.2) Using the Definition 3.1 we obtain Identifying the sum We re-write inequality (7.2) as (P + P T ) : [(S + S T )(P + P T )] ≤ 0 (7.5) But due the assumption of vertex-symmetry, it is easy to compute that (S + S T ) has strictly positive eigenvalues and thus equation (7.5) cannot hold, and the contradiction is thus shown to be false. The opposite case with µ K = 0 is analogous, however here the contradiction follows since for vertex-symmetric subcells for all u ∈ Π F V,s H T /R(Ω), and some finite constants θ ′ 2,s . We remark that corollary 7.3 states that the coercivity assumption is equivalent to verifying that the jumps in the elements of Π F V,s H T are bounded by the mechanical energy. This is trivially verified for hexagons, since the discretization is exact for deformations with constant gradients. A straightforward calculation also verifies the property for squares and parallelograms. However, we note that this property does not hold on equilateral triangles when the full set of quadrature points G σ s is used. 7.2. Simplex grids. As a consequence of section 7.1, we consider the reduced integration proposed in Section 6.2 for simplex grids.
Lemma 7.4. For simplex grids the local coercivity Assumption 4.5 with the reduced integration of Section 6.2 holds.
Proof. Due to equation (6.3) the discretization is symmetric, and the local coercivity Assumption 4.5 simplifies to showing that the jump terms are bounded, e.g. all However, the reduced integration imposes strong continuity at the points x β with β = G s σ , and the right-hand side of equation (7.6) is zero for all u ∈ Π F V,s H T . The above shows that the local coercivity Assumption 4.5 always holds for reduced integration on triangular grids. A similar situation was observed for the scalar MPFA method in [3]. However, this does not imply that the discretization is unconditionally convergent on triangular grids, since as pointed out in Section 6.2, the local problems (4.6)-(4.7) are not always well-posed. 7.3. Robustness with respect to Poisson ratio. It is of interest to understand the approximation quality of the method with respect to the incompressible limit. This is known to be a challenge for numerical methods, and is seen for e.g. lowest-order conforming finite element methods on simplexes [8].
The standard approach to understanding the issue locking is to recognize that in the limit λ → ∞, the solution must satisfy the equations (1.1) with the parameter choice λ = 0, subject to the constraint that This holds true both for the continuous and numerical solution. The phenomenon of locking occurs when equation (7.7) introduces more constraints in the discrete system then the available degrees of freedom [21].
For the methods discussed herein, we see that in the limit of λ → 0, each of the local problems (4.9)-(4.11) are constrained to satisfy equation (7.7). Thus if the number of vertexes in V, exceeds or equals the number of degrees of freedom of the global system, given by d times the number of elements in T , the finite volume discretization will lock. I.e., a locking phenomena will arise if card(V) > card(T ). This is the case for e.g. hexagons, where numerical locking has also been observed numerically [23].
To prove that locking will not occur for simplex grids, we need to establish the existence of a Fortin operator for the method [9]. In essence, we need to establish the existence of an operator with the following property.
Definition 7.5. We refer to an operator Π F : (H 1 ) d → Π F V,s H T as a Fortin operator if the following properties hold: The existence of a Fortin operator assures the robustness of the method with respect to locking, and we recall the following lemma: Proof. The proof uses classical arguments found in e.g. [9,14,19]. Lemma 7.6 guarantees that the approximation qualities of the space Π λ F V,s H T does not degenerate as λ → ∞. Our task remains to show the conditions under which a Fortin operator exists. To this end we will use ideas from linear solvers, where static condensation and multiscale finite volume methods provide the right partitioning of the grid.
The finite volume methods of the type discussed herein share the property that the discretization stencil is local, wherein "local" implies that the discrete conservation law for a cell K ∈ T depends only on cells L ∈ T s where s ∈ V K . This allows for a partitioning the tesselation T into a macromesh triplet D = {T, F, B} such that each K ∈ T belongs to either a macrocell k ∈ T, a macroedge z ∈ F or a makrovertex s ∈ B. We refer to e.g. [25] for the details of the macro-topology, but note the important property that any two macrocells k, l ∈ T are completely separated from the perspective of the discretization, such that static condensation can be performed solely dependent on the unknowns in the macroedges and macrovertexes. Furthermore, for all vertex s ∈ V it holds that at least one K ∈ T s lies in a macrocell k ∈ T.
Definition 7.7. A family D n of mesh triplets is locally underconstrained, if for each n there exists a macromesh triplet D n satisfying the following properties: 1. The maximum number of cells in any k ∈ T is bounded independent of n. 2. For each k ∈ T it holds that d card(T k ) > card(V k ). Note that due to Definition 7.7a) the macrocell diameter scales proportionally to the mesh diameter h. Regular family of meshes such as triangulations satisfy Definition 7.7.
Proof. Let the operator Π F : (H 1 ) d → Π F V,s H T be any operator such that (Π F u) K = u(x K ) for all K ∈ {F, B} (7.8) And such that (∇ · u) s K = 0 for all K ∈ T , s ∈ V K (7.9) Due to Definition 7.7b) such operators exist, and due to Definition 7.7a) and scaling, any such operator defined consistently across the grid sequence also satisfies Definition 7.7b). Thus Π F is a Fortin operator, and the numerical method is robust.
Regular Cartesian lattices with Cartesian macrocell provides a limiting case where d card(T k ) = card(V k ) − 1 independent of the coarsening ratio, and as such never satisfy Definition 7.7. Nevertheless, Cartesian grids appear to be locking free based on numerical investigation [23].
8. Concluding remarks. We have expanded the hybrid finite volume framework and consider a mixed discretization in the context of linear elasticity. We use ideas from the variational multi-scale method to re-write the discretization in a subspace of cell-centered variables. By exploiting a discrete Korn's inequality borrowed from the analysis of discontinuous Galerkin methods we obtain global coercivity of the method, dependent only on a locally computable coercivity condition. Convergence is then obtained by compactness and consistency arguments. The discretization is designed to be identical to the MPSA finite volume method recently proposed, and convergence of that method is therefore established by the results herein.
The local coercivity conditions are simplified and verified a priori in the setting of vertex-symmetric meshes in 2D and simplex meshes in 2D and 3D. We furthermore establish necessary conditions for the robustness of the scheme with respect to numerical locking using tools from mixed methods and static condensation. Finally, we also identify a (new) finite difference method based on the same framework, but for which no local coercivity assumption is needed.
The analysis presented herein is fully consistent with the numerical results presented in [23] on similar classes of grids, and provides a rigorous understanding of the method.
The current work has not considered the convergence rate of the method, and such results would necessarily require additional assumptions on the regularity of the physical parameters µ and λ and on the mesh sequence D n .