An engineering perspective to the virtual element method and its interplay with the standard finite element method

In this manuscript we thoroughly study the behavior of the virtual element method (VEM) in the context of two-dimensional linear elasticity problems for an engineering audience versed in standard FEM. Through detailed convergence studies we show the accuracy and the convergence rates recovered by VEM, and we compare them to those obtained by the h - and p -versions of the finite element method (FEM). We also demonstrate the mixability between FEM and VEM; in particular, applications of VEM for coupling non-conforming discretizations and for local refinement are presented, showing higher versatility compared to FEM. Computer implementation aspects in displacement-based finite element codes are thoroughly explained, remarking on the main differences with respect to standard FEM.


Introduction
The Virtual Element Method (VEM) is a numerical technique designed for polytopal discretizations (polygons in 2-D or polyhedra in 3-D).Its mathematical foundations have been well established and its benefits have been explored on numerous applications.However, some further insight into the inner-workings of the method is still required for the application-oriented user.The purpose of this paper is to give an engineering perspective to VEM, showing with practical applications not only its features but also the synergic advantages that arise when the method is blended with the standard finite element method (FEM).
Inspired by the ideas put forward by the Mimetic Finite Difference method [1,2], VEM was devised and put into a Galerkin framework [3], ultimately becoming a generalization of standard FEM that allows using arbitrary polytopal meshes.The main idea behind VEM is that the local shape function space in each element is defined implicitly-thus the name virtual.Furthermore, shape functions are constructed according to an integer k, which represents the interpolation order of the approximation.At variance with the p-version of the finite element method ( p-FEM) [4], however, VEM may contain non-polynomial functions.In spite of its short life span, VEM has already caught the attention of the computational mechanics community; the method has been investigated in the context of a wide number of problems, including small strain [5][6][7] and finite strain [8] linear elasticity, material inelasticity [9,10], topology optimization [11], plate bending [12], contact [13], geomechanics [14], and fracture mechanics [15].High-order VEM has been applied to scalar problems and has been developed theoretically for arbitrary polynomial interpolation k [16][17][18].In the present context of linear elastostatics, a C++ implementation and description for the lowest order case can be found in [19] and for an arbitrary order formulation see [20].
The method was designed to handle arbitrarily shaped polytopal discretizations, which are particularly convenient for meshing complex geometries.This endows VEM with unmatched flexibility in the pre-processing stage of the analysis.This is by far the most appealing feature of the method, not only because VEM has been shown to be robust to mesh distortion [18,21], but also because non-conforming discretizations in standard FEM-i.e., meshes containing hanging nodes-can be handled in VEM in a straightforward manner; VEM simply interprets a hanging node in a non-conforming mesh as a division of an edge at a 180 • angle, effectively increasing by one the number of edges in elements sharing this node.This requires no modification whatsoever to the VEM formulation and the element is treated normally, thereby resulting in a VEM-conforming mesh.Handling non-conforming discretizations gives VEM an edge over standard FEM, since creating conforming discretizations in the latter is a time-consuming process that may even result in inappropriate meshes for the analysis; remeshing may introduce elements with bad aspect ratios that can affect the accuracy of the solution and/or degrade the condition number of the system matrix [22,23].
We believe that the many innovative features of VEM can be buried beneath mathematical formality and bogged down by some of its technical complexity.For this reason, in this manuscript we thoroughly study VEM for linear elastostatics problems in the hopes of making the method more accessible to the FEM practitioner, who is interested in it but is otherwise discouraged by VEM's seemingly steep learning curve.We discuss the method's advantages and disadvantages and compare them to standard FEM.We further compare VEM to p-FEM in the context of high-order interpolations for a smooth problem.We provide detailed convergence analyses that show VEM is not only optimally convergent, but also converges faster than p-FEM.We investigate VEM as a procedure to facilitate standard FEM analyses.To this end, a hybrid VEM-FEM approach is used to take advantage of VEM's main virtue in dealing with non-conforming discretizations.We finally exploit VEM for an efficient refinement strategy in regions of high stress concentrations, where standard FEM would mandate at best for a complex strategy for handling hanging nodes, or at worst for a complete remeshing.

Problem statement
We are interested in solving the standard linear elasticity boundary value problem on a continuous domain Ω ⊂ R 2 , as shown in Fig. 1a.A body force per unit area f acts on Ω , and its domain boundary ∂Ω (with unit outward normal n) is decomposed into disjoint regions ∂Ω u and ∂Ω t ; these are, respectively, regions with prescribed essential (Dirichlet) and natural (Neumann) boundary conditions.For simplicity we assume homogeneous boundary condition ū = 0 on ∂Ω u .An arbitrary traction field t is prescribed on ∂Ω t .

VEM discrete equations
In VEM, the (Bubnov-)Galerkin or finite-dimensional form of (1) is obtained by choosing an appropriate discrete space V h ⊂ V for the trial solution u h and weight function v h .Unlike in standard FEM, polygonal elements of arbitrary shape can be used-both convex and non-convex-and thus the domain Ω is first partitioned into nonoverlapping polygonal regions ) , as illustrated in Fig. 1b.In order to deal with arbitrarily shaped polygonal elements, VEM uses general finite element spaces that may contain non-polynomial functions.However, as for standard FEM, numerical integration is required for the polynomial component and hence we assume that it is possible to compute their contribution by means of some quadrature procedure.One such technique, although by no means the only possibility, is to subdivide each polygon into triangles and use Gaussian integration.However, an appropriate choice of DOFs could drastically decrease explicit numerical integration, as explained later.

Projection operator. The local projection operator is the most crucial
] 2 indicates the local space of VEM shape functions on E (formally defined in Appendix A.1) and P k (E) ≡ [P k (E)] 2 , with P k (E) denoting the polynomial space of order less than or equal to k defined on element E. The projector is used to obtain the polynomial projection of VEM shape functions without losing accuracy in the process.This is done because VEM shape functions cannot be used directly-they are not explicitly known in general.Unless stated explicitly, it is implied hereafter that Π ≡ Π E,k to simplify the notation.
The operator is constructed locally in each polygon so that it satisfies the following orthogonality condition: where integration is performed on E (indicated by the subscript) and zero energy modes (infinitesimal rigid motions) are fixed for uniqueness of the projection.This equation manifests that the error between u h and Π u h in E, measured in the energy norm, cannot be captured by the polynomial space P k used, i.e., the error in the energy norm is orthogonal to P k .In VEM lexicon, this property is referred to as k-consistency.

Discrete variational form. The discrete bilinear form on an element
where the first term in the right hand side is known as the consistency part and is no different than what would be obtained by following standard procedures to obtain a finite-dimensional form of (1)-except for the use of the projection operator.The second term, which is scaled by a user-defined parameter τ h , stabilizes the method.An approximation L E (v h ) ≈ L E (v h ) of the local linear form is needed to compute the contribution of body forces (further explained in Section 2.5).The total contribution to the internal energy and the load term are obtained as respectively.With these considerations in mind, the discrete form of the problem expressed in abstract form is: Find A stabilization, which computes an approximation of the strain energy not taken into account in the projection, only needs to scale appropriately, i.e., it should be of the same order of magnitude as the consistency term regardless of the mesh.More specifically, the stabilization is required to preserve the coercivity of the problem so that â ( u h , v h ) is bounded below and therefore can be solved (and the solution is unique).In other words, the second term in ( 5) is needed because functions within each element are only polynomials when dealing with k = 1 and simplexestriangles (tetrahedra) in 2-D (3-D); for k > 1, a virtual element has more DOFs for triangles but less or equal for quadrilaterals.In general, for polygonally shaped elements, functions inside them may no longer be polynomials and the stability component controls the strain energy of terms of order greater than k.This, along with the fact the VEM local space contains all polynomials of order up to k, ensures that VEM passes the patch test of order ≤ k − 1 no matter the shape of elements used.In particular, virtual elements of any order k > 0 pass the zeroth order patch test, i.e., the local VEM shape functions are able to represent constant strain states.

Degrees of freedom (DOFs)
In each polygon of the discretization, a vector-valued test function in the local VEM space, v h ∈ V h (E), has the following properties: • components of v h are continuous functions along the polygon boundary ∂ E; • v h is a vector with polynomial components of degree k on every edge of e; • ∆v h ⏐ ⏐ E are polynomials of degree k − 2 in the interior of E.1 Although these properties allow us to recover the behavior of v h over the boundary ∂ E, there is no information about v h inside E-except for a condition on its Laplacian.In fact, DOFs associated to a polygon are defined in such a way that explicit knowledge about the shape functions is not actually needed-which is emphasized by the virtual word in the name of the method.Given a polynomial order k, DOFs associated to a given polygon are (schematically shown in Fig. 1b for k = 2): For k ≥ 2, the moments up to order k − 2 of v h in E: For DOFs of type (i) and (ii), an appropriate choice is the (k + 1)-point Gauss-Lobatto quadrature points on each edge as they simplify the computation of integrals.Another choice could be to replace DOFs of type (ii) with k − 1 "moments" on the edge.In fact, any way of univocally determining a polynomial of order k on an edge can be used, although the one presented here closely resembles the choice for standard FEM while the use of moments on an edge is closer to p-FEM.For a given polygon, the total number of DOFs along one spatial dimension is 1)/2; thus for 2-D elasticity a polygonal element has 2n d DOFs.At this point it is convenient to define operator dof i acting on v h ∈ V h (E), i.e., dof i (v h ), i = 1, . . ., 2n d , as the value of v h at the ith degree of freedom.Basis functions for the local VEM space along one dimension will be denoted by {ϕ i } i=1,...,n d and are defined so that the Kronecker-δ property holds, i.e., dof i ( ϕ j ) = δ i j .The vector that collects the shape functions is denoted by . A graphic representation of functions that associated with each DOF type is presented in Fig. 2, where functions that correspond to edge and vertex DOFs have a null mean value, i.e., their value at the internal DOF is null.The contrary happens to a function associated with an internal DOF, whose mean value is the area of the element and is identically zero on the boundary.These functions are defined implicitly by a local boundary value problem and therefore a numerical approximation was computed for illustration purposes.As we are dealing with problems in 2 dimensions, the bold symbol ϕ will stand for a vector-valued expression of the basis functions, such that The choice of DOFs univocally defines a polynomial of order k on any edge of the mesh.Therefore, contiguous elements sharing DOFs across a common edge define an inter-element continuous function so that conformity is preserved and the solution is globally C 0 -continuous across the domain.For C 1 VEM discretizations see [24,25].

Local stiffness matrix
In this section we show the construction of the element local stiffness matrix k E .In standard FEM, it is usually written as where C is the material constitutive tensor for either plane stress/strain and is the strain-displacement matrix, with ∂ denoting the differential operator (in Voigt notation) In VEM, because the functions are unknown within the element, constructing the stiffness matrix is possible thanks to the local DOFs and the aforementioned projection operator.The stiffness matrix is obtained by the contributions of consistency and stability components: where since only an approximation of the shape functions can be recovered using Π (ϕ i ≈ Πϕ i ).Nevertheless, this approximation will retain optimal convergence properties since the consistency part is exact up to order k, i.e., it can compute the exact strain energy associated with polynomial strains of order up to k − 1.While the consistency part has to be computed exactly, its expression differs from that of standard FEM, as explained next.
Consistency term.For a given admissible vector-valued displacement field ] ⊺ , the (engineering) strains associated with it can be computed using the strain operator: The projection operator Π , defined by (4), allows the computation of the consistency term of the stiffness matrix by first expressing the projection of a VEM function in the local polynomial base P k = { p α } α=1,...,n k ∈ P k (E) (see Appendix A.1).To wit, Next, the component-wise expression can be obtained by using linearity of the strain tensor: The matrix representation of the projection operator Π is given by where components of matrices B and G are given by respectively (see Appendices B.1 and B.3).The projection operator given by ( 14) is crucial to understanding the method; (15) and ( 16) can both be seen as internal virtual work produced on element E by a state of stress Cϵ( p β ) with respect to virtual strain fields ϵ(ϕ i ) and ϵ( p α ), respectively.Therefore, if we rewrite Eq. ( 14) as G Π = B, we see that the function of the projection is to ensure the internal virtual work produced by VEM functions is the same as that caused by the polynomial basis chosen.Finally, the consistency part of the stiffness matrix is computed as where G equals G but with the first three rows containing all zeros.In essence, we decompose the VEM shape functions into polynomials, and then use G.To better understand (17), we reformulate it into something more akin to (9) and thus more familiar to standard FEM; replacing G and using the strain operator (11) yields where B3×n k is defined as: Note that since p 1 , p 2 and p 3 are taken as infinitesimal rigid motions (A.3), there is no strain associated with them and thus the first 6 columns of B are zero.
It is important to note that G is computable while B is (a priori) not, since it involves integration of a VEM shape function inside the element.However, integration by parts allows to compute the internal energy when at least one of the terms is a polynomial, as it is the case here (see Appendix A.2).
Stability term.This term is not present in classical FEM discretizations and requires some preliminaries.The dof operator defined in 2.3 can also be applied to polynomials p ∈ P k (E) since P k (E) ⊂ V h (E).The matrix D 2n d ×n k collects the values of the polynomials p α at DOFs on E (see (A.1)), and is defined as Bear in mind that the value at the degree of freedom can be a simple evaluation of the function, in the case of point values, or an integration, in cases where a DOF is a moment on a length or an area.Note that for vector valued functions, the DOF operator can refer to different components of the function.The matrix representation of the projection operator acting from V h (E) into itself can then be written as The stability term can now be computed as where, as stated earlier, τ h is a user-defined parameter that can be taken as 1/2 for linear elasticity (see [20] for an analysis of this parameter in linear elasticity and [10] for inelastic solids), tr • denotes the trace operator, and I is the 2n d × 2n d identity matrix.An inspection of ( 22) reveals that the stability term is basically a rough approximation of the internal energy associated with the difference between a VEM shape function and its projection.Namely, the term scales the term with respect to the consistency part, which is the desired property used for proving convergence of the method [3].Since all polynomials of degree up to order k are included in the local VEM space, these terms in the stabilization are of higher order than the chosen discretization order k.For a similar presentation with explicit matrix construction and arbitrary interpolation order k see [20].

Local load vector
From (3), the local VEM loading term for a single polygonal element takes the form where in addition to the body force and boundary tractions, we added the contribution of concentrated forces F i widely used in engineering practice-and which could be mathematically described by means of Dirac delta functions.
The last two terms in (23) do not require any special treatment in VEM and thus follow the standard FEM approach.In particular, tractions on the Neumann boundary can be computed directly as in standard FEM since VEM local shape functions on elements' boundaries are polynomials of degree k and are thus known explicitly.Also, concentrated loads are associated with the displacement DOFs of the location where the load is applied, and since nodal DOFs for VEM and FEM are analogous, so is the handling of concentrated loads.In the case that this load is associated with one of the DOFs of type (iii) in ( 8), the same procedure as for generalized (non-Lagrangian) DOFs in FEM can be applied.Non-zero prescribed values at the Dirichlet boundary can be dealt with by obtaining the equivalent traction forces that the prescribed displacements generate, as in standard FEM.Alternatively, by using "lift functions" the problem can be solved as homogeneous with an added load term due to the prescribed displacements, which is again standard FEM practice.
The first term in (23) requires special consideration since it involves the VEM shape function in the interior of the element, and therefore the projection operator is needed to express the VEM displacement v h as a polynomial within E in order to perform the integration.We thus follow a similar procedure to that used for the construction of the stiffness matrix, and we distinguish between three cases based on the accuracy order k.For k > 2, we introduce the (orthogonal E) so that the loading term due to body forces can be approximated as Note that, although projection is done on polynomials of order k − 2, this is enough to retain optimal convergence properties as the approximation of the load vector and that of the discrete solution are of the same order2 .For k = 2 we use the L 2 projection Π 0 k onto the space of polynomials P k (E), so that the loading term becomes Finally, in the case of k = 1, the computation is done taking an average of f over the element and using a vertex-based integration rule for v h .More details are given in Appendix C.An alternative procedure can be found in [20], where the projection is taken on the body force f .

Assembly of local arrays
Once that local stiffness matrices and local load vectors are calculated for all elements in the discretization by using (10) and (23), respectively, these are assembled into the global stiffness matrix and load vector as where A denotes the standard finite element assembly operator.As it is commonly done in standard FEM, a system of linear equations K U = R is solved for the global degree of freedom vector U after prescribing boundary conditions.It is worth noting that the condition number of the prescribed stiffness matrix depends on the choice of the polynomial basis used; the monomial basis used here erodes the conditioning as the order k of the interpolation increases.This undesired behavior can be alleviated by choosing an orthogonal polynomial basis [26].

Implementation
In this section we discuss, with the aid of pseudo-code, the implementation of VEM in a standard displacementbased FE package.Implementing VEM requires the following modifications or additions: (i) Input/Output (I/O) Regarding input, VEM requires discretization formats that are not standard, as meshes can be comprised of polygons of an arbitrary number of vertices-and all polygons need not have the same number of vertices; this implies that the same data structure used for storing elements in standard FEM may not be reused.One way to accommodate for new mesh data is to set a separate data structure for the VEM discretization that uses an associative array (map) where the key is the number of polygon nodes and the value is the corresponding connectivity array (element freedom table).An important remark on polygonal meshes is that, in general, it is not straightforward to produce discretizations whose elements have a specific number of vertices.For instance, meshing procedures based on Voronoi cells usually result in discretizations where most of the elements are hexagons, with some pentagons or quadrangles near the boundaries.Typically, the simplest meshes to produce are those containing elements with 3, 4, or 6 vertices.However, depending on the nature of the problem, elements of different number of vertices may arise and this is where the mesh versatility offered by VEM becomes an asset.
Regarding output, some post-processing is needed to write the results in a format that is commensurate with the use of polygons of arbitrary shape.Many post-processors support the graphical representation of polygonal data by providing them with field data only at vertex locations; inside polygons the field is then interpolated from vertex data.This output type is straightforward to achieve in a VEM implementation: Because in VEM the DOFs associated with vertices correspond to the value of the field, ordering all vertex DOFs first would require minimal post-processing because the solution vector can be simply sliced.This approach may not be sufficient to properly display the field within the polygons, particularly for higher-order interpolations.A better representation of the field can then be obtained by creating an "output mesh", where the field value is not only given at polygon vertices, but also internally by evaluating Πu h (x) at various locations within the polygons.The post-processing here is not trivial and takes some time to produce [27].(ii) Local arrays The other main addition to a standard FEM package lies in the computation of local stiffness matrices and local load vectors.The routine does not only compute the arrays depending on the value of the interpolation order k, but also returns the corresponding DOF array so that the contributions can be readily assembled into the global sparse arrays.Pseudo-code for the assembly routine is given in Algorithm 1.
All other parts of a FEM package, from the data structures used to store the global stiffness matrix and the load vector, to the solver used to obtain the solution vector U, need no further modification.However, if a hybrid FEM-VEM code is sought, further modifications are needed.Detecting a non-conforming FEM interface and placing polygons to create a VEM conforming discretization requires a set of geometric operations that could be encapsulated into a geometric engine.Standard FEM elements that share nodes along the non-conforming interface could be removed altogether from the discretization, or simply masked so they are not taken into consideration in the numerical quadrature of local arrays.The procedure then detects non-conforming nodes in the mesh and embeds them into polygons that replace the removed/masked elements.Noteworthy, no new DOFs are added when dealing with linear interpolations and thus efficiency is not compromised when solving the global system matrix.
It is worth highlighting that there are some additional operations required for computing the stiffness matrix in VEM with respect to FEM.The computation of the projector operator is of course not present in FEM since shape functions are already polynomials; therefore, VEM requires a higher number of computations to obtain the local stiffness matrix.However, many of these added operations are not computationally intensive and numerical quadrature can be significantly reduced with a smart choice of DOFs.For instance, the computation of matrix G can be directly avoided for the computation of the projector.Matrix B can be obtained by 1-D integration on element boundaries simplified by the choice of Gauss-Lobatto integration points, and the internal integrations are directly obtained from the DOF definitions.Matrix D requires the evaluation of polynomials at specific locations as well as integration of polynomials over the element, which ultimately requires numerical quadrature (see the appendix for more details).Additionally, the stabilization part involves some matrix operations as well.The load vector also requires the computation of a projector for the body load for k ≥ 2, while boundary tractions and concentrated loads are handled as in standard FEM.Finally, the assembly process is basically the same, except for the fact that VEM elements of the same order do not necessarily have the same number of DOFs if they have different number of vertices.

Algorithm 1 Computation of local matrices and load vector due to body forces
Input: Interpolant order k, polygonal element E with nodal coordinates X = [ x i x j . . .] ⊺ , quadrature point weights γ and coordinates ξ , constitutive matrix C and body force vector f -get centroid, area, diameter and number of vertices (15) and Appendix B.1 -compute local load vector due to f (C.17)

Results
In this section we showcase VEM with the aid of numerical tests.We study the accuracy of VEM through a series of patch tests, and then investigate VEM both as a standalone method and as a helper method to standard FEM for the coupling of non-conforming discretizations.

Constant strain patch test
We first test whether VEM is able to recover a constant state of strain through a simple patch test, since it is a necessary condition for guaranteeing convergence [28].The problem, which is schematically shown in Fig. 3, consists of a unit square plate Ω = [0, 1] × [0, 1] mm 2 constrained as shown and loaded along its right edge with a traction t = 1 N/mm e 1 .The plate's material has Young modulus E = 1 MPa and Poisson ratio ν = 0.3.Plane strain conditions are considered.
Since the exact solution for this problem is a linear displacement field, first order VEM is enough to pass the test.We thus consider k = 1 on two discretizations.While in the first all elements are VEM polygons, in the second we further test the mixability between one concave VEM polygon and standard FEM elements-specifically six linear triangular (T3) and four bilinear quadrangular (Q4) elements.The results are shown in Figs.4a and 4b, respectively, for the discretizations thereof.The figures show, on the deformed configuration, that both discretizations recover

High-order patch test
The traditional patch test studied above aims at recovering a constant state of strain.However, this test does not shed any light into the behavior of VEM for higher-order interpolations.Here we study a high-order patch test, where we manufacture a solution with k > 1 and see whether we can capture it with error within numerical precision.In fact, it can be shown that an order k VEM discretization automatically satisfies patch tests of orders 1, 2, . . ., k − 1.This is due to the fact that a VEM space of order k contains all polynomials up to degree k.The constant strain patch test was kept for completeness.
Consider the following displacement field: with a 6th-order polynomial in the u 2 component, whose strain is therefore a 5th-order polynomial.We solve the problem on the square plate of Fig. 3 under plane stress conditions, and material properties E = 1 MPa and Poisson's ratio ν = 0.2.A pseudo-body force, obtained as f = −∇ • σ (u), is also applied.Fig. 5 shows the numerical solution for this test, where the magnitude of the displacement field is plotted on the deformed configuration.As expected, with k = 6 VEM was able to obtain this solution with an error in the L 2 norm O ( 10 −12 ) and an error in energy norm O ( 10 −11 ) (these are defined in Section 4.2).Noteworthy, the apparent curvature in the right edge of the square is the result of describing the deformed shape also considering the k − 1 internal points for each polygonal element edge in the discretization.

Non-conforming patch test
We now study a patch test where a non-conforming FEM discretization is resolved using VEM polygons.In standard FEM non-conforming discretizations can be resolved by either remeshing locally to create a conforming mesh, or by adding constraint equations to ensure the continuity of the field at hanging nodes.This process is usually not computationally efficient nor robust, since local remeshing could generate distorted elements that may degrade the accuracy of the method.Furthermore, if the mesh sizes at both sides of the non-conforming interface are considerably different, remeshing may no longer be localized to the interface due to the need for transitioning elements on the coarser side.
For this type of applications VEM has a clear edge over standard FEM because hanging nodes can be easily incorporated into VEM polygonal elements; these are only used along the non-conforming interface.As a result, VEM can be seen as a helper procedure for extending the capabilities of standard FEM.
The problem is schematically shown in Fig. 6a, where the domain Ω is discretized by T3 and Q4 elements, respectively, at the left and right of a non-conforming interface Γ .Tractions t = 1 N/mm e 1 and t = 1 N/mm e 2 are prescribed at the right and top edges of the domain, respectively.We use plane strain conditions, E = 1 MPa and Poisson's ratio ν = 3/8.The resulting von Mises stress field is shown on the deformed configuration in Fig. 6b, where it can be seen the stress is constant throughout the plate.

Convergence analysis on a smooth problem
We now focus our attention to the convergence behavior of VEM, i.e., we look at how fast the error with respect to an exact solution decreases as the mesh size h → 0 (h-convergence) or as the interpolant order increases for a fixed mesh size h ( p-convergence).Consequently, for each study, VEM will be compared to the hand p-versions of FEM, respectively.
Consider the smooth displacement field given by u = sin (4π x 1 ) sin (4π x 2 ) (e 1 + e 2 ) , on the unit square domain Ω shown in Fig. 3. Either component of ( 28) is illustrated in Fig. 7a, and the types of discretizations for VEM and p-FEM, respectively, are shown in Figs.7b-7c.The material constants for this plane stress problem are E = 1 MPa and ν = 3/8.Homogeneous essential boundary conditions are prescribed along the entire perimeter of the domain.In addition, a pseudo-body force f = −∇ • σ (u) is applied on the entire plate.The error of our VEM approximations with respect to (28) is quantified through two standard measures, although the projection of the discrete solution is used in the computations (as usually done when working with VEM):  • Relative error in L 2 norm: • Relative error in energy norm: The h-convergence results, with discretizations composed of 64, 256, 1024, and 4096 polygons, are summarized in Fig. 8; Figs.8a-8b show the relative error in L 2 norm as per (29) as a function of the total number of DOFs N dof 8(a) and mesh size h 8(b).Similar results are presented in Figs.8c-8d for the relative error in energy norm given by (30).As apparent from the figures, VEM attains optimal convergence rates for this smooth problem for all interpolation orders.It is worth noting that in order to get these convergence rates, the load vector is computed differently for k = 1 and k = 2, as explained in Section 2.5.
Figs. 9a-9d show the p-convergence results in L 2 and energy norms; we compare VEM results with those obtained by p-FEM on a structured mesh of Q4 elements.We therefore fix the mesh size h (see Figs. 7b-7c for the discretizations used) and increase the order of the interpolant.The figures show that, similarly to p-FEM, VEM recovers exponential rates of convergence of the form where C, γ , θ are constants.In Figs.9c-9d we show the results of applying the natural logarithm twice to (31), resulting in the equation of a line in ln N dof × ln (ln ∥e∥) space, i.e., ln (ln ∥e∥ E ) = C 1 − θ ln N dof .It is shown that VEM has a steeper slope θ than p-FEM; θ = 1.02 for VEM vs. θ = 0.9 for p-FEM in energy norm (theory predicts θ ≥ 1/2).VEM and p-FEM therefore perform similarly for smooth problems.The difference in accuracy observed is due to the lower accuracy of the polygonal element type used with respect to the number of DOFs.Moreover, the difference in slope reflects the type of approximation spaces used.For instance, VEM on triangles has more DOFs (except for the lowest order), so the slope is expected to be better for p-FEM.On the other hand, on quadrilateral meshes VEM uses less DOFs for k > 2 while still retaining same convergence properties, so the slope in this case should be higher for VEM.

Stress convergence
In this example we look at the classical problem of an infinite plate with a circular hole of radius a = 4 mm, as shown in Fig. 10a.The plate is subjected to a uniaxial tension of magnitude σ 0 = 1 N/mm at infinity.Due to symmetry, we only analyze the portion shown in the inset under plane strain conditions, where w = 20 mm.The exact traction t = σ • n is applied to top and right edges.The material of the plate is linearly elastic with E = 200 GPa and ν = 0.3.This problem has an analytical solution [29]: with associated stress field In this problem we study the convergence in stress using the L 2 norm of the stress error: which is homologous to the error measure defined in (30), and the stresses are computed from the projection of the discrete solution as σ h = Cϵ(Π u h ).The h-convergence is performed using meshes with 100, 400, 1600 and 6400 polygonal elements.We also compared the VEM solutions with standard FEM using constant strain triangular elements; uniform meshes composed of 167, 492, 1869 and 7497 T3 elements were used.Typical VEM and FEM discretizations used in this example are shown in Figs.10b and 10c, respectively.
The convergence results are summarized in Fig. 11, where we appreciate similar convergence behavior in stress for both VEM and FEM.Because of the stress concentration at the hole, better predictions could be obtained by using geometric meshes [30], where the mesh size is reduced in geometric progression towards the stress concentration.Later in this section we will explore an alternative approach where VEM is used to study stresses in regions with singularities.

Coupling non-matching meshes
Here we build on the non-conforming patch test of Section 4.1 and show how standard FEM can be coupled to VEM for solving a complex non-conforming mesh; as explained earlier, VEM polygons will be constructed only where a non-conforming node is detected.
We revisit the example of Fig. 10a for the 2-D stress analysis of an infinite plate containing a circular hole loaded in uniaxial tension.We choose x 2 = 0 to place our non-conforming interface and again we take advantage of symmetry to model half of the problem.The discretization is composed of T3 elements for x 2 ≥ 0 and Q4 elements for x 2 ≤ 0. We prescribe the exact value of the traction field at the domain's top, bottom, and right edges.
The results are shown in Fig. 12, where the stress in the e 1 direction is shown on the deformed configuration.A stress concentration factor K t = 2.98 is obtained, which approaches the well known theoretical value K t = 3 obtained by Kirsch [31].The stress concentration value obtained is, of course, closer to the theoretical prediction in the upper mesh due to the use of finer elements.In the next example we use this fact to improve the approximation properties close to a singularity.

VEM for localized mesh refinement
This section deals with the solution of the Flamant problem [32,33], a half-space loaded by a downward vertical force, as illustrated in Fig. 13a.Its interest lies in studying stresses in the neighborhood of the applied load considering that there is a singularity at the origin.This problem has applications in foundation design, where stresses in the soil (considered as a semi-infinite elastic medium) are sought.The exact solution for the displacement  field in the radial and angular directions is given by: where F y is the magnitude of the applied point load, G is the shear modulus, and r = √ x 2 1 + x 2 2 , θ = arctan (x 2 /x 1 ) are polar coordinates with origin located at the load's point of application.The exact stress field is given in polar coordinates by or in Cartesian coordinates by Inspecting (35) reveals that, for a fixed θ , displacements are directly proportional to log(r )-and are therefore undefined at r = 0-but also that displacements are unbounded for r → ∞.Finite element solutions of this problem show that indeed ∥u (0)∥ becomes increasingly large as h → 0. Furthermore, by simply integrating (36) in a neighborhood of the applied load we can infer that the solution to this problem does not have finite strain energy (the H 0 -norm of the stress is analogous to the energy norm); in fact, it can be easily shown that the solution is not in Despite all that, the finite element solution for this problem is well behaved if we exclude the element that contains the load, and we will investigate the stresses obtained by our VEM approximation.The computational domain consists of the square region shown in Fig. 13 with a = 1 m, so we take advantage of symmetry.At the bottom and right sides of the domain we prescribe the displacement field given by (35).We will then measure stresses at a horizontal line located at x 2 = 0.1 m (distance marked b in Fig. 13a).By taking full advantage of VEM's ability to handle arbitrary discretizations, a local refinement scheme can be devised to increase the number of elements only in the area of interest, in this case close to where the load is applied.The discretization we used is shown in Fig. 13b, where increasingly smaller polygons are placed close to the applied load.
The results are then summarized in Fig. 14, where the vertical stress σ 22 and the shear stress σ 12 obtained by the VEM approximation are compared to the exact values given by (37).VEM stresses were averaged in each element, as shown by the contour plots in the figure.Because stresses close to the load get very high values, the contour plot for σ 22 was capped to a minimum value of −5 Pa.Similarly, the stress contour for σ 12 was capped to a maximum value of 1 Pa.Along the line x 2 = b these averaged stress values are reported at the center of each element and compared to the exact values.The figures show good agreement between the two, and thus VEM can be used as a powerful technique for creating fine meshes close to points of interest in our domain.It is worth noting that all elements in the discretization chosen are VEM polygons.Nevertheless, as shown earlier, we could have used a standard FEM discretization and use VEM elements only close to the applied load-with appropriate transition VEM polygons.

Summary and conclusions
In this manuscript we studied high-order VEM in the context of 2-D elastostatics.We discussed the construction of the method, highlighting the main differences with respect to a standard displacement-based FEM formulation.Noteworthy, VEM defines a virtual space, which is composed of trial/test functions whose explicit knowledge is not required.By building a projection operator it is then possible to obtain the necessary information about those functions in the standard polynomial space needed for the computation of the stiffness matrix.The construction of the virtual space is performed element by element and requires only mesh geometry, allowing the framework to be used regardless of the type of polygons-including non-convex shapes.However, the computation of the projection operator can be computationally expensive, in particular when the displacement field is approximated with high-order functions.
VEM's accuracy and convergence properties were compared to those of hand p-versions of FEM.VEM guarantees the same performance but a higher robustness for arbitrary polygonal discretizations.With iso-parametric formulations, commonly used in standard FEM, a reference element is defined on which all computations are performed and then transformed to the actual position in space.Because the mapping between reference and actual elements involves the Jacobian of the transformation and badly shaped elements result in loss of accuracy, the performance of the method becomes dependent on mesh quality.Thus, for geometrically complicated problems, VEM provides a clear advantage in this regard.The numerical tests conducted showed not only that VEM's accuracy is on par with that of FEM, but also that optimal convergence rates are achieved for high-order approximations of the displacement field.VEM is therefore a viable alternative to FEM in problems with intricate geometries, where creating a matching mesh would require too much effort.
We showed how to exploit the advantages arising from this polygonal-mesh based technique for coupling nonconforming meshes.A fine VEM discretization was also used to improve the approximation properties near regions where the solution is no longer smooth.Therefore, integrating VEM in standard FEM packages results in a powerful framework for handling hanging nodes and local refinements.To that end, we gave detailed instructions for a straightforward computer implementation.At a glance, we have summarized the main differences and similarities of the two methods in Table 1.
VEM possesses features that, on the very least, put it on equal grounds to standard FEM and its variations for any problem, while for some problems it is clearly superior.The latter group includes problems with complex geometries for which a good quality mesh is difficult to obtain, solutions that require very local refinements, and situations for which two independently meshed domains need to be made compatible.These problems highlight VEM's robustness to mesh distortion.Further insight into this topic can be found in [18,21], but in general most works dealing with numerical results for VEM showcase the robustness of the method.

Mesh quality
Accuracy is highly dependent on mesh quality b It has been shown experimentally to be very robust to mesh distortion including small edges, small angles, collapsing nodes, etc.

Computational cost
Standard computation and assembly Requires some overhead computations for the projector and the stabilization term • Assembly of global matrices and vectors is standard a This holds for the choice of DOFs presented here because of the Kronecker-δ property.For a different choice of edge DOFs, a local problem for explicitly reconstructing the same function on the boundary is required.
b See [34] for a detailed survey on what constitutes a good finite element, including assessments of mesh quality.
The features of this space, as well as the corresponding DOFs and their total count were given in Section 2.3.It has been shown that these DOFs univocally define functions on V h (E) [6].As mentioned in the main text, a basis for this space will be denoted by {ϕ i } i=1,...,2n d ⊂ V h (E), where each basis function takes the value 1 at its corresponding DOF and 0 on all others.The global VEM space is defined analogously to standard FEM by Since shape functions shall be projected onto the space P k (E), we now introduce the basis P k = { p α } α=1,...,n k ⊂ P k (E), where n k = (k + 1)(k + 2) is the cardinality of P k (E) since they are polynomials in two variables for both spatial directions.
Denoting with x = ( x1 , x2 ) and h E the centroid and diameter of element E, respectively, let the scaled monomials ξ and η be defined as The cardinality of the basis used depends on the interpolation order k, such that all the monomials of Pascal's triangle of order less then or equal to k are included.For instance, • k = 0: ) ; (A.4) • k = 1: ) , ) ; (A.5) • k = 2: ) , ) , • Arbitrary k: The first three monomials (α = 1, 2, 3) in (A.5), for which the infinitesimal strain ϵ( p α ) = 0 , represent infinitesimal rigid body motions.

A.2. Projection operator Π
The definition of the local projection operator Π ≡ Π E,k : V h (E) → P k (E), stated in (4), can be rewritten as: This condition establishes that the internal elastic energy of the VEM function is computed exactly when tested against strains derived from polynomial displacement fields, such as those used in patch tests.In its matrix form, the projector Π can be thought of as a matrix that maps the coefficients of a decomposition of a shape function ϕ i in the base of V h (E) into the coefficients of a linear combination of polygons in the space P k (E).To compute Π we first have to express v h and p with respect to the bases {ϕ i } and P k , respectively.By expressing Π ϕ i as a linear combination of p β we recall (A.8) and use the linearity of a h E to obtain Taking α = 1, . . ., n k results in a linear system of n k equations with the n k unknowns s i,β that can be written in compact matrix form as Repeating for all basis functions ϕ i , we obtain the matrix B of vectors b i , and the matrix representation Π of the projection operator: The system (A.11) then becomes The definitions and computations of B and G are further explained in Appendices B.1 and B.3, respectively.Looking at the system above we can already point out that, for the monomials that represent rigid body motion (the first three), Eq. (A.10) becomes 0 = 0.These monomials constitute the kernel of the projector Π , so that for them we define a new projection operator, which can be thought of as the following Euclidean scalar product on the 2n v vertex DOFs: More intuitively, the conditions arising from (A.16) are nothing more than equality of average displacement in both spatial directions ( p 1 , p 2 ) and equality of (infinitesimal) rotation ( p 3 ) between v h and Π v h .Finally, decomposing a VEM function v h into a combination of basis functions ϕ i and using Eq.(A.9), we can solve system (A.15) for any VEM function, where the first three rows of B are taken from (A.16).

B.1. Matrix B
The computation of matrix B n k ×2n d in (A.13) is probably the most conceptually important one, since it involves integrating VEM shape functions inside an element, which is a priori not possible.This matrix is defined as • Kernel part: The first three rows are computed with (A.16) such that: Example: E is a triangle with vertices x 1 , x 2 , x 3 and order of accuracy k = 1, the kernel part of B is written as: When the interpolation order increases, the additional columns related to edge and internal DOFs are set to zero; this is because Eq. (A.16) uses only the values for the vertex DOFs.
• Non-kernel part: We can use integration by parts to transform the original computation into two separate integrals: where ne is the outward normal vector to the boundary.The first term on the right computes the product of a VEM function against a polynomial of degree k − 2 in the element's interior and the second one is taken on its boundary, where the VEM functions are explicitly known.Thus, both terms can be obtained only from the knowledge of the VEM function on the DOFs defined in Section 2.3.We examine them separately: (i) The first term is related to internal DOFs.Expressing the strain ϵ( p α ) in Voigt notation, and C as the usual 3 × 3 constitutive matrix in either plane stress/strain, the term ∇ • (Cϵ( p α )) can be written as a linear combination of the monomial basis p β of the space P k−2 (E) (since derivatives are taken twice): Substituting in (B.4), the first term becomes: Denoting by |E| the area of the element, and recalling the definition of the internal DOFs (8) as the moments up to order k − 2, the equation above can be expressed as: due to the property of the shape functions that dof (2kn v +β) (ϕ i ) = 1 when i = 2kn v + β and 0 otherwise.We are thus able to compute this integral without any kind of numerical integration by just using the definition of internal DOFs.

Example:
Consider an order of accuracy k = 2.The basis of the space )} , which means that each element E has 2 internal DOFs: ) , strain and divergence of the stress are, respectively, where d 7,1 is the first coefficient of p 7 , i.e., the constant part (B.5).This leads to: The second term of (B.4) is related to boundary DOFs.We can split the boundary integral as a sum of edge contributions, which are computed using a quadrature rule: where e j refers to the jth edge of ∂ E and ne j its corresponding outward unit normal vector.By choosing Gauss-Lobatto points ξ r (with associated weights w r ) for the definition of the edge DOFs, numerical integration over the boundary is simplified as the values of v h at these integration points are precisely the DOFs of v h .The same is true for ϕ i in particular.We point out that the choice of taking k + 1 Gauss-Lobatto points on each edge (including the vertices) guarantees an optimal choice; this is because the integrand of (B.11) is a polynomial of order 2k − 1, which is the order of integration of the k + 1 Gauss-Lobatto rule.For basis functions associated with internal DOFs, this term is of course zero.
Example: Let us consider the generic VEM shape function ϕ 2i = ( 0 ϕ i ) ⊺ such that its value at the node ξ r is ϕ since ξ r is the only integration point where ϕ 2i is nonzero.

B.2. Matrix D
This 2n d × n k matrix is used to express the projection of a VEM function as a linear combination of the VEM functions themselves, unlike before where the projection was expressed as a linear combination of polynomials.Therefore, through this matrix the projection can be applied from the VEM space onto itself.It is obtained by evaluating the polynomials of degree ≤ k at the DOFs defined for V h (E):

B.3. Matrix G
This n k × n k matrix, appearing in (A.15), is defined as: Its computation is straightforward as it consists on an integration of a product of known polynomials over a polygonal domain.However, it can also be expressed in terms of the already presented matrices B and D as G = B D. Therefore, it is not necessary to compute G if B and D were already obtained-although it is a useful debugging tool to verify that both procedures lead to the same result within numerical precision.

Appendix C. Load term
In this section we show how to compute the VEM loading term r f E due to body forces f = ( f 1 , f 2 ) as detailed in Section 2.3.As for the construction of the stiffness matrix, a projection operator is used to express the VEM displacement v h as a polynomial.We introduce the classical orthogonal L 2 projection Π k−2 : V h (E) → P k−2 (E) such that the loading term can be approximated.Π k−2 is defined with the following orthogonality condition: The computation of the local load vector due to body forces can be discriminated among different cases, depending on the interpolation order k, as follows: (i) Case k = 1: For this particular case the space P k−2 (E) is not defined, so we compute the loading term with an integration rule based on the value of v h at the vertices of the polygon: where x j is the spatial location of the jth vertex.Considering now the basis functions ϕ i to obtain the discrete load vector associated with the ith DOF, Eq. (C.2) becomes (ii) Case k = 2: In this case, a better choice can be made in order to retain optimal convergence rates [35].We use the L 2 projection Π 0 k onto the space of polynomials P k (E), defined as: where which is expressed in matrix form as Extending column-wise to all shape functions, the system can be rewritten as: While the matrix H 0 n k ×n k can be computed numerically since the integrand is a product of known polynomials, the matrix Q 0 n k ×2n d is computed as follows: where Q n k−2 ×2n d is computed as in (C.19).Once Π 0 k is obtained through (C.8), it can be substituted in (25) using Π 0 k instead of Π k−2 to construct the local loading term associated to shape function ϕ i as Extending column-wise to all shape functions, the system can be written as where Π k−2 is the matrix representation of the projection operator:  (24) to construct the local loading term associated to shape function ϕ i as: (C.17)For all cases presented, the complete local load vector due to body forces is obtained by collecting all base functions' contributions: So the generic terms of Q can be written as ∫ E p α • ϕ i dE = |E| dof (2kn v +α) (ϕ i ), which due to the Kronecker-δ property of the basis functions, is equal to |E| if 2kn v + α = i and 0 otherwise.

Fig. 1 .
Fig. 1.(a) Schematic of a solid domain Ω with boundary ∂Ω .Essential and natural boundary conditions, ū and t respectively, are prescribed on disjoint portions of the boundary ∂Ω u and ∂Ω t .The body force is denoted by f .(b) A polygonal discretization Ω h showing in the inset an element E with the corresponding degrees of freedom for order k = 2.

Fig. 2 .
Fig. 2. Visual representation of VEM functions obtained by solving a boundary value problem in a hexagonal domain.Referring back to Fig. 1b, the functions correspond to vertex (a), edge (b), and internal (c) DOFs (either component).

Fig. 4 .
Fig. 4. Stress component in the e 1 direction shown over the deformed configuration for VEM (a) and VEM-FEM (b) discretizations, whose original configuration is portrayed with dashed lines.

Fig. 5 .
Fig. 5. Deformed configuration for the high-order patch test, showing in colors the magnitude of the displacement field.

Fig. 6 .
Fig. 6.Non-conforming patch test.(a) Schematic showing a line Γ that is considered as a non-conforming interface in the discretization.von Mises Stress field σ vM (in MPa) plotted on the deformed configuration.

Fig. 8 .
Fig. 8. L 2 norm (a,b) and energy norm (c,d) of the relative error as a function of total DOFs N dof (a,c) and mesh size h (b,d).

Fig. 9 .
Fig. 9. p-convergence results for the error in L 2 and energy norms, (a) and (b), respectively, as a function of the total number of DOFs N dof .The natural logarithm of the corresponding error is applied again in Figures (c) and (d) to display the exponent θ in Eq. (31).

Fig. 10 .
Fig. 10.(a) Schematic for the infinite plate with a circular hole showing in the inset the computational domain chosen due to symmetry; (b,c) Typical VEM and FEM discretizations.

Fig. 11 .
Fig. 11.Stress convergence curves for the infinite plate with a hole under uniaxial tension and VEM discretizations with k = 1.The error in stresses is measured in L 2 norm by Eq. (34).

Fig. 12 .
Fig. 12. Stress plot for the plate with a hole under unidirectional tension σ 0 at x 1 = ∞.In this example VEM polygons are placed to resolve the non-conforming interface between the triangular and quadrangular meshes.The figure shows the normalized element stress (average value per element) in the e 1 direction, directly calculated from the finite element solution.The maximum value 2.98 approximates the well-known value of 3.

Fig. 13 .
Fig. 13.(a) Schematic for the Flamant problem, a half-space loaded with a point force as shown; (b) VEM discretization.

Fig. 14 .
Fig. 14.Comparison between the exact stresses and those resulting from the VEM approximation.The stress components σ 22 (top) and σ 12 (bottom) were measured along a horizontal line located at x 2 = 0.1 m below the applied load.

This n k− 2 ×
2n d matrix is obtained as

Table 1
Comparison between standard FEM and VEM.
(8) 1 ) dof 2n d ( p 2 ) ...dof 2n d ( p n k ).≤ i ≤ 2 k n v , so dof i ( p α ) is the value of p α in the node associated with the ith DOF.(ii) Internal DOFs: 2kn v + 1 ≤ i ≤ 2n d .Applying definition(8)to p α and performing numerical integration on the product of known polynomials, we are able to compute those terms as: dof (2kn v +β) p 1 ) a E ( p 1 , p 2 ) • • • a E ( p 1 , p n k ) a E ( p 2 , p 1 ) a E ( p 2 , p 2 ) • • • a E ( p 2 , p n k ) E ( p n k , p 1 ) a E ( p n k p 2 ) • • • a E ( p n k , p n k ) a r 1 r 2 . . .r 2n d ] .(C.16) Matrix Q n k−2 ×2n d is computable thanks to the internal DOFs, as shown in Appendix C.1, while matrix H n k−2 ×n k−2 can be computed numerically since the integrand is a product of known polynomials.Once Π k−2 is obtained through (C.15), it can be substituted in ϕ 2 dE . . .