p- and hp- virtual elements for the Stokes problem

We analyse the p- and hp-versions of the virtual element method (VEM) for the the Stokes problem on a polygonal domain. The key tool in the analysis is the existence of a bijection between Poisson-like and Stokes-like VE spaces for the velocities. This allows us to re-interpret the standard VEM for Stokes as a VEM, where the test and trial discrete velocities are sought in Poisson-like VE spaces. The upside of this fact is that we inherit from [7] an explicit analysis of best interpolation results in VE spaces, as well as stabilization estimates that are explicit in terms of the degree of accuracy of the method. We prove exponential convergence of the hp-VEM for Stokes problems with regular right-hand sides. We corroborate the theoretical estimates with numerical tests for both the p- and hp-versions of the method.


Introduction
The virtual element method (VEM) is an increasingly popular tool in the approximation to solutions of fluido-static and dynamic problems in polygonal/polyhedral meshes. In particular we recall: the very first paper on low-order VEM for Stokes [2]; its high-order conforming [11] and nonconforming versions [20,33]; conforming [12] and nonconforming VEM for the Navier-Stokes equation [32]; mixed VEM for the pseudo-stress-velocity formulation of the Stokes problem [17]; mixed VEM for quasi-Newtonian flows [19]; mixed VEM for the Navier-Stokes equation [24]; other variants of the VEM for the Darcy problem [18,44,45]; analysis of the Stokes complex in the VEM framework [9,13]; a stabilized VEM for the unsteady incompressible Navier-Stokes equations [29]; implementation details [23].
Notwithstanding, all the above articles refer to the h-version of the method (i.e., when the convergence is achieved by refinement of the underlying mesh while keeping the order of the approximation fixed) and the convergence analysis is performed assuming enough smoothness of the solutions to the problem under consideration. This is not the case when the domain of the equation is polygonal/polyhedral. In fact, even with smooth data, solutions are expected to have singularities at the corners of the domain; see, e.g., [26,35]. More precisely, it can be proven that they belong to Kondrat'ev spaces, i.e., weighted Sobolev spaces with weight given by a function of the distance from the corners of the domain; see definitions (4) and (5) below.
For this reason, employing hp spaces arises as a natural technique in order to construct methods, which lead to an exponential decay of the error. This approach has been investigated in a plethora of works, in the framework of conforming and nonconforming finite element methods. We recall the following works, which relate to the hp approximation of problems of Stokes and Navier-Stokes type: hp dG primal and mixed methods for the Stokes equation [40,42]; mixed discontinuous Galerkin (dG) finite element methods for the Navier-Stokes equation [34]; error indicator for the Stokes equation [14]; analysis of Stokes flows [25]; mixed hp-dG methods for incompressible flows [37][38][39] and their a posteriori version [28]; spectral elements for Stokes eigenvalue problems [43].
The main contribution of this paper is given by the development of the analysis of p-and hp-VEM for the approximation of solutions to the Stokes problem, building upon the analysis for p-and hp-VEM for the Poisson problem in [6,7]. The key tool in the analysis is the proof of the existence of a bijection between Poisson-like [5] and Stokes-like [11] VE spaces for the velocities. This allows us to re-interpret the standard VEM for Stokes [11] as a VEM, where the test and trial discrete velocities are sought in Poisson-like VE spaces. The upside of this fact is that we inherit from [7] an explicit analysis of best interpolation results in VE spaces, as well as stabilization estimates that are explicit in terms of the degree of accuracy of the method.
We prove that the hp-version of the method converges exponentially in terms of the cubic root of the number of degrees of freedom when the right-hand side of the Stokes problem in a polygonal domain is analytic. In addition, we also show that the p-version of the method converges algebraically if the solution is sufficiently regular, and exponentially in terms of the degree of accuracy when the solution is analytic.
In the remainder of this section, we introduce some notation, the continuous problem we are interested in, namely a Stokes problem in a two dimensional polygonal domain, and discuss the regularity of solutions to this kind of problems in polygonal domains. Finally, we conclude this section by presenting the structure of the paper.

Notation
We employ the standard notation for Sobolev spaces [1]. More precisely, given a domain D ⊂ R d , d = 1, 2, we denote the Sobolev space of integer order s ∈ N by H s (D). We endow H s with standard Sobolev inner products, seminorms and norms by (·, ·) s,D , | · | s,D , · s,D .
Fractional Sobolev spaces can be defined via interpolation theory. Moreover, we set P p (D) as the space of polynomials of total degree at most p over the domain D As customary, given two positive quantities a and b, we write a b meaning that there exists a positive constant c independent of the discretization parameters such that a ≤ c b. Moreover, we write a b if a b and b a at once.

The continuous problem
Let Ω ⊂ R 2 be a polygonal domain with boundary Γ and f ∈ [L 2 (Ω)] 2 . We want to approximate the solution to the following problem: find u and s such that Define the spaces and the bilinear forms The weak formulation of problem (1) reads Problem (3) is well-posed: see, e.g., [15].

Regularity of the solution
The regularity of the solution (u, s) to the Stokes problem (1) in the polygonal domain Ω depends on the shape of the domain. In particular, even if the right-hand side f is analytic, the corners of the domain give rise to corner singularities in the solution, which limit its regularity in the scale of classical Sobolev spaces. In order to properly characterize the solution to the Stokes problem, we resort to corner-weighted Sobolev spaces, of the kind firstly proposed in [30].
Assume that the polygon Ω has n c ∈ N corners, which we denote by C = {c i ∈ R 2 , i = 1 . . . , n c }. Set the amplitude of the internal angles at each corner c i ∈ C as φ ci ∈ (0, 2π)\{π} and the Euclidean norm in R 2 by | · |. Then, given the vector γ = {γ ci ∈ R, c i ∈ C} ∈ R nc and k ∈ N 0 , define the weight function For ∈ N 0 and γ ∈ R nc , introduce the seminorm and associated norm where we use the notation ∂ α = ∂ α1 x1 ∂ α2 x2 . We define the homogeneous Kondrat'ev space as Furthermore, we introduce the class of weighted analytic functions For each vertex c i ∈ C of Ω, λ ci denotes the smallest positive solution to the following equation: Observe that, for all φ ci ∈ (0, 2π) \ {π}, we have λ ci > 1/2. Furthermore, for all 0 < φ c < π, i.e., in presence of convex corners, we have λ ci = 1.
The following result is a finite regularity shift result in weighted Sobolev spaces for solutions to the Stokes problem; see [26,Theorem 5.7] and [31,Section 5]; see also [35,Proposition 1.8] for the case of homogeneous spaces. Theorem 1.1. Let ∈ N 0 and γ be such that 0 < γ ci − 1 < λ ci for all c i ∈ C. Assume that f ∈ K γ−2 (Ω) 2 and let (u, s) ∈ V × Q be the (unique) solution to (1) with right-hand side f . Then, there exists C > 0 such that Furthermore, if the right-hand side belongs to analytic weighted spaces, then also the solution to the Stokes problem belongs to the same spaces, as stated in the following result; see [26,Theorem 5.7]. Theorem 1.2. Let γ be such that 0 < γ c − 1 < λ c for all c ∈ C. Let f ∈ K γ−2 (Ω) 2 and (u, s) ∈ V × Q be the solution to (1) with right-hand side f . Then u ∈ K γ (Ω) 2 and s ∈ K γ−1 (Ω).

Structure of the paper.
In Section 2, we construct the VEM for the approximation of solutions to problem (3). Differently from the standard approach of [11], we show that the VEM for the Stokes equation can be reinterpreted as a VEM where the velocity space is Poisson-like [5]. Section 3 is concerned with the derivation of a priori estimates on velocities and pressures. Among the key points here, we prove the validity of the inf-sup condition and stabilization bounds, which are explicit in terms of the degree of accuracy of the method. The exponential convergence for the p-and hp-versions of the method are theoretically proven in Section 4, and numerically validated in Section 5. We draw some conclusions in Section 6.

Meshes and the virtual element method
In this section, we present the virtual element method for the approximation of solutions to (3). More precisely, we begin by introducing sequences of polygonal meshes partitioning the domain Ω and their properties in Section 2.1. Next, in Section 2.2, we recall the virtual element spaces introduced in [11], whereas, in Section 2.3, we construct computable bilinear forms and exhibit the method. We devote, then, Section 2.4 to recalling the standard virtual element method from [5]. Indeed, we show that the virtual element method for the Stokes equation can be re-interpreted as a method, where the velocity is sought, in Poisson-like virtual element spaces. This fact will play an important role in the analysis presented in Section 3 below.

Meshes
Here, we introduce the polygonal meshes upon which we will construct the virtual element method. Specifically, we consider sequences {T n } n∈N of meshes which partition the domain Ω into conforming, nonoverlapping polygons. Fix n ∈ N, i.e., fix one of the meshes in the sequence. We denote the set of vertices and edges in T n by V n and E n , respectively. Next, fix K ∈ T n . We denote its diameter and centroid by h K and x K , respectively. Moreover, E K represents its set of edges. We define h := max K∈Tn h K .
The set of vertices V n and edges E n can be decomposed into internal and boundary, i.e., contained in Γ = ∂Ω, ones. We write V I n , V B n , E I n , and E B n , respectively. We denote the length of each edge e ∈ E n by h e .
We state the following assumptions on the sequence of meshes: for all n ∈ N, there exists γ ∈ (0, 1) such that (A0-p) the mesh T n is quasi-uniform, i.e., for all K 1 and K 2 ∈ T n , there holds γh K1 ≤ h K2 ≤ γ −1 h K1 ; (A0-hp) the mesh T n is locally quasi-uniform, i.e., for all neighbouring K 1 and K 2 ∈ T n , there holds (A1) for all K ∈ T n , K is star-shaped with respect to a ball with radius larger than or equal to γh K ; (A2) for all K ∈ T n and for all e ∈ E K , there holds h K ≤ γh e .
The assumptions (A1) and (A2) will be used throughout the whole paper. Instead, assumptions (A0p) and (A0-hp) will be considered when dealing with the p-and hp-version of the method, respectively.
For the sake of exposition, we construct the method for uniform p only, and postpone to Section 4.2 the variable degree case. Remark 1. The forthcoming analysis can be also extended to more general geometries; see, e.g., [10,16,21]. For the sake of clarity, we stick to the setting detailed above.
We denote the space of piecewise discontinuous polynomials of degree p ∈ N over T n by P p (T n ).

The Stokes virtual element spaces
Here, we recall from [11] the virtual element spaces which we will use in the discretization of the Stokes problem (3). Henceforth, p ∈ N denotes the degree of accuracy of the method. Given K ∈ T n , set and introduce the subspace In [11], H p (K) is chosen as the L 2 (K)-orthogonal complement in [P p (K)] 2 of G p (K), denoted G ⊥ p (K). In practical computations, see [23], a convenient choice is provided by the space In what follows, we do not impose orthogonality in (8), but only require that H p (K) is such that (8) is a direct sum.
Recall that E K denotes the set of edges of the element K and introduce Define the local bilinear forms Consider the following local Stokes problem: Given q ⊕ p−2 ∈ H p (K) and q p−1 ∈ P p−1 (K)/R, Set the local Stokes-like virtual element space for the velocity as follows: V n (K) := v n ∈ [H 1 (K)] 2 | v n solves a problem of the form (9) .
We introduce the following linear functionals on V n (K): given v n ∈ V n (K), define the point values at the vertices of K; • Dv K 2 (v n ): the point values at the p − 1 Gauß-Lobatto points on each edge e ∈ E K ; • given {q ⊕ α } a basis of H p , the "complementary" moments • given {q α } p−1 α=1 a basis of P p−1 (K)/R, the "divergence" moments Lemma 2.1. The above linear functionals are a set of degrees of freedom for V n (K).

Proof. See [11, Proposition 3.2].
We define the H 1 -conforming global Stokes-like velocity space as follows: We endow this space with the set of degrees of freedom, which is obtained by a standard H 1conforming dof coupling of the local ones.
The above degrees of freedom allow us to compute two projection operators; see [11, Sections 3.2 and 3.3]. The first one is the We define the global projector These two operators are instrumental in the design of the virtual element methods; see Section 2.3 below.
For future convenience, introduce the broken Sobolev space and associate with it the broken Sobolev seminorm and norm Finally, set the pressure space as

The virtual element method
Here, we design computable discrete bilinear forms and right-hand side and introduce the virtual element method for the approximation of solutions to the Stokes problem (3).

Discrete bilinear forms.
We introduce the elementwise discrete bilinear form a K n given by where, for all K ∈ T n , S K : H 1 (K) × H 1 (K) → R is a computable local stabilizing bilinear form, which is computable from the degrees of freedom introduced in Section 2.2. We postpone the discussion about further properties of the stabilizing bilinear forms S K to Section 3.2 below. The global discrete bilinear form reads As for the discretization of the bilinear form b(·, ·) in (2), we observe that the divergence of functions in the space V n (K) is polynomial and can be expressed in closed form in terms of their degrees of freedom. Therefore, no approximation is necessary for the second bilinear form and we define b n (v n , q n ) := b(v n , q n ) ∀v n ∈ V n , ∀q n ∈ Q n .

Discrete right-hand side.
Define the global piecewise L 2 projector Π 0 p−2 as follows: The virtual element method.
The virtual element method for the Stokes problem (3) reads as follows:

An equivalent formulation in Poisson-like virtual element spaces
We recall the vector Poisson-like virtual element space, see [5], for this will allow us to reinterpret method (17) in a way that is more convenient for the sake of the analysis in Section 3 below. Given K ∈ T n , set The global H 1 standard Poisson-like virtual element space reads The operators Dv K i , i = 1, . . . , 4 introduced in Section 2.2 are unisolvent degrees of freedom for both V n (K) and V n (K), as stated in the following lemma, where we also prove that such degrees of freedom identify a bijection between the two virtual element spaces.

Lemma 2.2. For all K ∈ T n , there exists a Stokes-to-Poisson bijection
Proof. Given K ∈ T n , introduce the following auxiliary set of degrees of freedom: given v n ∈ V n (K), • Dv • given {q ⊕ α } the basis of H p−2 used in (10), the moments • given {p α } p−1 α=1 a basis of G p−2 such that p α = ∇q α , with q α defined in (11), the moments Since [P p−2 (K)] 2 = G p−2 ⊕H p−2 , this is indeed a set of degrees of freedom; see [5,Proposition 4.1].
Finally, for α = 1, . . . , p − 1, let {q α } α be the basis of P p−1 (K)/R used in (11). We require Using that dim(V n (K)) = dim( V n (K)), we get that T K StP is a bijection. As an immediate consequence, we have the following result.
The two next lemmata are instrumental in order to prove Proposition 2.6 below.

Lemma 2.4. Let T K
StP be the bijection introduced in Lemma 2.2. Then, the following identity is valid:

Lemma 2.5. Let T K
StP be the bijection introduced in Lemma 2.2. Then, we have Proof. Let v n ∈ V n (K) and denote v n = T K StP v n ∈ V n . An integration by parts yields Since v n|∂K = v n|∂K and using Lemma 2.4, we deduce The second identity in (24) is a direct consequence of Lemma 2.4.
Define the global bijection as (T StP v n ) |K = T K StP (v n|K ) for all v n ∈ V n and K ∈ T n . The following result is a direct consequence of Lemmata 2.4 and 2.5. 1 Here and in what follows (IBP) means 'integration by parts'. and Proof. The identities in (27) follow from Lemma 2.5 and the definitions of Π ∇ p and Π 0 p−2 directly. In order to show (26), remark that, due to Lemma 2.2, Then, (26) follows from summing up the contributions of each integral in K.
In words, Proposition 2.6 states that, given two functions in the virtual element spaces V n and V n sharing the same value of the degrees of freedom, their Π ∇ p and Π 0 p−2 projections, as well as their evaluations through b(·, q n ) for all q n ∈ Q n , are the same.

A priori estimates
In this section, we prove the well-posedness and provide an abstract error analysis of method (17). To this aim, we first prove that the bilinear form b(·, ·) satisfies a discrete inf-sup condition independently of the degree of accuracy of the method; see Section 3.1. Secondly, in Section 3.2, we analyse the discrete bilinear form a n (·, ·) and show that, under suitable assumptions on the stabilization terms, it is coercive and continuous. Notably, the coercivity and continuity constants are determined using Poisson-like spaces and are explicit in terms of the degree of accuracy p of the method. The abstract error analysis on the velocities and pressures is provided in Sections 3.3 and 3.4, respectively. The bounds herein proven are instrumental in deducing the rate of convergence of the error of the method, which is the topic of Section 4 below.

The discrete inf-sup condition
The discrete inf-sup stability of method (17) has been shown in [11] already. Here, we recall its proof, and show that the discrete inf-sup constant is independent of the degree of accuracy p.
We start by recalling a classical result on the inf-sup constant for star-shaped domains.
Proof. As is customary, we use Fortin's trick, i.e., we show the existence of an operator Π n : V → V n and a positive constant C independent of p such that This implies the validity of the inf-sup stability of the spaces V n and Q n ; see, e.g., [15]. We devote the remainder of the proof to showing the existence of such operator Π n and constant C.
Let W n be a low-order (p = 2) virtual element space for the velocity. By [11,Proposition 4.2], there exists v n ∈ W n such that b(v n , q n ) = b(v, q n ) for all q n ∈ P 0 (T ) v n 1,Ω ≤ C v 1,Ω .
In each element K, we introduce a bubble function w K n ∈ V n (K) such that In other words, we construct w K n such that, in each element K ∈ T n , Dv K i (w K n ) = 0, i = 1, 2, 3, By the standard well-posedness of the above Stokes problem, we claim that In order to show (30), first observe that Next, denote the inf-sup constant of the continuous Stokes problem in K with homogeneous Dirichlet boundary conditions by β(K). This gives whence (30) follows.
Next, consider w n ∈ V n defined as w K n in each element K ∈ T n and define Π n v = w n + v n . By construction, it follows that From (29) and (30), we deduce that Π n is H 1 (Ω)-stable, with stability constant independent of the degree of accuracy p.

Stabilization, coercivity, and continuity: well-posedness of the VEM
In this section, we analyse the properties of the discrete bilinear form a n (·, ·). Notably, we show that suitable choices of the stabilization forms yield to a coercive and continuous bilinear form. Furthermore, the coercivity and continuity constant are explicit in terms of the degree of accuracy of the method p. The main ingredient is given by the properties of the bijection T K StP ; see Lemma 2.2.
In order to investigate the stability of the method, we require an additional property on the stabilization bilinear forms: For all u n , v n ∈ V n (K) and u n , v n ∈ V n (K) such that Dv K i (v n ) = Dv K i ( v n ), i = 1, 2, 3, 4, i.e., v n = T StP v n with T StP defined in Lemma 2.2, Furthermore, we assume that, for all p ∈ N and K ∈ T n , there exist two positive constant α * (p) < α * (p), such that Following, e.g., [5], we can prove that α * (p) and α * (p) are the coercivity and continuity constants for the discrete bilinear form a n (·, ·). The actual dependence on p of the two constants hinges upon the definition of the stabilizing bilinear forms S K (·, ·) in (16); see Remark 2 below for an explicit choice of the stabilization together with the explicit dependence in terms of the degree of accuracy. As in [5], the properties of the discrete bilinear form a n (·, ·) entail that the method is stable and p-polynomially consistent. We have the following well-posedness result. Proof. The assertion follows from the continuity of the bilinear form a n and b n , the coercivity of a n , the discrete inf-sup condition (28), and standard argument as in [15]. Remark 2. An example of an explicit stabilization S K such that (31) and (32) are valid is as follows: All the terms on the right-hand side of (33) are computable via the degrees of freedom Dv K i , i = 1, . . . , 4 explicitly. Furthermore, (31) is valid thanks to Lemmata 2.2 and 2.5. On the other hand, the bounds in (32) can be proven as in [7,Theorem 2], with explicit stability constants for all > 0 and where ω K denotes the largest angle of K. In all fairness, the practical dependence of the stabilization constants in terms of p results to be much milder numerically; see [6, Section 4.6] and [7, Section 4.1].

Why did we assume (31)?
The reason we have introduced the auxiliary Poisson-like virtual element space V n in (18) and analysed its relation with the Stokes-like virtual element space V n in (12) is that we can exploit previous stability bounds that are explicit in terms of the degree of accuracy p; see [7,Section 4].

A priori estimate on the velocity
In this section, we prove some upper bounds, which will be instrumental in the analysis of the convergence for the error on the velocity.
Introduce the weakly divergence-free subspace of V n Z n := v n ∈ V n : b( v n , q n ) = 0 for all q n ∈ Q n .
For future use, we also introduce the weakly divergence-free subspace of V n Z n := {v n ∈ V n : b(v n , q n ) = 0 for all q n ∈ Q n } .
Moreover, let F n denote the smallest constant such that The first result is an upper bound on the error between the solution to the continuous problem and the discrete solution mapped through the bijection in (25). (3), u n ∈ V n be the virtual element solution to (17), and T StP be the bijection defined in (25). Then, the following bound is valid:

Lemma 3.4. Let u be the solution to
Proof. Introduce u n = T StP u n . Since b(u n , q n ) = 0 for all q n ∈ Q n , use (26) to get that u n ∈ Z n . Moreover, by (27) and (31), u n is the solution to the reduced problem find u n ∈ Z n such that a n ( u n , v n ) = (Π 0 p In fact, u n solves the Stokes-like counterpart find u n ∈ Z n such that a n (u n , v n ) = (Π 0 The analysis proceeds with classical tools for a priori estimates for virtual element methods; see, e.g., [5]. For any z n ∈ Z n , the triangle inequality yields |u − u n | 1,Tn ≤ |u − z n | 1,Tn + | z n − u n | 1,Tn . Denoting δ n = z n − u n ∈ Z n , we compute, for all u π ∈ [P p (T n )] 2 , where the last inequality follows from the definition of F n and from the Cauchy-Schwarz inequality. Dividing both sides by |δ n | 1,Tn gives | z n − u n | 1,Tn ≤ 1 α * (p) (F n + α * (p)| z n − u| 1,Tn + (α * (p) + 1)|u π − u| 1,Tn ) .
The next result is an upper bound on the error between the solution to the continuous problem and the H 1 projection of the discrete Stokes-like solution.
Lemma 3.5. Let u and u n ∈ V n be the solutions to (3) and (17), respectively. Then, we have Proof. Let u n = T StP u n ∈ V n . Use Proposition 2.6 to get For all u π ∈ [P p (T n )] 2 , we have Combining (35), (39), and (40), the assertion follows.
Next, we show an upper bound on the best error on the Poisson-like weakly divergence free subspace Z n in terms of a best error in terms of functions in the Poisson-like virtual element space V n .
Then, the following upper bound is valid: Proof. We begin by proving a discrete "switched inf-sup" condition. Introduce Z C n the complementary space of Z n defined in (34) in V n . In Lemma 3.2, we proved the existence of a surjective operator n : Q n → V n such that In particular, the discrete inf-sup condition (28) can be written as Thence, for all v n ∈ Z C n , thanks to the surjectivity of n , we can write ≤ sup qn∈Qn ( n q n , v n ) 1,Ω q n 0,Ω = sup qn∈Qn b(v n , q n ) q n 0,Ω .
For each v n ∈ V n , define w n ∈ Z C n as the solution of find w n ∈ Z C n such that b(w n , q n ) = b( v n , q n ) ∀q n ∈ Q n .
This problem has a unique solution due to the continuity and the discrete "switched inf-sup" stability in (44) of the bilinear form b(·, ·); see, e.g., [15]. Furthermore, the following a priori estimate is valid: Next, define where T StP is the bijection in (25). Thanks to (26), we get We deduce that z n ∈ Z n . Then, we have ≤ a n (T StP w n , T StP w n ) (27), (31) = a n (w n , w n ) (32) ≤ (1 + α * (p))|w n | 2 1,Ω .

A priori estimate on pressure
In this section, we prove some upper bounds which will be instrumental in the analysis of the convergence of the error on the pressure obtained by the VEM.
(Ω) and (u n , s n ) ∈ V n × Q n be the solutions to (3) and (17), respectively. Recall that the bijection T StP is defined in (25). Then, the following bound is valid: Proof. For all q n ∈ Q n , the triangle inequality yields s − s n 0,Ω ≤ s − q n 0,Ω + s n − q n 0,Ω .
By the discrete inf-sup condition (28), there exists v n ∈ V n such that β n s n − q n 0,Ω ≤ b(v n , s n − q n ) |v n | 1,Ω .
We have b(v n , s n − q n ) = b(v n , s − q n ) + b(v n , s n − s) and |b(v n , s − q n )| ≤ |v n | 1,Ω s − q n 0,Ω .
For u n = T StP u n , we deduce whence the assertion follows. Define Combining Lemmata 3.4, 3.5, and 3.7, we obtain the following result.

The convergence rate of the p-and hp-versions
In Section 3, we have established an abstract error analysis for method (17). Notably, we have proven that the error on the velocity and the pressure can be estimated from above in terms of best polynomial approximation and best interpolation results in virtual element spaces. With this at hand, in this section, we state the convergence of the p-and hp-versions of method (17) for analytic, weighted analytic, and finite Sobolev regularity solutions; see Sections 4.1 and 4.2, respectively.

p-VEM
Since all the necessary best approximation results have been proven in [6], we state the main convergence result only.
Furthermore, if u and s are the restrictions of suitable analytic functions over an extension of the domain 2 Ω, then there exist two positive constants C 1 and C 2 independent of the discretization parameters such that |u − Π ∇ p u n | 1,Tn + β n s − s n 0,Ω ≤ C 1 exp(−C 2 p).
Proof. Starting from the abstract error analysis in Theorem 3.8, it suffices to be able to show hand p-upper bounds on the four terms appearing on the right-hand side of (51 From Theorem 4.1, we have that the p-version of the method converges exponentially for analytic solutions and algebraically for solution with (sufficiently high) finite Sobolev regularity. However, since solutions to the Stokes problem are in general singular, as detailed in Theorem 1.1, we are also interested in analysing the convergence of the hp-version of the method. Indeed, it is known that such approach allows for exponential convergence with respect to a suitable root of the total number of degrees of freedom for singular solutions as well. We postpone the design of hp-virtual element spaces for the Stokes problem, as well as the convergence of the error, to Section 4.2 below.
Remark 4. An additional reason why the hp-version is more suited than the p-version for the approximation of singular solutions to the Stokes problem is that the algebraic rate of convergence in (52) contains the suboptimal term γ(p) due to the stabilization of the method.
which differ from those that are typically investigate in the VEM literature, i.e., |u − u n | 1,Tn + β n s − s n 0,Ω .
The reason for this is that we need to resort to Poisson-like spaces, when performing the theoretical analysis, and we know from Proposition 2.6 that functions in Poisson-like and Stokes-like virtual element spaces, sharing the same degrees of freedom, have the same Π ∇ p projection. In turn, we had to resort to Poisson-like virtual element spaces, because we are not able to construct a stabilization on Stokes-like virtual element spaces, with bounds on the stabilization constants, which are explicit in terms of the degree of accuracy of the method. On the positive side, the two errors that we bound are those that we actually compute in the numerical experiments presented in Section 5 below.

hp-VEM
In the present section, we construct hp-virtual element spaces for the approximation of nonsmooth solutions to the Stokes problem (3). The main idea of the construction hinges upon employing • geometric refinement of the mesh towards the singular points; • p-refinement in the elements where the solution is smooth.
For the sake of exposition, assume that the right-hand side f in (3) is smooth. Thanks to Theorems 1.1 and 1.2, the solution (u, s) to (3) consists of two functions that are smooth everywhere but at neighbourhoods of the vertices of the polygonal domain Ω. There, the Sobolev regularity is known a priori and depends on the amplitude of the angle.
The first step in the construction of hp-virtual element spaces resides in introducing the layer of the mesh associated with the set of vertices C. We assume that the mesh T n consists of n + 1 layers, where the first one is given by L 0 n := {K ∈ T n | there exists a unique c ∈ C such that c ∈ E K }, and the others are defined recursively as L j n := {K ∈ T n | K ∈ ∪ j−1 =0 L n ; ∃ K ∈ L j−1 n such that K ∩ K = ∅} ∀j = 1, . . . , n.
Further, for each K ∈ T n , we denote any of the closest corner of the domain to K, i.e., any of the c ∈ C such that dist(c, K) ≤ dist(c, K) for allc ∈ C \ {c}, by c K . For the sake of simplicity, we assume the uniqueness of such a vertex. With this at hand, we say that the sequence of meshes {T n } n∈N is geometrically refined towards c K if there exists a grading parameter σ ∈ (0, 1) such that, for all n ∈ N, and The conditions (55)-(56) asserts that the elements abutting the vertices in N are small, whereas the elements in the layers with large index j have fixed size asymptotically. Note that the assumption (A0-hp) is satisfied automatically. We require an additional assumption, which is necessary to show the exponential convergence result of Theorem 4.2 below; see [7,Assumption (D4)].
(A4-hp) For all n ∈ N, let T 1 n = T n \ L 0 n . There exist a collection of squares Q n such that • card(Q n ) = card(T 1 n ); for each K ∈ T 1 n , there exists Q = Q(K) ∈ Q n such that K ⊂ Q and h K h Q . Additionally, dist(c K , Q(K)) h K ; • every x ∈ Ω belongs at most to a fixed number of squares Q, uniformly in the discretization parameters.
In addition, for all K ∈ L 0 n , K is star shaped with respect to c K and the subtriangulation obtained by joining c K with the other vertices of K is shape regular.
Although necessary in the proof of Theorem 4.2, the condition (A4-hp) is not necessary in practice. For instance, the hp-version of the method converges exponentially also on meshes, as those depicted in Figure 2 (right); see Section 5.2 below.
Next, we introduce a distribution of degrees of accuracy, by picking a high degree on large elements, where the solution is smooth, and decrease such degree linearly while decreasing the size of the elements. More precisely, given a positive parameter µ, set n el := card(T n ) and introduce p ∈ N n el as follows: where K ∈ L j n ∀j = 0, . . . , n + 1.
The vector p represents the distribution of the degrees of accuracy over a mesh T n . Given n edge := card(E n ), we also introduce a vector p En ∈ N n edge , which represents the distribution of polynomial degrees over the skeleton of the mesh, and is defined as n and e ∈ E K for some K ∈ T n .
We can now define the hp-space for the velocities as the space of functions that are piecewise polynomials with distribution p En over the skeleton of the mesh and which solve problems of the form (9) with right-hand side being polynomials of degree p K −2 (vector) and p K −1, respectively, on K. On the other hand, we define the hp-virtual element space for the pressure as the space of piecewise polynomials of degree p K on K.
Using the abstract analysis in Theorem 3.8 together with the tools in [7, Section 5], we state the following result. Let the right-hand side f be analytic in Ω and let (u, s) and (u n , s n ) ∈ V n × Q n be the solutions to (3) and (17), respectively. For all n ∈ N, define N V := card(V n ) + card(Q n ). Then, there exist two positive constants C and b such that Proof. Starting from the abstract error analysis in Theorem 3.8, it suffices to be able to show hpupper bounds on the four terms appearing on the right-hand side of (51). More precisely, from [7, Lemmata 2 and 3], there exist constants C 1 and b 1 > 0 such that, for all n ∈ N, Furthermore, noting that the pressure s has the same regularity as the components of the gradient of the velocity u, with similar arguments, we deduce that there exist constants C 2 and b 2 > 0 such that, for all n ∈ N, inf qn∈Qn s − q n 0,Ω ≤ C 2 exp(−b 2 n).
Then, we deduce from [7, Lemmata 4 and 5] that there exist constants C 3 and b 3 > 0 such that, for all n ∈ N, Finally, the estimate inf vn∈ Vn for constants C 4 , b 4 > 0, independent of n is a consequence of [7, Lemmata 6 and 7]. We remark that (55), (56), and (57) imply card(T n ) n, see, e.g., [27,Equation (5.6)]. Since dim(V n (K)) p 2 K and dim(P p K (K)) p 2 K for each K ∈ T n , (57) gives N V n 3 . Since γ(p) grows at most algebraically in terms of p, we absorb the term γ(p) space V n (K), which is dual to the degrees of freedom {dof j (·)} dim(Vn(K)) j=1 introduced in Section 2.2. We define It is known [8,36] that stabilizations of this sort lead to effective performance of the method.
We highlight that we also tested the method with the stabilization (33), and this leads to results that are comparable to those that we present in the forthcoming sections.

Polynomial bases.
We refer to [23], as for the choice of the polynomial bases. We underline that this choice could be improved; see Remark 6 below.

Errors.
We are interested in the convergence rate of the two following quantities: Indeed, Theorems 4.1 and 4.2 provide upper bounds on such two quantities.

The p-version of the method
In this section, we present numerical results validating the theoretical predictions of Theorem 4.1 for the p-version of the method. We consider the exact solutions (u 1 , s 1 ) and (u 2 , s 2 ) in (58) and (60), respectively. We employ a coarse mesh of 2 × 2 uniform squares on the domain Ω 1 .  (58) and (60) in the left and right panel, respectively. We plot the errors s j − sn 0,Ω and |u j − un| 1,Tn , for j = 1, 2. We employ a coarse mesh of 2 × 2 uniform squares.
As expected from the theoretical predictions, in Figure 1, we observe exponential convergence for the test case with smooth solution, and only algebraic convergence for the singular solution case.
Remark 6. For the exact solution u 2 , the L 2 error on the pressure stagnates at around p = 4 and then grows. Similarly, the H 1 error stagnates starting from p = 9. This behaviour can be traced back to the ill-conditioning of the resulting linear system, which is mainly due to the choice of the polynomial bases in the definition of the degrees of freedom and in the expansion of the polynomial projectors. A possible remedy to this problem might be an orthogonalization process of the polynomial bases; see, e.g., [36]. For the sake of clarity, we avoid such investigation here.

The hp-version of the method
As predicted in Theorem 4.1 and observed in Figure 1 numerically, the method converges in terms of the degree of accuracy p algebraically, whenever the exact solution is not analytic. However, as discussed in Theorems 1.1 and 1.2, solutions to the Stokes problem (3) on polygonal domains with smooth data belong to the Kondrat'ev spaces K γ (Ω) in (5). In general, for solutions (u, s) to the Stokes problem in a nonconvex domain Ω, we can expect u ∈ H k (Ω) 2 and s ∈ H k−1 (Ω) for a given k < 2 only. Exponential convergence can be recovered for weighted analytic functions, by employing hpapproximation spaces, following the gospel of Babuška and collaborators, as proven in Theorem 4.2. See also, e.g., [3,4,41] and the references therein.
Thus, in this section, we validate the theoretical predictions of Theorem 4.2. To this aim, we consider the test case with exact solution (u 2 , s 2 ) in (60). We construct the distribution of the degrees of accuracy by picking µ = 1 in (57). Moreover, we employ hp-virtual element spaces based on geometric meshes as those depicted in Figure 2. There, we depict meshes with three layers, which are geometrically refined towards the re-entrant corner (0, 0) in three different ways. The numbers within the elements represent the local degrees of accuracy of the method. In Figures 3, 4, and 5, we depict the decay of the errors in (54) employing hp-virtual element spaces based on meshes as those in Figure 2. We pick different choices of the grading parameter σ. We employ meshes that are geometrically refined towards the re-entrant corner as those in Figure 2 (left). We pick three different choices of the parameter σ satisfying (55) and (56), namely σ = 1 2 , σ = √ 2 − 1, and σ = ( √ 2 − 1) 2 .
We observe exponential decay of the errors. The error saturation due to ill-conditioning manifests itself earlier for some of the meshes depicted in Figure 2. See Remark 6 for further details on this point.