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

TL;DR

This paper analyzes the ill-conditioning of the stiffness matrix in the 2D Virtual Element Method for Poisson problems, showing how basis modifications and stabilization choices affect numerical stability.

Contribution

It demonstrates how to improve the condition number by modifying internal moments and compares stabilization strategies for better numerical stability.

Findings

Ill-conditioning worsens with high-order methods and irregular polygons.

Modifying internal moments can improve the condition number.

Standard stabilization choices yield similar condition numbers.

Abstract

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.

Equations

$$|\cdot {|}_{\ell ,\omega },\Vert \cdot {\Vert }_{\ell ,\omega },(\cdot ,\cdot {)}_{0,\omega }.$$

$$\{\begin{array}{l}\text{find }u\in V\text{ such that}\\ a(u,v)=(f,v{)}_{0,\Omega }\forall v\in V\end{array},\text{where }V:=H_0^1(\Omega ),\text{and}a(\cdot ,\cdot ):=(\nabla \cdot ,\nabla \cdot {)}_{0,\Omega }.$$

$${\mathbb{B}}_p(\partial K):=\left\{v_p\in {\mathcal{C}}^0(\partial K)\mid v_p{|}_e\in {\mathbb{P}}_p(e)\forall e\text{ edge of }K\right\},$$

$$V_p(K):=\left\{v_p\in H^1(K)\mid v_p{|}_{\partial K}\in {\mathbb{B}}_p(\partial K),\Delta v_p\in {\mathbb{P}}_{p-2}(K)\right\}.$$

$$\frac{1}{|K|}{\int }_K{v}_p{q}_{\alpha },$$

$$\{{\phi }_i{\}}_{i=1}^{N_{\text{dof}}^K}$$

$${\text{dof}}_i({\phi }_j)={\delta }_{i,j}=\{\begin{array}{ll}1& \text{if }i=j\\ 0& \text{otherwise}\end{array}\forall i,j=1,\dots ,N_{\text{dof}}^K.$$

$$(q_{p-2},{\Pi }_{p-2}^0{v}_p-v_p{)}_{0,K}=0\forall q_{p-2}\in {\mathbb{P}}_{p-2}(K),\forall v_p\in V_p(K),$$

$$\{\begin{array}{l}{a}^K(q_p,{\Pi }_p^{\nabla }{v}_p-v_p)=0\\ P_0({\Pi }_p^{\nabla }{v}_p-v_p)=0\end{array}\forall q_p\in {\mathbb{P}}_p(K),\forall v_p\in V_p(K),$$

$$P_0(v_p)=\{\begin{array}{ll}\frac{1}{N_K}{\sum }_{j=1}^{N_K}{v}_p({\nu }_j)& \text{if }p=1\\ {\int }_K{v}_p& \text{otherwise}\end{array}\forall v_p\in V_p(K),$$

$$a^K(q_p,v_p)=-{\int }_K\Delta q_p{v}_p+{\int }_{\partial K}{\partial }_{\mathbf{n}}{q}_p{v}_p\forall q_p\in {\mathbb{P}}_p(K),\forall v_p\in V_p(K)$$

$$a^K(u_p,v_p)=a^K({\Pi }_p^{\nabla }{u}_p,{\Pi }_p^{\nabla }{v}_p)+a^K((I-{\Pi }_p^{\nabla })u_p,(I-{\Pi }_p^{\nabla })v_p)\forall u_p,v_p\in V_p(K).$$

$$\begin{array}{rl}{V}_p(K)& =\left\{v_p\in V_p(K)\mid v_p\in {\mathbb{P}}_p(K)\right\}\oplus \left\{v_p\in V_p(K)\mid {\Pi }_p^{\nabla }{v}_p=0\right\}\\ & ={\mathbb{P}}_p(K)\oplus \ker ({\Pi }_p^{\nabla })=:V_{p,1}(K)\oplus V_{p,2}(K),\end{array}$$

$$\begin{array}{rl}& \dim (V_{p,1}(K))=\dim ({\mathbb{P}}_p(K))=\frac{(p+1)(p+2)}{2},\\ \text{vskip}6.0ptplus2.0ptminus2.0pt& \begin{array}{rl}\dim (V_{p,2}(K))& =\dim (V_p(K))-\dim (V_{p,1}(K))\\ & =N_K\cdot p+\frac{(p-1)p}{2}-\frac{(p+1)(p+2)}{2}=(N_K-2)p-1.\end{array}\end{array}$$

$$a^K((I-{\Pi }_p^{\nabla })u_p,(I-{\Pi }_p^{\nabla })v_p)⟹S^K((I-{\Pi }_p^{\nabla })u_p,(I-{\Pi }_p^{\nabla })v_p),$$

$$c_{\ast }(p)|v_p{|}_{1,K}^2\le S^K(v_p,v_p)\le c^{\ast }(p)|v_p{|}_{1,K}^2\forall v_p\in \ker ({\Pi }_p^{\nabla }),$$

$$a_p^K(u_p,v_p)=a^K(u_p,v_p)+S^K((I-{\Pi }_p^{\nabla })u_p,(I-{\Pi }_p^{\nabla })v_p)\forall u_p,v_p\in V_p(K),$$

$$a^K(q_p,v_p)=a_p^K(q_p,v_p);$$

$${\alpha }_{\ast }(p)|v_p{|}^2\le a_p^K(v_p,v_p)\le {\alpha }^{\ast }(p)|v_p{|}^2,$$

$$⟨f_n,v_p{⟩}_K:=\{\begin{array}{ll}{\int }_K\left(f{\int }_{\partial K}{v}_p\right)& \text{if }p=1\\ {\int }_{K}f{\Pi }_{p-2}^0{v}_p& \text{otherwise}\end{array}\forall v_p\in V_p(K).$$

$$V_p:=\left\{v_p\in V\cap {\mathcal{C}}^0(\overline{\Omega })\mid v_p{|}_K\in V_p(K)\forall K\in {\mathcal{T}}_n\right\}$$

$$a_p(u_p,v_p)=\sum _{K\in {\mathcal{T}}_n}{a}_p^K(u_p,v_p),⟨f_n,v_p⟩=\sum _{K\in {\mathcal{T}}_n}⟨f_n,v_p{⟩}_K\forall u_p,v_p\in V_p.$$

$$\{\begin{array}{l}\text{find }{u}_p\in V_p\text{ such that }\\ a_p(u_p,v_p)=⟨f_n,v_p⟩\forall v_p\in V_p\end{array}.$$

$$|\cdot {|}_{1,{\mathcal{T}}_n;\Omega }^2=\sum _{K\in {\mathcal{T}}_n}|\cdot {|}_{1,K}^2.$$

$$|u-u_p{|}_{1,\Omega }\lesssim \frac{{\alpha }^{\ast }(p)}{{\alpha }_{\ast }(p)}\left\{|u-u_{\pi }{|}_{1,{\mathcal{T}}_n;\Omega }+|u-u_I{|}_{1,\Omega }+{\mathcal{F}}_n\right\}\forall u_{\pi }\in S^{p,-1}(\Omega ,{\mathcal{T}}_n),\forall u_I\in V_p,$$

$$|(f,v_p{)}_{0,\Omega }-⟨f_n,v_p⟩|\le {\mathcal{F}}_n|v_p{|}_{1,\Omega }\forall v_p\in V_p.$$

$$|u-u_p{|}_{1,\Omega }\lesssim \frac{{\alpha }^{\ast }(p)}{{\alpha }_{\ast }(p)}\frac{h^{\min (p,k)}}{p^k}|u{|}_{k+1,\Omega }$$

$$|u-u_p{|}_{1,\Omega }\lesssim \frac{{\alpha }^{\ast }(p)}{{\alpha }_{\ast }(p)}\exp (-bp),$$

$$\frac{{\alpha }^{\ast }(p)}{{\alpha }_{\ast }(p)}\lesssim p^5,$$

$$S_1^K(u_p,v_p)=\sum _{i=1}^{N_{\text{dof}}^K}{\text{dof}}_i(u_p){\text{dof}}_i(v_p),$$

Full text

1 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 tipically the solution of local PDEs and therefore are not known explicitely 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.

In a few years, thanks to high theoretical and practical flexibility of VEM, the size of associated literature has rapidly blown up. Among the other references, we recall the following: the $p$ and $hp$ version of the method [11, 12, 13, 14, 15], parabolic problems [16], Cahn-Hilliard, Stokes, Navier-Stokes and Helmoltz equations [17, 18, 19, 20, 21], linear and nonlinear elasticity problems [22, 23, 24], general elliptic problems [25], PDEs on surfaces [26], Domain Decomposition [27], application to discrete fracture networks [28], serendipity VEM [29], VEM on surfaces [26]. The implementation of the method is described in [30], whereas the basic principles of the 3D version of the method are the topic of [31, 32].

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 “bad-shaped” (i.e. nonuniformly star-shaped, with small edges…) polygons, as discussed [33].

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 employ 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 [32]. 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 [34, 35, 36] 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 bad-shaped polygons) VEM.

In particular, we show that, while the choice of the stabilization has not a deep impact on the ill-conditioning (at least for the 2D case, which is the focus of this paper), 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, i.e. low-to-moderate order VEM and VEM applied to shape-regular decompositions, the standard VEM (e.g. the one described in [10]) is preferable, since 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 paper 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 behaviour of the method when changing various stabilizations and VE bases.

In the remainder of the paper, we adopt the standard notation for Sobolev spaces, see e.g. [37, 38]. Given $\omega\subset\mathbb{R}^{2}$, we denote by $H^{\ell}(\omega)$, $\ell\in\mathbb{N}$, the Sobolev space of order $\ell$ over $\omega$; in the case $\ell=0$, we set $H^{0}(\omega)=L^{2}(\omega)$, where $L^{2}(\omega)$ is the Lebesgue space of square integrable functions over $\omega$. By $H^{1}_{0}(\Omega)$, we mean the Sobolev space $H^{1}$ of functions with zero trace. The Sobolev (semi)norms and inner products read, respectively:

[TABLE]

Further, given a Lipschitz domain $\omega$ with boundary $\Gamma$, we set $\mathbf{n}$ the normal versor associated with $\Gamma$ and $\partial_{\mathbf{n}}v=\nabla v\cdot\mathbf{n}$ the normal derivative of a sufficiently regular function $v$. By $\mathbb{P}_{\ell}(\omega)$, $\ell\in\mathbb{N}$, we denote the space of polynomials of degree $\ell$ over $\omega$. Finally, given two positive quantities $a$ and $b$, possibly depending on the discretization parameters $h$ and $p$, we write $a\lesssim b$ if there exists a positive and parameter-independent constant $c$ such that $a\leq c\,b$. Moreover, we write $a\approx b$ meaning that $a\lesssim b$ and $b\lesssim a$ simultaneously.

2 The Virtual Element Method: definition, choice of the stabilization and of the basis

Given $\Omega\subset\mathbb{R}^{2}$ a polygonal domain and $f\in L^{2}(\Omega)$, we consider the following 2D Poisson problem:

[TABLE]

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.

2.1 A family of VEM

Given $\Omega$ the computational domain of problem (2), we consider a family $\{\mathcal{T}_{n}\}_{n\in\mathbb{N}}$ of conforming polygonal decomposition of $\Omega$, where by conforming we mean that, given an edge in the skeleton of the decomposition which does not lie on $\partial\Omega$, then it is an edge of exactly two polygons. We also fix $p\in\mathbb{N}$, which will represent the degree of accuracy of the method. We associate to each $K\in\mathcal{T}_{n}$ its diameter $h_{K}$ and to decomposition $\mathcal{T}_{n}$ its mesh size function $h=\max_{K\in\mathcal{T}_{n}}h_{K}$.

Standard regularity assumptions on decomposition $\mathcal{T}_{n}$ that are usually required in VEM literature read:

• (D1)

For every $K\in\mathcal{T}_{n}$, $K$ is star-shaped (see [39]) with respect to a ball of radius greater or equal than $\gamma h_{K}$, $\gamma$ being a positive constant independent of the family of decompositions.

• (D2)

For every $K\in\mathcal{T}_{n}$ and for every $e$ edge of $K$, the length of $e$ is greater or equal than $\gamma h_{K}$, $\gamma$ being the same constant introduced in assumption (D1).

We point out that such assumptions can be relaxed, see [40].

We now define the local VE space on polygon $K\in\mathcal{T}_{n}$ following the standard definition given e.g. in [10]. Having set the space of piecewise continuous polynomials of degree $p$ over $\partial K$:

[TABLE]

we introduce the local VE space as:

[TABLE]

Space $V_{p}(K)$, which contains $\mathbb{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 $\mathcal{T}_{n}$.

We associate to space $V_{p}(K)$ the following set of linear functionals. Given $v_{p}\in 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:

[TABLE]

where $\{q_{\alpha}\}_{\alpha=1}^{\dim(\mathbb{P}_{p-2}(K))}$ is any basis of $\mathbb{P}_{p-2}(K)$.

It was proven in [10] that this is a set of unisolvent degrees of freedom. Let $N_{\text{dof}}^{K}$ be the number of such degrees of freedom, i.e. the dimension of space $V_{p}(K)$. Henceforth, we denote by $\text{dof}_{i}$ the $i$-th dof of space $V_{p}(K)$ and by:

[TABLE]

the local canonical basis of space $V_{p}(K)$, which is dual to the set of degrees of freedom $\{\text{dof}_{i}\}_{i=1}^{\dim(V_{p}(K))}$, i.e. the set of functions such that:

[TABLE]

Functions in the canonical basis associated with internal moments (5) are bubbles on element $K$ since 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_{\alpha}\}_{\alpha=1}^{\dim(\mathbb{P}_{p-2}(K))}$, 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_{\boldsymbol{\alpha}}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$.

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, 30] 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 $\Pi^{0}_{p-2}:V_{p}(K)\rightarrow\mathbb{P}_{p-2}(K)$ defined, for all $K\in\mathcal{T}_{n}$, by:

[TABLE]

which is clearly computable via internal moments (5).

The second one is a $H^{1}$ projector $\Pi^{\nabla}_{p}:V_{p}(K)\rightarrow\mathbb{P}_{p}(K)$ defined by:

[TABLE]

where $a^{K}(\cdot,\cdot)=(\nabla\cdot,\nabla\cdot)_{0,K}$ and where $P$ is an operator which fixes the constant part of the energy projector $\Pi^{\nabla}_{p}$ defined as:

[TABLE]

where $\{\nu_{j}\}_{j=1}^{N_{K}}$ is the set of vertices of polygon $K$. Operator $\Pi^{\nabla}_{p}$ is computable by means of the degrees of freedom, noting that:

[TABLE]

and recalling that $v_{p}\in 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 can not use their continuous counterparts since 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:

[TABLE]

If we split the local space $V_{p}(K)$ into a polynomial and a “pure virtual” part:

[TABLE]

then the first term in 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:

[TABLE]

This entails that the actual “pure virtual” part of space $V_{p}(K)$, i.e. $\ker(\Pi^{\nabla}_{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))\approx p$ whereas $\dim(V_{p,1}(K))\approx p^{2}$.

We emphasize that the first term in the right-hand side of (11) is explicitly computable but the second one is not. For this reason, one substitutes:

[TABLE]

where $S^{K}(\cdot,\cdot)$ is any bilinear form mimicking the continuous one, i.e. $a^{K}(\cdot,\cdot)$ on $\ker(\Pi^{\nabla}_{p})\times\ker(\Pi^{\nabla}_{p})$. In order to have a well-posed method, we demand the following stability assumption on the bilinear form $S^{K}$:

[TABLE]

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:

[TABLE]

one can prove the following properties for $a_{p}^{K}(\cdot,\cdot)$:

• (P1)

$p$-consistency: for every $q_{p}\in\mathbb{P}_{p}(K)$ and for every $v_{p}\in V_{p}(K)$, the local bilinear form $a_{p}^{K}$ satisfies:

[TABLE]

• (P2)

stability: for every $v_{p}\in V_{p}(K)$, the local bilinear form $a_{p}^{K}$ satisfies:

[TABLE]

where $\alpha_{*}(p)=\min(1,c_{*}(p))$ and $\alpha^{*}(p)=\max(1,c^{*}(p))$.

Regarding the local discrete right-hand side, we set:

[TABLE]

At this point, we are able to define the global VE space, which is obtained by a standard dof coupling of the local spaces:

[TABLE]

and the discrete global bilinear form and right-hand side:

[TABLE]

We now define a family of VEMs depending on the choice of the local stabilization $S^{K}$ defined in (14):

[TABLE]

Property (P1) guarantees that the VEM passes the patch test for piecewise-polynomial of degree $p$ solutions, i.e. if the solution of problem (2) is a polynomial of degree $p$ and the order of the method is also $p$, then the method returns as an output such a polynomial. For this a reason, $p$ is also called the degree of accuracy of the method.

On the other hand, 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 $\mathcal{T}_{n}$ as:

[TABLE]

and let us set $S^{p,-1}(\Omega,\mathcal{T}_{n})$ the space of piecewise discontinuous polynomials of degree $p$ over decomposition $\mathcal{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 of (2) and (21), respectively, the following holds true:

[TABLE]

where $\mathcal{F}_{n}$ is the smallest positive constant such that:

[TABLE]

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}(\Omega,\mathcal{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 of (2), in $H^{k+1}(\Omega)$, the following $hp$ a priori estimate holds true:

[TABLE]

and that for $u$ analytic on a proper extension of $\Omega$, the following pure $p$ approximation result holds true:

[TABLE]

where $b$ is a positive constant independent of the discretization parameters.

The dependence on $p$ of the pollution factor $\frac{\alpha^{*}(p)}{\alpha_{*}(p)}$ is clearly an issue of extreme importance in the $p$ analysis of the method. The following crude upper bound:

[TABLE]

employing a particular stabilization, was proved in [12, Theorem 4.1]. 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 $\frac{\alpha^{*}(p)}{\alpha_{*}(p)}\lesssim p^{\beta}$, for some $\beta\in(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, employing the $hp$ version of the method; see [12].

2.2 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_{\text{dof}}^{K}$ is the number of dofs of local space $V_{p}(K)$.

The first stabilization that we present is:

[TABLE]

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 [32, 31]. 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:

[TABLE]

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 [32]. If we denote by $\mathbf{K}_{C}^{K}$ the consistency part of the local stiffness matrix, i.e. the matrix counterpart of the first term in the right-hand side of (15), then we can define a stabilization $S^{K}_{3}$ through its matrix representation $\mathbf{S^{K}_{3}}$ associated with the second term in the right-hand side of (15) as follows. $\mathbf{S^{K}_{3}}$ is a diagonal matrix, whose $i$-th entry is given by:

[TABLE]

For the sake of clarity, the complete local stiffness matrix $\mathbf{K}_{n}^{K}$ reads:

[TABLE]

where $\mathbf{\Pi^{\nabla}}$ is the matrix representing the action of operator $\Pi^{\nabla}_{p}$ and $\mathbf{S^{K}_{3}}$ is defined in (28).

In [32], 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_{\text{dof}}^{\text{bndr}}$ be the number of boundary degrees of freedom of $V_{p}(K)$. Then we set:

[TABLE]

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 upon the choice of the polynomial basis dual to internal moments (5).

2.3 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_{\alpha}^{i}\}_{\alpha=1}^{\dim(\mathbb{P}_{p-2}(K))}$, $i=1,2,3$, of $\mathbb{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 employ, with an abuse of notation, the natural bijection $\mathbb{N}^{2}\leftrightarrow\mathbb{N}_{0}$ given by:

[TABLE]

We will also write occasionally:

[TABLE]

The first choice of the polynomial basis is the “standard” one, i.e. the one which is used in the majority of VEM literature, since, at the implementation level, is the most convenient. Given $\mathbf{x_{K}}=(x_{K},y_{K})$ the barycenter of $K$, we set:

[TABLE]

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 bad-shaped polygons, see [33] and the references therein.

For this reason, we suggest two possible modifications which rely on orthogonalizing processes of $\{q_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ with respect to the $L^{2}$ norm on polygon $K$.

The first modification, which allows to construct an $L^{2}(K)$ orthonormal basis $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$, is based on the stable Gram-Schmidt orthonormalization process presented in [41]. We point out that basis $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ 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_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ is well-known (see [30]), no details were given in [13] for the implementation of the method using basis $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$. Therefore, we address this issue in Appendix A.

The third choice of the polynomial basis is slightly inspired by an orthonormalizing procedure used in [33]. In order to present the third basis $\{q_{\boldsymbol{\alpha}}^{3}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$, we set matrix:

[TABLE]

Given $n_{\ell}=\dim(\mathbb{P}_{\ell}(K))$, we fix $q_{\boldsymbol{\alpha}}^{3}=q_{\boldsymbol{\alpha}}^{1}$ if $\alpha=1$ and we decompose matrix $\mathbf{H}$ into blocks as:

[TABLE]

we diagonalize matrix $\mathbf{H}_{2,2}$ and we get:

[TABLE]

This entails that matrix $\mathbf{V}\cdot\mathbf{D}^{-\frac{1}{2}}$ contains the coefficients which orthonormalize $\{q_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=1}^{p-2}$, the monomial basis of $\mathbb{P}_{p-2}(K)/\mathbb{R}$. Therefore, one has:

[TABLE]

It is worth to stress that here the orthonormalizing process is performed with a different target with respect to what done in [33]; in fact, there, the method is built employing the canonical bases computed taking moments with respect to scaled monomials $\{q_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$; however, the projectors onto polynomial spaces, i.e. $\Pi^{0}_{p-2}$ and $\Pi^{\nabla}_{p}$ defined in (8) and (9), respectively, are computed expanding with respect to polynomial basis $\{q_{\boldsymbol{\alpha}}^{3}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ defined in (32). Here, we define in addition internal moments with respect to basis $\{q_{\boldsymbol{\alpha}}^{3}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$.

The implementation issues regarding this new basis are not explicated in this paper, but are similar to those of basis $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$, which are shown in Appendix A.

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 “bad-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_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ 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 $\mathbb{P}_{p-2}(K)$; it is well-known in the theory of Spectral Elements, see e.g. [42, 43], that employing $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, employing $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_{\boldsymbol{\alpha}}^{i}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ the sets of dofs $\{\text{dof}^{i}_{j}\}_{j=1}^{N_{\text{dof}}^{K}}$ and the canonical bases $\{\varphi_{j}^{i}\}_{j=1}^{N_{\text{dof}}^{K}}$, for all $i=1,2,3$. This notation will be instrumental in Appendix A.

2.4 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 $\frac{\alpha^{*}(p)}{\alpha_{*}(p)}$, see (25) and especially (24). Secondly, since 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. Firstly, 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. Secondly, 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, since 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 $\frac{\alpha^{*}(p)}{\alpha_{*}(p)}$.

3 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 behaviour 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 behaviour 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.

3.1 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 three meshes, the basis $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ shows the best performances, whereas the standard monomial basis $\{q_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ shows the worst results.

From Figure 2, it is also clear that basis $\{q_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ entails an exponential growth of the condition number in terms of $p$.

Furthermore, employing $\{q_{\boldsymbol{\alpha}}^{i}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$, $i=2,3$, suggests instead an algebraic growth of the condition number in terms of $p$. A polynomial fitting yields:

[TABLE]

This behaviour is extremely interesting since it is well-known, see e.g. [42], that the growth in terms of $p$ of the condition number in triangular Spectral Elements with nodal bases is of the following sort:

[TABLE]

We want now to understand how much the ill-conditioning pollutes the convergence of the error:

[TABLE]

where we recall that $u$ is the solution of (2) and $u_{p}$ is its VEM approximation.

For the purpose, we consider a test case with analytic solution:

[TABLE]

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_{\boldsymbol{\alpha}}^{i}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$, $i=1,2,3$.

We observe that, due to the ill-conditioning of basis $\{q_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ for high values of $p$, the linear solver of the system (namely the one associated with the $\backslash$ command of MATLAB) does not work properly. For this reason, we highly recommend to use basis $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ in lieu of basis $\{q_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ when approximating with high-order VEM.

In order to understand better “how much the (linear) solver fails” when solving the system arising from (21), we consider as an exact solution:

[TABLE]

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, using $S^{K}_{1}$ (26) as a stabilization and the three bases $\{q_{\boldsymbol{\alpha}}^{i}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$, $i=1,2,3$.

The behaviour of basis $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ is again superior to the other two bases. More precisely, employing basis $\{q_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ has a large effect on the error for high degrees of accuracy $p$.

In order 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_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ 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).

3.2 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\in\mathbb{N}}$, sequence of “collapsing” hexagons, as those depicted in Figure 6.

In particular, the coordinates of $K_{i}$, the $i$-th element, are:

[TABLE]

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_{\boldsymbol{\alpha}}^{i}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$, $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_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$, which we recall is obtained by an orthonormalization of the standard monomial basis $\{q_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ via a stable Gram-Schmidt process. Basis $\{q_{\boldsymbol{\alpha}}^{3}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ , although behaves much better than the monomial basis, is not as good as $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$.

Next, in Figure 8, we compare the condition number of the stiffness matrix by fixing $p=6$ and the polynomial basis $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$, 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_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ as dual basis for internal moments (5). Again, we assume $p=6$.

3.3 Numerical results: hanging nodes

As a final set of numerical results, we study the behaviour of the condition number of the stiffness matrix employing 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\in\mathbb{N}}$ such that each $K_{i}$, $i\in\mathbb{N}$, has the following set of coordinates:

[TABLE]

In Figure 11, 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 (39) of the pentagons $K_{i}$. In particular, we compare such behaviour employing the three bases $\{q_{\boldsymbol{\alpha}}^{i}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$, $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 $\{q_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$, the condition number is higher.

Mimicking what done in Section 3.2, we compare in Figures 12 and 13 the condition number of the stiffness matrix by fixing $p=6$ and the polynomial bases $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ and $\{q_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$, respectively, and by considering the four stabilization discussed in Section 2.2.

We deduce from Figures 12 and 13 again that the behaviour of the method employing different stabilizations is basically the same.

4 Conclusions

In the present 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 presented in literature have almost the same effect on the condition number of the stiffness matrix, the choice of the internal degrees of freedom has a deep impact.

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 bad-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 [32, 15].

Appendix A Appendix: a hitchhikers guide for the implementation of VEM using basis $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$

We discuss here the implementation details for the construction of the method employing basis $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ as a dual basis for the internal moments (5). Henceforth, we adopt the notation of Section 2.3. Moreover, given a matrix $\mathbf{A}\in\mathbb{R}^{n\times m}$, we will denote by:

[TABLE]

the submatrix of $\mathbf{A}$ from row $i$ to $j$ and from column $\ell$ to $k$.

Let us define by $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p}$ the basis of $\mathbb{P}_{p}(K)$ obtained by an $L^{2}$ ortonormalization of basis $\{q_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ using the stable Gram-Schimdt process presented in [41]. In particular, we can write (always using with a little abuse of notation the bijection (30)):

[TABLE]

where $\mathbf{GS}$ is the lower triangular matrix containing the orthonormalizing coefficients.

In [30], the implementation details employing basis $\{q_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ 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:

[TABLE]

where $P_{0}$ is defined in (10).

The local stiffness matrix reads:

[TABLE]

where $\mathbf{S^{K}}$ denotes the matrix associated with any of the bilinear forms $S^{K}$ introduced in (14), where $\mathbf{\Pi^{\nabla}_{*}}$ denotes the matrix associated with operator $\Pi^{\nabla}_{p}$ introduced in (9) acting from the local VE space $V_{p}(K)$ to $\mathbb{P}_{p}(K)$ with respect to basis $\{q_{\boldsymbol{\alpha}}^{1}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ and where $\mathbf{\Pi^{\nabla}}$ denotes the matrix associated with the operator $\Pi^{\nabla}_{p}$ introduced in (9) acting from the local VE space $V_{p}(K)$ to $\mathbb{P}_{p}(K)$ with respect to basis $\{\varphi_{j}^{1}\}_{j=1}^{N_{\text{dof}}^{K}}$.

It was shown in [30] that:

[TABLE]

The aim of the present section is to write the local stiffness matrix employing the new canonical basis associated with polynomial basis $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ and expanding all the projectors with respect to the polynomial spaces on polynomial basis $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p}$.

For the sake of simplicity, we denote the counterpart of the VEM matrices in (40) and (41) associated with bases $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p}$ and $\{\varphi_{j}^{2}\}_{j=1}^{N_{\text{dof}}^{K}}$, with a bar at the top of each one of them.

We explain how to compute in the new setting such new matrices. We start with matrix $\overline{\widetilde{\mathbf{G}}}$, which is defined as:

[TABLE]

One simply has to compute:

[TABLE]

Now, we consider matrix $\overline{\mathbf{G}}$ defined as:

[TABLE]

where $P_{0}(\cdot)$ is defined in (10). We obviously have to take care only of the first line of the matrix since the remainder is inherited from $\overline{\widetilde{\mathbf{G}}}$. We distinguish two cases.

• ($p=1$)

$P_{0}(q_{\boldsymbol{\beta}}^{2})=\frac{1}{N_{K}}\sum_{\ell=1}^{N_{K}}q_{\boldsymbol{\beta}}^{2}(\nu_{\ell})$, where we recall that $\nu_{\ell}$ is the $\ell$-th vertex of $K$. We shall then write:

[TABLE]

• ($p\geq 2$)

In this case, we have:

[TABLE]

since basis $\{q_{\boldsymbol{\alpha}}^{2}\}_{\alpha=1}^{n_{p}}$ is $L^{2}(K)$ orthonormal by construction.

Next, we turn our attention to matrix $\overline{\mathbf{D}}$ which is defined as:

[TABLE]

and we distinguish two situations.

• •

Let us consider firstly the boundary dofs:

[TABLE]

where $\xi_{j}$ is a proper node on the boundary.

• •

Next, we deal with the internal dofs. One simply has, if $q_{\boldsymbol{\gamma}}^{2}$ is the polynomial associated with $\text{dof}_{j}$:

[TABLE]

owing again to the $L^{2}$ orthonormality of basis $\{q_{\boldsymbol{\alpha}}^{2}\}_{|\boldsymbol{\alpha}|=0}^{p}$.

Finally, we discuss the construction of matrix $\overline{\mathbf{B}}$ which is defined as:

[TABLE]

We firstly deal with the first line and we consider the two cases $p=1$ and $p\geq 2$.

• ($p=1$)

$P_{0}(\varphi^{2}_{j})=\frac{1}{N_{K}}\sum_{\ell=1}^{N_{K}}\varphi_{j}^{2}(\nu_{\ell})=\frac{1}{N_{K}}\sum_{\ell=1}^{N_{K}}\varphi_{j}^{1}(\nu_{\ell})$, where we recall that $\{\nu_{\ell}\}_{\ell=1}^{N_{K}}$ is the set of vertices of polygon $K$. Thus $\overline{\mathbf{B}}_{1,j}=\mathbf{B}_{1,j}$ $\forall j=1,\dots,N_{\text{dof}}^{K}$.

• ($p\geq 2$)

In this case, we can write:

[TABLE]

since $q_{1}^{2}=\mathbf{GS}_{1,1}q_{1}^{1}=\mathbf{GS}_{1,1}$. Thus $\overline{\mathbf{B}}_{1,j}=\mathbf{GS}_{1,1}^{-1}\mathbf{B}_{1,j}$ $\forall j=1,\dots,N_{\text{dof}}^{K}$.

Next, we treat all the other lines. We must compute $(\nabla q_{\boldsymbol{\alpha}}^{2},\nabla\varphi_{j}^{2})_{0,K}$. Again, we consider two different situations.

• •

If $\varphi_{j}^{2}$ is a boundary basis function, i.e. $j=1,\dots,pN_{K}$, where we recall that $N_{K}$ is the number of edges (and vertices) of $K$, then:

[TABLE]

where we used that it holds $\varphi_{j}^{1}|_{\partial K}=\varphi_{j}^{2}|_{\partial K}$.

• •

Assume now $\varphi_{j}^{2}$ is an internal basis function. This case is a bit more involved. We write:

[TABLE]

We expand $\Delta q_{\boldsymbol{\alpha}}^{2}$ into a combination of elements of the basis $\{q_{\boldsymbol{\beta}}^{2}\}_{\beta=0}^{n_{p-2}}$, since the Laplace operator eliminates the high ($p-1$ and $p$) polynomial degree contributions. We get:

[TABLE]

We only need to compute the entries of matrix $\overline{\mathbf{F}}$. For the purpose, we test (43) with $q_{\boldsymbol{\gamma}}^{2}$ thus obtaining:

[TABLE]

The first term is nothing but $-\overline{\widetilde{\mathbf{G}}}_{\alpha,\gamma}$. We wonder how to compute the second term. We note that matrix $\mathbf{L}$ defined as:

[TABLE]

can be computed exactly. For the sake of completeness, we explicitly write how. Given $\mathcal{E}(K)$ the set of edges of $K$:

[TABLE]

where $\omega_{k}^{e}$ and $\nu_{k}^{e}$, $k=0,\dots,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:

[TABLE]

then:

[TABLE]

As a consequence:

[TABLE]

Now, we plug (43) in (42) obtaining:

[TABLE]

where $\overline{\mathbf{C}}$ is a matrix defined as follows:

[TABLE]

One has:

[TABLE]