{'fraction_non_alphanumeric': 0.04867521879870519, 'fraction_numerical': 0.013067977460736123, 'mean_word_length': 4.1875, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A three-dimensional Hellinger-Reissner Virtual Element Method for linear elasticity problems

We present a Virtual Element Method for the 3D linear elasticity problems, based on Hellinger-Reissner variational principle. In the framework of the small strain theory, we propose a low-order scheme with a-priori symmetric stresses and continuous tractions across element interfaces. A convergence and stability analysis is developed and we confirm the theoretical predictions via some numerical tests.


Introduction
The Virtual Element Method (VEM), introduced in [11,13], is a recent technology for the approximation of partial differential equation problems. This method is a generalization of the Finite Element Method (FEM) which allows to deal with arbitrary polygonal/polyhedral meshes, also including non convex and distored elements. To garantee this flexibility, the virtual element method abandons the idea of the local polynomial approximation, typical of FEM, to use approximating functions which are solution of suitable local PDE. In general, these non-polynomial functions are not explicitly known. Therefore, the main idea of this method is to exploit the available information (the degrees of freedom) to compute the stiffness matrix and the righthand side of the discretized problem.
In the present paper we extend the study presented in [7] to the three dimensional case. More precisely, we design and analyze a low-order virtual element method for linear elasticity problems. Within the framework of small displacements and small deformations, we consider the Hellinger-Reissner variational principle as the basis of our discretization procedure. This mixed formulation describes the problem by means of both the displacement and the stress fields. In the Finite Element practice it is a difficult task designing a stable and accurate scheme that preserves both the symmetry of the stress tensor and the continuity of the tractions at the inter-elements, see for instance [20] and [4]. The fundamental reason behind this difficulty lies in the local polynomial approximation, which forces the introduction of nodal degrees of freedom for the stress unknown. The resulting finite element schemes are typically quite cumbersome, especially in the three dimensional case, see [2]. Furthermore, the presence of nodal degrees of freedom introduces an additional complication if one aims at using the hybridization procedure to solve the discrete linear system, see [3]. We exploit the flexibility of virtual element methods to avoid these drawbacks and to develop an optimal scheme which is reasonably cheap with respect to the delivered accuracy. Since the approximated stresses does not have nodal degrees of freedom (on the contrary, the degress of freedom are entirely local to each polyhedron face), the hybridization procedure could be easily applied to our VEM scheme. This aspects show that, even for tetrahedral or hexahedral meshes, the proposed VEM method is a valid alternatives to FEM schemes.
The paper is organized as follows. In Section 2 we briefly introduce the classical Hellinger-Reissner formulation of the 3D elasticity problem. Section 3 describes the Virtual Element approximation we propose, while Section 4 is about the convergence analysis of the method. The numerical experiments, which confirm the theoretical predictions, are detailed in Section 5. Finally, we draw some conclusions.
Space notation. Throughout the paper, we will make use of standard notations regarding Sobolev spaces, norms and seminorms (cf. [35] for example). In addition, given two quantities a and b, we write a b when there exists a constant C, indipendent of the meshsize, such that a ≤ C b. Finally, given any subset A ⊂ R n and an integer k ≥ 0, P k (A) denotes the space of polynomials up to degree k, defined on A; whereas, given a functional space X, we denote by [X] 3×3 s the 3 × 3 symmetric tensors whose components belong to the space X.
Mesh notation. Given a polyhedron E with n E f faces we denote its volume, diameter and barycenter by |E|, h E and x E , respectively. In a similar way we refer to the area, diameter and barycenter of a face f , while |e| denotes the length of the edge e. Given a polygonal face f , we use x andx to indicate the global and local coordinates of a generic point of f , respectively.
As usual, H(div; Ω) is the space of tensor in L 2 (Ω) 3×3 whose divergence is the vectorvalued operator in L 2 (Ω) 3 . The elasticity fourth-order symmetric tensor D is assumed to be uniformly bounded and positive-definite. It is well known that Problem (1) is well posed, cf. [20] for example. In particular, considering the natural norms where C is a constant depending on Ω and on the material tensor D, which, however, does not degenerate in the incompressible limit. We also remark that a possible interesting variant of the variational formulation (1) has been recently proposed and studied in [38] and [39]. When considering a polyhedral mesh T h of the domain Ω, the bilinear form a(·, ·) is split as

The Virtual Element Method
In this section we define our Virtual Element discretization of Problem (1). Let {T h } h be a sequence of decompositions of Ω into general polyhedral elements E with We suppose that for all h, each element E in T h fulfils the following standard assumptions (cf. [11]): • (A1) E is star-shaped with respect to a ball of radius ≥ γ h E , • (A2) for every face f ∈ ∂E we have h f ≥ γh E and f is star-shaped with respect to a disk of radius ≥ γ h f , where γ is a suitable positive constant. We also assume that the elasticity tensor D is piecewise constant with respect to the decomposition T h , i.e. D is constant on each polyhedron of the mesh [7].

The local spaces
To describe the local spaces employed in the VEM proposed scheme, we introduce two spaces: RM (E) and T h (f ).
Space RM (E). It is the space of local infinitesimal rigid body motions: whose dimension is 6. • The three dimensional tangent vector t is determined by a linear combination of two given linearly independent tangential vectors t 1 and t 2 , i.e.
• The rotational term a n f ∧ (x(x) − x f ) is determined by a single scalar value a ∈ R.
• The polynomial p 1 (x) ∈ P 1 (f ) is a two variable polynomial with respect to the local face coordinate system so it is determined by three parameters, for instance: Therefore, this space consists of vector functions whose tangential component is a 2D face rigid body motion (the first two terms of (3)), while the normal component is a linear two-variable polynomial (the last term of (3)).

Stress space.
We are now ready to introduce our local approximation space for the stress field: We have the following Proposition.
it holds allows to compute div τ τ τ h using the degrees of freedom of the space T h (f ). Testing (7) with constant functions r(x) = α α α and recalling (5), we have Hence, we obtain We then infer that ω i.e. ω ω ω E solves system (6). We also notice that from (8) we deduce that the linear operator is symmetric and positive definite.
From Proposition 3.1 and (3) we infer that the dimension of the space (4) is Displacement space. The local approximation space for the displacement field is defined by, see (2): and it follows that dim(U h (E)) = 6.

The local forms
In this section we introduce the VEM counterparts of the local forms associated with the continuous problem.
The local mixed term. Given E ∈ T h , we begin by noticing that the term For this reason, we do not need to introduce any approximation of the terms (div τ τ τ , u) and (div σ σ σ, v) in problem (1).
The local bilinear form a E (·, ·). The local bilinear form As it is standard in the VEM procedure (cf. [11]), we build a computable approximation of the bilinear form by defining a suitable projection operator onto local polynomial functions. In our case, we introduce This is a projection operator onto the constant symmetric tensor functions and it is computable. Indeed, we notice that each π π π 0 ∈ [P 0 (E)] 3×3 s can be written as the symmetric gradient of a linear vectorial function, i.e. π π π 0 = ε ε ε(p 1 ), with p 1 ∈ [P 1 (E)] 3 . Hence, using the divergence theorem, the right-hand side of (10) becomes which is clearly computable (see also Proposition 3.1). Then, the approximation of a E (·, ·) reads: where s E (·, ·) is a suitable stabilization term. In this paper we propose the following choice: Above, κ E is a positive constant to be chosen according to D. For instance, in the numerical examples of Section 5, κ E is set equal to 1 2 tr(D |E ). However, any norm of D |E can be used. A possible variant of (12) is The local loading term. We split the load term on each element and we have Since v h ∈ RM (E), the right-hand side is computable via quadrature rules for polyhedral domains.

Remark 1.
Since this integral involves a sufficiently regular function f , to get a "good" approximation we exploit a quadrature rule of high degree. For the numerical examples in Section 5, we use a quadrature rule of degree 4.

The discrete scheme
Starting from the local spaces and local terms introduced in the previous sections, we can set the global problem. More specifically, we introduce a global approximation space for the stress field, by glueing the local approximation spaces, see (4): For the global approximation of the displacement field, we take, see (9): Then, given a local approximation of a E (·, ·), see (11), we set The method we consider is then defined by

Stability and convergence analysis
Since some results of the analysis follows the guidelines of the theory developed in [7], in this section we do not provide full details of all the proofs. First, for all E ∈ T h , we introduce the space: The global space Σ is defined as In the sequel, given a measurable subset A ⊆ Ω and r > 2, we will use the following space equipped with the obvious norm.

An interpolation operator for stresses
We introduce the local interpolation operator where: We remark that if τ τ τ is not sufficiently regular, the integral in the right-hand side of (17) is intended as a duality between W − 1 r ,r (∂E) 3 and W 1 r ,r (∂E) 3 . Instead, if τ τ τ is a regular function, the above condition is equivalent to require: We now show that I E τ τ τ ∈ Σ h (E) is well-defined by conditions (17). Indeed this is an immediate consequence of the following Lemma.
Proof. Recalling (4), Proposition 3.1 and (18), it is sufficient to prove that, given a face imply (τ τ τ h n) |f = 0. To this end, we first set (cf. (3)) with t f an arbitrary constant vector tangent to the face f , a ∈ R, n f the outward normal vector and Choosing δ δ δ f = 0 and γ γ γ f arbitrary in (20), and considering (21), we infer t f = 0 and c 1 = 0. Hence, it holds: We now select γ γ γ f = 0 and δ δ δ f = n f in (20). From (22) we get We have: Above, we have used that n f ∧ (x f − x E ) is a constant vector and that n f ∧ (x(x) − x f ) has zero mean value over the face f , to infer that: From (23) we deduce a = 0, and therefore we get: We finally select γ γ γ f = 0 and δ δ δ f = t f in (20), with t f an arbitrary vector tangential to the face f . Thus, we obtain: We now choose t f such that the (tangential to the face) vector n f ∧ t f has components (c 2 , c 3 ) with respect to the local coordinate system (x,ỹ). Then, from (25) we get which implies c 2 = c 3 = 0. Therefore, (τ τ τ h n) |f = 0, see (24), and the proof is complete.
The global interpolation operator I h : W r (Ω) → Σ h is then defined by glueing the local contributions provided by I E . More precisely, for every τ τ τ ∈ W r (Ω) and E ∈ T h , we set Using the same steps detailed in [7], we have the following error estimates for the interpolation operator I h .

The ellipticity-on-the-kernel and inf-sup conditions
The proposed approach satisfies the compatibility conditions. First, we notice that (see (13), (4) and (14), (9)): Then, we introduce the discrete kernel K h ⊆ Σ h : and we infer from (26) that τ τ τ h ∈ K h implies div τ τ τ h = 0. Hence, it holds: This is essentially the property that leads to the following ellipticity-on-the-kernel condition.
For the discrete inf-sup condition, we need the following commuting diagram property.

Proposition 4.4. For the operator
where Π RM denotes the L 2 -projection operator onto the space of the rigid body motions, see (2).
Proof. It is sufficient to prove property (29) locally, in each element E ∈ T h . Fix now r ∈ RM (E) and τ τ τ ∈ W r (E). We have: From (30) and the definition of L 2 -projection operator, we get div( Using Proposition 4.4 the following discrete inf-sup condition follows from the theory developed in [7].

Error estimates
We denote with P 0 (T h ) the space of piecewise constant functions with respect to the given mesh T h . Using the techniques developed in [7], one can prove the following result.

Numerical results
In this section we numerically assess the proposed VEM approach through the study of the method accuracy on a selected number of test problems. The numerical results confirm the proved theoretical results. We consider the standard unit cube Ω = [0, 1] 3 as the domain of our problems and we take the following three types of mesh: • Cube, a mesh composed by standard structured cubes; • Tetra, a Delaunay tetrahedralization of the domain Ω; • CVT, a Voronoi tassellation obtained by the Lloyd algorithm [31]; • Random, a Voronoi tassellation achieved with random control points.
We remark that the meshes CVT and Random are very challenging. Indeed, they could have some elements with small faces and edges, and we remark that such case is not covered by the developed theory, i.e. assumptions (A1), (A2) and (A3). However, the numerical results show that the proposed methods are fairly robust with respect to this geometric situation. These two type of meshes are build via the voro++ library [40]. Moreover, the whole numerical scheme is developed inside the vem++ library, a c++ code realized at the Univeristy Milano -Bicocca during the CAVE project (https://sites.google.com/view/vembic/home).
In order to assess the convergence rate, for each type of mesh, we define the following mesh-size h: where we recall that N E is the number of elements in the mesh, and h E is the diameter of the polyhedron E. The accuracy and the convergence rate assessment is carried out using the following error norms: • L 2 error norm for the displacement field: • L 2 error on the divergence: • L 2 error on the projection:

Cube
Tetra CVT Random Figure 1: Overview of adopted meshes for convergence assessment numerical tests.
• Discrete error norms for the stress field: where F h is the set of faces for T h and κ = 1 2 tr(D) (the material is here homogeneous). We remark that the quantity above scales like the internal elastic energy, with respect to the size of the domain and of the elastic coefficients, i.e. ∼ h.
Example 2 (nearly incompressible material). We consider a problem with known analytical solution. A nearly incompressible material is chosen by selecting Lamé constants as λ = 10 5 , µ = 0.5. The test is designed by choosing a required solution for the displacement field and deriving the load f accordingly. The displacement solution is as follows:          u 1 = sin(2πx) 2 cos(2πy) sin(2πy)sin(2πz) 2 − cos(2πz) sin(2πz)sin(2πy) 2 u 2 = sin(2πy) 2 cos(2πz) sin(2πz)sin(2πx) 2 − cos(2πx) sin(2πx)sin(2πz) 2 u 3 = sin(2πz) 2 cos(2πx) sin(2πx)sin(2πy) 2 − cos(2πy) sin(2πy)sin(2πx) 2 In Figure 3 we report the convergence results for the proposed VEM approach. It can be clearly seen that our method shows the expected asymptotic rate of convergence for each kind of mesh. Moreover, also in this case the convergence lines are close to each other and this fact further confirms the robustness of the proposed scheme with respect to element shape.

Example 3 (unloaded body).
We consider a problem with polynomial solution, non-homogeneous Dirichlet boundary conditions and zero loading. We take a homogeneous and isotropic material with Lamé constants λ = 1 and µ = 1 (compressible case). As in the previous examples, the test is defined by choosing a required solution and deriving the correspondig body load f , as indicated in the following: We remark that this is a typical example where the displacement field is nontrivial, while stresses are divergence-free. As expected, this latter feature is numerically satisfied by our VEM scheme. Indeed, in Table 1 the errors E σ σ σ,div are close to the machine precision at each mesh refinement step. The other error behaviours are similar to the ones showed in the previous examples so we do not report such graphichs. Step