[2]Umberto Zerbinati

1]Dipartimento di Matematica, Uiniversità di Pavia, Via Ferrata 5, Pavia,
27100, Italy

2]Mathematical Institute, University of Oxford, Andrew Wiles Building, Oxford,
OX2 6GG, United Kingdom

# When rational functions meet virtual elements:
The lightning Virtual Element Method

Manuel Luigi Trezzi<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [

###### Abstract

We propose a lightning Virtual Element Method that eliminates the
stabilisation term by actually computing the virtual component of the local
VEM basis functions using a lightning approximation. In particular, the
lightning VEM approximates the virtual part of the basis functions using
rational functions with poles clustered exponentially close to the corners of
each element of the polygonal tessellation. This results in two great
advantages. First, the mathematical analysis of a priori error estimates is
much easier and essentially identical to the one for any other non-conforming
Galerkin discretisation. Second, the fact that the lightning VEM truly
computes the basis functions allows the user to access the point-wise value of
the numerical solution without needing any reconstruction techniques. The cost
of the local construction of the VEM basis is the implementation price that
one has to pay for the advantages of the lightning VEM method, but the
embarrassingly parallelizable nature of this operation will ultimately result
in a cost-efficient scheme almost comparable to standard VEM and FEM.

###### keywords:

virtual element method, lightning Laplace, rational functions, partial
differential equations

##  1 Introduction

Since its introduction in [1], the Virtual Element Method (VEM) has been
recognised as a valuable tool for the solution of partial differential
equations (PDEs). One of the key advantages of the VEM is that it allows any
type of polygonal meshes, thanks to the introduction of a virtual component in
the local basis functions. We will discuss this aspect in a later section.
Curved meshes can also be used more naturally than iso-parametric finite
elements [2]. Furthermore, given the rising interest in structure-preserving
discretisation, it is worth mentioning that the VEM allows the mimicry of many
interesting physical structures that arise at the level of the continuous PDE
[3]. For example, the VEM, unlike the standard finite element method (FEM),
can be used to construct arbitrary low-order divergence-free $H^{1}$
conforming discretisations for Stokes flow [3], without any restriction on the
mesh. The sophisticated mathematical infrastructure that allows the proof of
accurate a priori error estimates, which are a key advantage of the FEM over
other numerical schemes, can also be applied with some adjustments to the VEM.
Another feature of the VEM, perhaps the one of greatest importance, is the
ability to easily produce discretisations with a high order of continuity, a
property that is well known to be a weakness of the classical FEM. It is also
worth mentioning that VEM accuracy can be improved with the $p$ and $hp$
versions of VEM introduced in [4], clearly reflecting the intimate relation
between the FEM and VEM.

The advantages just mentioned lead to the application of the VEM to a wide
variety of problems from linear elasticity [5, 6, 7], fluid-dynamics [8, 9,
10], fourth-order problems [11, 12, 13], acoustic wave propagation and
Helmholtz problem [14, 15], and magnetostatics [16, 17]. While the
mathematical infrastructure is what makes the VEM shine and reveals the
advantages listed above, the challenges of the VEM are the practical aspects
of its implementation. The true novelty behind the VEM is the use of the so-
called projection operators that allow assembling stiffness and mass matrices
without the necessity of computing the virtual component of the basis
functions. Resorting to projection operators results in a lack of coercivity
that requires the introduction of an additional stabilisation term in the weak
formulation of the problem.

The lightning VEM proposed here eliminates the stabilisation term by actually
computing, in an extremely efficient manner, the virtual component of the
local VEM basis functions. In particular, the lightning VEM approximates the
virtual part of the basis functions using rational functions with poles
clustered exponentially close to the corners of each element of the polygonal
tessellation. It is worth mentioning that in the literature other ideas have
been proposed in order to get rid of the stabilisation term [18, 19, 20, 21].
The key difference between the lightning VEM and these methods is the fact
that the lightning VEM does not make use of any projection operators.
Therefore, the mathematical analysis of a priori error estimates is much
easier and essentially identical to that for any other non-conforming Galerkin
discretisation.

Furthermore, the fact that the lightning VEM truly computes the basis function
allows the user to access pointwise values of the numerical solution without
the need for any reconstruction technique. This is particularly appealing
considering that a common reconstruction technique is based on polynomial
interpolation and may require the triangulation of the polygonal mesh.

The local construction of the VEM basis is the implementation price that one
has to pay for the advantages of the lightning VEM method, but as we will see
the embarrassingly parallelizable nature of this operation will ultimately
result in a cost-efficient scheme compared to standard VEM and FEM.

Before diving in the core aspects of the paper, the authors would like to
stress the difference between the meaning of the word “virtual” in the context
of the VEM and in the context of the lightning VEM. In the first case“virtual”
refers to the fact that we have no point-wise knowledge of the basis functions
in the interior of an element. In the lightning VEM we compute a set of basis
functions, for which we can access point-wise values, and “virtual” refers to
the fact that we are approximating the standard virtual element space.

##  2 Virtual Element Method

For the sake of simplicity, we focus our attention on the Laplace problem. The
VEM has been applied to a large number of equations but we want to keep the
focus on the simplest scenario. Given a polygonal open, bounded domain
$\Omega\subset\mathbb{R}^{2}$ with boundary $\Gamma$, we consider the problem
of finding a function $u:\Omega\to\mathbb{R}$ such that

$\left\\{\begin{aligned} -\Delta u&=f\quad\text{in }\Omega\,,\\\
u&=0\quad\text{on }\Gamma\,,\end{aligned}\right.$ (1)

where $f\in L^{2}(\Omega)$ is the load term. It is well known that the
standard weak formulation of (1) reads as

$\left\\{\begin{aligned} &\text{find $u\in V\coloneqq H^{1}_{0}(\Omega)$ such
that: }\\\ &a(u,v)=\langle f,v\rangle\quad\forall v\in
H_{0}^{1}(\Omega)\,,\end{aligned}\right.$ (2)

where the bilinear form $a(\cdot,\cdot):V\times V\to\mathbb{R}$ is defined as

$\forall u,v\in V\,,\quad a(u,v)\coloneqq\int_{\Omega}\nabla u\cdot\nabla
v\,.$ (3)

###  2.1 The local space

We decompose the domain $\Omega$ in a tessellation $\mathcal{T}_{h}$ by a
finite number of non-overlapping convex polygons $K$. In particular, we assume
that there exists a positive constant $\rho$ such that every element
$K\in\mathcal{T}_{h}$ is star-shaped with respect to a ball of a radius
greater or equal than $\rho h_{K}$, where $h_{K}$ denotes the diameter of the
element $K$. Let $k\geq 1$ be the “order” of the method. For every element
$K\in\mathcal{T}_{h}$ we define the local virtual element space as

$V_{h}^{k}(K)\coloneqq\left\\{v_{h}\in H^{1}(K)\,\text{ such that }\,\Delta
v_{h}\in\mathbb{P}_{k-2}(K)\,\text{ and }\,v_{h}\in\mathbb{B}_{K}(\partial
K)\right\\},$ (4)

where

$\mathbb{B}_{k}(\partial K)\coloneqq\left\\{v_{h}\in\mathcal{C}^{0}(\partial
K)\,\text{ such that }\,\forall e\in\partial
K,\,v_{h}|_{e}\in\mathbb{P}_{k}(e)\right\\},$ (5)

where $e$ denotes an edge of $\partial K$ and we have denoted
$\mathbb{P}_{-1}\coloneqq\\{0\\}$. One can check that the following quantities
represent a set of degrees of freedom for the space $V_{h}^{k}(K)$:

  1. 1.

$\mathcal{V}^{K,k}$: the pointwise values of $v_{h}$ at the vertexes of the
polygon $K$,

  2. 2.

$\mathcal{E}^{K,k}$: the values of $v_{h}$ at $k-1$ internal points of a
Gauss-Lobatto quadrature for every edge $e\in\partial K$,

  3. 3.

$\mathcal{P}^{K,k}$: the moments
$\dfrac{1}{|K|}\int_{K}v_{h}\,m_{\alpha\beta}\,{\rm d}K\,$, $\forall
m_{\alpha\beta}\in\mathcal{M}_{k-2}(K)$ where $\mathcal{M}_{k-2}(K)$ is the
set of monomials defined as

$\mathcal{M}_{k-2}(K)\coloneqq\left\\{m_{\alpha\beta}\coloneqq\left(\dfrac{x-x_{K}}{h_{K}}\right)^{\alpha}\left(\dfrac{y-y_{K}}{h_{K}}\right)^{\beta}\
\alpha,\beta\in\mathbb{N}\,,\alpha+\beta\leq k-2\right\\}.$ (6)

Details of the proof are in [1]. A graphic representation of the degrees of
freedom is given in Figure 1.

$k=1$

$k=2$

$k=3$

Figure 1: Degrees of freedom on a pentagon.

Thanks to the definition of the degrees of freedom, it is possible to
construct the projection operator
$\Pi^{\nabla,K}_{k}:V_{h}^{k}(K)\to\mathbb{P}_{k}(K)$ defined as

$\left\\{\begin{aligned}
&\int_{K}\nabla\,p_{k}\cdot\nabla(v_{h}-\,{\Pi}_{k}^{\nabla,K}v_{h})\,{\rm
d}K=0\quad\text{for all $v_{h}\in V_{h}^{k}(K)$ and
$p_{k}\in\mathbb{P}_{k}(K)$,}\\\ &\int_{\partial
K}(v_{h}-{\Pi}_{k}^{\nabla,K}v_{h})\,{\rm d}s=0\,.\end{aligned}\right.$

Indeed, thanks to integration by parts, we note that

$\int_{K}\nabla\,p_{k}\cdot\nabla\,v_{h}=-\int_{K}v_{h}\,\Delta p_{k}\,{\rm
d}K+\int_{\partial K}\dfrac{p_{k}}{\partial n}\,v_{h}\,{\rm d}s.$ (7)

The first integral is known thanks to $\mathcal{P}^{K,k}$, for the second one
we use $\mathcal{E}^{K,k}$ and $\mathcal{V}^{K,k}$. The operator
$\Pi^{\nabla,K}$ is an orthogonal projection into the space of polynomials
with respect to the $H^{1}$ seminorm. This is the cornerstone around we can
construct a discretisation of the bilinear form $a(\cdot,\cdot)$. The global
space is obtained by gluing together the local spaces

$V_{h}^{k}(\mathcal{T}_{h})=\\{v_{h}\in V\quad\text{s.t.}\quad v_{h}|_{K}\in
V_{h}(K)\quad\text{for all $K\in\mathcal{T}_{h}$}\\}\,,$

therefore obtaining a space characterized by the following degrees of freedom:

  1. 1.

$\mathcal{V}^{k}$: the values of $v_{h}$ at the vertices;

  2. 2.

$\mathcal{E}^{k}$: the values of $v_{h}$ at $k-1$ points on each edge $e$;

  3. 3.

$\mathcal{P}^{k}$: the moments up to order $k-2$ for each element
$E\in\mathcal{T}_{h}$.

###  2.2 The discrete problem

First, we decompose the global bilinear form into local contributions

$a(u,v)=\sum_{K\in\mathcal{T}_{h}}a^{K}(u,v)\,.$ (8)

We point out that, except for a particular structures of the element $K$, we
don’t have an analytic expression for all the functions in $V_{h}^{k}(K)$.
Hence, given two generic virtual functions, we are not able to compute the
quantity

$a^{K}(u_{h},v_{h})=\int_{K}\nabla u_{h}\cdot\nabla v_{h}\,{\rm d}K\,,\qquad
u_{h},v_{h}\in V_{h}^{k}(\mathcal{T}_{h})$ (9)

We would like to construct a discrete bilinear form
$a_{h}^{K}:V_{h}^{k}(K)\times V_{h}^{k}(K)\to\mathbb{R}$ that is computable
for all the virtual functions and acts as a discrete counterpart of
$a^{K}(\cdot,\cdot)$. The idea is to split the virtual functions as

$v_{h}=\Pi^{\nabla,K}_{k}v_{h}+(I-\Pi^{\nabla,K}_{k})v_{h}\,.$ (10)

Thanks to this choice, the bilinear form is split as

$\displaystyle a^{K}(u_{h},v_{h})$
$\displaystyle=a^{K}(\Pi^{\nabla,K}_{k}u_{h},\Pi^{\nabla,K}_{k}v_{h})+a^{K}((I-\Pi^{\nabla,K}_{k})u_{h},(I-\Pi^{\nabla,K}_{k})v_{h})$
(11)
$\displaystyle+a^{K}(\Pi^{\nabla,K}_{k}u_{h},(I-\Pi^{\nabla,K}_{k})v_{h})+a^{K}((I-\Pi^{\nabla,K}_{k})u_{h},\Pi^{\nabla,K}_{k}v_{h})\,.$

Thanks to orthogonality of $\Pi^{\nabla,K}_{k}$, the last two terms are equal
to zero. The first term involves only polynomials hence is computable. This
term is known in the VEM literature as the consistency term. The only thing
that remains to be done is to handle the second term. The idea in the VEM
framework is to replace it with a computable bilinear form
$\mathcal{S}^{K}(\cdot,\cdot):V_{h}^{k}(K)\times V_{h}^{k}(K)\to\mathbb{R}$
that satisfies

$\alpha_{*}|v_{h}|_{1,E}^{2}\leq\mathcal{S}(v_{h},v_{h})\leq\alpha^{*}|v_{h}|_{1,E}^{2}\,,\qquad\text{for
all $v_{h}\in{\rm Ker}(\Pi^{\nabla,K}_{k})$}\,,$ (12)

for two positive uniform constants $\alpha_{*}$ and $\alpha^{*}$. This term is
called the stability term. We define the bilinear form

$a_{h}^{K}(u_{h},v_{h})\coloneqq
a^{K}(\Pi^{\nabla,K}_{k}u_{h},\Pi^{\nabla,K}_{k}v_{h})+\mathcal{S}^{K}((I-\Pi^{\nabla,K}_{k})u_{h},(I-\Pi^{\nabla,K}_{k})v_{h})\,.$
(13)

The global bilinear form is obtained summing all the local contributions

$a_{h}(u_{h},v_{h})\coloneqq\sum_{K\in\mathcal{T}_{h}}a^{K}_{h}(u_{h},v_{h})$
(14)

It remains to discretise the load term. A standard choice is the following
procedure:

  1. 1.

if $k\geq 2$, we replace $f$ with $f_{h}$ defined locally as the
$L^{2}$-orthogonal projection into the space $\mathbb{P}_{k-2}(K)$. In detail,
we consider

$\langle f_{h},v_{h}\rangle=\sum_{K\in\mathcal{T}_{h}}\int_{K}f_{h}\,v_{h}{\rm
d}K=\sum_{K\in\mathcal{T}_{h}}\int_{K}f\,\Pi^{0,E}_{k-2}v_{h}{\rm d}K\,;$ (15)

  2. 2.

if $k=1$, we replace $f$ with a piecewise constant and define

$\langle
f_{h},v_{h}\rangle=\sum_{K\in\mathcal{T}_{h}}\int_{K}P_{0}^{K}f\,\overline{v}_{h}{\rm
d}K\,$ (16)

where

$\overline{v}_{h}=\dfrac{1}{n}\sum_{i=1}^{n}v_{h}(V_{i})\,,\quad
V_{i}=\text{vertices of $K$}\,.$ (17)

The discrete problem reads as

$\left\\{\begin{aligned} &\text{find $u_{h}\in V_{h}^{k}(\Omega_{h})$ such
that: }\\\ &a_{h}(u_{h},v_{h})=\langle f_{h},v_{h}\rangle\quad\forall v_{h}\in
V_{h}^{k}(\Omega_{h})\,.\end{aligned}\right.$ (18)

######  Remark 1.

If we desire an $H^{1}$ conforming discretisation the most common choice of
virtual operator is the negative Laplacian as in (4). The well-posedness of
the local virtual element space is guaranteed even if we choose any other
elliptic operator with sufficiently smooth coefficients. Yet to recover the,
optimal approximation property, it will be convenient also to assume that
$\mathbb{P}_{k}$ is a subset of the co-domain of chosen the operator when
applied to the space $\mathbb{P}_{k-2}$ with the constraint of having boundary
data in $\mathbb{P}_{k}$. A common misunderstanding when first approaching the
VEM is to think in terms of Trefftz methods, i.e. to assume that the operator
appearing in (4) has to be the same as the one appearing in (1). We direct the
reader interested in Trefftz methods to [22, 23]. We should think of the
operator appearing in (4) only in terms of the conformity we desire.

######  Remark 2.

In general, if we require an $H^{k}$ conforming discretisation the most common
choice of operator in (4) would be the polylaplacian $-\Delta^{k}$. Of course,
we need to modify slightly our definition of local VEM space to accommodate
the fact that we have to prescribe not only the value of the solution at the
boundary but also the first $k-1$ normal derivatives. More detail on this
topic can be found in [11, 12, 13, 24, 25].

######  Remark 3.

An important observation is that the VEM only provides the value of the
degrees of freedom as result. If one is interested in the value of the
discrete solution outside of $\mathcal{V}^{k}$ and $\mathcal{E}^{k}$ a
reconstruction technique has to be used. A common reconstruction technique is
to triangulate each polygonal element and perform a linear interpolation using
the nodal degrees of freedom. We will see later on that the lightning virtual
element method avoids this issue.

##  3 The lightning approximation

The key idea behind the lightning VEM is to consider a different local
discrete space which is obtained by approximating $V_{h}^{k}(K)$ using
rational functions with poles clustered exponentially close to the corners of
$K$. For the sake of clarity, we will focus our discussion on the lowest order
case, $k=1$. We will discuss in a later section how to deal with the case
$k>1$ and the case where higher conformity is required. Thanks to the
linearity of the Laplacian, if $K$ has $N_{k}$ edges, to construct
$V_{h}^{k}(K)$ we are interested in solving $N_{K}$ problems of the form

$\displaystyle\Delta\phi_{i}=0\;\;$ $\displaystyle\text{in}\;K,$ (19)
$\displaystyle\phi_{i}=\varphi_{i}\;\;$ $\displaystyle\text{on}\;\partial K.$

where $\varphi_{i}$ are the functions in $\mathbb{B}_{1}(\partial K)$ that are
equal to one on the $i$-th vertex and zero on the others. If we are working
with a structured rectangular mesh we can solve this problem exactly. As the
geometry of $K$ gets more complicated we need to find an efficient way of
solving (19). To address this problem we turn to the lightning Laplace solver
presented in [26, 27, 28]. The idea behind the lightning Laplace method is to
search for an approximation $\hat{\phi}_{i}$ of $\phi_{i}$ of the form

$\hat{\phi}_{i}=\Re\bigg{\\{}\sum_{j=0}^{N_{P}}\frac{a_{j}}{z-z_{j}}+\sum_{j=0}^{N_{Z}}b_{j}(z-z_{*})^{j}\bigg{\\}},$
(20)

where $\\{z_{j}\\}_{j=1}^{N_{P}}$ and $z_{*}$ are points in the complex plane
and $\Re$ denotes the real part of a complex number. Finding the optimal
coefficients, $\\{z_{j}\\}_{j=1}^{N_{P}}$ and $z_{*}$ to minimize
$\lVert\hat{\phi}_{i}-\phi_{i}\rVert_{L^{\infty}(K)}$ in general is a highly
non-linear problem. To transform this non-linear problem into a linear one the
position of the points $\\{z_{j}\\}_{j=1}^{N_{P}}$ and $z_{*}$ are fixed only
based upon the geometry of $K$. In particular, the $N_{P}$ poles are clustered
exponentially closer to the corners of the polygon $K$. Under this hypothesis,
the following result was proven in [29]:

######  Lemma 1.

Let $K$ be a convex polygon with corners $w_{1},\dots,w_{N_{K}}$ and let $f$
be an analytic function on $K$ that is analytic on the interior of each side
segment. Furthermore, we assume that $f$ can be analytically continued to a
disk near any corner $w_{k}$ with a slit along the exterior bisector of the
corner. Lastly, we assume that at each corner $w_{k}$ the function $f$
satisfies $f(z)-f(w_{k})=\mathcal{O}(\lVert z-w_{k}\rVert^{\delta})$ as $z\to
w_{k}$, for some $\delta>0$. Under these assumptions, there exists a rational
function $r$ of the form

$r=\Re\bigg{\\{}\sum_{j=0}^{N_{P}}\frac{a_{j}}{z-z_{j}}+\sum_{j=0}^{N_{Z}}b_{j}(z-z_{*})^{j}\bigg{\\}},$
(21)

with $N_{P}$ poles $z_{j}$ only at points exponentially clustered along the
exterior bisectors at the corners, such that the following approximation bound
holds, for $C>0$:

$\lVert f-r\rVert_{L^{\infty}(K)}\leq C_{K}e^{-C\sqrt{N_{P}}}.$ (22)

The previous lemma is a strong indication that the lightning Laplace scheme
will be able to converge to an extremely accurate solution, yet it can not be
applied in any straightforward manner to derive a priori error estimates with
respect to the $H^{1}$ and $L^{2}$ norms. An idea might be to mimic the
reasoning presented in [30], which can be used to derive an a priori error
estimate on a least squares collocation method in terms of Lebesgue constant
and the best approximation estimate presented above. The reason why this is
not a viable path is that we are interested here in $H^{1}$ error estimates,
which can not be produced using the type of $L^{\infty}$ bound presented in
[30]. To produce a priori error estimates for the overall scheme, we decide to
proceed adaptively when it comes to the resolution of (19) using the lightning
Laplace method. The exact algorithm used to compute the solution
$\hat{\phi}_{i}$ to (19) is presented in Algorithm 1 and we once again direct
the reader interested in more detail to [26].

Algorithm 1 For computing the solution $\hat{\phi}_{i}^{(e)}$ of (19).

1:$n_{1},\;n_{2},\;\dots$ such that $\sqrt{n_{1}},\;\sqrt{n_{2}},\;\dots$ is
uniformly distributed, a tolerance $\varepsilon$.

2:for $n=n_{1},n_{2},\dots$ do

3: Fix $N_{P}\approx nN_{K}$ and cluster the poles outside $K$ and
exponentially close the corners;

4: Fix $z_{*}$ in the interior of $K$ and fix $N_{Z}\approx N_{K}$ elements of
the monomial basis;

5: Choose $M=6N_{P}+6N_{Z}+1$ sample points on the boundary and clustered near
the corners;

6: Evaluate at the sample points $\varphi_{i}$ to obtain a matrix
$A\in\mathbb{R}^{M\times(2N_{P}+2N_{Z}+1)}$ and $\mathbf{d}\in\mathbb{R}^{M}$;

7: Solve the least squares problem $A\mathbf{x}\approx\mathbf{d}$ for the
coefficients $\mathbf{x}=(\mathbf{a},\mathbf{b})$;

8: Construct the function
$\hat{\phi}_{i}=\Re\bigg{\\{}\sum_{j=0}^{N_{P}}\frac{a_{j}}{z-z_{j}}+\sum_{j=0}^{N_{Z}}b_{j}(z-z_{*})^{j}\bigg{\\}}$
using the coefficients $\mathbf{a}$ and $\mathbf{b}$;

9: if $\lVert\hat{\phi}_{i}-\varphi_{i}\rVert_{H^{1/2}({\partial
K})}\leq\varepsilon$ then

10: break;

11: end if

12:end forreturn $\hat{\phi}_{i}$.

We then observe that since
$\sum_{j=0}^{N_{P}}\frac{a_{j}}{z-z_{j}}+\sum_{j=0}^{N_{Z}}b_{j}(z-z_{*})^{j}$
is analytic in the convex hull described by the vertex of any polygonal
element $K$ of the tessellation $\mathcal{T}_{h}$ where $z_{j}$ are the
vertices of $K$, then its real part is a harmonic function. We direct the
reader unfamiliar with this notion to [31, 32]. Since $\hat{\phi}_{i}$ is
harmonic we know that the function defined as $\hat{\phi}_{i}-\phi_{i}$ is
harmonic and by applying boundary regularity estimates for harmonic functions
we obtain the following result.

######  Proposition 2.

Let $\phi_{i}$ be as in (19) and $\hat{\phi}_{i}$ be the outcome of Algorithm
1. Then

$\lVert\hat{\phi}_{i}-\phi_{i}\rVert_{H^{1}(K)}\leq
C_{1}\varepsilon,\qquad\lVert\hat{\phi}_{i}-\phi_{i}\rVert_{L^{\infty}(\partial
K)}\leq C_{2}\varepsilon,$ (23)

where the constants $C_{1}$ and $C_{2}$ only depend on $\Omega$ and
$\partial\Omega$.

###### Proof.

Since we assumed the elements $K$ of our tessellation $\mathcal{T}_{h}$ are
polygonal, proof for the boundary estimates used to produce the first
inequality in (23) can be found in [33, 34]. The second estimate comes from
the Morrey’s inequality on a parametrization of the boundary [35]. ∎

From now on we will adopt the notation $\hat{V}_{h,\varepsilon}(K)$ to express
the discrete function space that is obtained from the solution of (19) using
the lightning Lapalce method on the element $K$ of the tessellation and
$\hat{V}_{h,\varepsilon}(\mathcal{T}_{h})$ to denote the global space
constructed starting from the various $\hat{V}_{h,\varepsilon}(K)$.

######  Remark 4.

It is clear that a fundamental step in the lightning VEM is solving the least
squares problem

$A\mathbf{x}\approx\mathbf{d}$ (24)

that appears in Algorithm 1. In particular the least squares problem (24) is
observed to be terribly ill-conditioned. Yet is still possible to solve (24)
with great accuracy if standard regularisation are adopted [26]. In particular
in [36] it is discussed the effect that the pole at infinity in (20) have on
the $2$-norm of $\mathbf{x}$. In fact bounding the $2$-norm of $\mathbf{x}$
will guarantee the solvability of (24) to a high degree of accuracy in spite
its ill-conditioned nature. An other route to solve (24) to a great deal of
accuracy, inspite its ill-conditioned nature, is to use the Vandermonde with
Arnoldi algorithm, introduced in [37] as discussed in [30, 29].

##  4 The lightning virtual element method

We begin by observing that the degrees of freedom $\mathcal{V}_{h}$ and
$\mathcal{E}_{h}$ are not decoupled in adjacent elements. In fact if
$\hat{\phi}_{i}$ had been a polynomial on the edges of the element $K$ as in
the standard VEM space $V_{h}(K)$, then this would have ensured the continuity
of the function $\hat{\phi}_{i}$ across the elements interface. Instead, since
each $\hat{\phi}_{i}$ is a rational function on the edges of the element $K$
even if over adjacent elements $\hat{\phi}_{i}$ would have the same value on
the edges and vertex degrees of freedom we still might have a jump across the
elements interface, as depicted Figure 2.

Figure 2: We know that the basis function $\hat{\phi}_{i,K_{1}}$ and
$\hat{\phi}_{i,K_{2}}$ corresponding to the $i$-th vertex and constructed
respectively on $K_{1}$ and $K_{2}$, match at the degrees of freedom here
denoted in red. Yet we have no guarantee that $\hat{\phi}_{1}$ and
$\hat{\phi}_{2}$ are continuous along the blue edge.$K_{1}$$K_{2}$

It is important to notice that this jump is bounded by adaptively solving the
lightning approximation problem in Algorithm 1, i.e.
$\lVert\hat{\phi}_{i}-\phi_{i}\rVert_{H^{1/2}({\partial K})}\leq\varepsilon$.
Since the function in $\hat{V}_{h,\varepsilon}(\mathcal{T}_{h})$ are no longer
continuous we need to consider a broken version of the bilinear form $a$ in
(3), i.e.

$\hat{a}(\hat{u}_{h},\hat{v}_{h})=\sum_{K\in\mathcal{T}_{h}}a^{K}(\hat{u}_{h},\hat{v}_{h})$
(25)

This bilinear form on $H^{1}(\Omega)$ is still symmetric positive definite
with respect to the broken $H^{1}(\Omega)$ norm and for the continuous
solution to the problem we have that $\hat{a}(u,v)=a(u,v)$ since $u,v\in
H^{1}(\Omega)$. Furthermore, in $\hat{V}_{h,\varepsilon}(\mathcal{T}_{h})$ the
kernel of the bilinear form $\hat{a}(\cdot,\cdot)$ are only functions that are
constant on each element of the tessellation, using the fact that such
constants must have the same value on the vertices of neighboring elements we
have that the boundary conditions impose that $\hat{a}(\cdot,\cdot)$ has a
trivial kernel. We now derive an a priori error estimate using the following
Lemma from [38].

######  Lemma 3.

Let $\hat{a}(\cdot,\cdot)$ be a symmetric positive definite bilinear form on
$H^{1}(\Omega)+\hat{V}_{h,\varepsilon}(\mathcal{T}_{h})$ which reduces to
$a(\cdot,\cdot)$ on $H^{1}(\Omega)$. We will further assume that $u$ is the
solution of (2) and $\hat{u}_{h}$ is such that:

$\hat{a}(\hat{u}_{h},\hat{v}_{h})=\langle f,\hat{v}_{h}\rangle\textrm{ for all
}\hat{v}_{h}\in\hat{V}_{h,\varepsilon}(\mathcal{T}_{h}),$ (26)

where in both the previous variation equation and (2) we assume $f\in
L^{2}(\Omega)\cap\hat{V}_{h,\varepsilon}(\mathcal{T}_{h})^{*}$. Under this
hypothesis the following best approximation estimate holds,

$\lVert
u-\hat{u}_{h}\rVert_{h}\leq\underset{\hat{v}_{h}\in\hat{V}_{h,\varepsilon}}{\inf}\lVert
u-\hat{v}_{h}\rVert_{h}+\underset{\hat{w}_{h}\in\hat{V}_{h,\varepsilon}}{\sup}\frac{\lvert\hat{a}(u-\hat{u}_{h},\hat{w}_{h})\rvert}{\lVert\hat{w_{h}}\rVert_{h}},$
(27)

where $\lVert\cdot\rVert_{h}:=\sqrt{\hat{a}(\cdot,\cdot)}$.

To apply the previous lemma, in order to obtain an a priori error estimate, we
need first to provide an estimate for the second term in the right-hand side
of (27).

$\displaystyle\hat{a}(u-\hat{u}_{h},\hat{w}_{h})$
$\displaystyle=\sum_{K\in\mathcal{T}_{h}}(\nabla
u,\nabla\hat{w}_{h})_{L^{2}(\Omega)}-(f,\hat{w}_{h})_{L^{2}(\Omega)}$ (28)
$\displaystyle=\sum_{K\in\mathcal{T}_{h}}\int_{\partial K}(\nabla
u\cdot\mathbf{n})\hat{w}_{h}=\sum_{e\in\mathscr{E}_{h}}\int_{e}\nabla
u\cdot\llbracket\hat{w}_{h}\rrbracket\,,$ (29)

where $\mathscr{E}_{h}$ are the internal edges of the tessellation
$\mathcal{T}_{h}$ and $\llbracket\hat{w}_{h}\rrbracket$ denotes the vector
jump of $\hat{w}_{h}$ across the edge $e$, i.e.

$\llbracket\hat{w}_{h}\rrbracket=\hat{w}_{h}\big{\lvert}_{K_{1}}\mathbf{n}_{1}+\hat{w}_{h}\big{\lvert}_{K_{2}}\mathbf{n}_{2}\,.$
(30)

We now rewrite the last term using the basis functions
$\hat{\phi}_{i}\big{\lvert}_{e}$,

$\sum_{e\in\mathscr{E}_{h}}\int_{e}\nabla
u\cdot\llbracket\hat{w}_{h}\rrbracket\leq\sum_{e\in\mathscr{E}_{h}}\sum_{i=1}^{N_{K}}\int_{e}\hat{w}_{i}\big{\lvert}_{e}\nabla
u\cdot\llbracket\hat{\phi}_{i}\rrbracket\leq\sum_{K\in\mathcal{T}_{h}}\varepsilon\lVert\hat{w}_{h}\rVert_{h,K}\lvert
u\rvert_{H^{2}(K)},$ (31)

where the last inequality comes from the trace theorem combined with Cauchy-
Schwarz inequality and the fact that we have constructed the basis function
$\hat{\phi}^{(e)}_{i}$ in such a way that
$\llbracket\hat{\phi}^{(e)}_{i}\rrbracket\leq\hat{C}\varepsilon$, where
$\hat{C}$ only depends on the shape of all the elements in the tessellation
but not on their size. Therefore assuming that each element in the
tessellation is convex, so that we can use the usual elliptic regularity
result, we rewrite (27) as

$\lVert
u-\hat{u}_{h}\rVert_{h}\leq\\!\\!\underset{\hat{v}_{h}\in\hat{V}_{h,\varepsilon}}{\inf}\lVert
u-\hat{v}_{h}\rVert_{h}+\sum_{K\in\mathcal{T}_{h}}\hat{C}\varepsilon\lVert
f\rVert_{L^{2}(K)}\leq\\!\\!\underset{\hat{v}_{h}\in\hat{V}_{h,\varepsilon}}{\inf}\lVert
u-\hat{v}_{h}\rVert_{h}+\hat{C}\varepsilon\lVert f\rVert_{L^{2}(\Omega)}.$
(32)

We are therefore left estimating the best approximation property for objects
that are in $\hat{V}_{h,\varepsilon}(\mathcal{T}_{h})$. If we denote
$\Pi_{k}u$ the standard VEM interpolant, [1], then we can write $\Pi_{k}u$ in
terms of the basis functions $\phi_{i}$, i.e.
$\Pi_{k}u=\sum_{K\in\mathcal{T}_{h}}\sum_{i=1}^{N_{K}}u_{i}\phi_{i}$ and
introduce
$\hat{\Pi}_{k}u=\sum_{K\in\mathcal{T}_{h}}\sum_{i=1}^{N_{K}}u_{i}\hat{\phi}_{i}$.

$\displaystyle\lVert\Pi_{k}u-\hat{\Pi}_{j}u\rVert_{h}\leq\sum_{K\in\mathcal{T}_{h}}\sum_{i=1}^{N_{K}}\lvert
u_{i}\rvert\lVert\phi_{i}-\hat{\phi}_{i}\rVert_{H^{1}(K)}\leq\sum_{K\in\mathcal{T}_{h}}\sum_{i=1}^{N_{K}}\lvert
u_{i}\rvert\hat{C}\varepsilon\leq\lVert
f\rVert_{L^{2}(\Omega)}\hat{C}\varepsilon$ (33)

Therefore using the previous estimate together with the approximation
estimates of the standard VEM interpolant we obtain

$\displaystyle\lVert u-\hat{\Pi}_{k}u\rVert_{h}$ $\displaystyle\leq\lVert
u-{\Pi}_{k}u\rVert_{h}+\lVert\Pi_{k}u-\hat{\Pi}_{k}u\rVert_{h}$ (34)
$\displaystyle\leq C(\Omega)h^{\max\\{k,m-1\\}}\lvert
u\rvert_{H^{m}(\Omega)}+\lVert f\rVert_{L^{2}(\Omega)}\hat{C}\varepsilon.$
(35)

Substituting this last estimate in the infimum appearing in the right hand
side of (27), we get

$\lVert u-\hat{u}_{h}\rVert_{h}\leq C(\Omega)h^{\max\\{k,m-1\\}}\lvert
u\rvert_{H^{m}(\Omega)}+\lVert f\rVert_{L^{2}(\Omega)}\hat{C}\varepsilon.$
(36)

##  5 Numerical experiments

To investigate the behavior of the method, we consider a family of model
problems in the unit square $\Omega=[0,1]\times[0,1]$ with as analytic
solution the function

$u(x,y)\coloneqq\sin(\pi x)\sin(\pi y)+\log(1+xy)\,.$ (37)

The meshes $\mathcal{T}_{h}$ that we will consider are centroidal Voronoi
tessellations of the unit square. An example of a Voronoi tesselation is given
in Figure 3.

Figure 3: A Voronoi tesselation made by 64 polygons.

Before presenting the cases of interest and the numerical results, we point
out that in the usual VEM framework, the errors in the $L^{2}-$norm and
$H^{1}-$seminorm are replaced by the following quantities:

$e_{H^{1}}:=\sqrt{\sum_{E\in\mathcal{T}_{h}}\left\|\nabla(u-\Pi_{k}^{\nabla}u_{h})\right\|^{2}_{0,E}}\,\qquad
e_{L^{2}}:=\sqrt{\sum_{E\in\mathcal{T}_{h}}\left\|(u-\Pi_{k}^{\nabla}u_{h})\right\|^{2}_{0,E}}\,.$
(38)

This choice is due to not knowing the analytic expression of the virtual
functions. Using lighting approximation, we can access the pointwise values of
the approximated virtual element functions. This permits to compute the local
errors using a quadrature formula without projecting into the space of
polynomials. We therefore use the usual definitions of the $H^{1}$ and $L^{2}$
errors. The source code to reproduce all the numerical experiment presented in
this section can be found in [39].

#### The Laplace problem

As first model problem, we consider the PDE

$\displaystyle-\Delta u$ $\displaystyle=f\quad\text{in }\Omega\,,$ (39)
$\displaystyle u$ $\displaystyle\equiv g\quad\text{on }\partial\Omega\,.$

We start from this equation since it represents the simplest elliptic PDE that
one can consider. Considering a sequence of Voronoi tesselations that
quadruples the number of polygons, we obtain the orders of convergence
represented in Figure 4, with $k=1$.

Figure 4: Convergence results for the Laplace problem.

We observe that we have achieved the expected order of convergence for this
equation. In particular, the error in $L^{2}$-norm and $H^{1}-$seminorm decay
as $\mathcal{O}(h^{2})$ and $\mathcal{O}(h)$, respectively.

#### The diffusion-reaction problem

We add to the previous equation a reaction term. We obtain the following PDE

$\displaystyle-\epsilon\Delta u+\gamma u$ $\displaystyle=f\quad\text{in
}\Omega\,,$ (40) $\displaystyle u$ $\displaystyle\equiv g\quad\text{on
}\partial\Omega\,,$

where $\epsilon>0$ and $\gamma\in L^{\infty}(\Omega)$ is a bounded non-
negative function. This problem is of interest because the usual VEM approach
to discretising the $L^{2}-$scalar product is to construct the
$L^{2}-$orthogonal projection operator
$\Pi^{0,E}_{k}:V_{h}^{k}(E)\to\mathbb{P}_{k}(E)$. This is not possibile with
the standard definition of the virtual element space. To overcome this
difficulty, the idea is to use the enhanced virtual element space defined as

$\displaystyle\tilde{V}_{h}(E)=\bigl{\\{}v_{h}\in H^{1}(E)\cap C^{0}(\partial
E)\quad\text{:}\quad v_{h}|_{e}\in\mathbb{P}_{k}(e)$ $\displaystyle\text{for
all $e\in\partial E$,}\bigr{.}$ (41) $\displaystyle\bigl{.}\Delta
v_{h}\in\mathbb{P}_{k}(E)\,,\,\langle
v-\Pi^{\nabla,E}_{k}v,\,\widehat{p}_{k}\rangle=0$ $\displaystyle\text{for all
$\widehat{p}_{k}\in\mathbb{P}_{k}(E)/\mathbb{P}_{k-2}(E)$}\bigr{\\}}\,.$

Thanks to the lighting approximation, we do not need to project the virtual
functions into the space of polynomials. This implies that we do not have to
change the definition of the discrete space. This gives benefits also from the
theoretical point of view. For the numerical tests, we set
$\epsilon=\gamma=1$. The results are shown in Figure 5 and we observe that
also for this equation we recovered the expected order of convergence.

Figure 5: Convergence results for the diffusion-reaction problem.

#### The advection-diffusion-reaction problem

As the last model problem, we consider the advection-diffusion-reaction
problem given by

$\displaystyle-\epsilon\Delta u+\beta\cdot\nabla u+\gamma u$
$\displaystyle=f\quad\text{in }\Omega\,,$ (42) $\displaystyle u$
$\displaystyle\equiv g\quad\text{on }\partial\Omega\,.$

where $\epsilon$ and $\gamma$ are as in the previous case and
$\beta\in[W^{1,\infty}(\Omega)]^{2}$ with $\text{div}\beta=0$. We point out
that the bilinear form associated to the advective field

$b(u,v)\coloneqq\int_{\Omega}(\beta\cdot\nabla u)\,v{\rm d}\Omega\,,$ (43)

is skew-symmetric. When we discretise $b(\cdot,\cdot)$ by inserting the
projections onto the polynomial space, we lose this property. Instead, to
preserve the skew-symmetry property, we discrtize

$b^{\text{skew}}(u,v)\coloneqq\dfrac{1}{2}\bigl{(}b(u,v)-b(v,u)\bigr{)}\,.$
(44)

Using the lighting approximation, we avoid this difficulty and do not require
(44). We select the same solution of the previous case and we set

$\beta(x,\,y):=\left[\begin{array}[]{r}-2\,\pi\,\sin(\pi\,(x+2\,y))\\\
\pi\,\sin(\pi\,(x+2\,y))\end{array}\right]\,,$

with $\epsilon=\sigma=1$. To overcome problems related to the hyperbolicity of
the advection term, we have assumed that we are not in an advection-dominated
regime. The numerical results are represented in Figure 6. In Table 1 we
compare the average local assembly time between a vanilla VEM implementation
and the lightning VEM method, clearly as the mesh becomes finer the lightning
VEM method is outperformed by the standard VEM method because we need to solve
(19) to a greater accuracy.

Figure 6: Convergence results for the advection-diffusion-reaction problem Table 1: A comparison between a vanilla VEM implementation and the lightning VEM implementation, of the average time (in seconds) taken by the assembly of the local matrix for different number of elements. N | 4 | 16 | 64 | 256 | 1024
---|---|---|---|---|---
Vanilla VEM | 4.613250e-03 | 2.031250e-03 | 2.201594e-03 | 1.102801e-03 | 1.036618e-03
Lightning VEM | 3.679250e-03 | 3.221188e-03 | 6.074953e-03 | 9.153375e-03 | 1.845604e-02
|  |  |  |  |

##  6 Extensions of the lightning virtual element method

As seen in the literature about the VEM it is possible to generalise a vast
variety of ideas developed on the FEM framework also the lightning VEM. In
this section we would like to detail some of these extensions that we plan to
investigate in the near future. The authors would like to point out that most
of these extensions are only possible thanks to the large body of work on the
lightening and AAA approximation by Nick Trefethen, Yuji Nakatsukasa and
colleagues.

  1. 1.

High order discretisation, we have focused our attention on the lowest order
VEM, i.e. the one that has only degree of freedom the nodal evaluation in the
vertices of our mesh. We would like to point out that extending the lightning
VEM idea to higher order approximation is an easy task. In particular, it is
enough to consider as (4) a modified version similar to the lowest degree
case, i.e.

$V_{h}^{k}(K)=\Big{\\{}u\in
C^{2}(K)\;:\;u_{|e}\in\mathbb{P}_{k}(e)\;\text{and}\;\Delta(u-\pi)=0\Big{\\}},$
(45)

where $\pi\in\mathbb{P}_{k}$ is determined, in a non-unique manner by the
internal moment degrees of freedom.

  2. 2.

$C^{k}$ discretisation, as mentioned in Reamark 2 it is possible to obtain
discretisation with a high degree of conformity susbtituting (4) with

$\displaystyle V_{h}^{k}(K)=\Big{\\{}u\in
C^{k+2}(K)\;:\;\partial_{n}^{j}u_{|e}\in\mathbb{P}_{2k-j}(e),\;\;$
$\displaystyle u_{|e}\in\mathbb{P}_{2k-1}(e),\;\;\Delta^{k}u=0,\;\;$ (46)
$\displaystyle j=2,\dots,p\Big{\\}}.$ (47)

In this case we can still proceed constructing a set of basis functions in the
same spirit as the lightening VEM method but following [29].

  3. 3.

Curved elements, a careful reader might have notice that an essential
requirement for the lightning virtual element method is that the tessellation
$\mathcal{T}_{h}$ is made by polygonal elements. In fact in order to transform
the non-linear problem of fitting (20) in a linear one we made the a priori
choice to cluster the poles of (20) exponentially close to the corner of the
polygonal elements. How do we choose the position of the poles if we have a
curved element? This question has been addressed in [40] and we plan to use a
similar reasoning to extend the lightning virtual element method to
tessellations with curved elements.

  4. 4.

Eigenvalue problems, it has been observed in [41, 42] that the stabilisation
term plays a harmful role in the discretisation of eigenvalue problems using a
virtual element discretisation. Possible fixes have been proposed in [20, 43].
We plan to study the role that the absence of a stabilisation term plays when
discretising an eigenvalue problem using the lightning virtual element method.

## Acknowledgments

We would like to acknowledge Astrid Herremans, Nick Trefethen, Yuji
Nakatsukasa and Patrick E. Farrell for the invaluable conversation and
continuous help while writing this paper.

## References

  *   * Beirão da Veiga et al. [2013] Veiga, L., Brezzi, F., Cangiani, A., Manzini, G., Marini, L.D., Russo, A.: Basic principles of virtual element methods. Math. Models Methods Appl. Sci. 23(1), 199–214 (2013) https://doi.org/10.1142/S0218202512500492
  * Beirão da Veiga et al. [2019] Veiga, L., Russo, A., Vacca, G.: The virtual element method with curved edges. ESAIM Math. Model. Numer. Anal. 53(2), 375–404 (2019) https://doi.org/%****␣sn-article.bbl␣Line␣75␣****10.1051/m2an/2018052
  * Beirão da Veiga et al. [2022] Beirão da Veiga, L., Dassi, F., Di Pietro, D.A., Droniou, J.: Arbitrary-order pressure-robust DDR and VEM methods for the Stokes problem on polyhedral meshes. Computer Methods in Applied Mechanics and Engineering 397, 115061 (2022) https://doi.org/10.1016/j.cma.2022.115061
  * Čertík et al. [2020] Čertík, O., Gardini, F., Manzini, G., Mascotto, L., Vacca, G.: The $p$\- and $hp$-versions of the virtual element method for elliptic eigenvalue problems. Comput. Math. Appl. 79(7), 2035–2056 (2020) https://doi.org/10.1016/j.camwa.2019.10.018
  * Dassi et al. [2020] Dassi, F., Lovadina, C., Visinoni, M.: A three-dimensional Hellinger-Reissner virtual element method for linear elasticity problems. Comput. Methods Appl. Mech. Engrg. 364, 112910–17 (2020) https://doi.org/10.1016/j.cma.2020.112910
  * Beirão da Veiga et al. [2013] Veiga, L., Brezzi, F., Marini, L.D.: Virtual elements for linear elasticity problems. SIAM J. Numer. Anal. 51(2), 794–812 (2013) https://doi.org/10.1137/120874746
  * Artioli et al. [2017] Artioli, E., Miranda, S., Lovadina, C., Patruno, L.: A stress/displacement virtual element method for plane elasticity problems. Comput. Methods Appl. Mech. Engrg. 325, 155–174 (2017) https://doi.org/10.1016/j.cma.2017.06.036
  * Antonietti et al. [2014] Antonietti, P.F., Veiga, L., Mora, D., Verani, M.: A stream virtual element formulation of the Stokes problem on polygonal meshes. SIAM J. Numer. Anal. 52(1), 386–404 (2014) https://doi.org/%****␣sn-article.bbl␣Line␣175␣****10.1137/13091141X
  * Beirão da Veiga et al. [2018] Veiga, L., Lovadina, C., Vacca, G.: Virtual elements for the Navier-Stokes problem on polygonal meshes. SIAM J. Numer. Anal. 56(3), 1210–1242 (2018) https://doi.org/10.1137/17M1132811
  * Beirão da Veiga et al. [2017] Veiga, L., Lovadina, C., Vacca, G.: Divergence free virtual elements for the Stokes problem on polygonal meshes. ESAIM Math. Model. Numer. Anal. 51(2), 509–535 (2017) https://doi.org/10.1051/m2an/2016032
  * Antonietti et al. [2016] Antonietti, P.F., Veiga, L., Scacchi, S., Verani, M.: A $C^{1}$ virtual element method for the Cahn-Hilliard equation with polygonal meshes. SIAM J. Numer. Anal. 54(1), 34–56 (2016) https://doi.org/10.1137/15M1008117
  * Brezzi and Marini [2013] Brezzi, F., Marini, L.D.: Virtual element methods for plate bending problems. Comput. Methods Appl. Mech. Engrg. 253, 455–462 (2013) https://doi.org/10.1016/j.cma.2012.09.012
  * Zhao et al. [2016] Zhao, J., Chen, S., Zhang, B.: The nonconforming virtual element method for plate bending problems. Math. Models Methods Appl. Sci. 26(9), 1671–1687 (2016) https://doi.org/10.1142/S021820251650041X
  * Beirão da Veiga et al. [2017] Veiga, L., Mora, D., Rivera, G., Rodríguez, R.: A virtual element method for the acoustic vibration problem. Numer. Math. 136(3), 725–763 (2017) https://doi.org/10.1007/s00211-016-0855-5
  * Perugia et al. [2016] Perugia, I., Pietra, P., Russo, A.: A plane wave virtual element method for the Helmholtz problem. ESAIM Math. Model. Numer. Anal. 50(3), 783–808 (2016) https://doi.org/10.1051/m2an/2015066
  * Beirão da Veiga et al. [2018a] Veiga, L., Brezzi, F., Dassi, F., Marini, L.D., Russo, A.: Lowest order virtual element approximation of magnetostatic problems. Comput. Methods Appl. Mech. Engrg. 332, 343–362 (2018) https://doi.org/10.1016/j.cma.2017.12.028
  * Beirão da Veiga et al. [2018b] Veiga, L., Brezzi, F., Dassi, F., Marini, L.D., Russo, A.: A family of three-dimensional virtual elements with applications to magnetostatics. SIAM J. Numer. Anal. 56(5), 2940–2962 (2018) https://doi.org/10.1137/18M1169886
  * Berrone et al. [2023] Berrone, S., Borio, A., Marcon, F., Teora, G.: A first-order stabilization-free Virtual Element Method. Applied Mathematics Letters 142, 108641 (2023) https://doi.org/10.1016/j.aml.2023.108641
  * Berrone et al. [2022] Berrone, S., Borio, A., Marcon, F.: Comparison of standard and stabilization free Virtual Elements on anisotropic elliptic problems. Applied Mathematics Letters 129, 107971 (2022) https://doi.org/10.1016/j.aml.2022.107971
  * Meng et al. [2022] Meng, J., Wang, X., Bu, L., Mei, L.: A lowest-order free-stabilization Virtual Element Method for the Laplacian eigenvalue problem. Journal of Computational and Applied Mathematics 410, 114013 (2022) https://doi.org/10.1016/j.cam.2021.114013
  * Beirão da Veiga et al. [2023] Veiga, L., Canuto, C., Nochetto, R.H., Vacca, G., Verani, M.: Adaptive VEM: Stabilization-Free A Posteriori Error Analysis and Contraction Property. SIAM Journal on Numerical Analysis 61(2), 457–494 (2023) https://doi.org/10.1137/21M1458740
  * [22] Lehrenfeld, C., Stocker, P.: Embedded Trefftz discontinuous Galerkin methods. International Journal for Numerical Methods in Engineering n/a(n/a) https://doi.org/10.1002/nme.7258
  * Qin [2005] Qin, Q.-H.: Trefftz Finite Element Method and Its Applications. Applied Mechanics Reviews 58(5), 316–337 (2005) https://doi.org/10.1115/1.1995716 https://asmedigitalcollection.asme.org/appliedmechanicsreviews/article-pdf/58/5/316/5441315/316_1.pdf
  * Antonietti et al. [2020] Antonietti, P.F., Manzini, G., Verani, M.: The conforming virtual element method for polyharmonic problems. Comput. Math. Appl. 79(7), 2021–2034 (2020) https://doi.org/10.1016/j.camwa.2019.09.022
  * Antonietti et al. [[2022] ©2022] Antonietti, P.F., Manzini, G., Mazzieri, I., Scacchi, S., Verani, M.: The conforming virtual element method for polyharmonic and elastodynamics problems: a review. In: The Virtual Element Method and Its Applications. SEMA SIMAI Springer Ser., vol. 31, pp. 411–451. Springer, ??? ([2022] ©2022). https://doi.org/10.1007/978-3-030-95319-5_10
  * Gopal and Trefethen [2019] Gopal, A., Trefethen, L.N.: Solving Laplace problems with corner singularities via rational functions. SIAM J. Numer. Anal. 57(5), 2074–2094 (2019) https://doi.org/10.1137/19M125947X
  * Trefethen [2020] Trefethen, L.N.: Numerical conformal mapping with rational functions. Comput. Methods Funct. Theory 20(3-4), 369–387 (2020) https://doi.org/10.1007/s40315-020-00325-w
  * Gopal and Trefethen [2019] Gopal, A., Trefethen, L.: New laplace and helmholtz solvers. Proceedings of the National Academy of Sciences 116 (2019)
  * Brubeck and Trefethen [2022] Brubeck, P.D., Trefethen, L.N.: Lightning Stokes solver. SIAM J. Sci. Comput. 44(3), 1205–1226 (2022) https://doi.org/10.1137/21M1408579
  * Zhu and Nakatsukasa [2023] Zhu, W., Nakatsukasa, Y.: Convergence and Near-optimal Sampling for Multivariate Function Approximations in Irregular Domains via Vandermonde with Arnoldi. arXiv preprint arXiv:2301.12241 (2023)
  * Gilardi [2020] Gilardi, G.: Analisi 3. McGraw-Hill Education, ??? (2020)
  * Ahlfors and Collection [1979] Ahlfors, L.V., Collection, K.M.R.: Complex Analysis: An Introduction to The Theory of Analytic Functions of One Complex Variable. International series in pure and applied mathematics. McGraw-Hill Education, ??? (1979)
  * Grisvard [2011] Grisvard, P.: Elliptic Problems in Nonsmooth Domains. Classics in Applied Mathematics, vol. 69, p. 410. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, ??? (2011). https://doi.org/10.1137/1.9781611972030.ch1 . Reprint of the 1985 original [ MR0775683], With a foreword by Susanne C. Brenner
  * Gilbarg and Trudinger [2001] Gilbarg, D., Trudinger, N.S.: Elliptic Partial Differential Equations of Second Order. Classics in Mathematics, p. 517. Springer, ??? (2001). Reprint of the 1998 edition
  * Evans [2010] Evans, L.C.: Partial Differential Equations, 2nd edn. Graduate Studies in Mathematics, vol. 19, p. 749. American Mathematical Society, Providence, RI, ??? (2010). https://doi.org/10.1090/gsm/019
  * Herremans et al. [2023] Herremans, A., Huybrechs, D., Trefethen, L.N.: Resolution of singularities by rational functions (2023)
  * Brubeck et al. [2021] Brubeck, P.D., Nakatsukasa, Y., Trefethen, L.N.: Vandermonde with arnoldi. SIAM Review 63(2), 405–415 (2021) https://doi.org/10.1137/19M130100X https://doi.org/10.1137/19M130100X
  * Brenner and Scott [2007] Brenner, S., Scott, R.: The Mathematical Theory of Finite Element Methods vol. 15. Springer, ??? (2007)
  * Trezzi and Zerbinati [2023] Trezzi, M.L., Zerbinati, U.: LightningVEM. GitHub (2023). https://doi.org/10.5281/zenodo.8192892 . https://github.com/UZerbinati/LightningVEM
  * Xue et al. [2023] Xue, Y., Waters, S.L., Trefethen, L.N.: Computation of 2D Stokes flows via lightning and AAA rational approximation (2023)
  * Boffi et al. [2022] Boffi, D., Gardini, F., Gastaldi, L.: In: Antonietti, P.F., Veiga, L., Manzini, G. (eds.) Virtual Element Approximation of Eigenvalue Problems, pp. 275–320. Springer, Cham (2022). https://doi.org/10.1007/978-3-030-95319-5_7
  * Boffi et al. [2020] Boffi, D., Gardini, F., Gastaldi, L.: Approximation of PDE eigenvalue problems involving parameter dependent matrices. Calcolo 57(4), 41 (2020) https://doi.org/10.1007/s10092-020-00390-6
  * Gardini and Vacca [2017] Gardini, F., Vacca, G.: Virtual element method for second-order elliptic eigenvalue problems. IMA Journal of Numerical Analysis 38(4), 2026–2054 (2017) https://doi.org/10.1093/imanum/drx063