Virtual Element Method: an equilibriumbased stress recovery procedure
Abstract
Within the framework of the displacementbased Virtual Element Method (VEM) for plane elasticity a significant problem is represented by an accurate evaluation of the stress field. In particular, in the classical VEM formulation, a suitable operator which maps to the strain field is introduced in order to allow the calculation of the stiffness matrix. The stress field is then computed using that strian field, by using the constitutive law. Considering for example a firstorder formulation for a homogeneous material, strains are locally mapped onto constant functions, and stresses are accordingly piecewise constant. However, the virtual displacements might engender more complex strain fields for polygons which are not triangles. In this paper, Recovery by Compatibility in Patches is used in order to mitigate such an effect and, thus, enhance the accuracy of the recovered stress field. The procedure is simple, efficient and can be readily implemented in existing codes. Numerical tests confirm the soundness of the proposed approach.
keywords:
Virtual Element Method, Stress recovery, RCP[cor33]Corresponding author
1 Introduction
The Virtual Element Method (VEM) is a relatively new and very powerful discretisation scheme which is rapidly attracting the interest of the scientific community. The technique is wellknown for its formal elegance and flexibility, which allows to adopt meshes composed of general polygons/polyhedra as well as allowing the presence of hanging nodes and nonconforming grids (1); (2); (3); (4). The method, originally proposed in 2012 in its displacementbased formulation and presented for the Laplace operator in a twodimensional context (1), is rapidly developing and has been already applied to numerous physical problems (3); (5); (6); (7); (8); (9); (10); (11); (12); (13); (14); (15) as well as extended to threedimensional cases (16) and mixed formulations (17); (18).
Differently from the finite element framework, in the virtual element technique the local polynomial approximation concept is abandoned and the functions used for the discretization procedure are not known pointwise, in general. However, a careful selection of the degrees of freedom and assuming that the unknown fields satisfy appropriate differential equations, make it possible is to compute the stiffness matrix and the load vector.
In particular, for plane elasticity problems, classical VEM formulations assume that displacements are completely known only at the element boundaries while the internal field is unknown and, for highorder schemes, characterised only by means of internal degrees of freedom representing appropriate integral quantities.
Due to the virtual nature of the displacement field, strains (the symmetric gradient of the displacements) are not directly computable and need to be further approximated, in order to compute the stiffness matrix. The standard procedure consists in defining a strain field, polynomial inside the element, which can be determined starting from the displacement degrees od freedom. For example, in the case of the firstorder formulation, the aforementioned strategy leads to constant strains inside the elements, irrespective of the number of element edges, and thus irrespective of the number of displacement degrees of freedom. Therefore, applying the constitutive law, the stress field is essentially piecewise constant as well. Especially for meshes consisting of polygons with numerous edges, this often leads to not completely satisfactory results.
In this paper, a modified version of the Recovery by Compatibility in Patches (RCP) (19); (20) is used to devise an alternative procedure, which allows to partly alleviate the above drawback and compute accurate stresses in the displacementbased VEM framework. Patchbased stress recovery techniques, initially inspired by the Superconvergent Patch Recovery (SPR) procedure proposed by Zienkiewicz and Zhu (21); (22), have quickly developed and found numerous applications, ranging from mesh adaptivity (23); (24) to the recovery of accurate stresses for crack propagation in XFEM analyses (25). In contrast with classical approaches, the common characteristic of such techniques is to recover enhanced stress fields by considering groups of surrounding elements (i.e. patches). In the original SPR procedure this is obtained by calculating stresses as a leastsquare interpolation of values evaluated at superconvergent points. Later, the procedure has been refined (26); (27); (28); (29); (30); (31) introducing, for example, equilibrium constraints in order to lead to more accurate results. Alternative procedures, based on variational principles such as the Recovery by Equilibrium in Patches (REP) (32); (33); (34), have been proved to obtain extremely interesting results.
RCP has been originally proposed as a patchbased recovery technique which, by enforcing equilibrium and relaxing compatibility, allows to enhance the stresses obtained from standard displacementbased finite elements schemes.
The key idea underlying RCP is to minimise the complementary energy associated to a patch of elements treated as an isolated system. On each patch, stresses are obtained as a linear combination of selfequilibrated stress modes enriched by an appropriate particular solution, and the explicit knowledge of the displacements is required only along the patch boundaries. Therefore, the RCP approach is naturally wellsuited for the virtual element schemes.
In this paper, the application of RCP in the context of VEM is analysed by proposing two approaches. In the first one, a degenerate patch composed of a single polygon is considered. For triangles, this does not significantly improve the results obtained using the classical strain computation. However, the benefit becomes apparent when the number of edges increases. In the second form, RCP is adopted by considering patches of elements in analogy to its original formulation developed in the context of finite elements schemes.
The paper is organised as follows. In Section 2, the VEM formulation for plane elasticity problems is briefly recalled. Then, the RCP procedure is introduced in Section 3. Numerical results are presented for a wide selection of test cases in Section 4 and, finally, conclusions are drawn in Section 5.
2 Virtual Element Method
Consider a body occupying a region of the twodimensional space on which a reference system is introduced. Let us denote with the symbol the boundary of such a body and consider the case of homogeneous Dirichlet boundary conditions. Other types of boundary conditions can be treated in the same spirit of finite element schemes. Indicating as , the displacement field, the associated strains are defined as
(2.1) 
where represents the differential compatibility operator. The plane elasticity problem can be thus expressed as: find such that
(2.2) 
where is the elastic matrix, is the space of the kinematically admissible displacements, is the applied external load.
In order to build a VEM scheme, it is assumed that is subdivided into nonoverlapping polygons with straight edges. In the following, indicates one and each of such polygons, denotes its boundary while denotes one and each of the edges of . Moreover, we here assume a homogeneous material, but inhomogeneous bodies can be treated as well (cf. for instance (18); (35)) We introduce a local approximation space, , whose elements are typically not known pointwise, but are determined by means of a suitable set of degrees of freedom (thus, is said to be virtual).
In the following, reference is made to a firstorder VEM discretisation scheme. The space is chosen so that
(2.3) 
where is the formal adjoint of while indicates the space of linear functions on the edges of .
A set of degrees of freedom which can be used in order to describe such a space is represented by the values of the displacements at the element vertices. It must be noticed that, thanks to the requirements expressed in Eq. (2.3), displacements must be linear on edges so that the vertex values completely characterise the displacement field on . A sketch of the virtual field and of the considered degrees of freedom is provided in Fig. 1.
Due to the virtual nature of the approximating functions, it is not possible to directly compute the terms of Eq. (2.2), and especially the lefthand side. It is thus necessary to introduce an operator, , defined as
(2.4) 
Above, denotes the space of constant fields over the element. It can be observed that is an operator which maps the virtual displacements onto constant strain fields. A typical displacement and strain field can be formally represented as
(2.5) 
where the vector collects the degrees of freedom (i.e. the values of the displacements at the polygon vertices) while is a vector collecting the degrees of freedom of the approximated strain field. The term collects virtual displacement shape functions and, thus, is not explicitly known, while , still considering a first order scheme, takes the form
(2.6) 
Indicating as the operator expressed as a matrix, it is possible to write the discrete strains as
(2.7) 
The operator can be thus evaluated by means of Eq. (2.4). In particular, substituting Eq. (2.7) into Eq. (2.4), we obtain
(2.8) 
which, integrated by parts, leads to
(2.9) 
where, denoting the outward normal as , is
(2.10) 
Equation (2.11) allows to compute as the solution of a linear system while matrices and are computable because displacements appear only at the boundary where they are explicitly known. This leads to the evaluation of as
(2.13) 
and, thus, to approximate the l.h.s. of Eq. (2.2) as
(2.14) 
After algebraic manipulations, it is possible to identify the stiffness matrix, , as
(2.15) 
In order to ensure that has the correct rank, a stabilisation term (not here further discussed) must be added to so that the stiffness matrix is usually redefined as . For further details, we refer to (2); (7), for example.
It is important to notice that the operator leads to constant strain fields while the virtual displacements are assumed to be linear over the element edges, and lead to divergencefree stress distributions (cf. Eq. (2.3) and Fig. 1).
For triangular elements, the virtual (loworder) formulation is equivalent to standard (loworder) triangular finite elements. However, when polygons with a higher number of edges are considered, the degrees of freedom for the displacements could, in principle, give rise to stress/strain fields richer than the constants. A kind of loss of information is thus expected when using operator to compute the stresses.
In the following, RCP is introduced as a postprocessing technique for the VEM procedure, with the aim of improving the accuracy of the obtained stresses.
3 Recovery by Compatibility in Patches
The key idea of RCP is to minimise the complementary energy associated to a patch of elements treated as an isolated system (19). Such operation allows to enforce equilibrium, which might be desirable in some applications (see for instance (36)), as well as to increase the accuracy of the obtained stress field. In this work, RCP is presented for degenerate patches composed of a single virtual element. The generalisation of the procedure to the case of patches composed of more than one element does not introduce substantial modifications.
The complementary energy associated to can be expressed as
(3.1) 
where is the unknown stress field while is the known displacement field (i.e. the displacement field computed by the displacementbased VEM analysis). It must be noticed that the displacement field appears only in the boundary term. In order to minimise , the stress field is written as the sum of two contributions:
(3.2) 
which correspond to its divergencefree part and an appropriate particular solution which guarantees equilibrium with the applied external loads. The homogeneous part is approximated as
(3.3) 
where is a matrix of preselected selfequilibrated stress modes, while is the vector of the unknown parameters. In particular, if a linear approximation is adopted for , takes the form
(3.4) 
The particular solutions, , appearing in Eq. (3.6) can be calculated as
(3.5) 
Above, and are the body force components. Furthermore, is an antiderivative of with respect to the coordinate , i.e. it holds . Analogous definition applies to with respect to the coordinate . If the antiderivatives for or are not explicitly known, we first approximate the body force within the element by a constant field, then we take the corresponding antiderivatives.
Imposing null variations of the functional expressed in Eq. (3.7), the following is obtained
(3.7) 
where
(3.8) 
In the definition of , see Eq. (3.8), the first term can be integrated by parts and, exploiting that the stress modes are divergencefree, we get
(3.9) 
In the finite element framework, the vector can be computed using the lefthand side or the righthand side of Eq. (3.9), leading to identical results. On the contrary, when virtual element schemes are considered, the lefthand side should be used, since the displacements are pointwise known only at the element boundary.
The adoption of a patch of elements instead of a single element does not introduce additional complexity: simply, the integrals appearing in Eq. (3.8) must be computed over the entire patch. Figure 2 provides examples of the patches adopted in the following numerical tests. In particular, Patch 0 indicates the degenerate case of a single element, Patch 1 is a patch obtained by considering all surrounding elements sharing at least one node with a central element, while Patch 1B represents a boundary patch of the same type as Patch 1.
4 Numerical results
In this section, numerical tests are used in order to validate the proposed stress recovery procedure. In particular, a square domain characterised by edge length equal to unity is considered. The body is assumed to be homogeneous and isotropic and its mechanical parameters are assigned in terms of the Lamé constants and (37); (7). Plane strain regime is assumed.
Indicating as the analytical stress field, the following error norm is used in order to evaluate the accuracy of the obtained results
(4.1) 
As it can be seen in Fig. 3, eight meshes are considered. Four of them are structured and are indicated with the letter ’S’ in the following. Those are composed of triangles, quadrilaterals, hexagons and a mixture of convex and concave quadrangles. The remaining four are unstructured meshes composed of triangles, quadrilaterals, random polygons and concave hexagons. Those meshes are indicated with the letter ’U’. The average edge length is indicated as .
Prescribed displacements are assumed and the corresponding body forces are calculated. Such body forces are then applied to the body, together with Dirichlet conditions over the entire boundary. Three tests are considered:

Test a: ;

Test b: ;

Test c: .
In particular, Test a is characterised by polynomial solution and engenders null body forces with inhomogeneous Dirichlet boundary conditions, while in Test b the displacement solution is trigonometric and is characterised by homogeneous Dirichlet boundary conditions, see (18). Analogously to Test b, Test c is characterised by homogeneous Dirichlet boundary conditions, and is a product of polynomial and trigonometric functions, while vanishes.
Three stress recovery procedures are compared and denoted as VEM, RCP0 and RCP1. In particular, VEM corresponds to the standard procedure based on the adoption of the operator , together with the use of the constitutive law. Instead, RCP0 and RCP1 denote the use of RCP as presented in Eq. (3.8), for Patch 0 and Patch 1, respectively.
Figures 4, 5 and 6 report results obtained for Test a, Test b and Test c for all the aforementioned meshes, respectively. It can be clearly seen that RCP always has a favourable effect, leading to more accurate results if compared to the use of the operator and, in some cases, to a substantial increase in the convergence rate. Regarding RCP0, we notice that the improvement can be clearly appreciated when the Hex (S) mesh is considered while it tends to disappear for triangular elements. Furthermore, considering now results obtained using RCP1, it can be seen that the use of element patches has a significant beneficial effect, in agreement with what has been observed in the context of displacementbased finite elements (19); (20).
Finally, a comparison between the three techniques in terms of von Mises equivalent stress distributions obtained using Hex (S) and Poly (U) meshes is reported for Test a and Test b. The exact solutions are shown in Fig. 7 while numerical results are reported in Figs. 8 and 9, respectively. The improved compliance with the analytical solution, even for relatively coarse meshes, can be clearly observed, confirming again the effectiveness of the proposed approach.
5 Conclusions
In the framework of displacementbased virtual element methods, the standard approach to compute stresses is via the constitutive law, using the available approximated strains. When general polygons are considered, such a procedure does not fully exploit the adopted degrees of freedom, thus leading to a quite poor stress field.
In the present paper, the beneficial effect of RCP used in connection with VEM discretisation schemes has been presented. The RCP formulation allows to compute accurate patchwise equilibrated stress fields starting from approximated boundary displacements and from the knowledge of the applied loads. Such a formulation appears to be wellsuited for the recovery of stresses in the framework of displacementbased virtual elements, in which such quantities are known only on the elements edges. In its simplest version, the RCP is applied elementwise, so representing an efficient alternative to the standard stress recovery. It has been also shown that further improved effects might be obtained by considering patches of elements, in agreement with the results obtained in the context of finite elements schemes. We finally remark that RCP can be very easily integrated in existing codes without introducing significant computational costs.
Footnotes
 journal: Int. J. for Num. Methods in Engrg.
References
 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., vol. 23, pp. 119–214, 2013.
 L. Beirão da Veiga, F. Brezzi, L. D. Marini, and A. Russo, “The Hitchhikers Guide to the Virtual Element Method,” Math. Models Methods Appl. Sci., vol. 24, no. 8, pp. 1541–1573, 2014.
 L. Beirão da Veiga, F. Brezzi, and L. D. Marini, “Virtual Elements for linear elasticity problems,” Siam. J. Numer. Anal., vol. 51, pp. 794–812, 2013.
 B. Ayuso, K. Lipnikov, and G. Manzini, “The nonconforming virtual element method.” Preprint arXiv:1405.3741. Submitted for publication.
 P. Antonietti, L. Beirão da Veiga, D. Mora, and M. Verani, “A stream function formulation of the Stokes problem for the virtual element method,” Siam. J. Numer. Anal., vol. 52, pp. 386–404, 2014.
 F. Brezzi and L. Marini, “Virtual Element Method for plate bending problems,” Comput. Methods Appl. Mech. Engrg., vol. 253, pp. 455–462, 2012.
 E. Artioli, L. Beirão Da Veiga, C. Lovadina, and E. Sacco, “Arbitrary order 2D virtual elements for polygonal meshes: part I, elastic problem,” Comp. Mech., vol. 60, pp. 355–377, 2017.
 E. Artioli, L. Beirão Da Veiga, C. Lovadina, and E. Sacco, “Arbitrary order 2D virtual elements for polygonal meshes: part II, inelastic problem,” Comp. Mech., vol. 60, pp. 643–657, 2017.
 P. Wriggers, B. Reddy, W. Rust, and B. Hudobivnik, “Efficient virtual element formulations for compressible and incompressible finite deformations,” Computational Mechanics, vol. 60, no. 2, pp. 253–268, 2017.
 H. Chi, L. Beirão da Veiga, and G. H. Paulino, “Some basic formulations of the virtual element method (VEM) for finite deformations,” Comput. Meth. Appl. Mech. Engrg., vol. 318, pp. 148–192, 2017.
 T. C., P. G.H., P. A., and M. I.F.M., “Polygonal finite elements for topology optimization: A unifying paradigm.,” Int. J. for Num. Meth. Engnr., vol. 82, pp. 671–698, 2010.
 P. Antonietti, M. Bruggi, S. Scacchi, and M. Verani, “On the virtual element method for topology optimization on polygonal meshes: a numerical study,” Computers and Mathematics with Applications, vol. 74, pp. 1091–1109, 2017.
 M. Benedetto, S. Berrone, S. Pieraccini, and S. Scialò, “The Virtual Element Method for Discrete Fracture Network simulations.,” Comput. Meth. Appl. Mech. Engrg., vol. 280, pp. 135–156, 2014.
 I. Perugia, P. Pietra, and A. Russo, “The plane wave virtual element method for the Helmholtz problem.” Submitted for publication.
 D. Mora, G. Rivera, and R. Rodríguez, “A virtual element method for the Steklov eigenvalue problem.” CI2MA PrePublicación 201427, in press on Math. Mod. Meth. Appl. Math., 2015.
 A. Gain, C. Talischi, and G. Paulino, “On the Virtual Element Method for ThreeDimensional Elasticity Problems on Arbitrary Polyhedral Meshes,” Comp. Meth. Appl. Mech. Engrg., vol. 282, pp. 132–160, 2014.
 F. Brezzi, R. Falk, and L. Marini, “Basic principles of mixed Virtual Element Methods,” Math. Mod. Num. Anal., vol. 48, no. 4, pp. 1227–1240, 2014.
 E. Artioli, S. de Miranda, C. Lovadina, and L. Patruno, “A stress/displacement virtual element method for plane elasticity problems,” Comp. Meth. Appl. Mech. Engrg., vol. 325, pp. 155–ï¿½174, 2017.
 F. Ubertini, “Patch recovery based on complementary energy,” Int. J. for Num. Meth. Engnr., vol. 59, pp. 1501–1538, 2004.
 A. Benedetti, S. de Miranda, and F. Ubertini, “A posteriori error estimation based on the superconvergent recovery by compatibility in patches,” Int. J. for Num. Meth. Engnr., vol. 67, pp. 108–131, 2006.
 O. Zienkiewicz and J. Zhu, “The superconvergent patch recovery and a posteriori error estimates. Part I: the recovery technique,” Int. J. Numer. Methods. Eng., vol. 33, pp. 1331–1364, 1992.
 O. Zienkiewicz and J. Zhu, “The superconvergent patch recovery and a posteriori error estimates. Part II: error estimates and adaptivity,” Int. J. Numer. Methods. Eng., vol. 33, pp. 1365–1382, 1992.
 B. Boroomand, O. Zienkiewicz, and J. Zhu, “Recovery procedures in error estimation and adaptivity. Part II: adaptivity in nonlinear problems of elastoplasticity behaviour,” Comput. Methods Appl. Mech. Eng., vol. 176, pp. 127ï¿½–146, 1999.
 G. Castellazzi, S. de Miranda, and F. Ubertini, “Adaptivity based on the recovery by compatibility in patches,” Finite Elements in Analysis and Design, vol. 46, pp. 379–390, 2010.
 C. Prange, S. Loehnert, and P. Wriggers, “Error estimation for crack simulations using the XFEM,” Int. J. Numer. Methods Eng., vol. 91, pp. 1459ï¿½–1474, 2012.
 T. Blacker and T. Belytschko, “Superconvergent patch recovery with equilibrium and conjoint interpolant enhancements,” Int. J. Numer. Methods. Eng., vol. 37, pp. 517ï¿½–536, 1994.
 N. Wiberg and F. Abdulwahab, “Error estimation and adaptive procedures based on superconvergent patch recovery (SPR) techniques,” Arch. Comput. Methods Eng., vol. 4, p. 203ï¿½242, 1997.
 T. Lee, H. Park, and S. Lee, “A superconvergent stress recovery technique with equilibrium constraint,” Int. J. Numer. Methods Eng., vol. 40, pp. 1139–1160, 1997.
 H. Park and S. Shin, “Element stress recovery technique with equilibrium constraints for superconvergence boundary stress extraction,” AIAA J, vol. 37, pp. 623–628, 1999.
 T. Kvamsdal and K. Okstad, “Error estimation based on superconvergent patch recovery using statically admissible stress fields,” Int. J. Numer. Methods Eng., vol. 42, pp. 443ï¿½–472, 1998.
 K. Okstad, T. Kvamsdal, and K. Mathisen, “Superconvergent patch recovery for plate problems using statically admissible stress resultant fields,” Int. J. Numer. Methods Eng., vol. 44, pp. 697ï¿½–727, 1999.
 B. Boroomand and O. Zienkiewicz, “Recovery by equilibrium in patches (REP),” Int. J. Numer. Methods Eng., vol. 40, pp. 137ï¿½–164, 1997.
 B. Boroomand and O. Zienkiewicz, “An improved REP recovery and the effectivity robustness test,” Int. J. Numer. Methods Eng., vol. 40, pp. 3247–3277, 1997.
 B. Boroomand, O. Zienkiewicz, and J. Zhu, “Recovery procedures in error estimation and adaptivity. Part I: adaptivity in linear problems,” Comput. Methods Appl. Mech. Eng., vol. 176, pp. 111ï¿½–125, 1999.
 E. Artioli, S. de Miranda, C. Lovadina, and L. Patruno, “A family of virtual element methods for plane elasticity problems based on the hellingerreissner principle,” Comp. Meth. Appl. Mech. Engrg., submitted.
 S. de Miranda, L. Patruno, and F. Ubertini, “Transverse stress profiles reconstruction for finite element analysis of laminated plates,” Comp. Stru., vol. 94, pp. 2706–2715, 2012.
 L. Beirão Da Veiga, C. Lovadina, and D. Mora, “A Virtual Element Method for elastic and inelastic problems on polytope meshes,” Computer Methods in Applied Mechanics and Engineering, vol. 295, pp. 327 – 346, 2015.