The Stokes complex for Virtual Elements with application to Navier--Stokes flows
L. Beirão da Veiga D. Mora G. Vacca









































































0Π^0, E_k

Π^∇, E_k
Π^∇^2, E_k+1
0Π^0, E_k-1
𝐠_k-2^⊕
𝐠_k^⊕
𝒢_k-2^⊕(E)
𝒢_k^⊕(E)

sistema
{[ .
definition
definizioneDefinition[section]

remark
remarkRemark[section]

remark
testTest[section]

plain
theoremTheorem[section]
propositionProposition[section]
corollarioCorollary[section]
lemmaLemma[section]













citazione




1]L. Beirão da Veiga lourenco.beirao@unimib.it

2,3]D. Mora dmora@ubiobio.cl

1]G. Vacca giuseppe.vacca@unimib.it

[1]Dipartimento di Matematica e Applicazioni,
Università degli Studi di Milano Bicocca,
Via Roberto Cozzi 55 - 20125 Milano, Italy  

[2]Departamento de Matemática, 
Universidad del Bío-Bío,
Casilla 5-C, Concepción, Chile 

[3]CI^2MA, Universidad de Concepción, Concepción, Chile














 The Stokes complex for Virtual
Elements with application to Navier–Stokes flows
    [
    December 30, 2023
================================================================================


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 Φ_h ⊂ H^2(Ω) and
prove that the triad {Φ_h, _h, Q_h}
(with _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.







§ INTRODUCTION




The Virtual Element Method (VEM) was introduced in <cit.> as a generalization of the Finite Element Method that allows to use general polygonal and polyhedral meshes. The VEM enjoyed a good success in this recent years both in the mathematics and in the engineering communities, a very brief list of papers being 
<cit.>,

while for the specific framework of incompressible flows we refer to <cit.>.

Virtual Elements are part of a wider community of
polytopal methods, that include other technologies
(some representative papers being 
<cit.>).

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 <cit.>
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 <cit.>). 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 <cit.>,
also in order to get a deeper understanding of the underlying structure.
In <cit.> 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  formulation of the problem in a potential-like variable.
This matter is one of  interest in the FEM community, see for instance
<cit.>, 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 <cit.>
by introducing a Virtual Element space Φ_h ⊂ H^2(Ω) and proving that the triad {Φ_h, _h, Q_h}  (with _h and Q_h velocity and pressure spaces of <cit.>) 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 exact-sequence construction, we obtain a discrete  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 <cit.>, 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 <ref> we review the Navier–Stokes problem in strong and variational form, together with some basic theoretical facts. In Section <ref> (after introducing some preliminaries and definitions in Section <ref>) we review the Virtual scheme in <cit.>, propose a third option for the discretization of the convective term and extend the convergence results also to this case. In Section <ref> 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 <ref> 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 <cit.>.
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  denote the vector Laplacian,  the gradient operator
and the divergence  for vector fields.

Furthermore for a scalar function ϕ and a vector field := (v_1,   v_2) we set

    ϕ := ( ∂ϕ/∂ y,   -∂ϕ/∂ x) 
        and      := ∂ v_2/∂ x  -  ∂ v_1/∂ y ,
    
    
    ϕ× := ϕ( -v_2,   v_1)  .














§ THE NAVIER–STOKES EQUATION




We consider the steady Navier–Stokes equation on a polygonal
simply connected domain Ω⊆^2
(for more details, see for instance <cit.>) 



    {    find (,p) such that      
        -  ν  + ( )   -  ∇ p =           in Ω,
         = 0           in Ω,
        = 0            on ∂Ω,.



where , p are the velocity and the pressure fields, respectively, 
ν∈ℝ, ν > 0 is the viscosity of the fluid and ∈ [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 (<ref>) can be written in the equivalent rotational form

    {         
        -  ν  + ( ) × -  ∇ P =           in Ω,
         = 0           in Ω,
        = 0            on ∂Ω..


Systems (<ref>) and (<ref>)
are equivalent in the sense that the velocity solutions 
coincide and the rotational pressure solution P of Problem (<ref>),
the so-called Bernoulli pressure, and the convective
pressure solution p of Problem (<ref>) are jointed by the relation

    P := p - 1/2 ·  +  λ ,


where, for the time being, λ denotes a suitable constant.

Let us consider the spaces

    := [ H_0^1(Ω) ]^2,      Q:= L^2_0(Ω) = { q ∈ L^2(Ω)   s.t.  ∫_Ω q   dΩ = 0 }

endowed with natural norms

    _ := | |_ H^1(Ω) ,     q_Q :=  q_L^2(Ω).


Let us introduce the bilinear forms

    a(·, ·) ×→,         
    a (𝐮,  𝐯) := ∫_Ω ∇𝐮 : ∇𝐯  dΩ,     for all 𝐮,  𝐯∈𝐕,
    
    b(·, ·) 𝐕× Q →     b(𝐯, q) :=  ∫_Ωq  div 𝐯  dΩ    for all 𝐯∈𝐕, q ∈ Q,


and the trilinear forms


    (·;  ·, ·)    ××→        (;  , ) :=  ∫_Ω (  )  ·  dΩ ,
    
    (·;  ·, ·)    ××→        (;  , ) :=  
    1/2(;  , ) -
    1/2(;  , )  ,
    
    (·;  ·, ·)    ××→        (;  , ) :=  ∫_Ω ( × ) ·  dΩ .


Direct computations give that

    (;  , ) = (;  , )     for all , ∈ with= 0,
    (;  , ) = (;  , ) + 
    1/2∫_Ω (·)    dΩ    for all , ∈.


In the following we denote with c(·;  ·, ·) one of the trilinear forms listed above.

Then a standard variational formulation of Problem (<ref>) is:

    {   find (, p) ∈× Q, such that
       ν  a(, ) + c(;  , ) + b(, p) = (, )        for all ∈,
         b(, q) = 0        for all q ∈ Q,.


where

    (, ) := ∫_Ω·  dΩ .


From (<ref>) it is clear that if
c(·;  ·, ·) = (·;  ·, ·)
we recover instead the variational formulation of
system (<ref>) and the pressure solution
p is actually the  Bernoulli pressure P where λ
in (<ref>) is the mean value of 1/2 ·.



It is well known (see, for instance, <cit.>)
that in the diffusion dominated regime, i.e. under the assumption

    (𝐀0)    γ := C _H^-1/ν^2 < 1

where C denotes the continuity constant of
c(·;  ·, ·) with respect to the -norm,

Problem (<ref>) is well-posed and the
unique solution (𝐮, p) ∈𝐕× Q satisfies

    _≤_H^-1/ν .






































Finally, let us introduce the kernel of bilinear form b(·,·),

    := {∈  s.t.   b(, q) = 0   for all q ∈ Q} .


Then, Problem (<ref>) can be formulated in the equivalent kernel form:

    {   find ∈, such that
       ν  a(, ) + c(;  , ) = (, )        for all ∈..


In this case, from (<ref>) and (<ref>) it is straightforward to see that

    (;  , ) = (;  , ) = (;  , )
        for all , ∈.








 §.§ Curl and Stream Formulation of the Navier–Stokes Equations




If Ω is a simply connected domain, 
a well known result (see  <cit.> for the details) 
states that
a vector function ∈ if and only if there exists a scalar potential function ϕ∈ H^2(Ω), called stream function such that 

    =  ϕ .
 
Clearly the function ϕ is defined up to a constant.
Let us consider the space

    Φ := H^2_0(Ω) = {ϕ∈ H^2(Ω)   s.t  ϕ = 0,   ∂ϕ/∂ n = 0 
      on ∂Ω}
 

endowed with norm

    ϕ_Φ := |Φ|_H^2(Ω)    for all ϕ∈Φ.


Then, Problem (<ref>) can be
formulated in the following  formulation:


    {   find ψ∈Φ, such that
       ν  a(ψ,  ϕ) + c(ψ;  ψ, ϕ) = (, ϕ)
           for all ϕ∈Φ..





A different approach that makes use of the stream functions is the following.
Let ψ be the stream function of the velocity  solution  of (<ref>),  i.e. = ψ.
Then applying the  operator to the equation (<ref>), and 
using simple computations on the differential operators  
we obtain the equivalent following problem: 

    {         
       ν Δ^2 ψ -   ψ·(Δψ) =            in Ω,
       ψ = 0            on ∂Ω,
       ∂ψ/∂ n = 0            on ∂Ω..
 

This elliptic equation, can be reformulated in a variational
way as follows, obtaining the so-called stream formulation
(we refer again to <cit.>):


    {    find ψ∈Φ such that      
       ν a(ψ,  ϕ) 
      +   c(ψ;  ψ,  ϕ)  = ( ,  ϕ) 
              for all  ϕ∈Φ,.
 

where

    a(·, ·) Φ×Φ→,        a (ψ,  ϕ) := ∫_Ω Δ ψ Δ ϕ  dΩ,     for all ψ,  ϕ∈Φ,
    c(·;  ·, ·) Φ×Φ×Φ→  c(ζ;  ψ, ϕ) :=  ∫_ΩΔζ ψ·ϕ  dΩ  for all ζ, ψ, ϕ∈Φ.


Since the formulations (<ref>) and (<ref>)  are equivalent to Problem (<ref>) (in turn equivalent to Problem (<ref>)), the well-posedness of  and stream formulations  follows from assumption (𝐀0). Moreover from (<ref>) follows the stability estimate

    ψ_Φ≤_H^-1/ν .

 

























































































§ 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

    h_E :=  diameter(E) ,   
    h := sup_E ∈Ω_h h_E .


We suppose that for all h, each element E in Ω_h fulfils the following assumptions:


  (𝐀1) E is star-shaped with respect to a ball B_E of radius ≥ ρ  h_E, 

  (𝐀2) the distance between any two vertexes of E is ≥ρ  h_E, 

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  <cit.>. 


For any E ∈Ω_h, using standard VEM notations, for n ∈  let us introduce the spaces:



  * _n(E) the set of polynomials on E of degree ≤ n  (with the extended notation _-1(E)=∅),

  * _n(∂ E) := {v ∈ C^0(∂ E)   s.t   v_|e∈_n(e)   for all edge e ⊂∂ E}.






Note that

    [_n(E)]^2 = ∇(_n+1(E))  ⊕𝐱^⊥ _n-1(E)

 with 𝐱^⊥:= (x_2, -x_1).






















Note that (<ref>) implies that
the operator  is an isomorphism from 𝐱^⊥ _n-1(E) to the whole _n-1(E), i.e. for any q_n-1∈_n-1(E) there exists a unique p_n-1∈_n-1(E) such that

    q_n-1 =  (𝐱^⊥ p_n-1 )  .

 



We also have that

    (_n(E) ) = n(n+1)/2 ,
        (_n(∂ E)) = n_E   n


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 ∈ and E ∈Ω_h we introduce the following polynomial projections:




  * the L^2-projection Π_n^0, E L^2(E) →_n(E), given by

    ∫_Eq_n (v -  Π_n^0, E  v)   d E = 0     for all v ∈ L^2(E)  and for all q_n ∈_n(E),
 
with obvious extension for vector functions Π_n^0, E [L^2(E)]^2 → [_n(E)]^2, and tensor functions 
Π_n^0, E [L^2(E)]^2 × 2→ [_n(E)]^2 × 2,


  * the H^1-seminorm projection Π_n^∇,E H^1(E) →_n(E), defined by 

    {   ∫_E  q_n · ( v -  Π_n^∇,E   v)   d E = 0     for all v ∈ H^1(E) and for all q_n ∈_n(E),
       Π_0^0,E(v -  Π_n^∇, E  v) = 0   ,
    .
 
with obvious extension for vector functions Π_n^∇, E [H^1(E)]^2 → [_n(E)]^2,



  * the H^2-seminorm projection Π_n^∇^2,E H^2(E) →_n(E), defined by 

    {   ∫_E ∇^2 q_n : ∇^2 (v-  Π_n^∇^2,E   v)   d E = 0     for all v ∈ H^2(E) and for all q_n ∈_n(E),
       ∫_∂ Ev -  Π_n^∇^2, E  v    ds= 0      and    ∫_∂ E∂/∂ n(v -  Π_n^∇^2, E  v)    ds= 0  .
    .
 










Finally, let us recall a classical approximation result for _n
polynomials on star-shaped domains, see for instance <cit.>:



Under the assumptions (𝐀1) and (𝐀2),
for any v ∈ H^s+1(E), with s ≥ 0, it holds

    v - Π_n^0, E v_m, E≤ C   h^s +1 - m_E   |v|_s+1, E    for m ≤ s + 1 ≤ n+1,
    v - Π_n^∇, E v _m, E≤ C   h^s + 1- m_E   |v|_s+1, E    for m ≤ s + 1 ≤ n+1,

with C depending only on n and the shape  constant ρ in assumptions (𝐀1)
and (𝐀2).
 


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 (𝐀0), but which may depend on Ω,
the integer k (representing the “polynomial” order of the method) and
on the shape constant ρ in assumptions (𝐀1)
and (𝐀2). 
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 (<ref>). 
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
<cit.>.

We recall that in <cit.>  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 <cit.> has been presented an enhanced Virtual space
to be used in place of the original one, that, taking the inspiration from <cit.>,
allows the explicit knowledge of the L^2-projection onto
the  polynomial space _k (being k the order of the method). 

In <cit.> 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 (<ref>)
related to the trilinear form (·;  ·, ·)
defined in (<ref>) that was not treated in <cit.>. 
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



















    𝐔_h^E := {𝐯∈ [H^1(E)]^2   s.t  𝐯_|∂ E∈ [_k(∂ E)]^2   , .
    
    .
    {    - Δ𝐯  -  ∇ s ∈𝐱^⊥ _k-1(E),  
        div 𝐯∈_k-1(E),
    .      for some s ∈  L^2(E) ∖  }.

Now, we define the Virtual Element space 𝐕_h^E
as the restriction  of 𝐔_h^E given by (cf. <cit.>):

    𝐕_h^E := {𝐯∈𝐔_h^E   s.t.  (𝐯 - Π^∇,E_k 𝐯,  𝐱^⊥ p_k-1)_E = 0   for all p_k-1∈_k-1 ∖ k-3(E)} ,

where the symbol _k-1 ∖ k-3(E) := _k-1(E) ∖_k-3(E)
denotes the polynomials in _k-1(E) that are L^2-orthogonal
to all polynomials of _k-3(E).




We here summarize the main properties of the virtual
space _h^E  (we refer  <cit.> for a deeper analysis):



  * dimension:
 the dimension of _h^E  is
 
    ( 𝐕_h^E ) 
       = 2n_E   k + (k-1)(k-2)/2  + (k+1)k/2 - 1

where n_E is the number of vertexes of E;


  * degrees of freedom:
the following linear operators 𝐃_𝐕, split into four subsets (see Figure <ref>) constitute a set of DoFs for _h^E:




  * 𝐃_𝐕1:  the values of 𝐯 at the vertexes of the polygon E,


  * 𝐃_𝐕2: the values of 𝐯 at k-1 distinct points of every edge e ∈∂ E,


  * 𝐃_𝐕3: the moments of 𝐯

    ∫_E 𝐯·𝐱^⊥  p_k-3  dE     for all p_k-3∈_k-3(E),



  * 𝐃_𝐕4: the moments of div 𝐯 


    ∫_E ( div 𝐯)   q_k-1  dE     for all q_k-1∈_k-1(E) ∖;
 


  * projections:
the DoFs 𝐃_𝐕 allow us to compute exactly (c.f. (<ref>) and (<ref>))

    _h^E → [_k(E)]^2,     0 _h^E → [_k(E)]^2,     0 (_h^E) → [_k-1(E)]^2 × 2,

in the sense that, given any _h ∈_h^E, we are able to compute the polynomials
_h, 0 _h and 0∇_h only using, as unique information, the degree of freedom values 𝐃_𝐕 of _h.




The global virtual element space is obtained as usual by combining the local spaces _h^E
accordingly to the local degrees of freedom, as in standard  finite elements and considering the homogeneous boundary conditions. We obtain the global space


    𝐕_h := {𝐯∈ [H^1_0(Ω)]^2    s.t  𝐯_|E∈𝐕_h^E    for all E ∈Ω_h} .
 

The space _h has an optimal interpolation order of
accuracy with respect to k (see Theorem 4.1 in <cit.>).



Let ∈ H^s+1(Ω) ∩, for 0<s ≤ k.
Then, there exists _I ∈_h such that

    - _I _0 + h   - _I _≤ C   h^s+1  |  |_s+1,

where the constant C depends only on the degree k
and the shape regularity constant ρ (see assumptions (𝐀1)
and (𝐀2) of Section <ref>).


The pressure space is simply given by the piecewise polynomial functions

    Q_h := { q ∈ L_0^2(Ω)   s.t.   q_|E∈_k-1(E)   for all E ∈Ω_h},

with the obvious associated set of global degrees of freedom:



  * 𝐃_𝐐: the moments up to order k-1 of q, i.e.

    ∫_E q   p_k-1  dE     for all p_k-1∈_k-1(E)
    and for all E ∈Ω_h.



















































A crucial observation is that, by definitions (<ref>) and (<ref>), it holds 

    div 𝐕_h⊆ Q_h  .

Therefore the discrete kernel


    𝐙_h := {𝐯_h ∈𝐕_h   s.t.   b(𝐯_h, q_h) = 0   for all q_h ∈ Q_h},

is a subspace of the continuous kernel  defined in (<ref>), i.e.


    _h ⊆ .



This leads to a series of important advantages, as explored in 
<cit.>.

Finally, we remark that the kernel _h is obtained
by gluing the local kernels explicitly defined by

    𝐙_h^E := {𝐯∈ [H^1(E)]^2   s.t  (𝐯 - Π^∇,E_k 𝐯,  𝐱^⊥ p_k-1)_E = 0   for all p_k-1∈_k-1 ∖ k-3(E).
    
    .
    𝐯_|∂ E∈ [_k(∂ E)]^2   , 
      {    - Δ𝐯  -  ∇ s ∈𝐱^⊥ _k-1,  
        div 𝐯 = 0,
    .      for some s ∈ L^2(E) ∖  } .


The dimension of _h^E  is

    ( 𝐙_h^E ) 
       = 2n_E k + (k-1)(k-2)/2  - 1.

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 (<ref>)
and trilinear forms c(·; ·, ·)
(cf. (<ref>), (<ref>), (<ref>)).

We can follow in a rather slavish way the procedure initially introduced in <cit.> for the laplace problem and further developed in <cit.> for flow problems.

First, we decompose into local contributions the bilinear
form a(·, ·) and the trilinear forms c(·; ·, ·) by considering:


    a (𝐮,  𝐯) =: ∑_E ∈Ω_h a^E (𝐮,  𝐯)  ,
        
    c(;  , ) =: ∑_E ∈Ω_h c^E(;  , )      for all , , ∈.



Therefore, following a standard procedure in the VEM framework, we define a computable discrete local bilinear form


    a_h^E(·, ·) 𝐕_h^E ×𝐕_h^E →

approximating the continuous form a^E(·, ·), and defined by

    a_h^E(𝐮_h, 𝐯_h) := a^E (𝐮_h,  𝐯_h ) + 𝒮^E ((I - ) 𝐮_h,   (I -) 𝐯_h )

for all 𝐮_h, 𝐯_h ∈𝐕_h^E, where the (symmetric) stabilizing bilinear form 𝒮^E 𝐕_h^E ×𝐕_h^E → satisfies 

    α_* a^E(𝐯_h, 𝐯_h) ≤𝒮^E(𝐯_h,
    𝐯_h) ≤α^* a^E(𝐯_h, 𝐯_h)     ∀ 𝐯_h ∈𝐕_h^E,   Π_k^∇ ,E𝐯_h= 0

with α_* and α^*  positive  constants independent of the element E.

For instance, a standard choice is 
𝒮^E(𝐮_h, 𝐯_h) = ∑_i=1^N_DoFs𝐃_𝐕_i(𝐮_h)  𝐃_𝐕_i(𝐯_h)

where 𝐃_𝐕_i(𝐮_h) denotes the i-th  DoFs value of 
_h, opportunely scaled.

It is straightforward to check that the definition of Π^∇, E_k projection  (<ref>) and properties (<ref>) imply 
that the discrete form a_h^E(·, ·)
satisfies the consistency and stability properties.

The global approximated bilinear form a_h(·, ·) 𝐕_h ×𝐕_h → is defined by simply summing the local contributions:


    a_h(𝐮_h, 𝐯_h) := ∑_E ∈Ω_h  a_h^E(𝐮_h, 𝐯_h)     for all 𝐮_h, 𝐯_h ∈𝐕_h.



We now define discrete versions of the forms  c(·;  ·, ·). 
Referring to (<ref>),
(<ref>), (<ref>) 
we set
for all _h, _h, _h ∈_h^E:

    ^E(_h;  _h, _h)    := ∫_E [ (0  _h  )  (0 _h ) ] ·0 _h   dE 
    
    ^E(_h;  _h, _h)    := 1/2^E(;  , ) - 1/2^E(;  , )
    
    ^E(_h;  _h, _h)    := ∫_E [ (Π_k-1^0,E _h  ) ×(0 _h ) ] ·0 _h   dE

and note that all quantities in the previous formulas are computable. 

Let c_h^E(·;  ·, ·) be one of the discrete trilinear forms listed above.
As usual we define the global approximated trilinear form by adding the local contributions:

    c_h(_h;  _h, _h) := ∑_E ∈Ω_h  c_h^E(_h;  _h, _h),     for all _h, _h, _h ∈_h.


The forms c_h(·;  ·, ·) are immediately extendable to the whole  (simply apply the same definition for any , , ∈). 

Moreover, the trilinear forms c_h(·;  ·, ·) are continuous on , i.e. 
there exists a uniform bounded constant C_h such that

    |c_h(;  , )| ≤C_h  ___    for all , , ∈.

The proof of the continuity for the trilinear forms (·;  ·, ·) and
(·;  ·, ·) can be found in Proposition 3.3 in <cit.>. Analogous techniques can be used to prove  the continuity of the trilinear form (·;  ·, ·).


For what concerns b(·, ·), as noticed in <cit.> we do not need to introduce any approximation of the bilinear form, since it can be exactly computed by the DoFs 𝐃_.






The last step consists in constructing a computable
approximation of the right-hand side (𝐟,  𝐯)
in (<ref>).  We define the approximated load term 𝐟_h as 

    𝐟_h|_E := Π_k^0,E𝐟    for all E ∈Ω_h,

and consider:

    (𝐟_h, 𝐯_h)  = ∑_E ∈Ω_h∫_E 𝐟_h ·𝐯_h   dE = ∑_E ∈Ω_h∫_E Π_k^0,E𝐟·𝐯_h   dE = ∑_E ∈Ω_h∫_E 𝐟·Π_k^0,E𝐯_h   dE  .








 §.§ The discrete problem



Referring to (<ref>), (<ref>),   (<ref>),   (<ref>),  (<ref>)
and (<ref>), the virtual element approximation of the
Navier–Stokes equation in the velocity-pressure formulation is given by:


    {   find (_h, p_h) ∈_h × Q_h, such that
       ν  a_h(_h, _h) + c_h(_h;  _h, _h) + b(_h, p_h) = (_h, _h)        for all _h ∈_h,
         b(_h, q_h) = 0        for all q_h ∈ Q_h,.


with c_h(·;  ·, ·) given by (<ref>),
(<ref>) or (<ref>). 

Whenever the choice (<ref>) is adopted, the pressure
output in (<ref>) approximates the Bernoulli
pressure P in (<ref>) instead of the convective pressure p.

Recalling the kernel inclusion (<ref>),
Problem (<ref>) can be also formulated in the equivalent kernel form

    {   find _h ∈_h, such that
       ν  a_h(_h, _h) + c_h(_h;  _h, _h) = (_h, _h)        for all _h ∈_h..


The well-posedness of the discrete problems can be
stated in the following theorem (cf. <cit.>).



Under the assumption

    (𝐀0)_𝐡    γ_h := C_h  _h_H^-1/α_*^2  ν^2  < 1

Problem (<ref>) has a unique solution (_h, p_h) ∈_h × Q_h such that

    _h_≤_h_H^-1/α_*  ν .



We have the following approximation results (see Theorem 4.6
and Remark 4.2 in <cit.>
for the choices (<ref>) and (<ref>)). 



Under the assumptions  (𝐀0) and (𝐀0)_𝐡, 
let (, p) ∈× Q be the solution of Problem (<ref>) and (_h, p_h) ∈_h × Q_h be the solution of Problem (<ref>). 
Assuming moreover , ∈ [H^s+1(Ω)]^2 and p ∈ H^s(Ω), 0 < s ≤ k, then 

    - _h _≤  h^s ℱ(;  ν, γ, γ_h) +   h^s+2 ℋ(; ν, γ_h)
    
    p  - p_h_Q ≤ C   h^s  |p|_s +  h^s  𝒦(; ν, γ, γ_h) + C   h^s+2  ||_s+1

for a suitable functions 
ℱ, 
ℋ,
𝒦 
independent of h. 


Following the same steps as in <cit.>, Theorem <ref> easily extends also to the choice (<ref>). In such case 
we preliminary observe that if the velocity solution ∈ [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 (<ref>) by the Hölder inequality and Sobolev embedding H^s+1(Ω) ⊂ W^s_4(Ω), we recover

    P_s ≤p_s + 1/2 ·_s + λ_s ≤p_s + 1/2 _W^s_4^2 + λ_s ≤p_s + 1/2 ^2_s+1 + |λ| |Ω|^1/2 .



Now let (_h,  P_h) be the solution of the virtual problem (<ref>) with the trilinear form (<ref>) and (,   P) be the solution of the Navier–Stokes equation (<ref>). Then (<ref>) is substituted by


    P  - P_h_Q ≤ C   h^s  |P|_s +  h^s  𝒦(; ν, γ, γ_h) + C   h^s+2  ||_s+1 .


In such case the following computable approximation p_h of the convective pressure p is available:


    p_h_|E := P_h_|E +
    1/2 Π_k^0, E_h ·Π_k^0, E_h - 
    λ_h

where λ_h is the mean value of the piecewise polynomial function 1/2 Π_k^0, E_h ·Π_k^0, E_h.

The optimal order of accuracy for the convective pressure  can established as follows.
Definitions (<ref>) (taking λ as the mean value of 
1/2 ·)
and (<ref>) easily imply


    p - p_h_Q    ≤P - P_h_Q + 
    1/2(·  - 
      ∑_E ∈Ω_hΠ_k^0, E_h ·Π_k^0, E_h ) - (λ - λ_h) _Q
    
       ≤P - P_h_Q + 
    1/2·  - 
      ∑_E ∈Ω_hΠ_k^0, E_h ·Π_k^0, E_h _Q
    
        = 	P - P_h_Q + 1/2 ( ∑_E ∈Ω_h· - Π_k^0, E_h ·Π_k^0, E_h _Q, E^2 )^1/2
       =:  P - P_h_Q + 1/2 ( ∑_E ∈Ω_hμ_E^2 )^1/2 ,

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 (<ref>) is bounded by (<ref>).

Whereas for the terms μ^E the triangular inequality and  the Hölder inequality entail


    μ_E
       ≤· - Π_k^0, E·Π_k^0, E_E
    + Π_k^0, E·Π_k^0, E - Π_k^0, E_h ·Π_k^0, E_h_E
    
        =
     ( -  Π_k^0, E) ·  ( + Π_k^0, E)_E
    +  
    Π_k^0, E (  -_h) ·Π_k^0, E (  +_h)_E
    
       ≤ -  Π_k^0, E_L^4(E)  + Π_k^0, E_L^4(E) 
    +
    Π_k^0, E (  -_h)  _L^4(E) Π_k^0, E (  +_h) _L^4(E) .


Using the Sobolev embedding H^1(Ω) ⊂ L^4(Ω),
the continuity of the projection Π_k^0, E with respect
to the  L^4-norm and the H^1-norm  (see, for instance, <cit.>), 
and polynomial approximation estimates on star shaped polygons of Lemma <ref>,
from (<ref>) we infer

    μ_E
       ≤ -  Π_k^0, E_, E  + Π_k^0, E_, E  
    +   -_h _L^4(E)   +_h_L^4(E)
       ≤
    C   h^s   ||_s+1, E _, E + 
       -_h _, E   +_h_, E
       ≤
    C   h^s   ||_s+1, E _, E + 
       -_h _, E (_, E  +_h_, E)  .


Combining bound (<ref>) with the Hölder
inequality for sequences, the velocity error estimate
(<ref>) and with the stability estimates
(<ref>) and (<ref>), it follows

    ( ∑_E ∈Ω_hμ_E^2 )^1/2   ≤( 
    ∑_E ∈Ω_h C   h^2s  ||_s+1, E^2  _, E^2 + 
    ∑_E ∈Ω_h  -_h _, E^2     ( _, E^2  +_h_, E^2 ) 
    )^1/2
       ≤
    C   h^s   ||_s+1 _ + 
       -_h _ (_  +_h_)
    
       ≤
    C   h^s   ||_s+1 _H^-1/ν
    +  -_h _  ( _H^-1/ν + 
    _h_H^-1/α_*  ν)
    
       ≤  h^s  ℒ (, ; ν, γ, γ_h) +   h^s+2 ℐ(; ν, γ_h)

for a suitable functions ℒ and ℐ independent of h.

Finally, inserting estimates (<ref>) and (<ref>)
in (<ref>) we obtain the optimal convergence result
for the convective pressure also for choice (<ref>).



















We observe that, due to the divergence-free property, the estimate on the velocity error in Theorem <ref> does not depend on the continuous pressure, whereas the velocity errors of classical methods have a pressure contribution, see <cit.> for more details on this aspect.





A careful investigation of the proof of Theorem <ref> (see Lemma 4.3 in <cit.>) shows that if the exact velocity solution ∈ [_k(Ω)]^2 and the trilinear form (·;  ·, ·) or the trilinear form (·;  ·, ·) are adopted in  (<ref>), the corresponding schemes provide a higher order approximation errors, that are respectively

    - _h _≤ C   h^k+2  () _k+1 + C   h^k+2 _k+1 ,
    
     - _h _≤ C   h^k+2  ( ) ×_k+1 + C   h^k+2 _k+1 .


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.



Another interesting aspect related to method (<ref>)
is the so called “reduced” version. Noting that the discrete
solution satisfies div _h=0, one can immediately
set to zero all 𝐃_𝐕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
<cit.> for more details).





§ VIRTUAL ELEMENTS STOKES COMPLEX AND  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 (<ref>), such that 

    Φ_h = _h  .


In particular this will allow to exploit the kernel formulation (<ref>) in order to define an equivalent VEM approximation for the Navier–Stokes equation in   form (cf. (<ref>)).

Note that a related approach, but limited to a lowest order case and suitable  only for the Stokes problem, was presented in <cit.>.


In order to construct the space of the virtual stream functions Φ_h, we proceed step by step, following the enhanced technique used in Subsection <ref>.

On each element E ∈Ω_h 
we consider the enlarged local virtual space:











    Ψ_h^E := {ϕ∈ H^2(E)   s.t.  ϕ_|∂ E∈_k+1(∂ E)  ,  
    (∇ϕ)_|∂ E∈ [_k(∂ E)]^2  ,   Δ^2  ϕ∈_k-1(E) 
    } .


Then we define the enhanced space of the stream functions

    Φ_h^E := {ϕ∈Ψ_h^E   s.t.   (ϕ - Π^∇,E_k (ϕ),  𝐱^⊥ p_k-1)_E = 0   for all p_k-1∈_k-1 ∖ k-3(E)} .



It is straightforward to see that _k+1(E) ⊆Φ_h^E.

We are now ready to introduce a suitable set of degrees of freedom for the local space of virtual stream functions Φ_h^E.
Given a function ϕ∈Φ_h^E, 
we take the following linear operators 𝐃_Φ, split into five subsets (see Figure <ref>)



  * 𝐃_Φ1:  the values of ϕ at the vertexes of the polygon E,


  * 𝐃_Φ2:  the values of ϕ at the vertexes of the polygon E,


  * 𝐃_Φ3: the values of ϕ at k-2 distinct points of every edge e ∈∂ E,


  * 𝐃_Φ4: the values of ∂ϕ/∂ n at k-1 distinct points of every edge e ∈∂ E,


  * 𝐃_Φ5: the moments of ϕ 


    ∫_E  ϕ·𝐱^⊥  p_k-3  dE     for all p_k-3∈_k-3(E).
 





We note that the linear operators 𝐃_Φ1 and 𝐃_Φ2 are always
needed to enforce the C^1-continuity at the vertices.

Moreover it is immediate to check that for any stream function ϕ∈Φ_h^E
(the same holds for  Ψ_h^E), 
the linear operator evaluations of 𝐃_Φ1,
𝐃_Φ2, 𝐃_Φ3, 𝐃_Φ4 uniquely determine  
ϕ and its gradient on  ∂ E.



We now prove the following result.


The linear operators 𝐃_Φ are a unisolvent set of degrees of freedom for the
virtual space of stream functions Φ_h^E and the dimension of Φ_h^E is

    (Φ_h^E ) =  2   n_E   k + (k-1)(k-2)/2
 

where as usual n_E denotes the number of edges of the polygon E.



We preliminary prove that
the linear operators 𝐃_Φ plus 
the additional moments of ϕ 

    𝐃_Ψ5 :      ∫_E  ϕ·𝐱^⊥ p_k-1  dE     for all p_k-1∈_k-1 ∖ k-3(E)

constitute a set of degrees of freedom for Ψ_h^E. 

An integration by parts and recalling Remark <ref> imply that the linear operators 𝐃_Ψ5 + 𝐃_Φ5 are equivalent to prescribe the moments
∫_E ϕ  q_k-1  dE for all q_k-1∈_k-1(E).
Indeed, Remark <ref> and simple computations give

    ∫_E ϕ  q_k-1  dE    = 
    ∫_E ϕ  ( 𝐱^⊥  p_k-1)   dE = 
    ∫_E ( ϕ) ·𝐱^⊥  p_k-1  dE 
    - ∑_e ∈∂ E∫_eϕ 𝐱^⊥  p_k-1·𝐭_e  ds  ,

where the boundary integral is computable using the DoFs values.

Now the assertion easily follows by a direct application of Proposition 4.1 in
<cit.>.
In particular, from (<ref>) it holds that

    ( Ψ_h^E ) = 2   n_E   k + k(k+1)/2 .


The next step is to prove that the linear operators 𝐃_Φ 
are unisolvent for Φ_h^E.

From (<ref>) it holds

    ( _k-1 ∖ k-3(E)  ) =
    ( _k-1(E) ) -
    ( _k-3(E)) 
    = 2   k - 1  .



Hence, neglecting the independence of the additional (2 k - 1)
conditions in (<ref>), 
the dimension of Φ_h^E is bounded from below by

    ( Φ_h^E ) ≥( Ψ_h^E ) - (2 k - 1) = 2   n_E   k 
    + (k-1)(k-2)/2 = number of DoFs  𝐃_Φ.


Due to (<ref>), the proof is concluded if we
 show that a function ϕ∈Φ_h^E
such that 𝐃_Φ(ϕ) = 0 is identically zero.

In such case, 
ϕ = 0 and  ∇ϕ = 0 on the skeleton ∂ E and this entails 
ϕ = 0 on ∂ E. 

Moreover we note that in this case the 
Π^∇,E_k ( ϕ) = 0; as a matter of fact, by definition (<ref>), we get

    ∫_E  ( ϕ) : 𝐩_k   dE 
    = - ∫_E  ϕ·𝐩_k   dE 
    + ∫_∂ E ϕ·𝐩_k 	 𝐧_E   dS  .

The boundary integral is zero being ϕ = 0 on the skeleton ∂ E. 

For the internal integral, in the light of (<ref>), let us set 
𝐩_k =  q_k-1 + 𝐱^⊥  p_k-3 with
q_k-1∈_k-1(E) ∖.

Then we infer 

    ∫_E  ( ϕ) : 𝐩_k   dE 
       = 
    -  ∫_E  ϕ · q_k-1  dE 
    -  ∫_E  ϕ·𝐱^⊥  p_k-3  dE 
    
       =
    ∑_e ∈∂ E∫_e ϕ ∂ q_k-1/∂ t  dS
    -  ∫_E  ϕ·𝐱^⊥  p_k-3  dE 
    = 0   ,


where the boundary integral is zero since ϕ = 0 on ∂ E,
whereas the second term is zero since   𝐃_Φ5(ϕ)= 0.

In particular we proved that, since Π_k^∇, E(  ϕ) =0, 
recalling (<ref>) also the moments 𝐃_Ψ5 of ϕ are zero. 

Since 𝐃_Φ(ϕ) = 0 by assumption, recalling that 
ϕ∈Φ_h^E ⊂Ψ_h^E and that 
{𝐃_Φ,  𝐃_Ψ5} are a set of
degrees of freedom for Ψ_h^E, it follows ϕ = 0.











An alternative way to define a unisolvent set of DoFs for the
space Φ_h^E is to provide in the place of 𝐃_Φ5 
the following  operators <cit.>


  * 𝐃_Φ5 : the moments of ϕ
against the polynomial of degree up to degree k-3


    ∫_E ϕ  p_k-3  dE     for all p_k-3∈_k-3(E),


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
Φ_h^E accordingly to the local degrees of freedom,
taking into account the boundary conditions:

    Φ_h := {ϕ∈Φ  s.t.  ϕ_| E∈Φ_h^E    for all E∈Ω_h} .

The dimension is
(Φ_h ) =  3   n_V + (2   k - 3) n_e + n_P  (k-1)(k-2)/2 ,
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 <cit.>:

    0   
    H_0^2(Ω)   
    [H_0^1(Ω)]^2   
    L_0^2(Ω)   
    0  ,

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 Φ_h^E as the
space of the stream functions associated to the discrete kernel _h^E.




For any E ∈Ω_h let _h^E and Φ_h^E be the
spaces defined in (<ref>) and (<ref>), respectively.
Then, it holds that

    Φ_h^E = _h^E  .





For any ϕ_h ∈Φ_h^E, we show that the function _h : =  ϕ_h ∈_h^E.

First of all, it is straightforward to check that

    _h = ( ϕ_h ) = 0  .


Concerning the condition on the  skeleton ∂ E, we observe that 
ϕ_h_| ∂ E∈_k+1(∂ E) and (ϕ_h)_| ∂ E∈ [_k(∂ E)]^2 easily imply that 


    _h_| ∂ E  = ( ϕ_h)_| ∂ E∈ [_k(∂ E)]^2  .


Inside the element, by simple calculations and by definition (<ref>), we infer
_h  =
  ( ϕ_h) =  
Δ ^2ϕ_h
∈_k-1(E).
In the light of Remark <ref>, the previous relation is equivalent to

    _h  
     ∈ ( 𝐱^⊥ _k-1(E) )  .

Therefore, there exists p_k-1∈_k-1(E)
such that ( _h - 𝐱^⊥  p_k-1) = 0. 
Since E is simply connected, there exists s such that
_h - 𝐱^⊥  p_k-1 = ∇ s.
Thus we have shown that

    (_h -  ∇  s)  ∈ 𝐱^⊥ _k-1(E) .

Moreover, by definition (<ref>), for all
p_k-1∈_k-1 ∖ k-3(E) it holds

    (_h - Π^∇,E_k  _h,  𝐱^⊥ p_k-1)_E = 
    (ϕ_h - Π^∇,E_k (ϕ_h),  𝐱^⊥ p_k-1)_E = 0  .


At this point is clear that (<ref>), (<ref>), (<ref>),
(<ref>) and definition (<ref>), imply _h =  ϕ_h ∈_h^E for any scalar potential ϕ_h ∈Φ_h^E,
i.e. 
Φ_h^E ⊆_h^E.

The proof now follows by a dimensional argument.
In fact, from (<ref>) and (<ref>) easily follows that 

    ((Φ_h) ) =
     (Φ_h^E ) - 1 =
    2    n_E   k + (k-2)(k-1)/2   - 1 = ( _h^E )  .


Therefore we can conclude that Φ_h^E = _h^E for all E ∈Ω_h.






Given any ϕ_h ∈Φ_h^E, from the degrees of freedom values 𝐃_Φ of ϕ_h,
we are able to compute the DoFs values 𝐃_𝐕 of ϕ_h.
In particular it holds that 

    𝐃_𝐕3( ϕ_h) = 𝐃_Φ5(ϕ_h)
        and    𝐃_𝐕4( ϕ_h) = 0  .


Therefore, for any ϕ_h ∈Φ_h^E, the DoFs 𝐃_Φ allow to compute
the polynomial projections
( ϕ_h), 0 ( ϕ_h) and 0∇( ϕ_h).


As a consequence of Proposition <ref> we have the following
Stokes exact sequence for our discrete VEM spaces and its reduced
version (see also Figures <ref> and <ref>).



Let  Φ_h^E and _h^E  the spaces defined in  (<ref>)
and (<ref>), respectively,
and let _h^E denote the reduced velocity space, see Remark <ref>.
Then, the following sequences are exact 

    Φ_h^E   _h^E _k-1(E)  
    0  ,


    Φ_h^E   _h^E  _0(E)   
    0  .















In terms of degrees of freedom, our lowest order element
(when restricted to triangles!) can be compared with the
Zienkiewicz element <cit.>,
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 <ref> and Figure <ref>) 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 <cit.>. 

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 <ref> and Corollary <ref>.



Let _h and Φ_h be the spaces defined in (<ref>) and
(<ref>), respectively.
Then, it holds that

    Φ_h = _h  .




We note that Proposition <ref> endowed with the boundary
condition in the definitions  (<ref>) and (<ref>)
imply Φ_h ⊆_h.
The proof now follows by a dimensional argument using the Euler formula.




Let Φ_h^E, _h and Q_h be the spaces defined in
(<ref>), (<ref>) and (<ref>), respectively,
and let _h and Q_h denote the reduced
velocity space and the piecewise constant pressures, respectively,
see Remark <ref>.
Then, the following sequences are exact 

    0   Φ_h   _h Q_h  
    0  ,


    0   Φ_h   _h  Q_h   
    0  .

The case without boundary conditions follows analogously.





 §.§ The discrete problem



In the light of Proposition <ref>, referring to
(<ref>) and (<ref>)
we can set the virtual element approximation of the
Navier–Stokes equation in the  formulation:

    {   find ψ_h ∈Φ_h, such that
       ν  a_h( ψ_h,   ϕ_h) + c_h( ψ_h;   ψ_h,   ϕ_h) = (_h,   ϕ_h)     for all ϕ_h ∈Φ_h..


Due to Proposition <ref>, Problem (<ref>)
is equivalent to (<ref>).

We remark that all forms in (<ref>) are exactly
computable by the DoFs 𝐃_Φ. In fact, recalling
Remark <ref>, the polynomials
( ϕ_h), 0 ( ϕ_h) and 0∇( ϕ_h)
are computable on the basis of 𝐃_Φ, so that,
referring to (<ref>), (<ref>) and (<ref>), we infer that

    a_h( ·,   ·)  ,      
    c_h( ·;   ·,   ·)  ,      
    (_h,   ·)

are exactly computable from DoFs 𝐃_Φ, see also Remark <ref>.

The well-posedness of Problem (<ref>) under
the assumption (𝐀0)_𝐡  
immediately follows from Corollary <ref>. 
In fact  the 
operator on the global space Φ_h is an isomorphism into _h. 
Owing to this isomorphism, the existence and uniqueness result of Theorem <ref> carry over to Problem (<ref>), yielding the following theorem.




Under
the assumption (𝐀0)_𝐡  
Problem (<ref>) has a unique solution ψ_h ∈Φ_h such that

    ψ_h_Φ≤_h_H^-1/α_*  ν .



The convergence
of the discrete solution ψ_h of (<ref>) 
to the continuous solution ϕ of (<ref>)
follows immediately from Theorem <ref>, taking =  ϕ and
_h =  ϕ_h.






















Clearly Problem (<ref>) 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 (<ref>) and (<ref>),  starting from the solution ψ_h of Problem (<ref>), we infer the following problem

    {   find p_h ∈ Q_h, such that
        b(_h,   p_h) = 
    - ν  a_h( ψ_h,  _h) 
    - c_h( ψ_h;   ψ_h,  _h) +  (_h,  _h) 
        for all _h ∈_h. .

Since (_h) > (Q_h), the previous system, is actually an  overdetermined system, i.e. there are more equations than unknowns.
Nevertheless the well-posedness  of Problem (<ref>) is guaranteed by Theorem <ref>.
We refer to Section <ref> for a deeper analysis and computational aspects of (<ref>).

We stress that the  virtual formulation (<ref>) exhibits important differences
from the computational point of view compared with the velocity-pressure formulation (<ref>).

First of all the linear system associated to Problem (<ref>) has
2(n_P - 1) less DoFs than Problem  (<ref>), even if considering its equivalent reduced form (see Remark <ref>).

Moreover the first iteration of the Newton method applied to the 
the non-linear virtual stream formulation
(<ref>) results in a linear system which is symmetric and positive definite,
whereas applied to the virtual element method (<ref>) 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.





Simple integration by parts gives
(,   ϕ) = ( ,  ϕ).
By Remark <ref>, the DoFs 𝐃_Φ allow us to compute  the L^2-projection Π^0,E_k-1Φ_h^E →_k-1(E), so that we can consider a  new computable right hand-side

    (( )_h, ϕ_h )  := 
    
    ∑_E ∈Ω_h∫_E Π_k-1^0, E ( )  ϕ_h   dE = 
    ∑_E ∈Ω_h∫_E ( )  Π_k-1^0,Eϕ_h   dE  .

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, we improve the error estimate in (<ref>) by removing the dependence of the error by the load.
More generally, with the choice (<ref>), 
we completely remove the influence in the error stemming  from the irrotational part in the Helmholtz decomposition of the load.

Clearly (<ref>) 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 (<ref>), also compared with a direct C^1 VEM discretization of the stream formulation (<ref>) described in the Appendix, see equation (<ref>). 

For the scheme (<ref>), in all tests we investigate the three possible options for the trilinear form  in (<ref>), (<ref>), (<ref>).


In Test <ref>  we study the convergence of the proposed
virtual element schemes (<ref>) and (<ref>) for the  discretization of the Navier–Stokes equation
in  formulation and stream formulation respectively.

A comparison of (<ref>) (in terms of errors, number of DoFs, condition number of the resulting linear systems) with the equivalent virtual element scheme (<ref>) for the Navier–Stokes equation in  velocity-pressure formulation is also performed. 


In Test <ref> we consider a benchmark problem for the Navier–Stokes equation (<ref>) 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 <ref>, estimate (<ref>) and Remark <ref>). 

In order to compute the VEM errors, we consider the computable error quantities:

    (𝐮, H^1) := ( ∑_E ∈Ω_h∇  -  Π_k-1^0, E (∇ _h) _0,E^2 )^1/2

for the velocity-pressure formulation (<ref>) and 

    (ψ, H^2) := ( ∑_E ∈Ω_h∇ ψ -  Π_k-1^0, E (∇ ψ_h) _0,E^2 )^1/2

for the  and stream formulations (see (<ref>)
and (<ref>), respectively).

For what concerns the pressures we simply compute the standard L^2 error
(p, L^2)  := p - p_h_0.

For the computation of the discrete pressure for the virtual element scheme (<ref>) 
we follow (<ref>) and solve the overdetermined system by means of the least squares method.
We briefly sketch the construction of the least square formula.
Let {_j}_j=1^ dim(_h) be  the canonical basis functions of _h
and let us denote with  r_h the vector with component

    r_h, j :=
    - ν  a_h( ψ_h,  _j) 
    - c_h( ψ_h;   ψ_h,  _j) +  (_h,  _j) 
        for j=1, …,  dim(_h).

i.e. r_h contains the values of the degrees of freedom associated
to the right hand side of (<ref>) with respect to the basis {_j}_j=1^ dim(_h).

Similarly, let {q_i}_i=1^ dim(Q_h) 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 (<ref>) is

    B  B^T  p_h = B  r_h
  
where the matrix B ∈^ dim(Q_h) × dim(_h) is defined by

    B_i,j = b(_j,   q_i) 
        for i=1, …,  dim(Q_h) and j=1, …,  dim(_h).



We note that, if a code for formulation (<ref>) is available, all is needed in order to implement (<ref>) is the construction of the rectangular matrixes that represent, in terms of the degrees of freedom, the local operator : Φ_h^E →_h^E. Once this matrixes are built, one just needs to combine them with the local stiffness matrixes of scheme (<ref>) and finally assemble the global systems for the space Φ_h. 
 

The polynomial degree of accuracy for the numerical tests is k=2.

In the experiments we consider the computational domains Ω_ Q := [0,1]^2 and Ω_ D := {𝐱∈^2  s.t.  |𝐱| ≤ 1 }. The square domain Ω_ Q is partitioned using the following sequences of polygonal meshes:


  * {𝒱_h}_h: sequence of CVT (Centroidal Voronoi tessellation) with h=1/8, 1/16, 1/32, 1/64,

  * {𝒬_h}_h: sequence of distorted quadrilateral meshes with h=1/10, 1/20, 1/40, 1/80.

An example of the adopted meshes is shown in Figure <ref>.  

For the generation of the Voronoi meshes we use the code Polymesher <cit.>.

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. 



For what concerns the disk Ω_ D we consider the sequences of polygonal meshes:


  * {𝒯_h}_h: sequence of triangular meshes with h=1/5, 1/10, 1/20, 1/40,

  * {𝒲_h}_h: sequence of mapped CVT with h= 1/4, 1/8, 1/16, 1/20.


The meshes 𝒲_h are obtained by mapping a CVT on the square 
[-1, 1]^2 on the disk Ω_ D through the map

    Σ [-1, 1]^2 →Ω_ D     with    Σ (x, y) ↦( x √(1 - y^2/2),   y √(1 - x^2/2))  .
 

Figure <ref> displays an example of the adopted meshes.






In this test we solve the Navier–Stokes equation 

on the square domain Ω_ Q with  viscosity ν = 1 and with
the load term 𝐟 chosen such that the analytical velocity-pressure solution 
and the corresponding stream function solution are respectively 


    𝐮(x,y) = 1/2 [   sin(2 π x)^2  sin(2 π y)  cos(2 π y); - sin(2 π y)^2  sin(2 π x)  cos(2 π x) ]     
    p(x,y) = π^2  sin(2 π x)  cos(2 π y)  ,
    
    ψ(x,y) = 
    1/8 π sin(2 π x)^2  sin(2 π y)^2  .



We test the virtual element scheme (<ref>). 

In Figures <ref>
and <ref> we show the results obtained with  the sequences of Voronoi meshes 𝒱_h and  quadrilateral meshes 𝒬_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 (·; ·, ·) we consider the definition (<ref>)).  





We notice that the theoretical predictions of Section <ref> 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 <ref> (left) and <ref> (left) we also depict the error for the direct C^1 discretization of the stream formulation (<ref>), that follows a similar behaviour to (<ref>). Note that we do not compute a pressure error for scheme (<ref>) since the computation of a discrete pressure is a more complex issue in this case, see Remark <ref>.

Finally we test the corresponding virtual element method (<ref>) with the same sequences of polygonal meshes 𝒱_h, 𝒬_h. 
Table <ref> shows the results obtained respectively  with VEM (<ref>) and (<ref>) obtained considering the trilinear form (·; ·, ·).
The results are analogous also for the other two proposed trilinear forms (not shown).



In Table <ref> we compare the number of DoFs and the condition number of the resulting linear systems (stemming from the fist iteration of the Newton method) 
for both formulations (<ref>)  and (<ref>).

As observed in Section <ref>, the scheme (<ref>) has the advantage of having (2  n_p - 2) less of unknowns, even when considering the reduced version (see Remark <ref>) for formulation (<ref>).

The drawback is that the condition number of the system resulting from the velocity-pressure scheme (<ref>) behaves as h^-2, while the asymptotic rate of the condition number of the linear system resulting from the scheme (<ref>) formulation is h^-4.







In this test, inspired by <cit.>,
we consider a problem  for the Navier–Stokes
with the property of having the velocity and
stream solution in the corresponding  discrete spaces.  

In particular we consider the Navier–Stokes
(<ref>) on the disk Ω_ D  
with  viscosity ν = 1 and with
the load 
𝐟 chosen such that the analytical velocity-pressure solution 
and the corresponding stream function solution are respectively 

    𝐮(x,y) = [ x^2 + y^2;  -2   x y ]     
    p(x,y) = x^3 y^3 - 1/16 ,
    
    ψ(x, y) = x^2 y + 1/3 y^3  .
  
We notice that the stream function solution ψ belongs to the discrete spaces Φ_h and Φ_h (cf. (<ref>) and (<ref>) respectively) 
and the corresponding velocity =  ψ belongs to [_k(Ω)]^2 ⊂_h.
In the light of Remark <ref>, the methods (<ref>) obtained considering the trilinear forms (·; ·, ·) or (·; ·, ·) provide a higher order approximation error, that is h^k+2 (instead of h^k as a standard inf-sup stable FEM would do). Note that instead an exactly divergence-free FEM would acquire machine precision velocity error on this test, but this comes at the price of the additional complications of this methods.

Figures <ref> and <ref> plot the results obtained with  the sequence of triangular meshes 𝒯_h and  mapped Voronoi meshes 𝒲_h considering the scheme (<ref>) (with all possible choices of the trilinear forms)  and the scheme (<ref>).

The error for the C^1 discretization of the stream formulation (<ref>) provides the same higher order of convergence h^k+2. The error for scheme (<ref>) in the case of the finest Voronoi mesh is not depicted since the Newton iterations did not converge in that case, probably an indication on the reduced robustness of such scheme when compared to (<ref>).












§ APPENDIX: VIRTUAL ELEMENTS STREAM FORMULATION




In the present section we  briefly present a stream-function formulation approach, that is based on a direct approximation with C^1 Virtual Elements <cit.> of the continuous stream function formulation of the problem (<ref>). This approach is not associated to a discrete Stokes complex and therefore lacks an equivalent discrete velocity-pressure formulation.

Recalling (<ref>), we consider on each element
E ∈Ω_h the finite dimensional local virtual space

    Φ_h^E := {ϕ∈Ψ_h^E    s.t.  ( ϕ - Π^∇^2,E_k  ϕ,  p_k-1)_E = 0   for all p_k-1∈_k-1 ∖ k-3(E)} .



We here mention only the main properties of the virtual space
Φ_h^E and refer to <cit.>
for a deeper description:



  * dimension:
it holds _k+1(E) ⊂Φ_h^E and 
( Φ_h^E ) = 2n_E k + (k-1)(k-2)/2,
that is the same dimension of Φ_h^E (cf. (<ref>));


  * degrees of freedom:
the linear operators 𝐃_Φ:
𝐃_Φ1, 𝐃_Φ2, 𝐃_Φ3,
𝐃_Φ4, 𝐃_Φ5
(see Remark <ref>) constitute a set of
DoFs for Φ_h^E;

  * projections:
the DoFs 𝐃_Φ allow us to compute exactly 
(c.f. (<ref>) and (<ref>))

    Π_k+1^∇^2, EΦ_h^E →_k+1(E)  ,     Π_k-1^0, EΔ (Φ_h^E) →_k-1(E)  ,     Π_k-1^0, EΦ_h^E →_k-1(E)  ,
    
    Π_k-1^0, E (Φ_h^E) → [_k-1(E)]^2  ,     Π_k-1^0, E (Φ_h^E) → [_k-1(E)]^2  ,

in the sense that, given any ϕ_h ∈Φ_h^E,
we are able to compute the polynomials
Π_k+1^∇^2, Eϕ_h, 
Π_k-1^0, EΔϕ_h,
Π_k-1^0, Eϕ_h,
Π_k-1^0, Eϕ_h,
and 
Π_k-1^0, Eϕ_h
only using, as unique information, the degree of
freedom values 𝐃_Φ of ϕ_h.


The global virtual element space is obtained as usual by combining the local spaces Φ_h^E
accordingly to the local degrees of freedom 


    Φ_h := {ϕ∈ H^2_0(Ω)   s.t  ϕ_|E∈Φ_h^E    for all E ∈Ω_h} .
 










Now we briefly describe the
construction of a discrete computable version of the bilinear
form a(·, ·)  given in (<ref>)
and trilinear form c(·; ·, ·) (cf. (<ref>)).


First, we decompose into local contributions
a(·, ·) and c(·; ·, ·)
by considering:


    a (ψ,  ϕ) =: 
    ∑_E ∈Ω_ha^E (ψ,  ϕ) ,     c(ζ;  ψ, ϕ) =: ∑_E ∈Ω_hc^E(ζ;  ψ, ϕ)     for all ζ, ψ, ϕ∈Φ.


Following a standard procedure in the VEM framework,
we define a computable discrete local bilinear form



a_h^E(·, ·) 
Φ_h^E ×Φ_h^E →

approximating the continuous form a_h^E(·, ·), and defined by

    a_h^E(ψ,  ϕ) := 
    (Π_k-1^0, EΔψ,  Π_k-1^0, EΔϕ)_E + 
    𝒮^E ( (I - ) ψ,  (I -) ϕ)

for all ψ, ϕ∈Φ_h^E, where the (symmetric)
stabilizing bilinear form
𝒮^E Φ_h^E ×Φ_h^E →,
satisfies the stability condition

    α_* ϕ^2_Φ,E≤𝒮^E(ϕ,  ϕ) ≤α_*  ϕ^2_Φ,E    for all ϕ∈Φ_h^E such that ϕ = 0

with α_* and α^*  positive  constants independent of the element E.

We notice that the bilinear form  a_h^E(·,  ·)
is not coercive, since its kernel consists of all harmonic polynomial functions.

For what concerns the approximation of the local trilinear form c^E(·;  ·, ·), we simply set

    c_h^E(ζ_h;  ψ_h,  ϕ_h) := 
    ∫_E [ (Π^0, E_k-1 Δζ_h  )  (Π_k^0, E ψ_h ) ] ·(Π^0, E_kψ_h )   dE

for all ζ_h, ψ_h, ϕ_h ∈Φ_h, which is computable by the DoFs 𝐃_Φ.

We define the global approximated forms a_h(·, ·) 
and c_h(·;  ·, ·)
 by simply summing the local contributions:


    a_h(ψ_h,  ϕ_h) := 
    ∑_E ∈Ω_ha_h^E(ψ_h,  ϕ_h)      c_h(ζ_h;  ψ_h,  ϕ_h) := 
    ∑_E ∈Ω_hc_h^E(ζ_h;  ψ_h,  ϕ_h) ,

for all ζ_h, ψ_h, ϕ_h ∈Φ_h. 
Finally the right hand side is computed as explored in (<ref>).

Referring to (<ref>),  (<ref>) and (<ref>), the virtual element stream formulation of the Navier–Stokes equation is given by


    {         
       ν a_h(ψ_h,  ϕ_h) 
      +   c_h(ψ_h;  ψ_h,  ϕ_h)  = (( )_h,  ϕ_h) 
              for all  ϕ_h ∈Φ_h..
 



Note that the C^1 discretization of Navier–Stokes (<ref>)
does not allow to reconstruct an approximation of the continuous pressure
p analogous to that given in (<ref>).
Indeed the core idea behind (<ref>) is to exploit the
Stokes complex associated to Φ_h and _h (cf. (<ref>)
and (<ref>), respectively).






§ ACKNOWLEDGEMENTS

The authors L. Beirão da Veiga and G. Vacca were
partially supported by the European Research Council through
the H2020 Consolidator Grant (grant no. 681162) CAVE,
Challenges and Advancements in Virtual Elements.
This support is gratefully acknowledged.
The author D. Mora was partially supported by
CONICYT-Chile through project FONDECYT 1180913 and
project AFB170001 of the PIA Program: Concurso Apoyo
a Centros Científicos y Tecnológicos de Excelencia con
Financiamiento Basal.

tocsection
plain ]