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

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.


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 structurepreserving 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 socalled 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 costefficient 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 pointwise values, and "virtual" refers to the fact that we are approximating the standard virtual element space.

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 Ω ⊂ R 2 with boundary Γ, we consider the problem of finding a function u : Ω → R such that where f ∈ L 2 (Ω) is the load term.It is well known that the standard weak formulation of (1) reads as find u ∈ where the bilinear form a(•,

The local space
We decompose the domain Ω in a tessellation T h by a finite number of non-overlapping convex polygons K.In particular, we assume that there exists a positive constant ρ such that every element K ∈ T h is star-shaped with respect to a ball of a radius greater or equal than ρh K , where h K denotes the diameter of the element K. Let k ≥ 1 be the "order" of the method.For every element K ∈ T h we define the local virtual element space as where where e denotes an edge of ∂K and we have denoted P −1 := {0}.One can check that the following quantities represent a set of degrees of freedom for the space V k h (K): 1. V K,k : the pointwise values of v h at the vertexes of the polygon K, 2. E K,k : the values of v h at k − 1 internal points of a Gauss-Lobatto quadrature for every edge e ∈ ∂K, is the set of monomials defined as Details of the proof are in [1].A graphic representation of the degrees of freedom is given in Figure 2.1.Thanks to the definition of the degrees of freedom, it is possible to construct the projection operator Π ∇,K Indeed, thanks to integration by parts, we note that The first integral is known thanks to P K,k , for the second one we use E K,k and V K,k .The operator Π ∇,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(•, •).The global space is obtained by gluing together the local spaces therefore obtaining a space characterized by the following degrees of freedom: 1. V k : the values of v h at the vertices; 2. E k : the values of v h at k − 1 points on each edge e; 3. P k : the moments up to order k − 2 for each element E ∈ T h .

The discrete problem
First, we decompose the global bilinear form into local contributions 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 k h (K).Hence, given two generic virtual functions, we are not able to compute the quantity We would like to construct a discrete bilinear form a K h : that is computable for all the virtual functions and acts as a discrete counterpart of a K (•, •).The idea is to split the virtual functions as Thanks to this choice, the bilinear form is split as Thanks to orthogonality of Π ∇,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 for two positive uniform constants α * and α * .This term is called the stability term.We define the bilinear form The global bilinear form is obtained summing all the local contributions It remains to discretise the load term.A standard choice is the following procedure: 1. if k ≥ 2, we replace f with f h defined locally as the L 2 -orthogonal projection into the space P k−2 (K).In detail, we consider 2. if k = 1, we replace f with a piecewise constant and define where The discrete problem reads as 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 P k is a subset of the co-domain of chosen the operator when applied to the space P k−2 with the constraint of having boundary data in 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 −∆ 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-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 V k and 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.

The lightning approximation
The key idea behind the lightning VEM is to consider a different local discrete space which is obtained by approximating V k h (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 k h (K) we are interested in solving N K problems of the form where φ i are the functions in B 1 (∂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 φi of ϕ i of the form φi = ℜ where {z j } N P j=1 and z * are points in the complex plane and ℜ denotes the real part of a complex number.Finding the optimal coefficients, {z j } N P j=1 and z * to minimize ∥ φi − ϕ i ∥ L ∞ (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 } N P j=1 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 , . . ., 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 Under these assumptions, there exists a rational function r of the form 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: 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 ∞ 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 φi to (19) is presented in Algorithm 1 and we once again direct the reader interested in more detail to [26].We then observe that Algorithm 1 For computing the solution φ(e) i of (19).
Fix N P ≈ nN K and cluster the poles outside K and exponentially close the corners;

3:
Fix z * in the interior of K and fix N Z ≈ N K elements of the monomial basis;

4:
Choose M = 6N P + 6N Z + 1 sample points on the boundary and clustered near the corners;

5:
Evaluate at the sample points φ i to obtain a matrix A ∈ R M ×(2N P +2N Z +1) and d ∈ R M ; N Z j=0 b j (z − z * ) j is analytic in the convex hull described by the vertex of any polygonal element K of the tessellation 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 φi is harmonic we know that the function defined as φi − ϕ i is harmonic and by applying boundary regularity estimates for harmonic functions we obtain the following result.Proposition 2. Let ϕ i be as in (19) and φi be the outcome of Algorithm 1. Then where the constants C 1 and C 2 only depend on Ω and ∂Ω.
Fig. 2 We know that the basis function φi,K 1 and φ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 φ1 and φ2 are continuous along the blue edge.
Proof.Since we assumed the elements K of our tessellation 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 Vh,ε (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 Vh,ε (T h ) to denote the global space constructed starting from the various Vh,ε (K).Remark 4. It is clear that a fundamental step in the lightning VEM is solving the least squares problem Ax ≈ d 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 x.In fact bounding the 2-norm of x will guarantee the solvability of (24) to a high degree of accuracy in spite its illconditioned 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 [29,30].

The lightning virtual element method
We begin by observing that the degrees of freedom V h and E h are not decoupled in adjacent elements.In fact if φ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 φi across the elements interface.Instead, since each φi is a rational function on the edges of the element K even if over adjacent elements φ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 This bilinear form on H 1 (Ω) is still symmetric positive definite with respect to the broken H 1 (Ω) norm and for the continuous solution to the problem we have that â(u, v) = a(u, v) since u, v ∈ H 1 (Ω).Furthermore, in Vh,ε (T h ) the kernel of the bilinear form â(•, •) 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 â(•, •) has a trivial kernel.We now derive an a priori error estimate using the following Lemma from [38].Lemma 3. Let â(•, •) be a symmetric positive definite bilinear form on H 1 (Ω) + Vh,ε (T h ) which reduces to a(•, •) on H 1 (Ω).We will further assume that u is the solution of (2) and ûh is such that: where in both the previous variation equation and (2) Under this hypothesis the following best approximation estimate holds, where 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). = where E h are the internal edges of the tessellation T h and ŵh denotes the vector jump of ŵh across the edge e, i.e.
We now rewrite the last term using the basis functions φi e , where the last inequality comes from the trace theorem combined with Cauchy-Schwarz inequality and the fact that we have constructed the basis function φ(e) i in such a way that φ(e) i ≤ Ĉε, where Ĉ 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 We are therefore left estimating the best approximation property for objects that are in Vh,ε (T h ).If we denote Π k u the standard VEM interpolant, [1], then we can write Π k u in terms of the basis functions ϕ Therefore using the previous estimate together with the approximation estimates of the standard VEM interpolant we obtain Substituting this last estimate in the infimum appearing in the right hand side of ( 27), we get
The meshes T h that we will consider are centroidal Voronoi tessellations of the unit square.An example of a Voronoi tesselation is given in Figure 5. 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: 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 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.
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 O(h 2 ) and O(h), respectively.

The diffusion-reaction problem
We add to the previous equation a reaction term.We obtain the following PDE where ϵ > 0 and γ ∈ L ∞ (Ω) 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 Π 0,E k : V k h (E) → 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 (41) 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 ϵ = γ = 1.The results are shown in Figure 5 and we observe that also for this equation we recovered the expected order of convergence.

The advection-diffusion-reaction problem
As the last model problem, we consider the advection-diffusion-reaction problem given by −ϵ∆u 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 β(x, y) := −2 π sin(π (x + 2 y)) π sin(π (x + 2 y)) , with ϵ = σ = 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.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. 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.
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. Curved elements, a careful reader might have notice that an essential requirement for the lightning virtual element method is that the tessellation 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. 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.

Fig. 1
Fig. 1 Degrees of freedom on a pentagon.

Fig. 4
Fig. 4 Convergence results for the Laplace problem.

Fig. 5
Fig. 5 Convergence results for the diffusion-reaction problem.whereϵ and γ are as in the previous case and β ∈ [W 1,∞ (Ω)] 2 with divβ = 0. We point out that the bilinear form associated to the advective field
. It is important to notice that this jump is bounded by adaptively solving the lightning approximation problem in Algorithm 1, i.e. ∥ φi − ϕ i ∥ H 1/2 (∂K) ≤ ε.Since the function in Vh,ε (T h ) are no longer continuous we need to consider a broken version of the bilinear form a in (3), i.e. â(û h , vh ) = K (û h , vh )

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.