Curvilinear Virtual Elements for 2D solid mechanics applications
CCurvilinear Virtual Elements for 2D solidmechanics applications
E. Artioli ∗ , L. Beirão da Veiga † , and F. Dassi ‡ Department of Civil Engeneering and Computer science,University of Rome Tor Vergata, Via del Politecnico 1 - 00133 Roma, Italy Department of Mathematics and Applications,University Milano - Bicocca, Via R. Cozzi 55 - 20125 Milano, ItalyOctober 3, 2019
Abstract
In the present work we generalize the curvilinear Virtual Element technology, intro-duced for a simple linear scalar problem in a previous work, to generic 2D solid mechanicproblems in small deformations. Such generalization also includes the development of anovel Virtual Element space for displacements that contains rigid body motions. Ourapproach can accept a generic black-box (elastic or inelastic) constitutive algorithm and,in addition, can make use of curved edges thus leading to an exact approximation of thegeometry. Rigorous theoretical interpolation properties for the new space on curvilinearelements are derived. We undergo an extensive numerical test campaign, both on elasticand inelastic problems, to assess the behavior of the scheme. The results are very promis-ing and underline the advantages of the curved VEM approach over the standard one inthe presence of curved geometries.
The Virtual Element Method (VEM) is a recent technology introduced in [13, 10] (see also[11]) for the discretization of problems in partial differential equations, generalizing the FiniteElement Method to meshes using general polygonal/polyhedral elements. In the framework ofStructural Mechanics, after the initial papers [26, 30, 15, 43], many applications and develop-ments were pursued by the community. Among the many publications, we here refer to a longbut non exhaustive list [20, 29, 33, 6, 5, 21, 3, 42, 44, 4, 8, 7, 16, 41, 34, 27, 22]. Finally, itmust be noted that other polytopal methods exist in the Structural Mechanics literature; wehere limit ourselves in citing [39, 18, 23, 28] and references therein.Since standard Virtual and Finite Elements make use of meshes with straight edges, thedomain of interest is approximated by a linear interpolation of the boundary, a procedure whichintroduces a geometry approximation error that can dominate the analysis for second-to-higherorder schemes. In the literature there has been introduced many approaches and variantsin order to have a better (or even exact) approximation of the geometry. Among the manypossibilities, we here cite two among the most known, that are isoparametric Finite Elements(see for instance [24]) and Isogeometric Analysis (see for instance [32, 25]). As shown in [12], itturns out that for Virtual Elements in 2D the extension to cases with curved edges (and thus ableto approximate exactly the geometry of interest) is quite natural. By changing the definitionof the edge-wise space, the method is able to accommodate for general boundaries described ∗ [email protected] † [email protected] ( corresponding author ) ‡ [email protected] a r X i v : . [ m a t h . NA ] O c t y a given parametrization (for instance obtained by CAD). The approximation properties ofthe space and the convergence of the method turns out to be optimal, and there is no need tobuild a volumetric parametrization of the object as in standard IGA or to carefully positionthe “mapping points” as in isoparametric FEM. Currently, a limitation of this VEM approachis that extension to 3D is not immediate and needs more investigation. Furthermore, we citealso [2, 17, 19] that represent other examples of approaches that allow for curved polygonalelements, still on a simple linear elliptic model problem as in [12].In the present work we start from the basic construction introduced in [12], that was on asimple (scalar) linear model problem, and generalize it to the more complex situation of smalldeformation problems in structural mechanics. Our approach, in the spirit of [15, 6, 5], canaccept a generic black-box (elastic or inelastic) constitutive algorithm but, in addition, canmake use of curved edges thus leading to an exact approximation of the geometry. When theVirtual Spaces of [12] are applied in the present context of structural mechanics, a potentialdrawback appears. Indeed, such spaces (although holding optimal approximation properties)contain constant functions but do not contain linear ones. As a consequence, the ensuing vectordisplacement space does not contain rigid body rotations. In order to deal with this issue, inthe present paper we propose also an alternative novel VEM space, that has the same degreesof freedom as the original one but the added property of containing all rigid body motions.Such space is based on a non-conventional construction that is partially in physical space andpartially mapped from a parameter space. We furthermore prove the optimal approximationproperties of this new space, the proof being interesting due to the particular nature of theproposed space.After presenting the Virtual space, the associated method and the theoretical results, weundergo an extensive numerical test campaign, both on elastic and inelastic problems, to assessthe behavior of the scheme and make some comparison. The results are very promising andunderline the advantages of the curved VEM approach over the standard one for elements ofsecond (or higher) order. It is notable that, from the pure coding standpoint, the extensionfrom standard to curved elements is essentially only a change in the quadrature rules. Thismakes the proposed approach a viable choice, since modifying a standard VEM scheme into acurved approach implies a modification (and not a revolution) in the existing codes.The paper is organized as follows. After introducing some preliminaries in Sec. 2, we describethe two VEM spaces (original and novel) in Sec. 3. Afterwards, we derive the theoreticalinterpolation properties of the new space in Sec. 4, and describe the proposed method in Sec. 5.Finally, Sec. 6 is devoted to the numerical tests, divided into elastic and inelastic cases. We assume a bounded 2D domain Ω describing the body of interest, such that its boundary ∂ Ω is a finite union of regular components Γ i , i = 1 , , .., m . Each of such regular componentsis parametrized by an invertible mapping γ i from an interval in the real line into Γ i (this isa typical situation, for instance, for domains described by CAD). We assume that we have aconforming mesh of polygons Ω h , that may have curved edges on the boundary in order todescribe exactly the domain Ω (see Fig. 1). It is not restrictive to assume that each curvedface is a subset of only one Γ i and therefore regular. In order to symplify the notation in thefollowing we will simply use Γ and γ : [0 , L ] −→ Γ . to indicate a generic curved part of the boundary, and its associated parametrization.In the following we will denote with E a generic element of the mesh and with e a genericedge. As usual the symbol h will be associated to the diameter of objects, for instance h E will denote the diameter of the element E and h e the (curvilinear) length of the edge e . An h without indexes denotes as usual the maximum mesh element size.We will adopt the standard notation for H s Sobolev norms and semi-norms (with s ∈ R non-negative) on open measurable sets. In the case of a curved edge, the norm is intended withrespect to the arc-length parametrization of the edge. We will denote by (cid:46) , (cid:38) an inequality of2igure 1: Sample meshed domain Ω.two real numbers that holds up to a costant that is independent of the particular mesh Ω h or,in case of a local estimate, of the element E .In this paper we focus on the equilibrium problem for a two-dimensional medium Ω witheither elastic or inelastic constitutive material behavior, in small deformation regime. Morespecifically, in the elastic case the variational formulation of the problem reads Find u ∈ V such that Z Ω σ ( x , ε ( u )) : ε ( v )d x = F ( v ) ∀ v ∈ δ V , (1)where V (and δ V ) represents the space of kinematically admissible displacements (and its vari-ations), ε denotes the usual linear strain operator, σ represents the elastic constitutive lawrelating strains to stresses, F indicates a linear operator representing volume and surface forcesacting on Ω.In the case of an inelastic material behaviour, the constitutive law σ will depend on somehistory variables and the problem will be written in a pseudo-time evolution formulation, seefor instance [35]. Remark . We note that the method here proposed trivially extends also to the case wherethere are (fixed) interfaces also inside Ω, for instance in the presence of material coefficientjumps. In such case also internal elements may have curved edges in order to adapt to thematerial interface. For the sake of simplicity, in the exposition we will refer only to the boundaryof Ω as a source of curved edges.
In the present section we briefly review the space proposed in [12] and introduce also analternative space that has the additional advantage of including rigid body motions. Indeed,the space proposed for the diffusion problem in [12] has good approximation properties and istherefore suitable also as a displacement space for problems in structural mechanics. On theother hand, such space does contain constant polynomials (that is translations of the body)but not (linearized) rigid body rotations. Therefore we here propose an alternative space thatincludes also rigid body rotations without increasing the number of degrees of freedom.
We here quickly review the space of [12], here extended to the vector-valued version. Asusual, we define the space element by element. Let therefore E ∈ Ω h . Note that E may havesome curved edge, laying on the curved boundary Γ. For each curved edge e , let γ e : [ a, b ] → e denote the restriction of the parametrization describing Γ to the edge e . Then, we indicate the3pace of mapped polynomials (living on e ) as e P k ( e ) = n p ◦ γ − e : p ∈ ( P k [ a, b ]) o . The local virtual element space on E is then defined as V h ( E ) = n v ∈ [ H ( E ) ∩ C ( E )] : v | e ∈ [ P k ( e )] if e is straight , v | e ∈ [ e P k ( e )] if e is curved , − ∆ v ∈ [ P k − ( E )] o . (2)The associated degrees of freedom are (see [12] for the simple proof) • D1 : pointwise evaluation at each vertex of E , • D2 : pointwise evaluation at k − E • D3 : moments R E v · p k − for all p k − ∈ [ P k − ( E )] .As usual, the global space is obtained by a standard gluing procedure V h = n v ∈ H (Ω) : v | E ∈ V h ( E ) ∀ E ∈ Ω h o with the obvious extension of the local degrees of freedom to global ones. Note that, clearly,Remark 3.3 below applies also to the space above. Let e be a (possibly curved) edge of the polygon E , with endpoints ν , ν . In the followingwe denote by h e the (curvilinear) length of the edge e . Note that, as a consequence of the factthat the boundary Γ is fixed once and for all (see [12]) it holds h e ’ k ν − ν k .We will now present a preliminary linear space of (two component) vector fields living on e , with associated degrees of freedom given by the vector evaluation at the two endpoints of e .Let therefore u and u in R represent vector values associated to ν and ν , respectively. Webuild an associated vector field on e as follows. Let F : R → R be the unique mapping thatsatisfies ( F ( x ) = a + A b a , b ∈ R ,F ( ν ) = u , F ( ν ) = u , (3)where the matrix field A = A ( x ) depends on the position x = ( x , x ) A = (cid:18) x − ν − ( x − ν ) x − ν x − ν (cid:19) . A mapping F as above is a composition of a translation, an homotopy and a (linearized)rotation; therefore mappings of this type include, in particular, linearized rigid body motions(that is the kernel of the symmetric gradient operator). Since the map F depends on u , u wewill denote it also by F u , u when we need to be explicit on such dependence. We will give anexplicit expression of F in Lemma 3.1 below.Then, the space R h ( e ) of two-component vector fields living on e is defined as the restrictionto e of all possible maps F : R h ( e ) = n F u , u | e : u , u ∈ R o . The following lemma follows by direct computations (the simple proof is not shown).
Lemma 3.1.
Given u , u ∈ R , the coefficients a , b of the map F = F u , u in (3) are a = u , b = B − ( u − u ) , B = (cid:18) ν − ν − ( ν − ν ) ν − ν ν − ν (cid:19) , B − = 1 k ν − ν k (cid:18) ν − ν ν − ν − ( ν − ν ) ν − ν (cid:19) .
4e recall γ e : [ a, b ] → e denotes the restriction of the parametrization describing Γ to theedge e . We introduce another space of vector fields living on e , defined as the push forward ofthe P k polynomials living on e that vanish on a, b : B h ( e ) = [ e P k ( e ) ∩ H ( e )] = n p ◦ γ − e : p ∈ ( P k [ a, b ]) , p ( a ) = p ( b ) = o . Finally, the virtual edge space is V h ( e ) = R h ( e ) ⊕ B h ( e ) . (4) Remark . A set of degrees of freedom for the space V h | e is simply given by k + 1 distinctpointwise evaluations along the edge, including the two extrema. This is easy to check, exploit-ing the fact that the evaluation operators at the two endpoints constitutes a set of degrees offreedom for R h ( e ) that vanish when computed on B h ( e ).As a byproduct of the above remark, the two spaces in (4) constitute a direct sum and thedimension of the space V h ( e ) is equal to k + 1. We are now ready to define the virtual space V h ( E ) on the (possibly curved) polygon E . V h ( E ) = n v ∈ [ H ( E ) ∩ C ( E )] : v | e ∈ [ P k ( e )] if e is straight , v | e ∈ V h ( e ) if e is curved , − ∆ v ∈ [ P k − ( E )] o . (5)Note that, due to Remark 3.2 below, the definition above is robust to the limit when a curvededge becomes straight. The standard virtual element theory [13], combined with Remark 3.1,yields that D1 , D2 and D3 of Sec. 3.1 is also a set of degrees of freedom for the novel space V h ( E ). As usual, the global space is obtained by a standard gluing procedure V h = n v ∈ H (Ω) : v | E ∈ V h ( E ) ∀ E ∈ Ω h o with the obvious extension of the local degrees of freedom to global ones. Remark . It can be checked that, whenever the edge e is straight and the associatedparametrization γ e is affine, one has V h ( e ) = P k ( e ). Therefore in such case we correctly recoverthe standard virtual space on straight edges. Remark . Also an enhanced version (see for instance [1]) of this space could be introduced.Indeed, since the enhancement of the space does not involve the boundary of the elements, it istrivial to build the enhanced version of the spaces V h ( E ) above combining the present resultswith those in [1, 14, 10]. Optimal approximation properties for the original space in the H norm were proved inTheorem 3.1 of [12]. These immediately extend (component-wise) to the vector valued version(2) here considered. On the other hand, the space (5) differs from (2). The only differencebetween the space (2) and (5) is the definition of the boundary space on curved edges (thatis V h ( e ) instead of [ e P k ( e )] ). Therefore we can directly apply the theory of [12] also in thisnew case, provided we can prove approximation estimates for the new curved edge space V h ( e ).That is, we need to prove a new counterpart of Lemma 3.2 in [12] for the novel space V h ( e ), thechallenge being that V h ( e ) is not a mapped polynomial space but stems instead from a “hybrid”definition (4). This is the topic of the present section.Let now a generic vector field u ∈ [ H ( e )] . Let the unique interpolant u RI ∈ R h ( e ) , u RI ( ν ) = u ( ν ) , u RI ( ν ) = u ( ν ) . (6) Lemma 4.1.
It exists a positive constant C, depending only on k and the parametrization γ ,such that | u RI | H m ( e ) ≤ C | u | H ( e ) ∀ m ∈ { , , .., k + 1 } . roof. By definition u RI = F | e , where F is given in (3) and coefficients determined by Lemma3.1 with u = u ( ν ) and u = u ( ν ). Let t represent the arc-length parametrization along theedge e , and % : [0 , h e ] → e the corresponding parametrization. Since all norms on R × areequivalent, it holds k B − k (cid:46) (cid:0) B − T : B − (cid:1) / = k ν − ν k − (cid:46) h − e . Then, by an Holder inequality, the euclidean norm of b satisfies k b k ≤ k B − k k u − u k = k B − k k u ( ν ) − u ( ν ) k (cid:46) h − e k Z ν ν ∂ t u ( t )d t k (cid:46) h − / e | u | H ( e ) . (7)The expression of u RI in term of the arc-length parametrization is, by definition, u RI ( t ) = a + A b with A = (cid:18) % ( t ) − ν − ( % ( t ) − ν ) % ( t ) − ν % ( t ) − ν (cid:19) . Therefore a direct derivation and (7) yield, for all t ∈ (0 , h e ) and m ∈ { , , .., k + 1 } , k ∂ mt u RI ( t ) k = k (cid:18) ∂ mt % ( t ) − ∂ mt % ( t ) ∂ mt % ( t ) ∂ mt % ( t ) (cid:19) b k (cid:46) k ∂ mt % ( t ) kk b k (cid:46) h − / e | u | H ( e ) . (8)Note that the term above k ∂ mt % ( t ) k (cid:46)
1, independently of the mesh, because the curve Γ isfixed once and for all. The Holder inequality together with bound (8) immediately grants | u RI ( t ) | H m ( e ) = k ∂ mt u RI ( t ) k L ( e ) (cid:46) h / e h − / e | u | H ( e ) = | u | H ( e ) . We now need to introduce some notation. Given a curved edge e , we denote by { x ei } k +1 i =1 the k +1 distinct points on the edge associated to the DoFs of V h ( e ) = V h | e (the cases i = 1 , i = k +1representing the two extrema of e ). Given any sufficiently regular vector field v living on e ,we denote by v II its interpolant in B h ( e ), that is the unique function in B h ( e ) that satisfies v II ( x ei ) = v ( x ei ) for i = 2 , .., k . We moreover denote by v I the interpolant of v in the space V h ( e ), that is the unique function in V h ( e ) that satisfies v I ( x ei ) = v ( x ei ) for i = 1 , .., k + 1. It isthen easy to check that, for any sufficiently regular vector field u living on e , one has u I = u RI + ( u − u RI ) II , (9)where we recall u RI was defined in (6).We now prove our counterpart of Lemma 3.2 in [12]. Proposition 4.1.
Let u ∈ [ H s ( e )] with ≤ s ≤ k + 1 . Then it holds | u − u I | H m ( e ) (cid:46) h s − mE k u k H s ( e ) for all real numbers ≤ m ≤ s .Proof. Recalling (9) we have | u − u I | H m ( e ) = | ( u − u RI ) − ( u − u RI ) II | H m ( e ) . (10)Let now ˜ p ∈ e P k ( e ) be the unique interpolant of ( u − u RI ) at the points { x ei } k +1 i =1 . It was provedin Lemma 3.2 of [12] that the interpolation into mapped polynomials e P k ( e ) enjoys standardapproximation properties; therefore | ( u − u RI ) − ˜ p | H m ( e ) (cid:46) h s − me k u − u RI k H s ( e ) . (11)Since the interpolation on mapped polynomials preserves constant functions (that is, the in-terpolation of a constant is the same constant), it is not restrictive to assume that u − u RI has6ero average. Therefore a standard Poincaré inequality for zero average functions yields that(11) can be extended to the slightly sharper bound (where the L part of the norm in the righthand side does not appear) | ( u − u RI ) − ˜ p | H m ( e ) (cid:46) h s − me s X j =1 | u − u RI | H j ( e ) . (12)We now note that, since ( u − u RI ) vanishes at the endpoints of e , it holds ( u − u RI ) II = ˜ p .Therefore, combining this identity with (10) and (12), we obtain | u − u I | H m ( e ) (cid:46) h s − me s X j =1 | u − u RI | H j ( e ) . (13)From (13), we conclude by the triangle inequality and using Lemma 4.1 | u − u I | H m ( e ) (cid:46) h s − me (cid:16) m X j =1 | u | H j ( e ) + | u RI | H j ( e ) (cid:17) (cid:46) h s − me m X j =1 | u | H j ( e ) (cid:46) h s − mE k u k H s ( e ) . Combining the results in [12] with Proposition 4.1 (that substitutes Lemma 3.2 in [12]) weimmediately have the following approximation property for the global space V h associated to(5). Proposition 4.2.
Let u ∈ H s (Ω) and < s ≤ k + 1 . Let all elements E ∈ Ω h be star shapedwith respect to a ball of radius uniformly comparable to h E . Then, there exists u I ∈ V h suchthat | u − u I | H (Ω) ≤ Ch s − k u k H s (Ω) where the constant C depends on the degree k , the “chunkiness” constant associated to the shaperegularity condition above and the parametrization γ . The condition s > above could be generalized to s >
1, but we here prefer to follow thesame choice as in [12]. See Remark 3.8 in [12].
The discretization of the problem is a combination of the scheme proposed in [15, 6, 5] forthe case with standard straight edges and the curved-edge technology introduced in [12] for amodel linear diffusion problem. Clearly, we here have the possibility of using the (vector valuedversion) of the space in [12] (briefly reviewed in Sec. 3.1 here above) or the novel space proposedin Sec. 3.2. The construction below stays the same for both spaces.We start by introducing the following projection operator that is used to compute, on eachmesh element E , a polynomial strain function. Let [ P k − ( E )] × denote the set of polynomialsymmetric tensors of degree k − E . Given v h ∈ V h and E ∈ Ω h , the operatorΠ ε : V h → [ P k − ( E )] × is defined by Π ε ( v h ) ∈ [ P k − ( E )] × Z E Π ε ( v h ) : p k − = Z E ε ( v h ) : p k − ∀ p k − ∈ [ P k − ( E )] × , where ε ( v h ) denotes as usual the symmetric gradient of v h . On each element E , the polynomialtensor Π ε ( v h ) represents the local approximation of the strains. Note that the above operatoris computable. Indeed an integration by parts shows that Z E ε ( v h ) : p k − = − Z E v h · (div p k − ) + Z ∂E v h · ( p k − n E ) . p k − is a vector polyno-mial of degree k − v h . The second term onthe right hand side can be computed since we have complete knowledge of v h on the boundaryof E . Note that all this computations clearly require the integration of known functions on acurved element and a curved boundary; those can be accomplished as shown for instance in[12, 36, 37, 38], see also Sec. 6.We can now describe the proposed numerical method. For simplicity of exposition wewill focus on the (nonlinear) elastic case (1), the inelastic case being treated analogously ascommented in the notes below. We start by introducing the local discrete“energy” form (forany E ∈ Ω h and u h , v h ∈ V h ) divided into a consistency and a stability part a Eh ( u h , v h ) = Z E σ (Π ε ( u h )) : Π ε ( v h ) + s E ( u h ; ( I − Π) u h , ( I − Π) u h ) , (14)where the above quantities and operators are described below. The operator σ : R × → R × represents the “black box” constitutive law associating strains to stresses (note that this rulemay depend on the position x ∈ Ω). The operator Π : V h ( E ) → [ P k ( E )] can be chosen as anyprojection operator on vector-valued polynomials of degree k , for instance one that minimizesthe distance of the euclidean norm of the degree of freedom values (such particular choice hasthe advantage of being very simple to code, see for instance [6]). The stabilization form can betaken, for instance, as s E ( u h ; ( I − Π) u h , ( I − Π) u h ) = α ( u h ) dofs X i =1 (cid:16) dof i ( u h − Π u h ) (cid:17)(cid:16) dof i ( v h − Π v h ) (cid:17) (15)where the dof i ( · ) symbol denotes evaluation at the i th local degree of freedom. The positivescalar α depending on u h is important in order to scale the discrete form correctly with respectto the constitutive law [15]. It can, for instance, be taken as α ( u h ) = (cid:13)(cid:13)(cid:13)(cid:13) ∂ σ ∂ ε (Π u h ( x b )) (cid:13)(cid:13)(cid:13)(cid:13) (16)with k · k denoting any (chosen) norm on the tangent tensor, that is here computed for Π u h evaluated at the baricenter x b of the element. Note that the above stabilization, that is quiteawkward to write on paper, is instead very simple to code since it is directly based on the degreeof freedom values (that is what the code operates on).In the numerical test of Sec. 6 we adopt a more advanced stabilization, inspired from the socalled “D-recipe” in [11], to whom we refer for more details on the underlying idea. We take s E ( u h ; ( I − Π) u h , ( I − Π) u h ) = dofs X i =1 α i ( u h ) (cid:16) dof i ( u h − Π u h ) (cid:17)(cid:16) dof i ( v h − Π v h ) (cid:17) , (17)with the positive scalars α i ( u h ) = max { α ( u h ) , M ii } where α ( u h ), defined in (16), is computed at the previous converged incremental step, and M ii denotes the i th diagonal entry of the consistent tangent matrix (which is the consistent tangentmatrix associated to the first term in (14)) computed in the centroid of the element at theprevious converged incremental step. Finally, we note that other options for s E could also beadopted and that the method is quite robust in this respect.Let V h represent the subspace of V h including the Dirichlet-type boundary conditions forthe problem under study, that we consider homogeneous for simplicity of exposition. Let f h denote the piecewise polynomial function of degree k − E ,by the L ( E ) projection of the volume loading f on [ P k − ( E )] (in the case k = 1 we takethe projection on [ P ( E )] ). Assuming an incremental loading procedure with M steps, thenumerical method at each incremental step i = 1 , .., M is written as Find u i ∈ V h such that X E ∈ Ω h a Eh ( u ih , v h ) = iM Z E f h · v h ∀ v h ∈ V h (18)8nd the final solution to the problem taken as u h = u Mh . At each incremental step, a standardNewton-Raphson procedure is applied to solve the nonlinear problem.Also the right hand side above is computable for k ≥ k = 1 the above integral can be approximated by avertex-based quadrature rule [13, 10, 15, 6] (see also the notes below).An analogous construction is applied in the presence on non-homogeneous boundary condi-tions (such as boundary tractions or enforced displacements), that are also treated incremen-tally. A series of remarks are in order. Alternative calculation of α . For M > α in (15) can also be computed at the previousconverged incremental step. Indeed, since the method is not sensible to such parameter, thiswill anyway allow for a satisfactory update of the material scaling factor. This choice carries theadvantage of leading to a simpler tangent matrix in the Newton iterations, since the derivativesof the α term disappear in the calculation. More accurate loading . By applying an enhanced construction for the space V h (see Remark 3.3)one can also compute a more accurate loading, using f h piecewise in [ P k ( E )] . The modificationis exactly the same as for standard straight-edged VEM and thus not detailed here. Inelastic case . The extension to the inelastic case is done exactly as in the Finite ElementMethod. Indeed, one simply needs to keep track of the history variables at the chosen integrationpoints (within each element) and apply at each point the constitutive law, exactly as in FEM(see for instance [15, 5]).
Simplicity of the curved case . Note that (at the practical level) the above construction followsthe same logic and structure as for the straight-edge case [15, 6, 5]. In the code, the maindifference is only the need to integrate along curved edges and on curved domains (that can behandled following the literature given above).
The present section has twofold purposes. First we give numerical evidence for error analysison linear elasticity benchmarks with available analytic solution. Then we present a coupleof classical benchmarks from computational mechanics which entail solution of equilibriumboundary value problems on curved domains in conjunction with nonlinear inelastic materialbehavior. Such examples numerically show the efficiency of the proposed curved methodologyin real scale structural simulations.In the remainder of the section, three virtual element approaches will be invoked accordingto the following acronyms: VEM co i.e. the original curvilinear VEM of Sec. 3.1, VEM cv i.e.the variant curvilinear VEM of Sec. 3.2, and VEM s i.e. the straight-edge VEM presented andvalidated in [6]. As already mentioned, we use the stabilization in Equation (17). Remark . Integration rules . In the following numerical tests, in order to compute thevolume integrals on (possibly) curved elements, we adopt the same integration algorithm usedin [12] (that in turn was taken from [36, 37, 38]) combined with a rotation of the cartesiancoordinates. Such algorithm allows, for any given order n ∈ N and element E , to generate a setof points and weights yielding an integration rule that is exact for polynomials of degree up to n . In the presented tests, when combined with a suitable coordinate rotation, such algorithmyielded (in all cases considered) points that are internal to the element and weights that arepositive. Nevertheless, such algorithm has the drawback of using many integration points andthus needs to be combined with a compression rule for the sake of efficiency. In order to keep thepositivity of the weights during the compression, we follow the compression rule of [9] insteadof the one used in [12]. The final number of integration points is equal to ( n + 1)( n + 2) /
2. Aninvestigation on the number n to be used in the VEM analysis is shown in Example 1. Finally,we point out that for all edge integrals (that are standard mapped one dimensional integralson the reference interval), one can use a classical Gauß-Lobatto integration rule. The minimalrequired number of points (in analogy with the straight edge case) for a scheme of order k are k + 1 points, including the two extrema, but a higher number of points can be used to obtaina safer rule. 9 .1 Numerical tests for elastic materials In the present section we consider 2D elasticity problems with known solution. The firstnumerical test investigates the convergence of the scheme and develops some comparisons. Thesecond test gives numerical evidence that the novel virtual element space with curved edgesproposed in Sec. 3.2 contains rigid body motions.In both examples we consider as domain Ω the circle centered at the origin with unit radius,and the following meshes: • quadN : a mesh composed by squares and a polygon at the center, see Fig. 2 (a); • rhexN : a mesh composed by regular hexagons inside the domain and arbitrary polygonsclose to the boundary, see Fig. 2 (b); • voroN : a mesh composed by Voronoi cells, see Fig. 2 (c).All of the above polygons abutting the circle boundary have curved edges. Here N refers to thenumber of polygons in the mesh. The last two meshes are generated with Polymesher [40].(a) (b) (c)Figure 2: Meshes used in Examples 1 and 2: (a) quad80 , (b) rhex70 and (c) voro128 .In order to develop a convergence analysis, in the numerical tests we re-define the mesh-size h as the mean value of element diameters, i.e. h = 1 N X E ∈ Ω h h E , where h E is the diameter of element E ∈ Ω h . For each mesh type, we consider a set ofuniformly refined meshes and adopt the following error measures on displacement and strainfield, respectively: • l ∞ -error on displacement field: e u := k u − u h k l ∞ ( N h ) k u k l ∞ ( N h ) , with the norm k u − u h k l ∞ ( N h ) := max x ∈N h k u ( x ) − u h ( x ) k defined through the values of the solution on the mesh skeleton, i.e. at the vertex/edgenodes x ∈ N h , and where k · k is the Euclidean norm. Although we do not have a rigoroustheoretical analysis for this norm, in analogy with Finite Elements the expected trend ofsuch error is O ( h k +1 ) for VEM elements of order k and sufficiently regular solutions.10 L -error on strain field: in this case we exploit the projection operator Π ε and wecompute e ε := k ε ( u ) − Π ε u h k L (Ω h ) k ε ( u ) k L (Ω h ) . On the basis of the interpolation results of Proposition 4.2 the expected trend of sucherror is O ( h k ) for VEM elements of order k and sufficiently regular solutions. Example 1. Convergence analysis.
In this example we develop a convergence analysisfor the proposed VEM curved elements, also comparing integration rules of different order.Moreover, we compare the results obtained via the proposed approach (which can deal withcurved elements and thus reproduce exactly the geometry of interest) and the standard VEMthat uses straight polygonal meshes (and thus the geometry boundary is approximated by apiecewise linear curve). In the present tests we do not distinguish between the scheme VEM co of Sec. 3.1 and the scheme VEM cv of Sec. 3.2 since they yielded almost identical results.We consider a nonlinear elastic material, characterized by the Hencky-von Mises constitutivelaw [31] σ ( x, ∇ u ) := ˜ λ (dev( ε ( u ))) tr( ε ( u )) I + ˜ µ (dev( ε ( u ))) ε ( u )where˜ µ ( % ) := 34 (1 + (1 + % ) − / ) · MPa and ˜ λ ( % ) := 34 (1 − µ ( % )) · MPa ∀ % ∈ R + are the non-linear Lamé functions, ε ( u ) := ( ∇ u + ( ∇ u ) t ) is the small deformation straintensor, tr( τ ) is the trace operator anddev( τ ) := (cid:13)(cid:13)(cid:13)(cid:13)(cid:18) τ −
12 tr( τ ) I (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) F , is the Frobenius norm of the deviatoric part of the tensor τ .The volume force density field and the homogeneous Dirichlet boundary conditions arechosen in such a way that the solution of Problem 1 u ( x, y ) = (cid:18) sin( π ( x + y ))2 cos( π ( x + y )) (cid:19) . The problem is solved with two VEM approaches: the VEM co for curved geometry describedin Sec. 3.1 and the standard VEM s displacement-based VEM method [13, 15, 6, 5] with straightedges (that approximates the domain of interest by a straight edge polygonal). We reportnumerical evidence for method degrees k = 1 , , minimal quadrature rule we use an integration ruleof order n = 2 k − higher-order quadrature rule we instead use an integration rule of order n = 2 k , that is moreaccurate but more expensive (expecially in view of using the scheme for inelastic problems).Analogously, for edge integrals we use a Gauß-Lobatto rule with k + 1 points in the minimal quadrature case, and k + 2 points in the higher-order case.In Figs. 3 and 4 we show the convergence plots in terms of the above introduced error normsfor the curved VEM co elements, comparing the minimal and the higher-order rules. We cansee that in all cases the convergence slopes are the expected optimal ones ( O ( k k +1 ) for nodaldisplacement error and O ( h k ) for the strain L error). When considering the displacementerror, the higher-order rule yields slightly more accurate (although of the same convergencerate) results. Our conclusion is that, in general, the minimal rule seems the better choice sinceit provides similar results with fewer quadrature points.In Fig. 5 we depict the error curves for the straight VEM s elements, that is without an exactgeometry approximation. We use the high-order quadrature rule. Due to the error introducedby the geometry approximation, for k ≥ k andfiner meshes. -2 -1 -7 -6 -5 -4 -3 -2 -1 -2 -1 -7 -6 -5 -4 -3 -2 -1 (a) (b)Figure 3: Example 1 (curved VEM co elements): convergence plot for e u with minimal quadra-ture (a), and higher order quadrature, (b). -2 -1 -4 -3 -2 -1 -2 -1 -4 -3 -2 -1 (a) (b)Figure 4: Example 1 (curved VEM co elements): convergence plot for e ε with minimal quadra-ture (a), and higher order quadrature, (b). Example 2. Rigid body motions.
In this paragraph we compare the virtual elementapproach VEM co of Sec. 3.1, and its variant VEM cv of Sec. 3.2. The main difference betweenthese two methods relies on the inclusion of the rigid body motion in the virtual displacementspace. Indeed, both of them are defined for curved elements and they will present an optimalconvergence rate, but only the latter method VEM cv will guarantee exact rigid body motionrepresentation in the solution.To give a numerical evidence of this feature, we consider a linear elastic problem, where thevolume force density is zero. As in the previous test, we consider a unit radius circle centeredat the origin, with Dirichlet boundary conditions set in such a way that the exact solution is12 -2 -1 -3 -2 -1 -2 -1 -3 -2 -1 (a) (b)Figure 5: Example 1 (straight VEM s elements): convergence with high order quadrature, error e u (a) and error plot e ε (b).the vector field (representing a rigid body motion) u ( x, y ) = (cid:18) − yx (cid:19) . (19)All volume and edge integrals are computed with a very high quadrature rule in this test, inorder to eliminate the effects of numerical quadrature in the comparison.In Table 1 we report the e u error, k = 1 , ,
3, for the two curvilinear VEM variants and somesample meshes. While it is evident that the VEM cv is exact up to machine precision, the errorlevels associated with the VEM co are much higher. This shows that the virtual element spacein the new variant (VEM cv ) contains the solution of the problem at hand, see Equation (19),while the standard space (VEM co ) does not contains such solution. quad481 rhex626 voro1800 VEM co VEM cv VEM co VEM cv VEM co VEM cvk = 1 2.2018e-03 1.0847e-15 4.0851e-04 3.2530e-13 4.4464e-04 1.0139e-15 k = 2 1.0379e-06 5.4964e-15 4.1290e-08 1.6785e-13 4.3888e-08 6.9502e-15 k = 3 5.6965e-08 2.4372e-14 1.9404e-09 9.8653e-14 1.6537e-09 1.0776e-14Table 1: Example 2 (rigid body motions): error e u for different mesh types and VEM approxi-mation degree. In this section we consider two classical benchmark examples involving inelastic materials.Since we do not have a reference solution for such problems, we use the software FEAP on avery fine mesh to validate the results.
Example 3. Thick-walled viscoelastic cylinder subjected to internal pressure.
Weconsider a classical numerical test regarding a thick-walled cylinder characterized by a viscoelas-tic constitutive response [45]. The cylinder has inner [resp. outer] radius R i = 2 [ R o = 4] andis subjected to uniform pressure p on the inner surface, see Fig. 6. The material is isotropicand obeys a viscoelastic constitutive law as outlined Sec. 3.1 of reference [5]. The materialproperties are set assuming M = 1 and λ ≡ λ = 1, i.e. a standard linear solid is considered[46]. Young’s modulus and Poisson’s ratio are E = 1000, ν = 0 .
3, respectively. Two sets of13iscoelastic parameters are adopted for the present analysis, i.e. ( µ , µ ) ve1 = (0 . , .
99) and,( µ , µ ) ve2 = (0 . , . K/G (0) = 2 .
167 and for long time loading,say at t = 8, by K/G (8) = 216 .
7, which indicates a nearly incompressible behavior for sus-tained loading (for instance, at t = ∞ the Poisson ratio results 0 . t = ∞ after loading is applied. Plane strain assumptionis applied in this numerical simulation.For symmetry, only a quarter of the cylinder cross section is meshed, as reported in Fig. 7,where structured convex curvilinear quadrilaterals (a), and unstructured convex quadrilaterals(b) meshes are portrayed, respectively. Zero normal displacement is enforced along straightedges. The structural response for an internal pressure p = 10 applied at t = 0 and kept constantuntil t = 20 is computed through 20 equally spaced time instants by using the generalizedMaxwell model in Prony series form, which, albeit intrinsically three-dimensional, translatesnaturally in the present two-dimentional context [45].A displacement accuracy analysis is drawn between the proposed curved VEM co / cv for-mulations, for k = 1 , ,
3, with a reference solution reported for comparison. Such referencesolution is obtained with quadratic quadrilateral displacement based finite elements with ninenodes Q ve1 and ve2 introduced above. We here show only the results for the(more demanding) distorted meshes (b), the ones for other case being analogous. It is observedthat compared VEM formulations present a response in excellent agreement with the referencesolution also for non structured and distorted discretizations of the domain, thus proving theefficiency of the proposed formulation.We consider a quadrilateral mesh and we take a subset of quadrature points deployed alonga radial direction, see the quadrature points high-lighted with crosses in Fig. 9 (a). Accordingto the theory proposed in [45], the radial stress component, σ % , has a regular trend from -10to 0, see Fig. 4.5 (b) of [45]. In such example, standard displacement based finite elementsoften display an irregular trend characterized by spurious oscillations [45]. On the contrary, theproposed virtual element method results free from spurious oscillations and it is better alignedwith the reference trend, compare Fig. 9 (b) of the present paper and Fig. 4.5 (b) of [45].We finally compare again the accuracy of the novel VEM formulation with curved edgesand the original straight edge formulation. In both cases we adopt quite coarse meshes: a meshmade by a single quadrilateral element ( quad1 ), and meshes composed by 2 ×
2, 4 × × quad4 , quad16 and quad64 respectively). It is emphasized that, in thecurved case, the elements of the mesh have a curved boundary and thus no error is introducedin reproducing the geometry, while in the straight case the elements are straight quadrilateralshence introducing a rectification error for the domain boundary approximation. The problemis the same described at the start of this section and the reference solution was obtained withFEAP as described above.In Table 2 we report the error in the radial displacement at points A and B (for the finaltime instant t = 20), comparing, for various meshes and polynomial orders k , the straightcase with the curved case. The advantage of the curved formulation (also considering that thecomputational cost is essentially equivalent to that of the straight formulation) appears clearly,in particular for higher values of k and finer meshes. Example 4. Perforated plastic plate.
A rectangular strip with width of 2 L = 200 mmwidth and length of 2 H = 360 mm with a circular hole of 2 R = 100 mm diameter in itscenter is considered (see Fig. 10). Material response here follows classical von Mises plasticconstitutive model, with material parameters: E = 7000 kg / mm , ν = 0 .
3, and yield stress σ y , = 24 . / mm [46]. Plane strain assumption is assumed, and a standard backwardEuler scheme with return map projection is used for stress and material moduli computationat the quadrature point level [35]. Displacement boundary restraints are prescribed for normal14igure 6: Example 3. Thick-walled viscoelastic cylinder subjected to internal pressure. Geom-etry, boundary conditions, applied load.(a) (b)Figure 7: Example 3. Thick-walled viscoelastic cylinder subjected to internal pressure. (a)Structured quadrilateral mesh composed by 272 elements. (b) Unstructured quadrilateral meshcomposed by 324 elements.components on symmetry boundaries and on top and lateral boundaries. Loading is appliedby a uniform normal displacement δ = 2 mm with 400 equal increments on the upper edge, seeFig. 10.Owing to symmetry, only one quadrant of the perforated strip is discretized, as shown inFig. 11, with structured quadrilaterals, and a centroid based Voronoi tessellation, respectively.The simulation refers to virtual elements of order k = 2 ,
3, in the straight and curved variant,15 (a) (b)Figure 8: Example 3. Thick-walled viscoelastic cylinder with internal pressure. Integrationstep vs. radial displacement curves for control points A (higher curve), and B (lower curve).(a) case ( µ , µ ) ve1 = (0 . , . µ , µ ) ve2 = (0 . , . (a) (b)Figure 9: Example 3. Thick-walled viscoelastic cylinder with internal pressure. (a) Location ofradial quadrature nodes. (b) Radial plot for σ % stress component.16 rror on the radial displacement for point A k = 1 k = 2 k = 3VEM s VEM cv VEM s VEM cv VEM s VEM cv quad1 quad4 quad16 quad64 Error on the radial displacement for point B k = 1 k = 2 k = 3VEM s VEM cv VEM s VEM cv VEM s VEM cv quad1 quad4 quad16 quad64 A and B .respectively.Accuracy and robustness is assessed by plotting the structural response of the structure,in terms of force reaction sum at the imposed displacement top edge vs. integration step,which seems correct for all compared methods and mesh types, showing no significant spuriouslocking phenomena (at least for the presented problem and polynomial degrees). The presentbenchmark confirms the proposed novel curved VEM methodology as a powerful tool to beimplemented in a standard nonlinear FEM structural analysis platform.Figure 10: Example 4. Perforated plastic plate. Geometry, boundary conditions, loading.17a) (b)Figure 11: Example 4. Perforated plastic plate. (a) Structured quadrilateral mesh. (b) Voronoimesh. We generalized the curvilinear Virtual Element technology to generic 2D solid mechanicsproblems in the small strain regime. The proposed approach can accept a generic black-box(elastic or inelastic) constitutive algorithm and, in addition, can make use of polygonal mesheswith curved edges thus leading to an exact approximation of the geometry. We also introduceda novel Virtual Element space for displacements on curvilinear elements that, differently fromthe original one, contains also rigid body motions. After presenting rigorous theoretical inter-polation properties for the new space, we developed an extensive numerical test campaign toassess the behavior of the scheme. The campaign included a convergence analysis on a nonlinearelastic problem, a patch test to verify rigid body motions, standard benchmarks with inelasticmaterials, a study on the integration points rule to be adopted, and a comparison among curvi-linear and standard Virtual Elements. It is noted that a curvilinear implementation does notimply any significant complication with respect to standard VEM technology and can hence beimplemented into existing codes at very limited expense. Numerical results are very promisingand indicate the proposed curvilinear VEM as a viable and competitive technology.
Acknowledgements
E. Artioli gratefully acknowledges the partial financial support of the University of RomeTor Vergata Mission Sustainability Programme through project SPY-E81I18000540005; and thepartial financial support of PRIN 2017 project ”3D PRINTING: A BRIDGE TO THE FUTURE(3DP_Future). Computational methods, innovative applications, experimental validations ofnew materials and technologies.”, grant 2017L7X3CS_004.L. Beirão da Veiga and F. Dassi were partially supported by the European Research Councilthrough the H2020 Consolidator Grant (grant no. 681162) CAVE, Challenges and Advance-ments in Virtual Elements. This support is gratefully acknowledged.
References [1] B. Ahmad, A. Alsaedi, F. Brezzi, L. D. Marini, and A. Russo. Equivalent projectors forvirtual element methods.
Comput. Math. Appl. , 66(3):376–391, 2013.18
100 200 300 40000.511.522.5 10 (a) (b)Figure 12: Example 4. Perforated plastic plate. Structural response. (a) Quadrilateral mesh.(b) Voronoi mesh.[2] A. Anand, J. S. Ovall, S. Reynolds, and S. Weiß er. Trefftz Finite Elements on CurvilinearPolygons. Preprint arXiv:1906.09015.[3] P. F. Antonietti, G. Manzini, and M. Verani. The fully nonconforming virtual elementmethod for biharmonic problems. Math. Models Methods Appl. Sci. , 28(02):387–407, 2017.[4] E. Artioli. Asymptotic homogenization of fibre-reinforced composites: a virtual elementmethod approach.
Meccanica , 53(6):1187–1201, 2018.[5] E. Artioli, L. Beirão da Veiga, C. Lovadina, and E. Sacco. Arbitrary order 2D virtualelements for polygonal meshes: Part II, inelastic problem.
Comput. Mech. , 60(4):643–657,2017.[6] E. Artioli, L. Beirão da Veiga, C. Lovadina, and E. Sacco. Arbitrary order 2D virtualelements for polygonal meshes: Part I, elastic problem.
Comput. Mech. , 60(3):355–377,2017.[7] E. Artioli, S. de Miranda, C. Lovadina, and L. Patruno. A family of virtual elementmethods for plane elasticity problems based on the Hellinger-Reissner principle.
Comput.Methods Appl. Mech. Engrg. , 340:978 – 999, 2018.[8] E. Artioli, S. Marfia, and E. Sacco. High-order virtual element method for the homoge-nization of long fiber nonlinear composites.
Comput. Methods Appl. Mech. Engrg. , 341:571– 585, 2018.[9] E. Artioli, A. Sommariva, and M. Vianello. Algebraic cubature on polygonal elements witha circular edge.
Submitted , 2019.[10] L. Beirão da Veiga, F. Brezzi, L. D. Marini, and A. Russo. The hitchhiker’s guide to thevirtual element method.
Math. Models Methods Appl. Sci. , 24(08):1541–1573, 2014.[11] L. Beirão da Veiga, F. Dassi, and A. Russo. High-order virtual element method on poly-hedral meshes.
Comput. Math. Appl. , 74(5):1110–1122, 2017.[12] L. Beirão da Veiga, A. Russo, and G. Vacca. The virtual element method with curvededges.
ESAIM: M2AN , 53(2):375–404, 2019.[13] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo. Basicprinciples of virtual element methods.
Math. Models Methods Appl. Sci. , 23(1):199–214,2013. 1914] L. Beirão da Veiga, F. Brezzi, L. D. Marini, and A. Russo. Virtual element method forgeneral second-order elliptic problems on polygonal meshes.
Math. Models Methods Appl.Sci. , 26(04):729–750, 2016.[15] L. Beirão da Veiga, C. Lovadina, and D. Mora. A virtual element method for elastic andinelastic problems on polytope meshes.
Comput. Methods Appl. Mech. Eng. , 295:327–346,2015.[16] M. F. Benedetto, A. Caggiano, and G. Etse. Virtual elements and zero thickness interface-based approach for fracture analysis of heterogeneous materials.
Comput. Methods Appl.Mech. Engrg. , 338:41–67, 2018.[17] S. Bertoluzza, M. Pennacchio, and D. Prada. High order VEM on curved domains. PreprintarXiv:1811.04755.[18] S.O.R. Biabanaki, A.R. Khoei, and P. Wriggers. Polygonal finite element methods forcontact-impact problems on non-conformal meshes.
Comput. Methods Appl. Mech. Engrg. ,269:198–221, 2014.[19] L. Botti and D. A. Di Pietro. Assessment of Hybrid High-Order methods on curved meshesand comparison with discontinuous Galerkin methods.
J. Comput. Phys. , 370:58–84, 2018.[20] F. Brezzi and L. D. Marini. Virtual element methods for plate bending problems.
Comput.Methods Appl. Mech. Engrg. , 253:455–462, 2013.[21] H. Chi, L. Beirão da Veiga, and G.H. Paulino. Some basic formulations of the virtualelement method (VEM) for finite deformations.
Comput. Methods Appl. Mech. Engrg. ,318:148–192, 2017.[22] H. Chi, Beirão da Veiga L., and Paulino G.H. A simple and effective gradient recoveryscheme and a posteriori error estimator for the Virtual Element Method (VEM).
Comput.Methods Appl. Mech. Engrg. , 347:21 – 58, 2019.[23] H. Chi, C. Talischi, O. Lopez-Pamies, and G. H. Paulino. Polygonal finite elements forfinite elasticity.
Int. J. Numer. Meth. Eng. , 101(4):305–328, 2015.[24] P. G. Ciarlet.
The finite element method for elliptic problems , volume 40. Siam, 2002.[25] J. A. Cottrell, T. J. R. Hughes, and Y. Bazilevs.
Isogeometric analysis: toward integrationof CAD and FEA . John Wiley & Sons, 2009.[26] L. Beirão da Veiga, F. Brezzi, and L. D. Marini. Virtual elements for linear elasticityproblems.
SIAM J. Numer. Anal. , 51(2):794–812, 2013.[27] F. Dassi, C. Lovadina, and M. Visinoni. A three-dimensional Hellinger-Reissner virtualelement method for linear elasticity problems. Preprint arXiv:1906.06119.[28] D. A. Di Pietro and A. Ern. A hybrid high-order locking-free method for linear elasticityon general meshes.
Comput. Methods Appl. Mech. Engrg. , 283:1–21, 2015.[29] A. L. Gain, G. H. Paulino, L. S. Duarte, and I.F.M. Menezes. Topology optimization usingpolytopes.
Comput. Methods Appl. Mech. Engrg. , 293:411–430, 2015.[30] A. L. Gain, C. Talischi, and G. H. Paulino. On the virtual element method for three-dimensional linear elasticity problems on arbitrary polyhedral meshes.
Comput. MethodsAppl. Mech. Engrg. , 282:132–160, 2014.[31] G. N. Gatica, A. Márquez, and W. Rudolph. A priori and a posteriori error analyses ofaugmented twofold saddle point formulations for nonlinear elasticity problems.
Comput.Methods Appl. Mech. Engrg. , 264:23–48, sep 2013.2032] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs. Isogeometric analysis: CAD, finiteelements, NURBS, exact geometry and mesh refinement.
Comput. Methods Appl. Mech.Engrg. , 194(39-41):4135–4195, 2005.[33] G. H. Paulino and A. L. Gain. Bridging art and engineering using Escher-based virtualelements.
Struct. Multidiscip. O. , 51(4):867–883, 2015.[34] M. Pingaro, E. Reccia, and P. Trovalusci. Homogenization of random porous materials withlow-order virtual elements.
ASCE-ASME Journal of Risk and Uncertainty in EngineeringSystems, Part B: Mechanical Engineering , 5(3):030905, 2019.[35] J. C. Simo and T. J. R. Hughes.
Computational Inelasticity . Springer, 1998.[36] A. Sommariva and M. Vianello. Product Gauss cubature over polygons based on Green’sintegration formula.
BIT , 47(2):441–453, 2007.[37] A. Sommariva and M. Vianello. Gauss-green cubature and moment computation overarbitrary geometries.
J. Comput. App. Math. , 231:886–896, 2009.[38] A. Sommariva and M. Vianello. Compression of multivariate discrete measures and appli-cations.
Numer. Funct. Anal. Optim. , 36(9):1198–1223, 2015.[39] N. Sukumar and A. Tabarraei. Conforming polygonal finite elements.
Int. J. Numer. Meth.Eng. , 61(12):2045–2066, 2004.[40] C. Talischi, G. H. Paulino, A. Pereira, and I. F. M. Menezes. Polymesher: a general-purpose mesh generator for polygonal elements written in Matlab.
Struct. Multidiscip. O. ,45(3):309–328, Mar 2012.[41] R. L. Taylor and E. Artioli. VEM for inelastic solids.
Comput. Methods Appl. Sci. , 46:381–394, 2018.[42] P. Wriggers and B. Hudobivnik. A low order virtual element formulation for finite elasto-plastic deformations.
Comput. Methods Appl. Mech. Engrg. , 327:459–477, 2017.[43] P. Wriggers, W.T. Rust, and B.D. Reddy. A virtual element method for contact.
Comput.Mech. , 58(6):1039–1050, 2016.[44] P. Wriggers, W.T. Rust, B.D. Reddy, and B. Hudobivnik. Efficient virtual element formu-lations for compressible and incompressible finite deformations.
Comput. Mech. , 60(2):253–268, 2017.[45] O. C. Zienkiewicz, R. L. Taylor, and D. D. Fox.
The Finite Element Method for Solid andStructural Mechanics . Butterworth Heinemann, 2014.[46] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu.