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

1.   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\geq3$ in 2D and $n\geq4$ 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.

2.   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.

2.1.   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 $\Omega\subset \mathbb{R}^2$ with boundary $\partial\Omega$ by $L^2(\Omega)$, and the Hilbert space of functions in $L^2(\Omega)$ with all partial derivatives up to a positive integer $m$ also in $L^2(\Omega)$ by $H^{{m}}(\Omega)$, see ^[28]. We endow $H^{{m}}(\Omega)$ with the usual norm and seminorm that we denote as $||{\, \cdot\, }||_{{m, \Omega}}$ and $|{\, \cdot\, }|_{{m, \Omega}}$, respectively.

The virtual element methods that we consider in the next sections are formulated on the mesh family $\big\{ \Omega_{h}\big\}_{h}$, where each mesh $\Omega_{h}$ is a partition of the computational domain $\Omega$ into nonoverlapping star-shaped polygonal elements $\text{E}$ with boundary $\partial \text{E}$, area $\left|{\text{E}}\right|$, center of gravity $\mathbf{x}_{\text{E}}$, and diameter $h_{\text{E}} = \sup_{ \mathbf{x}, \mathbf{y}\in \text{E}}\vert \mathbf{x}- \mathbf{y}\vert$. The mesh elements of $\Omega_{h}$ form a finite cover of $\Omega$ such that $\overline{\Omega} = \cup_{ \text{E}\in \Omega_{h}} \text{E}$ and the mesh size labeling each mesh $\Omega_{h}$ is defined by $h = \max_{ \text{E}\in \Omega_{h}} h_{\text{E}}$. A mesh edge $\text{e}$ has center $\mathbf{x}_{\text{e}}$ and length $h_{\text{e}}$; a mesh vertex $\text{v}$ has position vector $\mathbf{x}_{\text{v}}$. We denote the edges of the polygonal boundary $\partial \text{E}$ by ${ \mathcal{E}_{{ \text{E}}}}$.

For any integer number $\ell\geq0$, we let $\mathbb{P}_{{\ell}}(\text{E})$ and $\mathbb{P}_{{\ell}}(\text{e})$ denote the space of polynomials defined on the element $\text{E}$ and the edge $\text{e}$, respectively; $\mathbb{P}_{{\ell}}(\Omega_{h})$ denotes the space of piecewise polynomials of degree $\ell$ on the mesh $\Omega_{h}$. For the convenience of exposition, we also use the notation. $\mathbb{P}_{{-2}}(\text{E}) = \mathbb{P}_{{-1}}(\text{E}) = \{0\}$. Accordingly, it holds that ${{ q}}_{|{{ \text{E}}}}\in\mathbb{P}_{{\ell}}(\text{E})$ if $\text{E}\in \Omega_{h}$ for all $q\in\mathbb{P}_{{\ell}}(\Omega_{h})$. Throughout the paper, we use the multi-index notation, so that $\nu = (\nu_1, \nu_2)$ is a two-dimensional index defined by the two integer numbers $\nu_1, \nu_2\geq0$. Moreover, $D^{\nu} w = \partial^{|\nu|} w/{\partial x_1^{\nu_1}\partial x_2^{\nu_2}}$ denotes the partial derivative of order $|\nu| = \nu_1+\nu_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 $\nu = (0, 0)$. We also denote the partial derivatives of $w$ versus $x$ and $y$ by the shortcuts $\partial_{x} w$ and $\partial_{y} w$, and the normal and tangential derivatives with respect to a given edge by $\partial_{n} w$ and $\partial_{t} w$. Finally, we define the trace operator $\text{tr} \left({ A} \right) = \sum_{i} a_{ii}$ for any matrix (two-index tensor) $A = \big(a_{ij}\big)$. Given two tensors $A = (a_{ij})$ and $B = (b_{ij})$, we also define $A:B: = \sum_{i}\sum_{j}a_{ij}b_{ij}.$

2.2.   Governing equations in strong form

2.2.1.   Linear momentum balance equation

The linear momentum balance equation for the displacement $\boldsymbol{{u}}:\Omega\times(0, T]\to \mathbb{R}^2$ on the 2-D computational domain $\Omega$ reads as

where $\rho:\Omega\to \mathbb{R}$ is a scalar density field, $\ddot{ \boldsymbol{{u}}}$ denotes the second-order time derivative of $\boldsymbol{{u}}$, $\boldsymbol{{\sigma}}:\Omega\times(0, T]\to\mathbb{R}_{sym}^{{{2}}\times {{2}}}$ is the (time-dependent) stress tensor ($\mathbb{R}_{sym}^{{{2}}\times {{2}}}$ being the set of symmetric $2\times2$-sized real-valued tensors), and $\boldsymbol{{f}}:\Omega\times(0, T]\to \mathbb{R}^2$ is a suitable vector-valued forcing term. In this paper, we assume the density field is constant. The corresponding stress tensor $\boldsymbol{{\sigma}}$ is given by

where $g_{ \boldsymbol{{u}}}(\mathbf{x})\, : \mathbb{R}^2\to \mathbb{R}$ is a given function modeling different material behaviors such as hardening or softening effects, for instance, due to damage or fracture, and $\boldsymbol{{\varepsilon}}^{\text{vol}}$ and $\boldsymbol{{\varepsilon}}^{\text{dev}}$ are the volumetric and deviatoric strain tensors.

Tensors $\boldsymbol{{\varepsilon}}^{\text{vol}}$ and $\boldsymbol{{\varepsilon}}^{\text{dev}}$ are respectively given by

where $\boldsymbol{{I}}$ the $2\times2$ identity tensor, and the volumetric strain tensor can also be written as $\boldsymbol{{\varepsilon}}^{\text{vol}} = \big(\nabla\cdot \boldsymbol{{u}}\big)\boldsymbol{{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

Using (2.4)-(2.5) allows us to define the stress tensor $\sigma$ at $t = 0$ through definitions (2.2)-(2.3). To introduce the boundary conditions, we first assume that the domain boundary $\Gamma$ is decomposed into a Dirichlet boundary, e.g., $\Gamma_{D, \boldsymbol{{u}}}$ that is a closed subset of $\Gamma$ in the Euclidean topology, and Neumann boundary, e.g., $\Gamma_{N, \boldsymbol{{u}}}$, such that $\Gamma = \Gamma_{D, \boldsymbol{{u}}}\cup\Gamma_{N, \boldsymbol{{u}}}$ and $\big|\Gamma_{D, \boldsymbol{{u}}}\cap\Gamma_{N, \boldsymbol{{u}}} \big| = 0$. Then the boundary conditions read as:

Herein, $\boldsymbol{{u}}_0, \boldsymbol{{v}}_0, \boldsymbol{{u}}$, and $\boldsymbol{{g}}_N$ are given data; $\mathbf{n} = \big(\text{n}_x, \text{n}_y\big)^T$ is the unit vector orthogonal to $\Gamma$ and pointing out of $\Omega$.

2.2.2.   Fourth-order phase-field equation

The fourth-order equation for the phase field $d$ on the 2-D computational domain $\Omega$ reads as

where $\alpha_2, \alpha_1$, and $\alpha_0$ are model parameters, $g_{d}(d): \mathbb{R}\to \mathbb{R}$ is an affine function, and the history variable $\mathcal{H}_t:\Omega\to \mathbb{R}$ is an additional multiplicative term. We assume that $\alpha_2$ and $\alpha_1$ are strictly positive real numbers, and $\alpha_0$ is a non-negative real number. We remark that $g_{d}(d) : = \alpha_3 d + \alpha_4$, where $\alpha_3$ and $\alpha_4$ are constants. For simplicity, we let $\alpha_4 = 0$ in this work. We refer to ^[12] for a detailed discussion about the physical meaning of $g_{d}(d)$ and $\mathcal{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 $\Gamma = \Gamma_{D, d}\cup\Gamma_{N, d}$, where $\Gamma_{D, d}$ and $\Gamma_{N, d}$ represents the Dirichlet and Neumann parts of the domain boundary, respectively. Subsets $\Gamma_{D, d}$ and $\Gamma_{N, d}$ are disjoint in the sense that their intersection has zero 2D measure, i.e., $\big|\Gamma_{D, d}\cap\Gamma_{N, d}\big| = 0$. Additionally, we split the Dirichlet boundary as $\Gamma_{D, d} = \Gamma_{ D_1, d}\cup\Gamma_{ D_2, d}$ and the Neumann boundary as $\Gamma_{N, d}$ as $\Gamma_{N, d} = \Gamma_{ N_1, d}\cup\Gamma_{ N_2, d}$. We complete the mathematical formulation of the model problem (2.6) with the following boundary conditions

2.3.   Variational formulations

2.3.1.   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}}(\Omega)$, and have an assigned value on the boundary Dirichlet boundary $\Gamma_{D, \boldsymbol{{u}}}$:

Denoting $\boldsymbol{{V}}_{ \boldsymbol{{u}}, \boldsymbol{{0}}}$ if $\overline{{ \boldsymbol{{u}}}} = \boldsymbol{{0}}$ and testing Eq (2.1) with elements in $\boldsymbol{{V}}_{ \boldsymbol{{u}}, \boldsymbol{{0}}}$, we arrive at the variational formulation of the elastodynamics equation as

with the definitions:

and

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 = \sum_{i, j = 1, 2} a_{ij} b_{ij}$. As shown, for example, by Raviart and Thomas (see [13, Theorem 8-3.1]) the variational problem (2.8) is well-posed and its unique solution satisfies

assuming that the external load $\mathbf{f}\in L^2\big((0, T);[ L^2(\Omega)]^2\big)$, the boundary function $\mathbf{g}_N\in C^{{1}}\big((0, T);[H^{\frac{1}{2}}_{0, \Gamma_N}]^2\big)$, and the initial functions $\mathit{\boldsymbol{u}}_{0}\in[H^{1}_{0, \Gamma_D}(\Omega)]^2$, $\mathit{\boldsymbol{v}}_{0}\in[L^2(\Omega)]^2$.

2.3.2.   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 $\overline{{ {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.

Find $d\in V_{ {d}}$, such that

where

The wellposedness of the variational formulation, i.e., existence and uniqueness of ${d}\in H^2(\Omega)$ solving (2.10) with definitions (2.11)-(2.12), follows from applications of the Lax-Milgram theorem ^[14], since the bilinear form $\mathcal{A}_d(\cdot, \cdot)$ is coercive and continuous.

3.   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.

3.1.   VEM for the second-order momentum balance equation

Let $\boldsymbol{{V}}_{ \boldsymbol{{u}}, k}^h$ be the virtual element approximation of the affine functional space $\boldsymbol{{V}}_{ \boldsymbol{{u}}}$ defined in (2.7). The semi-discrete virtual element approximation of (2.8) reads as:

For all $t\in(0, T]$, find $\mathit{\boldsymbol{u}}_{h}(t)\in\boldsymbol{{V}}_{ \boldsymbol{{u}}, k}^h$, such that for $t = 0$ it holds that

Here, $\mathit{\boldsymbol{u}}_{h}(t)$ is the virtual element approximation of $\mathit{\boldsymbol{u}}$ and $\mathit{\boldsymbol{v}}_{h}$ the generic test function in $\mathbf{V}^{h}_{k}$, while $(\mathit{\boldsymbol{u}}_{0})_{ {I}}$ and $(\mathit{\boldsymbol{u}}_{1})_{ {I}}$ are the virtual element interpolants of the initial solution functions $\mathit{\boldsymbol{u}}(0)$ and $\dot{ \mathit{\boldsymbol{u}}}(0)$, respectively. The bilinear forms $\mathcal{M}^h_{ \boldsymbol{{u}}}(\cdot, \cdot)$, $\mathcal{A}^h_{ \boldsymbol{{u}}}(\cdot, \cdot)$, and the linear functional $\mathcal{F}^h_{ \boldsymbol{{u}}}(\cdot)$ are the virtual element approximations of $\mathcal{M}_{ \boldsymbol{{u}}}(\cdot, \cdot)$, $\mathcal{A}_{ \boldsymbol{{u}}}(\cdot, \cdot)$, and $\mathcal{F}_{ \boldsymbol{{u}}}(\cdot)$, in (2.8).

The total time interval $[0, T]$ is divided into $N_T$ subintervals with time step $\Delta t = \frac{T}{N_T}$ and time instants $t^{(n)} = n \Delta 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:

where $\beta$ in (0, 0.5] and $\gamma \in (0, 1]$ are Newmark parameters.

3.1.1.   Degrees of freedom, projection operators, local and global spaces

The vector-valued virtual element space is such that $\boldsymbol{{V}}_{ \boldsymbol{{u}}, k}^h(\text{E}) = V_{ \boldsymbol{{u}}, k}^h(\text{E}) \times V_{ \boldsymbol{{u}}, k}^h(\text{E})$, where $V_{ \boldsymbol{{u}}, k}^h(\text{E})$ denotes the local scalar conforming virtual element space (which will be introduced below). Accordingly, the vector-valued displacement field $\boldsymbol{{u}}_h$ is provided by the two components $\boldsymbol{{u}}_h = \left(u_{x, h}, u_{y, h} \right)$, where each component belongs to $V_{ \boldsymbol{{u}}, k}^h(\text{E})$. Both components $u_{i, h}$, for $i\in\{x, y\}$, are uniquely identified by the following degrees of freedom:

(U1) the values of $u_{i, h}$ at the vertices of $\text{E}$;

(U2) the edge polynomial moments of $u_{i, h}$ of order up to $k-2$ on each one-dimensional edge $e\in{ \mathcal{E}_{{ \text{E}}}}$:

(U3) the cell polynomial moments of $u_{i, h}$ of order up to $k-2$ on $\text{E}$:

In this definition, we use the symbols $M_{k-2}(\text{e})$ and $M_{k-2}(\text{E})$ to denote a basis for the polynomial spaces $\mathbb{P}_{{k}}(\text{e})$ and $\mathbb{P}_{{k-2}}(\text{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 $\Pi^{\nabla, \text{E}}_{ \boldsymbol{{u}}, {k}}: V_{ \boldsymbol{{u}}, k}^h(\text{E})\to\mathbb{P}_{{k}}(\text{E})$, which is defined as

where the second condition is needed to fix the kernel of the gradient operator. The local scalar virtual element space of order $k\ge 1$ is defined by

where $\mathbb{P}_k(\text{E}) \backslash \mathbb{P}_{k-2}(\text{E})$ denotes as the space spanned by the monomials of degree $k$ and $k-1$ on $\text{E}$, and

According to this definition, the orthogonal projection $\Pi^{0, \text{E}}_{ \boldsymbol{{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\geq1$, i.e., the finite-dimensional space ${\bf{V}}_{\boldsymbol{u}, k}^h$, is obtained by combining all the elemental spaces $\big[ V_{ \boldsymbol{{u}}, k}^h(\text{E})\big]^2$. Building upon the local spaces $V_{ \boldsymbol{{u}}, k}^h(\text{E})$, $k\ge1$ for all $\text{E}\in \Omega_{h}$, we define it as

3.1.2.   Virtual element bilinear forms and linear functional

In the virtual element setting, we define the bilinear forms $\mathcal{M}^h_{ \boldsymbol{{u}}}(\cdot, \cdot)$ and $\mathcal{A}^h_{ \boldsymbol{{u}}}(\cdot, \cdot)$ as the sum of elemental contributions, which are denoted by $\mathcal{M}^{ \text{E}}_{ \boldsymbol{{u}}}(\cdot, \cdot)$ and $\mathcal{A}^{ \text{E}}_{ \boldsymbol{{u}}}(\cdot, \cdot)$, respectively:

The local bilinear form $\mathcal{M}^{h, \text{E}}_{ \boldsymbol{{u}}}(\cdot, \cdot)$ is given by

where $S^{ \text{E}}_{{m}}(\cdot, \cdot)$ is the local stabilization term for $\mathcal{M}^{h, \text{E}}_{ \boldsymbol{{u}}}$. The bilinear form $\mathcal{M}^{ \text{E}}_{ \boldsymbol{{u}}}$ depends on the orthogonal projections $\Pi^{0}_{\boldsymbol{{u}}, {k}} \mathit{\boldsymbol{u}}_{h}$ and $\Pi^{0}_{\boldsymbol{{u}}, {k}} \mathit{\boldsymbol{v}}_{h}$, which are computable from the degrees of freedom of $\mathit{\boldsymbol{u}}_{h}$ and $\mathit{\boldsymbol{v}}_{h}$, respectively, see the previous section. The local bilinear form $\mathcal{A}^{h, \text{E}}_{ \boldsymbol{{u}}}(\cdot, \cdot)$ is given by

where $S^{ \text{E}}_{{a}}(\cdot, \cdot)$ is the local stabilization term for $\mathcal{A}^{ \text{E}}_{ \boldsymbol{{u}}}$. The bilinear form $\mathcal{A}^{ \text{E}}_{ \boldsymbol{{u}}}$ depends on the orthogonal projections $\Pi^{0, \text{E}}_{ \boldsymbol{{u}}, {k-1}} \boldsymbol{{\sigma}}(\mathit{\boldsymbol{u}}_{h})$ and $\Pi^{0, \text{E}}_{ \boldsymbol{{u}}, {k-1}} {\boldsymbol \varepsilon}(\mathit{\boldsymbol{v}}_{h})$, which are computable from the degrees of freedom of $\mathit{\boldsymbol{u}}_{h}$ and $\mathit{\boldsymbol{v}}_{h}$, respectively, see the previous section.

The local stabilization terms $S^{ \text{E}}_{{m}}(\cdot, \cdot), \, S^{ \text{E}}_{{a}}(\cdot, \cdot)\, :\, {\bf{V}}_{\boldsymbol{u}, k}^h\times {\bf{V}}_{\boldsymbol{u}, k}^h\to \mathbb{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 $\mathcal{F}^h_{ \boldsymbol{{u}}}(\cdot)$ is clearly computable since the edge trace $\mathit{\boldsymbol{v}}_{h}|_ \text{e}$ is a polynomial and $\Pi^{0, \text{E}}_{ \boldsymbol{{u}}, {k}}(\mathit{\boldsymbol{v}}_{h})$ is computable from the degrees of freedom of $\mathit{\boldsymbol{v}}_{h}$.

3.2.   VEM for the fourth-order phase-field equation

The virtual element method approximates the variational formulation (2.10) as follows

In this formulation, $V_{d, r}^h$ is the $H^{{2}}$-conforming approximation of the space $H^{{2}}(\Omega)$ provided by the VEM, $d_{h}$ and $c_{h}$ are the trial and test functions from this space.

3.2.1.   Degrees of freedom, projection operators, local and global spaces

On every mesh element $\text{E}\in \Omega_{h}$, the following set of values comprises the degrees of freedom of a virtual element function $c_{h}$ as shown in Figure 2:

($\mathbf{{D1}}$) for $r\geq2$, $c_{h}(\mathbf{x}_{\text{v}})$, $\partial_x c_{h}(\mathbf{x}_{\text{v}})$, $\partial_y c_{h}(\mathbf{x}_{\text{v}})$ for any vertex $\text{v}$ of $\partial \text{E}$;

($\mathbf{{D2}}$) for $r\geq4$, $\frac{1}{ \left|{\text{e}}\right|}\int_{ \text{e}} q c_{h} \, \text{d}S$ for any $q\in\mathbb{P}_{{ r-4}}(\text{e})$, and any edge $\text{e}\in{ \mathcal{E}_{{ \text{E}}}}$;

($\mathbf{{D3}}$) for $r\geq3$, $\int_{ \text{e}} q\partial_n c_{h} \, \text{d}S$ for any $q\in\mathbb{P}_{{ r-3}}(\text{e})$, and any edge $\text{e}\in{ \mathcal{E}_{{ \text{E}}}}$;

($\mathbf{{D4}}$) for $r\geq2$, $\frac{1}{ \left|{\text{E}}\right|}\int_{ \text{E}} q v \, \text{d}V$ for any $q\in\mathbb{P}_{{ r-2}}(\text{E})$.

Consider the integer $r\geq2$ and the bilinear form $\mathcal{B}(u, v) = \alpha_2 \mathcal{A}_{d, 2}(u, v)+\alpha_1 \mathcal{A}_{d, 1}(u, v)$ with $\alpha_2, \alpha_1 > 0$. We define the elliptic projection operator $\Pi^{ \mathcal{L}, \text{E}}_{d, { r}}:H^{{2}}(\text{E})\to\mathbb{P}_{{ r}}(\text{E})$ such that for every $v\in H^{{2}}(\text{E})$, the $r$-degree polynomial $\Pi^{ \mathcal{L}, \text{E}}_{d, { r}} v$ is the solution to the variational problem:

The elliptic projection operator $\Pi^{\mathcal{L}, \text{E}}_{ {d}, {r}}:H^{{2}}(\text{E})\cap C^{{1}}(\bar{ \text{E}})\to\mathbb{P}_{{r}}(\text{E})$ for any integer number $r\geq k$ is such that $\Pi^{\mathcal{L}, \text{E}}_{ {d}, {r}} {d}$ for all ${d}\in H^{{2}}(\text{E})\cap C^{{1}}(\bar{ \text{E}})$ is the the solution of the finite-dimensional variational problem

with the following additional conditions

and for $i = \left\{ x, y \right\}$:

For $r\ge 3$, the local virtual element space on element $\text{E}$ is defined by

where

For $r = 2$, we have the special case of the low-order virtual element space:

Note that $\mathbb{P}_{{r}}(\text{E})\subset V_{d, r}^h(\text{E})$. Building upon the local spaces $V_{d, r}^h(\text{E})$, $r\ge 2$, for all $\text{E}\in \Omega_{h}$ the global conforming virtual element space $V_{d, r}^h$ is defined as

This implies that $V_{d, r}^h \subset C^{{1}}(\Omega)$ because ${d}_h\in H^{{2}}(\text{E})$ in each element and $C^{{1}}$-regularity across the internal mesh faces.

3.2.2.   Virtual element bilinear forms and linear functional

Let $\text{E}\in \Omega_{h}$ be a mesh element, and consider the bilinear forms $\mathcal{A}^{ \text{E}}_{d, 2}, \mathcal{A}^{ \text{E}}_{d, 1}, \mathcal{A}^{ \text{E}}_{d, 0}: V_{d, r}^h(\text{E})\times V_{d, r}^h(\text{E})\to \mathbb{R}$ given by integrating on $\text{E}$ instead of $\Omega$ in the corresponding bilinear forms in (2.11)-(2.12). Let $\mathcal{A}^{ \text{E}}_{d}(\cdot, \cdot) = \alpha_2 \mathcal{A}^{ \text{E}}_{d, 2}(\cdot, \cdot)+\alpha_1 \mathcal{A}^{ \text{E}}_{d, 1}(\cdot, \cdot)+\alpha_0 \mathcal{A}^{ \text{E}}_{d, 0}(\cdot, \cdot)$. We use the elliptic projection $\Pi^{ \mathcal{L}, \text{E}}_{ {d}, { r}}$ and the $L^{{2}}$-orthogonal projection $\Pi^{0, \text{E}}_{ {d}, { r}}$ to define the virtual element bilinear form $\mathcal{A}_d^{ \text{E}}: V_{d, r}^h(\text{E})\times V_{d, r}^h(\text{E})\to \mathbb{R}$:

Herein, the stabilization term is also built by using the projection $\Pi^{ \mathcal{L}, \text{E}}_{ {d}, { r}}(c_{h})$, and the usual formula, so that the bilinear form $S^{\text{E}}_{h}: V_{d, r}^h(\text{E})\times V_{d, r}^h(\text{E})\to \mathbb{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.

4.   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.

4.1.   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_{ \boldsymbol{{u}}, k}^h$ and $H^2$-conforming $V_{d, r}^h$. For conciseness, we use the virtual element space $V_{ \boldsymbol{{u}}, k}^h$ as an example to introduce vector and matrix notations. Such notations can be easily extended to $V_{d, r}^h$.

We consider the following compact notation. For all element $\text{E}\in \Omega_{h}$, we locally number the degrees of freedom ($\mathbf{{U1}}$), ($\mathbf{{U2}}$), and ($\mathbf{{U3}}$) from $1$ to $N^\text{dofs}$. Then, we introduce the bounded, linear functionals $\text{dof}_i: V_{{ \boldsymbol{{u}}, k}}^{h}(\text{E})\rightarrow \mathbb{R}$, $i = 1, \dots, N^\text{dofs}$, such that

for $v_{h}\in V_{{ \boldsymbol{{u}}, k}}^{h}(\text{E})$. Let $\Lambda_{ \text{E}} = \big\{\text{dof}_{{i}}\left({\cdot} \right)\big\}_{i}$ denote the set of such functionals and collect the degrees of freedom of $v_{h}$ in the vector $\mathit{\boldsymbol{v}}_{h} = \big(\text{dof}_{{1}}\left({ v_{h}} \right), \ldots, \text{dof}_{{ N^\text{dofs}}}\left({ v_{h}} \right)\big)^T$. Since the degrees of freedom ($\mathbf{{U1}}$), ($\mathbf{{U2}}$), and ($\mathbf{{U3}}$) are unisolvent in $V_{{ \boldsymbol{{u}}, k}}^{h}(\text{E})$, the triplet $\big(\text{E}, V_{{ \boldsymbol{{u}}, k}}^{h}(\text{E}), \Lambda_{ \text{E}}\big)$ is a finite element in the sense of Ciarlet ^[16]. This property implies the existence of a Lagrangian basis set $\big\{\varphi^{ \boldsymbol{{u}}}_i\big\}_i$, with $\varphi^{ \boldsymbol{{u}}}_i\in V_{{ \boldsymbol{{u}}, k}}^{h}(\text{E})$, $i = 1, \dots, N^\text{dofs}$, which satisfies

We refer to the basis function set $\big\{\varphi^{ \boldsymbol{{u}}}_i\big\}_i$ as the "canonical" basis of $V_{{ \boldsymbol{{u}}, k}}^{h}(\text{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 $\mathbb{P}_{{ k}}(\text{E})\subset V_{{ \boldsymbol{{u}}, k}}^{h}(\text{E})$, which reads as

where $n_k$ is the cardinality of $\mathbb{P}_{{ k}}(\text{E})$. Since the polynomials $m^{ \boldsymbol{{u}}}_\alpha(\mathbf{x})$ are also virtual element functions, we can expand them on the canonical basis $\pmb{\varphi}$. We express such expansions as

where matrix $\pmb{D}_{ \boldsymbol{{u}}}$ has size $N^\text{dofs}\times n_k$ and collects all the expansion coefficients

so that

Following this notation, we also express the action of a differential operator $\mathcal{D}$, e.g., $\mathcal{D} = \Delta$ or $\mathcal{D} = \nabla$, in a entry-wise way, so that

and

Similarly, we express the action of the projectors $\Pi^{\nabla, \text{E}}_{\boldsymbol{{u}}, {k}}$, and $\Pi^{0, \text{E}}_{\boldsymbol{{u}}, {k}}$ on the canonical basis functions $\pmb{\varphi}_{ \boldsymbol{{u}}}$ and their expansion on the polynomial basis $\pmb{m}_{ \boldsymbol{{u}}}$ as follows:

where $\pmb{\Pi}^{\nabla, \text{E}}_{\boldsymbol{{u}}, {k}}$ and $\pmb{\Pi}^{0, \text{E}}_{\boldsymbol{{u}}, {k}}$ are the matrix representation of $\Pi^{\nabla, \text{E}}_{\boldsymbol{{u}}, {k}}$ and $\Pi^{0, \text{E}}_{\boldsymbol{{u}}, {k}}$, respectively. The expansion coefficients for the three projection operators applied to the basis function $\varphi_j$ are collected along the $j$-th column of the projection matrices $\pmb{\Pi}^{\nabla, \text{E}}_{\boldsymbol{{u}}, {k}}$, and $\pmb{\Pi}^{0, \text{E}}_{\boldsymbol{{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_{{d, r}}^{h}(\text{E})$, degree of freedoms ($\mathbf{{D1}}$)–($\mathbf{{D4}}$), canonical basis $\pmb{\varphi}_{d}$, the polynomial basis $\pmb{m}_d$, and more.

4.2.   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.

4.2.1.   Elliptic projector

The elliptic projector $\Pi^{\nabla, \text{E}}_{\boldsymbol{{u}}, {k}}$ is defined as, for $v_h \in V_{{ \boldsymbol{{u}}, k}}^{h}(\text{E})$,

Using vector and matrix notation as introduced in Section 4.1, we have

Substituting $\Pi^{\nabla, \text{E}}_{\boldsymbol{{u}}, {k}} \pmb{\varphi}_{ \boldsymbol{{u}}}^\textbf{T}$ with $\pmb{m}_{ \boldsymbol{{u}}}^\textbf{T} \pmb{\Pi}^{\nabla, \text{E}}_{ \boldsymbol{{u}}, k}$, we obtain

As a consequence, we have the following matrix equation

where

However, $\widetilde{ \pmb{G}}_{\boldsymbol{{u}}}$ is singular. We additionally define

Collecting Eqs (4.4) and (4.5), we obtain

It is worth noting that $\pmb{G}_{\boldsymbol{{u}}}$ is nonsingular by construction, and we can formally state that $\pmb{\Pi}^{\nabla, \text{E}}_{\boldsymbol{{u}}, { k}} = \pmb{G}_{\boldsymbol{{u}}}^{-1} \pmb{B}_{\boldsymbol{{u}}}$.

4.2.2.   $L^{{2}}$ orthogonal projectors

Recall the definition of the $L^2$-orthogonal projector $\Pi^{0, \text{E}}_{\boldsymbol{{u}}, {k}}$ acting on a scalar $v_h\in V_{{ \boldsymbol{{u}}, k}}^{h}(\text{E})$, is

Rewrite Eq (4.7) using vector and matrix notation as

We replace $\Pi^{0, \text{E}}_{\boldsymbol{{u}}, {k}} \pmb{\varphi}_{ \boldsymbol{{u}}}^\textbf{T}$ with $\pmb{m}_{ \boldsymbol{{u}}}^\textbf{T} \pmb{\Pi}^{0, \text{E}}_{ \boldsymbol{{u}}, k}$, and obtain

The above equation can be rewritten as

where

and

We remark that $\mathbf{H}_{ \boldsymbol{{u}}}$ is nonsingular, therefore $\pmb{\Pi}^{0, E}_{ \boldsymbol{{u}}, k} = \mathbf{H}_{ \boldsymbol{{u}}}^{-1} \mathbf{C}_{ \boldsymbol{{u}}}$.

In addition, the $L^2$-orthogonal projector acting on $\nabla v_h$, for $v_h \in V_{{ \boldsymbol{{u}}, k}}^{h}(\text{E})$ is defined as

We can choose $\boldsymbol{{q}} = \left[q^{\text{up}}, 0 \right]^T$ or $\boldsymbol{{q}} = \left[0, q^{\text{down}} \right]^T$ where $q^{\text{up}}, q^{\text{down}} \in \mathbb{P}_{{ k-1}}(\text{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 $\pmb{\Pi}^{0, y, \text{E}}_{ \boldsymbol{{u}}, k-1}$

where

It is easy to see that $\mathbf{H}_{ \boldsymbol{{u}}, k-1}$ is invertible, therefore,

4.2.3.   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.,

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, $\pmb{M}_ \text{E}^{ \boldsymbol{{u}}, \text{stab}} = \boldsymbol{{0}}$. Therefore, the local mass matrix is given by

Using the $L^2$-projector $\Pi^{0, \text{E}}_{\boldsymbol{{u}}, {k}}$ and adopting the vector and matrix notations, we have

where we have assumed the density of the material, $\rho$, 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 $\boldsymbol{{\varepsilon}}(\boldsymbol{{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 $\Pi^{0, x, \text{E}}_{\boldsymbol{{u}}, {k-1}}$, $\Pi^{0, y, \text{E}}_{\boldsymbol{{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 $\pmb{K}_ \text{E}^{ \boldsymbol{{u}}, \text{stab}}$ should scale with the consistency matrix $\pmb{K}_ \text{E}^{ \boldsymbol{{u}}, \text{cons}}$, which depends on the Lame parameters. Therefore, we let $\pmb{K}_ \text{E}^{ \boldsymbol{{u}}, \text{stab}}$ scale with $\max(2 \mu, \lambda)$. 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.

4.2.4.   Right-hand side approximation

Using the compact notation again, we rewrite the virtual element approximation of the first term of $\mathcal{F}^h_{ \boldsymbol{{u}}}$, see Eq (3.11), as

4.3.   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.

4.3.1.   Elliptic projectors

In this section, we first construct $\Pi^{\nabla, \text{E}}_{d, {r}}$, followed by $\Pi^{ \mathcal{L}, \text{E}}_{d, {r}}$. The elliptic projector $\Pi^{\nabla, \text{E}}_{d, {r}}$ is defined as, for $v_h \in V_{{d, k}}^{h}(\text{E})$

The above equation can be recast into the following matrix equation:

where

To fix the kernel of the differential operator $\nabla$, we additionally require

Gathering Eqs (4.29) and (4.30), we arrive at

Next, from the definition of $\Pi^{ \mathcal{L}, \text{E}}_{d, {r}}$,

we obtain

where

A (double) integration by parts of the right-hand-side of the definition of $\widetilde{ \pmb{B}}_d$ yields

Collecting Eqs (4.30) and (4.32), we have

It is worth noting that Eq (4.30) also fixes the kernel of $\mathcal{L}$, namely, $\pmb{G}^0_d \, \pmb{\Pi}^{ \mathcal{L}, \text{E}}_{d, { r}} = \pmb{B}^0_d$. Therefore, $\pmb{G}_{d, 1}$ and $\pmb{G}_d$ are nonsingular by construction, and we can formally state that

4.3.2.   $L^{{2}}$-orthogonal projector

Recall that $\Pi^{0, \text{E}}_{d, {r}}$ is defined as, for $v_h \in V_{{d, k}}^{h}(\text{E})$

The construction of the $L^{{2}}$ projection matrix $\pmb{\Pi}^{0, \text{E}}_{d, {r}}$ is similar to the derivation of Eq (4.10), therefore

where

4.3.3.   Local matrices

The construction of the local mass matrices uses the $L^{{2}}$-orthogonal projection, $\Pi^{0, \text{E}}_{d, {r}}$. Following a similar procedure to obtain Eq (4.22), we have

and

where we have assumed the linear function $g_d(d) = d$ and

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:

where

and

It is worth noting that we can also use other stabilizations as in ^[31].

5.   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 second-order 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.

Acknowledgments

The authors gratefully acknowledge the support of the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory under project number 20220129ER. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).

Conflict of interest

All authors declare no conflicts of interest in this paper.