The conforming virtual element method for polyharmonic and elastodynamics problems: a review
Paola Francesca Antonietti, Gianmarco Manzini, Ilario Mazzieri, Simone Scacchi, Marco Verani
TThe conforming virtual element method for polyharmonic andelastodynamics problems: a review
P. F. Antonietti a , G. Manzini b , I. Mazzieri a , S. Scacchi a , M. Verani a a MOX, Dipartimento di Matematica, Politecnico di Milano, Italy b IMATI, Consiglio Nazionale delle Ricerche, Pavia, Italy
Abstract
In this paper we review recent results on the conforming virtual element approximation of polyharmonic andelastodynamics problems. The structure and the content of this review is motivated by three paradigmaticexamples of applications: classical and anisotropic Cahn-Hilliard equation and phase field models for brittlefracture, that are briefly discussed in the first part of the paper. We present and discuss the mathematical detailsof the conforming virtual element approximation of linear polyharmonic problems, the classical Cahn-Hilliardequation and linear elastodynamics problems.
1. Introduction
In the recent years, there has been a tremendous interest to numerical methods that approximate partialdifferential equations (PDEs) on computational meshes with arbitrarily-shaped polytopal elements. One of themost successful method is the virtual element method (VEM), originally proposed in [14] for second-orderelliptic problems and then extended to a wide range of applications. The VEM was originally developedas a variational reformulation of the nodal mimetic finite difference (MFD) method [21, 34, 73] for solvingdiffusion problems on unstructured polygonal meshes. A survey on the MFD method can be found in the reviewpaper [70] and the research monograph [22]. The VEM inherits the flexibility of the MFD method with respectto the admissible meshes and this feature is well reflected in the many significant applications using polytopalmeshes that have been developed so far, see, for example, [51, 5, 40, 18, 17, 83, 93, 25, 76, 30, 82, 24, 26, 45, 46].Meanwhile, the mixed VEM for elliptic problems were introduced in setting a la‘
Raviart-Thomas in [16]and in a BDM-like setting in [35]. The nonconforming formulation for diffusion problems was proposedin [12] as the finite element reformulation of [69] and later extended to general elliptic problems [44, 29],Stokes problem [41], eigenvalue problems [60], and the biharmonic equation [8, 94]. equation [8]. Moreover,the connection between the VEM and the finite elements on polygonal/polyhedral meshes is thoroughlyinvestigated in [74, 43], between VEM and discontinuous skeletal gradient discretizations in [52], andbetween the VEM and the BEM-based FEM method in [42]. The VEM was originally formulated in [14]as a conforming FEM for the Poisson problem. Then, it was later extended to convection-reaction-diffusionproblems with variable coefficients in [17].The virtual element method combines a great flexibility in using polytopal meshes with a great versatilityand easiness in designing approximation spaces with high-order continuity properties on general polytopalmeshes. These two features turn out to be essential in the numerical treatment of the classical plate bendingproblem, for which a 𝐶 -regular conforming virtual element approximation has been introduced in [36, 50].Virtual elements with 𝐶 - regularity have been proposed to solve elliptic problems on polygonal meshes [24]and polyedral meshes in [20], the transmission eigenvalue problem in [78], the vibration problem of Kirchhoffplates in [77], the buckling problem of Kirchhoff-Love plates in [79]. The use of 𝐶 -virtual elements hasalso been employed in the conforming approximation of the Cahn-Hilliard problem [5] and the von Kármán Preprint submitted to Elsevier February 3, 2021 a r X i v : . [ m a t h . NA ] F e b quations [72], and in the context of residual based a posteriori error estimators for second-order ellipticproblems [25].Higher-order of regularity of the numerical approximation is also required when addressing PDEs withdifferential operators of order higher than two as the already mentioned biharmonic problem and the moregeneral case of the polyharmonic equations. An example of the latter is found in the work of Reference [9].In this paper we consider three paradigmatic examples of applications where the conforming discretizationrequires highly regular approximation spaces. The first two examples are the the classical and the anisotropicCahn-Hilliard equations, that are used in modeling a wide range of problems such as the tumor growth,the origin of the Saturn rings, the separation of di-block copolymers, population dynamics, crystal growth,image processing and even the clustering of mussels, see [5] and the references therein. The third examplehighlights the importance of coupling phase field equations with the elastodynamic equation in the context ofmodeling fracture propagation (see also [3] for a phase-field based VEM and the references therein). Thesethree examples motivate the structure of this review, where we consider the conforming virtual approximationof the polyharmonic equation, the classical Cahn-Hilliard equation and the time-dependent elastodynamicsequation.Historically, the numerical approximation of polyharmonic problems dates back to the eighties [32], andmore recently, this problem has been addressed in the context of the finite element method by [13, 64, 91,88, 59]. The conforming virtual element approximation of the biharmonic problem has been addressed in[36, 50]. while a non-conforming approximation has been proposed in [94, 8, 95]. In Section 2, we reviewthe conforming virtual element approximation of polyharmonic problems following [9]. A nonconformingapproximation is studied in [49].The Cahn-Hilliard equation involves fourth-order spatial derivatives and the conforming finite elementmethod is not really popular approach because primal variational formulations of fourth-order operatorsrequires the use of finite element basis functions that are piecewise-smooth and globally 𝐶 - continuous.Only a few finite element formulations exists with the 𝐶 -continuity property, see for example [57, 53], butin general, these methods are not simple and easy to implement. This high-regularity issue has successfullybeen addressed in the framework of isogeometric analysis [62]. The virtual element method provides a veryeffective framework for the design and development of highly regular conforming approximation, and inSection 3 we review the method proposed in [5].Alternative approaches are offered by nonconforming methods [54] or discontinuous methods [92]), butthese methods do not provide 𝐶 -regular approximations. Another common strategy employed in practice to solve the Cahn-Hilliard equation by finite elements resorts to mixed methods; see, e.g., [55, 56] and [65]for the continuous and discontinuous setting, respectively. Recently, mixed based discretization schemes onpolytopal meshes have been addressed in [47] in the context of the Hybrid High Order Method, and in [71] inthe context of the mixed Virtual Element Method. However, mixed finite element methods requires a biggernumber of degrees of freedom, which implies, as a drawback, an increased computational cost.Very popular strategies for numerically solving the time-dependent elastodynamics equations in the dis-placement formulation are based on spectral elements [66, 58], discontinuous Galerkin and discontinuousGalerkin spectral elements [86, 4, 11]. High-order DG methods for elastic and elasto-acoustic wave propaga-tion problems have been extended to arbitrarily-shaped polygonal/polyhedral grids [10, 6] to further enhancethe geometrical flexibility of the discontinuous Galerkin approach while guaranteeing low dissipation anddispersion errors. Recently, the lowest-order Virtual Element Method has been applied for the solution of theelastodynamics equation on nonconvex polygonal meshes [80, 81]. See also [15] for the approximation of thelinear elastic problem, [23] for elastic and inelastic problems on polytopal meshes, [90] for virtual elementapproximation of hyperbolic problems. In Section 4, we review the conforming virtual element method ofarbitrary order of accuracy proposed in [7]. 2 .1. Notation and technicalities Throughout the paper, we consider the usual multi-index notation. In particular, if 𝑣 is a sufficientlyregular bivariate function and 𝛼 = ( 𝛼 , 𝛼 ) a multi-index with 𝛼 , 𝛼 nonnegative integer numbers, thefunction 𝐷 𝛼 𝑣 = 𝜕 | 𝛼 | 𝑣 / 𝜕𝑥 𝛼 𝜕𝑥 𝛼 is the partial derivative of 𝑣 of order | 𝛼 | = 𝛼 + 𝛼 >
0. For 𝛼 = ( , ) ,we adopt the convention that 𝐷 𝛼 𝑣 coincides with 𝑣 . Also, for the sake of exposition, we may use theshortcut notation 𝜕 𝑥 𝑣 , 𝜕 𝑦 𝑣 , 𝜕 𝑥𝑥 𝑣 , 𝜕 𝑥𝑦 𝑣 , 𝜕 𝑦𝑦 𝑣 , to denote the first- and second-order partial derivatives alongthe coordinate directions 𝑥 and 𝑦 ; 𝜕 𝑛 𝑣 , 𝜕 𝑡 𝑣 , 𝜕 𝑛𝑛 𝑣 , 𝜕 𝑛𝑡 𝑣 , 𝜕 𝑡𝑡 𝑣 to denote the first- and second-order normaland tangential derivatives along a given mesh edge; and 𝜕 𝑚𝑛 𝑣 and 𝜕 𝑚𝑡 𝑣 to denote the normal and tangentialderivative of 𝑣 of order 𝑚 along a given mesh edge. Finally, let n = ( 𝑛 𝑥 , 𝑛 𝑦 ) and 𝝉 = ( 𝜏 𝑥 , 𝜏 𝑦 ) be the unitnormal and tangential vectors to a given edge 𝑒 of an arbitrary polygon P , respectively. We recall the followingrelations between the first derivatives of 𝑣 : 𝜕 𝑛 𝑣 ∇ 𝑣 𝑇 n = 𝑛 𝑥 𝜕 𝑥 𝑣 + 𝑛 𝑦 𝜕 𝑦 𝑣, 𝜕 𝜏 𝑣 = ∇ 𝑣 𝑇 𝝉 = 𝜏 𝑥 𝜕 𝑥 𝑣 + 𝜏 𝑦 𝜕 𝑦 𝑣, (1)and the second derivatives of 𝑣 : 𝜕 𝑛𝑛 𝑣 = n 𝑇 H ( 𝑣 ) n , 𝜕 𝑛𝜏 𝑣 = n 𝑇 H ( 𝑣 ) 𝝉 , 𝜕 𝜏𝜏 𝑣 = 𝝉 𝑇 H ( 𝑣 ) 𝝉 , (2)respectively, where the matrix H ( 𝑣 ) is the Hessian of 𝑣 , i.e., H ( 𝑣 ) = 𝜕 𝑥𝑥 𝑣 , H ( 𝑣 ) = H ( 𝑣 ) = 𝜕 𝑥𝑦 𝑣 , H ( 𝑣 ) = 𝜕 𝑦𝑦 𝑣 .We use the standard definitions and notation of Sobolev spaces, norms and seminorms [1]. Let 𝑘 be anonnegative integer number. The Sobolev space 𝐻 𝑘 ( 𝜔 ) consists of all square integrable functions with allsquare integrable weak derivatives up to order 𝑘 that are defined on the open bounded connected subset 𝜔 of R . As usual, if 𝑘 =
0, we prefer the notation 𝐿 ( 𝜔 ) . Norm and seminorm in 𝐻 𝑘 ( 𝜔 ) are denoted by || · || 𝑘,𝜔 and | · | 𝑘,𝜔 , respectively, and (· , ·) 𝜔 denote the 𝐿 -inner product. We omit the subscript 𝜔 when 𝜔 is thewhole computational domain Ω .Given the mesh partitioning Ω ℎ = { P } of the domain Ω into elements P , we define the broken (scalar)Sobolev space for any integer 𝑘 > 𝐻 𝑘 ( Ω ℎ ) = (cid:214) P ∈ Ω ℎ 𝐻 𝑘 ( P ) = (cid:8) 𝑣 ∈ 𝐿 ( Ω ) : 𝑣 | P ∈ 𝐻 𝑘 ( P ) (cid:9) , which we endow with the broken 𝐻 𝑘 -norm || 𝑣 || 𝑘,ℎ = ∑︁ P ∈ Ω ℎ || 𝑣 || 𝑘, P ∀ 𝑣 ∈ 𝐻 𝑘 ( Ω ℎ ) , (3)and, for 𝑘 =
1, with the broken 𝐻 -seminorm | 𝑣 | ,ℎ = ∑︁ P ∈ Ω ℎ ||∇ 𝑣 || , P ∀ 𝑣 ∈ 𝐻 ( Ω ℎ ) . (4)We denote the linear space of polynomials of degree up to ℓ defined on 𝜔 by P ℓ ( 𝜔 ) , with the usefulconventional notation that P − ( 𝜔 ) = { } . We denote the space of two-dimensional vector polynomials ofdegree up to ℓ on 𝜔 by (cid:2) P ℓ ( 𝜔 ) (cid:3) ; the space of symmetric 2 × ℓ on 𝜔 by P × ℓ, sym ( 𝜔 ) . Space P ℓ ( 𝜔 ) is the span of the finite set of scaled monomials of degree up to ℓ , that are3iven by M ℓ ( 𝜔 ) = (cid:26) (cid:18) x − x 𝜔 ℎ 𝜔 (cid:19) 𝛼 with | 𝛼 | ≤ ℓ (cid:27) , where • x 𝜔 denotes the center of gravity of 𝜔 and ℎ 𝜔 its characteristic length, as, for instance, the edge lengthor the cell diameter for 𝑑 = , • 𝛼 = ( 𝛼 , 𝛼 ) is the two-dimensional multi-index of nonnegative integers 𝛼 𝑖 with degree | 𝛼 | = 𝛼 + 𝛼 ≤ ℓ and such that x 𝛼 = 𝑥 𝛼 𝑥 𝛼 for any x ∈ R .We will also use the set of scaled monomials of degree exactly equal to ℓ , denoted by M ∗ ℓ ( 𝜔 ) and obtainedby setting | 𝛼 | = ℓ in the definition above.Finally, we use the letter 𝐶 in the estimates below to denote a strictly positive constant whose value canchange at any instance but that is independent of the discretization parameters such as the mesh size ℎ . Notethat 𝐶 may depend on the the polynomial order, on the constants of the model equations or the variationalproblem, like the coercivity and continuity constants, or even constants that are uniformly defined for thefamily of meshes of the approximation while ℎ →
0, such as the mesh regularity constant, the stabilityconstants of the discrete bilinear forms, etc. Whenever it is convenient, we will simplify the notation by usingexpressions like 𝑥 (cid:46) 𝑦 and 𝑥 (cid:38) 𝑦 to mean that 𝑥 ≤ 𝐶 𝑦 and 𝑥 ≥ 𝐶 𝑦 , respectively, 𝐶 being the generic constantin the sense defined above. Throughout the paper we assume that T = (cid:8) Ω ℎ (cid:9) ℎ is a family of decompositions of the computationaldomain Ω , where each mesh Ω ℎ is a collection of nonoverlapping polygonal elements P with boundary 𝜕 P ,such that ¯ Ω = (cid:100) P ∈ Ω ℎ ¯ P . Each mesh is labeled by the mesh size ℎ , the diameter of the mesh, defined as usualby ℎ = max P ∈ Ω ℎ ℎ P , where ℎ P = sup x , y ∈ P | x − y | . We assume the mesh sizes of family T form a countablesubset of H = ( , ∞) having zero as its unique accumulation point. We denote the set of mesh vertices v by V ℎ and the set of mesh edges 𝑒 by E ℎ Moreover, the symbol ℎ v is a characteristic length associated with eachvertex; more precisely, ℎ v is the average of the diameters of the polygons sharing vertex v . We consider thefollowing mesh regularity assumptions: (M) There exists a positive constant 𝛾 , mesh regularity constant , which is independent of ℎ (and P ) and suchthat for 𝐾 ∈ Ω ℎ there hold: • (M1) P is star-shaped with respect to every point of a ball of radius 𝛾ℎ P , where ℎ P is the diameterof P ; • (M2) for every edge 𝑒 of the cell boundary 𝜕 P of every cell P of Ω ℎ , it holds that ℎ 𝑒 ≥ 𝛾ℎ P , where ℎ 𝑒 denotes the length of 𝑒 .All the results contained in the rest of the paper are obtained under assumptions (M1) - (M2) .4 . The virtual element method for the polyharmonic problem Let Ω ⊂ R be a open, bounded, convex domain with polygonal boundary Γ . For any integer 𝑝 ≥
1, weintroduce the conforming virtual element method for the approximation of the following problem: (− Δ ) 𝑝 𝑢 = 𝑓 in Ω , (5a) 𝜕 𝑗𝑛 𝑢 = 𝑗 = , . . . , 𝑝 − Γ , (5b)(recall the conventional notation 𝜕 𝑛 𝑢 = 𝑢 ). Let 𝑉 ≡ 𝐻 𝑝 ( Ω ) = (cid:8) 𝑣 ∈ 𝐻 𝑝 ( Ω ) : 𝜕 𝑗𝑛 𝑣 = Γ , 𝑗 = , . . . , 𝑝 − (cid:9) . Denoting the duality pairing between 𝑉 and its dual 𝑉 (cid:48) by (cid:10) · , · (cid:11) , the variational formulation of the polyharmonicproblem (5) reads as: Find 𝑢 ∈ 𝑉 such that 𝑎 ( 𝑢, 𝑣 ) = (cid:10) 𝑓 , 𝑣 (cid:11) ∀ 𝑣 ∈ 𝑉, (6)where, for any nonnegative integer ℓ , the bilinear form is given by: 𝑎 ( 𝑢, 𝑣 ) = ∫ Ω ∇ Δ ℓ 𝑢 · ∇ Δ ℓ 𝑣 𝑑 x for 𝑝 = ℓ + , ∫ Ω Δ ℓ 𝑢 Δ ℓ 𝑣 𝑑 x for 𝑝 = ℓ. (7)Whenever 𝑓 ∈ 𝐿 ( Ω ) we have (cid:10) 𝑓 , 𝑣 (cid:11) = ( 𝑓 , 𝑣 ) = ∫ Ω 𝑓 𝑣 𝑑𝑉 𝑑 x . (8)where (· , ·) denotes the 𝐿 -inner product. The existence and uniqueness of the solution to (6) follows fromthe Lax-Milgram Theorem because of the continuity and coercivity of the bilinear form 𝑎 (· , ·) with respect to (cid:107) · (cid:107) 𝑉 = | · | 𝑝, Ω which is a norm on 𝐻 𝑝 ( Ω ) . Moreover, since Ω is a convex polygon, from [61] we know that 𝑢 ∈ 𝐻 𝑝 − 𝑚 ( Ω ) ∩ 𝐻 𝑝 ( Ω ) if 𝑓 ∈ 𝐻 − 𝑚 ( Ω ) , 𝑚 ≤ 𝑝 and it holds that || 𝑢 || 𝑝 − 𝑚 ≤ 𝐶 || 𝑓 || − 𝑚 . In the following, wedenote the coercivity and continuity constants of 𝑎 (· , ·) by 𝛼 and 𝑀 , respectively.Let P be a polygonal element and set 𝑎 P ( 𝑢, 𝑣 ) = ∫ P ∇ Δ ℓ 𝑢 · ∇ Δ ℓ 𝑣 𝑑 x for 𝑝 = ℓ + , ∫ P Δ ℓ 𝑢 Δ ℓ 𝑣 𝑑 x for 𝑝 = ℓ. For an odd 𝑝 , i.e., 𝑝 = ℓ +
1, a repeated application of the integration by parts formula yields 𝑎 P ( 𝑢, 𝑣 ) = − ∫ P Δ 𝑝 𝑢 𝑣 𝑑 x + ∫ 𝜕 P 𝜕 𝑛 ( Δ ℓ 𝑢 ) Δ ℓ 𝑣 𝑑𝑠 + ℓ ∑︁ 𝑖 = (cid:18)∫ 𝜕 P 𝜕 𝑛 ( Δ 𝑝 − 𝑖 𝑢 ) Δ 𝑖 − 𝑣 𝑑𝑠 − ∫ 𝜕 P Δ 𝑝 − 𝑖 𝑢 𝜕 𝑛 ( Δ 𝑖 − 𝑣 ) 𝑑𝑠 (cid:19) , (9)5hile, for an even 𝑝 , i.e., 𝑝 = ℓ , we have 𝑎 P ( 𝑢, 𝑣 ) = ∫ P Δ 𝑝 𝑢 𝑣 𝑑 x + ℓ ∑︁ 𝑖 = (cid:18)∫ 𝜕 P 𝜕 𝑛 ( Δ 𝑝 − 𝑖 𝑢 ) Δ 𝑖 − 𝑣 𝑑𝑠 − ∫ 𝜕 P Δ 𝑝 − 𝑖 𝑢 𝜕 𝑛 ( Δ 𝑖 − 𝑣 ) 𝑑𝑠 (cid:19) . (10) The conforming virtual element discretization of problem (6) hinges upon three mathematical objects: (1)the finite dimensional conforming virtual element space 𝑉 𝑝ℎ,𝑟 ⊂ 𝑉 ; (2) the continuous and coercive discretebilinear form 𝑎 ℎ (· , ·) ; (3) the linear functional (cid:10) 𝑓 ℎ , · (cid:11) .Using such objects, we formulate the virtual element method as: Find 𝑢 ℎ ∈ 𝑉 𝑝ℎ,𝑟 such that 𝑎 ℎ ( 𝑢 ℎ , 𝑣 ℎ ) = (cid:10) 𝑓 ℎ , 𝑣 ℎ (cid:11) ∀ 𝑣 ℎ ∈ 𝑉 𝑝ℎ,𝑟 . (11)The existence and uniqueness of the solution 𝑢 ℎ is again a consequence of the Lax-Milgram theorem.[33,Theorem 2.7.7, page 62]. For 𝑝 ≥ 𝑟 ≥ 𝑝 −
1, the local Virtual Element space on element P is defined by 𝑉 𝑝ℎ,𝑟 ( P ) = (cid:110) 𝑣 ℎ ∈ 𝐻 𝑝 ( P ) : Δ 𝑝 𝑣 ℎ ∈ P 𝑟 − 𝑝 ( P ) , 𝑣 ℎ ∈ P 𝑟 ( 𝑒 ) , 𝜕 𝑖𝑛 𝑣 ℎ ∈ P 𝑟 − 𝑖 ( 𝑒 ) , 𝑖 = , . . . , 𝑝 − ∀ 𝑒 ∈ 𝜕 P (cid:111) , with the conventional notation that P − ( P ) = { } . The virtual element space 𝑉 𝑝ℎ,𝑟 ( P ) contains the spaceof polynomials P 𝑟 ( P ) , for 𝑟 ≥ 𝑝 −
1. Moreover, for 𝑝 =
1, it coincides with the conforming virtualelement space for the Poisson equation [14], and for 𝑝 =
2, it coincides with the conforming virtual elementspace for the biharmonic equation [36]. The requirement 𝑣 ℎ ∈ 𝐻 𝑝 ( P ) implies that suitable compatibilityconditions for 𝑣 ℎ and its derivatives up to order 𝑝 − 𝑉 𝑝ℎ,𝑟 ( P ) through the following set of degrees of freedom : (D1) ℎ | 𝜈 | v 𝐷 𝜈 𝑣 ℎ ( v ) , | 𝜈 | ≤ 𝑝 − v of the polygonal boundary 𝜕 P ; (D2) ℎ − 𝑒 ∫ 𝑒 𝑞𝑣 ℎ 𝑑𝑠 for any 𝑞 ∈ P 𝑟 − 𝑝 ( 𝑒 ) and any edge 𝑒 of the polygonal boundary 𝜕 P ; (D3) ℎ − + 𝑗𝑒 ∫ 𝑒 𝑞𝜕 𝑗𝑛 𝑣 ℎ 𝑑𝑠 for any 𝑞 ∈ P 𝑟 − 𝑝 + 𝑗 ( 𝑒 ) , 𝑗 = , . . . , 𝑝 − 𝑒 of 𝜕 P ; (D4) ℎ − P ∫ P 𝑞 ℎ 𝑣 ℎ 𝑑 x for any 𝑞 ∈ P 𝑟 − 𝑝 ( P ) . 6 = , 𝑟 = 𝑝 = , 𝑟 = 𝑝 = , 𝑟 = 𝑝 = , 𝑟 = 𝑝 = , 𝑟 = 𝑝 = , 𝑟 = Figure 1: [Taken from [9]]. Edge degrees of freedom of the virtual element space 𝑉 𝑝ℎ,𝑟 ( P ) for the polyharmonic problem with 𝑝 = 𝑝 = 𝑝 = 𝑝 is theorder of the partial differential operator; 𝑟 = , , . . . , P 𝑟 ( P ) of the VEM space 𝑉 ℎ,𝑟 ( P ) . The (green) dots at the vertices represent the vertex values and each (red) vertex circle represents anorder of derivation. The (black) dot on the edge represents the moment of 𝑣 ℎ | 𝑒 ; the arrows represent the moments of 𝜕 𝑛 𝑣 ℎ | 𝑒 ; the doublearrows represent the moments of 𝜕 𝑛𝑛 𝑣 ℎ | 𝑒 . The corresponding internal degrees of freedom (D4) are absent in the case 𝑟 = 𝑝 −
1, whilereduce to a single one in the case 𝑟 = 𝑝 . Here, as usual, we assume that P − 𝑛 (·) = { } for 𝑛 ≥
1. Figure 1 illustrates the degrees of freedom on a givenedge 𝑒 for 𝑝 = , , 𝑟 = 𝑝 − , 𝑝 ; the correspondinginternal degrees of freedom (D4) are absent in the case 𝑟 = 𝑝 −
1, while reduce to a single one in the case 𝑟 = 𝑝 . Finally, we note that in general the internal degrees of freedom (D4) make it possible to define the 𝐿 -orthogonal polynomial projection of 𝑣 ℎ onto the space of polynomial of degree 𝑟 − 𝑝 .The dimension of 𝑉 𝑝ℎ,𝑟 ( P ) is 𝑑 ( 𝑉 𝑝ℎ,𝑟 ( P )) = 𝑝 ( 𝑝 + ) 𝑁 P + 𝑁 P 𝑝 − ∑︁ 𝑗 = ( 𝑟 − 𝑝 + 𝑗 + ) + ( 𝑟 − 𝑝 + )( 𝑟 − 𝑝 + ) , where 𝑁 P is the number of vertices, which equals the number of edges, of P .In [9], it is proved that the above choice of degrees of freedom is unisolvent in 𝑉 𝑝ℎ,𝑟 ( P ) .Building upon the local spaces 𝑉 𝑝ℎ,𝑟 ( P ) for all P ∈ Ω ℎ , the global conforming virtual element space 𝑉 𝑝ℎ,𝑟 is defined on Ω as 𝑉 𝑝ℎ,𝑟 = (cid:110) 𝑣 ℎ ∈ 𝐻 𝑝 ( Ω ) : 𝑣 ℎ | P ∈ 𝑉 𝑝ℎ,𝑟 ( P ) ∀ P ∈ Ω ℎ (cid:111) . (12)We remark that the associated global space is made of 𝐻 𝑝 ( Ω ) functions. Indeed, the restriction of a virtualelement function 𝑣 ℎ to each element P belongs to 𝐻 𝑝 ( P ) and glues with 𝐶 𝑝 − -regularity across the internalmesh faces. The set of global degrees of freedom inherited by the local degrees of freedom are: • ℎ | 𝜈 | v 𝐷 𝜈 𝑣 ℎ ( v ) , | 𝜈 | ≤ 𝑝 − v of Ω ℎ ; • ℎ − 𝑒 ∫ 𝑒 𝑞𝑣 ℎ 𝑑𝑠 for any 𝑞 ∈ P 𝑟 − 𝑝 ( 𝑒 ) and every interior edge 𝑒 ∈ E ℎ ;7 ℎ − + 𝑗𝑒 ∫ 𝑒 𝑞𝜕 𝑗𝑛 𝑣 ℎ 𝑑𝑠 for any 𝑞 ∈ P 𝑟 − 𝑝 + 𝑗 ( 𝑒 ) 𝑗 = , . . . , 𝑝 − 𝑒 ∈ E ℎ ; • ℎ − P ∫ P 𝑞𝑣 ℎ 𝑑 x for any 𝑞 ∈ P 𝑟 − 𝑝 ( P ) and every P ∈ Ω ℎ . In this section, we briefly discuss the possibility of introducing modified lowest order virtual elementspaces with a reduced number of degrees of freedom with respect to the corresponding lowest order ones thatwere introduced previously. The price we pay is a reduced order of accuracy since the polynomial functionsincluded in such modified spaces has a lower degree.For the sake of presentation we start from the case 𝑝 =
3, while we refer the reader to [36] for the case of 𝑝 = 𝑉 ℎ, ( P ) = (cid:110) 𝑣 ℎ ∈ 𝐻 ( P ) : Δ 𝑣 ℎ = , 𝑣 ℎ ∈ P ( 𝑒 ) , 𝜕 𝑛 𝑣 ℎ ∈ P ( 𝑒 ) , 𝜕 𝑛𝑛 𝑣 ℎ ∈ P ( 𝑒 ) ∀ 𝑒 ∈ 𝜕 P (cid:111) with associated degrees of freedom: (D1’) ℎ | 𝜈 | v 𝐷 𝜈 𝑣 ℎ ( v ) , | 𝜈 | ≤ v of 𝜕 P ; (D2’) ℎ 𝑒 ∫ 𝑒 𝜕 𝑛𝑛 𝑣 ℎ 𝑑𝑠 for any edge 𝑒 of 𝜕 P .In Ref. [9], we proved that the degrees of freedom (D1’) and (D2’) are unisolvent in ˜ 𝑉 ℎ, ( P ) and this spacecontains the linear subspace of polynomials of degree up to 4. Moreover, the associated global space obtainedby gluing together all the elemental spaces ˜ 𝑉 ℎ, ( P ) reads as:˜ 𝑉 ℎ, = (cid:110) 𝑣 ℎ ∈ 𝐻 ( Ω ) : 𝑣 ℎ | P ∈ ˜ 𝑉 ℎ, ( P ) ∀ P ∈ Ω ℎ (cid:111) , (13)is made of 𝐻 ( Ω ) functions.Analogously, in the general case we can build the modified lowest order spaces containing the space ofpolynomials of degree up to 2 𝑝 − 𝑉 𝑝ℎ, 𝑝 − ( P ) = (cid:110) 𝑣 ℎ ∈ 𝐻 𝑝 ( P ) : Δ 𝑝 𝑣 ℎ = , 𝑣 ℎ ∈ P 𝑝 − ( 𝑒 ) , 𝜕 𝑖𝑛 𝑣 ℎ ∈ P 𝑝 − − 𝑖 ( 𝑒 ) ,𝑖 = , . . . , 𝑝 − ∀ 𝑒 ∈ 𝜕 P (cid:111) , with associated degrees of freedom: (D1’) ℎ | 𝜈 | v 𝐷 𝜈 𝑣 ℎ ( v ) , | 𝜈 | ≤ 𝑝 − v of 𝜕 P ; (D2’) ℎ − + 𝑗𝑒 ∫ 𝑒 𝑞𝜕 𝑖𝑛 𝑣 ℎ 𝑑𝑠 for any 𝑞 ∈ P 𝑗 − ( 𝑒 ) and edge 𝑒 of 𝜕 P , 𝑗 = , . . . , 𝑝 − .2.3. Discrete bilinear form To define the elliptic projection Π ∇ , P 𝑟 : 𝑉 𝑝ℎ,𝑟 ( P ) → P 𝑟 ( P ) , we first need to introduce the vertex averageprojector (cid:98) Π P : 𝑉 𝑝ℎ,𝑟 ( P ) → P ( P ) , which projects any smooth enough function defined on P onto the spaceof constant polynomials. To this end, consider the continuous function 𝜓 defined on P . The vertex averageprojection of 𝜓 onto the constant polynomial space is given by: (cid:98) Π P 𝜓 = 𝑁 P ∑︁ v ∈ 𝜕 P 𝜓 ( v ) . (14)Finally, we define the elliptic projection Π ∇ , P 𝑟 : 𝑉 𝑝ℎ,𝑟 ( P ) → P 𝑟 ( P ) as the solution of the following finitedimensional variational problem 𝑎 P ( Π ∇ , P 𝑟 𝑣 ℎ , 𝑞 ) = 𝑎 P ( 𝑣 ℎ , 𝑞 ) ∀ 𝑞 ∈ P 𝑟 ( P ) , (15) (cid:98) Π P 𝐷 𝜈 Π ∇ , P 𝑟 𝑣 ℎ = (cid:98) Π P 𝐷 𝜈 𝑣 ℎ | 𝜈 | ≤ 𝑝 − . (16)According to Reference [9], such operator has two important properties: ( 𝑖 ) it is a polynomial-preserving operator in the sense that Π ∇ , P 𝑟 𝑞 = 𝑞 for every 𝑞 ∈ P 𝑟 ( P ) ; ( 𝑖𝑖 ) Π ∇ , P 𝑟 𝑣 ℎ is computable using only the degrees of freedom of 𝑣 ℎ .We write the symmetric bilinear form 𝑎 ℎ : 𝑉 𝑝ℎ,𝑟 × 𝑉 𝑝ℎ,𝑟 → R as the sum of local terms 𝑎 ℎ ( 𝑢 ℎ , 𝑣 ℎ ) = ∑︁ P ∈ Ω ℎ 𝑎 ℎ, P ( 𝑢 ℎ , 𝑣 ℎ ) , (17)where each local term 𝑎 ℎ, P : 𝑉 𝑝ℎ,𝑟 ( P ) × 𝑉 𝑝ℎ,𝑟 ( P ) → R is a symmetric bilinear form. We set 𝑎 ℎ, P ( 𝑢 ℎ , 𝑣 ℎ ) = 𝑎 P ( Π ∇ , P 𝑟 𝑢 ℎ , Π ∇ , P 𝑟 𝑣 ℎ ) + 𝑆 P ( 𝑢 ℎ − Π ∇ , P 𝑟 𝑢 ℎ , 𝑣 ℎ − Π ∇ , P 𝑟 𝑣 ℎ ) , (18)where 𝑆 P : 𝑉 𝑝ℎ,𝑟 ( P ) × 𝑉 𝑝ℎ,𝑟 ( P ) → R is a symmetric positive definite bilinear form such that 𝜎 ∗ 𝑎 P ( 𝑣 ℎ , 𝑣 ℎ ) ≤ 𝑆 P ( 𝑣 ℎ , 𝑣 ℎ ) ≤ 𝜎 ∗ 𝑎 P ( 𝑣 ℎ , 𝑣 ℎ ) ∀ 𝑣 ℎ ∈ 𝑉 𝑝ℎ,𝑟 ( P ) with Π ∇ , P 𝑟 𝑣 ℎ = , (19)for two some positive constants 𝜎 ∗ , 𝜎 ∗ independent of ℎ and P . The bilinear form 𝑎 ℎ, P (· , ·) has the twofundamental properties of 𝑟 - consistency and stability [9]: ( 𝑖 ) 𝑟 - Consistency : for every polynomial 𝑞 ∈ P 𝑟 ( P ) and function 𝑉 𝑝ℎ,𝑟 ( P ) we have: 𝑎 ℎ, P ( 𝑣 ℎ , 𝑞 ) = 𝑎 P ( 𝑣 ℎ , 𝑞 ) ; (20) ( 𝑖𝑖 ) Stability : there exist two positive constants 𝛼 ∗ , 𝛼 ∗ independent of ℎ and P such that for every 𝑣 ℎ ∈ 𝑉 𝑝ℎ,𝑟 ( P ) it holds: 𝛼 ∗ 𝑎 P ( 𝑣 ℎ , 𝑣 ℎ ) ≤ 𝑎 ℎ, P ( 𝑣 ℎ , 𝑣 ℎ ) ≤ 𝛼 ∗ 𝑎 P ( 𝑣 ℎ , 𝑣 ℎ ) . (21) We denote by 𝑓 ℎ the piecewise polynomial approximation of 𝑓 on Ω ℎ given by 𝑓 ℎ | P = Π , P 𝑟 − 𝑝 𝑓 , (22)9or 𝑟 ≥ 𝑝 − P ∈ Ω ℎ . Then, we set (cid:10) 𝑓 ℎ , 𝑣 ℎ (cid:11) = ∑︁ P ∈ Ω ℎ ∫ P 𝑓 ℎ 𝑣 ℎ 𝑑 x (23)which implies, using the 𝐿 -orthogonal projection, that (cid:10) 𝑓 ℎ , 𝑣 ℎ (cid:11) = ∑︁ P ∈ Ω ℎ ∫ P Π , P 𝑟 − 𝑝 𝑓 Π , P 𝑟 − 𝑝 𝑣 ℎ 𝑑 x = ∑︁ P ∈ Ω ℎ ∫ P 𝑓 Π , P 𝑟 − 𝑝 𝑣 ℎ 𝑑 x . (24)The right-hand side of (24) is computable by a combined use of the degrees of freedom (D1) - (D4) and theenhanced approach of Reference [2]. In this section we briefly sketch the construction of global virtual element spaces with arbitrary high orderof continuity. More precisely, we consider the local virtual element space defined as before, for 𝑟 ≥ 𝑝 − 𝑉 𝑝ℎ,𝑟 ( P ) = (cid:110) 𝑣 ℎ ∈ 𝐻 𝑝 ( P ) : Δ 𝑝 𝑣 ℎ ∈ P 𝑟 − 𝑝 ( P ) , 𝑣 ℎ ∈ P 𝑟 ( 𝑒 ) , 𝜕 𝑗𝑛 𝑣 ℎ ∈ P 𝑟 − 𝑗 ( 𝑒 ) ,𝑗 = , . . . , 𝑝 − ∀ 𝑒 ∈ 𝜕 P (cid:9) . Differently from the previous section, we make the degrees of freedom depend on a given parameter 𝑡 with 0 ≤ 𝑡 ≤ 𝑝 −
1. For a given value of 𝑡 we choose the degrees of freedom as follows (D1) ℎ | 𝜈 | v 𝐷 𝜈 𝑣 ℎ ( v ) , | 𝜈 | ≤ 𝑝 − v of P ; (D2) ℎ − 𝑒 ∫ 𝑒 𝑣 ℎ 𝑞 𝑑𝑠 for any 𝑞 ∈ P 𝑟 − 𝑝 ( 𝑒 ) , for any edge 𝑒 of 𝜕 P ; (D3) ℎ − + 𝑗𝑒 ∫ 𝑒 𝜕 𝑗𝑛 𝑣 ℎ 𝑞 𝑑𝑠 for any 𝑞 ∈ P 𝑟 − 𝑝 + 𝑗 ( 𝑒 ) and edge 𝑒 ∈ 𝜕 P , 𝑗 = , . . . , 𝑝 − (D4’) ℎ − P ∫ P 𝑞𝑣 ℎ 𝑑 x for any 𝑞 ∈ P 𝑟 − ( 𝑝 − 𝑡 ) ( P ) ;where as usual we assume P − 𝑛 (·) = { } for 𝑛 = , , , . . . .This set of degrees is still unisolvent, cf. [9]. Moreover, for 𝑟 ≥ 𝑝 − P 𝑟 ( P ) ⊂ 𝑉 𝑝ℎ,𝑟 ( P ) .Finally, it is worth noting that the choice (D4’) , if compared with (D4) , still guarantees that the associatedglobal space is made of 𝐶 𝑝 − (D1) - (D4’) to solve a differential probleminvolving the Δ 𝑝 − 𝑡 operator and 𝐶 𝑝 − ( Ω ) basis functions. For the sake of exposition, let us consider thefollowing two examples, in the context of the Laplacian and the Bilaplacian problem.1. Choosing 𝑝 and 𝑡 such that 𝑝 − 𝑡 = 𝐶 𝑝 − -conforming virtual element method for thesolution of the Laplacian problem. For example, for 𝑝 = 𝑡 = 𝑟 =
5, the local space 𝑉 ℎ, ( P ) endowed with the corresponding degrees of freedom (D1) - (D4’) can be employed to build a globalspace made of 𝐶 functions. It is also worth mentioning that the new choice (D4’) , differently from10he original choice (D4) , is essential for the computability of the elliptic projection, see (15)-(16), withrespect to the bilinear form 𝑎 P (· , ·) = ∫ P ∇(·)∇(·) 𝑑 x .2. Choosing 𝑝 and 𝑡 such that 𝑝 − 𝑡 = 𝐶 𝑝 − -conforming virtual element method for the solutionof the Bilaplacian problem. For example, for 𝑝 = 𝑡 = 𝑟 =
5, similarly to the previous case, thespace 𝑉 ℎ, ( P ) together with (D1) - (D4’) provides a global space of 𝐶 functions that can be employedfor the solution of the biharmonic problem.It is worth remembering that 𝐶 -regular virtual element basis function has been employed, e.g., in [25]to study residual based a posteriori error estimators for the virtual element approximation of second orderelliptic problems. Moreover, the solution of coupled elliptic problems of different order can take advantagefrom this flexibility of the degree of continuity of the basis functions. Indeed, for the sake of clarity considerthe conforming virtual element approximation of the following simplified situation: − Δ 𝑢 = 𝑓 in Ω , Δ 𝑢 = 𝑓 in Ω ,𝑢 = 𝑢 on Γ = Ω ∩ Ω ,𝜕 𝑛 𝑢 = 𝜕 𝑛 𝑢 on Γ ,𝑢 = 𝜕 Ω \ Γ ,𝑢 = 𝜕 Ω \ Γ ,𝜕 𝑛 𝑢 = 𝜕 Ω \ Γ . Handling the coupling conditions on Γ asks for the use of 𝐶 -regular virtual basis functions not only in Ω where the bilaplacian problem is defined, but also in Ω , where the second order elliptic problem is defined.Indeed, a simple use of 𝐶 -basis functions in Ω , which would be natural given the second order of theproblem, would not allow the imposition (or at least a simple imposition) of the gluing condition on thenormal derivatives. The following convergence result in the energy norm holds (see [9] for the proof).
Theorem 2.1.
Let 𝑓 ∈ 𝐻 𝑟 − 𝑝 + ( Ω ) be the forcing term at the right-hand side, 𝑢 the solution of the variationalproblem (6) and 𝑢 ℎ ∈ 𝑉 𝑝ℎ,𝑟 the solution of the virtual element method (11) . Then, it holds that || 𝑢 − 𝑢 ℎ || v ≤ 𝐶ℎ 𝑟 −( 𝑝 − ) (cid:0) | 𝑢 | 𝑟 + + | 𝑓 | 𝑟 − 𝑝 + (cid:1) . (25)Moreover, the following convergence results in lower order norms can established [9]. Theorem 2.2 (Even 𝑝 , even norms). Let 𝑓 ∈ 𝐻 𝑟 − 𝑝 + ( Ω ) , 𝑢 the solution of the variational problem (6) with 𝑝 = ℓ and 𝑣 ℎ ∈ 𝑉 𝑝ℎ,𝑟 the solution of the virtual element method (11) . Then, there exists a positive constant 𝐶 independent of ℎ such that | 𝑢 − 𝑢 ℎ | 𝑖 ≤ 𝐶ℎ 𝑟 + − 𝑖 (cid:16) | 𝑢 | 𝑟 + + | 𝑓 | 𝑟 −( 𝑝 − ) (cid:17) , (26) for every integer 𝑖 = , . . . , ℓ − . heorem 2.3 (Even 𝑝 , odd norms). Let 𝑓 ∈ 𝐻 𝑟 − 𝑝 + ( Ω ) , and 𝑢 the solution of the variational problem (6) with 𝑝 = ℓ and 𝑢 ℎ ∈ 𝑉 𝑝ℎ,𝑟 the solution of the virtual element method (11) . Then, there exists a positiveconstant 𝐶 independent of ℎ such that | 𝑢 − 𝑢 ℎ | 𝑖 + ≤ 𝐶ℎ ( 𝑟 + )−( 𝑖 + ) (cid:16) | 𝑢 | 𝑟 + + | 𝑓 | 𝑟 −( 𝑝 − ) (cid:17) , (27) for every integer 𝑖 = , . . . , ℓ − . Theorem 2.4 (Odd 𝑝 , even norms). Let 𝑢 be the solution of the variational problem (6) and 𝑢 ℎ ∈ 𝑉 𝑝ℎ,𝑟 thesolution of the virtual element method (11) . Then, there exists a positive constant 𝐶 independent of ℎ suchthat | 𝑢 − 𝑢 ℎ | 𝑖 ≤ 𝐶ℎ ( 𝑟 + )− 𝑖 (cid:16) | 𝑢 | 𝑟 + + | 𝑓 | 𝑟 −( 𝑝 − ) (cid:17) , (28) for every integer 𝑖 = , . . . , ℓ − . Theorem 2.5 (Odd 𝑝 , odd norms). Let 𝑢 be the solution of the variational problem (6) and 𝑢 ℎ ∈ 𝑉 𝑝ℎ,𝑟 thesolution of the virtual element method (11) . Then, there exists a positive constant 𝐶 independent of ℎ suchthat | 𝑢 − 𝑢 ℎ | 𝑖 + ≤ 𝐶ℎ ( 𝑟 + )−( 𝑖 + ) (cid:16) | 𝑢 | 𝑟 + + | 𝑓 | 𝑟 −( 𝑝 − ) (cid:17) , (29) for every integer 𝑖 = , . . . , ℓ − .
3. The virtual element method for the Cahn-Hilliard problem
Let Ω ⊂ R be an open, bounded domain with polygonal boundary Γ , 𝜓 : R → R with 𝜓 ( 𝑥 ) = ( − 𝑥 ) / 𝜙 ( 𝑥 ) = 𝜓 (cid:48) ( 𝑥 ) . We consider the Cahn-Hilliard problem: Find 𝑢 ( 𝑥, 𝑡 ) : Ω × [ , 𝑇 ] → R such that: (cid:164) 𝑢 − Δ (cid:0) 𝜙 ( 𝑢 ) − 𝛾 Δ 𝑢 (cid:1) = Ω × [ , 𝑇 ] , (30a) 𝑢 (· , ) = 𝑢 (·) in Ω , (30b) 𝜕 𝑛 𝑢 = 𝜕 𝑛 (cid:0) 𝜙 ( 𝑢 ) − 𝛾 Δ 𝑢 (cid:1) = 𝜕 Ω × [ , 𝑇 ] , (30c)where 𝜕 𝑛 denotes the (outward) normal derivative and 𝛾 ∈ R + , 0 < 𝛾 (cid:28)
1, represents the interface parameter.On the domain boundary we impose a no flux-type condition on 𝑢 and the chemical potential 𝜙 ( 𝑢 ) − 𝛾 Δ 𝑢 .To define the variational formulation of problem (30a)-(30c) we introduce the three bilinear forms: 𝑎 Δ ( 𝑣, 𝑤 ) = ∫ Ω (∇ 𝑣 ) : (∇ 𝑤 ) 𝑑 x ∀ 𝑣, 𝑤 ∈ 𝐻 ( Ω ) ,𝑎 ∇ ( 𝑣, 𝑤 ) = ∫ Ω ∇ 𝑣 · ∇ 𝑤 𝑑 x ∀ 𝑣, 𝑤 ∈ 𝐻 ( Ω ) ,𝑎 ( 𝑣, 𝑤 ) = ∫ Ω 𝑣 𝑤 𝑑 x ∀ 𝑣, 𝑤 ∈ 𝐿 ( Ω ) , ∇ being the Hessian operator) and the semi-linear form 𝑟 ( 𝑧 ; 𝑣, 𝑤 ) = ∫ Ω 𝜙 (cid:48) ( 𝑧 )∇ 𝑣 · ∇ 𝑤 𝑑 x ∀ 𝑧, 𝑣, 𝑤 ∈ 𝐻 ( Ω ) . Finally, introducing the functional space 𝑉 = (cid:8) 𝑣 ∈ 𝐻 ( Ω ) : 𝜕 𝑛 𝑣 = Γ (cid:9) , (31)which is a subspace of 𝐻 ( Ω ) .The weak formulation of problem (30a)-(30c) reads as: Find 𝑢 (· , 𝑡 ) ∈ 𝑉 such that 𝑎 ( (cid:164) 𝑢, 𝑣 ) + 𝛾 𝑎 Δ ( 𝑢, 𝑣 ) + 𝑟 ( 𝑢 ; 𝑢, 𝑣 ) = ∀ 𝑣 ∈ 𝑉, (32a) 𝑢 (· , ) = 𝑢 . (32b) In this section, we introduce the main building blocks for the conforming virtual discretization of the Cahn-Hilliard equation, report a convergence result and collect some numerical results assessing the theoreticalproperties of the proposed scheme. 𝐶 Virtual Element space
We briefly recall the construction of the virtual element space 𝑊 ℎ ⊂ 𝐻 ( Ω ) that we use to discretize(32a)-(32b); see [5] for more details.Given an element P ∈ Ω ℎ , the augmented local space (cid:101) 𝑉 ℎ | P is defined by (cid:101) 𝑉 ℎ | P = (cid:110) 𝑣 ∈ 𝐻 ( P ) : Δ 𝑣 ∈ P ( P ) , 𝑣 | 𝜕 P ∈ 𝐶 ( 𝜕 P ) , 𝑣 | 𝑒 ∈ P ( 𝑒 ) ∀ 𝑒 ∈ 𝜕 P , ∇ 𝑣 | 𝜕 P ∈ (cid:2) 𝐶 ( 𝜕 P ) (cid:3) , 𝜕 𝑛 𝑣 | 𝑒 ∈ P ( 𝑒 ) ∀ 𝑒 ∈ 𝜕 P (cid:111) , (33)with 𝜕 𝑛 denoting the (outward) normal derivative. Remark 3.1.
The space (cid:101) 𝑉 ℎ | P corresponds to the space ˜ 𝑉 𝑝ℎ, 𝑝 − ( P ) with 𝑝 = introduced in Section 2.2.2. We consider the two sets of linear operators from (cid:101) 𝑉 ℎ | P into R denoted by (D1) and (D2) and defined asfollows: (D1) contains linear operators evaluating 𝑣 ℎ at the 𝑛 = 𝑛 ( P ) vertices of P ; (D2) contains linear operators evaluating ∇ 𝑣 ℎ at the 𝑛 = 𝑛 ( P ) vertices of P .The output values of the two sets of operators (D1) and (D2) are sufficient to uniquely determine 𝑣 ℎ and ∇ 𝑣 ℎ on the boundary of P (cf. Section 2.2.2).We use of the following local bilinear forms for all P ∈ Ω ℎ 𝑎 Δ P ( 𝑣, 𝑤 ) = ∫ P (∇ 𝑣 ) : (∇ 𝑤 ) 𝑑 x ∀ 𝑣, 𝑤 ∈ 𝐻 ( P ) , (34) 𝑎 ∇ P ( 𝑣, 𝑤 ) = ∫ P ∇ 𝑣 · ∇ 𝑤 𝑑 x ∀ 𝑣, 𝑤 ∈ 𝐻 ( P ) , (35) 𝑎 P ( 𝑣, 𝑤 ) = ∫ P 𝑣 𝑤 𝑑 x ∀ 𝑣, 𝑤 ∈ 𝐿 ( P ) . (36)13ow, we introduce the elliptic projection operator Π Δ , P : (cid:101) 𝑉 ℎ | P → P ( P ) defined by 𝑎 Δ P ( Π Δ , P 𝑣 ℎ , 𝑞 ) = 𝑎 Δ P ( 𝑣 ℎ , 𝑞 ) ∀ 𝑞 ∈ P ( P ) , (37) (( Π Δ , P 𝑣 ℎ , 𝑞 )) P = (( 𝑣 ℎ , 𝑞 )) P ∀ 𝑞 ∈ P ( P ) , (38)for all 𝑣 ℎ ∈ (cid:101) 𝑉 ℎ | P where ((· , ·)) P is the Euclidean scalar product acting on the vectors that collect the vertexfunction values, i.e. (( 𝑣 ℎ , 𝑤 ℎ )) P = ∑︁ v ∈V P 𝑣 ℎ ( v ) 𝑤 ℎ ( v ) ∀ 𝑣 ℎ , 𝑤 ℎ ∈ 𝐶 ( P ) . As shown in [5], the operator Π Δ , P : (cid:101) 𝑉 ℎ | P → P ( P ) is well defined and uniquely determined on the basis ofthe information carried by the linear operators in (D1) and (D2) .Hinging upon the augmented space (cid:101) 𝑉 ℎ | P and employing the projector Π Δ , P we define our virtual localspace 𝑊 ℎ | P = (cid:8) 𝑣 ∈ (cid:101) 𝑉 ℎ | P : 𝑎 P ( Π Δ , P ( 𝑣 ) , 𝑞 ) = 𝑎 P ( 𝑣, 𝑞 ) ∀ 𝑞 ∈ P ( P ) (cid:9) . (39)Since 𝑊 ℎ | P ⊂ (cid:101) 𝑉 ℎ | P , operator Π Δ , P is well defined on 𝑊 ℎ | P and computable by using the values provided by (D1) and (D2) . Moreover, the set of operators (D1) and (D2) constitutes a set of degrees of freedom for thespace 𝑊 ℎ | P . Finally, there holds P ( P ) ⊆ 𝑊 ℎ | P .We now introduce two further projectors on the local space 𝑊 ℎ | P , namely Π , P and Π ∇ , P , that will beemployed together with the above projector Π Δ , P to build the discrete counterparts of the bilinear formsin (34). Operator Π , P : 𝑊 ℎ | P → P ( P ) is the standard 𝐿 projector on the space of quadratic polynomials in P . This is computable by means of the values of the degrees of freedom (D1) and (D2) (cf. [5]). To define Π ∇ , P : 𝑊 ℎ | P → P ( P ) we need the additional bilinear form 𝑎 ∇ (· , ·) : 𝑊 ℎ | P × 𝑊 ℎ | P → R that is given by 𝑎 ∇ ( 𝑣, 𝑤 ) = ∫ Ω ∇ 𝑣 · ∇ 𝑤 𝑑 x ∀ 𝑣, 𝑤 ∈ 𝐻 ( Ω ) . Operator Π ∇ , P is the elliptic projection defined with respect to 𝑎 ∇ (· , ·) : 𝑎 ∇ P ( Π ∇ , P 𝑣 ℎ , 𝑞 ) = 𝑎 ∇ P ( 𝑣 ℎ , 𝑞 ) ∀ 𝑞 ∈ P ( P ) , (40a) ∫ P Π ∇ , P 𝑣 ℎ 𝑑 x = ∫ P 𝑣 ℎ 𝑑 x . (40b)Such operator is well defined and uniquely determined by the values of (D1) and (D2) [5].We are now ready to introduce the global virtual element space, which defined as follows 𝑊 ℎ = (cid:8) 𝑣 ∈ 𝑉 : 𝑣 | P ∈ 𝑊 ℎ | P ∀ P ∈ Ω ℎ (cid:9) . The virtual element functions in 𝑊 ℎ and their gradients are continuous fields on Ω , so this functional spaceis a conforming subspace of 𝐻 ( Ω ) . The global degrees of freedom of 𝑊 ℎ are obtained by collecting theelemental degrees of freedom, so the dimension of 𝑊 ℎ is three times the number of the mesh vertices, andevery virtual element function 𝑣 ℎ defined on Ω is uniquely determined by ( 𝑖 ) its values at the mesh vertices; ( 𝑖𝑖 ) its gradient values at the mesh vertices. 14inally, we recommended to scale the degrees of freedom (D2) by some local characteristic mesh size ℎ v inorder to obtain a better condition number of the final system. We start by introducing the discrete versions of the elemental bilinear form forms in (34). Let P ∈ Ω ℎ bea generic mesh element and 𝑠 P (· , ·) : 𝑊 ℎ | P × 𝑊 ℎ | P → R the positive definite bilinear form given by: 𝑠 P ( 𝑣 ℎ , 𝑤 ℎ ) = ∑︁ v ∈V P (cid:16) 𝑣 ℎ ( v ) 𝑤 ℎ ( 𝜈 ) + ℎ v ∇ 𝑣 ℎ ( v ) · ∇ 𝑤 ℎ ( v ) (cid:17) ∀ 𝑣 ℎ , 𝑤 ℎ ∈ 𝑊 ℎ | P , where ℎ v is a characteristic mesh size length associated with node v , e.g., the maximum diameter among theelements having v as a vertex.Recalling (34), we consider the virtual element bilinear forms: 𝑎 Δ ℎ, P ( 𝑣 ℎ , 𝑤 ℎ ) = 𝑎 Δ P ( Π Δ , P 𝑣 ℎ , Π Δ , P 𝑤 ℎ ) + ℎ − P 𝑠 P (cid:0) 𝑣 ℎ − Π Δ , P 𝑣 ℎ , 𝑤 ℎ − Π Δ , P 𝑤 ℎ (cid:1) , (41) 𝑎 ∇ ℎ, P ( 𝑣 ℎ , 𝑤 ℎ ) = 𝑎 ∇ P ( Π ∇ , P 𝑣 ℎ , Π ∇ , P 𝑤 ℎ ) + 𝑠 P (cid:0) 𝑣 ℎ − Π ∇ , P 𝑣 ℎ , 𝑤 ℎ − Π ∇ , P 𝑤 ℎ ) , (42) 𝑎 ℎ, P ( 𝑣 ℎ , 𝑤 ℎ ) = 𝑎 P ( Π , P 𝑣 ℎ , Π , P 𝑤 ℎ ) + ℎ P 𝑠 P (cid:0) 𝑣 ℎ − Π , P 𝑣 ℎ , 𝑤 ℎ − Π , P 𝑤 ℎ (cid:1) (43)for all 𝑣 ℎ , 𝑤 ℎ ∈ 𝑊 ℎ | P . Under the mesh regularity conditions of Section 1.2, we can prove the consistency andstability of the discrete bilinear forms. Let the symbol † stands for “ Δ ”, “ ∇ ” or “0”. We have: (A) (polynomial consistency) 𝑎 † ℎ, P ( 𝑝, 𝑣 ℎ ) = 𝑎 † P ( 𝑝, 𝑣 ℎ ) ∀ 𝑝 ∈ P ( P ) , 𝑣 ℎ ∈ 𝑊 ℎ | P ; (B) (stability) there exist two positive constants 𝑐 ∗ and 𝑐 ∗ independent of ℎ and the element P ∈ Ω ℎ such that 𝑐 ∗ 𝑎 † P ( 𝑣 ℎ , 𝑣 ℎ ) ≤ 𝑎 † ℎ, P ( 𝑣 ℎ , 𝑣 ℎ ) ≤ 𝑐 ∗ 𝑎 † P ( 𝑣 ℎ , 𝑣 ℎ ) ∀ 𝑣 ℎ ∈ 𝑊 ℎ | P . A consequence of the above properties is that the bilinear form 𝑎 † ℎ, P (· , ·) is continuous with respect to therelevant norm, which is 𝐻 for (41), 𝐻 for (42), and 𝐿 for (43). For every choice of † , the correspondingglobal bilinear form is 𝑎 † ℎ ( 𝑣 ℎ , 𝑤 ℎ ) = ∑︁ P ∈ Ω ℎ 𝑎 † ℎ, P ( 𝑣 ℎ , 𝑤 ℎ ) ∀ 𝑣 ℎ , 𝑤 ℎ ∈ 𝑊 ℎ . (30a)-(30c)We now turn our attention to the semilinear form 𝑟 (· ; · , ·) , which we can also write as the sum of elementalcontributions: 𝑟 ( 𝑧 ; 𝑣, 𝑤 ) = ∑︁ P ∈ Ω ℎ 𝑟 P ( 𝑧 ; 𝑣, 𝑤 ) ∀ 𝑧, 𝑣𝑤 ∈ 𝐻 ( Ω ) where 𝑟 P ( 𝑧 ; 𝑣, 𝑤 ) = ∫ P ( 𝑧 − )∇ 𝑣 · ∇ 𝑤𝑑 x ∀ P ∈ Ω ℎ . On each element P , we approximate the term 𝑧 ( 𝑥 ) by means of its cell average, which we compute using15he 𝐿 ( P ) bilinear form 𝑎 ℎ, P (· , ·) : 𝑧 ℎ | P ≈ | P | − 𝑎 ℎ, P ( 𝑧 ℎ , 𝑧 ℎ ) , where we recall that | P | is the area of element P . This approach has the correct approximation properties andpreserves the positivity of 𝑧 .We therefore propose the following approximation of the local nonlinear forms 𝑟 ℎ, P ( 𝑧 ℎ ; 𝑣 ℎ , 𝑤 ℎ ) = (cid:155) 𝜙 (cid:48) ( 𝑧 ℎ ) | P 𝑎 ∇ ℎ, P ( 𝑣 ℎ , 𝑤 ℎ ) ∀ 𝑧 ℎ , 𝑣 ℎ , 𝑤 ℎ ∈ 𝑊 ℎ | P , where (cid:155) 𝜙 (cid:48) ( 𝑧 ℎ ) | P = | P | − 𝑎 ℎ, P ( 𝑧 ℎ , 𝑧 ℎ ) −
1. The global form is then assembled as 𝑟 ℎ ( 𝑧 ℎ ; 𝑣 ℎ , 𝑤 ℎ ) = ∑︁ P ∈ Ω ℎ 𝑟 ℎ, P ( 𝑧 ℎ ; 𝑣 ℎ , 𝑤 ℎ ) ∀ 𝑧 ℎ , 𝑣 ℎ 𝑤 ℎ ∈ 𝑊 ℎ . The virtual element discretization of problem (32a)(32b) follows a Galerkin approach in space combinedwith a backward Euler time-stepping scheme. Consider the functional space 𝑊 ℎ = 𝑊 ℎ ∩ 𝑉 = (cid:8) 𝑣 ∈ 𝑊 ℎ : 𝜕 𝑛 𝑣 = 𝜕 Ω (cid:9) , which includes the boundary conditions. Then, we introduce the the semi-discrete approximation: Find 𝑢 ℎ (· , 𝑡 ) in 𝑊 ℎ such that 𝑎 ℎ ( (cid:164) 𝑢 ℎ , 𝑣 ℎ ) + 𝛾 𝑎 Δ ℎ ( 𝑢 ℎ , 𝑣 ℎ ) + 𝑟 ℎ ( 𝑢 ℎ ; 𝑢 ℎ , 𝑣 ℎ ) = ∀ 𝑣 ℎ ∈ 𝑊 ℎ , (44) 𝑢 ℎ ( , ·) = 𝑢 ,ℎ (·) , (45)where 𝑢 ,ℎ is a suitable approximation of 𝑢 in 𝑊 ℎ and 𝑎 ℎ (· , ·) , 𝑎 Δ ℎ and 𝑟 ℎ are the virtual element bilinearforms defined in the previous section.To formulate the fully discrete scheme, we subdivide the time interval [ , 𝑇 ] into 𝑁 uniform sub-intervalsof length 𝑘 = 𝑇 / 𝑁 by means of the time nodes 0 = 𝑡 < 𝑡 < . . . < 𝑡 𝑁 − < 𝑡 𝑁 = 𝑇 , and denote the virtualelement approximation of the solution 𝑢 (· , 𝑡 ) at 𝑢 (· , 𝑡 𝑖 ) in 𝑊 ℎ by 𝑢 𝑖ℎ,𝑘 . The fully discrete problem reads as: Given 𝑢 ℎ𝑘 = 𝑢 ,ℎ ∈ 𝑊 ℎ , find 𝑢 𝑖ℎ𝑘 ∈ 𝑊 ℎ , 𝑖 = , . . . , 𝑁 such that 𝑘 − 𝑎 ℎ ( 𝑢 𝑖ℎ𝑘 − 𝑢 𝑖 − ℎ𝑘 , 𝑣 ℎ ) + 𝛾 𝑎 Δ ℎ ( 𝑢 𝑖ℎ𝑘 , 𝑣 ℎ ) + 𝑟 ℎ ( 𝑢 𝑖ℎ𝑘 , 𝑢 𝑖ℎ𝑘 ; 𝑣 ℎ ) = ∀ 𝑣 ℎ ∈ 𝑊 ℎ . (46)The semidiscrete Virtual Element formulation given in (44)-(45) converges to the exact solution of problem(32a)-(32b) according to the result stated in this theorem and proved in [5]. Theorem 3.2.
Let 𝑢 be the solution of problem (32a) - (32b) . Let 𝑢 ℎ be the virtual element approximationprovided by (44) - (45) and assume that || 𝑢 ℎ || 𝐿 ∞ ( Ω ) ≤ 𝐶 for all 𝑡 ∈ ( , 𝑇 ] and some positive constant 𝐶 independent of ℎ . Then, it holds that || 𝑢 − 𝑢 ℎ || 𝐿 ( Ω ) (cid:46) ℎ for every 𝑡 ∈ [ , 𝑇 ] .3.3. Numerical results In this test, taken from [5] we study the convergence of our VEM discretization applied to the Cahn-Hilliard problem with a load term 𝑓 obtained by enforcing as exact solution 𝑢 ( 𝑥, 𝑦, 𝑡 ) = 𝑡 cos ( 𝜋𝑥 ) cos ( 𝜋𝑦 ) .The parameter 𝛾 is set to 1 /
10 and the time step size Δ 𝑡 is 1 𝑒 −
7. The 𝐻 , 𝐻 and 𝐿 errors are computedat 𝑡 = . igure 2: Spinoidal decomposition on the unit square at three temporal frames for a Voronoi polygonal mesh of 4096 elements. method, using the 𝑙 norm of the relative residual as a stopping criterion. The tolerance for convergence is1 𝑒 − Table 1: 𝐻 , 𝐻 and 𝐿 errors and convergence rates 𝛼 computed on four quadrilateral meshes discretizing the unit square [5]. ℎ | 𝑢 − 𝑢 ℎ | 𝐻 ( Ω ) 𝛼 | 𝑢 − 𝑢 ℎ | 𝐻 ( Ω ) 𝛼 || 𝑢 − 𝑢 ℎ || 𝐿 ( Ω ) 𝛼 𝐿 norm as expected from Theorem 3.2. In the 𝐻 and 𝐻 seminorms, the method convergeswith order 1 and 2 respectively, as we can expect from the FEM theory and the approximation propertiesof the virtual element space. Finally, in Figure 2 we report the results of a spinoidal decomposition. Forcompleteness, we recall that spinoidal decomposition is a physical phenomenon consisting of the separationof a mixture of two or more components to bulk regions of each, which occurs when a high-temperaturemixture of different components is rapidly cooled. We employ an initial datum 𝑢 chosen to be a uniformlydistributed random perturbation between −
4. The virtual element method for the elastodynamics problem
We consider an elastic body occupying the open, bounded polygonal domain Ω ⊂ R with Lipschitzboundary Γ . We assume that boundary Γ can be split into the two disjoint subsets Γ 𝐷 and Γ 𝑁 , so that Γ = Γ 𝐷 ∪ Γ 𝑁 and with the one-dimensional Lebesgue measure (length) | Γ 𝐷 ∩ Γ 𝑁 | =
0. For the well-posedness of the mathematical model, we further require length of Γ 𝐷 is nonzero, i.e., | Γ 𝐷 | >
0. Let
𝑇 > f ∈ 𝐿 (cid:0) , 𝑇 ; [ 𝐿 ( Ω )] (cid:1) , the boundary function g 𝑁 ∈ 𝐶 (cid:0) , 𝑇 ; [ 𝐻 / , Γ 𝑁 ] (cid:1) , and the initial functions u ∈ [ 𝐻 , Γ 𝐷 ( Ω )] , u ∈ [ 𝐿 ( Ω )] . For such time-dependent vector fields, we may indicate the dependence on time explicitly, e.g., f ( 𝑡 ) : = f (· , 𝑡 ) ∈ [ 𝐿 ( Ω )] ,or drop it out to ease the notation when it is obvious from the context.17he equations governing the two-dimensional initial/boundary-value problem of linear elastodynamicsfor the displacement vector u : Ω × [ , 𝑇 ] → R are: 𝜌 (cid:165) u − ∇ · 𝝈 ( u ) = f in Ω × ( , 𝑇 ] , (47) u = on Γ 𝐷 × ( , 𝑇 ] , (48) 𝝈 ( u ) n = g 𝑁 on Γ 𝑁 × ( , 𝑇 ] , (49) u = u in Ω × { } , (50) (cid:164) u = u in Ω × { } . (51)Here, 𝜌 is the mass density, which we suppose to be a strictly positive and uniformly bounded functionand 𝝈 ( u ) is the stress tensor. In (48) we assume homogeneous Dirichlet boundary conditions on Γ 𝐷 . Thisassumption is made only to ease the exposition and the analysis, as our numerical method is easily extendableto nonhomogeneous Dirichlet boundary conditions.We denote the space of the symmetric, 2 × R × and assume that the stresstensor 𝝈 : Ω × [ , 𝑇 ] → R × is expressed, according to Hooke’s law, by 𝝈 ( u ) = D 𝜺 ( u ) , where, 𝜺 ( u ) denotesthe symmetric gradient of u , i.e., 𝜺 ( u ) = (cid:0) ∇ u + (∇ u ) 𝑇 (cid:1) /
2, and D = D ( x ) : R × −→ R × is the stiffness tensor D 𝝉 = 𝜇 𝝉 + 𝜆 tr ( 𝝉 ) I (52)for all 𝝉 ∈ R × . In this definition, I and tr (·) are the identity matrix and the trace operator; 𝜆 and 𝜇 are the firstand second Lamé coefficients, which we assume to be in 𝐿 ∞ ( Ω ) and nonnegative. The compressional (P) andshear (S) wave velocities of the medium are respectively obtained through the relations 𝑐 𝑃 = √︁ ( 𝜆 + 𝜇 )/ 𝜌 and 𝑐 𝑆 = √︁ 𝜇 / 𝜌 .Let V = (cid:2) 𝐻 Γ 𝐷 ( Ω ) (cid:3) be the space of 𝐻 vector-valued functions with null trace on Γ 𝐷 . We consider the twobilinear forms 𝑚 (· , ·) , 𝑎 (· , ·) : V × V → R defined as 𝑚 ( w , v ) = ∫ Ω 𝜌 w · v 𝑑 x ∀ w , v ∈ V , (53) 𝑎 ( w , v ) = ∫ Ω 𝝈 ( w ) : 𝜺 ( v ) 𝑑 x ∀ w , v ∈ V , (54)and the linear functional 𝐹 (·) : V → R defined as 𝐹 ( v ) = ∫ Ω f · v 𝑑 x + ∫ Γ 𝑁 g 𝑁 · v 𝑑𝑠 ∀ v ∈ V . (55)The variational formulation of the linear elastodynamics equations reads as: For all 𝑡 ∈ ( , 𝑇 ] find u ( 𝑡 ) ∈ V such that for 𝑡 = it holds that u ( ) = u and (cid:164) u ( ) = u and 𝑚 ( (cid:165) u , v ) + 𝑎 ( u , v ) = 𝐹 ( v ) ∀ v ∈ V . (56)As shown, for example, by Raviart and Thomas (see Theorem 8-3.1 [85]) the variational problem (56) is wellposed and its unique solution satisfies u ∈ 𝐶 (cid:0) , 𝑇 ; V (cid:1) ∩ 𝐶 (cid:0) , 𝑇 ; [ 𝐿 ( Ω )] (cid:1) . In this section we introduce the main building blocks for the conforming virtual element discretizationof the elastodynamics equation, report stability and convergence results and collect some numerical resultsassessing the theoretical properties of the proposed scheme.18 .2.1. Virtual element spaces
Let 𝑘 ≥ V ℎ𝑘 : = (cid:110) v ∈ V : v | P ∈ V ℎ𝑘 ( P ) for every P ∈ Ω ℎ (cid:111) (57)where V ℎ𝑘 ( P ) = (cid:2) 𝑉 ℎ𝑘 ( P ) (cid:3) , with 𝑉 ℎ𝑘 ( P ) : = (cid:110) 𝑣 ℎ ∈ 𝐻 ( P ) : 𝑣 ℎ | 𝜕 P ∈ 𝐶 ( 𝜕 P ) , 𝑣 ℎ | 𝑒 ∈ P 𝑘 ( 𝑒 ) ∀ 𝑒 ∈ 𝜕 P , Δ 𝑣 ℎ ∈ P 𝑘 ( P ) , (cid:0) 𝑣 ℎ − Π ∇ 𝑘 𝑣 ℎ , 𝜇 ℎ (cid:1) P = ∀ 𝜇 ℎ ∈ P 𝑘 ( P )\ P 𝑘 − ( P ) (cid:111) , (58)where Π ∇ 𝑘 : 𝐻 ( P ) ∩ 𝐶 ( P ) → P 𝑘 ( P ) is the usual elliptic projection of a function 𝑣 ℎ on the space ofpolynomials of degree 𝑘 , cf. (15)-(16).Each virtual element function 𝑣 ℎ ∈ 𝑉 ℎ𝑘 ( P ) is uniquely characterized by (C1) the values of 𝑣 ℎ at the vertices of P ; (C2) the moments of 𝑣 ℎ of order up to 𝑘 − 𝑒 ∈ 𝜕 P :1 | 𝑒 | ∫ 𝑒 𝑣 ℎ 𝑚 𝑑𝑠, ∀ 𝑚 ∈ M 𝑘 − ( 𝑒 ) , ∀ 𝑒 ∈ 𝜕 P ; (59) (C3) the moments of 𝑣 ℎ of order up to 𝑘 − P :1 | P | ∫ P 𝑣 ℎ 𝑚 𝑑 x , ∀ 𝑚 ∈ M 𝑘 − ( P ) . (60)As usual, the degrees of freedom of the global space V ℎ𝑘 are provided by collecting all the local degrees offreedom (which allow the computation of the elliptic projection Π ∇ 𝑘 ), and their unisolvence is an immediateconsequence of the unisolvence of the local degrees of freedom for the elemental spaces 𝑉 ℎ𝑘 ( P ) . In the virtual element setting, we define the bilinear forms 𝑚 ℎ (· , ·) and 𝑎 ℎ (· , ·) as the sum of elementalcontributions, which are respectively denoted by 𝑚 ℎ, P (· , ·) and 𝑎 ℎ, P (· , ·) : 𝑚 ℎ (· , ·) : V ℎ𝑘 × V ℎ𝑘 → R , with 𝑚 ℎ ( v ℎ , w ℎ ) = ∑︁ P ∈ Ω ℎ 𝑚 ℎ, P ( v ℎ , w ℎ ) ,𝑎 ℎ (· , ·) : V ℎ𝑘 × V ℎ𝑘 → R , with 𝑎 ℎ ( v ℎ , w ℎ ) = ∑︁ P ∈ Ω ℎ 𝑎 ℎ, P ( v ℎ , w ℎ ) . The local bilinear form 𝑚 ℎ, P (· , ·) is given by 𝑚 ℎ, P ( v ℎ , w ℎ ) = ∫ P 𝜌 Π 𝑘 v ℎ · Π 𝑘 w ℎ 𝑑𝑉 + 𝑆 P 𝑚 ( v ℎ , w ℎ ) , (61)where 𝑆 P 𝑚 (· , ·) is the local stabilization term. The bilinear form 𝑚 ℎ, P depends on the orthogonal projections Π 𝑘 v ℎ and Π 𝑘 w ℎ , which are computable from the degrees of freedom of v ℎ and w ℎ . The local form 𝑆 P 𝑚 (· , ·) : V ℎ𝑘 × V ℎ𝑘 → R can be any symmetric and coercive bilinear form that is computable from thedegrees of freedom and for which there exist two strictly positive real constants 𝜎 ∗ and 𝜎 ∗ such that 𝜎 ∗ 𝑚 P ( v ℎ , v ℎ ) ≤ 𝑆 P 𝑚 ( v ℎ , v ℎ ) ≤ 𝜎 ∗ 𝑚 P ( v ℎ , v ℎ ) v ℎ ∈ ker (cid:0) Π 𝑘 (cid:1) ∩ V ℎ𝑘 ( P ) . (62)Computable stabilizations 𝑆 P 𝑚 (· , ·) are provided by resorting to the two-dimensional stabilizations of theeffective choices for the scalar case proposed in the literature[75, 51]. The local bilinear form 𝑎 ℎ, P is given by 𝑎 ℎ, P ( v ℎ , w ℎ ) = ∫ P D Π 𝑘 − ( 𝜺 ( v ℎ )) : Π 𝑘 − ( 𝜺 ( w ℎ )) 𝑑𝑉 + 𝑆 P 𝑎 ( v ℎ , w ℎ ) , (63)19here 𝑆 P 𝑎 (· , ·) is the local stabilization term. The bilinear form 𝑎 ℎ, P depends on the orthogonal projections Π 𝑘 − ∇ v ℎ and Π 𝑘 − ∇ w ℎ , which are computable from the degrees of freedom of v ℎ and w ℎ . On its turn, 𝑆 P 𝑎 (· , ·) : V ℎ𝑘 × V ℎ𝑘 → R can be any symmetric and coercive bilinear form that is computable from thedegrees of freedom and for which there exist two strictly positive real constants 𝜎 ∗ and 𝜎 ∗ such that 𝜎 ∗ 𝑎 P ( v ℎ , v ℎ ) ≤ 𝑆 P 𝑚 ( v ℎ , v ℎ ) ≤ 𝜎 ∗ 𝑎 P ( v ℎ , v ℎ ) v ℎ ∈ ker (cid:0) Π 𝑘 (cid:1) ∩ V ℎ𝑘 ( P ) . (64)Moreover, the bilinear form 𝑆 P 𝑎 (· , ·) must scale with respect to ℎ like 𝑎 P (· , ·) , i.e., as O ( ) . As before, wecan define computable stabilizations 𝑆 P 𝑎 (· , ·) by resorting to the two-dimensional stabilizations for the scalarcase proposed in the literature [75, 51]. As usual, the discrete bilinear forms 𝑎 ℎ, P (· , ·) and 𝑚 ℎ, P (· , ·) satisfythe 𝑘 -consistency and stability properties. The stability constants may depend on physical parameters and thepolynomial degree 𝑘 [18, 7]. We approximate the right-hand side (67) of the variational formulation by means of the linear functional 𝐹 ℎ (·) : V ℎ𝑘 → R given by 𝐹 ℎ ( v ℎ ) = ∫ Ω f · Π 𝑘 − ( v ℎ ) 𝑑𝑉 + ∑︁ 𝑒 ∈ Γ 𝑁 ∫ 𝑒 g 𝑁 · v ℎ 𝑑𝑠 ∀ v ℎ ∈ V ℎ𝑘 . (65)The linear functional 𝐹 ℎ (·) is clearly computable since the edge trace v ℎ | 𝑒 is a known polynomial and Π 𝑘 ( v ℎ ) is computable from the degrees of freedom of v ℎ . Moreover, 𝐹 ℎ (·) is a bounded functional. In fact, when g 𝑁 = | 𝐹 ℎ ( v ℎ )| ≤ (cid:12)(cid:12)(cid:12)(cid:12)∫ Ω f ( 𝑡 ) · Π 𝑘 − ( v ℎ ) 𝑑𝑉 (cid:12)(cid:12)(cid:12)(cid:12) ≤ || f ( 𝑡 )|| (cid:12)(cid:12)(cid:12)(cid:12) Π 𝑘 − ( v ℎ ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ || f ( 𝑡 )|| || v ℎ || ∀ 𝑡 ∈ [ , 𝑇 ] . (66)This estimate is used in the proof of the stability of the semi-discrete virtual element approximation (seeTheorem 4.1). The semi-discrete virtual element approximation of (56) reads as:
For all 𝑡 ∈ ( , 𝑇 ] find u ℎ ( 𝑡 ) ∈ V ℎ𝑘 suchthat for 𝑡 = it holds that u ℎ ( ) = ( u ) 𝐼 and (cid:164) u ℎ ( ) = ( u ) 𝐼 and 𝑚 ℎ ( (cid:165) u ℎ , v ℎ ) + 𝑎 ℎ ( u ℎ , v ℎ ) = 𝐹 ℎ ( v ℎ ) ∀ v ℎ ∈ V ℎ𝑘 . (67)Here, u ℎ ( 𝑡 ) is the virtual element approximation of u and v ℎ is the generic test function in V ℎ𝑘 , while ( u ) 𝐼 and ( u ) 𝐼 are the virtual element interpolants of the initial solution functions u ( ) and (cid:164) u ( ) .We carry out the time integration by applying the leap-frog time marching scheme [84] to the secondderivative in time (cid:165) u ℎ . To this end, we subdivide the interval ( , 𝑇 ] into 𝑁 𝑇 subintervals of amplitude Δ 𝑡 = 𝑇 / 𝑁 𝑇 and at every time level 𝑡 𝑛 = 𝑛 Δ 𝑡 we consider the variational problem for 𝑛 ≥ 𝑚 ℎ ( u 𝑛 + ℎ , v ℎ ) − 𝑚 ℎ ( u 𝑛ℎ , v ℎ ) + 𝑚 ℎ ( u 𝑛 − ℎ , v ℎ ) + Δ 𝑡 𝑎 ℎ ( u 𝑛ℎ , v ℎ ) = Δ 𝑡 𝐹 𝑛ℎ ( v ℎ ) ∀ v ℎ ∈ V ℎ𝑘 , (68)and initial step 𝑚 ℎ ( u ℎ , v ℎ ) − 𝑚 ℎ ( u , v ℎ ) − Δ 𝑡𝑚 ℎ ( u , v ℎ ) + Δ 𝑡 𝑎 ℎ ( u , v ℎ ) = Δ 𝑡 𝐹 ℎ ( v ℎ ) ∀ v ℎ ∈ 𝑉 ℎ𝑘 . We employ the energy norm ||| v ℎ ( 𝑡 )||| = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) 𝜌 (cid:164) v ℎ ( 𝑡 ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + | v ℎ ( 𝑡 )| , 𝑡 ∈ [ , 𝑇 ] , (69)which is defined for all v ℎ ∈ V ℎ𝑘 . The local stability property of the bilinear forms 𝑚 ℎ (· , ·) and 𝑎 ℎ (· , ·) impliesthe equivalence relation 𝑚 ℎ ( (cid:164) v ℎ , (cid:164) v ℎ ) + 𝑎 ℎ ( v ℎ , v ℎ ) (cid:46) ||| v ℎ ( 𝑡 )||| (cid:46) 𝑚 ℎ ( (cid:164) v ℎ , (cid:164) v ℎ ) + 𝑎 ℎ ( v ℎ , v ℎ ) (70)for all time-dependent virtual element functions v ℎ ( 𝑡 ) with square integrable derivative (cid:164) v ℎ ( 𝑡 ) .The hidden constants in (70) are independent of the mesh size parameter ℎ [7]. However, they may dependon the stability parameters, the physical parameters and the polynomial degree 𝑘 [19]. It is worth noting thatthe dependence on 𝑘 does not seem to have a relevant impact on the optimality of the convergence rates in thenumerical experiments of Section 4.3. The following stability result has been proved in [7]. Theorem 4.1.
Let f ∈ 𝐿 (cid:0) , 𝑇 ; [ 𝐿 ( Ω )] (cid:1) and let u ℎ ∈ 𝐶 (cid:0) , 𝑇 ; V ℎ𝑘 (cid:1) be the solution of (67) . Then, it holds ||| u ℎ ( 𝑡 )||| (cid:46) |||( u ) 𝐼 ||| + ∫ 𝑡 || f ( 𝜏 )|| , Ω 𝑑𝜏. (71) The hidden constant in (cid:46) is independent of ℎ , but may depend on the model parameters and approximationconstants and the polynomial degree 𝑘 . We point out that in the case of f null external force, i.e. f = , the above bound reduces to ||| u ℎ ( 𝑡 )||| (cid:46) |||( u ) 𝐼 ||| that is the virtual element approximation is dissipative.Now, we recall [7] the convergence of the semi-discrete virtual element approximation in the energy norm(69). Theorem 4.2.
Let u ∈ 𝐶 (cid:0) , 𝑇 ; [ 𝐻 𝑚 + ( Ω )] (cid:1) , 𝑚 ∈ N , be the exact solution of problem (56) . Let u ℎ ∈ V ℎ𝑘 bethe solution of the semi-discrete problem (67) . For f ∈ 𝐿 (cid:0) ( , 𝑇 ) ; (cid:2) 𝐻 𝑚 − ( Ω ) (cid:3) (cid:1) we have that sup <𝑡 ≤ 𝑇 ||| u ( 𝑡 ) − u ℎ ( 𝑡 )||| (cid:46) ℎ 𝜇 𝑘 𝑚 sup <𝑡 ≤ 𝑇 (cid:16) || (cid:164) u ( 𝑡 )|| 𝑚 + + || u ( 𝑡 )|| 𝑚 + (cid:17) + ∫ 𝑇 (cid:18) ℎ 𝜇 + 𝑘 𝑚 (|| (cid:165) u ( 𝜏 )|| 𝑚 + + || (cid:164) u ( 𝜏 )|| 𝑚 + ) + ℎ 𝜇 𝑘 𝑚 (cid:0) || (cid:165) u ( 𝜏 )|| 𝑚 + + || (cid:164) u ( 𝜏 )|| 𝑚 + (cid:1)(cid:19) 𝑑𝜏 + ∫ 𝑇 ℎ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:0) 𝐼 − Π 𝑘 − (cid:1) f ( 𝜏 ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) 𝑑𝜏, (72) where 𝜇 = min ( 𝑘, 𝑚 ) . The hidden constant in “ (cid:46) “ is independent of ℎ , but may depend on the modelparameters and approximation constants, the polynomial degree 𝑘 , and the final observation time 𝑇 . Finally, we state the convergence result in the 𝐿 norm, whose proof is again found in [7]. Theorem 4.3.
Let u be the exact solution of problem (56) under the assumption that domain Ω is 𝐻 -regularand u ℎ ∈ V ℎ𝑘 the solution of the virtual element method stated in (67) . If u , (cid:164) u , (cid:165) u ∈ 𝐿 (cid:0) , 𝑇 ; (cid:2) 𝐻 𝑚 + ( Ω ) ∩ 𝐻 ( Ω ) (cid:3) (cid:1) , with integer 𝑚 ≥ , then the following estimate holds for almost every 𝑡 ∈ [ , 𝑇 ] by setting esh 1 Mesh 2 Mesh 3 Figure 3: Base meshes (top row) and first refined meshes (bottom row) of the following mesh families from left to right: randomizedquadrilateral mesh; mainly hexagonal mesh; nonconvex octagonal mesh [7]. 𝜇 = min ( 𝑚, 𝑘 ) : || u ( 𝑡 ) − u ℎ ( 𝑡 )|| (cid:46) || u ℎ ( ) − u || + || (cid:164) u ℎ ( ) − u || + ℎ 𝜇 + 𝑘 𝑚 + (cid:16) || (cid:165) u || 𝐿 ( ,𝑇 ; [ 𝐻 𝑚 + ( Ω ) ] ) + || (cid:164) u || 𝐿 ( ,𝑇 ; [ 𝐻 𝑚 + ( Ω ) ] ) + || u || 𝐿 ( ,𝑇 ; [ 𝐻 𝑚 + ( Ω ) ] ) (cid:17) + ∫ 𝑇 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:0) − Π 𝑘 − (cid:1) f ( 𝜏 ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) 𝑑𝜏. (73) The hidden constant in “ (cid:46) “ is independent of ℎ , but may depend on the model parameters and approximationconstants 𝜚 , 𝜇 ∗ , and the polynomial degree 𝑘 , and the final observation time 𝑇 .4.3. Numerical Results In this section, we report from [7] a set of numerical results assessing the convergence properties of thevirtual element discretization by using a manufactured solution on three different mesh families, each onepossessing some special feature.In particular, we let
Ω = ( , ) for 𝑡 ∈ [ , 𝑇 ] , 𝑇 =
1, and consider initial condition u , boundary condition g and forcing term f determined from the exact solution: u ( 𝑥, 𝑦, 𝑡 ) = cos (cid:18) 𝜋 𝑡𝑇 (cid:19) (cid:32) sin ( 𝜋𝑥 ) sin ( 𝜋𝑦 ) sin ( 𝜋𝑥 ) sin ( 𝜋𝑦 ) (cid:33) . (74)To this end, we consider three different mesh partitionings, denoted by: • Mesh 1 , randomized quadrilateral mesh; • Mesh 2 , mainly hexagonal mesh with continuously distorted cells;22 .020.040.070.10.20.410 -9 -8 -7 -6 -5 -4 -3 -2 -1 k=1k=2k=3k=4 L re l a t i v e a pp r ox i m a t i o n err o r Mesh size h 2345 -7 -6 -5 -4 -3 -2 -1 k=1k=2k=3k=4 H re l a t i v e a pp r ox i m a t i o n err o r Mesh size h 1234
Figure 4: Convergence plots for the virtual element approximation of Problem (47)-(51) with exact solution (74) using family
Mesh 1 of randomized quadrilateral meshes. Error curves are computed using the 𝐿 norm (left panels) and 𝐻 norm (right panels) and are plotversus the number of degrees of freedom [7]. -9 -8 -7 -6 -5 -4 -3 -2 -1 k=1k=2k=3k=4 L re l a t i v e a pp r ox i m a t i o n err o r Mesh size h 2345 -7 -6 -5 -4 -3 -2 -1 k=1k=2k=3k=4 H re l a t i v e a pp r ox i m a t i o n err o r Mesh size h 1234
Figure 5: Convergence plots for the virtual element approximation of Problem (47)-(51) with exact solution (74) using family
Mesh 2 of mainly hexagonal meshes. Error curves are computed using the 𝐿 norm (left panels) and 𝐻 norm (right panels) and are plot versusthe number of degrees of freedom [7]. .020.040.070.10.20.410 -7 -6 -5 -4 -3 -2 -1 k=1k=2k=3k=4 L re l a t i v e a pp r ox i m a t i o n err o r Mesh size h 2345 -7 -6 -5 -4 -3 -2 -1 k=1k=2k=3k=4 H re l a t i v e a pp r ox i m a t i o n err o r Mesh size h 1234
Figure 6: Convergence plots for the virtual element approximation of Problem (47)-(51) with exact solution (74) using family
Mesh 3 of nonconvex octagonal meshes. Error curves are computed using the 𝐿 norm (left panels) and 𝐻 norm (right panels) and are plotversus the number of degrees of freedom [7]. -9 -8 -7 -6 -5 -4 -3 -2 -1 monomialsorthogonal L re l a t i v e a pp r ox i m a t i o n err o r -9 -8 -7 -6 -5 -4 -3 -2 -1 monomialsorthogonal H re l a t i v e a pp r ox i m a t i o n err o r Figure 7: Convergence plots for the virtual element approximation of Problem (47)-(51) with exact solution (74) using family
Mesh 1 of randomized quadrilateral meshes. Error curves are computed using k-refinement the 𝐿 norm (left panel) and 𝐻 norm (right panel)and are plot versus the number of degrees of freedom by performing a refinement of type “p” on a 5 × Mesh 3 , nonconvex octagonal mesh.The base mesh and the first refined mesh of each mesh sequence are shown in Figure 3.The discretization in time is given by applying the leap-frog method with Δ 𝑡 = − and carried out for10 time cycles in order to reach time 𝑇 = 𝑉 ℎ𝑘 with 𝑘 = , , , × 𝐼 the order of the virtual element space is increased from 𝑘 = 𝑘 = 𝑘 isgenerated by the standard scaled monomials, while in the second one we consider an orthogonal polynomialbasis. The behavior of the VEM when using nonorthogonal and orthogonal polynomial basis shown inFigure 7 is in accordance with the literature, see, e.g., [28, 75]. Acknowledgement
PFA and MV acknowledge the financial support of PRIN research grant number 201744KLJL “
VirtualElement Methods: Analysis and Applications ” funded by MIUR. PFA, IM, and MV, and SS acknowledges thefinancial support of INdAM-GNCS. GM acknowledges the financial support of the ERC Project CHANGE,which has received funding from the European Research Council under the European Union’s Horizon 2020research and innovation program (grant agreement no. 694515).
References [1] R. A. Adams and J. J. F. Fournier.
Sobolev spaces . Pure and Applied Mathematics. Academic Press, 2edition, 2003.[2] B. Ahmad, A. Alsaedi, F. Brezzi, L. D. Marini, and A. Russo. Equivalent projectors for virtual elementmethods.
Comput. Math. Appl. , 66(3):376–391, 2013.[3] F. Aldakheel, B. Hudobivnik, A. Hussein, and P. Wriggers. Phase-field modeling of brittle fracture usingan efficient virtual element scheme.
Comput. Methods Appl. Mech. Engrg. , 341:443–466, 2018.[4] P. F. Antonietti, B. Ayuso de Dios, I. Mazzieri, and A. Quarteroni. Stability analysis of discontinuousGalerkin approximations to the elastodynamics problem.
J. Sci. Comput. , 68(1):143–170, 2016.[5] P. F. Antonietti, L. Beirão da Veiga, S. Scacchi, and M. Verani. A 𝐶 virtual element method for theCahn-Hilliard equation with polygonal meshes. SIAM J. Numer. Anal. , 54(1):34–56, 2016.[6] P. F. Antonietti, F. Bonaldi, and I. Mazzieri. A high-order discontinuous Galerkin approach to theelasto-acoustic problem.
Comput. Methods Appl. Mech. Engrg. , 358:112634, 29, 2020.[7] P. F. Antonietti, G. Manzini, I. Mazzieri, H. Mourad, and M Verani. The virtual element method for linearelastodynamics models. Convergence, stability and dissipation-dispersion analysis. arXiv:1912.07122,2020. accepted on
International Journal for Numerical Methods in Engineering .[8] P. F. Antonietti, G. Manzini, and M. Verani. The fully nonconforming virtual element method forbiharmonic problems.
Math. Models Methods Appl. Sci. , 28(2):387–407, 2018.259] P. F. Antonietti, G. Manzini, and M. Verani. The conforming virtual element method for polyharmonicproblems.
Comput. Math. Appl. , 79(7):2021–2034, 2020.[10] P. F. Antonietti and I. Mazzieri. High-order discontinuous Galerkin methods for the elastodynamicsequation on polygonal and polyhedral meshes.
Comput. Methods Appl. Mech. Engrg. , 342:414–437,2018.[11] P. F. Antonietti, I. Mazzieri, A. Quarteroni, and F. Rapetti. Non-conforming high order approximationsof the elastodynamics equation.
Comput. Methods Appl. Mech. Engrg. , 209/212:212–238, 2012.[12] B. Ayuso de Dios, K. Lipnikov, and G. Manzini. The non-conforming virtual element method.
ESAIMMath. Model. Numer. , 50(3):879–904, 2016.[13] J. W. Barrett, S. Langdon, and R. Nürnberg. Finite element approximation of a sixth order nonlineardegenerate parabolic equation.
Numer. Math. , 96(3):401–434, 2004.[14] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo. Basic principles ofvirtual element methods.
Math. Models Methods Appl. Sci. , 23(1):199–214, 2013.[15] L. Beirão da Veiga, F. Brezzi, and L. D. Marini. Virtual elements for linear elasticity problems.
SIAMJ. Numer. Anal. , 51(2):794–812, 2013.[16] L. Beirão da Veiga, F. Brezzi, L. D. Marini, and A. Russo. Mixed virtual element methods for generalsecond order elliptic problems on polygonal meshes.
ESAIM: Mathematical Modelling and NumericalAnalysis , 50(3):727–747, 2016.[17] L. Beirão da Veiga, F. Brezzi, L. D. Marini, and A. Russo. Virtual element methods for general secondorder elliptic problems on polygonal meshes.
Math. Models Methods Appl. Sci. , 26(4):729–750, 2016.[18] L. Beirão da Veiga, A. Chernov, L. Mascotto, and A. Russo. Basic principles of ℎ𝑝 virtual elements onquasiuniform meshes. Math. Models Methods Appl. Sci. , 26(8):1567–1598, 2016.[19] L. Beirão da Veiga, A. Chernov, L. Mascotto, and A. Russo. Exponential convergence of the ℎ𝑝 virtualelement method in presence of corner singularities. Numer. Math. , 138(3):581–613, 2018.[20] L. Beirão da Veiga, F. Dassi, and A. Russo. A 𝐶 virtual element method on polyhedral meshes. Comput.Math. Appl. , 79(7):1936–1955, 2020.[21] L. Beirão da Veiga, K. Lipnikov, and G. Manzini. Arbitrary order nodal mimetic discretizations ofelliptic problems on polygonal meshes.
SIAM J. Numer. Anal. , 49(5):1737–1760, 2011.[22] L. Beirão da Veiga, K. Lipnikov, and G. Manzini.
The Mimetic Finite Difference Method , volume 11 of
MS&A. Modeling, Simulations and Applications . Springer, I edition, 2014.[23] L. Beirão da Veiga, C. Lovadina, and D. Mora. A virtual element method for elastic and inelasticproblems on polytope meshes.
Comput. Methods Appl. Mech. Engrg. , 295:327–346, 2015.[24] L. Beirão da Veiga and G. Manzini. A virtual element method with arbitrary regularity.
IMA J. Numer.Anal., , 34(2):782–799, 2014. DOI: 10.1093/imanum/drt018, (first published online 2013).[25] L. Beirão da Veiga and G. Manzini. Residual a posteriori error estimation for the virtual element methodfor elliptic problems.
ESAIM Math. Model. Numer. Anal. , 49(2):577–599, 2015.2626] E. Benvenuti, A. Chiozzi, G. Manzini, and N. Sukumar. Extended virtual element method for the Laplaceproblem with singularities and discontinuities.
Comput. Methods Appl. Mech. Engrg. , 356:571 – 597,2019.[27] C. Bernardi, M. Dauge, and Y. Maday. Polynomials in the Sobolev world. Technical report, HAL, 2007.hal-00153795,.[28] S. Berrone and A. Borio. Orthogonal polynomials in badly shaped polygonal elements for the virtualelement method.
Finite Elem. Anal. Des. , 129:14–31, 2017.[29] S. Berrone, A. Borio, and Manzini. SUPG stabilization for the nonconforming virtual element methodfor advection–diffusion–reaction equations.
Comput. Methods Appl. Mech. Engrg. , 340:500–529, 2018.[30] S. Berrone, S. Pieraccini, S. Scialò, and F. Vicini. A parallel solver for large scale DFN flow simulations.
SIAM J. Sci. Comput. , 37(3):C285–C306, 2015.[31] M. J. Borden, T. J. R. Hughes, C. M. Landis, and C. V. Verhoosel. A higher-order phase-field model forbrittle fracture: formulation and analysis within the isogeometric analysis framework.
Comput. MethodsAppl. Mech. Engrg. , 273:100–118, 2014.[32] J. H. Bramble and R. S. Falk. A mixed-Lagrange multiplier finite element method for the polyharmonicequation.
RAIRO Modél. Math. Anal. Numér. , 19(4):519–557, 1985.[33] S. C. Brenner and R. Scott.
The mathematical theory of finite element methods , volume 15. SpringerScience & Business Media, 2008.[34] F. Brezzi, A. Buffa, and K. Lipnikov. Mimetic finite differences for elliptic problems.
M2AN Math.Model. Numer. Anal. , 43:277–295, 2009.[35] F. Brezzi, R. S. Falk, and L. D. Marini. Basic principles of mixed virtual element methods.
ESAIMMath. Model. Numer. Anal. , 48(4):1227–1240, 2014.[36] F. Brezzi and L. D. Marini. Virtual element methods for plate bending problems.
Comput. MethodsAppl. Mech. Engrg. , 253:455–462, 2013.[37] J. W. Cahn. On spinodal decomposition.
Acta Metall. , 9:795–801, 1961.[38] J. W. Cahn and J. E. Hilliard. Free energy of a nonuniform system. I. Interfacial free energy.
The Journalof Chemical Physics , 28:258–0, 1958.[39] J. W. Cahn and J. E. Hilliard. Free energy of a nonuniform system. III. nucleation in a two-componentincompressible fluid.
The Journal of Chemical Physics , 31:688–, 1959.[40] A. Cangiani, E. H. Georgoulis, T. Pryer, and O. J. Sutton. A posteriori error estimates for the virtualelement method.
Numer. Math. , 137:857–893, 2017.[41] A. Cangiani, V. Gyrya, and G. Manzini. The non-conforming virtual element method for the Stokesequations.
SIAM J. Numer. Anal. , 54(6):3411–3435, 2016.[42] A. Cangiani, V. Gyya, G. Manzini, and Sutton. O. Chapter 14: Virtual element methods for ellipticproblems on polygonal meshes. In K. Hormann and N. Sukumar, editors,
Generalized BarycentricCoordinates in Computer Graphics and Computational Mechanics , pages 1–20. CRC Press, Taylor &Francis Group, 2017. 2743] A. Cangiani, G. Manzini, A. Russo, and N. Sukumar. Hourglass stabilization of the virtual elementmethod.
Internat. J. Numer. Methods Engrg. , 102(3-4):404–436, 2015.[44] A. Cangiani, G. Manzini, and O. Sutton. Conforming and nonconforming virtual element methods forelliptic problems.
IMA J. Numer. Anal., , 37:1317–1354, 2017. (online August 2016).[45] O. Certik, F. Gardini, G. Manzini, L. Mascotto, and G. Vacca. The p- and hp-versions of the virtualelement method for elliptic eigenvalue problems.
Comput. Math. Appl. , 79(7):2035–2056, 2020.[46] O. Certik, F. Gardini, G. Manzini, and G. Vacca. The virtual element method for eigenvalue problemswith potential terms on polytopic meshes.
Applications of Mathematics , 63(3):333–365, 2018.[47] F. Chave, D. A. Di Pietro, F. Marche, and F. Pigeonneau. A hybrid high-order method for the Cahn-Hilliard problem in mixed form.
SIAM J. Numer. Anal. , 54(3):1873–1898, 2016.[48] F. Chen and J. Shen. Efficient energy stable schemes with spectral discretization in space for anisotropicCahn-Hilliard systems.
Communications in Computational Physics , 13(5):1189–1208, 2013.[49] L. Chen and X. Huang. Nonconforming virtual element method for 2 𝑚 th order partial differentialequations in R 𝑛 . Math. Comp. , 89(324):1711–1744, 2020.[50] C. Chinosi and L. D. Marini. Virtual element method for fourth order problems: 𝐿 -estimates. Comput.Math. Appl. , 72(8):1959–1967, 2016.[51] F. Dassi and L. Mascotto. Exploring high-order three dimensional virtual elements: bases and stabiliza-tions.
Comput. Math. Appl. , 75(9):3379–3401, 2018.[52] D. A. Di Pietro, J. Droniou, and G. Manzini. Discontinuous skeletal gradient discretisation methods onpolytopal meshes.
J. Comput. Phys. , 355:397–425, 2018.[53] C. M. Elliott and D. A. French. Numerical studies of the Cahn-Hilliard equation for phase separation.
IMA J. Appl. Math. , 38(2):97–128, 1987.[54] C. M. Elliott and D. A. French. A nonconforming finite-element method for the two-dimensionalCahn-Hilliard equation.
SIAM J. Numer. Anal. , 26(4):884–903, 1989.[55] C. M. Elliott, D. A. French, and F. A. Milner. A second-order splitting method for the Cahn-Hilliardequation.
Numer. Math. , 54(5):575–590, 1989.[56] C. M. Elliott and S. Larsson. Error estimates with smooth and nonsmooth data for a finite elementmethod for the Cahn-Hilliard equation.
Math. Comp. , 58(198):603–630, S33–S36, 1992.[57] C. M. Elliott and Z. Songmu. On the Cahn-Hilliard equation.
Arch. Rational Mech. Anal. , 96(4):339–357,1986.[58] E. Faccioli, F. Maggio, A. Quarteroni, and A. Taghan. Spectral-domain decomposition methods forthe solution of acoustic and elastic wave equations.
The Leading Edge , 61:1160–1174, 07 1996.Faccioli1996.[59] D. Gallistl. Stable splitting of polyharmonic operators by generalized Stokes systems.
Math. Comp. ,86(308):2555–2577, 2017.[60] F. Gardini, G. Manzini, and G. Vacca. The nonconforming virtual element method for eigenvalueproblems.
ESAIM Math. Model. Numer. , 53:749–774, 2019. Accepted for publication: 29 November2018. DOI-DUMMY: 10.1051/m2an/2018074. 2861] F. Gazzola, H.-C. Grunau, and G. Sweers.
Polyharmonic boundary value problems , volume 1991 of
Lecture Notes in Mathematics . Springer-Verlag, Berlin, 2010. Positivity preserving and nonlinear higherorder elliptic equations in bounded domains.[62] H. Gómez, V. M. Calo, Y. Bazilevs, and T. J. R. Hughes. Isogeometric analysis of the Cahn-Hilliardphase-field model.
Comput. Methods Appl. Mech. Engrg. , 197(49-50):4333–4352, 2008.[63] P. Grisvard.
Elliptic problems in nonsmooth domains , volume 24 of
Monographs and Studies inMathematics . Pitman (Advanced Publishing Program), Boston, MA, 1985.[64] T. Gudi and M. Neilan. An interior penalty method for a sixth-order elliptic equation.
IMA J. Numer.Anal. , 31(4):1734–1753, 2011.[65] D. Kay, V. Styles, and E. Süli. Discontinuous Galerkin finite element approximation of the Cahn-Hilliardequation with convection.
SIAM J. Numer. Anal. , 47(4):2660–2685, 2009.[66] D. Komatitsch and J. Tromp. Introduction to the spectral element method for three-dimensional seismicwave propagation.
Geophysical Journal International , 139(3):806–822, 1999.[67] D. J. Korteweg. Sur la forme que prenent les équations du mouvements des fluides si l’on tient comptedes forces capilaires causées par des variations de densité considérables mains continues et sur la théoriede la capillarité dans l’hypothése d’une varation continue de la densité. Arch. Néerl Sci. Exactes Nat.Ser. II, 1901.[68] L. D. Landau. On the theory of superconductivity. In D. ter Haar, editor,
Collected papers of L. D.Landau , pages 546–568. Pergamon, 1965.[69] K. Lipnikov and G. Manzini. A high-order mimetic method for unstructured polyhedral meshes.
J.Comput. Phys. , 272:360–385, 2014.[70] K. Lipnikov, G. Manzini, and M. Shashkov. Mimetic finite difference method.
J. Comput. Phys. , 257 –Part B:1163–1227, 2014.[71] X. Liu and Z. Chen. A virtual element method for the Cahn-Hilliard problem in mixed form.
Appl.Math. Lett. , 87:115–124, 2019.[72] C. Lovadina, D. Mora, and I. Velásquez. A virtual element method for the von Kármán equations.Technical report, Preprint CI2MA:2019-36, 2019.[73] G. Manzini, K. Lipnikov, J. D. Moulton, and M. Shashkov. Convergence analysis of the mimetic finitedifference method for elliptic problems with staggered discretizations of diffusion coefficients.
SIAM J.Numer. Anal. , 55(6):2956–2981, 2017.[74] G. Manzini, A. Russo, and N. Sukumar. New perspectives on polygonal and polyhedral finite elementmethods.
Math. Models Methods Appl. Sci , 24(8):1621–1663, 2014.[75] L. Mascotto. Ill-conditioning in the virtual element method: stabilizations and bases.
Numer. MethodsPartial Differential Equations , 34(4):1258–1281, 2018.[76] D. Mora, G. Rivera, and R. Rodríguez. A virtual element method for the Steklov eigenvalue problem.
Math. Methods Appl. Sci. , 25(08):1421–1445, 2015.[77] D. Mora, G. Rivera, and I. Velásquez. A virtual element method for the vibration problem of Kirchhoffplates.
ESAIM Math. Model. Numer. Anal. , 52(4):1437–1456, 2018.2978] D. Mora and I. Velásquez. A virtual element method for the transmission eigenvalue problem.
Math.Models Methods Appl. Sci. , 28(14):2803–2831, 2018.[79] D. Mora and I. Velásquez. Virtual element for the buckling problem of Kirchhoff-Love plates.
Comput.Methods Appl. Mech. Engrg. , 360:112687, 22, 2020.[80] K. Park, H. Chi, and G. H. Paulino. On nonconvex meshes for elastodynamics using virtual elementmethods with explicit time integration.
Comput. Methods Appl. Mech. Engrg. , 356:669–684, 2019.[81] K. Park, H. Chi, and G. H. Paulino. Numerical recipes for elastodynamic virtual element methods withexplicit time integration.
Internat. J. Numer. Methods Engrg. , 121(1):1–31, 2020.[82] G. H. Paulino and A. L. Gain. Bridging art and engineering using Escher-based virtual elements.
Struct.and Multidisciplinary Optim. , 51(4):867–883, 2015.[83] I. Perugia, P. Pietra, and A. Russo. A plane wave virtual element method for the Helmholtz problem.
ESAIM Math. Model. Num. , 50(3):783–808, 2016.[84] A. Quarteroni, R. Sacco, and F. Saleri.
Numerical Mathematics , volume Vol. 37 of
Texts in AppliedMathematics . Springer, 2007.[85] P.-A. Raviart and J.-M. Thomas.
Introduction à l’analyse numérique des équations aux dérivées par-tielles . Collection Mathématiques Appliquées pour la Maîtrise. [Collection of Applied Mathematics forthe Master’s Degree]. Masson, Paris, 1983.[86] B. Rivière and M. F. Wheeler. Discontinuous finite element methods for acoustic and elastic waveproblems. In
Current trends in scientific computing (Xi’an, 2002) , volume 329 of
Contemp. Math. ,pages 271–282. Amer. Math. Soc., Providence, RI, 2003.[87] J. S. Rowlinson. Translation of J. D. van der Waals’ “The thermodynamic theory of capillarity underthe hypothesis of a continuous variation of density”.
J. Statist. Phys. , 20(2):197–244, 1979.[88] M. Schedensack. A new discretization for 𝑚 th-Laplace equations with arbitrary polynomial degrees. SIAM J. Numer. Anal. , 54(4):2138–2162, 2016.[89] S. Torabi, J. Lowengrub, A. Voigt, and S. Wise. A new phase-field model for strongly anisotropic systems.
Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. , 465(2105):1337–1359, 2009. With supplementarymaterial available online.[90] G. Vacca. Virtual element methods for hyperbolic problems on polygonal meshes.
Comput. Math. Appl. ,74(5):882–898, 2017.[91] M. Wang and J. Xu. Minimal finite element spaces for 2 𝑚 -th-order partial differential equations in 𝑅 𝑛 . Math. Comp. , 82(281):25–43, 2013.[92] G. N. Wells, E. Kuhl, and K. Garikipati. A discontinuous Galerkin method for the Cahn-Hilliardequation.
J. Comput. Phys. , 218(2):860–877, 2006.[93] P. Wriggers, W. T. Rust, and B. D. Reddy. A virtual element method for contact.
Comput. Mech. ,58(6):1039–1050, 2016.[94] J. Zhao, S. Chen, and B. Zhang. The nonconforming virtual element method for plate bending problems.
Math. Models Methods Appl. Sci. , 26(9):1671–1687, 2016.[95] J. Zhao, B. Zhang, S. Chen, and S. Mao. The Morley-type virtual element for plate bending problems.