A guide to the design of the virtual element methods for second-and fourth-order partial di ﬀ erential equations

: We discuss the design and implementation details of two conforming virtual element methods for the numerical approximation of two partial di ﬀ erential equations that emerge in phase-ﬁeld modeling of fracture propagation in elastic material. The two partial di ﬀ erential equations are: (i) a linear hyperbolic equation describing the momentum balance and (ii) a fourth-order elliptic equation modeling the damage of the material. Inspired by [1–3], we develop a new conforming VEM for the discretization of the two equations, which is implementation-friendly, i


Introduction
The virtual element method (VEM) is a generalization of the finite element method (FEM) that enables computational simulations using n-sided polygonal/polyhedral elements (with n ≥ 3 in 2D and n ≥ 4 in 3D), significantly reducing the difficulty of meshing geometrically complex domains.In addition, the VEM provides a means of designing approximations of an arbitrarily high order of accuracy and regularity [4], thus making it ideal for problems admitting C 1 -continuous solutions.The VEM is also known to be accurate, stable, and robust on highly deformed meshes, as shown in numerous theoretical [5][6][7] and numerical [8][9][10] studies on the treatment of large-deformation problems.The method preserves these desirable performance characteristics when the mesh includes elements with obtuse or re-entrant angles (i.e., nonconvexity) and hanging nodes, facilitating the development of adaptive mesh refinement (AMR) schemes to avoid excessive mesh refinement [11].
Despite these excellent properties and an established and solid mathematical foundation, which is well documented in the literature [12][13][14][15][16][17][18], the VEM lags behind the FEM's adoption in practical applications mainly because it might be perceived as being conceptually more complex and, therefore, more challenging to implement and apply.Several works in the literature help in addressing this issue (e.g., see [2,19]) sometimes by presenting simple virtual element formulations and detailing their implementations (e.g., see [20][21][22]).In this work, we address the implementation design of two distinct virtual element methods for approximating a time-dependent, second-order momentum balance equation and a stationary, fourth-order, general elliptic equation, which includes second-and zero-order terms.The VEM is developed as a numerical method to approximate different model problems on polygonal elements and could be seen as a generalization of the FEM.Unlike the FEM, the discrete functions are formally the solutions of a local PDE that is defined on every mesh element Instead of directly evaluating the discrete functions, which are never explicitly constructed, we employ suitable projection operators, which are computable from the degrees of freedom, and project such functions onto a polynomial space.Consequently, the matrices associated with the discrete-bilinear forms involve projection operators.The VEM allows us to compute the local matrices on different polygonal elements in a unified way and the scheme's formulation and implementation do not depend on the shape of the elements.These features make the VEM different from the existing techniques.However, some issues may exist, such as the stability related to the stiffness and mass matrices.If the polygonal element contains very small edges compared to the diameter of the element, the dofi-dofi stabilization will not work, and we need boundary integration of the trace of the basis functions, which is more technical to compute [5,23].In this work, we have assumed the meshes do not contain small edges and all polygonal elements are star-shaped.
The momentum balance equation governs the linear elastodynamic initial boundary-value problem (IBVP); the second equation mathematically models crack propagation in materials.When coupled, these two PDEs constitute fourth-order phase-field models of quasi-brittle fractures, similar to those presented in [24,25].Part of the authors has previously presented VEMs for the solution of the first equation in [26], while the virtual element discretization of the second equation is currently under study, see, for instance, the work in progress of [27].These previous works focused on both methods' theoretical aspects and convergence analysis, providing convergence proofs, error estimates, and the numerical investigation of their convergence behavior in representative benchmark test case applications.Instead, herein we focus on the implementation design for each equation separately, emphasizing similar computational aspects.At the same time, a forthcoming publication will address the VEM-based coupled approach to solving the fracture model.
The remainder of this paper is as follows.We present the mathematical formulation of these two model problems in the strong and weak form in Section 2 and their numerical discretizations using the VEM in Section 3. We discuss the implementation design in Section 4 and offer our concluding remarks in Section 5.

Model problems and variational formulations
In this section, we introduce the model problems that we consider in this paper, and their variational formulation, which we will discretize using the virtual element method.In Subsection 2.1, we introduce the notation that is used throughout the paper; in Subsection 2.2, we present the model problems in strong form; in Subsection 2.3, we present the model problems in variational form.

Notation and technicalities
Throughout this paper, we adopt the notation of Sobolev spaces as in [28].Accordingly, we denote the space of square-integrable functions defined on any open, bounded, connected domain Ω ⊂ R 2 with boundary ∂Ω by L 2 (Ω), and the Hilbert space of functions in L 2 (Ω) with all partial derivatives up to a positive integer m also in L 2 (Ω) by H m (Ω), see [28].We endow H m (Ω) with the usual norm and seminorm that we denote as The virtual element methods that we consider in the next sections are formulated on the mesh family Ω h h , where each mesh Ω h is a partition of the computational domain Ω into nonoverlapping star-shaped polygonal elements E with boundary ∂E, area |E|, center of gravity x E , and diameter h E = sup x,y∈E |x − y|.The mesh elements of Ω h form a finite cover of Ω such that Ω = ∪ E∈Ω h E and the mesh size labeling each mesh Ω h is defined by h = max E∈Ω h h E .A mesh edge e has center x e and length h e ; a mesh vertex v has position vector x v .We denote the edges of the polygonal boundary ∂E by E E .
For any integer number ≥ 0, we let P (E) and P (e) denote the space of polynomials defined on the element E and the edge e, respectively; P (Ω h ) denotes the space of piecewise polynomials of degree on the mesh Ω h .For the convenience of exposition, we also use the notation.P −2 (E) = P −1 (E) = {0}.Accordingly, it holds that q |E ∈ P (E) if E ∈ Ω h for all q ∈ P (Ω h ).Throughout the paper, we use the multi-index notation, so that ν = (ν 1 , ν 2 ) is a two-dimensional index defined by the two integer numbers ν 1 , ν 2 ≥ 0.Moreover, D ν w = ∂ |ν| w/∂x ν 1 1 ∂x ν 2 2 denotes the partial derivative of order |ν| = ν 1 + ν 2 > 0 of a sufficiently regular function w(x 1 , x 2 ), and we use the conventional notation that D (0,0) w = w for ν = (0, 0).We also denote the partial derivatives of w versus x and y by the shortcuts ∂ x w and ∂ y w, and the normal and tangential derivatives with respect to a given edge by ∂ n w and ∂ t w.Finally, we define the trace operator tr (A) = i a ii for any matrix (two-index tensor) A = a i j .Given two tensors A = (a i j ) and B = (b i j ), we also define A : B := i j a i j b i j .

Linear momentum balance equation
The linear momentum balance equation for the displacement u : Ω × (0, T ] → R 2 on the 2-D computational domain Ω reads as where ρ : Ω → R is a scalar density field, ü denotes the second-order time derivative of u, σ : sym is the (time-dependent) stress tensor (R 2×2 sym being the set of symmetric 2 × 2-sized real-valued tensors), and f : Ω × (0, T ] → R 2 is a suitable vector-valued forcing term.In this paper, we assume the density field is constant.The corresponding stress tensor σ is given by Mathematics in Engineering Volume 5, Issue 6, 1-22. where g u (x) : R 2 → R is a given function modeling different material behaviors such as hardening or softening effects, for instance, due to damage or fracture, and ε vol and ε dev are the volumetric and deviatoric strain tensors.Tensors ε vol and ε dev are respectively given by tr (ε(u)) I, and where I the 2 × 2 identity tensor, and the volumetric strain tensor can also be written as ε vol = ∇ • u I.
To complete the mathematical formulation of this model problem, we supplement Eq (2.1) with the initial and boundary conditions.The initial conditions on the displacement and its derivative at t = 0 are given by Herein, u 0 , v 0 , u, and g N are given data; n = n x , n y T is the unit vector orthogonal to Γ and pointing out of Ω.

Fourth-order phase-field equation
The fourth-order equation for the phase field d on the 2-D computational domain Ω reads as where α 2 , α 1 , and α 0 are model parameters, g d (d) : R → R is an affine function, and the history variable H t : Ω → R is an additional multiplicative term.We assume that α 2 and α 1 are strictly positive real numbers, and α 0 is a non-negative real number.We remark that g d (d) := α 3 d + α 4 , where α 3 and α 4 are constants.For simplicity, we let α 4 = 0 in this work.We refer to [12] for a detailed discussion about the physical meaning of g d

Linear momentum balance equation
Consider the affine functional space of the two-dimensional vector-valued fields whose components belong to the Sobolev space H 1 (Ω), and have an assigned value on the boundary Dirichlet boundary Γ D,u : Denoting V u,0 if u = 0 and testing Eq (2.1) with elements in V u,0 , we arrive at the variational formulation of the elastodynamics equation as with the definitions: where we denote the scalar product of two (matrix) tensors by ":".In particular, for any tensor fields A = (a il ) i, j=1,2 , and B = (b il ) i, j=1,2 , we consider the standard scalar product of 2 × 2-matrices: A : B = i, j=1,2 a i j b i j .As shown, for example, by Raviart and Thomas (see [13,) the variational problem (2.8) is well-posed and its unique solution satisfies

Fourth-order phase-field equation
Similar to the previous section, we consider the affine functional space and a special case of V d , namely, V d,0 , given by setting d = 0 in the previous definition.On testing (2.6) with elements in V d,0 , we obtain the variational formulation of the fourth-order phase-field equation. where 12) The wellposedness of the variational formulation, i.e., existence and uniqueness of d ∈ H 2 (Ω) solving (2.10) with definitions (2.11)-(2.12),follows from applications of the Lax-Milgram theorem [14], since the bilinear form A d (•, •) is coercive and continuous.

Numerical discretizations
In this section, we briefly review the formulation of the virtual element method that approximates the second-order momentum balance equation and the high-order phase field equation introduced in the previous section.

VEM for the second-order momentum balance equation
Let V h u,k be the virtual element approximation of the affine functional space V u defined in (2.7).The semi-discrete virtual element approximation of (2.8) reads as: Here, u h (t) is the virtual element approximation of u and v h the generic test function in V h k , while (u 0 ) I and (u 1 ) I are the virtual element interpolants of the initial solution functions u(0) and u(0), respectively.The bilinear forms M h u (•, •), A h u (•, •), and the linear functional The total time interval [0, T ] is divided into N T subintervals with time step ∆t = T N T and time instants t (n) = n∆t.The Newmark method [15] is widely used to solve elastodynamic problems because it allows us to solve the second-order equation directly without splitting the equation into two first-order equations.Therefore, we employ the Newmark method to advance the system in time.In this implicit scheme, accelerations and velocities at time t (n+1) are approximated as follows: .. .. . ..

Degrees of freedom, projection operators, local and global spaces
The vector-valued virtual element space is such that , where V h u,k (E) denotes the local scalar conforming virtual element space (which will be introduced below).Accordingly, the vector-valued displacement field u h is provided by the two components u h = u x,h , u y,h , where each component belongs to V h u,k (E).Both components u i,h , for i ∈ {x, y}, are uniquely identified by the following degrees of freedom: (U1) the values of u i,h at the vertices of E; (U2) the edge polynomial moments of u i,h of order up to k − 2 on each one-dimensional edge e ∈ E E : (U3) the cell polynomial moments of u i,h of order up to k − 2 on E: 1 In this definition, we use the symbols M k−2 (e) and M k−2 (E) to denote a basis for the polynomial spaces P k (e) and P k−2 (E), respectively.Suitably scaled monomials or orthogonal polynomials can provide this basis.
Figure 1 shows the degrees-of-freedom of the three scalar conforming virtual element spaces (k = 1, 2, 3) defined on a pentagonal cell.Using these degrees of freedom, we compute the elliptic projection operator where the second condition is needed to fix the kernel of the gradient operator.The local scalar virtual element space of order k ≥ 1 is defined by where P k (E)\P k−2 (E) denotes as the space spanned by the monomials of degree k and k − 1 on E, and The degrees-of-freedom of the scalar conforming virtual element spaces V h u,k (E), (k = 1, 2, 3) on a pentagonal cell, which approximates the components of the vector-valued field u, solving the momentum balance equation.Green circles, red squares, and blue crosses are (U1), (U2), and (U3), respectively.
According to this definition, the orthogonal projection Π 0,E u,k c h of a virtual element function c h onto the polynomial space of degree k, defined as, is also computable.The global conforming space of vector-valued virtual element field of order k ≥ 1, i.e., the finite-dimensional space V h u,k , is obtained by combining all the elemental spaces V h u,k (E) 2 .
Building upon the local spaces where S E m (•, •) is the local stabilization term for M h,E u .The bilinear form M E u depends on the orthogonal projections Π 0 u,k u h and Π 0 u,k v h , which are computable from the degrees of freedom of u h and v h , respectively, see the previous section.The local bilinear form A h,E u (•, •) is given by where S E a (•, •) is the local stabilization term for A E u .The bilinear form A E u depends on the orthogonal projections Π 0,E u,k−1 σ(u h ) and Π 0,E u,k−1 ε(v h ), which are computable from the degrees of freedom of u h and v h , respectively, see the previous section.
The local stabilization terms u,k → R can be any symmetric and coercive bilinear forms that are computable from the degrees of freedom, and suitably provides the k-consistency and stability property.The role of the stabilization in the VEM is discussed in detail in [23].The stabilization used in our implementation is discussed in Subsection 4.2.Finally, we approximate the right-hand side (3.1) of the semi-discrete formulation as follows: where The linear functional F h u (•) is clearly computable since the edge trace v h | e is a polynomial and Π 0,E u,k (v h ) is computable from the degrees of freedom of v h .

VEM for the fourth-order phase-field equation
The virtual element method approximates the variational formulation (2.10) as follows In this formulation, V h d,r is the H 2 -conforming approximation of the space H 2 (Ω) provided by the VEM, d h and c h are the trial and test functions from this space.

Degrees of freedom, projection operators, local and global spaces
On every mesh element E ∈ Ω h , the following set of values comprises the degrees of freedom of a virtual element function c h as shown in Figure 2: 1 |e| e qc h dS for any q ∈ P r−4 (e), and any edge e ∈ E E ; (D3) for r ≥ 3, e q∂ n c h dS for any q ∈ P r−3 (e), and any edge e ∈ E E ; (D4) for r ≥ 2, 1 |E| E qv dV for any q ∈ P r−2 (E).
Consider the integer r ≥ 2 and the bilinear form We define the elliptic projection operator Π L,E d,r : H 2 (E) → P r (E) such that for every v ∈ H 2 (E), the r-degree polynomial Π L,E d,r v is the solution to the variational problem: The degrees-of-freedom of the scalar conforming virtual element spaces V h d,r (E), (r = 1, 2, 3) on a pentagonal cell, which approximates the scalar field solving the high-order phase-field equation.The solid green circles and empty red cirles are (D1), solid red squares are (D2), empty red squares are (D3), and blue crosses are (D4).

The elliptic projection operator
is the the solution of the finite-dimensional variational problem B(Π L,E d,r d, q) = B(d, q), q ∈ P r (E), (3.16) with the following additional conditions and for i = {x, y}: For r ≥ 3, the local virtual element space on element E is defined by where For r = 2, we have the special case of the low-order virtual element space: Note that P r (E) ⊂ V h d,r (E).Building upon the local spaces V h d,r (E), r ≥ 2, for all E ∈ Ω h the global conforming virtual element space V h d,r is defined as This implies that V h d,r ⊂ C 1 (Ω) because d h ∈ H 2 (E) in each element and C 1 -regularity across the internal mesh faces.

Virtual element bilinear forms and linear functional
Let E ∈ Ω h be a mesh element, and consider the bilinear forms R given by integrating on E instead of Ω in the corresponding bilinear forms in (2.11)-(2.12).Let We use the elliptic projection Π L,E d,r and the L 2orthogonal projection Π 0,E d,r to define the virtual element bilinear form Herein, the stabilization term is also built by using the projection Π L,E d,r (c h ), and the usual formula, so that the bilinear form S E h : V h d,r (E) × V h d,r (E) → R can be any symmetric, positive definite, bilinear form that suitably provides the k-consistency and stability properties.The stabilization used in our implementation is discussed in Subsection 4.3.

Implementation
In this section, we follow the general guidelines of [19] and briefly describe how we implemented the VEM for the second-order elastodynamics equation and the fourth-order phase-field equation, respectively.

Vector and matrix notation
In this section, we introduce VEM vector and matrix notations that will be used in the following sections.We have introduced two virtual element spaces, H 1 -conforming V h u,k and H 2 -conforming V h d,r .For conciseness, we use the virtual element space V h u,k as an example to introduce vector and matrix notations.Such notations can be easily extended to V h d,r .We consider the following compact notation.For all element E ∈ Ω h , we locally number the degrees of freedom (U1), (U2), and (U3) from 1 to N dofs .Then, we introduce the bounded, linear functionals i denote the set of such functionals and collect the degrees of freedom of v h in the vector v h = dof 1 (v h ) , . . ., dof N dofs (v h ) T .Since the degrees of freedom (U1), (U2), and (U3) are unisolvent in V h u,k (E), the triplet E, V h u,k (E), Λ E is a finite element in the sense of Ciarlet [16].This property implies the existence of a Lagrangian basis set ϕ u i i , with We refer to the basis function set ϕ u i i as the "canonical" basis of V h u,k (E).We introduce the compact notation and write the expansion of a virtual element function v h on such a basis set as We also introduce a compact notation for the basis of the polynomial subspace where n k is the cardinality of P k (E).Since the polynomials m u α (x) are also virtual element functions, we can expand them on the canonical basis ϕ ϕ ϕ.We express such expansions as
where matrix D D D u has size N dofs × n k and collects all the expansion coefficients Following this notation, we also express the action of a differential operator D, e.g., D = ∆ or D = ∇, in a entry-wise way, so that and Similarly, we express the action of the projectors Π ∇,E u,k , and Π 0,E u,k on the canonical basis functions ϕ ϕ ϕ u and their expansion on the polynomial basis m m m u as follows: where Π Π Π ∇,E u,k and Π Π Π 0,E u,k are the matrix representation of Π ∇,E u,k and Π 0,E u,k , respectively.The expansion coefficients for the three projection operators applied to the basis function ϕ j are collected along the j-th column of the projection matrices Π Π Π ∇,E u,k , and Π Π Π 0,E u,k .Remark 1.The vector and matrix notation for the fourth-order phase-field equation can be written in a similar fashion, such as, the virtual element space V h d,r (E), degree of freedoms (D1)-(D4), canonical basis ϕ ϕ ϕ d , the polynomial basis m m m d , and more.

VEM implementation for the linear momentum equation
In this section, we start with constructing elliptic and L 2 projectors, followed by assembling local mass and stiffness matrices.

Elliptic projector
The elliptic projector Π ∇,E u,k is defined as, for (4.1) Using vector and matrix notation as introduced in Section 4.1, we have
As a consequence, we have the following matrix equation where However, G G G u is singular.We additionally define Collecting Eqs (4.4) and (4.5), we obtain It is worth noting that G G G u is nonsingular by construction, and we can formally state that

L 2 orthogonal projectors
Recall the definition of the L 2 -orthogonal projector Π 0,E u,k acting on a scalar Rewrite Eq (4.7) using vector and matrix notation as We replace The above equation can be rewritten as where and We remark that H u is nonsingular, therefore Π Π Π 0,E u,k = H −1 u C u .In addition, the L 2 -orthogonal projector acting on ∇v h , for We can choose q = q up , 0 T or q = 0, q down T where q up , q down ∈ P k−1 (E).The above equation can be split into the following two equations, We start with the first equation, As a result, the above equation can be cast into where and Following a similar procedure, we have the matrix equation for where It is easy to see that H u,k−1 is invertible, therefore,

Local matrices
After constructing projection matrices, we can assemble the local mass and stiffness matrices.We follow two major ideas in this section.The first idea is that we split the vector-valued function into two scalar-valued components, i.e., 2 , and u = u up , 0 T + 0, u down T , where u up , Second, it is well-known that we can split the local matrix into consistency and stability terms in VEM.
The consistency term is calculated using polynomials, while the stability term is approximated using the local degrees of freedom.
Using the splitting idea Eq (4.19), we can rewrite the mass term as in the variational form Eq (2.8) as The local mass matrix is defined as We neglect the stability term of the local mass matrix, namely, M M M u,stab E = 0. Therefore, the local mass matrix is given by Using the L 2 -projector Π 0,E u,k and adopting the vector and matrix notations, we have where we have assumed the density of the material, ρ, is a constant.Next, substituting Eqs (2.2) and (2.3) into the stiffness term in the variational formulation Eq (2.8), we arrive at where Moreover, using Eq (4.19) and expanding ε(u), we further obtain
+ µ and Now, we are ready to assemble the stiffness matrix, which consists of consistency and stability terms as From (4.23), we obtain the consistency term We first address the deviatoric terms, and the volumetric components follow similarly.Using Π 0,x,E u,k−1 , Π 0,y,E u,k−1 , and Eq (4.23), we obtain where It is evident that Then, the components of the volumetric term of the local stiffness matrix can be derived from Eq (4.24) as, Finally, we follow [26] and adopt the stability term of the stiffness matrix as The stability matrix K K K u,stab E should scale with the consistency matrix K K K u,cons E , which depends on the Lame parameters.Therefore, we let K K K u,stab E scale with max(2µ, λ).We remark that more options to design the stability term can be found in [17,18], and [29,30] for nearly incompressible material, where the stabilization term depends on the Lame parameters.

Right-hand side approximation
Using the compact notation again, we rewrite the virtual element approximation of the first term of F h u , see Eq (3.11), as (4.28)

VEM implementation for the fourth-order phase field equation
In this section, we follow a similar procedure as in Section 4.2 and discuss how to construct elliptic and L 2 -orthogonal projectors, and the fourth-order equation's local mass and stiffness matrices.

Elliptic projectors
In this section, we first construct Π ∇,E d,r , followed by Π L,E d,r .The elliptic projector Π ∇,E d,r is defined as, for v h ∈ V h d,k (E) The above equation can be recast into the following matrix equation: where To fix the kernel of the differential operator ∇, we additionally require Next, from the definition of Π L,E d,r , where Recall that Π 0,E d,r is defined as, for v h ∈ V h d,k (E) The construction of the L 2 projection matrix Π Π Π 0,E d,r is similar to the derivation of Eq (4.10), therefore Note that we do not need to specify stabilization terms in the local mass matrices.The stiffness matrix is given by the sum of two terms: a rank-deficient consistency term, which guarantees the consistency of the approximation, and a stability term, which fixes the correct rank: It is worth noting that we can also use other stabilizations as in [31].

Conclusions
In this work, we present, in detail, using vector and matrix notation, how to implement two distinct virtual element methods for approximating a time-dependent, second-order momentum balance equation and a fourth-order elliptic equation.Specifically for the momentum equation, we discuss the implementation details by adopting a deviatoric and volumetric strain split.Such split is essential to model various material behavior such as hardening and softening.
The momentum balance equation governs the linear elastodynamic IBVP; the second equation mathematically models crack propagation in materials.Future works include coupling the secondorder elastodynamic equation with the fourth-order phase-field equation.When connected, these two PDEs constitute fourth-order phase-field models of dynamic-brittle fractures, similar to those presented in [24,25].We plan to solve the coupled equations using virtual element methods discussed in the work and study dynamic brittle fracture problems.

Use of AI tools declaration
The authors declare they have not used any Artificial Intelligence (AI) tools in the creation of this article.
(d) and H t when modeling crack initiation and propagation phenomena.To set the boundary conditions that establish the mathematical model, we first split the domain boundary as Γ = Γ D,d ∪ Γ N,d , where Γ D,d and Γ N,d represents the Dirichlet and Neumann parts of the domain boundary, respectively.Subsets Γ D,d and Γ N,d are disjoint in the sense that their intersection has zero 2D measure, i.e., Γ D,d ∩ Γ N,d = 0. Additionally, we split the Dirichlet boundary as Γ D,d = Γ D 1 ,d ∪Γ D 2 ,d and the Neumann boundary as Γ N,d as Γ N,d = Γ N 1 ,d ∪Γ N 2 ,d .We complete the mathematical formulation of the model problem (2.6) with the following boundary conditions d = d on Γ D 1 ,d , and ∇d = 0 on Γ D 2 ,d , ∆d = 0 on Γ N 1 ,d , and ∇∆d = 0 on Γ N 2 ,d .

2 E ∆m m m d ∆ϕ ϕ ϕ T d dV + α 1 EEα 2 ∆
double) integration by parts of the right-hand-side of the definition of B B B d yieldsα 2 B B B d,2 + α 1 B B B d,1 = α ∇m m m d • ∇ϕ ϕ ϕ T d dV = 2 m m m d − α 1 ∆m m m d ϕ ϕ ϕ d dV + ∂E ∂ n (−α 2 m m m d + α 1 m m m d ) ϕ ϕ ϕ T d + α 2 ∆m m m d ∂ n ϕ ϕ ϕ T d dS .Collecting Eqs (4.30) and (4.32), we haveG G G d Π Π Π L,E d,r = B B B d , with G G G d = G G G d + G G G 0 d , B B B d = B B B d + B B B 0 d .(4.33)It is worth noting that Eq (4.30) also fixes the kernel of L, namely,G G G 0 d Π Π Π L,E d,r = B B B 0 d .Therefore, G G G d,1and G G G d are nonsingular by construction, and we can formally state thatΠ Π Π ∇,E d,r = G G G −1 d,1 B B B d,1 , and Π Π Π L,E d,r = G G G −1 d B B B d .(4.34) 4.3.2.L 2 -orthogonal projector .35) where H d := E m m m d m m m T d dV, and C d := E m m m d ϕ ϕ ϕ T d dV.4.3.3.Local matricesThe construction of the local mass matrices uses the L 2 -orthogonal projection, Π 0,E d,r .Following a similar procedure to obtain Eq (4.22), we haveM M M d E := M M M d,cons E = α 0 A d,0 Π 0,E d,k ϕ ϕ ϕ d , Π 0,E d,k ϕ ϕ ϕ T d = α 0 Π Π Π 0,E d,k T H d Π Π Π 0,E d,k, andM M M d,H t E := M M M d,H t ,cons E = A g d Π 0,E d,k ϕ ϕ ϕ d , Π 0,E d,k ϕ ϕ ϕ T d = Π Π Π 0,E d,k T H d,H t Π Π Π 0,E d,k , where we have assumed the linear function g d (d) = d and H d,H t := E H t (x) m m m d m m m T d dV.

= α 2
A d,2 Π L,E d,r ϕ ϕ ϕ d , Π L,E d,r ϕ ϕ ϕ d + α 1 A d,1 Π L,E d,r ϕ ϕ ϕ d , Π L,E d,r ϕ ϕ ϕ d = Π Π Π L,E d,r T G G G d Π Π Π L,Ed,r and K K K stab E = I I I − D D D d Π Π Π L,E d,r T I I I − D D D d Π Π Π L,E d,r .
Γ D,u that is a closed subset of Γ in the Euclidean topology, and Neumann boundary, e.g., Γ N,u , such that Γ = Γ D,u ∪ Γ N,u and Γ D,u ∩ Γ N,u = 0. Then the boundary conditions read as: .8) 3.1.2.Virtual element bilinear forms and linear functional In the virtual element setting, we define the bilinear forms M h u (•, •) and A h u (•, •) as the sum of elemental contributions, which are denoted by M E u (•, •) and A E u (•, •), respectively: