# Edge-averaged virtual element methods for convection-diffusion and
convection-dominated problems111Submitted to the editors in 2024.

Shuhao Cao School of Science and Engineering, University of Missouri-Kansas
City, Kansas City, MO 64110<EMAIL_ADDRESS>Long Chen Department of
Mathematics, University of California at Irvine, Irvine, CA 92697
<EMAIL_ADDRESS>Seulip Lee Department of Mathematics, University of Georgia,
Athens, GA 30602<EMAIL_ADDRESS>
###### Abstract

This manuscript develops edge-averaged virtual element (EAVE) methodologies to
address convection diffusion problems effectively in the convection-dominated
regime. It introduces a variant of EAVE that ensures monotonicity (producing
an $M$-matrix) on Voronoi polygonal meshes, provided their duals are Delaunay
triangulations with acute angles. Furthermore, the study outlines a
comprehensive framework for EAVE methodologies, introducing another variant
that integrates with the stiffness matrix derived from the lowest-order
virtual element method for the Poisson equation. Numerical experiments confirm
the theoretical advantages of the monotonicity property and demonstrate an
optimal convergence rate across various mesh configurations. Keywords: virtual
elements; edge-averaged finite element; steady-state convection-diffusion
equation; monotone schemes; convection-dominated problems; polygonal meshes.
AMS Subject Classification: 65N30, 65N12

##  1 Introduction

We consider the convection-diffusion equation in a bounded domain
$\Omega\subset\mathbb{R}^{2}$ with the simply-connected polygonal Lipschitz
boundary $\partial\Omega$:

$\displaystyle-\nabla\cdot\left(\alpha(\mathbf{x})\nabla
u+\bm{\beta}(\mathbf{x})u\right)$ $\displaystyle=f$
$\displaystyle\quad\text{in }\Omega,$ (1.1a) $\displaystyle u$
$\displaystyle=0$ $\displaystyle\quad\text{on }\partial\Omega,$ (1.1b)

where $\alpha\in C^{0}(\bar{\Omega})$ with
$0<\alpha_{\text{min}}\leq\alpha(\mathbf{x})\leq\alpha_{\text{max}}$ for every
$\mathbf{x}\in\Omega$, $\bm{\beta}\in(C^{0}(\bar{\Omega}))^{2}$, and $f\in
L^{2}(\Omega)$. The convection-diffusion problem (1.1) expresses the
convective and molecular transport along a stream moving at the given velocity
$\bm{\beta}$ and with the diffusive effect from $\alpha$. The convection-
dominated regime means the situation that

$\alpha(\mathbf{x})\ll\left|\bm{\beta}(\mathbf{x})\right|,\quad\forall\mathbf{x}\in\Omega,$

and the convection-diffusion problem in this regime is important in describing
viscous fluid models with high Reynolds numbers governed by the Navier-Stokes
equations.

It has long been known that the standard-conforming finite element method
(FEM) fails to provide accurate numerical solutions for convection-dominated
problems. The finite element solutions based on standard Galerkin projections
contain spurious oscillations that deteriorate the solutions’ quality (e.g.,
see [31, 25, 40, 30]). Hence, various stabilized FEMs have been developed in
two primary directions. The first direction is based on presenting a provable
order of convergence in appropriate norms independent of small $\alpha$ by
adding stabilization terms to the weak formulation. This category includes the
streamline-upwind Petrov-Galerkin (SUPG) methods [17], continuous interior
penalty methods [19], and local projection stabilization (LPS) methods [32].
This stabilization technique mitigates such unexpected oscillations but does
not eliminate them. The second direction is the development of stabilized
methods satisfying the discrete maximum principle (DMP). When solving a
second-order elliptic equation numerically, it is known that a numerical
scheme satisfying the DMP guarantees discrete solutions without spurious
oscillations (see, e.g., [23, 18, 38]). Among many examples, the algebraic
flux corrections (AFC) method [4] is an upwind-based discretization that
directly satisfies the DMP with a reduced artificial diffusion effect, but it
requires solving a nonlinear system. The monotonicity property [27] is a
sufficient condition for the DMP, so monotone schemes have received attention
despite their first-order convergence when $\alpha=0$ [28]. Moreover, getting
a non-singular $M$-matrix in a discretization implies the monotonicity
property. Such an $M$-matrix is characterized by nonpositive off-diagonal
entries and diagonal dominance, including at least one strictly diagonally
dominant column. We note that this $M$-matrix condition makes it easier to
validate the DMP.

The edge-averaged finite element (EAFE) method [43, 33] is a well-known linear
monotone FEM, and it has been recently applied to space-time discretization
[3] and nearly inviscid incompressible flows [36]. The EAFE method is monotone
if the stiffness matrix of the linear conforming FEM for the Poisson equation
is an $M$-matrix, achieved by a Delaunay triangulation. The EAFE method has
been motivated by prior studies of the Scharfetter-Gummel method [39] (e.g.,
the finite volume Scharfetter-Gummel method [2], the exponential fitting
method [15, 16, 24], and the inverse-average-type finite element [37]). This
idea has recently been extended to a general framework for scalar and vector
convection-diffusion problems, called the simplex-averaged finite element
(SAFE) methods [42]. A stable mimetic finite difference (MFD) method [1] has
also been presented in the SAFE framework. However, none of these methods work
for general polygonal meshes because it may not be straightforward to
translate the conformity of the polynomial elements on simplices to polygons.
For this reason, our work proposes a generalization of the EAFE method
compatible with virtual element methods on polygonal meshes.

Virtual element methods (VEMs), a generalization of FEMs, are novel numerical
PDE methods on general polygonal and polyhedral meshes. The VEMs use a non-
polynomial approximation in a polygonal element, and exact pointwise values of
local basis functions are not needed inside the element. Instead the degrees
of freedom (DoFs) can be used to compute all necessary quantities to build a
stable and accurate discretization. More details can be found in [5, 6] and
Section 3 in this paper.

This paper specializes in developing stabilized virtual element methods for
convection dominated problems. The SUPG stabilization and local projection
stabilization (LPS) have been successfully applied to the VEMs in [11, 12, 13]
and [35], respectively. However, these approaches do not guarantee the DMP,
which may yield undesired oscillations in discrete solutions. Therefore, we
focus on generalizing the EAFE stabilization to the VEMs and call the
resulting stabilized methods edge-averaged virtual element (EAVE) methods. We
present two variants of the EAVE methods: (1) a monotone EAVE method on
Voronoi meshes with dual Delaunay and (2) a general EAVE framework working
with general polygonal meshes.

(1) The monotone EAVE method guarantees the DMP on Voronoi meshes with dual
Delaunay triangulations consisting of acute triangles. The basic idea of
deriving this method is utilizing flux approximations and dual edge patches
for mass lumping inspired by [43] and [14], respectively.

(2) The general EAVE framework presents another EAVE method on general
polygonal meshes. Its monotonicity property holds if the stiffness matrix of
the linear VEM for the Poisson equation is an $M$-matrix. To derive a bilinear
form in this framework, we define a flux approximation in the lowest-order
$H(\textbf{curl})$-conforming virtual element space and observe the
relationship between the lowest-order $H(\textbf{curl})$\- and linear
$H^{1}$-conforming spaces.

Through numerical experiments, we see that the monotone EAVE method produces a
more stable and effective solution than the general framework on the same
Voronoi meshes. The numerical results also demonstrate that all the EAVE
methods have first-order convergence on various polygonal meshes. Finally, we
compare the numerical performance of the EAVE methods with other stabilized
numerical methods and show the robustness and first order accuracy of EAVE.

The remaining sections are structured: Section 2 introduces necessary
notations and observes the EAFE stabilization’s main idea. The essential
definitions and theoretical aspects of the VEMs are introduced in Section 3.
In Section 4, we present a finite-volume bilinear form for the Poisson
equation compatible with a virtual element space and stabilization-free,
providing fundamental ideas of using dual edge patches for mass lumping.
Section 5 proposes a monotone EAVE method on Voronoi meshes with dual Delaunay
triangulations consisting of acute triangles. Section 6 presents a general
framework for the EAVE methods and demonstrates a sufficient condition for the
monotonicity property. In Section 7, our numerical experiments show the effect
of edge-averaged stabilization on VEMs as well as the optimal order of
convergence on different mesh types. Finally, we summarize our contribution in
this paper and discuss related research in Section 8.

##  2 Preliminaries

To begin with, we introduce some notation used throughout this paper. For a
bounded Lipschitz domain $\mathcal{D}\in\mathbb{R}^{2}$, we denote the Sobolev
space as $W^{m,p}(\mathcal{D})$ for $m\geq 0$ and $p\geq 1$. Its norm and
seminorm are denoted by $\|\cdot\|_{W^{k,p}(\mathcal{D})}$ and
$|\cdot|_{W^{k,p}(\mathcal{D})}$, respectively. When $p=2$, the Sobolev space
is written as $H^{m}(\mathcal{D})$, and the space $H^{0}(\mathcal{D})$
coincides with $L^{2}(\mathcal{D})$. The notation $H_{0}^{1}(\mathcal{D})$
means the space of $v\in H^{1}(\mathcal{D})$ such that $v=0$ on
$\partial\mathcal{D}$ in the trace sense. The polynomial spaces of degrees
less than or equal to $k$ are denoted as $\mathcal{P}_{k}(\mathcal{D})$, and
the space of continuous functions is denoted as $C^{0}(\mathcal{D})$.

###  2.1 Model problem

In the convection-diffusion equation (1.1a), we denote the flux as

$J(u):=\alpha\nabla u+\bm{\beta}u.$ (2.1)

The weak formulation of the problem (1.1) is to find $u\in H_{0}^{1}(\Omega)$
such that

$B(u,v)=F(v),\quad\forall v\in H_{0}^{1}(\Omega),$ (2.2)

where the bilinear form and the functional on the right-hand side are defined
as

$B(w,v):=\int_{\Omega}J(w)\cdot\nabla v\;{\rm d}\mathbf{x}\quad\text{and}\quad
F(v):=\int_{\Omega}fv\;{\rm d}\mathbf{x}.$

The well-posedness of the continuous problem (2.2) can be found in [27].

###  2.2 Exponential fitting of the flux

If there exists a potential function $\psi$ such that
$\nabla\psi=\alpha^{-1}\bm{\beta},$ then the flux $J(u)$ can be expressed as a
diffusion flux with a variable coefficient,

$J(u)=\kappa\nabla(e^{\psi}u)\quad\text{with}\quad\nabla\psi=\alpha^{-1}\bm{\beta},\quad\kappa({\bf
x})=\alpha({\bf x})e^{-\psi({\bf x})}.$

In terms of differential forms, $e^{\psi}u$ is a 0-form, $\nabla(e^{\psi}u)$
is a 1-form, $J(u)$ is a 2-form, and $\kappa$ is a Hodge star mapping a 1-form
to a 2-form.

If $\alpha$ and $\bm{\beta}$ are constant, the function $\psi$ is defined as
$\psi(\mathbf{x})=\alpha^{-1}\bm{\beta}\cdot\mathbf{x}$ and it is unique up to
a constant. However, such a potential function $\psi$ may not exist for
general $\alpha$ and $\bm{\beta}$ if
$\text{curl}\;(\alpha^{-1}\bm{\beta})\not=0$.

A proper flux approximation [43, 33] is defined in the lowest order Nédélec
space on triangular meshes, and it has been successfully applied to develop
the edge-averaged finite element (EAFE) scheme. In what follows, we briefly
explain the main idea of the EAFE scheme on a triangulation.

###  2.3 Flux approximations

Let $\mathcal{T}_{h}$ be a conforming triangulation for $\Omega$ and
$V_{h}=\mathcal{P}_{1}(\mathcal{T}_{h})$ be the piecewise linear
$H^{1}$-conforming finite element space. We consider the local flux
$J(u_{h})|_{T}=(\alpha\nabla u_{h}+\bm{\beta}u_{h})|_{T}$ for $u_{h}\in V_{h}$
and $T\in\mathcal{T}_{h}$. We define $\bm{\tau}_{E}$ as a scaled tangent
vector with $|\bm{\tau}_{E}|=|E|$. Suppose ${\bf x}(t)$ is a parametrization
of $E$ such that ${\bf x}(0)={\bf x}_{i}$ and ${\bf x}(1)={\bf x}_{j}$ are two
endpoints of $E$. Let us define

$\psi_{E}({\bf
x}(t))=\int_{0}^{t}\frac{1}{|E|}\alpha^{-1}(\bm{\beta}\cdot\bm{\tau}_{E})\;{\rm
d}s,\quad\quad\kappa_{E}=\alpha e^{-\psi_{E}}>0.$ (2.3)

Then, it follows

$\left.J(u_{h})\cdot\bm{\tau}_{E}\right|_{E}=\left.\kappa_{E}\nabla(e^{\psi_{E}}u_{h})\cdot\bm{\tau}_{E}\right|_{E}.$

That is, along any line segment, the flux in that direction is a diffusion
flux with variable coefficients $\kappa_{E}$. Note that $\psi_{E}$ is unique
up to a constant, and $\left.J(u_{h})\cdot\bm{\tau}_{E}\right|_{E}$ is
invariant to the constant change.

The local Nédélec space of the lowest order is defined as

$\mathcal{N}_{0}(T):=(\mathcal{P}_{0}(T))^{2}+\mathbf{x}^{\perp}\mathcal{P}_{0}(T),$

where $\mathbf{x}^{\perp}$ is a $90^{\circ}$ rotation of $\mathbf{x}$
satisfying $\mathbf{x}\cdot\mathbf{x}^{\perp}=0$. Its DoFs for any ${\bf
v}\in\mathcal{N}_{0}(T)$ are defined as

$\text{dof}_{E}(\mathbf{v})=\frac{1}{|E|}\int_{E}\mathbf{v}\cdot\bm{\tau}_{E}\;{\rm
d}s=\left.(\mathbf{v}\cdot\bm{\tau}_{E})\right|_{E}.$

Its canonical basis functions
$\\{\bm{\chi}_{E_{j}}\\}_{j=1}^{3}\subset\mathcal{N}_{0}(T)$ are obtained from
$\text{dof}_{E_{i}}(\bm{\chi}_{E_{j}})=\delta_{ij}$ for $1\leq i,j\leq 3$. If
we apply $\text{dof}_{E}$ to $\kappa_{E}^{-1}J(u_{h})|_{T}$, then we have

$\text{dof}_{E}(\kappa_{E}^{-1}J(u_{h}))=\text{dof}_{E}(\nabla(e^{\psi_{E}}u_{h}))=\delta_{E}(e^{\psi_{E}}u_{h}):=(e^{\psi_{E}}u_{h})({\bf
x}_{j})-(e^{\psi_{E}}u_{h})({\bf x}_{i}),$ (2.4)

where $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ are the endpoints of $E$, that is,
$\bm{\tau}_{E}={\bf x}_{j}-{\bf x}_{i}$. In the terminology of differential
forms, the flux $J(u_{h})$ is a $2$-form, and $\kappa_{E}^{-1}J(u_{h})$ is a
$1$-form. The functional $\text{dof}_{E}$ is applied to the $1$-form.

Define $\bar{\kappa}_{E}$ as the harmonic average of $\kappa_{E}=\alpha
e^{-\psi_{E}}$ on $E$, i.e.,

$\bar{\kappa}_{E}=\left(\frac{1}{|E|}\int_{E}\kappa_{E}^{-1}\;{\rm
d}s\right)^{-1}.$ (2.5)

We look for a flux approximation $J_{T}(u_{h})\in\mathcal{N}_{0}(T)$ written
as

$J_{T}(u_{h})=\sum_{E\subset\partial
T}\text{dof}_{E}(\kappa_{E}^{-1}J(u_{h}))\bar{\kappa}_{E}\bm{\chi}_{E}=\sum_{E\subset\partial
T}\delta_{E}(e^{\psi_{E}}u_{h})\bar{\kappa}_{E}\bm{\chi}_{E},$

where $\bar{\kappa}_{E}$ can be considered an edgewise Hodge star mapping a
1-form $\bm{\chi}_{E}$ to a 2-form. The flux approximation $J_{T}(u_{h})$ is
also invariant to the constant shift of $\psi_{E}$ to $\psi_{E}+c_{E}$.

Moreover, we can express $\nabla v_{h}$ in $\mathcal{N}_{0}(T)$ as

$\nabla v_{h}=\sum_{E\subset\partial T}(\nabla
v_{h}\cdot\bm{\tau}_{E})\bm{\chi}_{E}=\sum_{E\subset\partial
T}\delta_{E}(v_{h})\bm{\chi}_{E}.$

As a result, the following auxiliary bilinear form is used in the
discretization

$\displaystyle B^{T}_{\mathcal{N}}(u_{h},v_{h}):=$
$\displaystyle\int_{T}J_{T}(u_{h})\cdot\nabla v_{h}\;{\rm d}\mathbf{x}$
$\displaystyle=$ $\displaystyle\sum_{E_{i},E_{j}\subset\partial
T}\bar{\kappa}_{E_{i}}\delta_{E_{i}}(e^{\psi_{E_{i}}}u_{h})\delta_{E_{j}}(v_{h})\int_{T}\bm{\chi}_{E_{i}}\cdot\bm{\chi}_{E_{j}}\;{\rm
d}\mathbf{x}.$ (2.6)

The mass matrix $(\int_{T}\bm{\chi}_{E_{i}}\cdot\bm{\chi}_{E_{j}}\;{\rm
d}\mathbf{x})$ for $\mathcal{N}_{0}(T)$ is in general a full matrix. We shall
apply a mass lumping technique to make it diagonal.

###  2.4 Edge-averaged finite element scheme

To derive the EAFE scheme [43, 33], we approximate the bilinear form (2.6)
using mass lumping. We let
$\\{\lambda_{i}(\mathbf{x})\\}_{i=1}^{3}\subset\mathcal{P}_{1}(T)$ be the
barycentric coordinates to the vertices $\\{\mathbf{x}_{i}\\}_{i=1}^{3}$ in
$T\in\mathcal{T}_{h}$ and define

$\omega_{E}^{T}:=\frac{1}{2}\cot{\theta_{E}^{T}}=-\int_{T}\nabla\lambda_{i}\cdot\nabla\lambda_{j}\;{\rm
d}\mathbf{x}\quad\text{when}\quad\bm{\tau}_{E}=\mathbf{x}_{j}-\mathbf{x}_{i}.$

Then, the mass lumping for the local mass matrix in $\mathcal{N}_{0}(T)$ [29]
is

$\int_{T}\mathbf{v}\cdot\mathbf{w}\;{\rm
d}\mathbf{x}\approx\sum_{E\subset\partial
T}\omega_{E}^{T}\left.(\mathbf{v}\cdot\bm{\tau}_{E})\right|_{E}\left.(\mathbf{w}\cdot\bm{\tau}_{E})\right|_{E},\quad\forall{\bf
v},{\bf w}\in\mathcal{N}_{0}(T).$ (2.7)

That is, we use the diagonal matrix ${\rm
diag}(\omega_{E_{1}}^{T},\omega_{E_{2}}^{T},\omega_{E_{3}}^{T})$ to
approximate the $3\times 3$ dense matrix
$(\int_{T}\bm{\chi}_{E_{i}}\cdot\bm{\chi}_{E_{j}}\;{\rm d}\mathbf{x})$. The
mass lumping is exact when
$\mathbf{v},\mathbf{w}\in\nabla\mathcal{P}_{1}(T)\subset\mathcal{N}_{0}(T)$.
Consequently, we have

$\displaystyle B^{T}_{\mathcal{N}}(u_{h},v_{h})$
$\displaystyle=\int_{T}J_{T}(u_{h})\cdot\nabla v_{h}\;{\rm d}\mathbf{x}$
$\displaystyle\approx\sum_{E\subset\partial
T}\omega_{E}^{T}\bar{\kappa}_{E}\delta_{E}(e^{\psi_{E}}u_{h})\delta_{E}(v_{h})=:B_{h}^{T}(u_{h},v_{h}).$

The EAFE discrete problem is to find $u_{h}\in V_{h}$ such that

$\sum_{T\in\mathcal{T}_{h}}B_{h}^{T}(u_{h},v_{h})=F(v_{h}),\quad\forall
v_{h}\in V_{h}.$ (2.8)

In the EAFE bilinear form, it is clear to see $\bar{\kappa}_{E}>0$ and
$\delta_{E}(e^{\psi_{E}}\lambda_{i})\delta_{E}(\lambda_{j})=-e^{\psi_{E}(\mathbf{x}_{i})}<0$.
Therefore, the off-diagonal entries in the stiffness matrix of the EAFE
bilinear form in (2.8) are nonpositive if and only if
$\omega_{E}^{T}+\omega_{E}^{T^{\prime}}\geq 0$ for any interior edge
$E=\partial T\cap\partial T^{\prime}$, which is also equivalent to the
Delaunay condition of $\mathcal{T}_{h}$ in two dimensions. Hence, the
stiffness matrix corresponding to the EAFE bilinear form is an $M$-matrix if
and only if the stiffness matrix for the Poisson equation (computed by
$-\omega_{E}^{T}$) is an $M$-matrix.

In implementation, the coefficients $\alpha$ and $\bm{\beta}$ are approximated
as a constant $\alpha_{E}$ and a constant vector $\bm{\beta}_{E}$ on each edge
$E\subset\partial T$, respectively. For example,
$\alpha_{E}=(\alpha(\mathbf{x}_{i})+\alpha(\mathbf{x}_{j}))/2$ and
$\bm{\beta}_{E}=(\bm{\beta}(\mathbf{x}_{i})+\bm{\beta}(\mathbf{x}_{j}))/2$.
The edgewise potential function is given as
$\psi_{E}(\mathbf{x})={\alpha_{E}}^{-1}\bm{\beta}_{E}\cdot\mathbf{x}$. Thus,
the DoFs are explicitly expressed as

$\bar{\kappa}_{E}\delta_{E}(e^{\psi_{E}}u_{h})=\alpha_{E}\mathbb{B}({\alpha_{E}}^{-1}\bm{\beta}_{E}\cdot(\mathbf{x}_{i}-\mathbf{x}_{j}))u_{h}(\mathbf{x}_{j})-\alpha_{E}\mathbb{B}({\alpha_{E}}^{-1}\bm{\beta}_{E}\cdot(\mathbf{x}_{j}-\mathbf{x}_{i}))u_{h}({\bf
x}_{i}),$ (2.9)

where $\mathbb{B}(z)$ is the Bernoulli function,

$\mathbb{B}(z)=\left\\{\begin{array}[]{cl}\displaystyle\frac{z}{e^{z}-1}&z\not=0,\\\
1&z=0.\end{array}\right.$ (2.10)

##  3 Virtual element methods

This section introduces the essential definitions and fundamental theories of
the VEMs. Let $\mathcal{K}_{h}$ denote a decomposition of $\Omega$ into
polygonal elements $K$ with a mesh size $h=\max_{K\in\mathcal{K}_{h}}h_{K}$,
where $h_{K}$ is the diameter of $K$. In a polygonal element $K$,
$\mathbf{x}_{i}$ for $1\leq i\leq N_{V}$ denote vertices and $E_{j}$ for
$1\leq j\leq N_{E}$ denote edges, where $N_{V}$ is the number of vertices and
$N_{E}$ is the number of edges. For a simple polygon, $N_{V}=N_{E}$. We label
the indices of the vertices or edges counterclockwise. We denote scaled
tangent vectors on the edges by
$\bm{\tau}_{j}=\mathbf{x}_{j+1}-\mathbf{x}_{j}$ for $1\leq j\leq N_{E}-1$, and
$\bm{\tau}_{N_{E}}=\mathbf{x}_{1}-\mathbf{x}_{N_{V}}$.

###  3.1 The lowest order nodal and edge virtual spaces

We first introduce the lowest order local nodal space in [5, 6]. Let the
boundary space be

$B_{1}(\partial K):=\left\\{v_{h}\in C^{0}(\partial
K):v_{h}|_{E}\in\mathcal{P}_{1}(E),\ \forall E\subset\partial K\right\\}.$

Then, the local nodal space [5, 6] is defined as

$\mathcal{V}_{1}(K):=\left\\{v_{h}\in H^{1}(K):v_{h}|_{\partial K}\in
B_{1}(\partial K),\ \Delta v_{h}|_{K}=0\right\\},$

and its DoFs are given as the values of $v_{h}$ at the vertices of $K$, i.e.,

$\text{dof}_{i}(v_{h})=v_{h}(\mathbf{x}_{i}),\quad 1\leq i\leq N_{V}.$ (3.1)

The dimension of $\mathcal{V}_{1}(K)$ is $N_{V}$, and its canonical basis,
$\\{\phi_{j}\\}_{j=1}^{N_{V}}\subset\mathcal{V}_{1}(K)$, is defined by
$\text{dof}_{i}(\phi_{j})=\delta_{ij}$. It has been proved in [5] that the
DoFs (3.1) are unisolvent for $\mathcal{V}_{1}(K)$. We also highlight that
$\mathcal{P}_{1}(K)\subseteq\mathcal{V}_{1}(K)$, and the local
$H^{1}$-projection is computable using the DoFs (3.1), that is,
$\Pi^{\nabla}_{1}:\mathcal{V}_{1}(K)\rightarrow\mathcal{P}_{1}(K)$ with

$\displaystyle\int_{K}\left(\nabla\Pi_{1}^{\nabla}v_{h}-\nabla
v_{h}\right)\cdot\nabla p_{1}\;{\rm d}\mathbf{x}$
$\displaystyle=0,\quad\forall p_{1}\in\mathcal{P}_{1}(K),$
$\displaystyle\frac{1}{N_{V}}\sum_{i=1}^{N_{V}}\left(\Pi_{1}^{\nabla}v_{h}(\mathbf{x}_{i})-v_{h}(\mathbf{x}_{i})\right)$
$\displaystyle=0.$

The global nodal space is defined as

$\mathcal{V}_{1}(\mathcal{K}_{h}):=\left\\{v_{h}\in
H_{0}^{1}(\Omega):\left.v_{h}\right|_{K}\in\mathcal{V}_{1}(K),\ \forall
K\in\mathcal{K}_{h}\right\\},$

and its global DoFs are all the values of $v_{h}$ at the internal vertices of
$\mathcal{K}_{h}$.

We introduce the lowest order edge space with reduced DoFs [7, 9],

$\displaystyle\mathcal{N}_{0}(K)=\left\\{\mathbf{v}_{h}\in(L^{2}(K))^{2}\right.:$
$\displaystyle\text{div}\;\mathbf{v}_{h}=0,\
\text{rot}\;\mathbf{v}_{h}\in\mathcal{P}_{0}(K),$
$\displaystyle\left.\mathbf{v}_{h}\cdot\bm{\tau}_{E}\in\mathcal{P}_{0}(E),\
\forall E\subset\partial K\right\\},$

where $\text{rot}\;{\bf v}=\partial v_{2}/\partial x-\partial v_{1}/\partial
y$ for a vector ${\bf v}=\langle v_{1},v_{2}\rangle$. The DoFs for
$\mathcal{N}_{0}(K)$ [7, 8, 9, 10] use the lowest order edge moments, similar
to those of $\mathcal{N}_{0}(T)$,

$\text{dof}_{E}(\mathbf{v}_{h})=\frac{1}{|E|}\int_{E}\mathbf{v}_{h}\cdot\bm{\tau}_{E}\;{\rm
d}s=\left.(\mathbf{v}_{h}\cdot\bm{\tau}_{E})\right|_{E},$ (3.2)

for each edge $E\subset\partial K$, and the canonical basis,
$\\{\bm{\chi}_{j}\\}_{j=1}^{N_{E}}\subset\mathcal{N}_{0}(K)$, is defined by
$\text{dof}_{E_{i}}(\bm{\chi}_{j})=\delta_{ij}$ for $1\leq i\leq N_{E}$. Note
that, unlike the finite element spaces, explicit forms of
$\\{\bm{\chi}_{j}\\}_{j=1}^{N_{E}}\subset\mathcal{N}_{0}(K)$ are unknown.
Instead, only the DoFs are used to construct a discretization.

By integration by parts, $\text{rot}\;\mathbf{v}_{h}$ is determined by the
DoFs (3.2), so the dimension of $\mathcal{N}_{0}(K)$ is $N_{E}$. It is easy to
verify the unisolvence for $\mathcal{N}_{0}(K)$ because the divergence-free
condition implies
$\mathbf{v}_{h}=\textbf{rot}\;\vartheta=\langle\partial\vartheta/\partial
y,-\partial\vartheta/\partial x\rangle$ for some differentiable function
$\vartheta$, and integration by parts leads to

$\int_{K}|\mathbf{v}_{h}|^{2}\;{\rm
d}\mathbf{x}=\int_{K}\mathbf{v}_{h}\cdot\textbf{rot}\;\vartheta\;{\rm
d}\mathbf{x}=\int_{K}(\text{rot}\;\mathbf{v}_{h})\vartheta\;{\rm
d}\mathbf{x}-\sum_{E\subset\partial
K}\frac{1}{|E|}\int_{E}(\mathbf{v}_{h}\cdot\bm{\tau}_{E})\vartheta\;{\rm d}s.$

We also emphasize that $\mathcal{N}_{0}(T)\subset\mathcal{N}_{0}(K)$, and the
local $L^{2}$-projection
$\bm{\Pi}_{0}^{0}:\mathcal{N}_{0}(K)\rightarrow(\mathcal{P}_{0}(K))^{2}$ is
computable by the DoFs (3.2),

$\int_{K}\mathbf{v}_{h}\cdot\mathbf{p}_{0}\;{\rm
d}\mathbf{x}=\int_{K}\mathbf{v}_{h}\cdot\textbf{rot}\;p_{1}\;{\rm
d}\mathbf{x}=\int_{K}(\text{rot}\;\mathbf{v}_{h})p_{1}\;{\rm
d}\mathbf{x}-\sum_{E\subset\partial
K}\frac{1}{|E|}\int_{E}(\mathbf{v}_{h}\cdot\bm{\tau}_{E})p_{1}\;{\rm d}s,$

for some $p_{1}\in\mathcal{P}_{1}(K)$ with
$\textbf{rot}\;p_{1}=\mathbf{p}_{0}$.

The global edge virtual element space is required to preserve the
$H(\textbf{curl})$ conformity [8, 9, 10] with respect to the space

$H_{0}(\text{rot};\Omega)=\\{{\bf v}\in(L^{2}(\Omega))^{2}:\text{rot}\;{\bf
v}\in L^{2}(\Omega),\ {\bf v}\cdot\mathbf{t}=0\ \text{on}\ \partial\Omega\\},$

where $\mathbf{t}$ denotes the unit tangent vector. Hence, the two-dimensional
global space is defined as

$\mathcal{N}_{0}(\mathcal{K}_{h}):=\left\\{\mathbf{v}_{h}\in
H_{0}(\text{rot};\Omega):\left.\mathbf{v}_{h}\right|_{K}\in\mathcal{N}_{0}(K),\
\forall K\in\mathcal{K}_{h}\right\\},$

The global DoFs are $\mathbf{v}_{h}\cdot\bm{\tau}_{E}$ for all interior edges
of the decomposition $\mathcal{K}_{h}$.

We present an important relation between the local nodal space
$\mathcal{V}_{1}(K)$ and the local edge space $\mathcal{N}_{0}(K)$. The
detailed proof can be found in [10].

######  Proposition 3.1.

$\nabla\mathcal{V}_{1}(K)$ is a subset of $\mathcal{N}_{0}(K)$, and moreover

$\nabla\mathcal{V}_{1}(K)=\left\\{\mathbf{v}_{h}\in\mathcal{N}_{0}(K)\mid\textnormal{rot}\;\mathbf{v}_{h}=0\right\\}.$

In addition, we show the relation between the canonical basis functions
$\\{\phi_{i}\\}_{i=1}^{N_{V}}$ for $\mathcal{V}_{1}(K)$ and
$\left\\{\bm{\chi}_{j}\right\\}_{j=1}^{N_{E}}$ for $\mathcal{N}_{0}(K)$.

######  Lemma 3.2.

For the canonical basis $\bm{\chi}_{m},\bm{\chi}_{n}\in\mathcal{N}_{0}(K)$
with $1\leq m<n\leq N_{E}$, we have

$\bm{\chi}_{n}-\bm{\chi}_{m}=-\sum_{l=m+1}^{n}\nabla\phi_{l},$

where $\phi_{l}$ is the canonical basis of $\mathcal{V}_{1}(K)$.

###### Proof.

Let $\mathbf{w}_{h}=-\sum_{l=m+1}^{n}\nabla\phi_{l}$. Then,
$\mathbf{w}_{h}\in\nabla\mathcal{V}_{1}(K)\subset\mathcal{N}_{0}(K)$ by
Proposition 3.1. It is clear to see that
$\text{dof}_{E}(\mathbf{w}_{h})=\text{dof}_{E}(\bm{\chi}_{n}-\bm{\chi}_{m})$,
for all $E\subset\partial K$. Therefore,
$\mathbf{w}_{h}=\bm{\chi}_{n}-\bm{\chi}_{m}$ of $\mathcal{N}_{0}(K)$ by the
unisolvent condition. ∎

###  3.2 Bilinear forms for the Poisson equation

We introduce bilinear forms appearing in virtual element methods for the
Poisson equation. The local continuous bilinear form for the Poisson equation
is denoted as

$a^{K}(u,v)=\int_{K}\nabla u\cdot\nabla v\;{\rm d}\mathbf{x}.$ (3.3)

Since we do not know the explicit forms of $u_{h},v_{h}\in\mathcal{V}_{1}(K)$
inside $K$, it would not be possible to compute $a^{K}(u_{h},v_{h})$ exactly.
Therefore, we use the equivalent bilinear form [5, 6],

$a_{h}^{K}(u_{h},v_{h})=\int_{K}\nabla\Pi_{1}^{\nabla}u_{h}\cdot\nabla\Pi_{1}^{\nabla}v_{h}\;{\rm
d}\mathbf{x}+S^{K}(u_{h}-\Pi_{1}^{\nabla}u_{h},v_{h}-\Pi_{1}^{\nabla}v_{h})$
(3.4)

for all $u_{h},v_{h}\in\mathcal{V}_{1}(K)$. In this case, the stabilization
term $S^{K}$ needs to satisfy that there are positive constants $c_{0}$ and
$c_{1}$ independent of $h_{K}$ such that

$c_{0}a^{K}(v_{h},v_{h})\leq S^{K}(v_{h},v_{h})\leq c_{1}a^{K}(v_{h},v_{h}),$

for every $v_{h}\in\mathcal{V}_{1}(K)$ with $\Pi_{1}^{\nabla}v_{h}=0$. These
inequalities directly imply the norm equivalence,

$\gamma_{*}a^{K}(v_{h},v_{h})\leq
a_{h}^{K}(v_{h},v_{h})\leq\gamma^{*}a^{K}(v_{h},v_{h}),\quad\forall
v_{h}\in\mathcal{V}_{1}(K),$ (3.5)

for some positive constants $\gamma_{*}$ and $\gamma^{*}$ independent of
$h_{K}$ (see [5, 6, 22] for details). Here, we introduce the most popular
choice for $S^{K}$,

$S_{\mathcal{V}}^{K}(u_{h},v_{h}):=\sum_{i=1}^{N_{V}}\text{dof}_{i}(u_{h})\text{dof}_{i}(v_{h}).$
(3.6)

Another choice of the stabilization term is the $H^{\frac{1}{2}}$-inner
product [20],

$\displaystyle S^{K}_{\mathcal{E}}(u_{h},v_{h}):=\sum_{E\subset\partial
K}\delta_{E}(u_{h})\delta_{E}(v_{h}),$ (3.7)

where $\delta_{E}$ is defined in (2.4). The verification of norm equivalence
(3.5) for (3.6) can be found in [22] and for (3.7) in [20], respectively.

In Section 4, we will present a discretization of VEMs based on the finite
volume formulation requiring no stabilization.

###  3.3 Construction of the right-hand side

We use the approximation technique [5] for the right-hand side. Let
$f_{0}^{K}$ be a piecewise constant approximation of $f$. Then, we define

$F_{h}(v_{h}):=\sum_{K\in\mathcal{K}_{h}}\int_{K}f_{0}^{K}\bar{v}_{h}\;{\rm
d}\mathbf{x}=\sum_{K\in\mathcal{K}_{h}}|K|f_{0}^{K}\bar{v}_{h},\quad\bar{v}_{h}:=\frac{1}{N_{V}}\sum_{i=1}^{N_{V}}v_{h}(\mathbf{x}_{i}).$

This approximation yields an optimal estimate [5],

$F_{h}(v_{h})-F(v_{h})\leq
Ch\left(\sum_{K\in\mathcal{T}_{h}}\left|f\right|_{H^{1}(K)}^{2}\right)^{1/2}\left|v_{h}\right|_{H^{1}(\Omega)}.$
(3.8)

##  4 A finite volume bilinear form with virtual element

In this section, we present a finite volume bilinear form for the Poisson
equation corresponding to $a^{K}(\cdot,\cdot)$ in (3.3) on Voronoi meshes with
respect to dual Delaunay triangulations consisting of acute triangles. We
employ piecewise constant approximations on dual edge patches inspired by [14]
and show that the finite volume bilinear form is compatible with the virtual
element space $\mathcal{V}_{1}(\mathcal{K}_{h})$ and stabilization-free.

###  4.1 Voronoi meshes with dual Delaunay triangulations

In two dimensions, a Delaunay triangulation consists of triangles whose
circumcircles include no other vertex of the triangulation, which means that
the sum of the two angles opposite to any internal edge is less than or equal
to $\pi$. On the other hand, we call a partition a Voronoi tessellation if
each seed point has a corresponding region of all vertices closer to that seed
than any other seed using the Euclidean distance. The Delaunay and Voronoi
partitions are dual in a convex domain because there are one-to-one
relationships between vertices on one mesh and elements on the other. Indeed,
each seed point of a Voronoi element is a vertex of a Delaunay triangle, while
each circumcenter of a Delaunay triangle is a vertex of a Voronoi polygon.
Moreover, the duality between Delaunay and Voronoi means that the edges on a
mesh are orthogonal to the ones on the other. This dual configuration plays an
essential role in finite volume methods [26].

In this section, we assume all the triangles determined by interior nodes are
acute; each triangle includes its circumcenter. Then, we generate a Voronoi
mesh by connecting the neighboring circumcenters; each Voronoi element
contains a vertex of a Delaunay triangle. We consider the Voronoi mesh primary
for cell-centered finite volume methods, while the Delaunay triangulation is
dual (see the left figure in Figure 1). This primary mesh is considered a sub-
class of general Voronoi meshes whose dual Delaunay may include obtuse
triangles.

On the boundary of the domain, we obtain Voronoi nodes by extending
perpendicular bisectors of triangles’ edges around the boundary, where
additional corner Voronoi nodes may help to express the boundary more
accurately. Boundary Delaunay nodes are located between two Voronoi nodes
along the boundary, satisfying the orthogonality of the triangles’ edges and
the boundary. Hence, the control volumes (or triangles) in the dual Delaunay
triangulation can cover the whole domain. However, since considering the zero
boundary condition (1.1b), we exclude the control volumes containing the
boundary Deluanay nodes and employ the one-to-one relationships between a
control volume and a primary (Voronoi) vertex (see the right figure in Figure
1).

Figure 1: A primary Voronoi mesh with a dual Delaunay triangulation consisting
of acute triangles; black squares are Voronoi nodes, and blue circles are seed
points (or Delaunay nodes).

We let $\mathcal{K}_{h}$ denote a primary Voronoi mesh with a dual Delaunay
triangulation $\mathcal{T}_{h}=\mathcal{K}_{h}^{*}$ consisting of acute
triangles with interior nodes. Each polygon $K\in\mathcal{K}_{h}$ contains
only one Delaunay vertex $\mathbf{x}_{K}^{*}$ in the interior while each
interior Voronoi vertex is the circumcenter of a triangle
$T\in\mathcal{T}_{h}$. Then, an edge $E$ of $K$ and $\mathbf{x}_{K}^{*}$
determine a sub-triangle in $K$. Any polygon $K$ is split into sub-triangles,
and the number of sub-triangles is the number of edges of $K$. Every internal
edge in $\mathcal{K}_{h}$ is shared by two such sub-triangles. We let $D_{E}$
denote the union of the two sub-triangles, and $E^{*}$ denote the dual edge in
$\mathcal{T}_{h}$ included in $D_{E}$ with $E^{*}\perp E$. For a boundary edge
$E\subset\partial\Omega$, $D_{E}$ is defined as one triangle appearing in the
domain $\Omega$, and $E^{*}$ is defined as the line segment from ${\bf
x}_{K}^{*}$ to $E$. Then, all the dual edge patches $D_{E}$ form a polygonal
mesh of $\Omega$ and will be denoted by $\mathcal{D}_{h}$ (see Figure 1 for an
illustration). Moreover, $\mathcal{E}_{h}^{o}$ denotes the collection of all
interior edges in $\mathcal{K}_{h}$. We also note that if ${\bf n}_{E^{*}}$ is
a unit normal vector on $E^{*}$ and $\bm{\tau}_{E}$ is a scaled tangent vector
on $E$, then ${\bf n}_{E^{*}}=\bm{\tau}_{E}/|E|$ because $E^{*}\perp E$. By
the construction, their intersection $\mathbf{m}_{E^{*}}=E^{*}\cap E$ is the
midpoint of $E^{*}$ but may not be the midpoint of $E$.

Finally, we recall the virtual element spaces on the primary Voronoi mesh and
define function spaces on the dual triangulation and patches.

  * •

$\mathcal{V}_{1}(\mathcal{K}_{h})\subset H_{0}^{1}(\Omega)$: the linear nodal
virtual element space on the polygons $K\in\mathcal{K}_{h}$.

  * •

$\mathcal{N}_{0}(\mathcal{K}_{h})\subset H_{0}(\text{rot};\Omega)$: the lowest
order edge virtual element space on the polygons $K\in\mathcal{K}_{h}$.

  * •

$V_{0}(\mathcal{T}_{h})$: the set of piecewise constant functions $v_{h}^{*}$
on the triangles $T\in\mathcal{T}_{h}$, excluding the control volumes
containing the boundary Delaunay nodes.

  * •

$W_{0}(\mathcal{D}_{h})$: the set of piecewise constant functions $w_{h}$ on
the edge patches $D_{E}\in\mathcal{D}_{h}$, where $w_{h}|_{D_{E}}\equiv 0$ if
$E\subset\partial\Omega$.

###  4.2 A finite volume bilinear form on Voronoi meshes

We define a one-to-one map between $\mathcal{V}_{1}(\mathcal{K}_{h})$ and
$V_{0}(\mathcal{T}_{h})$ by the DoFs. In each interior triangle
$T\in\mathcal{T}_{h}$, there is only one interior Voronoi vertex ${\bf x}_{T}$
because of the duality between $\mathcal{K}_{h}$ and $\mathcal{T}_{h}$, so the
DoFs of $v_{h}\in\mathcal{V}_{1}(\mathcal{K}_{h})$ at the vertex, $v_{h}({\bf
x}_{T})$, maps to $v_{h}^{*}|_{T}$ of $v_{h}^{*}\in V_{0}(\mathcal{T}_{h})$.
Through this one-to-one mapping, we can extend the definition of $\delta_{E}$
to functions in $V_{0}(\mathcal{T}_{h})$, i.e.,
$\delta_{E}(v_{h}^{*})=\delta_{E}(v_{h})$ is the difference of function values
of $v_{h}$ on the vertices of $E$ following the orientation of $E$.

We consider a balance equation for the Poisson equation over each
$T\in\mathcal{T}_{h}$,

$-\int_{\partial{T}}\nabla u\cdot{\bf n}\;{\rm d}s=\int_{T}f\;{\rm d}{\bf x}.$
(4.1)

We define a normal flux approximation $\mathcal{G}_{{\bf n}}(u_{h})\in
W_{0}(\mathcal{D}_{h})$ as

$\left.\mathcal{G}_{{\bf n}}(u_{h})\right|_{D_{E}}:=\nabla
u_{h}(\mathbf{m}_{E^{*}})\cdot{\bf n}_{E^{*}}=\nabla
u_{h}(\mathbf{m}_{E^{*}})\cdot\frac{\bm{\tau}_{E}}{|E|}=\frac{\delta_{E}(u_{h})}{|E|}.$
(4.2)

A cell-centered finite volume method uses the discrete normal flux on
$E^{*}\subset D_{E}$

$\int_{E^{*}}\nabla u_{h}\cdot{\bf n}_{E^{*}}\;{\rm
d}s\approx\int_{E^{*}}\mathcal{G}_{{\bf n}}(u_{h})\;{\rm
d}s=\frac{|E^{*}|}{|E|}\delta_{E}(u_{h}).$

For $u_{h}$ in the virtual element space $\mathcal{V}_{1}(\mathcal{K}_{h})$,
$\nabla u_{h}\cdot{\bf n}_{E^{*}}$ is not computable but can be evaluated at
the midpoint $\mathbf{m}_{E^{*}}$. By the duality, it equals to the tangential
derivative along $E$, which is computable as $u_{h}|_{E}$ is linear. Then, if
$\phi_{T}^{*}\in V_{0}(\mathcal{T}_{h})$ is the characteristic function of a
triangle $T\in\mathcal{T}_{h}$, i.e. $\phi_{T}^{*}\equiv 1$ in $T$ and
$\phi_{T}^{*}\equiv 0$ in $\Omega\setminus T$, the balance equation (4.1)
implies a discrete problem

$\sum_{E\in\mathcal{E}_{h}^{o}}\frac{|E^{*}|}{|E|}\delta_{E}(u_{h})\delta_{E}(\phi_{T}^{*})=\int_{T}f\;{\rm
d}{\bf x}=\int_{\Omega}f\phi_{T}^{*}\;{\rm d}{\bf x}=F(\phi_{T}^{*}).$

By using linear combinations in $V_{0}(\mathcal{T}_{h})$, we define a local
bilinear form

$\mathbf{a}_{h}^{K}(u_{h},v_{h}^{*}):=\sum_{E\subset\partial
K}\omega_{E}^{K}\delta_{E}(u_{h})\delta_{E}(v_{h}^{*}),\quad\forall
v_{h}^{*}\in V_{0}(\mathcal{T}_{h}),$

where $\omega_{E}^{K}:=|E^{*}|/(2|E|)$. Therefore, the corresponding discrete
problem is to find $u_{h}\in\mathcal{V}_{1}(\mathcal{K}_{h})$ such that

$\sum_{K\in\mathcal{K}_{h}}\mathbf{a}_{h}^{K}(u_{h},v_{h}^{*})=F(v_{h}^{*}),\quad\forall
v_{h}^{*}\in V_{0}(\mathcal{T}_{h}),$ (4.3)

which is a Petrov-Galerkin formulation of the Poisson equation. By changing
$\delta_{E}(v_{h}^{*})$ to $\delta_{E}(v_{h})$, we obtain a Galerkin
formulation $\mathbf{a}_{h}^{K}(u_{h},v_{h})$.

In the virtual element framework, the bilinear form can be explained in terms
of mass lumping in the edge virtual element space $\mathcal{N}_{0}(K)$ for
$K\in\mathcal{K}_{h}$ (see also (2.7)),

$a^{K}(u_{h},v_{h})=\int_{K}\nabla u_{h}\cdot\nabla v_{h}\;{\rm d}{\bf
x}\approx\sum_{E\subset\partial K}\omega_{E}^{K}\left.(\nabla
u_{h}\cdot\bm{\tau}_{E})\right|_{E}\left.(\nabla
v_{h}\cdot\bm{\tau}_{E})\right|_{E}=\mathbf{a}_{h}^{K}(u_{h},v_{h})$ (4.4)

because $\left.(\nabla v_{h}\cdot\bm{\tau}_{E})\right|_{E}=\delta_{E}(v_{h})$.
Moreover, we emphasize that the DoFs of
$v_{h}\in\mathcal{V}_{1}(\mathcal{K}_{h})$ compute the bilinear form
$\mathbf{a}_{h}^{K}(\cdot,\cdot)$ without any stabilization term, so the
discretization (4.3) can be viewed as a stabilization-free virtual element
method.

For the right-hand side, we apply the mass lumping with the vertices of a
triangle $T\in\mathcal{T}_{h}$ by using the constant approximation $f_{0}^{K}$
on $K$ introduced in Section 3.3. More precisely, when the three vertices of
$T_{j}\in\mathcal{T}_{h}$ are in polygons $K_{1}$, $K_{2}$, and $K_{3}$ in
$\mathcal{K}_{h}$, the right hand side is calculated by

$F(\phi_{j}^{*})=\int_{\Omega}f\phi_{j}^{*}\;{\rm d}{\bf
x}\approx\frac{|T_{j}|}{3}\left(f_{0}^{K_{1}}+f_{0}^{K_{2}}+f_{0}^{K_{3}}\right)=:\mathcal{F}_{h}(\phi_{j}^{*}).$
(4.5)

Moreover, if the polygons $K_{1}$, $K_{2}$, and $K_{3}$ are regular, then
$\mathcal{F}_{h}(\phi_{j}^{*})=F_{h}(\phi_{j})$ in Section 3.3.

###  4.3 $M$-matrix property

We stress that the stiffness matrix corresponding to ${\bf
a}_{h}^{K}(\cdot,\cdot)$ is an $M$-matrix through the following lemma.

######  Lemma 4.1.

If $\mathcal{K}_{h}$ is a primary Voronoi mesh with a dual Delaunay
triangulation consisting of acute triangles, then the stiffness matrix
corresponding to ${\bf a}_{h}^{K}(\cdot,\cdot)$ is an $M$-matrix.

###### Proof.

Let $\zeta_{j}$ be a global nodal virtual function in
$\mathcal{V}_{1}(\mathcal{K}_{h})$ such that
$\zeta_{j}(\mathbf{x}_{i})=\delta_{ij}$ for any internal vertex
$\mathbf{x}_{i}$ in $\mathcal{K}_{h}$. We consider the stiffness matrix whose
entries are

$\mathbf{A}_{ij}=\sum_{K\in\mathcal{K}_{h}}{\bf
a}_{h}^{K}(\zeta_{j},\zeta_{i}).$

If $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ ($i\not=j$) are the endpoints of an
edge $E\in\mathcal{E}_{h}^{o}$, we have

${\bf a}_{h}^{K}(\zeta_{j},\zeta_{i})=-\omega_{E}^{K}=-|E^{*}|/(2|E|)<0.$

The partition of unity implies that

$\sum_{i}\mathbf{A}_{ij}=\sum_{K\in\mathcal{K}_{h}}\mathbf{a}_{h}^{K}(\zeta_{j},1)=0$

if $\mathbf{x}_{j}$ has no neighbor on the boundary. If $\mathbf{x}_{j}$ has a
neighbor on the boundary, then $\sum_{i}\mathbf{A}_{ij}>0$. Hence, the
stiffness matrix from the bilinear form $\mathbf{a}_{h}^{K}(\cdot,\cdot)$ is
an invertible $M$-matrix. ∎

###  4.4 Convergence analysis for Poisson equation

We use an alternative form with the normal flux approximation (4.2) on the
edge patches $D_{E}$ for convergence analysis,

$\sum_{K\in\mathcal{K}_{h}}\mathbf{a}_{h}^{K}(u_{h},v_{h})=\sum_{E\in\mathcal{E}_{h}^{o}}\frac{|E^{*}|}{|E|}\delta_{E}(u_{h})\delta_{E}(v_{h})=2\sum_{E\in\mathcal{E}_{h}^{o}}\int_{D_{E}}\mathcal{G}_{{\bf
n}}(u_{h})\mathcal{G}_{{\bf n}}(v_{h})\;{\rm d}\mathbf{x}$ (4.6)

for $u_{h},v_{h}\in\mathcal{V}_{1}(\mathcal{K}_{h})$ because
$|E^{*}||E|=2|D_{E}|$.

######  Theorem 4.2.

Assume that $\mathcal{K}_{h}$ is a primary Voronoi mesh with a dual Delaunay
triangulation consisting of acute triangles, $\nabla u\cdot{\bf n}_{E^{*}}\in
H^{1}(D_{E})$ for all $E\in\mathcal{E}_{h}^{o}$, and $u\in
C^{0}(\bar{\Omega})$. Then, we have for
$u_{h}\in\mathcal{V}_{1}(\mathcal{K}_{h})$ and $\mathcal{G}_{\bf n}(u_{h})$ in
(4.2),

$\left(\sum_{E\in\mathcal{E}_{h}^{o}}\lVert\nabla
u\cdot{\mathbf{n}}_{E^{*}}-\mathcal{G}_{{\bf
n}}(u_{h})\rVert^{2}_{L^{2}(D_{E})}\right)^{1/2}\leq
Ch\left(\sum_{E\in\mathcal{E}_{h}^{o}}|\nabla
u\cdot{\mathbf{n}}_{E^{*}}|^{2}_{H^{1}(D_{E})}\right)^{1/2}.$

###### Proof.

We adopt the idea of the analysis in [14]. We define two different
approximations of the exact normal flux $\nabla u\cdot{\bf n}_{E^{*}}$. The
first one, $\bar{\mathcal{G}}_{{\bf n}}(u)\in W_{0}(\mathcal{D}_{h})$, is the
average of the flux on $E^{*}\subset D_{E}$,

$\left.\bar{\mathcal{G}}_{{\bf
n}}(u)\right|_{D_{E}}:=\frac{1}{|E^{*}|}\int_{E^{*}}\nabla u\cdot{\bf
n}_{E^{*}}\;{\rm d}s$

for each $D_{E}\in\mathcal{D}_{h}$. Then, it follows from the balance equation
(4.1), the definitions of $\bar{\mathcal{G}}_{{\bf n}}(u)$ and
$\mathcal{G}_{{\bf n}}(v_{h})$, and the fact $|E^{*}||E|=2|D_{E}|$ that

$\displaystyle
2\sum_{E\in\mathcal{E}_{h}^{o}}\int_{D_{E}}\bar{\mathcal{G}}_{{\bf
n}}(u)\mathcal{G}_{{\bf n}}(v_{h})\;{\rm
d}\mathbf{x}=\sum_{E\in\mathcal{E}_{h}^{o}}|E^{*}|(\bar{\mathcal{G}}_{{\bf
n}}(u)|_{D_{E}})\delta_{E}(v_{h}^{*})=F(v_{h}^{*}).$

Hence, if we compare this equation with the discrete problem (4.3) applying
(4.6), then we have the orthogonality

$\sum_{E\in\mathcal{E}_{h}^{o}}\int_{D_{E}}(\mathcal{G}_{{\bf
n}}(u_{h})-\bar{\mathcal{G}}_{{\bf n}}(u))\mathcal{G}_{{\bf n}}(v_{h})\;{\rm
d}\mathbf{x}=0,\quad\forall v_{h}\in\mathcal{V}_{1}(\mathcal{K}_{h}).$

The second approximation, $u_{I}\in\mathcal{V}_{1}(\mathcal{K}_{h})$, is the
nodal value interpolant of $u$ onto the virtual element space. Then, we obtain

$\displaystyle\sum_{E\in\mathcal{E}_{h}^{o}}\lVert\mathcal{G}_{{\bf
n}}(u_{h})-\bar{\mathcal{G}}_{{\bf n}}(u)\rVert$
${}^{2}_{L^{2}(D_{E})}=\sum_{E\in\mathcal{E}_{h}^{o}}\int_{D_{E}}(\mathcal{G}_{{\bf
n}}(u_{h})-\bar{\mathcal{G}}_{{\bf n}}(u))(\mathcal{G}_{{\bf
n}}(u_{I})-\bar{\mathcal{G}}_{{\bf n}}(u))\;{\rm d}{\bf x}$
$\displaystyle\leq\sum_{E\in\mathcal{E}_{h}^{o}}\lVert\mathcal{G}_{{\bf
n}}(u_{h})-\bar{\mathcal{G}}_{{\bf
n}}(u)\rVert_{L^{2}(D_{E})}\lVert\mathcal{G}_{{\bf
n}}(u_{I})-\bar{\mathcal{G}}_{{\bf n}}(u)\rVert_{L^{2}(D_{E})}.$

Moreover, the Cauchy-Schwarz inequality implies that

$\sum_{E\in\mathcal{E}_{h}^{o}}\lVert\mathcal{G}_{{\bf
n}}(u_{h})-\bar{\mathcal{G}}_{{\bf
n}}(u)\rVert^{2}_{L^{2}(D_{E})}\leq\sum_{E\in\mathcal{E}_{h}^{o}}\lVert\mathcal{G}_{{\bf
n}}(u_{I})-\bar{\mathcal{G}}_{{\bf n}}(u)\rVert^{2}_{L^{2}(D_{E})}.$

Using the triangle inequality, we have

$\lVert\mathcal{G}_{{\bf n}}(u_{I})-\bar{\mathcal{G}}_{{\bf
n}}(u)\rVert_{L^{2}(D_{E})}\leq\lVert\nabla u\cdot{\bf
n}_{E^{*}}-\mathcal{G}_{{\bf n}}(u_{I})\rVert_{L^{2}(D_{E})}+\lVert\nabla
u\cdot{\bf n}_{E^{*}}-\bar{\mathcal{G}}_{{\bf n}}(u)\rVert_{L^{2}(D_{E})}.$

Each term means the first order approximation because $\mathcal{G}_{{\bf
n}}(u_{I})$ equals $\nabla u\cdot{\bf n}_{E^{*}}$ at a point on $E\subset
D_{E}$, and $\bar{\mathcal{G}}_{{\bf n}}(u)$ equals $\nabla u\cdot{\bf
n}_{E^{*}}$ at a point on $E^{*}\subset D_{E}$. ∎

##  5 A monotone virtual element method

We present a monotone virtual element scheme by utilizing the finite volume
approach in Section 4.

###  5.1 A monotone edge-averaged virtual element method

Let $\mathcal{K}_{h}$ be a primary Voronoi mesh with a dual Delaunay
triangulation $\mathcal{T}_{h}=\mathcal{K}_{h}^{*}$ consisting of acute
triangles (see Figure 1). We consider the same one-to-one map between
$\mathcal{V}_{1}(\mathcal{K}_{h})$ and $V_{0}(\mathcal{T}_{h})$, and the
constant space $W_{0}(\mathcal{D}_{h})$ on the edge patches
$D_{E}\in\mathcal{D}_{h}$ as in Section 4. For
$J_{K}(u_{h})\in\mathcal{N}_{0}(K)$, $J_{K}(u_{h})\cdot{\bf
n}_{E^{*}}=(J_{K}(u_{h})\cdot\bm{\tau}_{E})/|E|$ on each interior edge
$E\in\mathcal{E}_{h}^{o}$. Thus, we define a flux approximation
$\mathcal{J}_{{\bf n}}(u_{h})\in W_{0}(\mathcal{D}_{h})$ in light of (2.4),

$\left.\mathcal{J}_{{\bf
n}}(u_{h})\right|_{D_{E}}:=J_{K}(u_{h})(\mathbf{m}_{E^{*}})\cdot{\bf
n}_{E^{*}}=J_{K}(u_{h})(\mathbf{m}_{E^{*}})\cdot\frac{\bm{\tau}_{E}}{|E|}=\frac{\bar{\kappa}_{E}\delta_{E}(e^{\psi_{E}}u_{h})}{|E|}.$
(5.1)

The normal flux $J_{K}(u_{h})\cdot{\bf n}_{E^{*}}$ is computable at the
intersection point $\mathbf{m}_{E^{*}}$. In a balance equation corresponding
to the convection-diffusion equation (1.1a) over each $T\in\mathcal{T}_{h}$,

$-\int_{\partial T}J(u)\cdot{\bf n}\;{\rm d}s=\int_{T}f\;{\rm d}{\bf x},$
(5.2)

the discrete flux $\mathcal{J}_{{\bf n}}(u_{h})$ approximates the flux
$J(u_{h})\cdot{\bf n}_{E^{*}}$ on each $E^{*}\subset D_{E}$,

$\int_{E^{*}}J(u_{h})\cdot{\bf n}_{E^{*}}\;{\rm
d}s\approx\int_{E^{*}}\mathcal{J}_{{\bf n}}(u_{h})\;{\rm
d}s=\frac{|E^{*}|}{|E|}\bar{\kappa}_{E}\delta_{E}(e^{\psi_{E}}u_{h}).$

By introducing $\phi_{T}^{*}\in V_{0}(\mathcal{T}_{h})$ such that
$\phi_{T}^{*}\equiv 1$ in $T$ and $\phi_{T}^{*}\equiv 0$ in $\Omega\setminus
T$, a cell-centered finite volume method for the balance equation (5.2) is
represented as

$\sum_{E\in\mathcal{E}_{h}^{o}}\frac{|E^{*}|}{|E|}\bar{\kappa}_{E}\delta_{E}(e^{\psi_{E}}u_{h})\delta_{E}(\phi_{T}^{*})=\int_{T}f\;{\rm
d}{\bf x}=F(\phi_{T}^{*}).$ (5.3)

Thus, using linear combinations of such $\phi_{T}^{*}$ in
$V_{0}(\mathcal{T}_{h})$, we define a local bilinear form on
$K\in\mathcal{K}_{h}$,

$\mathcal{B}_{h}^{K}(u_{h},v_{h}^{*}):=\sum_{E\subset\partial
K}\omega_{E}^{K}\bar{\kappa}_{E}\delta_{E}(e^{\psi_{E}}u_{h})\delta_{E}(v_{h}^{*}),$

where $\omega_{E}^{K}=|E^{*}|/(2|E|)$. Therefore, the corresponding discrete
problem is to seek $u_{h}\in\mathcal{V}_{1}(\mathcal{K}_{h})$ such that

$\sum_{K\in\mathcal{K}_{h}}\mathcal{B}_{h}^{K}(u_{h},v_{h}^{*})=F(v_{h}^{*}),\quad\forall
v_{h}^{*}\in V_{0}(\mathcal{T}_{h}),$ (5.4)

which is a Petrov-Galerkin discretization for the convection-diffusion
equation (1.1a).

We can also apply the mass lumping to derive a Galerkin formulation

$\int_{K}J_{K}(u_{h})\cdot\nabla v_{h}\;{\rm d}{\bf x}\approx\sum_{E\subset
K}\omega_{E}^{K}\left.(J_{K}(u_{h})\cdot\bm{\tau}_{E})\right|_{E}\left.(\nabla
v_{h}\cdot\bm{\tau}_{E})\right|_{E}=:\mathcal{B}_{h}^{K}(u_{h},v_{h}).$

This approach can be considered as a cell-centered finite volume Scharfetter-
Gummel method. One can find a vertex-centered finite volume Scharfetter-Gummel
method [2] when the primary mesh is a Delaunay triangulation and the dual one
is a Voronoi mesh. The vertex-centered method is identical to the EAFE method
[33, 43]. Also, a stable mimetic finite difference (MFD) method [1] is
considered as a vertex-centered method.

######  Corollary 5.1.

Provided that a polygonal mesh $\mathcal{K}_{h}$ is a Voronoi mesh with a dual
Delaunay triangulation consisting of acute triangles. For $\alpha\in
C^{0}(\bar{\Omega})$ and $\bm{\beta}\in(C^{0}(\bar{\Omega}))^{2}$, the EAVE
method (5.4) is monotone, and consequently (5.4) is well-posed.

###### Proof.

The proof is similar to that of Lemma 4.1, so we omit the proof. ∎

##  6 General framework for EAVE methods

This section presents a general framework for edge-averaged virtual element
(EAVE) methods. In this framework, we derive another EAVE bilinear form
working with general polygonal meshes. The monotone EAVE method in Section 5
is a particular case of this general framework by applying the mass lumping in
$\mathcal{N}_{0}(K)$ on special polygonal meshes.

###  6.1 EAVE method on general polygonal meshes

We begin with the bilinear form with the flux approximation
$J_{K}(u_{h})\in\mathcal{N}_{0}(K)$ in light of (2.4),

$\int_{K}J_{K}(u_{h})\cdot\nabla v_{h}\;{\rm
d}\mathbf{x}=\sum_{E\subset\partial
K}\bar{\kappa}_{E}\delta_{E}(e^{\psi_{E}}u_{h})\int_{K}\bm{\chi}_{E}\cdot\nabla
v_{h}\;{\rm d}\mathbf{x}.$ (6.1)

We then apply a mass lumping technique and Lemma 3.2 (see [34] for details).
Let $E_{ij}$ be the line segment joining $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$
for $i<j$ in $K$. Then, the EAVE bilinear form in a general polygon $K$ is
defined as

$B_{h}^{K}(u_{h},v_{h}):=\sum_{1\leq i<j\leq
N_{V}}\omega_{ij}^{K}\bar{\kappa}_{E_{ij}}\delta_{ij}^{K}(e^{\psi_{E_{ij}}}u_{h})\delta_{ij}^{K}(v_{h}),$
(6.2)

where $\omega_{ij}^{K}=-a_{h}^{K}(\phi_{i},\phi_{j})$ in (3.4),
$\bar{\kappa}_{E_{ij}}$ is the harmonic average of $\kappa_{E_{ij}}$, and
$\delta_{ij}^{K}(v_{h})=v_{h}({\bf x}_{j})-v_{h}({\bf x}_{i})$. Thus, the
corresponding discrete problem is to find
$u_{h}\in\mathcal{V}_{1}(\mathcal{K}_{h})$ such that

$\sum_{K\in\mathcal{K}_{h}}B_{h}^{K}(u_{h},v_{h})=F_{h}(v_{h}),\quad\forall
v_{h}\in\mathcal{V}_{1}(\mathcal{K}_{h}).$ (6.3)

######  Remark 6.1.

We note that $\omega_{ij}^{K}=-a_{h}^{K}(\phi_{i},\phi_{j})$ in (6.2) is
computable and already given in the virtual element methods for the Poisson
equation. If we choose
$\omega_{ij}^{K}=-\mathbf{a}_{h}^{K}(\phi_{i},\phi_{j})$ as in (4.3) on a
Voronoi mesh corresponding to a Delaunay triangulation consisting of acute
triangles, the general bilinear form (6.2) becomes the one in the monotone
EAVE method (5.4). In practice, constant approximations $\alpha_{E_{ij}}$ and
$\bm{\beta}_{E_{ij}}$ on $E_{ij}$ give a function
$\psi_{E_{ij}}(\mathbf{x})={\alpha_{E_{ij}}}^{-1}\bm{\beta}_{E_{ij}}\cdot\mathbf{x}$
locally defined on $E_{ij}$. Therefore, the bilinear form is computed as

$\displaystyle\bar{\kappa}_{E_{ij}}\delta_{E_{ij}}^{K}(e^{\psi_{E_{ij}}}u_{h})$
$\displaystyle=\alpha_{E_{ij}}\mathbb{B}({\alpha_{E_{ij}}}^{-1}\bm{\beta}_{E_{ij}}\cdot(\mathbf{x}_{i}-\mathbf{x}_{j}))u_{h}(\mathbf{x}_{j})$
$\displaystyle\qquad\qquad-\alpha_{E_{ij}}\mathbb{B}({\alpha_{E_{ij}}}^{-1}\bm{\beta}_{E_{ij}}\cdot(\mathbf{x}_{j}-\mathbf{x}_{i}))u_{h}(\mathbf{x}_{i}),$

where $\mathbb{B}(\cdot)$ is the Bernoulli function defined in (2.10).

###  6.2 A sufficient condition for monotonicity property

In addition, the following lemma presents a sufficient condition for the
monotonicity property of the general framework of the EAVE methods.

######  Lemma 6.2.

If the stiffness matrix for the Poisson equation with the virtual bilinear
form (3.4) is an $M$-matrix, then the stiffness matrix from the EAVE method
(6.3) is an $M$-matrix for $\alpha\in C^{0}(\bar{\Omega})$ and
$\bm{\beta}\in\left(C^{0}(\bar{\Omega})\right)^{2}$.

###### Proof.

It suffices to show that the off-diagonal components of the stiffness matrix
are nonpositive. For $i\not=j$, it is clear to see that $\mathbf{B}_{ij}=0$ if
${\bf x}_{i}$ and ${\bf x}_{j}$ are not the vertices of a polygon. Otherwise,
we let $\mathcal{K}_{ij}$ be the set of the polygons having both
$\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ as vertices. Then, each component
becomes

$\mathbf{B}_{ij}=\sum_{K\in\mathcal{K}_{ij}}B_{h}^{K}(\zeta_{j},\zeta_{i}).$

For any polygon $K\in\mathcal{K}_{ij}$, we have

$B_{h}^{K}(\zeta_{j},\zeta_{i})=-\omega_{ij}^{K}\bar{\kappa}_{E_{ij}}e^{\psi_{E_{ij}}(\mathbf{x}_{j})}\quad\text{or}\quad-\omega_{ji}^{K}\bar{\kappa}_{E_{ji}}e^{\psi_{E_{ji}}(\mathbf{x}_{j})}.$

Note that $B_{h}^{K}(\zeta_{j},\zeta_{i})$ depends on the local indices of $i$
and $j$ which are determined counterclockwise in $K$. Hence, since
$\omega_{ij}^{K}=-a^{K}_{h}(\zeta_{i},\zeta_{j})\geq 0$ and
$\bar{\kappa}_{E_{ij}}>0$, we get $\mathbf{B}_{ij}\leq 0$. ∎

##  7 Numerical experiments

In this section, we present numerical experiments to demonstrate the effect of
the monotonicity property and the optimal order of convergence on different
mesh types. First, in the virtual element space
$\mathcal{V}_{1}(\mathcal{K}_{h})$, we compute a mesh-dependent norm defined
by the bilinear form (3.4), and we denote it as follows:

$\left\|v_{h}\right\|_{A}:=\left(\sum_{K\in\mathcal{K}_{h}}a_{h}^{K}(v_{h},v_{h})\right)^{1/2},\quad\forall
v_{h}\in\mathcal{V}_{1}(\mathcal{K}_{h}).$ (7.1)

This norm, which we call the $A$-norm, is equivalent to the $H^{1}$-norm from
the property of the bilinear form in (3.5). We also compute the mesh-dependent
$L^{\infty}$-norm,

$\left\|v_{h}\right\|_{\infty}:=\max_{i}\left(\text{dof}_{i}(|v_{h}|)\right),\quad\forall
v_{h}\in\mathcal{V}_{1}(\mathcal{K}_{h}),$

to check the presence of spurious oscillations on numerical solutions. The
numerical experiments are implemented by authors’ codes developed based on
$i$FEM [21]. Furthermore, we consider the following example problem:

######  Example 7.1.

Let the computational domain be $\Omega=(0,1)\times(0,1)$. The convection
dominated problem in [25] is governed by,

$\displaystyle-\nabla\cdot\left(\epsilon\nabla u+\bm{\beta}u\right)$
$\displaystyle=f$ $\displaystyle\quad\text{in }\Omega,$ (7.2a) $\displaystyle
u$ $\displaystyle=g$ $\displaystyle\quad\text{on }\partial\Omega,$ (7.2b)

where $\epsilon\ll 1$ is the diffusion coefficient, $\bm{\beta}=\langle
0,-1\rangle$, and $f\equiv 0$. The exact solution to the problem is given as

$u(x,y)=x\left(\frac{1-e^{(y-1)/\epsilon}}{1-e^{-2/\epsilon}}\right),$ (7.3)

so it has sharper boundary layer around $y=1$ when $\epsilon$ gets smaller.

###  7.1 Monotonicity test on Voronoi meshes with dual Delaunay

In this test, we compare the edge-averaged virtual element methods (5.4) and
(6.3) on Voronoi meshes with dual Delaunay triangulations consisting of acute
triangles:

  * •

M-EAVE: Monotone edge-averaged virtual element method (5.4);

  * •

EAVE: Edge-averaged virtual element method (6.3).

hexa-dual |  | voro-dual
---|---|---
|  |
Figure 2: Voronoi meshes corresponding to dual Delaunay triangulations with
$h=2^{-3}$.

We use hexagonal and unstructured Voronoi meshes with dual Delaunay
triangulations consisting of acute triangles. In Figure 2, the blue dots mean
vertices in a Delaunay triangulation, and it is easy to see that each triangle
with the blue dots contains an internal Voronoi vertex.

Table 1: A mesh refinement study for EAVE and M-EAVE with $\epsilon=10^{-2}$ ($e_{h}=u_{I}-u_{h}$). | hexa-dual
---|---
| EAVE | M-EAVE
$h$ |  $\lVert e_{h}\rVert_{A}$ | Order |  $\lVert e_{h}\rVert_{\infty}$ | Order |  $\lVert e_{h}\rVert_{A}$ | Order |  $\lVert e_{h}\rVert_{\infty}$ | Order
$2^{-3}$ | 1.631e-1 | - | 4.549e-2 | - | 1.728e-1 | - | 4.358e-2 | -
$2^{-4}$ | 1.058e-1 | 0.62 | 2.113e-2 | 1.11 | 9.922e-2 | 0.80 | 1.912e-2 | 1.19
$2^{-5}$ | 5.062e-2 | 1.06 | 7.153e-3 | 1.56 | 4.265e-2 | 1.22 | 5.871e-3 | 1.70
$2^{-6}$ | 1.737e-2 | 1.54 | 1.649e-3 | 2.12 | 1.430e-2 | 1.58 | 1.295e-3 | 2.18
| voro-dual
| EAVE | M-EAVE
$h$ | $\lVert e_{h}\rVert_{A}$ | Order | $\lVert e_{h}\rVert_{\infty}$ | Order | $\lVert e_{h}\rVert_{A}$ | Order | $\lVert e_{h}\rVert_{\infty}$ | Order
$2^{-2}$ | 1.127e-1 | - | 3.291e-2 | - | 2.425e-1 | - | 7.086e-2 | -
$2^{-3}$ | 9.554e-2 | 0.24 | 1.829e-2 | 0.85 | 1.682e-1 | 0.53 | 2.922e-2 | 1.28
$2^{-4}$ | 8.379e-2 | 0.19 | 9.731e-3 | 0.91 | 1.005e-1 | 0.74 | 8.397e-3 | 1.80
$2^{-5}$ | 4.399e-2 | 0.93 | 3.972e-3 | 1.29 | 4.647e-2 | 1.11 | 2.351e-3 | 1.84

Hence, we solve the convection-dominated problem in Example 7.3 and perform a
mesh refinement study for the M-EAVE and EAVE methods on the two Voronoi
meshes with varying mesh sizes $h$ and $\epsilon=10^{-2}$. Our test checks the
$A$-norm error, $\lVert e_{h}\rVert_{A}$ where $e_{h}=u_{I}-u_{h}$, equivalent
to the $H^{1}$-error. Table 1 shows that the M-EAVE and EAVE methods have
first-order convergence in the $A$-norm, and the $L^{\infty}$-errors tend to
converge to zero of almost second-order. We see from Table 1 that the M-EAVE
method performs better in the orders of convergence and magnitude of the
errors in most cases. Thus, we conclude that the M-EAVE method satisfying the
monotonicity property is beneficial for solving convection-dominated problems
on such Voronoi meshes.

###  7.2 Accuracy test on general polygonal meshes

This test confirms that the EAVE method (6.3) works with general polygonal
meshes by checking the accuracy. In order to confirm the accuracy of the EAVE
method, we conduct a mesh refinement study with different mesh sizes $h$ and
the fixed $\epsilon=10^{-2}$ on various polygonal meshes. We also solve the
problem in Example 7.3 on the following mesh types (see also Figure 3 as an
example):

  * •

voro: a Voronoi tessellation;

  * •

opti: a Voronoi tessellation where the cell shapes are optimized via a Lloyd
algorithm [41];

  * •

ncvx: a mesh composed by structured non-convex polygons and triangles.

voro | opti | ncvx
---|---|---
|  |
Figure 3: Various polygonal meshes with $h=2^{-3}$.

Table 2 shows the EAVE method’s $A$-norm error and $L^{\infty}$-error for the
voro, opti, and ncvx meshes.

Table 2: A mesh refinement study for EAVE on various meshes with $\epsilon=10^{-2}$. | voro | opti | ncvx
---|---|---|---
$h$ | $\lVert u_{I}-u_{h}\rVert_{A}$ | Order | $\lVert u_{I}-u_{h}\rVert_{A}$ | Order | $\lVert u_{I}-u_{h}\rVert_{A}$ | Order
$2^{-4}$ | 2.383e-1 | - | 1.078e-1 | - | 1.355e-1 | -
$2^{-5}$ | 1.586e-1 | 0.59 | 6.250e-2 | 0.79 | 2.488e-2 | 2.45
$2^{-6}$ | 8.057e-2 | 0.98 | 2.914e-2 | 1.10 | 4.987e-3 | 2.32
$2^{-7}$ | 4.263e-2 | 0.92 | 1.414e-2 | 1.04 | 9.253e-4 | 2.43
$2^{-8}$ | 2.126e-2 | 1.00 | 6.902e-3 | 1.03 | 1.665e-4 | 2.47
$h$ | $\lVert u_{I}-u_{h}\rVert_{\infty}$ | Order | $\lVert u_{I}-u_{h}\rVert_{\infty}$ | Order | $\lVert u_{I}-u_{h}\rVert_{\infty}$ | Order
$2^{-4}$ | 3.714e-2 | - | 1.139e-2 | - | 2.563e-2 | -
$2^{-5}$ | 1.494e-2 | 1.31 | 5.218e-3 | 1.13 | 3.551e-3 | 2.85
$2^{-6}$ | 4.183e-3 | 1.84 | 1.797e-3 | 1.54 | 5.052e-4 | 2.81
$2^{-7}$ | 1.224e-3 | 1.77 | 4.414e-4 | 2.03 | 6.625e-5 | 2.93
$2^{-8}$ | 3.123e-4 | 1.97 | 1.140e-4 | 1.95 | 8.424e-6 | 2.98

On all the given meshes, the orders of convergence for the $A$-norm errors are
at least first-order. Also, the $L^{\infty}$-errors converge to zero of almost
second-order. For the ncvx mesh, both $A$-norm error and $L^{\infty}$-error
show super-convergence because the ncvx mesh consists of structured grid
points and the points are aligned with the boundary layer. Therefore, these
numerical results confirm that the EAVE method produces accurate numerical
solutions on general polygonal meshes.

###  7.3 Test for convection-dominated case

In this test, we check the performance of the M-EAVE and EAVE methods compared
to other stabilized numerical methods for convection-dominated problems.
Various numerical methods having the same order of convergence are employed to
solve the problem in Example 7.3, and they are denoted as follows:

  * •

FE: Conforming finite element method with piecewise linear functions;

  * •

SUPG: Streamline-upwind Petrov-Galerkin method [17] with piecewise linear
functions and the stability parameter

$s_{E}=\frac{0.25h^{2}_{T}}{\epsilon
P_{e}}\left(1-\frac{1}{P_{e}}\right)\quad\text{and}\quad
P_{e}=\frac{\beta_{T}h_{T}}{2\epsilon},$

where $\beta_{T}=\|\bm{\beta}\|_{L^{\infty}(T)}$ and $h_{T}=|T|^{1/2}$;

  * •

EAFE: Edge-averaged finite element scheme (2.8);

  * •

M-EAVE: Monotone edge-averaged virtual element method (5.4);

  * •

EAVE: Edge-averaged virtual element method (6.3).

The performance of the stabilized numerical methods can be explained by the
presence of unexpected oscillations in their numerical solutions. According to
the exact solution (7.3), numerical solutions are expected to capture the
boundary layer without oscillation. However, even though numerical methods
imply small $H^{1}$\- or $L^{2}$-errors, nonphysical oscillations may locally
occur on their numerical solutions because they may not satisfy the discrete
maximum principle with the given mesh size. The $L^{\infty}$-error, $\lVert
u_{I}-u_{h}\rVert_{\infty}$, can identify such nonphysical oscillations, where
$u_{I}$ is the nodal value interpolant of $u$. For comparisons across the
stabilized numerical methods, we solved the problem in Example 7.3 with
varying $\epsilon$ values and mesh size $h$, respectively.

Figure 4: Error profiles of the numerical methods with varying $\epsilon$
(left) and $h$ (right).

Figure 4 shows the $L^{\infty}$-errors of the numerical methods in two ways;
(i) with varying $\epsilon$ values from $10^{-2}$ to $10^{-6}$ and the fixed
mesh size $h=2^{-3}$, (ii) with different mesh size from $2^{-3}$ to $2^{-7}$
and the fixed $\epsilon=10^{-2}$. As $\epsilon$ decreases, the
$L^{\infty}$-errors of the FE method rapidly increase, which shows the need
for stabilization techniques for convection-dominated cases. On the other
hand, the $L^{\infty}$-errors of the stabilized methods (SUPG, EAFE, M-EAVE,
and EAVE) seem independent of $\epsilon$. As $h$ decreases, all the methods
produce the $L^{\infty}$-errors decreasing in a similar order, but the errors
of edge-averaged methods (EAFE, M-EAVE, and EAVE) are almost a hundred times
smaller than the SUPG method.

In order to perform a quantitative comparison, we display the numerical
solutions obtained by the FE, SUPG, EAFE, and M-EAVE methods with
$\epsilon=10^{-2}$ and $h=2^{-4}$ in Figure 5. Since the solution of the EAVE
method is almost the same as the M-EAVE, we omit EAVE here.

FE | SUPG | EAFE |  M-EAVE
---|---|---|---
|  |  |

(a) Numerical solutions for the methods
FE SUPG EAFE M-EAVE
(b) Cross sections of the numerical solutions from $x=0$ to $x=0.9$

Figure 5: Comparison of the numerical solutions with $\epsilon=10^{-2}$ and
$h=2^{-4}$.

As expected, the numerical solution for the FE method contains spurious
oscillations around the boundary layer, and the magnitude of the biggest
oscillation is over 20% of the maximum of the exact solution. The SUPG, a
stabilized method, provides a more stable numerical solution than the FE
method. However, the solution still contains a nonphysical oscillation whose
magnitude is over 10% of the exact solution. These spurious oscillations may
deteriorate the quality of the numerical solutions and cause difficulty in
capturing hidden boundary layers. On the other hand, the EAFE and M-EAVE
methods satisfy the monotonicity property, so they produce numerical solutions
without such oscillations and nicely capture the sharp boundary layer.
Therefore, we confirm that the M-EAVE and EAVE methods produce a quality
approximation on polygonal meshes.

##  8 Conclusion

We have developed edge-averaged virtual element methods for convection-
diffusion problems by generalizing the edge-averaged finite element
stabilization [43] to the virtual element methods. A proper flux approximation
and mass lumping allow us to obtain a monotone EAVE method on a Voronoi mesh
with a dual Delaunay triangulation consisting of acute triangles. The
Bernoulli function and the geometric data of such a Voronoi mesh readily
compute the bilinear form in the monotone method, while the monotone method
guarantees stable numerical solutions. We also presented a general framework
for EAVE methods, proposing another EAVE bilinear form working with general
polygonal meshes. Our numerical experiments demonstrated the effect of the
monotonicity property and the optimal order of convergence on different mesh
types.

Extending the EAVE methods to three-dimensional cases will be one of our
future research directions. The EAVE methods can be also used for application
problems that include convection dominance and require polygonal meshes.

## References

  * Adler et al. [2023] James H Adler, Casey Cavanaugh, Xiaozhe Hu, Andy Huang, and Nathaniel Trask.  A stable mimetic finite-difference method for convection-dominated diffusion equations.  _SIAM Journal on Scientific Computing_ , 45(6):A2973–A3000, 2023.
  * Bank et al. [1998] Randolph E. Bank, W.M. Coughran, Jr., and Lawrence C. Cowsar.  The finite volume Scharfetter-Gummel method for steady convection diffusion equations.  _Computing and Visualization in Science_ , 1(3):123–136, 1998.
  * Bank et al. [2017] Randolph E Bank, Panayot S Vassilevski, and Ludmil T Zikatanov.  Arbitrary dimension convection–diffusion schemes for space–time discretizations.  _Journal of Computational and Applied Mathematics_ , 310:19–31, 2017.
  * Barrenechea et al. [2018] Gabriel R. Barrenechea, Volker John, Petr Knobloch, and Richard Rankin.  A unified analysis of algebraic flux correction schemes for convection-diffusion equations.  _SeMA Journal_ , 75(4):655–685, 2018.
  * Beirão da Veiga et al. [2013] Lourenco Beirão da Veiga, Franco Brezzi, Andrea Cangiani, Gianmarco Manzini, Luisa Donatella Marini, and Alessandro Russo.  Basic principles of virtual element methods.  _Mathematical Models and Methods in Applied Sciences_ , 23(01):199–214, 2013.
  * Beirão da Veiga et al. [2014] Lourenco Beirão da Veiga, Franco Brezzi, Luisa Donatella Marini, and Alessandro Russo.  The hitchhiker’s guide to the virtual element method.  _Mathematical models and methods in applied sciences_ , 24(08):1541–1573, 2014.
  * Beirão da Veiga et al. [2016] Lourenco Beirão da Veiga, Franco Brezzi, Luisa Donatella Marini, and Alessandro Russo.  ${H}(\text{div})$ and ${H}(\textbf{curl})$-conforming virtual element methods.  _Numerische Mathematik_ , 133(2):303–332, 2016\.
  * Beirão da Veiga et al. [2017a] Lourenco Beirão da Veiga, Franco Brezzi, Franco Dassi, Luisa Donatella Marini, and Alessandro Russo.  Virtual element approximation of 2D magnetostatic problems.  _Computer Methods in Applied Mechanics and Engineering_ , 327:173–195, 2017a.
  * Beirão da Veiga et al. [2017b] Lourenco Beirão da Veiga, Franco Brezzi, Luisa Donatella Marini, and Alessandro Russo.  Serendipity face and edge VEM spaces.  _Rendiconti Lincei_ , 28(1):143–180, 2017b.
  * Beirão da Veiga et al. [2018] Lourenco Beirão da Veiga, Franco Brezzi, Franco Dassi, Luisa Donatella Marini, and Alessandro Russo.  Lowest order virtual element approximation of magnetostatic problems.  _Computer Methods in Applied Mechanics and Engineering_ , 332:343–362, 2018.
  * Beirão da Veiga et al. [2021] Lourenco Beirão da Veiga, Franco Dassi, Carlo Lovadina, and Giuseppe Vacca.  SUPG-stabilized virtual elements for diffusion-convection problems: a robustness analysis.  2021\.
  * Benedetto et al. [2016] Matias Fernando Benedetto, Stefano Berrone, Andrea Borio, Sandra Pieraccini, and Stefano Scialó.  Order preserving SUPG stabilization for the virtual element formulation of advection–diffusion problems.  _Computer Methods in Applied Mechanics and Engineering_ , 311:18–40, 2016.
  * Berrone et al. [2018] Stefano Berrone, Andrea Borio, and Gianmarco Manzini.  SUPG stabilization for the nonconforming virtual element method for advection–diffusion–reaction equations.  _Computer Methods in Applied Mechanics and Engineering_ , 340:500–529, 2018.
  * Brezzi et al. [2006] F Brezzi, M Fortin, and LD Marini.  Error analysis of piecewise constant pressure approximations of darcy’s law.  _Computer methods in applied mechanics and engineering_ , 195(13-16):1547–1559, 2006.
  * Brezzi et al. [1989a] Franco Brezzi, Luisa Donatella Marini, and Paola Pietra.  Numerical simulation of semiconductor devices.  _Computer methods in applied mechanics and engineering_ , 75(1-3):493–514, 1989a.
  * Brezzi et al. [1989b] Franco Brezzi, Luisa Donatella Marini, and Paola Pietra.  Two-dimensional exponential fitting and applications to drift-diffusion models.  _SIAM Journal on Numerical Analysis_ , 26(6):1342–1355, 1989b.
  * Brooks and Hughes [1982] Alexander N. Brooks and Thomas J.R. Hughes.  Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations.  _Computer Methods in Applied Mechanics and Engineering_ , 32(1-3):199–259, 1982.
  * Burman and Ern [2005] Erik Burman and Alexandre Ern.  Stabilized galerkin approximation of convection-diffusion-reaction equations: discrete maximum principle and convergence.  _Mathematics of computation_ , 74(252):1637–1652, 2005.
  * Burman and Hansbo [2004] Erik Burman and Peter Hansbo.  Edge stabilization for Galerkin approximations of convection-diffusion-reaction problems.  _Computer Methods in Applied Mechanics and Engineering_ , 193(15-16):1437–1453, 2004.
  * Cao and Chen [2018] Shuhao Cao and Long Chen.  Anisotropic error estimates of the linear virtual element method on polygonal meshes.  _SIAM Journal on Numerical Analysis_ , 56(5):2913–2939, 2018.
  * Chen [2009] Long Chen.  iFEM: An Integrated Finite Element Methods Package in MATLAB.  _Technical Report, University of California at Irvine_ , 2009.
  * Chen and Huang [2018] Long Chen and Jianguo Huang.  Some error analysis on virtual element methods.  _Calcolo_ , 55(1):1–23, 2018.
  * Ciarlet and Raviart [1973] Philippe G Ciarlet and P-A Raviart.  Maximum principle and uniform convergence for the finite element method.  _Computer methods in applied mechanics and engineering_ , 2(1):17–31, 1973.
  * Dörfler [1999] Willy Dörfler.  Uniform error estimates for an exponentially fitted finite element method for singularly perturbed elliptic equations.  _SIAM journal on numerical analysis_ , 36(6):1709–1738, 1999.
  * Elman et al. [2014] Howard C. Elman, David J. Silvester, and Andrew J. Wathen.  _Finite Elements and Fast Iterative Solvers: with Applications in Incompressible Fluid Dynamics_.  Numerical Mathematics and Science, 2014.
  * Eymard et al. [2000] Robert Eymard, Thierry Gallouët, and Raphaèle Herbin.  Finite volume methods.  _Handbook of numerical analysis_ , 7:713–1018, 2000.
  * Gilbarg and Trudinger [1977] David Gilbarg and Neil S. Trudinger.  _Elliptic Partial Differential Equations of Second Order_ , volume 224.  Springer, 1977.
  * Godunov [1959] Sergei K. Godunov.  A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics.  _Matematicheskii Sbornik_ , 47(89):271–306, 1959\.
  * Haugazeau and Lacoste [1993] Yves Haugazeau and Patrick Lacoste.  Condensation de la matrice de masse pour les éléments finis mixtes de h (rot).  _Comptes rendus de l’Académie des sciences. Série I, Mathématique_ , 316(5):509–512, 1993.
  * John and Knobloch [2007] Volker John and Petr Knobloch.  On spurious oscillations at layers diminishing (SOLD) methods for convection–diffusion equations: Part I–A review.  _Computer methods in applied mechanics and engineering_ , 196(17-20):2197–2215, 2007.
  * Kikuchi [1977] Fumio Kikuchi.  Discrete maximum principle and artificial viscosity in finite element approximations to convective diffusion equations.  _Tokyo, University, Institute of Space and Aeronautical Science, Report no. 550_ , 42:153–166, 1977.
  * Knobloch [2010] Petr Knobloch.  A generalization of the local projection stabilization for convection-diffusion-reaction equations.  _SIAM Journal on Numerical Analysis_ , 48(2):659–680, 2010.
  * Lazarov and Zikatanov [2012] Raytcho D. Lazarov and Ludmil T. Zikatanov.  An exponential fitting scheme for general convection-diffusion equations on tetrahedral meshes.  _arXiv:1211.0869_ , 2012.
  * Lee [2021] Seulip Lee.  _Edge-averaged virtual element methods for convection-diffusion problems_.  University of California, Irvine, 2021.
  * Li and Feng [2021] Yang Li and Minfu Feng.  A local projection stabilization virtual element method for convection-diffusion-reaction equation.  _Applied Mathematics and Computation_ , 411:126536, 2021\.
  * Li and Zikatanov [2022] Yuwen Li and Ludmil T Zikatanov.  New stabilized ${P}1\times{P}0$ finite element methods for nearly inviscid and incompressible flows.  _Computer Methods in Applied Mechanics and Engineering_ , 393:114815, 2022.
  * Markowich and Zlámal [1988] Peter A. Markowich and Miloš A. Zlámal.  Inverse-average-type finite element discretizations of selfadjoint second-order elliptic problems.  _Mathematics of computation_ , 51(184):431–449, 1988.
  * Roos [2008] Hans-Görg Roos.  _Robust numerical methods for singularly perturbed differential equations_.  Springer, 2008.
  * Scharfetter and Gummel [1969] Donald L. Scharfetter and Hermann K. Gummel.  Large-signal analysis of a silicon read diode oscillator.  _IEEE Transactions on electron devices_ , 16(1):64–77, 1969.
  * Stynes [2005] Martin Stynes.  Steady-state convection-diffusion problems.  _Acta Numerica_ , 14:445–508, 2005.
  * Talischi et al. [2012] Cameron Talischi, Glaucio H. Paulino, Anderson Pereira, and Ivan F.M. Menezes.  PolyMesher: A general-purpose mesh generator for polygonal elements written in MATLAB.  _Structural and Multidisciplinary Optimization_ , 45(3):309–328, 2012.
  * Wu and Xu [2020] Shuonan Wu and Jinchao Xu.  Simplex-averaged finite element methods for ${H}(\text{grad})$, ${H}(\text{curl})$, and ${H}(\text{div})$ convection-diffusion problems.  _SIAM Journal on Numerical Analysis_ , 58(1):884–906, 2020.
  * Xu and Zikatanov [1999] Jinchao Xu and Ludmil T. Zikatanov.  A monotone finite element scheme for convection-diffusion equations.  _Mathematics of Computation_ , 68(228):1429–1446, 1999.