A Stress/Displacement Virtual Element Method for Plane Elasticity Problems
AA Stress/Displacement Virtual Element Method for PlaneElasticity Problems
E. Artioli ∗ , S. de Miranda † , C. Lovadina ‡ , L. Patruno § Abstract
The numerical approximation of 2D elasticity problems is considered, in theframework of the small strain theory and in connection with the mixed Hellinger-Reissner variational formulation. A low-order Virtual Element Method (VEM)with a-priori symmetric stresses is proposed. Several numerical tests are provided,along with a rigorous stability and convergence analysis.
The Virtual Element Method (VEM) is a new technology for the approximation ofpartial differential equation problems. VEM was born in 2012, see [7], as an evolution ofmodern mimetic schemes (see for instance [18, 5, 15]), which shares the same variationalbackground of the Finite Element Method (FEM). The initial motivation of VEM isthe need to construct an accurate conforming
Galerkin scheme with the capability todeal with highly general polygonal/polyhedral meshes, including “hanging vertexes”and non-convex shapes. The virtual element method reaches this goal by abandoningthe local polynomial approximation concept , and uses, instead, approximating functionswhich are solution to suitable local partial differential equations (of course, connectedwith the original problem to solve). Therefore, in general, the discrete functions arenot known pointwise, but a limited information of them is at disposal. The key pointis that the available information are indeed sufficient to implement the stiffness matrixand the right-hand side. We remark that VEM is not the only available technology fordealing with polytopal meshes: a brief representative sample of the increasing list oftechnologies that make use of polygonal/polyhedral meshes can be found in [17, 8, 15,9, 11, 13, 28, 30, 32, 35, 39, 37, 40, 41, 42, 23, 34, 43, 20]. We here recall, in particular,the polygonal finite elements and the mimetic discretisation schemes. However, VEMis experiencing a growing interest towards Structural Mechanics problems, also in theengineering community. We here cite the recent works [24, 22, 2, 3, 44, 21, 1] and[4, 19], for instance.In the present paper we apply the VEM concept to two-dimensional elasticity prob-lems in the framework of small displacements and small deformations. More precisely, ∗ Department of Civil Engineering and Computer Science, University of Rome Tor Vergata, Via delPolitecnico 1, 00133 Rome, Italy, [email protected] † DICAM, University of Bologna, Viale Risorgimento 2, 40136 Bologna, Italy, [email protected] ‡ Dipartimento di Matematica, Università di Milano, Via Saldini 50, 20133 Milano, and IMATI delCNR, Via Ferrata 1, 27100 Pavia, Italy, [email protected] § DICAM, University of Bologna, Viale Risorgimento 2, 40136 Bologna, Italy, [email protected] a r X i v : . [ m a t h . NA ] F e b e consider the (mixed) Hellinger-Reissner functional (see, for instance, [12, 14]) as thestarting point of the discretization procedure. Thus, the numerical scheme approxi-mates both the stress and the displacement fields.It is well-known that in the Finite Element practice, designing a stable and accurateelement for the Hellinger-Reissner functional, is not at all a trivial task. Essentially,one is led either to consider quite cumbersome schemes, or to relax the symmetry of theCauchy stress field, or to employ composite elements (a discussion about this issue canbe found in [12], for instance). We here exploit the flexibility of the VEM approach topropose and study a low-order scheme, with a-priori symmetric Cauchy stresses, thatcan be used for general polygons, from triangular shapes on . Furthermore, the methodis robust with respect to the compressibility parameter, and therefore can be used fornearly incompressible situations. Our scheme approximates the stress field by usingtraction degrees of freedom (three per each edge), while the displacement field insideeach polygon is essentially a rigid body motion. The VEM concept is then appliedessentially for the stress field. We also remark that the construction of the discretestress field is somehow similar to the construction of the discrete velocity field used forthe Stokes problem in [10]. Instead, the displacement field is modelled with polynomialfunctions, in accordance with the classical Finite Element procedure.An outline of the paper is as follows. In Section 2 we briefly introduce the Hellinger-Reissner variational formulation of the elasticity problem. Section 3 concerns with thediscrete problem: all the bilinear and linear forms are introduced and detailed. Numeri-cal experiments are reported in Section 4, where suitable error measures are considered.These numerical tests are supported by the stability and convergence analysis devel-oped in Section 5. Finally, Section 6 draws some conclusions, including possible futureextensions of the present study.Throughout the paper, given two quantities a and b , we use the notation a (cid:46) b to mean: there exists a constant C , independent of the mesh-size, such that a ≤ C b .Moreover, we use standard notations for Sobolev spaces, norms and semi-norms (cf.[27], for example).
In this section we briefly present the elasticity problem as it stems from the Hellinger-Reissner principle. More details can be found in [12, 14].
Find ( σσσ, u ) such that − div σσσ = f in Ω σσσ = C εεε ( u ) in Ω u | ∂ Ω = (1)Defining ( · , · ) as the scalar product in L , and a ( σσσ, τττ ) := ( D σσσ, τττ ), a mixed variationalformulation of the problem reads: Find ( σσσ, u ) ∈ Σ × U such that a ( σσσ, τττ ) + ( div τττ , u ) = 0 ∀ τττ ∈ Σ( div σσσ, v ) = − ( f , v ) ∀ v ∈ U (2)where Ω ⊂ R is a polygonal domain, Σ = H ( div ; Ω), U = × L (Ω) , and the loading f ∈ L (Ω) . We recall that div is the vector-valued divergence operator, acting on a2econd order tensor field. Thus, div τττ is, in Cartesian components: ∂τ ij ∂x j (Einstein’ssummation convention is here adopted). The elasticity fourth-order symmetric tensor D := C − is assumed to be uniformly bounded and positive-definite. It is well knownthat problem (2) is well-posed (see [12], for instance). in particular, it holds: || σσσ || Σ + || u || U ≤ C || f || , (3)where C is a constant depending on Ω and on the material tensor C .Note also that the bilinear form a ( · , · ) in (2) can obviously be split as a ( σσσ, τττ ) = X E ∈T h a E ( σσσ, τττ ) with a E ( σσσ, τττ ) := Z E D σσσ : τττ (4)for all σσσ, τττ ∈ Σ. Above, T h is a polygonal mesh of meshsize h .Similarly, it holds( div τττ , v ) = X E ∈T h ( div τττ , v ) E with ( div τττ , v ) E := Z E div τττ · v , (5)for all ( τττ , v ) ∈ Σ × U . Remark 1.
As discussed in [12], estimate (3) does not break down for nearly incom-pressible materials. More precisely, considering the constitutive law: C εεε = 2 µεεε + λ tr( εεε ) Id ∀ symmetric tensor εεε, (6) with λ, µ > the Lame’s parameters and tr( · ) the trace operator, the constant C in (3) can be chosen independent of λ . The key point is that it is sufficient to check the Σ -coercivity of the bilinear form a ( · , · ) in (2) for the subspace: K = { τττ ∈ Σ : ( div τττ , v ) = 0 ∀ v ∈ U } . (7) In fact, there exists a positive constant α such that (see [12]): a ( τττ , τττ ) ≥ α || τττ || ∀ τττ ∈ K, (8) with α independent of λ . We outline the Virtual Element discretization of problem (2). Let {T h } h be a sequenceof decompositions of Ω into general polygonal elements E with h E := diameter( E ) , h := sup E ∈T h h E . In what follows, | E | and | e | = h e will denote the area of E and the length of the side e ∈ ∂E , respectively.We suppose that for all h , each element E in T h fulfils the following assumptions: • ( A1 ) E is star-shaped with respect to a ball of radius ≥ γ h E , • ( A2 ) the distance between any two vertexes of E is ≥ c h E ,where γ and c are positive constants. We remark that the hypotheses above, thoughnot too restrictive in many practical cases, can be further relaxed, as noted in [7].In addition, we suppose that the tensor D is piecewise constant with respect to theunderlying mesh T h . 3 .1 The local spaces Given a polygon E ∈ T h with n E edges, we first introduce the space of local infinitesimalrigid body motions: RM ( E ) = n r ( x ) = a + b ( x − x C ) ⊥ a ∈ R , b ∈ R o . (9)Here above, given c = ( c , c ) T ∈ R , c ⊥ is the counterclock-wise rotated vector c ⊥ =( c , − c ) T , and x C is the baricenter of E . For each edge e of ∂E , we introduce thespace R ( e ) = n t ( s ) = c + d s n c ∈ R , d ∈ R , s ∈ [ − / , / o . (10)Here above, s is a local linear coordinate on e , such that s = 0 corresponds to the edgemidpoint. Furthermore, n is the outward normal to the edge e . Hence, R ( e ) consists ofvectorial functions which have the edge tangential component constant, and the edgenormal component linear along the edge. Our local approximation space for the stressfield is then defined byΣ h ( E ) = n τττ h ∈ H ( div ; E ) : ∃ w ∗ ∈ H ( E ) such that τττ h = C εεε ( w ∗ );( τττ h n ) | e ∈ R ( e ) ∀ e ∈ ∂E ; div τττ h ∈ RM ( E ) o . (11) Remark 2.
Alternatively, the space (11) can be defined as follows. Σ h ( E ) = n τττ h ∈ H ( div ; E ) : τττ h = τττ Th ; curl curl ( D τττ h ) = 0;( τττ h n ) | e ∈ R ( e ) ∀ e ∈ ∂E ; div τττ h ∈ RM ( E ) o . (12) Here above, the equation curl curl ( D τττ h ) = 0 is to be intended in the distribution sense. We remark that, once ( τττ h n ) | e = c e + d e s n is given for all e ∈ ∂E , cf. (10), thequantity div τττ h ∈ RM ( E ) is determined. Indeed, denoting with ϕϕϕ : ∂E → R thefunction such that ϕϕϕ | e := c e + d e s n , the obvious compatibility condition Z E div τττ h · r = Z ∂E ϕϕϕ · r ∀ r ∈ RM ( E ) , (13)allows to compute div τττ h using the c e ’s and the d e ’s. More precisely, setting (cf (9)) div τττ h = ααα E + β E ( x − x C ) ⊥ , (14)from (13) we infer ααα E = 1 | E | Z ∂E ϕϕϕ = 1 | E | X e ∈ ∂E Z e c e β E = 1 R E | x − x C | Z ∂E ϕϕϕ · ( x − x C ) ⊥ = 1 R E | x − x C | X e ∈ ∂E Z e ( c e + d e s n ) · ( x − x C ) ⊥ . (15)The local approximation space for the displacement field is simply defined by, see(9): U h ( E ) = n v h ∈ L ( E ) : v h ∈ RM ( E ) o . (16)We notice that dim(Σ h ( E )) = 3 n E , while dim( U h ( E )) = 3.4 .2 The local bilinear forms Given E ∈ T h , we first notice that, for every τττ h ∈ Σ h ( E ) and v h ∈ U h ( E ), the term Z E div τττ h · v h (17)is computable from the knowledge of the degrees of freedom. Therefore, there is no needto introduce any approximation in the structure of the terms ( div τττ , u ) and ( div σσσ, v )in problem (2). Instead, the term a E ( σσσ h , τττ h ) = Z E D σσσ h : τττ h (18)is not computable for a general couple ( σσσ h , τττ h ) ∈ Σ h ( E ) × Σ h ( E ). As usual in the VEMapproach (see [7], for instance), we then need to introduce a suitable approximation a hE ( · , · ) of a E ( · , · ). To this end, we first define the projection operator Π E : Σ h ( E ) → P ( E ) × s τττ h Π E τττ h a E (Π E τττ h , πππ ) = a E ( τττ h , πππ ) ∀ πππ ∈ P ( E ) × s (19)Above and in the sequel, given a domain ω and an integer k ≥
0, the space P k ( ω )denotes the polynomials up to degree k , defined on ω . Furthermore, given a functionalspace X , X × s denotes the 2 × X .Therefore, the operator in (19) is a projection onto the piecewise constant symmetrictensors.We then set a hE ( σσσ h , τττ h ) = a E (Π E σσσ h , Π E τττ h ) + s E (( Id − Π E ) σσσ h , ( Id − Π E ) τττ h )= Z E D (Π E σσσ h ) : (Π E τττ h ) + s E (( Id − Π E ) σσσ h , ( Id − Π E ) τττ h ) , (20)where s E ( · , · ) is a suitable stabilization term. We propose the following choice: s E ( σσσ h , τττ h ) := κ E h E Z ∂E σσσ h n · τττ h n , (21)where κ E is a positive constant to be chosen. For instance, in the numerical examplesof Section 4, κ E is set equal to tr( D | E ); however, any norm of D | E can be used. Apossible variant of (21) is provided by s E ( σσσ h , τττ h ) := κ E X e ∈ ∂E h e Z e σσσ h n · τττ h n . (22) We need to consider the term, see (2):( f , v h ) = Z Ω f · v h = X E ∈T h Z E f · v h . (23)We remark that, since v h ∈ RM ( E ), computing (23) is possible once a suitablequadrature rule is available for polygonal domains. For such an issue, see for instance[31, 36, 30]. 5 .4 The discrete scheme We are now ready to introduce the discrete scheme. We introduce a global approxima-tion space for the stress field, by glueing the local approximation spaces, see (11):Σ h = n τττ h ∈ H ( div ; Ω) : τττ h | E ∈ Σ h ( E ) ∀ E ∈ T h o . (24)For the global approximation of the displacement field, we take, see (16): U h = n v h ∈ L (Ω) : v h | E ∈ U h ( E ) ∀ E ∈ T h o . (25)Furthermore, given a local approximation of a E ( · , · ), see (20), we set a h ( σσσ h , τττ h ) := X E ∈T h a hE ( σσσ h , τττ h ) . (26)The method we consider is then defined by Find ( σσσ h , u h ) ∈ Σ h × U h such that a h ( σσσ h , τττ h ) + ( div τττ h , u h ) = 0 ∀ τττ h ∈ Σ h ( div σσσ h , v h ) = − ( f , v h ) ∀ v h ∈ U h . (27)Introducing the bilinear form A h : (Σ h × U H ) × (Σ h × U h ) → R defined by A h ( σσσ h , u h ; τττ h , v h ) := a h ( σσσ h , τττ h ) + ( div τττ h , u h ) + ( div σσσ h , v h ) , (28)problem (27) can be written as ( Find ( σσσ h , u h ) ∈ Σ h × U h such that A h ( σσσ h , u h ; τττ h , v h ) = − ( f , v h ) ∀ ( τττ h , v h ) ∈ Σ h × U h . (29)We will prove in Section 5 that our method is first order convergent with respect tothe natural norms, see in particular Theorem 5.8. More precisely, the following errorestimate holds true. || σσσ − σσσ h || Σ + || u − u h || U (cid:46) C h, (30)where C = C (Ω , σσσ, u ) is independent of h but depends on the domain Ω and on theSobolev regularity of σσσ and u . The present section is devoted to the validation of the proposed methodology throughthe assessment of accuracy on a selected number of test problems. Applicability tostructural analysis is then demonstrated through a classical benchmark.
We consider two boundary value problems on the unit square domain Ω = [0 , , withknown analytical solution, discussed in [22, 2]. The material obeys to a homogeneousisotropic constitutive law, see (6), with material parameters assigned in terms of theLamé constants, here set as λ = 1 and µ = 1. Plane strain regime is invoked throughout.The tests are defined by choosing a required solution and deriving the correspondingbody load f , as synthetically indicated in the following:6 Test a u = x − xy u = y − x y f = (31) • Test b ( u = u = sin( πx ) sin( πy ) f = f = − π [ − (3 µ + λ ) sin( πx ) sin( πy ) + ( µ + λ ) cos( πx ) cos( πy )] (32)As it can be observed, Test a is a problem with Dirichlet non-homogeneous boundaryconditions, zero loading and a polynomial solution; whereas Test b has homogeneousDirichlet boundary conditions, trigonometric distributed loads with a trigonometricsolution.In order to test the robustness of the proposed procedure with respect to elementtopology and mesh distortion, six different meshes are considered, as can be inspectedin Fig. 1. Three are structured meshes composed of triangles, quadrilaterals and aset of quads, pentagons and hexagons. In the following, such meshes are denoted bythe letter "S". Three unstructured meshes are considered as well, comprising triangles,quadrilaterals and random polygons; these are denoted by the letter "U". In the nu-merical campaign the mesh size parameter is chosen to be the average edge length,denoted with ¯ h e . We remark that, under mesh assumptions ( A1 ) and ( A2 ) and for aquasi-uniform family of mesh, ¯ h e is indeed equivalent to both h E and h . The accuracyand the convergence rate assessment is carried out using the following error norms: • Discrete error norms for the stress field:
Eσσσ := X e ∈E h κ e | e | Z e | ( σσσ − σσσ h ) n | / , (33)where κ e = κ = tr( D ) (the material is here homogeneous). We remark that thequantity above scales like the internal elastic energy, with respect to the size ofthe domain and of the elastic coefficients.We make also use of the L error on the divergence: Eσσσ , div := X E ∈T h Z E | div ( σσσ − σσσ h ) | / . (34) • L error norm for the displacement field: E u := X E ∈T h Z E | u − u h | / = || u − u h || . (35)Figure 2 reports the ¯ h e − convergence of the proposed method for Test a . As ex-pected, the asymptotic convergence rate is approximately equal to 1 for all the con-sidered error norms and meshes. It is noted that, in this case, the Eσσσ , div plots arenot reported because such a quantity is captured up to machine precision for all theconsidered computational grids. 7 Tri (S)
Quad (S)
Hex (S)
Tri (U)
Quad (U)
Poly (U)
Figure 1: Overview of adopted meshes for convergence assessment numerical tests. log( h e ) -2 -1.5 -1 -0.5 l og ( E < % ) -2-10 1 Tri (S)Quad (S)Hex (S) (a) log( h e ) -2 -1.5 -1 l og ( E < % ) -2-1.5-1-0.5 1 Tri (U)Quad (U)Poly (U) (b) log( h e ) -2 -1.5 -1 -0.5 l og ( E u % ) -2-10 1 Tri (S)Quad (S)Hex (S) (c) log( h e ) -2 -1.5 -1 l og ( E u % ) -1.6-1.4-1.2-1-0.8 1 Tri (U)Quad (U)Poly (U) (d)
Figure 2: ¯ h e − convergence results for Test a on structured and unstructured meshes:(a) and (b) Eσσσ error norm plots, (c) and (d) E u error norm plots.Figure 3 reports ¯ h e − convergence for Test b . Asymptotic converge rate is approxi-mately equal to 1 for all investigated mesh types and error measures, including Eσσσ , div .These results highlight the expected optimal performance of the proposed VEM ap-proach and its robustness with respect to the adopted computational grid.8 og( h e ) -2 -1.5 -1 -0.5 l og ( E < % ) -2-10 1 Tri (S)Quad (S)Hex (S) (a) log( h e ) -2 -1.5 -1 l og ( E < % ) -1.5-1-0.5 1 Tri (U)Quad (U)Poly (U) (b) log( h e ) -2 -1.5 -1 -0.5 l og ( E < ; d i v % ) -2-10 1 Tri (S)Quad (S)Hex (S) (c) log( h e ) -2 -1.5 -1 l og ( E < ; d i v % ) -1.5-1-0.5 1 Tri (U)Quad (U)Poly (U) (d) log( h e ) -2 -1.5 -1 -0.5 l og ( E u % ) -2-10 1 Tri (S)Quad (S)Hex (S) (e) log( h e ) -2 -1.5 -1 l og ( E u % ) -1.5-1-0.5 1 Tri (U)Quad (U)Poly (U) (f)
Figure 3: ¯ h e − convergence results for Test b on structured and unstructured meshes:(a) and (b) Eσσσ error norm plots, (c) and (d)
Eσσσ , div error norm plots, (e) and (f) E u error norm plots. A problem on the unit square domain Ω = [0 , , with known analytical solution, isconsidered. A nearly incompressible material is chosen by selecting Lamé constantsas λ = 10 , µ = 0 .
5. The test is designed by choosing a required solution for thedisplacement field and deriving the load f accordingly. The displacement solution is asfollows: ( u = 0 . πx )) sin(2 πy ) cos(2 πy ) u = − . πy )) sin(2 πx ) cos(2 πx ) . (36)Figure 4 reports the results obtained for both structured and unstructured meshes.In can be clearly seen that the proposed method shows the expected asymptotic rateof convergence also in this case. The present section deals with the classical Cook’s membrane 2D problem [45]. Thegeometry of the domain Ω is presented in Fig. 5 with length data H = 44, H = 16,9 og( h e ) -2 -1.5 -1 -0.5 l og ( E < % ) -1.5-1-0.50 1 Tri (S)Quad (S)Hex (S) (a) log( h e ) -2 -1.5 -1 l og ( E < % ) -1-0.50 1 Tri (U)Quad (U)Poly (U) (b) log( h e ) -2 -1.5 -1 -0.5 l og ( E < ; d i v % ) -1.5-1-0.50 1 Tri (S)Quad (S)Hex (S) (c) log( h e ) -2 -1.5 -1 l og ( E < ; d i v % ) -1.2-1-0.8-0.6-0.4 1 Tri (U)Quad (U)Poly (U) (d) log( h e ) -2 -1.5 -1 -0.5 l og ( E u % ) -1.5-1-0.500.5 1 Tri (S)Quad (S)Hex (S) (e) log( h e ) -2 -1.5 -1 l og ( E u % ) -1-0.50 1 Tri (U)Quad (U)Poly (U) (f)
Figure 4: Results for the nearly incompressible test on structured and unstructuredmeshes: (a) and (b) convergence of
Eσσσ , (c) and (d) convergence of
Eσσσ , div , (e) and (f)convergence for E u . L = 48. The loading is given by a constant tangential traction q = 6 .
25 on theright edge of the domain. The Young modulus, E , is set equal to 70 and two Poissonratios are considered, one corresponding to ν = 1 / ν = 0 . v A , thevertical displacement of point A (see Fig. 5), approximated as the vertical displacementat the centroid of the closest polygon. In particular, Fig. 7(a) corresponds to the case inwhich ν = 1 / Quad
CVor
RVor
Figure 6: Cook’s membrane. Examples of the adopted meshes.ported in Fig. 7. We remark that, inside the polygons, the stress distribution σσσ h is notknown, but its projection Π E σσσ h onto the constant tensors is (cf. (19)). Thus, we haveused this latter quantity to compute the von Mises equivalent stress displayed in Fig.7. Finally, the results refer to the case ν = 1 /
3, being the nearly incompressible caseextremely similar.
In this section, we provide a rigorous analysis of the proposed VEM method. For all E ∈ T h , we first introduce the space: e Σ( E ) := n τττ ∈ H ( div ; E ) : ∃ w ∈ H ( E ) such that τττ = C εεε ( w ) o . (37)The global space e Σ is defined as 11 = h e v A [ mm ] QuadCVorRVorReference (a) = h e v A [ mm ] QuadCVorRVorReference (b)
Figure 7: Convergence of the tip vertical displacement v A : (a) ν = 1 / ν =0 . (a) (b) (c) Figure 8: Contours representing the von Mises equivalent stress distributions for ν =1 /
3: (a) Quad, (b) CVor, (c) RVor. e Σ := n τττ ∈ H ( div ; Ω) : ∃ w ∈ H (Ω) such that τττ = C εεε ( w ) o . (38)In the sequel, given a measurable subset ω ⊆ Ω and r >
2, we will use the space W r ( ω ) := n τττ : τττ ∈ L r ( ω ) × s , div τττ ∈ L ( ω ) o , (39)equipped with the obvious norm. Under our assumptions on the mesh, we recall thefollowing version of the Korn’s inequality:inf r ∈ RM ( E ) (cid:16) h − E || v − r || ,E + | v − r | ,E (cid:17) (cid:46) || εεε ( v ) || ,E ∀ v ∈ H ( E ) . (40)Given v ∈ H ( E ) , the above inequality can be derived by classical results (see [33],for instance), and by choosing r v ∈ RM ( E ) such that R E ( v − r v ) = .We will also use the following result. Lemma 5.1.
Suppose that assumptions ( A1 ) and ( A2 ) are fulfilled. Given E ∈ T h ,let w ∈ H ( E ) be a solution of the problem: − div ( C εεε ( w )) = g in E ( C εεε ( w )) n = h on ∂E, (41) where g ∈ L ( E ) and h ∈ L ( ∂E ) satisfy the compatibility condition Z E g · r + Z ∂E h · r = 0 ∀ r ∈ RM ( E ) . (42) Then it holds: || C εεε ( w ) || ,E (cid:46) h E || g || ,E + h / E || h || ,∂E . (43) Proof.
For every r ∈ RM ( E ), we have || C εεε ( w ) || ,E (cid:46) Z E C εεε ( w ) : εεε ( w )= Z E C εεε ( w ) : εεε ( w − r ) = Z E g · ( w − r ) + Z ∂E h · ( w − r ) , (44)by which we get || C εεε ( w ) || ,E (cid:46) || g || ,E || w − r || ,E + || h || ,∂E || w − r || ,∂E . (45)Under assumptions ( A1 ) and ( A2 ), the Agmon’s inequality then gives || C εεε ( w ) || ,E (cid:46) || g || ,E || w − r || ,E + || h || ,∂E (cid:16) h − / E || w − r || L ( E ) + h / E | w − r | H ( E ) (cid:17) . (46)Estimate (43) now follows from (40). We now introduce the local interpolation operator I E : W r ( E ) → Σ h ( E ), defined asfollows. Given τττ ∈ W r ( E ), I E τττ ∈ Σ h ( E ) is determined by: Z ∂E ( I E τττ ) n · ϕϕϕ ∗ = Z ∂E τττ n · ϕϕϕ ∗ ∀ ϕϕϕ ∗ ∈ R ∗ ( ∂E ) , (47)where R ∗ ( ∂E ) = n ϕϕϕ ∗ ∈ L ( ∂E ) : ϕϕϕ ∗| e = γγγ e + δ e ( x − x C ) ⊥ γγγ e ∈ R , δ e ∈ R , ∀ e ∈ ∂E o . (48)If τττ is not sufficiently regular, the integral in the right-hsnd side of (47) is intended asa duality between W − r ,r ( ∂E ) and W r ,r ( ∂E ) . If τττ is a regular function, the abovecondition is equivalent to require: Z e ( I E τττ ) n = Z e τττ n ∀ e ∈ ∂E ; Z e ( I E τττ ) n · ( x − x C ) ⊥ = Z e τττ n · ( x − x C ) ⊥ ∀ e ∈ ∂E. (49)The following result shows, in particular, that I E τττ ∈ Σ h ( E ) is well-defined byconditions (47). 13 emma 5.2. If τττ h ∈ Σ h ( E ) , then Z ∂E τττ h n · ϕϕϕ ∗ = 0 ∀ ϕϕϕ ∗ ∈ R ∗ ( ∂E ) (50) imply τττ h = .Proof. First, recall that for τττ h ∈ Σ h ( E ) it holds ( τττ h n ) | e = c e + d e s n for each edge e ∈ ∂E , cf. (10) and (11). By (50), choosing ϕϕϕ ∗ such that ϕϕϕ ∗| e = γγγ e for each e ∈ ∂E , itfollows that ( τττ h n ) | e = d e s n . Choosing now ϕϕϕ ∗| e = δ e ( x − x C ) ⊥ , conditions (50) thengive d e Z e s n · ( x − x C ) ⊥ = 0 ∀ e ∈ ∂E. (51)A direct computation (for instance by using the Cavalieri-Simpson rule) shows that(51) is equivalent to d e | e | n · ( q e − p e ) ⊥ = 0 ∀ e ∈ ∂E. (52)Above, p e and q e denote the endpoints of e . From (52) we infer d e = 0 for each e ∈ ∂E ,which concludes the proof.The global interpolation operator I h : W r (Ω) → Σ h is then defined by simplyglueing the local contributions provided by I E . More precisely, we set ( I h τ ) | E := I E τττ | E for every E ∈ T h and τττ ∈ W r (Ω). Proposition 5.3.
Under assumptions ( A1 ) and ( A2 ) , for the interpolation operator I E defined in (49) , the following estimates hold: || τττ − I E τττ || ,E (cid:46) h E | τττ | ,E ∀ τττ ∈ e Σ( E ) ∩ H ( E ) × s . (53) || div ( τττ − I E τττ ) || ,E (cid:46) h E | div τττ | ,E ∀ τττ ∈ e Σ( E ) ∩ H ( E ) × s s.t. div τττ ∈ H ( E ) . (54) Proof.
Let τττ ∈ e Σ( E ) ∩ H ( E ) × s , and let w ∈ H ( E ) be such that τττ = C εεε ( w ), see(37). Furthermore, consider I E τττ ∈ Σ h ( E ) and w ∗ ∈ H ( E ) such that I E τττ = C εεε ( w ∗ ),see (11). Hence, setting δδδ := ( w − w ∗ ) ∈ H ( E ) , it holds: τττ − I E τττ = C εεε ( δδδ ) . (55)Furthermore, using (49), (14) and (15), we infer that δδδ ∈ H ( E ) satisfies: div ( C εεε ( δδδ )) = div τττ − | E | X e ∈ ∂E Z e τττ n − ( x − x C ) ⊥ R E | x − x C | X e ∈ ∂E Z e τττ n · ( x − x C ) ⊥ in E ( C εεε ( δδδ )) n = X e ∈ ∂E (cid:18) τττ n − | e | Z e τττ n (cid:19) χ e on ∂E, (56)where χ e denotes the characteristic function of the edge e . Applying Lemma 5.1 with:14 g := 1 | E | X e ∈ ∂E Z e τττ n + ( x − x C ) ⊥ R E | x − x C | X e ∈ ∂E Z e τττ n · ( x − x C ) ⊥ − div τττ h := X e ∈ ∂E (cid:18) τττ n − | e | Z e τττ n (cid:19) χ e , (57)we get || τττ − I E τττ || ,E = || C εεε ( δδδ ) || ,E (cid:46) h E || g || ,E + h / E || h || ,∂E . (58)We now estimate g and h . We denote respectively with Π ,E , Π RM,E and Π ,∂E the L -projection operators onto the constant functions on E , onto the space RM ( E ) (see (9)),and on the piecewise constant functions on ∂E (with respect to the edge subdivisionof ∂E ).The divergence theorem and a direct computation show that:1 | E | X e ∈ ∂E Z e τττ n + ( x − x C ) ⊥ R E | x − x C | X e ∈ ∂E Z e τττ n · ( x − x C ) ⊥ = Π RM,E div τττ . (59)Therefore, from the first equation of (57), we have g = Π RM,E div τττ − div τττ . (60)Noting that P ( E ) ⊂ RM ( E ), from the properties of the L projection operator, wethen get || g || ,E = || Π RM,E div τττ − div τττ || ,E ≤ || Π ,E div τττ − div τττ || ,E (cid:46) || div τττ || ,E (61)and || g || ,E = || Π RM,E div τττ − div τττ || ,E ≤ || Π ,E div τττ − div τττ || ,E (cid:46) h E | div τττ | ,E . (62)For the second equation of (57), we remark that: h = X e ∈ ∂E (cid:18) τττ n − | e | Z e τττ n (cid:19) χ e = X e ∈ ∂E (cid:18) τττ − | e | Z e τττ (cid:19) n χ e = ( τττ − Π ,∂E τττ ) n . (63)Hence, using a standard approximation estimate and a trace inequality, we get || h || ,∂E = || ( τττ − Π ,∂E τττ ) n || ,∂E ≤ || τττ − Π ,∂E τττ || ,∂E (cid:46) h / E | τττ | / ,∂E (cid:46) h / E | τττ | ,E . (64)Taking into account (61) and (64), from (58) we obtain estimate (53).We now notice that from (55), (56) and (57), we have: div ( τττ − I E τττ ) = − g . (65)Then, using (62), we immediately get (54).15 .3 Proving the el lipticity-on-the-kernel condition We first notice that by (19), (20) and (21), using the techniques of [7, 16], one has: || τττ h || ,E (cid:46) a hE ( τττ h , τττ h ) (cid:46) || τττ h || ,E ∀ τττ h ∈ Σ h ( E ) . (66)We also notice that (see (24), (11) and (25), (16)): div (Σ h ) ⊆ U h . (67)As a consequence, introducing the discrete kernel K h ⊆ Σ h : K h = { τττ h ∈ Σ h : ( div τττ h , v h ) = 0 ∀ v h ∈ U h } , (68)we infer that τττ h ∈ K h implies div τττ h = . Hence, it holds: || τττ h || Σ = || τττ h || ∀ τττ h ∈ K h . (69)We are now ready to prove the following ellipticity-on-the-kernel condition. Proposition 5.4.
For the method described in Section 3, there exists a constant α > such that a h ( τττ h , τττ h ) ≥ α || τττ h || ∀ τττ h ∈ K h . (70) Proof.
By recalling (26), from (66) we get the existence of α > a h ( τττ h , τττ h ) ≥ α || τττ h || ∀ τττ h ∈ Σ h . (71)Estimate (70) now follows by recalling (69). Remark 3.
Notice that for our method it holds K h ⊂ K , where K is defined by (7) .Considering an isotropic material, see (6) , from Remark 1 we infer that the coercivityconstant α can be chosen independent of λ . Therefore, our numerical method does notsuffer from volumetric locking (see [26], for instance) and can be used also for nearlyincompressible materials. This feature is confirmed by the numerical tests presented inSection 4. inf-sup condition We start by stating the following proposition, which can be derived by regularity resultsfor the elasticity problem on Lipschitz domains (see [25], for example).
Proposition 5.5.
Given the polygonal domain Ω , there exist s > and β ∗ > suchthat sup τττ ∈ W s (Ω) ( div τττ , v ) || τττ || W s (Ω) ≥ β ∗ || v || , Ω ∀ v ∈ L (Ω) , (72) where W s (Ω) is the Banach space defined by (39) . We are now ready to prove the discrete inf-sup condition for our choice of theapproximation spaces. 16 roposition 5.6.
Suppose that assumptions ( A1 ) and ( A2 ) are fulfilled. There exists β > such that sup τττ h ∈ Σ h ( div τττ h , v h ) || τττ h || Σ ≥ β || v h || , Ω ∀ v h ∈ U h . (73) Proof.
We will apply Fortin’s criterion (see [12]), using the operator I h : W s (Ω) → Σ h ,see (49) for the definition of the local contributions. More precisely, we will show thatit holds: Z Ω div ( I h τττ ) · v h = Z Ω div τττ · v h ∀ v h ∈ U h , ∀ τττ ∈ W s (Ω) , ||I h τττ || Σ (cid:46) || τττ || W s (Ω) ∀ τττ ∈ W s (Ω) . (74)Together with (72), conditions (74) imply (73), see [12].To prove the first condition in (74), recalling that v h | E ∈ RM ( E ), it is sufficient toshow that: Z E div ( I E τττ ) · r = Z E div τττ · r ∀ r ∈ RM ( E ) , ∀ E ∈ T h . (75)The above equation directly follows from the divergence theorem and definition (47).We now prove the continuity estimate (i.e. the second equation in (74)). We willexploit again Lemma 5.1. More precisely, we take w ∗ ∈ H ( E ) such that I E τττ = C εεε ( w ∗ ). It follows that w ∗ solves, cf. (59): div ( C εεε ( w ∗ )) = Π RM,E div τττ ( C εεε ( w ∗ )) n = X e ∈ ∂E c e χ e on ∂E, (76)where the c e ’s are given by the dualities for the couple < W − s ,s ( ∂E ) , W s ,s ( ∂E ) > : c e := 1 | e | ( < τττ n , χ e t > t + < τττ n , χ e n > n ) . (77)From (76) we obviously deduce || div ( I E τττ ) || ,E = || div ( C εεε ( w ∗ )) || ,E = || Π RM,E div τττ || ,E ≤ || div τττ || ,E . (78)We now apply Lemma 5.1 with: g := − Π RM,E div τττ , h := X e ∈ ∂E c e χ e , (79)and estimate || h || ,∂E . We start by noting that: || h || ,∂E = X e ∈ ∂E | c e | | e | / (cid:46) h / E X e ∈ ∂E | c e | / . (80)A duality estimate and a trace bound shows that < τττ n , χ e t > (cid:46) || τττ n || W − s ,s ( ∂E ) || χ e t || W s ,s ( ∂E ) (cid:46) || τττ n || W − s ,s ( ∂E ) (cid:46) || τττ || W s ( E ) . (81)17imilarly, it holds: < τττ n , χ e n > (cid:46) || τττ || W s ( E ) . (82)From (77), (81) and (82) we get | c e | (cid:46) h − E || τττ || W s ( E ) , (83)by which we deduce, see (80): || h || ,∂E (cid:46) h − / E || τττ || W s ( E ) . (84)Lemma 5.1 thus gives ||I E τττ || ,E = || C εεε ( w ∗ ) || ,E (cid:46) || τττ || W s ( E ) . (85)The continuity estimate in (74) now follows by collecting all the local estimates (85). We denote with P ( T h ) the space of piecewise constant functions with respect to thegiven mesh T h . We can prove the Proposition: Proposition 5.7.
Suppose that assumptions ( A1 ) and ( A2 ) are fulfilled. For every ( σσσ I , u I ) ∈ Σ h × U h and every σσσ π ∈ P ( T h ) × s , the following error equation holds: || σσσ − σσσ h || Σ + || u − u h || U (cid:46) || σσσ − σσσ I || Σ + || u − u I || U ++ h || div σσσ I || , Ω + || σσσ − σσσ π || , Ω . (86) Proof.
Given ( σσσ I , u I ) ∈ Σ h × U h , we form ( σσσ h − σσσ I , u h − u I ) ∈ Σ h × U h . Then, usingthe ellipticity-on-the-kernel condition of Proposition 5.4 and the inf-sup condition ofProposition 5.6, there exists ( τττ h , v h ) ∈ Σ h × U h such that (see [12] and [14], for instance): || τττ h || Σ + || v h || U (cid:46) || σσσ h − σσσ I || Σ + || u h − u I || U (cid:46) A h ( σσσ h − σσσ I , u h − u I ; τττ h , v h ) . (88)We have A h ( σσσ h − σσσ I , u h − u I ; τττ h , v h ) = A h ( σσσ h , u h ; τττ h , v h ) − A h ( σσσ I , u I ; τττ h , v h )= − ( f , v h ) − A h ( σσσ I , u I ; τττ h , v h )= A ( σσσ, u ; τττ h , v h ) − A h ( σσσ I , u I ; τττ h , v h )= [ a ( σσσ, τττ h ) − a h ( σσσ I , τττ h )] + ( div τττ h , u − u I ) + ( div ( σσσ − σσσ I ) , v h )= T + T + T (89)Concerning T , it holds: 18 = X E ∈T h h a E ( σσσ, τττ h ) − a hE ( σσσ I , τττ h ) i = X E ∈T h (cid:2) a E ( σσσ, τττ h ) − a E (Π E σσσ I , Π E τττ h ) − κ E h E Z ∂E [( Id − Π E ) σσσ I n ] · [( Id − Π E ) τττ h n ] (cid:3) = X E ∈T h (cid:2) a E ( σσσ − σσσ π , τττ h ) − a E (Π E ( σσσ I − σσσ π ) , Π E τττ h ) − κ E h E Z ∂E [( Id − Π E ) σσσ I n ] · [( Id − Π E ) τττ h n ] (cid:3) . (90)We have, using the continuity of a E ( · , · ) and of Π E : X E ∈T h (cid:2) a E ( σσσ − σσσ π , τττ h ) − a E (Π E ( σσσ I − σσσ π ) , Π E τττ h ) (cid:3) (cid:46) ( || σσσ − σσσ π || , Ω + || σσσ I − σσσ π || , Ω ) || τττ h || , Ω (cid:46) ( || σσσ − σσσ π || , Ω + || σσσ I − σσσ || , Ω + || σσσ − σσσ π || , Ω ) || τττ h || , Ω (cid:46) ( || σσσ − σσσ π || , Ω + || σσσ − σσσ I || , Ω ) || τττ h || Σ . (91)Furthermore, it holds: X E ∈T h κ E h E Z ∂E [( Id − Π E ) σσσ I n ] · [( Id − Π E ) τττ h n ] (cid:3) (cid:46) X E ∈T h h / E || ( Id − Π E ) σσσ I n || ,∂E h / E || ( Id − Π E ) τττ h n || ,∂E (92)Under assumptions ( A1 ) and ( A2 ), we notice that, given τττ h ∈ Σ h ( E ), we have the 1Dinverse estimate on ∂E : h / E || τττ h n || ,∂E (cid:46) || τττ h n || − / ,∂E ∀ τττ h ∈ Σ h ( E ) . (93)Using the techniques developed in [6], we deduce the scaled trace estimate: || τττ h n || − / ,∂E (cid:46) || τττ h || ,E + h E || div τττ h || ,E ∀ τττ h ∈ Σ h ( E ) . (94)Hence, we get: h / E || τττ h n || ,∂E (cid:46) || τττ h || ,E + h E || div τττ h || ,E ∀ τττ h ∈ Σ h ( E ) . (95)From (92) and (95) we then deduce X E ∈T h κ E h E Z ∂E [( Id − Π E ) σσσ I n ] · [( Id − Π E ) τττ h n ] (cid:3) (cid:46) X E ∈T h (cid:16) || ( Id − Π E ) σσσ I || ,E + h E || div σσσ I || ,E (cid:17) / || τττ h || Σ . (96)Since it holds, using also the L continuity of Π E : || ( Id − Π E ) σσσ I || ,E = || ( σσσ I − σσσ π ) + Π E ( σσσ π − σσσ I ) || ,E (cid:46) || σσσ I − σσσ π || ,E (cid:46) || σσσ I − σσσ || ,E + || σσσ − σσσ π || ,E . (97)19herefore, we get: X E ∈T h κ E h E Z ∂E [( Id − Π E ) σσσ I n ] · [( Id − Π E ) τττ h n ] (cid:3) (cid:46) ( || σσσ I − σσσ || , Ω + || σσσ − σσσ π || , Ω + h || div σσσ I || , Ω ) || τττ h || Σ . (98)Combining (90), (91) and (98), we infer T (cid:46) ( || σσσ − σσσ I || , Ω + || σσσ − σσσ π ) || , Ω + h || div σσσ I || , Ω ) || τττ h || Σ . (99)Regarding T , T and T , one obviously have: ( T (cid:46) || u − u I || U || τττ h || Σ T (cid:46) || σσσ − σσσ I || Σ || v h || U . . (100)From (88), (89), (99) and (100), we get: || σσσ h − σσσ I || Σ + || u h − u I || U (cid:46) (cid:16) || σσσ − σσσ I || Σ + || σσσ − σσσ π ) || , Ω + h || div σσσ I || , Ω + || u − u I || U (cid:17) ( || τττ h || Σ + || v h || U ) (101)Estimate (86) follows from the triangle inequality, estimate (101) and bound (87).We are now ready to state and prove our main convergence result. Theorem 5.8.
Let ( σσσ, u ) ∈ Σ × U be the solution of Problem (2) , and let ( σσσ h , bbu h ) ∈ Σ h × U h be the solution of the discrete problem (27) . Suppose that assumptions ( A1 ) and ( A2 ) are fulfilled. Assuming σσσ | E ∈ H ( E ) × s and ( div σσσ ) | E ∈ H ( E ) , the followingestimate holds true: || σσσ − σσσ h || Σ + || u − u h || U (cid:46) C (Ω , σσσ, u ) h, (102) where C (Ω , σσσ, u ) is independent of h but depends on the domain Ω and on the Sobolevregularity of σσσ and u .Proof. In Proposition 5.7 let us choose σσσ I = I h σσσ ∈ Σ h as detailed in Section 5.1, u I = P u ∈ U h and σσσ π = P σσσ ∈ Σ h . Estimate (102) easily follows from Proposition5.3 and standard approximation results. Remark 4.
An alternative way to develop the stability and error analysis might be theuse of suitable mesh-dependent norms, as detailed in [29] for the Poisson problem inmixed form.
We have proposed, numerically tested and analysed a new Virtual Element Method forthe Hellinger-Reissner formulation of two-dimensional elasticity problems. Our schemeis low-order, it has a-priori symmetric stresses and it optimally converges. Possiblefuture developments of the present study include the design of higher-order schemes inthe framework of the same variational principle. In addition, accurate post-processeddisplacements might be considered and used for mesh adaptive strategies, based onsuitable a-posteriori error estimators. 20 knowledgements
EA gratefully acknowledges the partial financial support of the Italian Minister ofUniversity and Research, MIUR (Program: Consolidate the Foundations 2015; Project:BIOART; Grant number (CUP): E82F16000850005).
References [1] O. Andersen, H. M. Nilsen, and X. Raynaud,
On the use of thevirtual element method for geomechanics on reservoir grids , online: https://arxiv.org/abs/1606.09508 .[2] E. Artioli, L. Beirão Da Veiga, C. Lovadina, and E. Sacco,
Arbitrary order 2Dvirtual elements for polygonal meshes: Part I, elastic problem , submitted for pub-lication, and online on: https://arxiv.org/abs/1701.06670 .[3] E. Artioli, L. Beirão Da Veiga, C. Lovadina, and E. Sacco,
Arbitrary order 2Dvirtual elements for polygonal meshes: Part II, inelastic problems , submitted forpublication, and online on: https://arxiv.org/abs/1701.06676 .[4] L. Beirão da Veiga, F. Brezzi, and L. D. Marini,
Virtual Elements for linear elas-ticity problems , Siam. J. Numer. Anal. (2013), 794–812.[5] L. Beirão da Veiga, K. Lipnikov, and G. Manzini, The mimetic finite differencemethod for elliptic problems , Springer, series MS&A (vol. 11), 2014.[6] L. Beirão da Veiga, C. Lovadina, and A. Russo,
Stability analysis for the VirtualElement Methods , Preprint arXiv:1607.05988. Submitted for publication.[7] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo,
Basic principles of virtual element methods , Math. Models Methods Appl. Sci. (2013), no. 1, 199–214.[8] L. Beirão da Veiga, K. Lipnikov, and G. Manzini, Arbitrary-order nodal mimeticdiscretizations of elliptic problems on polygonal meshes , SIAM J. Numer. Anal. (2011), no. 5, 1737–1760.[9] L. Beirão da Veiga, K. Lipnikov, and G. Manzini, The mimetic finite differ-ence method for elliptic problems , MS&A. Modeling, Simulation and Applications,vol. 11, Springer, 2014. MR 3135418[10] L. Beirão Da Veiga, C. Lovadina, and G. Vacca,
Divergence free virtual elementsfor the Stokes problem on polygonal meshes , to appear on ESAIM: M2AN.[11] J. E. Bishop,
A displacement-based finite element formulation for general polyhedrausing harmonic shape functions , Internat. J. Numer. Methods Engrg. (2014),no. 1, 1–31. MR 3146670[12] D. Boffi, F. Brezzi, and M. Fortin, Mixed finite element methods and applications ,Springer Series in Computational Mathematics, vol. 44, Springer, Heidelberg, 2013.MR 3097958 2113] J. Bonelle and A. Ern,
Analysis of compatible discrete operator schemes for ellipticproblems on polyhedral meshes , ESAIM Math. Model. Numer. Anal. (2014),no. 2, 553–581. MR 3177857[14] D. Braess, Finite elements. Theory, fast solvers, and applications in elasticitytheory. , third ed., Cambridge University Press, 2007.[15] F. Brezzi, A. Buffa, and K. Lipnikov,
Mimetic finite differences for elliptic prob-lems , Math. Mod. Numer. Anal. (2009), 277–295.[16] F. Brezzi, R. S. Falk, and L. D. Marini, Basic principles of mixed virtual elementmethods , ESAIM Math. Model. Numer. Anal.[17] F. Brezzi, K. Lipnikov, and M. Shashkov,
Convergence of the mimetic finite differ-ence method for diffusion problems on polyhedral meshes , SIAM J. Numer. Anal. (2005), no. 5, 1872–1896.[18] F. Brezzi, K. Lipnikov, M Shashkov, and V. Simoncini, A new discretizationmethodology for diffusion problems on generalized polyhedral meshes , Comput.Methods Appl. Mech. Engrg. (2007), 3682–3692.[19] F. Brezzi and L.D. Marini,
Virtual Element Method for plate bending problems ,Comput. Methods Appl. Mech. Engrg. (2012), 455–462.[20] A. Cangiani, E.H. Georgoulis, and P. Houston, hp-version discontinuous Galerkinmethods on polygonal and polyhedral meshes , Math. Mod. Meth. Appl. Sci. (2014), no. 10, 2009–2041.[21] H. Chi, L. Beirão da Veiga, and G. H. Paulino, Some basic formulations of thevirtual element method (vem) for finite deformations , Comput. Meth. Appl. Mech.Engrg., in press.[22] L. Beirão da Veiga, C. Lovadina, and D. Mora,
A virtual element method forelastic and inelastic problems on polytope meshes , Computer Methods in AppliedMechanics and Engineering (2015), 327 – 346.[23] D. Di Pietro and A. Alexandre Ern,
A hybrid high-order locking-free method forlinear elasticity on general meshes , Comput. Methods Appl. Mech. Engrg. (2015), no. 0, 1–21.[24] A. L. Gain, C. Talischi, and G. H. Paulino,
On the virtual element method for three-dimensional linear elasticity problems on arbitrary polyhedral meshes , Comput.Methods Appl. Mech. Engrg. (2014), 132–160. MR 3269894[25] R. Herzog, C. Meyer, and G. Wachsmuth,
Integrability of displacement and stressesin linear and nonlinear elasticity with mixed boundary conditions , J. Math. Anal.Appl. (2011), 802–813.[26] T. J. R. Hughes,
The finite element method. linear static and dynamic finite ele-ment analysis. , second ed., Dover, 2000.[27] J.-L. Lions and E. Magenes,
Problèmes aux limites non homogènes et applications.Vol. 1 , Travaux et Recherches Mathématiques, No. 17, Dunod, Paris, 1968.2228] K. Lipnikov, G. Manzini, and M. Shashkov,
Mimetic finite difference method , J.Comput. Phys. (2014), 1163–1227.[29] C. Lovadina and R. Stenberg,
Energy norm a posteriori error estimates for mixedfinite element methods , Math. Comp. (2006), no. 256, 1659–1674 (electronic).MR 2240629 (2007h:65129)[30] S. E. Mousavi and N. Sukumar, Numerical integration of polynomials and discon-tinuous functions on irregular convex polygons and polyhedrons , Comput. Mech. (2011), no. 5, 535–554.[31] S. E. Mousavi, H. Xiao, and N. Sukumar, Generalized Gaussian quadrature ruleson arbitrary polygons , International Journal for Numerical Methods in Engineering (2010), no. 1, 99–113.[32] S. Natarajan, S. Bordas, and D. R. Mahapatra, Numerical integration over arbi-trary polygonal domains based on Schwarz–Christoffel conformal mapping , Int. J.Numer. Meth. Eng. (2009), no. 1, 103–134.[33] O. Oleinik and V. Kondratiev, On Korn’s inequalities , C.R. Acad. Sci. Paris (1989), 483–487.[34] A. Rand, A. Gillette, and C. Bajaj,
Interpolation error estimates for mean valuecoordinates over convex polygons , Advances in Computational Mathematics (2013), no. 2, 327–347.[35] S. Rjasanow and S. Weißer, Higher order BEM-based FEM on polygonal meshes ,SIAM J. Numer. Anal. (2012), no. 5, 2357–2378.[36] A. Sommariva and M. Vianello, Product Gauss cubature over polygons based ongreen’s integration formula , BIT Numerical Mathematics (2007), no. 2, 441–453.[37] N. Sukumar and A. Tabarraei, Conforming polygonal finite elements , Internat. J.Numer. Methods Engrg. (2004), no. 12, 2045–2066.[38] Dassault Systèmes, Abaqus documentation , Providence, RI, 2011.[39] C. Talischi and G. H. Paulino,
Addressing integration error for polygonal finiteelements through polynomial projections: a patch test connection , Math. ModelsMethods Appl. Sci. (2014), no. 8, 1701–1727.[40] C. Talischi, G. H. Paulino, A. Pereira, and I. F. M. Menezes, Polygonal finite ele-ments for topology optimization: A unifying paradigm , Internat. J. Numer. Meth-ods Engrg. (2010), no. 6, 671–698.[41] M. Vohralik and B. I. Wohlmuth, Mixed finite element methods: implementationwith one unknown per element, local flux expressions, positivity, polygonal meshes,and relations to other methods , Math. Models Methods Appl. Sci. (2013), no. 5,803–838.[42] E. Wachspress, Rational bases for convex polyhedra , Comput. Math. Appl. (2010), no. 6, 1953–1956.[43] J. Wang and X. Ye, A weak Galerkin finite element method for second-order ellipticproblems , J. Comput. Appl. Math. (2013), 103–115.2344] P. Wriggers, W.T. Rust, and B.D. Reddy,