Introduction of virtual element method (VEM)

Last updated on Sep 7, 2021

VEM is a recently proposed numerical discretization technique in Galerkin framework which is inspired by mimic finite difference (MFD) method, and since the date of birth it has been extensively developed and applied to a wide range of engineering problems. The core of the method is to construct a projector which projects the function in local virtual element space (or say, local shape function space) to a polynomial space with prescribed order. In this way VEM avoids the explicit construction of shape function and its integral over the element domain; and thus VEM is able to handle arbitrary polygonal mesh.

In order to give a simple and general illustration of how VEM works, here I take a linear elasticity boundary value problem in 2-dimensional domain

Ω

as example. The domain is partitioned into arbitrary polygonal meshes $T_{h}$ , and suppose one of the elements is pentagon as shown in Fig. 1.The weak form of the boundary value problem of this pentagon element K is given as

\begin{matrix} (1) & a^{K} (u^{K}, v^{K}) = \int_{Ω^{K}} σ (u^{K}) : ϵ (v) d Ω^{K} = f^{K} (v^{K}) = \int_{\partial Ω^{K}} v^{K} \cdot t d \partial Ω^{K}, \forall v^{K} \in V (K) \times V (K) \end{matrix}

where in Eq.(1)

a^{K}

and

f^{K}

are the bilinear and linear form defined in the element K;

σ

and

ϵ

are the stress and strain tensors, respectively, and

v

is the test vector-value function which belongs to the local Sobolev space

V (K) \times V (K)

defined within the element K;

u

is the solution of vector-value function to be found, typically the displacement vector in x- and y-directions;

Ω^{K}

is the element domain and

\partial Ω^{K}

is the corresponding Neumann boundary of element K.

FIG.1: DOFs of second-order VEM on pentagon element K

Generally the solution

u

is approximated by a function in the subspace of

V_{h} (K) \times V_{h} (K) \subseteq V (K) \times V (K)

called local virtual element space (or say, local shape function space), and the number of basis functions (shape functions) that span

V_{h} (K) \times V_{h} (K)

is equal to the number of degrees of freedom (DOFs) of element K defined in the following. Let the element in

V_{h} (K) \times V_{h} (K)

denote as

v_{h}

, the the vector-value function

b o l d s y m b o l v_{h}

in second-order VEM has the following properties:

$v_{h}$ is continuous on the boundary of element K
$v_{h}$ is a vector with second-order polynomial components on each edge of element K
$Δ v_{h}$ is a vector with second-order polynomial components on each edge of element K

where

Δ

is the Laplace operator. The corresponding DOFs of element K can be classified into the following three types

The values of $v_{h}$ at the 5 vertices of element K
The values of $v_{h}$ at the 5 midpoints of 5 edges of element K
The two moments of $v_{h}$ with respect to the constant vectors $p_{1} = [1, 0]^{T}$ and $p_{2} = [0, 1]^{T}$ in element K, that is

\begin{matrix} (2) & \frac{1}{| Ω^{K} |} \int_{Ω^{K}} v_{h} \cdot p_{1} d Ω^{K}; \frac{1}{| Ω^{K} |} \int_{Ω^{K}} v_{h} \cdot p_{2} d Ω^{K} \end{matrix}

where

| Ω^{K} |

is the area of element K. Based on the DOFs of types (i)-(iii), the total number of DOFs of element K for the second-order VEM is 10+10+2=22. Therefore, the function

v_{h}

is obtained by a linear combination of the basis functions in

V_{h} (K) \times V_{h} (K)

as follow

\begin{matrix} (3) & v_{h} = \sum_{i = 1}^{22} {d o f}_{i} (v_{h}) φ_{i} \end{matrix}

where

{d o f}_{i} (v_{h})

represents the ith DOF of

v_{h}

as given in types (i)-(iii) above, and

φ_{i}

is the ith vector-valued basis function of the local virtual element space

V_{h} (K) \times V_{h} (K)

of element K with the following two properties

Property 1: $φ_{i} (1, 2, . . ., 22)$ is a second-order polynomial on the element edge;
Property 2: $φ_{i} (1, 2, . . ., 22)$ satisfies the Kronecker-delta property, that is, for the jth DOF of element K we have ${d o f}_{j} (φ_{i}) = δ_{i j}, i, j = 1, 2, . . ., 22$

Similar to classical Galerkin finite element method (FEM), the entry of local stiffness matrix

k^{K}

of element K is calculated as

\begin{matrix} (4) & {(k^{K})}_{i j} = a^{K} (φ_{i}, φ_{j}) = \int_{Ω^{K}} σ (φ_{i}) : ϵ (φ_{j}) d Ω^{K}, f o r i, j = 1, 2, . . ., 22 \end{matrix}

where the subscripts ij indicate the location of entry in

k^{K}

. Now comes the most important part of VEM: define a projector $Π^{\nabla} : V_{h} (K) \times V_{h} (K) \to P^{2} \times P^{2}$ which maps the function $v_{h}$ in the space $V_{h} (K) \times V_{h} (K)$ onto the second-order polynomial space $P^{2} \times P^{2}$ satisfying the following orthogonality condition:

\begin{matrix} (5) & a^{K} (p_{α}, Π^{\nabla} v_{h} - v_{h}) = 0, f o r \forall p_{α} \in P_{2} \times P_{2} a n d v_{h} \in V_{h} \times V_{h} \end{matrix}

where

p_{α}

are the polynomial basis functions that span the second-order polynomial space

P_{2} \times P_{2}

. Then Eq.(4) is reformulated by using the projector

Π^{\nabla}

\begin{matrix} (6) & \begin{aligned} {(k^{K})}_{i j} & = a^{K} (Π^{\nabla} φ_{i} + (φ_{i} - Π^{\nabla} φ_{i}), Π^{\nabla} φ_{j} + (φ_{j} - Π^{\nabla} φ_{j})) \\ = \underset{c o n s i s t e n c y}{\underset{⏟}{a^{K} (Π^{\nabla} φ_{i}, Π^{\nabla} φ_{j})}} + \underset{s t a b i l i t y}{\underset{⏟}{a^{K} (φ_{i} - Π^{\nabla} φ_{i}, φ_{j} - Π^{\nabla} φ_{j})}} = {(k_{c}^{K})}_{i j} + {(k_{s}^{K})}_{i j} \end{aligned} \end{matrix}

where

Π^{\nabla} φ_{i} = \sum_{α = 1}^{12} S_{i, α} p_{α}, p_{α} \in P^{2} \times P^{2}

is the image of basis function

φ_{i}

projected on the second-order polynomial space which is a linear combination of polynomial basis functions

p_{α}

with coefficients

S_{i, α}

, and

{(k^{K})}_{i j}

can be thus considered as the sum of consistency term

{(k_{c}^{K})}_{i j}

and stability term

{(k_{s}^{K})}_{i j}

as shown in Eq.(6).

Based on the defination of projector

Π^{\nabla}

and linearity of the stain tensor, the component-wise consistency term

(k_{c}^{K})

can be obtained explicitly as

\begin{matrix} (7) & \begin{aligned} (k_{c}^{K}) & = \int_{Ω^{K}} σ (Π^{\nabla} φ_{i}) : ϵ (Π^{\nabla} φ_{j}) d Ω^{K} = \int_{Ω^{K}} σ (\sum_{α = 1}^{12} S_{i, α} p_{α}) : ϵ (\sum_{β = 1}^{12} S_{j, β} p_{β}) d Ω^{K} \\ = \sum_{α = 1}^{12} \sum_{β = 1}^{12} S_{i, α} S_{j, β} \int_{Ω^{K}} σ (p_{α}) : ϵ (p_{β}) d Ω^{K} == \sum_{α = 1}^{12} \sum_{β = 1}^{12} S_{i, α} S_{j, β} a^{K} (p_{α}, p_{β}) \\ = \sum_{α = 1}^{12} \sum_{β = 1}^{12} Π_{i α} Π_{β j} G_{α β} = {(Π^{T} G Π)}_{i j}, f o r i, j = 1, 2, . . ., 12 \end{aligned} \end{matrix}

where

Π

is the matrix representation of projector

Π^{\nabla}

, and matrix

G

is defined with entry

G_{α β} = a^{K} (p_{α}, p_{β})

.

As for the stability term

{(k_{s}^{K})}_{i j}

of element K, fistly define a 22-by-12 matrix D with entry as

\begin{matrix} (8) & D_{i α} = {d o f}_{i} (p_{α}) f o r i = 1, 2, . . ., 22 a n d α = 1, 2, . . ., 12 \end{matrix}

Next, by decomposing the image of projection

Π^{\nabla} φ_{i}

by the vector-valued functions in the local virtual element space

V_{h} \times V_{h}

, we can have the following matrix representation of another projection

{\overset{―}{Π}}^{\nabla} : V_{h} \times V_{h} \to V_{h} \times V_{h}

which maps the function in space

V_{h} \times V_{h}

onto itself

\begin{matrix} (9) & {\overset{―}{Π}}_{i j} = {d o f}_{i} (Π \nabla φ_{j}) = {d o f}_{i} (\sum_{β = 1}^{12} S_{j, β} p_{β}) = \sum_{β = 1}^{12} D_{i β} Π_{β j} = {(D Π)}_{i j} f o r i, j = 1, 2, . . ., 22 \end{matrix}

where

\overset{―}{Π}

is the matrix representation of projection

{\overset{―}{Π}}^{\nabla}

. Based on Eq.(9), the stability term

{(k_{s}^{K})}_{i j}

is calculated as

\begin{matrix} (10) & (k_{s}^{K}) = γ τ^{*} (I - {\overset{―}{Π}}_{i}^{T}) (I - {\overset{―}{Π}}_{j}) f o r i, j = 1, 2, . . ., 22 \end{matrix}

where

I

is the 22-by-22 identity matrix;

γ

is the user-defined parameter which can be chosen as 1 for elastic problem, and

τ *

can be calculated by various formulas, such as

τ^{*} = t r a c e (k_{c}^{K})

. Hence, through Eqs.(1) to (10), one can calculate all the entries in local stiffness matrix

k^{K}

of element K, and the global stiffness matrix can be obtained by assembling all the local stiffness matrices in the same manner as the classical Galerkin FEM.

Introduction of virtual element method (VEM)

Wei Shen

Phd student