A virtual element formulation for general element shapes

The virtual element method is a lively field of research, in which considerable progress has been made during the last decade and applied to many problems in physics and engineering. The method allows ansatz function of arbitrary polynomial degree. However, one of the prerequisite of the formulation is that the element edges have to be straight. In the literature there are several new formulations that introduce curved element edges. These virtual elements allow for specific geometrical forms of the course of the curve at the edges. In this contribution a new methodology is proposed that allows to use general mappings for virtual elements which can model arbitrary geometries.


Introduction
The many different approaches to the approximate solution of problems involving partial differential equations include finite difference schemes, finite elements, finite volume techniques, boundary elements, and particle methods. Within the finite element method there have been various significant developments, including for example classical isoparametric mapping, see e.g. Hughes [15] and Zienkiewicz and Taylor [35], but also isogeometric schemes, see Cottrell et al. [10].
Research continues to be motivated by the goal of developing stable, efficient and robust discretization schemes for finite deformation applications in solid mechanics. Within this line also the virtual element technology is further refined and applied to nonlinear problems in mechanics, see e.g. Beirão da Veiga et al. [6], Chi et al. [9], Wriggers et al. [34], Artioli et al. [3], De Bellis et al. [11], Wriggers and Hudobivnik [33], Aldakheel et al. [1], Hussein et al. [17] and De Bellis et al. [12]. So far the assumptions-even for higher order virtual elements-contain the restriction to straight edges of the elements that are directly defined in the physical space, see e.g. Beirão da Veiga et al. [5]. This makes the definition of virtual elements having a general geometric shape more complicated which is due to the fact that mappings like the isoparametric map for finite elements or NURBS type B P. Wriggers wriggers@ikm.uni-hannover.de 1 Institute for Continuum Mechanics, Leibniz Universität, Hannover, Germany maps in isogeometric analysis are not easily applicable. New formulations that introduce curved element edges allow specific geometrical forms of the course of the curve at the edges were discussed in Beirão da Veiga et al. [7], Artioli et al. [4] and Aldakheel et al. [2]. This paper introduces a possibility to employ general mappings also within the virtual element formulation. The idea is to map a virtual element, defined at a reference configuration, to a general shape in the initial (physical) configuration. With such a map general shapes of virtual elements can be created in the initial configuration while preserving the straight edges in the reference configuration. This type of mapping was investigated for shell problems in e.g. Pimenta and Campello [27], but so far has not been used in two-and three-dimensional solid mechanics and within the virtual element methodology.
The method is developed for hyperelastic materials undergoing finite strains. However it can as well be applied for small strain cases as for material nonlinearities such as plasticity. Here we use the neo-Hookean strain energy as a model. The mapping from the reference to the initial configuration is performed either with an isoparametric map based on a quadratic ansatz for the geometry or a Bezier type map using NURBS: Both formulations perform extremely well in a series of benchmark tests involving regular and Voronoi meshes.
After presenting the governing equations for nonlinear elasticity in Sect. 2. Section 3 is devoted to the general mapping procedure, followed by the formulation of a low order

Governing equations for finite elasticity
Consider an elastic body that occupies the bounded domain 0 ⊂ R 2 . The body 0 has a boundary which comprises non-overlapping sections D and N such that D ∪ N = (Fig. 1).
The position x of a material point initially at X is given by the motion where u is the displacement field that depends generally on the time t. We also define the deformation gradient F by the gradient being evaluated with respect to X. The body satisfies for the static case, on 0 the equation of equilibrium with the body forcef and the first Piola-Kirchhoff stress P. The Dirichlet and Neumann boundary conditions are respectively with N the outward unit normal vector,ū the prescribed displacement, andt the prescribed surface traction on N .
By introducing a strain energy function Ψ (u) for elastic problems the first Piola-Kirchhoff stresses follow from For a homogeneous compressible isotropic hyperelastic material we adopt the neo-Hookean strain energy function for the two-dimensional case in which λ and μ are the Lamé constants. This strain energy is known as Neo-Hookean model. The right Cauchy-Green tensor C(u) is defined as C = F T F and the Jacobian J (u) of the deformation is given as J = det F > 0.
In case of a hyperelastic material the development of a virtual element can start from the potential energy function directly instead of using the weak form, see e.g. Wriggers et al. [34]. In that case the potential energy can be written with (7) as

General mapping
The idea for the construction of general shaped virtual elements is to use a function that maps the coordinates X r of a virtual element from a reference configuration r to the coordinates X 0 of the initial configuration 0 , see Fig. 2. The deformation is then given by the vector x which describes the current configuration d that a solid assumes under loading. The mapping from the reference to the initial configuration can be defined in an arbitrary way using either an isoparametric map or e.g. NURBS functions. For isoparametric maps we can use shape function N I that describe a polynomial basis in the reference configuration. The map is given by In this paper we use a quadratic function for the isoparametric map, explicit expressions for the ansatz functions can be found in e.g. Hughes [15], Wriggers [32] and Onate [25]. The  1 A more general mapping that includes also mapping onto domains that cannot be exactly defined by polynomials is provided by NURBS, see e.g. Farin [13] and Piegl and Tiller [26]. These Non-Uniform Rational B-Splines are obtained for a two-dimensional map by a projection of a B-spline quantity from three-dimensional space which allows exact constructions of geometrical objects such as circles and ellipses. The NURBS projection maps from a parameter space r with 0 ≤ η 1 ≤ 1 and 0 ≤ η 2 ≤ 1 to the initial configuration 0 follows as where P 0 I J are the control points that define the domain in the initial configuration. The functions R I J (η 1 , η 2 ) are the rational B-spline basis functions. They are explicitly given by with the weights w I J . For more details related to the order of the B-spline basic functions we refer to Piegl  [26], Hughes et al. [16] and Cottrell et al. [10]. The latter used this parametrization to develop the discretization technique known as isogemetrical analysis. We note that the NURBS functions can also be linearly mapped to . Thus we will write the mapping for both, isoparametric and NURBS, mapping as being a function of the coordinates X r = (X r , Y r ). According to Fig. 2 the deformation gradient in (2) can be computed via with following from the mapping functions (9) and (10). The Jacobian of this mapping is given by J 0 = det F 0 . Furthermore the volume element in the initial configuration is connected to a volume element in the reference configuration by d 0 = J 0 d r . Depending on the used mapping we either have to differentiate in (13) with respect to (X r , Y r ) or (η 1 , η 2 ). The deformation gradient F u that describes the deformation from the initial to the current configuration yields Based on these general kinematical quantities we can compute the right Cauchy Green tensor in terms of the deformation gradient F u as Now the potential (8) can be reformulated using (15) where in the general three-dimensional case is the mapping of the area elements of the reference to the initial configuration, either with respect to the coordinates of the isoparametric or the NURBS mapping, respectively. Furthermore, the Jacobian J describing the volume change in the strain energy (7) is given by J u = J J −1 0 . In the two-dimensional case one can compute the gradient of the mapping from the reference to the initial configuration directly from the isoparametric or the NURBS map as The differentiations in (18) have to be performed with respect to the reference coordinates. Note that the computation of the derivatives of the rational B-splines is more involved due to the definition of the B-spline in (11). Note that description of a solid may require more than one reference configurations. In this case it is possible to combine several configurations to model the solid which will be shown in the numerical example section.

Formulation of the virtual element method
The main idea of the virtual element method (VEM) is a projection of the deformation onto a specific ansatz space. In classical analysis using virtual elements a domain is partitioned into non-overlapping polygonal elements which need not be convex. Here we extend this possibility by mapping the polygons into other geometrical objects using a nonlinear map as discussed in the last section. In this contribution we restrict ourselves to two dimensional solids and to a low-order approach using linear ansatz functions. Due to the construction of the mappings the developed methodology works without any changes for higher order ansatz spaces of virtual elements.
The construction of a low order virtual element is based on a linear ansatz space u h , see e.g. Beirão da Veiga et al. [5], that has unknown displacements u k at the vertices k of the polygon, a linear ansatz for the displacement field u h at the edges of the polygon and the property that Div(∇u h ) = 0. The polygonal elements can have arbitrary number of vertices where n V is the total number. Furthermore, the element shape is restricted to straight edges, see Fig. 3, in the reference configuration e r .

Ansatz functions
Generally the virtual element method relies on the split of the ansatz space for the displacement field u h into a part u π and a remainder Classically the projection u π is defined at element level by a function which is directly formulated in the coordinates that describe the real geometry of the domain. Hence the ansatz for a linear polynomial function yields for the two dimensional case in the initial configuration

Computation of the projection
There exist several methods to obtain the parameters a i in terms of the unknown nodal values of a virtual element e as depicted in Fig. 3. Here we use a Galerkin projection of the gradients to project u π onto the ansatz u h . Since this formulation does not rely on the mechanical weak form it leads to the same projection for linear and nonlinear elasticity problems. The Galerkin projection of the gradient ∇u π is carried out in the initial configuration e 0 of the virtual element. It can be formulated as Here the index "0" in ∇ refers to the differentiation with respect to the initial coordinates X 0 . N π is the vector containing a weighting polynomial that has the same order as u π . For the chosen linear ansatz the gradient ∇ 0 N π is constant when we select N π = {1, X 0 , Y 0 }. Thus we can write (21) as Since ∇ 0 u π is constant in e 0 for the linear ansatz in (20) we can write for the left side of the integral above as Using now the Gauss theorem the second integral in (22) can be transformed to a boundary integral with N 0 being the outward normal to the element area 0 in the initial configuration.
The integral on the right hand side in (24) can be transformed to the reference domain directly. The question is how to compute the normal vector N 0 at an edge γ ∈ 0 where γ stands for an edge between two vertices and 0 is the entire boundary of the virtual element. Two different possibilities are available for the transformation: • We can introduce along an edge γ the coordinate ξ , see Fig. 4. The normal vector N 0 can then be discretized directly using the map (think of convective coordinates) Where the edge γ 0 is defined through the reference configuration by which yields the coordinate at the edge X 0 (X(ξ )) as a function of ξ . Now we define a linear ansatz for u h at the edge γ 0 as With d γ 0 = X γ 0,ξ dξ we can write the integral as a sum over all n γ edges γ ∈ 0 of the virtual element Note that the coordinate ξ is related to a line in r that is straight, see right side in Fig. 4.
For a specific map X 0 the integration of the right hand side in (28) follows by using an integration rule, e.g. Gauss integration, where w g are the weights and ξ g the integration points with respect to the reference configuration. This integration holds for isoparametric and NURBS mappings. The integration order in (29) depends on the mapping function, e.g for a quadratic isoparametric map the integrand in (29) is represented by a second order polynomial, thus a two point Gauss integration is sufficient. But for a Bezier type mapping the number of Gauss points has to be increased. • Another way to integrated the right hand side in (24) is to use Nanson's formula. Then the integral can be transferred to the reference configuration by where now N r is the normal related to the straight edges of the virtual element in the reference configuration. The inverse of this gradient is given for the isoparametric map in the two-dimensional case simply by One can combine (23) and (33) to finally obtain the projection of the gradient or use (23) and (28) which yields the alternative form In (34) and (35) the area 0 in the initial configuration is computed using the reference configuration of the virtual element. The integral can be evaluated over the edges using partial integration and Gauss integration All quantities are related to the reference configuration, see Fig. 4. These are the coordinates X r , the normal N γ r at edge γ , the Gauss points ξ g and associated weights w g and the length of an edge l γ r . Note that also the first formulation leading to (35) can be applied to compute 0 .
Finally the projection in (34) can be evaluated by inserting ansatz (27)

into this equation. Again Gauss integration is employed
which yields the parameters a 3 to a 6 in (20) that now depend on the displacements u k at the vertices k of a virtual element. Again (35) can also be used to compute the projected gradient.
In the projection in (38) only the gradient plays a role. This allows to compute the parameters a 3 to a 6 of (20) directly. Hence the constant parts of the ansatz disappear. Thus the parameters a 1 and a 2 , related to the constant parts, have to be obtained in a different way. Here the average displacement of the virtual element in the initial configuration along the edge 0 is used. The equivalence includes the parameters a 1 and a 2 and from (39) all parameters in (20) can be computed. A simplified expression follows, if the integrals are evaluated by the Gauss-Lobatto rule. In that case the nodal values of u π and u h appear in a sum over all vertices. This leads to Ansatz (20) can now be inserted into (40) and evaluated at X 0 k = n iso I =1 N I (X rk , Y rk )X 0 I for the left hand side, while on the right hand side only the nodal values of u h have to be used. This finally leads to two equations for the parameters a 1 and a 2 .

Construction of the virtual element
In this contribution a linear ansatz is employed for the discretization using virtual elements. Thus the gradient of the displacement field is approximated by a constant part as a result of the projection in the last section. A construction of a virtual element which is based only on this projection would lead to a rank deficient element once the number of vertices is greater than 3. Thus the formulation has to be stabilized, as in the case of the classical one-point integrated elements developed by Flanagan and Belytschko [14], Belytschko and Bindeman [8], Reese et al. [30], Reese and Wriggers [29], Reese [28], Mueller-Hoeppe et al. [24], Korelc et al. [21] and Krysl [22], to mention some key contributions.
The formulation is provided here for finite deformation formulations. The potential U is a split into a constant part of the deformation gradient and an associated stabilization term. In this contribution we employ the hyperelastic potential function (16) as basis for the virtual element. By summing up all element contributions for the n e virtual elements that discretize a domain 0 one obtains In the following we will first discuss the formulation of the consistency part U e c (u π ) that stems from the projection, see last section. Furthermore, a possibility for the stabilization U e stab of the virtual element method will be considered.

Constant part due to projection
The consistency part in the potential (41) can be evaluated by inserting the results obtained in the last section. This yields with respect to the initial configuration for a virtual element Due to the fact that we are able to compute the projection of the gradient ∇ 0 u π in the initial configuration directly in (38) we obtain the deformation gradient describing the map from the initial to the current configuration as which is constant. Note that here also (35) could have been inserted for ∇ 0 u π which leads to exactly the same results. Now the right Cauchy-Green tensor is given by which is also constant within the virtual element e but depends in a nonlinear fashion on the nodal displacements. The strain energy of the hyperelastic material can now be computed in terms of the projected deformation measures which are the right Cauchy-Green tensor C u and the Jacobian J u . The latter is given by Hence all quantities in (45) depend only the projection ∇ 0 u π . The first integral in (42) can easily be computed since Ψ (C u ) is constant in (X 0 , Y 0 ) and e 0 known from (37). The two loading terms have to be transformed for evaluation to the reference configuration, see also (28), where γ denotes the edges of the virtual element that are loaded by surface tractiont γ . Once the integration over the reference coordinate is finalized, the derivations of the potential (47) with respect to the n V unknown nodal displacements at the vertices of the virtual element u e = {u 1 , u 2 , . . . , u n V } can be carried out. This approach exploits the fact that ∇ 0 u π depends directly  (20) and (40) to the nodal displacements as well.
All derivations, leading to the residual vector R c e and the tangent matrix K c T e , were performed with the symbolic tool AceGen, see Korelc and Wriggers [20]. This yields for (47) where u e are the nodal displacements of the virtual element e at the vertices.

Stabilization techniques for nonlinear virtual elements
In the literature on virtual element technologies basically two different stabilization techniques were discussed that work well for classical solid mechanics problems. The first stabilization is directly based on the degrees of freedom. It introduces a point wise error measure between the nodal values u k and the approximation function u π evaluated at the vertices X k , see e.g. Beirão da Veiga et al. [5], Beirão da Veiga et al. [6] and Chi et al. [9].
Here we use the stabilization approach developed in Wriggers et al. [34] for virtual elements. The essence of the approach is to introduce a new, positive definite strain energŷ Ψ and to define the stabilization contribution to the strain energy by Such a stabilization was also used in Krysl [22] for stabilized mean strain formulations of finite elements. The second term on the right side ensures the consistency of the total potential energy, which is now given by ActuallyÛ can be selected differently from the strain energy that describes the constitutive behavior of the solid, see e.g. Wriggers et al. [34], since it is only for stabilization. Note thatÛ (u h ) −Û (u π ) disappears for element sizes e 0 → 0. In this paper we use the same energy Ψ as in the consistency part U c , however with different Lamé parametersλ andμ. The stabilization energy follows then by assembly over all virtual elementŝ The computation of the consistency part U c (u π ) was already discussed above. The last part of the potential energŷ U (u π ) depends also on u π , thus it can be evaluated as the consistency part, see Sec. 4.3.1.
The question is now how to compute the stabilization term U (u h ) since u h is not known within the virtual element. The idea is to use an internal triangular mesh. This inscribed triangular finite element mesh, see Fig. 5, consists of n int linear three-noded triangles that are connected to the nodes of the virtual element. Hence this internal mesh does not creat new nodal unknowns. The ansatz functions for the triangles are then assumed to approximate u h in (19).
For this discretization of the potentialÛ e (u h ) within a virtual element e 0 one has to add all contributions of the internal triangles i which can be directly formulated with respect to the initial configuration 0m .
To fulfill the condition thatÛ (u h )−Û (u π ) → 0 for small virtual elements e 0 → 0 one has to compute the poten-tialŪ (u π ) using the area related to the triangular elements e OT = n int m=1 i 0m which leads tō Due to the straight edges of the inscribed triangular mesh, see Fig. 5, this slightly differs from the way the element area is computed in the consistency part in (47). However such approximation can be admitted since (53) is only the stabilization energy which approaches zero for fine meshes, see above.
All further derivations leading to the residual vector R s e and the tangent matrix K s T e were performed in the same manner as in (48) with the symbolic tool AceGen, see Korelc and Wriggers [20]. This yields for (49) Thus the final residual and tangent matrix of the virtual element are given by the sum of expressions (48) and (54): R e = R c e + R s e and K T e = K c T e + K s T e . The values of the Lamé parameters in the strain energy (52) have to be selected in a proper way. Krysl [23] suggested for cuboid finite elements a procedure that is based on a comparison of the bending energy of a thick beam and the bending energy of the finite element. This yields Lamé parametersλ andμ that enhance the bending behavior of the element. Since this procedure is not directly applicable to arbitrary virtual elements we suggest a simple and computationally efficient way to compute the parameters from the basic geometric data of the virtual element in the initial configuration e 0 , see Fig. 6. The procedure was described in detail in Wriggers et al. [34] and will not be reported here. We just state the final results. The Young's modulus is corrected bŷ while the Poisson ratioν is kept constant asν = 0.3 since the Poisson ratio does not influence the convergence behavior of the element and avoids locking in the stabilization term for incompressible problems. The factor β depends on the height to length ratio of the brick element. For a virtual element with arbitrary shape this ratio can be approximated by using the inner and outer radii, R 2 i , R 2 a respectively, see Fig. 6, to obtain

Fig. 6 Inner and outer radius of a virtual element
The inner radius is computed by using the distance from the geometrical centre to the convex hull of the virtual element while the outer radius is defined by the maximum distance of nodes related to the virtual element. Another way to obtain β can be found in van Huyssteen and Reddy [31] which is based on employing ellipses instead of the circles in Fig. 6.
Once the corrected Young's modulus and Poisson ratio are determined the Lamé parameters follow from Remark From (55) the Lamé parameters can be written aŝ λ = φ λ andμ = φ μ, and with that also the strain energy aŝ Ψ = φ Ψ . This means that we can re-write the strain energy parts in (41) for a virtual element e 0 with (47) and (53) as where U e (u h ) is the stabilization energy from (52) computed for the Lamé parameters λ and μ. Observe that for φ = 0 only the energy related to the consistency part acts which would lead to a singular element. On the other hand for φ = e 0 / e 0T the consistency part disappears and a pure finite element formulation based on the internal triangles is left over. In Wriggers et al. [34] it was shown that the latter case does not have the same good solution behavior as a virtual element with a factor φ computed from (55).

Numerical examples
In this section we compare the new mapping procedure for virtual element discretizations with existing formulations.
The main goal is to show that this formulation blends perfectly into the idea of virtual elements and produces the same results for higher order isoparametric or NURBS mappings. The examples are subjected to loads that lead to finite deformation strain states. For an efficient solution all numerical simulations of the nonlinear problems have to be performed by using a Newton-Raphson algorithm with load stepping when necessary. Due to the fact that all formulations are linearized in a consistent manner, using AceGen, see Korelc [18], quadratic convergence is achieved. All geometric and constitutive data in the examples are given in consistent units.
In all examples either NURBS (denoted by a MB) or isoparametric functions (denoted by a MN) were used for the map from the reference to the initial configuration. Simulations were conducted with different virtual elements using the software AceFEM, see Korelc [19] and Korelc and Wriggers [20]: The results of these discretizations were compared with linear Q1 and quadratic Q2 quadrilateral finite elements.

Squeezing of an ellipse
The solid structure shown in Fig. 7 is subjected to a traction loadt y = −300 at the top. The ellipse is defined by its principal axes R 2)R y . In Fig. 7b, c the structured and Voronoi meshes are shown, respectively, with respect to the initial configuration 0 . The associated meshes in the reference configuration r are displayed in Fig. 7d, e.
The von Mises stress distribution is depicted in Fig. 8 for a relatively coarse mesh of 8 × 12 elements. The discretizations in Fig. 8a, d show results that were computed with a standard finite element and virtual element discretization without mapping. The results depicted in Fig. 8b, c, e, f were obtained with the new mapping techniques from reference to initial configuration using either the isoparametric (M0) map or the NURBS (MB) approach. The results are almost equal for the standard virtual element formulation with 8 node elements (d) and the approach using M0 or MB mapping in (b) and (e) which shows that the mapping formulation works correctly. Here we do not expect much improvement of the results since there is for the linear ansatz only a small difference between the curved geometry and the straight edge approach of the classical virtual element scheme.

Square plate with a hole
A plate with a hole, is depicted in Fig. 9, is subjected to an extension. The plate in Fig. 9a has a height of H = 2, a width of L = 2. The radius of the internal hole is R = 0.5. The displacement is fixed at the bottom in y-direction. In the mid of the plate, at x = L 2 the displacement in x-direction was fixed. This example was considered in De Bellis et al. [12].
A prescribed displacementū y = 1.5 H is applied at the top of the plate. The Lamé parameters, that govern the hyperelastic response of the solid, are given as λ = 71.0719 and μ = 36.6128. In Fig. 9b, c the structured and Voronoi meshes are shown, respectively, with respect to the initial configuration. In Fig. 9d, e the associated meshes in the reference configuration are depicted. As well a quadratic isoparametric map as a NURBS map were employed to map the structure from the reference to the initial configuration. Note that here a combination of 8 mappings was used to model the problem. The different mappings were tied together in a suitable way such that the displacement constraints related to the boundary value problem were enforced correctly.
In Fig. 10 the von Mises stress is depicted for different discretizations and mappings at the deformed configuration. It is clearly visible that the new scheme leads to results that cannot be distinguished from existing virtual element formulations, see Wriggers et al. [34], and a finite element solution using quadrilateral elements with quadratic shape functions.
In order to compare the different approaches a convergence study was performed for the displacement component u x at the point (x = L, y = H 2 ) and the total force R y , needed to obtain the resulting deformation state. Both results are depicted in Fig. 11 where a range from 30 to 9600 elements was used. Since there is no analytical solution available for this finite strain example an overkill solution based on the quadratic finite element with a mesh of around 0.5 Mio. element is used in this convergence study as reference solution for the displacement component u x,re f and the total reaction R y,re f .
As expected the solution with the quadratic finite element Q2 is superior. However it seems that the virtual element based on the NURBS mapping performs slightly better than the other formulations, especially for the Voronoi discretization. All virtual element formulations yield a better response than the Q1 finite element.

Bending of a C-profile
In the last example a C-profil is considered, see Fig. 12. On the right side the web of the profile is cut. The geometrical data are H = 3, L = 1 and R = 0.5. The profile is subjected at its flanges to a traction loadt x = 2 in opposite directions, see Fig. 12a which varies from 0 to 2 in magnitude. The Lamé constitutive parameters of the hyperelastic strain energy function are λ = 71.0719 and μ = 36.6128. The profile is fixed in the middle of the left web at y = H 2 . As in the previous example, the total mesh of the boundary value problem is obtained by using 10 mappings from different reference configurations, shown in Fig. 12d, e. These generate the initial configuration 0 in Fig. 12b, c for the isoparametric and NURBS map.
Due to the loading the C-profil bends and at the largest load level it assumes the deformed state that is depicted in Fig. 13. Again the results are plotted for different discretizations. The von Mises stresses are displayed in Fig. 13a for a finite element solution using Q2 finite element quadrilaterals and in (b) for a regular mesh with virtual elements that have 8 nodes. Furthermore, Fig. 13 includes results for the new mapping procedure. In (c) and (e) are the results for the isoparametric mapping for regular and Voronoi meshes, respectively, while (d) and (f) relate to the solutions obtained with the NURBS mapping. All discretizations produce basically the same answers. Difference can be seen in the convergence study that is shown in Fig. 14. A solution based on quadratic finite elements with a very fine discretization of about 0.65 Mio. elements is used as a reference since analytical solutions are not available.
All discretization schemes produce the same load displacement curve, see the left part of Fig. 14. Differences occur in the convergence of the displacement u y that is measured at the point (x = L, y = H − L). All formulations converge.
Again the virtual element formulation using a Voronoi discretizations depicts good coarse mesh accuracy.

Summary and conclusions
In this contribution a virtual element method for finite strain elasticity was derived. The key novel feature is a general mapping which allows to deviate from the assumption of straight edges of virtual elements. It was shown that general mappings like the isoparametric map and the NURBS approach, related to computer aided design, can be introduced to model exact geometries of complex structures with virtual elements. The formulation is simple and has been presented for both isoparametric and NURBS maps. It seems that the projected formulation is better in the case of Voronoi meshes when compared to solutions obtained with standard virtual elements using Voronoi meshes. For the special case of structured meshes the difference is very small. Here the virtual element solutions are basically reproduced in an exact manner. The method proposed here is amenable to extensions of various kinds: for example to higher-order VEM formulations, like serendipity elements, to problems in three dimensions, and to other nonlinear problems such as those involving inelastic material behavior. First results indicate that the same mapping procedure can be employed for virtual elements with higher order ansatz functions. Acknowledgements Open Access funding provided by Projekt DEAL.
The first author gratefully acknowledges support for this research by the "German Research Foundation" (DFG) in the collaborative research center CRC 1153 "Process chain for the production of hybrid highperformance components through tailored forming" while the third author acknowledges support for this research by the Priority Program SPP 2020 under the project WR 19/58-2.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.