The Stokes complex for Virtual Elements with application to Navier--Stokes flows

In the present paper, we investigate the underlying Stokes complex structure of the Virtual Element Method for Stokes and Navier--Stokes introduced in previous papers by the same authors, restricting our attention to the two dimensional case. We introduce a Virtual Element space $\Phi_h \subset H^2(\Omega)$ and prove that the triad $\{\Phi_h, V_h, Q_h\}$ (with $V_h$ and $Q_h$ denoting the discrete velocity and pressure spaces) is an exact Stokes complex. Furthermore, we show the computability of the associated differential operators in terms of the adopted degrees of freedom and explore also a different discretization of the convective trilinear form. The theoretical findings are supported by numerical tests.

It was soon recognized that the more general construction of VEM, that is not limited to polynomial functions on the elements, may allow for additional interesting features in additional to polytopal meshing. An example can be found in [10,11] where the authors developed a (conforming) Virtual Element Method for the Stokes and Navier-Stokes problems that guarantees a divergence free velocity, a property that yields advantages with respect to standard inf-sup stable schemes (see for instance [33]). And, most importantly, the proposed approach fitted quite naturally in the virtual element setting, so that the ensuing element is not particularly complicated to code or to handle.
Our aim is to further develop the idea in [10,11], also in order to get a deeper understanding of the underlying structure. In [26] the term "Stokes exact complex" was introduced; in that paper the authors underline that, if a given velocity/pressure FE scheme is associated to a discrete Stokes exact complex, than not only the existence of an unique solution is guaranteed, but also the divergence-free character of the discrete velocity. In addition, this allows to construct an equivalent curl formulation of the problem in a potential-like variable. This matter is one of interest in the FEM community, see for instance [36,33], also due to the difficulty in deriving exact Stokes complexes for Finite Elements, that often yield quite "cumbersome" schemes.
In the present paper, we unveil the underlying 2D Stokes complex structure of the VEM in [10,11] by introducing a Virtual Element space Φ h ⊂ H 2 (Ω) and proving that the triad {Φ h , V h , Q h } (with V h and Q h velocity and pressure spaces of [11]) is an exact Stokes complex. Furthermore, we show the computability of the associated differential operators in terms of the adopted degrees of freedom (a key aspect in VEM discretizations) and we explore also a different discretization of the convective trilinear form. As a byproduct of the above exactsequence construction, we obtain a discrete curl formulation of the Navier-Stokes problem (set in the smaller space Φ h ) that yields the same velocity as the original method (while the pressure needs to be recovered by solving a global rectangular system). For completeness, we also briefly present and compare a stream-function formulation approach, that is based on a direct discretization (with C 1 Virtual Elements) of the continuous stream function formulation of the problem. Some numerical tests are developed at the end of the paper, in order to show the performance of the methods, also comparing aspects such as condition number and size of the linear system. We note that a related study was developed in [3], but only for the lowest order case without enhancements (that is, suitable for Stokes but not for Navier-Stokes).
The paper is organized as follows. In Section 2 we review the Navier-Stokes problem in strong and variational form, together with some basic theoretical facts. In Section 4 (after introducing some preliminaries and definitions in Section 3) we review the Virtual scheme in [11], propose a third option for the discretization of the convective term and extend the convergence results also to this case. In Section 5 we introduce the space Φ h together with the associated degrees of freedom, prove the exact Stokes complex property and state the alternative curl formulation for the discrete problem. In Section 6 we present a set of numerical tests, that also compare the proposed method with a direct C 1 discretization of the stream-function problem, briefly described in the Appendix, that is not associated to a Stokes complex.
Throughout the paper, we will follow the usual notation for Sobolev spaces and norms [1]. Hence, for an open bounded domain ω, the norms in the spaces W s p (ω) and L p (ω) are denoted by · W s p (ω) and · L p (ω) respectively. Norm and seminorm in H s (ω) are denoted respectively by · s,ω and |·| s,ω , while (·, ·) ω and · ω denote the L 2 -inner product and the L 2 -norm (the subscript ω may be omitted when ω is the whole computational domain Ω). Moreover with a usual notation, the symbols ∇, ∆ and ∇ 2 denote the gradient, Laplacian and Hessian matrix for scalar functions, while ∆, ∇, and div denote the vector Laplacian, the gradient operator and the divergence for vector fields. Furthermore for a scalar function ϕ and a vector field v := (v 1 , v 2 ) we set

The Navier-Stokes equation
We consider the steady Navier-Stokes equation on a polygonal simply connected domain Ω ⊆ R 2 (for more details, see for instance [31]) where u, p are the velocity and the pressure fields, respectively, ν ∈ R, ν > 0 is the viscosity of the fluid and f ∈ [L 2 (Ω)] 2 represents the external force. For sake of simplicity we here consider Dirichlet homogeneous boundary conditions, different boundary conditions can be treated as well. Problem (1) can be written in the equivalent rotational form Systems (1) and (2) are equivalent in the sense that the velocity solutions u coincide and the rotational pressure solution P of Problem (2), the so-called Bernoulli pressure, and the convective pressure solution p of Problem (1) are jointed by the relation where, for the time being, λ denotes a suitable constant. Let us consider the spaces Let us introduce the bilinear forms It is well known (see, for instance, [31]) that in the diffusion dominated regime, i.e. under the assumption where C denotes the continuity constant of c(·; ·, ·) with respect to the V-norm, Problem (11) is well-posed and the unique solution (u, p) ∈ V × Q satisfies Finally, let us introduce the kernel of bilinear form b(·, ·), Then, Problem (11) can be formulated in the equivalent kernel form: In this case, from (9) and (10) it is straightforward to see that

Curl and Stream Formulation of the Navier-Stokes Equations
If Ω is a simply connected domain, a well known result (see [31] for the details) states that a vector function v ∈ Z if and only if there exists a scalar potential function ϕ ∈ H 2 (Ω), called stream function such that v = curl ϕ .
Clearly the function ϕ is defined up to a constant. Let us consider the space endowed with norm Then, Problem (15) can be formulated in the following curl formulation: A different approach that makes use of the stream functions is the following. Let ψ be the stream function of the velocity solution u of (1), i.e. u = curlψ. Then applying the curl operator to the equation (2), and using simple computations on the differential operators we obtain the equivalent following problem: This elliptic equation, can be reformulated in a variational way as follows, obtaining the socalled stream formulation (we refer again to [31]): where Since the formulations (17) and (19) are equivalent to Problem (15) (in turn equivalent to Problem (11)), the well-posedness of curl and stream formulations follows from assumption (A0). Moreover from (13) follows the stability estimate

Definitions and preliminaries
In the present section we introduce some basic tools and notations useful in the construction and theoretical analysis of Virtual Element Methods. Let { Ω h } h be a sequence of decompositions of Ω into general polygonal elements E with We suppose that for all h, each element E in Ω h fulfils the following assumptions: where is a uniform positive constant. We remark that the hypotheses above, though not too restrictive in many practical cases, can be further relaxed, as investigated in [9]. For any E ∈ Ω h , using standard VEM notations, for n ∈ N let us introduce the spaces: • P n (E) the set of polynomials on E of degree ≤ n (with the extended notation P −1 (E) = ∅), • B n (∂E) := {v ∈ C 0 (∂E) s.t v |e ∈ P n (e) for all edge e ⊂ ∂E}.
Remark 3.1. Note that (23) implies that the operator curl is an isomorphism from x ⊥ P n−1 (E) to the whole P n−1 (E), i.e. for any q n−1 ∈ P n−1 (E) there exists a unique p n−1 ∈ P n−1 (E) such that q n−1 = curl(x ⊥ p n−1 ) .
We also have that where n E is the number of edges (or the number of vertexes) of the polygon E. One core idea in the VEM construction is to define suitable (computable) polynomial projections. For any n ∈ N and E ∈ Ω h we introduce the following polynomial projections: • the L 2 -projection Π 0,E n : L 2 (E) → P n (E), given by for all v ∈ L 2 (E) and for all q n ∈ P n (E), with obvious extension for vector functions Π 0,E n : [L 2 (E)] 2 → [P n (E)] 2 , and tensor func- and for all q n ∈ P n (E), with obvious extension for vector functions Π ∇,E and for all q n ∈ P n (E), Finally, let us recall a classical approximation result for P n polynomials on star-shaped domains, see for instance [17]: with C depending only on n and the shape constant in assumptions (A1) and (A2).
In the following the symbol C will indicate a generic positive constant, independent of the mesh size h, the viscosity ν and the constant γ appearing in (A0), but which may depend on Ω, the integer k (representing the "polynomial" order of the method) and on the shape constant in assumptions (A1) and (A2). Furthermore, C may vary at each occurrence.

Virtual elements velocity-pressure formulation
In the present section we outline a short overview of the Virtual Element discretization of Navier-Stokes Problem (11). We will make use of various tools from the virtual element technology, that will be described briefly; we refer the interested reader to the papers [3,10,44,11]. We recall that in [10] a new family of Virtual Elements for the Stokes Equation has been introduced. The core idea is to define suitable Virtual space of velocities, associated to a Stokes-like variational problem on each element, such that the discrete velocity kernel is pointwise divergence-free. In [44] has been presented an enhanced Virtual space to be used in place of the original one, that, taking the inspiration from [2], allows the explicit knowledge of the L 2 -projection onto the polynomial space P k (being k the order of the method). In [11] a Virtual Element scheme for the Navier-Stokes equation in classical velocity-pressure formulation has been proposed. In the following we give some basic tools and a brief overview of such scheme. We focus particularly on the virtual element discretization of Navier-Stokes equation in rotation form (2) related to the trilinear form c rot (·; ·, ·) defined in (21) that was not treated in [11]. Specifically for the resulting discrete scheme we develop the convergence analysis for both the Bernoulli and the related convective pressure.

Virtual elements spaces
Let k ≥ 2 the polynomial degree of accuracy of the method. We consider on each element E ∈ Ω h the (enlarged) finite dimensional local virtual space Now, we define the Virtual Element space V E h as the restriction of U E h given by (cf. [11]): We here summarize the main properties of the virtual space V E h (we refer [44,11] for a deeper analysis): where n E is the number of vertexes of E; • degrees of freedom: the following linear operators D V , split into four subsets (see Figure 1) constitute a set of DoFs for V E h : -D V 1: the values of v at the vertexes of the polygon E, • projections: the DoFs D V allow us to compute exactly (c.f. (26) and (25)) in the sense that, given any v h ∈ V E h , we are able to compute the polynomials Π ∇,E k v h , Π 0,E k v h and Π 0,E k−1 ∇v h only using, as unique information, the degree of freedom values The global virtual element space is obtained as usual by combining the local spaces V E h accordingly to the local degrees of freedom, as in standard finite elements and considering the homogeneous boundary conditions. We obtain the global space The space V h has an optimal interpolation order of accuracy with respect to k (see Theorem 4.1 in [11]).
where the constant C depends only on the degree k and the shape regularity constant (see assumptions (A1) and (A2) of Section 3).
The pressure space is simply given by the piecewise polynomial functions with the obvious associated set of global degrees of freedom: A crucial observation is that, by definitions (30) and (31), it holds Therefore the discrete kernel is a subspace of the continuous kernel Z defined in (14), i.e.
This leads to a series of important advantages, as explored in [34,33,11,44]. Finally, we remark that the kernel Z h is obtained by gluing the local kernels explicitly defined by where n E is the number of vertexes of E.

Discrete bilinear forms and load term approximation
In this subsection we briefly describe the construction of a discrete version of the bilinear form a(·, ·) given in (4) and trilinear forms c(·; ·, ·) (cf. (6), (7), (8)). We can follow in a rather slavish way the procedure initially introduced in [7] for the laplace problem and further developed in [10,44,11] for flow problems. First, we decompose into local contributions the bilinear form a(·, ·) and the trilinear forms c(·; ·, ·) by considering: Therefore, following a standard procedure in the VEM framework, we define a computable discrete local bilinear form approximating the continuous form a E (·, ·), and defined by with α * and α * positive constants independent of the element E. For instance, a standard choice (26) and properties (39) imply that the discrete form a E h (·, ·) satisfies the consistency and stability properties. The global approximated bilinear form a h (·, ·) : V h × V h → R is defined by simply summing the local contributions: We now define discrete versions of the forms c(·; ·, ·). Referring to (6), (7), (8) we set for all and note that all quantities in the previous formulas are computable. Let c E h (·; ·, ·) be one of the discrete trilinear forms listed above. As usual we define the global approximated trilinear form by adding the local contributions: The forms c h (·; ·, ·) are immediately extendable to the whole V (simply apply the same definition for any w, u, v ∈ V). Moreover, the trilinear forms c h (·; ·, ·) are continuous on V, i.e. there exists a uniform bounded constant C h such that The proof of the continuity for the trilinear forms c conv,h (·; ·, ·) and c skew,h (·; ·, ·) can be found in Proposition 3.3 in [11]. Analogous techniques can be used to prove the continuity of the trilinear form c rot,h (·; ·, ·). For what concerns b(·, ·), as noticed in [10] we do not need to introduce any approximation of the bilinear form, since it can be exactly computed by the DoFs D V . The last step consists in constructing a computable approximation of the right-hand side (f , v) in (11). We define the approximated load term f h as and consider:

The discrete problem
Referring to (30), (31), (40), (44), (5) and (47), the virtual element approximation of the Navier-Stokes equation in the velocity-pressure formulation is given by: with c h (·; ·, ·) given by (41), (42) or (43). Whenever the choice (43) is adopted, the pressure output in (48) approximates the Bernoulli pressure P in (3) instead of the convective pressure p. Recalling the kernel inclusion (34), Problem (48) can be also formulated in the equivalent kernel form The well-posedness of the discrete problems can be stated in the following theorem (cf. [11]).

Theorem 4.2. Under the assumption
We have the following approximation results (see Theorem 4.6 and Remark 4.2 in [11] for the choices (41) and (42)).
for a suitable functions F, H, K independent of h.
Following the same steps as in [11], Theorem 4.3 easily extends also to the choice (43). In such case we preliminary observe that if the velocity solution u ∈ [H s+1 (Ω)] 2 and the convective pressure p ∈ H s (Ω) then the Bernoulli pressure P is in H s (Ω). As a matter of fact, recalling (3) by the Hölder inequality and Sobolev embedding H s+1 (Ω) ⊂ W s 4 (Ω), we recover Now let (u h , P h ) be the solution of the virtual problem (48) with the trilinear form (43) and (u, P ) be the solution of the Navier-Stokes equation (2). Then (53) is substituted by In such case the following computable approximation p h of the convective pressure p is available: where λ h is the mean value of the piecewise polynomial function 1 2 Π 0,E k u h ·Π 0,E k u h . The optimal order of accuracy for the convective pressure can established as follows. Definitions (3) (taking λ as the mean value of 1 2 u · u) and (55) easily imply where in the second inequality we have used the fact that all terms inside the norms are zero averaged. The first term in the right hand side of (56) is bounded by (54). Whereas for the terms µ E the triangular inequality and the Hölder inequality entail Using the Sobolev embedding H 1 (Ω) ⊂ L 4 (Ω), the continuity of the projection Π 0,E k with respect to the L 4 -norm and the H 1 -norm (see, for instance, [11]), and polynomial approximation estimates on star shaped polygons of Lemma 3.1, from (57) we infer Combining bound (58) with the Hölder inequality for sequences, the velocity error estimate (52) and with the stability estimates (13) and (51), it follows for a suitable functions L and I independent of h. Finally, inserting estimates (54) and (59) in (56) we obtain the optimal convergence result for the convective pressure also for choice (43). Remark 4.1. We observe that, due to the divergence-free property, the estimate on the velocity error in Theorem 4.3 does not depend on the continuous pressure, whereas the velocity errors of classical methods have a pressure contribution, see [11] for more details on this aspect.  [11]) shows that if the exact velocity solution u ∈ [P k (Ω)] 2 and the trilinear form c conv,h (·; ·, ·) or the trilinear form c rot,h (·; ·, ·) are adopted in (48), the corresponding schemes provide a higher order approximation errors, that are respectively These are to be compared with the error of standard inf-sup stable Finite Elements, that in the same situations would be O(h k ) due to the pressure contribution to the velocity error. Remark 4.3. Another interesting aspect related to method (48) is the so called "reduced" version. Noting that the discrete solution satisfies div u h = 0, one can immediately set to zero all D V 4 degrees of freedom, and correspondingly eliminate also the associated (locally zero average) discrete pressures. The resulting equivalent scheme has much less internal-to-element velocity DoFs and only piecewise constant pressures (we refer to [10,11] for more details).

Virtual elements Stokes complex and curl formulation
In the present part we present the VEM stream function space and the associated Stokes Complex.

Virtual element space of the stream functions
The aim of the present section is to introduce a suitable virtual space Φ h approximating the continuous space of the stream functions Φ defined in (16), such that In particular this will allow to exploit the kernel formulation (49) in order to define an equivalent VEM approximation for the Navier-Stokes equation in curl form (cf. (17)). Note that a related approach, but limited to a lowest order case and suitable only for the Stokes problem, was presented in [3]. In order to construct the space of the virtual stream functions Φ h , we proceed step by step, following the enhanced technique used in Subsection 4.2. On each element E ∈ Ω h we consider the enlarged local virtual space: (61) Then we define the enhanced space of the stream functions (62) It is straightforward to see that P k+1 (E) ⊆ Φ E h . We are now ready to introduce a suitable set of degrees of freedom for the local space of virtual stream functions Φ E h . Given a function ϕ ∈ Φ E h , we take the following linear operators D Φ , split into five subsets (see We note that the linear operators D Φ 1 and D Φ 2 are always needed to enforce the C 1continuity at the vertices. Moreover it is immediate to check that for any stream function ϕ ∈ Φ E h (the same holds for Ψ E h ), the linear operator evaluations of D Φ 1, D Φ 2, D Φ 3, D Φ 4 uniquely determine ϕ and its gradient on ∂E. We now prove the following result.
where as usual n E denotes the number of edges of the polygon E. Proof. We preliminary prove that the linear operators D Φ plus the additional moments of curl ϕ constitute a set of degrees of freedom for Ψ E h . An integration by parts and recalling Remark 3.1 imply that the linear operators D Ψ 5 + D Φ 5 are equivalent to prescribe the moments E ϕ q k−1 dE for all q k−1 ∈ P k−1 (E). Indeed, Remark 3.1 and simple computations give where the boundary integral is computable using the DoFs values. Now the assertion easily follows by a direct application of Proposition 4.1 in [18]. In particular, from (24) it holds that The next step is to prove that the linear operators D Φ are unisolvent for Φ E h . From (24) it holds dim P k−1\k−3 (E) = dim (P k−1 (E)) − dim (P k−3 (E)) = 2 k − 1 .
Hence, neglecting the independence of the additional (2 k − 1) conditions in (62), the dimension of Φ E h is bounded from below by Due to (65), the proof is concluded if we show that a function ϕ ∈ Φ E h such that D Φ (ϕ) = 0 is identically zero. In such case, ϕ = 0 and ∇ϕ = 0 on the skeleton ∂E and this entails curl ϕ = 0 on ∂E. Moreover we note that in this case the Π ∇,E k (curl ϕ) = 0; as a matter of fact, by definition (26), we get The boundary integral is zero being curl ϕ = 0 on the skeleton ∂E. For the internal integral, in the light of (23), let us set where the boundary integral is zero since ϕ = 0 on ∂E, whereas the second term is zero since D Φ 5(ϕ) = 0. In particular we proved that, since Π ∇,E k (curl ϕ) = 0, recalling (62) also the moments D Ψ 5 of ϕ are zero. Since D Φ (ϕ) = 0 by assumption, recalling that ϕ ∈ Φ E h ⊂ Ψ E h and that {D Φ , D Ψ 5} are a set of degrees of freedom for Ψ E h , it follows ϕ = 0.
1. An alternative way to define a unisolvent set of DoFs for the space Φ E h is to provide in the place of D Φ 5 the following operators [18] • D Φ 5 : the moments of ϕ against the polynomial of degree up to degree k − 3 but such choice is less suitable for the exact sequence construction of the present work.
The global virtual space Φ h is obtained by combining the local spaces Φ E h accordingly to the local degrees of freedom, taking into account the boundary conditions: , where n P (resp., n e and n V ) is the number of elements (resp., internal edges and vertexes) in the decomposition Ω h .

Virtual element Stokes complex
The aim of the present subsection is to provide a virtual element counterpart of the continuous Stokes complex [32]: where i denotes the mapping that to every real number r associates the constant function identically equal to r and we recall that a sequence is exact if the image of each operator coincides with the kernel of the following one. The case without boundary conditions is handled analogously. We start by characterizing the space Φ E h as the space of the stream functions associated to the discrete kernel Z E h . Proposition 5.2. For any E ∈ Ω h let Z E h and Φ E h be the spaces defined in (35) and (62), respectively. Then, it holds that Concerning the condition on the skeleton ∂E, we observe that ϕ h|∂E ∈ B k+1 (∂E) and Inside the element, by simple calculations and by definition (62), we infer curl ∆v h = curl ∆(curl ϕ h ) = ∆ 2 ϕ h ∈ P k−1 (E). In the light of Remark 3.1, the previous relation is equivalent to Therefore, there exists p k−1 ∈ P k−1 (E) such that curl(∆v h − x ⊥ p k−1 ) = 0. Since E is simply connected, there exists s such that ∆v h − x ⊥ p k−1 = ∇s. Thus we have shown that Moreover, by definition (62), for all At this point is clear that (69), (70), (71), (72) and definition (35), The proof now follows by a dimensional argument. In fact, from (63) and (36) easily follows that Therefore we can conclude that curl h , from the degrees of freedom values D Φ of ϕ h , we are able to compute the DoFs values D V of curl ϕ h . In particular it holds that and Therefore, for any ϕ h ∈ Φ E h , the DoFs D Φ allow to compute the polynomial projections . As a consequence of Proposition 5.2 we have the following Stokes exact sequence for our discrete VEM spaces and its reduced version (see also Figures 3 and 4). (62) and (28), respectively, and let V E h denote the reduced velocity space, see Remark 4.3. Then, the following sequences are exact Remark 5.3. In terms of degrees of freedom, our lowest order element (when restricted to triangles!) can be compared with the Zienkiewicz element [23,32], all other FEM elements in the literature being either higher order or needing a sub-element partition and more DoFs. The reduced version of our VEM element for k = 2 (see Remark 4.3 and Figure 4) has piecewise constant pressures and no internal degrees of freedom for velocities, and thus in terms of degrees of freedom exactly corresponds to the above finite element (the difference is that we use the VE approach instead of introducing rational basis functions). But note that the element here presented yields O(h 2 ) convergence rate for velocities and also for the local pressure average (full O(h 2 ) pressure convergence can be recovered by a local post-processing), instead of linear convergence as [32]. In addition, we avoid integration of rational functions. Clearly, this comes at the price of having a virtual formulation and thus the absence of an explicit expression of the shape functions. The following results are the global counterpart of Proposition 5.2 and Corollary 5.1.

Proposition 5.3.
Let Z h and Φ h be the spaces defined in (33) and (67), respectively. Then, it holds that Proof. We note that Proposition 5.2 endowed with the boundary condition in the definitions (33) and (67) imply curl Φ h ⊆ Z h . The proof now follows by a dimensional argument using the Euler formula. (67), (30) and (31), respectively, and let V h and Q h denote the reduced velocity space and the piecewise constant pressures, respectively, see Remark 4.3. Then, the following sequences are exact

Corollary 5.2. Let Φ E h , V h and Q h be the spaces defined in
The case without boundary conditions follows analogously.

The discrete problem
In the light of Proposition 5.2, referring to (49) and (67) we can set the virtual element approximation of the Navier-Stokes equation in the curl formulation: (77) Due to Proposition 5.3, Problem (77) is equivalent to (49). We remark that all forms in (77) are exactly computable by the DoFs D Φ . In fact, recalling Remark 5.2, the polynomials are computable on the basis of D Φ , so that, referring to (40), (44) and (47), we infer that are exactly computable from DoFs D Φ , see also The convergence of the discrete solution curl ψ h of (77) to the continuous solution curl ϕ of (17) follows immediately from Theorem 4.3, taking u = curl ϕ and u h = curl ϕ h .
Clearly Problem (77) does not provide any information on the pressure p. Nevertheless, the Stokes complex associated to the proposed scheme turns out to be very helpful if we are interested in computing suitable approximation p h of p. Indeed referring to (30) and (48), starting from the solution ψ h of Problem (77), we infer the following problem ), the previous system, is actually an overdetermined system, i.e. there are more equations than unknowns. Nevertheless the well-posedness of Problem (79) is guaranteed by Theorem 4.2. We refer to Section 6 for a deeper analysis and computational aspects of (79).
We stress that the curl virtual formulation (77) exhibits important differences from the computational point of view compared with the velocity-pressure formulation (48). First of all the linear system associated to Problem (77) has 2(n P − 1) less DoFs than Problem (48), even if considering its equivalent reduced form (see Remark 4.3). Moreover the first iteration of the Newton method applied to the the non-linear virtual stream formulation (77) results in a linear system which is symmetric and positive definite, whereas applied to the virtual element method (48) in velocity-pressure formulation leads to an indefinite linear system. These advantages come at the price of a higher condition number of the involved linear systems. Remark 5.4. Simple integration by parts gives (f , curl ϕ) = (curl f , ϕ). By Remark 5.1, the DoFs D Φ allow us to compute the L 2 -projection Π 0,E k−1 : Φ E h → P k−1 (E), so that we can consider a new computable right hand-side This new formulation of the right-hand side gets the same order of accuracy of the original one.
In particular if the external force is irrotational, i.e. f = ∇f , we improve the error estimate in (52) by removing the dependence of the error by the load. More generally, with the choice (80), we completely remove the influence in the error stemming from the irrotational part in the Helmholtz decomposition of the load. Clearly (80) can be applied only when f is given as an explicit function.

Numerical Tests
In this section we present two sets of numerical experiments to test the practical performance of the proposed virtual element methods (77), also compared with a direct C 1 VEM discretization of the stream formulation (18) described in the Appendix, see equation (87). For the scheme (77), in all tests we investigate the three possible options for the trilinear form in (41), (42), (43). In Test 6.1 we study the convergence of the proposed virtual element schemes (77) and (87) for the discretization of the Navier-Stokes equation in curl formulation and stream formulation respectively. A comparison of (77) (in terms of errors, number of DoFs, condition number of the resulting linear systems) with the equivalent virtual element scheme (48) for the Navier-Stokes equation in velocity-pressure formulation is also performed. In Test 6.2 we consider a benchmark problem for the Navier-Stokes equation (18) with the property of having the velocity and stream solution in the corresponding discrete spaces. It is well known that classical mixed finite element methods lead to significant velocity errors, stemming from the velocity/pressure coupling in the error estimates. This effect is greatly reduced by the presented methods (cf. Theorem 4.3, estimate (52) and Remark 4.2). In order to compute the VEM errors, we consider the computable error quantities: for the velocity-pressure formulation (48) and for the curl and stream formulations (see (77) and (87), respectively). For what concerns the pressures we simply compute the standard L 2 error error(p, L 2 ) := p − p h 0 . For the computation of the discrete pressure for the virtual element scheme (77) we follow (79) and solve the overdetermined system by means of the least squares method. We briefly sketch the construction of the least square formula.
be the canonical basis functions of V h and let us denote with r h the vector with component i.e. r h contains the values of the degrees of freedom associated to the right hand side of (79) with respect to the basis be the canonical basis functions of Q h and for any piecewise polynomial p k−1 ∈ Q h we denote with p h the vector containing the values of the coefficients with respect to the basis {q i } associated to p k−1 . Then the least squares formula associated to (79) is An example of the adopted meshes is shown in Figure 5. For the generation of the Voronoi meshes we use the code Polymesher [42]. The distorted quadrilateral meshes are obtained starting from the uniform square meshes and displacing the internal vertexes with a proportional "distortion amplitude" of 0. 3     We test the virtual element scheme (77). In Figures 7 and 8 we show the results obtained with the sequences of Voronoi meshes V h and quadrilateral meshes Q h , by considering the three possible choices of the trilinear forms. We stress that in all cases considered we compare the discrete convective pressure p h (for the trilinear form in c rot,h (·; ·, ·) we consider the definition (55)).
We notice that the theoretical predictions of Section 5 are confirmed. Moreover, we observe that the virtual element methods obtained with the three different trilinear forms exhibit almost identical results, at least for this example and with the adopted meshes. In Figures 7 (left) and 8 (left) we also depict the error for the direct C 1 discretization of the stream formulation (87), that follows a similar behaviour to (77). Note that we do not compute a pressure error for scheme (87) since the computation of a discrete pressure is a more complex issue in this case, see Remark 7.1. Finally we test the corresponding virtual element method (48) with the same sequences of polygonal meshes V h , Q h . Table 1 shows the results obtained respectively with VEM (48) and (77) obtained considering the trilinear form c rot,h (·; ·, ·). The results are analogous also for the other two proposed trilinear forms (not shown). In Table  2 we compare the number of DoFs and the condition number of the resulting linear systems  V h 1/8 3.704032467e-1 3.891840615e-1 3.704032467e-1 3.891840615e-1 1/16 9.153568669e-2 8.875084726e-2 9.153568669e-2 8.875084726e-2 1/32 2.308710367e-2 1.994452869e-2 2.308710367e-2 1.994452869e-2 1/64 5.791512013e-3 4.602515029e-3 5.791512013e-3 4.602515029e-3 Q h 1/10 3.047752518e-1 3.714633884e-1 3.047752518e-1 3.714633884e-1 1/20 8.709526360e-2 8.363888240e-2 8.709526360e-2 8.363888240e-2 1/40 2.188243443e-2 1.945853612e-2 2.188243443e-2 1.945853612e-2 1/80 5.523374104e-3 4.762632907e-3 5.523374104e-3 4.762632907e-3 (stemming from the fist iteration of the Newton method) for both formulations (48) and (77). As observed in Section 5, the scheme (77) has the advantage of having (2 n p − 2) less of unknowns, even when considering the reduced version (see Remark 4.3) for formulation (48). The drawback is that the condition number of the system resulting from the velocity-pressure scheme (48) behaves as h −2 , while the asymptotic rate of the condition number of the linear system resulting from the scheme (77)     We here mention only the main properties of the virtual space Φ E h and refer to [18,4] for a deeper description: , that is the same dimension of Φ E h (cf. (63)); • degrees of freedom: the linear operators D Φ : D Φ 1, D Φ 2, D Φ 3, D Φ 4, D Φ 5 (see Remark 5.1) constitute a set of DoFs for Φ E h ; • projections: the DoFs D Φ allow us to compute exactly (c.f. (27) and (25)) in the sense that, given any ϕ h ∈ Φ E h , we are able to compute the polynomials Π ∇ 2 ,E k+1 ϕ h , Π 0,E k−1 ∆ϕ h , Π 0,E k−1 ϕ h , Π 0,E k−1 ∇ϕ h , and Π 0,E k−1 curlϕ h only using, as unique information, the degree of freedom values D Φ of ϕ h .
The global virtual element space is obtained as usual by combining the local spaces Φ E