Ill-conditioning in the Virtual Element Method: stabilizations and bases

In this paper we investigate the behavior of the condition number of the stiffness matrix resulting from the approximation of a 2D Poisson problem by means of the Virtual Element Method. It turns out that ill-conditioning appears when considering high-order methods or in presence of"bad-shaped"(for instance nonuniformly star-shaped, with small edges...) sequences of polygons. We show that in order to improve such condition number one can modify the definition of the internal moments by choosing proper polynomial functions that are not the standard monomials. We also give numerical evidence that, at least for a 2D problem, standard choices for the stabilization give similar results in terms of condition number.


INTRODUCTION
The interest in numerical methods for the approximation of partial differential equations (PDEs in short) based on polytopic grids has grown in the last decade, due to the high flexibility that polygonal/polyhedral meshes allow. Among the other methods, we here recall only a short list including: mimetic finite differences [1,2], discontinuous Galerkin-finite element method (DG-FEM) [3,4], hybridizable and hybrid high-order methods [5,6], weak Galerkin method [7], BEM-based FEM [8], and polygonal FEM [9].
An alternative approach is offered by the virtual element method (VEM in short), recently introduced in [10]. VEM are a generalization of the finite element method (FEM in short) enabling the employment of polytopal meshes and the possibility of building high-order methods. In addition to polynomials, local VE spaces consist of other functions instrumental for constructing global H 1 conforming approximation space; such functions are typically the solution to local PDEs and therefore are not known explicitly in a closed-form. For this reason, VEM can be considered, for all practical purposes, Trefftz methods.
In VEM, the bilinear forms and the right-hand sides are not computed exactly due to the fact that the functions in VE spaces are not fully explicit. The VEM gospel states that such bilinear forms and right-hand sides are then approximated by means of discrete counterparts that are exactly computable only through the degrees of freedom and which scale like the continuous ones.
It was observed that the VEM stiffness matrix, similarly as FEM, can become ill-conditioned in various situations. This is the case for instance of the p and the hp version of VEM, as discussed in [11,13], and in presence of "badly-shaped" (i.e., nonuniformly star-shaped, with small edges...) polygons, as discussed [34].
Among the possible reasons of this ill-conditioning we highlight two of them. The first one is related to the fact that in the VEM framework one does not use the exact bilinear form but an approximated one; the choice of the discrete bilinear form, and, in particular, of one of its two components, namely the stabilization of the method, may have an impact on the 3D version of VEM as observed in [33]. The second one is the choice of the basis. This is also the case for FEM, where the choice of the basis has an important role on the ill-conditioning of the system, see [35][36][37] and the references therein.
The aim of the present paper is to discuss various choices for both the discrete bilinear forms and the VE bases and check numerically that particular choices can cure the ill-conditioning which arises in high-order (or in presence of badly-shaped polygons) VEM.
In particular, we show that while the choice of the stabilization has not a deep impact on the illconditioning (at least for the 2D case, which is the focus of this article), a proper choice of the basis can actually improve the condition number of the stiffness matrix. However, it is worth to mention that in "standard" situations, that is, low-to-moderate order VEM and VEM applied to shape-regular decompositions, the standard VEM (e.g., the one described in [10]) is preferable, as the implementation of that version of the method turns out to be much simpler than those we are going to present in this article.
The outline of the article is the following. We firstly discuss the model problem and its VEM approximation in Section 2, while in Appendix A, we give a hint on the implementation details of the method using one of the new VEM basis. In Section 3, we present a number of numerical experiments comparing the behavior of the method when changing various stabilizations and VE bases.
In the remainder of the article, we adopt the standard notation for Sobolev spaces, see for example, [38,39]. Given ω ⊂ R 2 , we denote by H (ω), ∈ N, the Sobolev space of order over ω; in the case = 0, we set H 0 (ω) = L 2 (ω), where L 2 (ω) is the Lebesgue space of square integrable functions over ω. By H 1 0 (Ω), we mean the Sobolev space H 1 of functions with zero trace. The Sobolev (semi)norms and inner products read, respectively: Further, given a Lipschitz domain ω with boundary Γ, we set n the normal versor associated with Γ and ∂ n v = ∇v · n the normal derivative of a sufficiently regular function v. By P (ω), ∈ N, we denote the space of polynomials of degree over ω. Finally, given two positive quantities a and b, possibly depending on the discretization parameters h and p, we write a b if there exists a positive and parameter-independent constant c such that a ≤ cb. Moreover, we write a ≈ b meaning that a b and b a simultaneously.

THE VIRTUAL ELEMENT METHOD: DEFINITION, CHOICE OF THE STABILIZATION AND OF THE BASIS
Given Ω ⊂ R 2 a polygonal domain and f ∈ L 2 (Ω), we consider the following 2D Poisson problem: The outline of the present section is the following. In Section 2.1, we introduce a family of VEMs for the approximation of problem (2); here, both the stabilization, typical of VEM, and the choice of the basis of local VE space are kept at a very general level. Explicit choices for the stabilization are investigated in Section 2.2, whereas explicit choices for the local basis elements are the topic of Section 2.3. Finally, in Section 2.4, we highlight the influence of the choices of the stabilizations and of the local basis elements on the performances of the method.

A family of VEM
Given Ω the computational domain of problem (2), we consider a family {T n } n∈N of conforming polygonal decomposition of Ω, where by conforming we mean that, given an edge in the skeleton of the decomposition which does not lie on ∂Ω, then it is an edge of exactly two polygons. We also fix p ∈ N, which will represent the degree of accuracy of the method. We associate to each K ∈ T n its diameter h K and to decomposition T n its mesh size function h = max K∈Tn h K . Standard regularity assumptions on decomposition T n that are usually required in VEM literature read: • For every K ∈ T n , K is star-shaped (see [40]) with respect to a ball of radius greater or equal than γh K , γ being a positive constant independent of the family of decompositions.
• For every K ∈ T n and for every e edge of K, the length of e is greater or equal than γh K , γ being the same constant introduced in assumption (D1).
We point out that such assumptions can be relaxed, see [41]. We now define the local VE space on polygon K ∈ T n following the standard definition given, for example, in [10]. Having set the space of piecewise continuous polynomials of degree p over ∂K: we introduce the local VE space as: Space V p (K), which contains P p (K), contains also other functions that in general are not known pointwise (hence, the name virtual), but which are added to polynomials in order to guarantee the possibility of building a H 1 conforming method over T n .
We associate to space V p (K) the following set of linear functionals. Given v p ∈ V p (K): • the point values of v p at the vertices of K; • for every e edge of K, the point values of v p at the p − 1 internal Gauß-Lobatto nodes of e; • the (scaled) internal moments: where is any basis of P p−2 (K).
It was proven in [10] that this is a set of unisolvent degrees of freedom. Let N K dof be the number of such degrees of freedom, that is, the dimension of space V p (K). Henceforth, we denote by dof i the i-th dof of space V p (K) and by: , that is, the set of functions such that: Functions in the canonical basis associated with internal moments (5) are bubbles on element K as they vanish on the boundary. It is fundamental to observe that the definition of space V p (K) is completely independent of the choice of polynomial basis {q α } dim(P p−2 (K)) α=1 , which is, so far, only instrumental for the definition of the internal moments (5). Nonetheless, such a choice plays a crucial role in the behaviour of the condition number of the stiffness matrix of the method as discussed and shown in Section 3. In Subsection 2.3, we present many choices of basis {q α } p−2 |α|=0 . As already stressed, not all the functions in space V p (K) are known explicitly; nevertheless, it is possible to compute local H 1 and L 2 projections as shown in [10,31] via the degrees of freedom above described. More precisely, it is possible to compute the following two operators.
The first one is the L 2 projector Π 0 p−2 : V p (K) → P p−2 (K) defined, for all K ∈ T n , by: which is clearly computable via internal moments (5). The second one is a H 1 projector Π ∇ p : V p (K) → P p (K) defined by: where a K (·, ·) = (∇·, ∇·) 0,K and where P is an operator which fixes the constant part of the energy projector Π ∇ p defined as: where ν j N K j=1 is the set of vertices of polygon K. Operator Π ∇ p is computable by means of the degrees of freedom, noting that: and recalling that v p ∈ V p (K).
At this point, we turn our attention to the computation of the local stiffness matrix and of the local right-hand side. In both cases, we cannot use their continuous counterparts as we do not know pointwise functions in space V p (K). Therefore, we follow the VEM gospel and we note that Pythagoras Theorem for Hilbert spaces asserts: If we split the local space V p (K) into a polynomial and a "pure virtual" part: then the first term on the right-hand side of (11) is identically zero on V p,2 (K), while the second one annihilates on V p,1 (K). We point out that, given N K the number of vertices of K, one has: This entails that the actual "pure virtual" part of space V p (K), that is, ker(Π ∇ p ), is asymptotically smaller than its polynomial counterpart, if the number of vertices of K remains uniformly bounded. More precisely, dim(V p,2 (K)) ≈ p whereas dim(V p,1 (K)) ≈ p 2 .
We emphasize that the first term on the right-hand side of (11) is explicitly computable but the second one is not. For this reason, one substitutes: where S K (·, ·) is any bilinear form mimicking the continuous one, that is, a K (·, ·) on ker(Π ∇ p )×ker(Π ∇ p ). To have a well-posed method, we demand the following stability assumption on the bilinear form S K : where c * (p) and c * (p) are two positive constants possibly depending on p. Various possible stabilizations S K are presented in Section 2.2. By defining the local discrete bilinear form as: one can prove the following properties for a K p (·, ·): (P1) p-consistency: for every q p ∈ P p (K) and for every v p ∈ V p (K), the local bilinear form a K p satisfies: (P2) stability: for every v p ∈ V p (K), the local bilinear form a K p satisfies: where α * (p) = min(1, c * (p)) and α * (p) = max(1, c * (p)).
Regarding the local discrete right-hand side, we set: At this point, we are able to define the global VE space, which is obtained by a standard dof coupling of the local spaces: and the discrete global bilinear form and right-hand side: We now define a family of VEMs depending on the choice of the local stabilization S K defined in (14): Property (P1) guarantees that the patch test is satisfied by any solution with a polynomial degree smaller than or equal to p (the polynomial accuracy of the method). Conversely, property (P2) implies the coercivity and the continuity of the discrete bilinear form and thus the well-posedness of methods (21). It is worth to stress that the coercivity and continuity constants may depend on p, owing to (17).
Let us set the H 1 broken Sobolev seminorm associated with decomposition T n as: and let us set S p,−1 (Ω, T n ) the space of piecewise discontinuous polynomials of degree p over decomposition T n . Having properties (P1) and (P2), one can prove as in [11,12] an abstract error analysis result. Given u and u p the solutions to (2) and (21), respectively, the following holds true: where F n is the smallest positive constant such that Thus, being able to study the convergence of the method is equivalent to being able to prove best approximation results by piecewise discontinuous polynomials in S p,−1 (Ω, T n ) and by functions in V p , the global VE space defined in (19). In [11,12], it was proven that for any u, solution to (2), in H k+1 (Ω), the following hp a priori estimate holds true: and that for u analytic on a proper extension of Ω, the following pure p approximation result holds true: where b is a positive constant independent of the discretization parameters. The dependence on p of the pollution factor α * (p) α * (p) is clearly an issue of extreme importance in the p analysis of the method. The following crude upper bound: using a particular stabilization, was proved in [12,Theorem 2]. However, numerical experiments on different stabilizations highlight that the dependence on p is much milder; typically, even on nonconvex and nonstar-shaped polygons, one has α * (p) α * (p) p β , for some β ∈ (0, 1); see [12,11]. It is worth to stress that such algebraic dependence on the degree of accuracy does not affect the exponential convergence (25) when approximating analytic solution. Importantly, it is possible to get rid of the dependence on the pollution factor also when approximating solutions presenting corner singularities, using the hp version of the method; see [12].

Choices for the stabilization
Here, we address the issue of choosing a proper stabilization S K (14) for method (21). In particular, we define four possible stabilizations which can be found in literature and we discuss their properties. We recall that N K dof is the number of dofs of local space V p (K).
The first stabilization that we present is: which can be regarded as the "standard" stabilization of VEM. It was firstly introduced in [10] and it can be proved that, for a fixed degree of accuracy p, S K 1 scales like the H 1 seminorm. We highlight that this holds true only in two dimension; for a three dimensional problem one should put a proper scaling factor in front of S K 1 , see [33,32]. The advantage of picking S K 1 as a stabilization is that its implementation is extremely easy.
If one takes into account also variable p, then one has to prove the dependence on p of the stability constant in (14). This was done in [12], where the authors introduced the following computable stabilization: At the current theoretical level, the stabilization constants c * (p) and c * (p) which appear in (14) have a dreadful dependence on p; nonetheless, at the practical level, it was observed in [12] an extremely mild dependence on p of c * (p) and c * (p). Next, we recall another stabilization which was introduced in [33]. If we denote by K K C the consistency part of the local stiffness matrix, that is, the matrix counterpart of the first term on the right-hand side of (15), then we can define a stabilization S K 3 through its matrix representation S K 3 associated with the second term on the right-hand side of (15) as follows. S K 3 is a diagonal matrix, whose i-th entry is given by: For the sake of clarity, the complete local stiffness matrix K K n reads: where Π ∇ is the matrix representing the action of operator Π ∇ p and S K 3 is defined in (28). In [33], numerical experiments, along with a heuristic motivation, show that choice (28) entails a better behaviour of the method for high p, when approximating a 3D Poisson problem.
Finally, we also consider a fourth stabilization. Let N bndr dof be the number of boundary degrees of freedom of V p (K). Then, we set: which is basically stabilization S K 1 without the contribution of internal dofs. We point out that other choices of stabilization can be found in literature, but we prefer to investigate the four presented here in order to avoid an overwhelming set of numerical experiments.
We also emphasize that the first three stabilizations above hinge on the choice of the polynomial basis dual to internal moments (5).

Choices for the basis
In this section, we discuss three possible choices of the local VE basis. More precisely, we consider internal moments (5) taken with respect to three different polynomial bases q i α dim(P p−2 (K)) α=1 , i = 1, 2, 3, of P p−2 (K), thus modifying the definition of the VE bubble functions.
The hope is that a proper choice of the polynomial basis dual to internal moments (5), and therefore, of internal VEM basis elements, entails a better conditioning for the stiffness matrix.
Henceforth, we will use, with an abuse of notation, the natural bijection N 2 ↔ N 0 given by: We will also write occasionally: , for i = 1, 2, 3.
The first choice of the polynomial basis is the "standard" one, that is, the one which is used in the majority of VEM literature, since, at the implementation level, is the most convenient. Given x K = (x K , y K ) the barycenter of K, we set: Although choice (31) is very suitable from the computational point of view, it has bad effects on the condition number of the stiffness matrix for high local degrees of accuracy p, see [11], and in presence of badly-shaped polygons, see [34] and the references therein. For this reason, we suggest two possible modifications which rely on orthogonalizing processes of q 1 α p−2 |α|=0 with respect to the L 2 norm on polygon K. The first modification, which allows to construct an L 2 (K) orthonormal basis q 2 α p−2 |α|=0 , is based on the stable Gram-Schmidt orthonormalization process presented in [42]. We point out that basis q 2 α p−2 |α|=0 was firstly introduced in the context of p-VEM multigrid algorithm for the construction of the multrigrid scheme, see [13].
While the implementation of VEM with basis q 1 α p−2 |α|=0 is well-known (see [31]), no details were given in [13] for the implementation of the method using basis q 2 α p−2 |α|=0 . Therefore, we address this issue in Appendix A.
The third choice of the polynomial basis is slightly inspired by an orthonormalizing procedure used in [34]. To present the third basis q 3 α p−2 |α|=0 , we set matrix: Given n = dim(P (K)), we fix q 3 α = q 1 α if α = 0 and we decompose matrix H into blocks as: we diagonalize matrix H 2,2 and we get: This entails that matrix V·D − 1 2 contains the coefficients which orthonormalize q 1 α p−2 |α|=1 , the monomial basis of P p−2 (K)/R. Therefore, one has: Importantly, the basis in (32) is not completely orthogonal, as the orthogonalization process is performed on the block H 2,2 only. It is worth to stress that here the orthonormalizing process is performed with a different target with respect to what done in [34]; in fact, there, the method is built using the canonical bases computed taking moments with respect to scaled monomials q 1 In summary, we presented three choices for the polynomial basis dual to internal moments (5). One is of easy implementation, but it may be the cause of a high condition number for the stiffness matrix when using high-order methods or in presence of "badly-shaped" polygons. The other two bases are obtained by two distinct orthonormalization processes; their performances with respect to the condition number are investigated in Section 3. What we can anticipate is that they outclass the performances of their counterpart using the standard basis q 1 α p−2 |α|=0 in the two situations above mentioned. A heuristic reason for this fact is the following. If, for each element K, the local virtual element space V p (K) were a space consisting of polynomials only, then picking internal moments (5) with respect to an L 2 (K) orthonormal polynomial basis would automatically entail that the local canonical basis is made of polynomials and contain a subset of L 2 (K) orthonormal polynomials spanning P p−2 (K); it is well-known in the theory of Spectral Elements, see for example, [43,44], that using L 2 orthonormal canonical basis damp the condition number of the stiffness matrix when increasing the polynomial degree.
Nonetheless, local virtual element spaces does not contain polynomials only, but also other functions needed for prescribing H 1 conformity. As stated in (13), the dimension of the subspace of nonpolynomial functions is, in terms of the degree of accuracy, asymptotically smaller than the dimension of the subspace of polynomial functions. Therefore, in a very rough sense, using L 2 orthonormal polynomials in the definition of internal moments entails a sort of partial "orthonormalization" of the local canonical basis.
Before concluding this section, we associate to bases q i , for all i = 1, 2, 3. This notation will be instrumental in Appendix A.

Stabilizations and bases: the effects on the method
Having presented in the two foregoing sections various choices of stabilizations and canonical bases, we want here to highlight the effects of such choices on the method and on the ill-conditioning of the stiffness matrix.

Stabilization
The choice of the stabilization has two effects. The first one is related to the convergence of the method since nonproperly tailored choices of the stabilization automatically entail higher pollution factor α * (p) α * (p) ,  (25) and especially (24). Second, as the stabilization appears in the discrete bilinear form of the method, see (15), there is also an effect on the condition number of the stiffness matrix.

Canonical basis
The choice of the canonical basis has also two effects. First, it has an impact on the condition number of the global stiffness matrix, simply because by changing the basis automatically the entries of the stiffness matrix modify. Second, by picking different canonical bases, one also changes the definition of the stabilization; as an example, if we fix stabilization S K 1 defined in (26) and we apply it to functions in the virtual element space, then we get in general different values, as the definition of the internal degrees of freedom vary depending on the choice of the basis; in particular, one also modify the behaviour of the pollution factor α * (p) α * (p) .

NUMERICAL RESULTS
In this section, we present some numerical experiments in which we compare the performances of the stabilizations and polynomial bases introduced in Sections 2.2 and 2.3, respectively. More precisely, we investigate the behaviour and the related effects of the condition number in two critical situations.
In Section 3.1, we investigate the behavior of the condition number of the p version of VEM and the effects on the linear solver used for the resolution of the associated linear system. We are also interested in the behavior of the condition number when varying the stabilization and the polynomial basis dual to internal moments (5) in presence of a sequence of "bad shaped" polygons (collapsing bulk, collapsing edges...); this is probed in Sections 3.2 and 3.3. The condition number is computed as the ratio between the maximum and the minimum (nonzero) eigenvalue of the stiffness matrix.

Numerical results: the p version of VEM
Let us consider the three meshes depicted in Figure 1. We investigate in this section the behaviour of the condition number of the stiffness matrix associated with method (21) by keeping fixed the meshes of Figure 1 and by increasing p. For the purpose, we modify the choice of the stabilization and the choice of the polynomial basis dual to internal moments (5).
In Figure 2, we depict the behaviour of the condition number by fixing the stabilization to be S K 1 presented in (26) and we consider the three polynomial bases introduced in Section 2.3. For all the This behavior is extremely interesting since it is well-known, see for example, [43], that the growth in terms of p of the condition number in triangular Spectral Elements with nodal bases is of the following sort: We want now to understand how much the ill-conditioning pollutes the convergence of the error: where we recall that u is the solution to (2) and u p is its VEM approximation. For the purpose, we consider a test case with analytic solution: u(x, y) = sin(πx) sin(πy), for which we know that the method converges exponentially, see (25). In Figure 3, we compare the errors (35) using the three meshes in Figure 1 (always using S K 1 (26) as a stabilization) and comparing the three bases q i α p−2 |α|=0 , i = 1, 2, 3. We observe that, due to the ill-conditioning of basis q 1  To understand better "how much the (linear) solver fails" when solving the system arising from (21), we consider as an exact solution: which, owing to polynomial consistency assumption (16), should be approximated exactly by the VEM.
In exact-arithmetic one would expect the error (35) to vanish, while in floating-point arithmetic the error is not zero but grows along with the condition number of the stiffness matrix.
In Figure 4, we compare error (35) using the three meshes in Figure 1 To conclude this section, we present in Figure 5 a numerical test where we fix the polynomial basis dual to the internal moments (5) to be q 2 α p−2 |α|=0 and we consider the four different stabilizations of Section 2.2. The stabilizations, as in the case of "bad" geometrical deformation presented in Section 3.2, have almost the same impact on the condition number (stabilization S K 2 seems to perform slightly worse than the other stabilizations).

Numerical results: collapsing polygons
It is also interesting to understand which is the impact of the choice of the stabilization and of the polynomial basis dual to internal moments (5) in presence of a sequence of "bad shaped" polygons (i.e., with collapsing bulk) on the condition number of the local stiffness matrix. In this way, we also test the robustness of the method when assumption (D1) is not valid.
For the purpose, we here present a quite limited and preliminary study. More precisely, we consider {K i } i∈N , sequence of "collapsing" hexagons, as those depicted in Figure 6. In particular, the coordinates of K i , the i-th element, are: Needless to say, sequence K i does not satisfy the star-shapedeness assumption (D1). In Figure 7, we depict the behaviour of the condition number of the local stiffness matrix in terms of i, parameter used in the definition of the coordinates (38) of the pentagons K i . In particular, we compare such behaviour employing the three bases q i α p−2 |α|=0 , i = 1, 2, 3, discussed in Section 2.3 and choosing p = 3 and p = 6, respectively. The stabilization is fixed to be S K 1 defined in (26). From Figure  7, we deduce that the standard choice for the polynomial basis (31) leads to a dramatic growth of the condition number. It turns out that the safest choice, in terms of ill-conditioning, is the one associated with basis q 2 α p−2 |α|=0 , which we recall is obtained by an orthonormalization of the standard monomial basis q 1 α p−2 |α|=0 via a stable Gram-Schmidt process. Basis q 3 α p−2 |α|=0 , although behaves much better than the monomial basis, is not as good as q 2 α p−2 |α|=0 . Next, in Figure 8, we compare the condition number of the stiffness matrix by fixing p = 6 and the polynomial basis q 2 α p−2 |α|=0 , which, from the previous tests, seems to be the best for the conditioning  of VEM, and by modifying the choice of the stabilizations; more precisely, we will consider the four stabilization discussed in Section 2.2. We deduce from Figure 8 that the choice of the stabilization does not have evident effects on the condition number, at least for the presented tests.
As a byproduct, in Figure 9, we consider a comparison between the four stabilizations by employing the standard monomial basis q 1 α p−2 |α|=0 as dual basis for internal moments (5). Again, we assume p = 6.

Numerical results: hanging nodes
As a final set of numerical results, we study the behavior of the condition number of the stiffness matrix using various bases and stabilizations in presence of hanging nodes collapsing on a vertex, checking thus the robustness of the method when assumption (D2) is not fulfilled.  Again, we present here only a quite limited and preliminary study. In particular, we present a sequence of "squared pentagons," that is a sequence of squares with a hanging node on a prescribed edge. More precisely, see Figure 10, we consider a sequence {K i } i∈N such that each K i , i ∈ N, has the following set of coordinates: In Figure 11, we depict the behavior of the condition number of the local stiffness matrix in terms of i, parameter used in the definition of the coordinates (39) of the pentagons K i . In particular, we compare such behavior using the three bases q i α p−2 |α|=0 , i = 1, 2, 3, discussed in Section 2.3 and choosing p = 3 and p = 6, respectively. The stabilization is fixed to be S K 1 defined in (26). The condition number is almost independent of parameter i for all choices of the canonical basis. This is not surprising since the bulk of the elements in the sequence remains the same for all i. However, when employing basis   We deduce from Figures 12 and 13 again that the behaviour of the method using different stabilizations is basically the same.

CONCLUSIONS
In this work, we addressed and suggested possible cures to the problem of the ill-conditioning of the virtual element method, arising from high values of the polynomial degree p and in presence of highly anisotropic elements. In particular, we focused our attention on the effects of the stabilization of the method and the choice of internal degrees of freedom. It turned out that, whereas various stabilizations We suggested two practical modifications of such internal degrees of freedom which greatly improve the behaviour of high-order VEM and VEM in presence of badly-shaped polygons.
It is worth to mention that we focused our attention to a simple 2D Poisson problem only. If one turns for instance to 3D problems, then the choice of stabilization (14) plays a major role, as shown in [33,15].

ACKNOWLEDGMENTS
The author acknowledges the support of the Austrian Science Fund (FWF) project F 65. We discuss here the implementation details for the construction of the method employing basis q 2 α p−2 |α|=0 as a dual basis for the internal moments (5). Henceforth, we adopt the notation of Section 2.3. Moreover, given a matrix A ∈ R n×m , we will denote by: the submatrix of A from row i to j and from column to k.
Let us define by q 2 α p |α|=0 the basis of P p (K) obtained by an L 2 ortonormalization of basis q 1 α p−2 |α|=0 using the stable Gram-Schimdt process presented in [42]. In particular, we can write (always using with a little abuse of notation the bijection [30]): where GS is the lower triangular matrix containing the orthonormalizing coefficients.
In [31], the implementation details employing basis q 1 α p−2 |α|=0 defined in (31) were discussed. In particular, it was proven that the local stiffness matrix was defined through some auxiliary matrices which we recall here: where P 0 is defined in (10). The local stiffness matrix reads: where S K denotes the matrix associated with any of the bilinear forms S K introduced in (14), where Π ∇ * denotes the matrix associated with operator Π ∇ p introduced in (9) acting from the local VE space V p (K) to P p (K) with respect to basis q 1 α p−2 |α|=0 and where Π ∇ denotes the matrix associated with the operator Π ∇ p introduced in (9) acting from the local VE space V p (K) to P p (K) with respect to basis It was shown in [31] that: The aim of this section is to write the local stiffness matrix using the new canonical basis associated with polynomial basis q 2 α p−2 |α|=0 and expanding all the projectors with respect to the polynomial spaces on polynomial basis q 2 α p |α|=0 . For the sake of simplicity, we denote the counterpart of the VEM matrices in (40) and (41)  We explain how to compute in the new setting such new matrices. We start with matrixG, which is defined as:G α,β = (∇q 2 α , ∇q 2 β ) 0,K ∀α, β ∈ N 2 , |α|, |β| = 0, . . . , p.
One simply has to compute:G = GS ·G · GS T . Now, we consider matrix G defined as: where P 0 (·) is defined in (10). We obviously have to take care only of the first line of the matrix since the remainder is inherited fromG. We distinguish two cases.
where ξ j is a proper node on the boundary.
We first, deal with the first line and we consider the two cases p = 1 and p ≥ 2.
• If ϕ 2 j is a boundary basis function, that is, j = 1, . . . , pN K , where we recall that N K is the number of edges (and vertices) of K, then: where we used that it holds ϕ 1 j | ∂K = ϕ 2 j | ∂K . • Assume now ϕ 2 j is an internal basis function. This case is a bit more involved. We write: We expand Δq 2 α into a combination of elements of the basis q 2 β n p−2 β=0 , since the Laplace operator eliminates the high (p − 1 and p) polynomial degree contributions. We get: We only need to compute the entries of matrix F. For the purpose, we test (43) with q 2 γ thus obtaining: F α,γ = (Δq 2 α , q 2 γ ) 0,K = −(∇q 2 α , ∇q 2 γ ) 0,K + ∂ n q 2 α , q 2 γ 0,∂K .
The first term is nothing but −G α,γ . We wonder how to compute the second term. We note that matrix L defined as: can be computed exactly. For the sake of completeness, we explicitly write how. Given E(K) the set of edges of K: where ω e k and ν e k , k = 0, . . . , p, are the k-th weights and nodes of the Gauß-Lobatto quadrature over edge e. It is easy to check that if we set: L α,β = ∂K (∂ n q 2 α )q 2 β , then: As a consequence: Now, we plug (43) in (42)  where C is a matrix defined as follows: C α,j = (q 2 α , ϕ 2 j ) 0,K ∀α ∈ N 2 , |α| = 0, . . . , p − 2, ∀j = 1, . . . N K dof .