An H^1-conforming Virtual Element Method for Darcy equations and Brinkman equations
Giuseppe Vacca































remark
remarkRemark[section]

remark
testTest[section]

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





















sistema
{[ .

remark































An H^1-conforming  Virtual Element Methods for  Darcy equations and  Brinkman equations
    Giuseppe Vacca  
  Dipartimento di Matematica e Applicazioni,  Università degli Studi di Milano Bicocca, Via Roberto Cozzi 55 - 20125 Milano, Italy; E-mail: giuseppe.vacca@unimib.it. 
  
    
==============================================================================================================================================================================================

 The focus of the present paper is on developing a Virtual Element Method for Darcy and Brinkman equations.

In <cit.>  we presented a family of Virtual Elements for Stokes equations and we defined a new Virtual Element space of velocities  such that the associated discrete kernel is pointwise divergence-free.  We use a slightly different Virtual Element space having two fundamental properties: the L^2-projection onto ℙ_k is exactly computable on the basis of the degrees of freedom, and the
associated discrete kernel is still pointwise divergence-free. 

The resulting numerical scheme for the Darcy equation has optimal order of convergence and H^1 conforming velocity solution. 

We can apply the same approach to develop a robust virtual element method for the Brinkman equation that is stable for both the Stokes and Darcy limit case.

We provide a rigorous error analysis of the method and several numerical tests.
 









An H^1-conforming  Virtual Element Methods for  Darcy equations and  Brinkman equations
    Giuseppe Vacca  
  Dipartimento di Matematica e Applicazioni,  Università degli Studi di Milano Bicocca, Via Roberto Cozzi 55 - 20125 Milano, Italy; E-mail: giuseppe.vacca@unimib.it. 
  
    
==============================================================================================================================================================================================






§ INTRODUCTION






The Virtual Element Methods (in short, VEM or VEMs) is a recent technique for solving PDEs. VEMs were recently introduced in <cit.> as a generalization of the finite element method on polyhedral or polygonal meshes. 

In the numerical analysis and engineering literature there has been a recent growth of interest 
in developing numerical methods that 
can make use of general polygonal and polyhedral meshes, as opposed to more standard triangular/quadrilateral (tetrahedral/hexahedral) 
grids. Indeed, making use of polygonal meshes brings forth a range of advantages, including for instance automatic hanging node treatment, 

more efficient approximation of geometric data features, better domain meshing 
capabilities, more efficient and easier 
adaptivity, more robustness to mesh deformation, and others. This interest in the literature is also reflected in commercial codes, 
such as CD-Adapco, that have recently included polytopal meshes.

We refer to the recent papers and monographs 
<cit.> 
as a brief representative sample of the increasing list of technologies that make use of polygonal/polyhedral meshes. We mention 
here in particular the polygonal finite elements, that generalize finite elements to polygons/polyhedrons by making use of 
generalized non-polynomial shape functions, and the mimetic discretisation schemes <cit.>, that combine ideas from the finite difference and 
finite element methods.

 


The principal idea behind VEM is to use approximated discrete bilinear forms that require only integration of polynomials on the (polytopal) 
element in order to be computed. The resulting discrete solution is conforming and the accuracy granted by such discrete bilinear forms 
turns out to be sufficient to achieve the correct order of convergence. 

Following this approach, VEM is able to make use of very general polygonal/polyhedral meshes without the need to integrate 
complex non-polynomial functions on the elements and without loss of accuracy. Moreover, 
VEM is not restricted to low order converge and can be easily applied to three dimensions and use non convex (even non simply connected) elements.

The Virtual Element Method has been developed successfully for a large range of problems, see for instance <cit.>. 
A helpful paper for the implementation of the method is <cit.>.


The focus of this paper is on developing a new Virtual Element Method for the Darcy equation that is suitable for a robust extension to the (more complex) Brinkman problem. 

For such a problem, other VEM numerical schemes  have been proposed, see for example <cit.>. 




In <cit.> the authors developed a new Virtual Element Method for  Stokes problems by exploiting  the flexibility of the Virtual Element construction in a new way. In particular, they define a new Virtual Element space of velocities carefully designed to solve the Stokes problem. In connection with a suitable pressure space, the new Virtual Element space leads to an exactly divergence-free discrete velocity, a favorable property when more complex problems, such as the Navier-Stokes problem, are considered. We highlight that this feature is not shared by the method defined in <cit.> or by most of the standard mixed Finite Element methods, where the divergence-free constraint is imposed only in a weak (relaxed) sense. 


In the present contribution we develop  the Virtual Element Method for Darcy equations by introducing a slightly different virtual space for the velocities such that the local L^2 orthogonal  projection onto the space of polynomials of degree less or equal than k (where k is the polynomial degree of accuracy of the method) can be computed using the local degrees of freedom. The resulting Virtual Elements family inherits the advantages on the scheme proposed in <cit.>, in particular it yields an exactly divergence-free discrete kernel. Thus we obtain a stable Darcy element that is also uniformly stable for the Stokes problem. A sample of uniformly stable methods for Darcy-Stokes model is for instance  <cit.>.

The last part of the paper deals with the analysis of a new mixed finite element method for  Brinkman equations that stems from the above scheme for the Darcy problem. Mathematically, the Brinkman problem resembles both the Stokes problem for fluid flow and the Darcy problem for flow in porous media (see <cit.>).
Constructing finite element methods to solve the Brinkman equation that are robust for both (Stokes and Darcy) limits is challenging. We will see how the above Virtual Element approach offers a natural and straightforward
framework for constructing stable numerical algorithms for the Brinkman
equations.

We remark that the proposed scheme belongs to the class of the pressure-robust method, i.e. delivers a velocity error independent of the continuous pressure.









The paper is organized as follows. 

In Section <ref> we introduce the model continuous Darcy problem. 

In Section <ref> we present its VEM discretisation. 

In Section <ref> we detail the theoretical features and the convergence analysis of the problem. 


In Section <ref> we  develop a stable numerical methods for  Brinkman equations.

In Section <ref> we show the numerical tests.

Finally in the Appendix we present the theoretical analysis of the extension to the Darcy equation of the scheme of  <cit.>. Even though this latter method is not recommended for the Darcy problem, the numerical  experiments showed an unexpected optimal convergence rate for the pressure. We theoretically prove this behaviour, developing an inverse inequality for the VEM space, which is interesting on its own. 








§ THE CONTINUOUS PROBLEM



We consider the classical Darcy equation that describes the flow of a fluid through a porous medium.

Let Ω⊆^2 be a bounded polygon then the Darcy equation in mixed form is


    {   find (𝐮, p) such that
       ^-1𝐮  + ∇ p = 0          in Ω,
        div 𝐮 = f           in Ω,
       𝐮·𝐧 = 0            on ∂Ω,.


where 𝐮 and p are respectively the velocity and the pressure fields, f ∈ L^2(Ω) is the source term and  is a uniformly symmetric, positive definite tensor that represents the permeability of the medium. 







From (<ref>), since we have assumed no flux boundary
conditions all over ∂Ω, the external force f has zero mean value on Ω. 

We consider the spaces

    𝐕:= {𝐮∈ H( div, Ω),   s.t  𝐮·𝐧 = 0   on ∂Ω},      Q:= L^2_0(Ω) = { q ∈ L^2(Ω)   s.t.  ∫_Ω q   dΩ = 0 }

equipped with the natural norms

    𝐯_𝐕^2 := 𝐯_[L^2(Ω) ]^2^2   +   div 𝐯_L^2(Ω)^2    ,     q_Q :=  q_L^2(Ω),


and the bilinear forms a(·, ·) 𝐕×𝐕→ and b(·, ·) 𝐕× Q → defined by:

    a (𝐮,  𝐯) := ∫_Ω^-1 𝐮·𝐯  dΩ,     for all 𝐮,  𝐯∈𝐕



    b(𝐯, q) :=  ∫_Ω div 𝐯  q   dΩ    for all 𝐯∈𝐕, q ∈ Q.


Then the variational formulation of Problem (<ref>) is

    {   find (𝐮, p) ∈𝐕× Q, such that
        a(𝐮, 𝐯) + b(𝐯, p) = 0       for all 𝐯∈𝐕,
         b(𝐮, q) = (f, q)        for all q ∈ Q,.

where

    (f, q) := ∫_Ω f  q   dΩ    for all q ∈ Q.


Let us introduce the kernel

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

then it is straightforward to see that

    𝐯_𝐕 := 𝐯_[L^2(Ω) ]^2    for all 𝐯∈𝐙.

It is well known that (see for instance <cit.>):



  * a(·, ·) and b(·, ·) are continuous, i.e.

    |a(𝐮, 𝐯)| ≤a𝐮_𝐕𝐯_𝐕    for all 𝐮, 𝐯∈𝐕,
 

    |b(𝐯, q)| ≤b𝐯_𝐕q_Q     for all 𝐯∈𝐕 and q ∈ Q;


  * a(·, ·) is coercive on the kernel 𝐙, i.e. there exists a positive constant α depending on  such that 

    a(𝐯, 𝐯) ≥α𝐯^2_𝐕    for all 𝐯∈𝐙;


  * b(·,·) satisfies the inf-sup condition, i.e.   

    ∃ β >0   such that  sup_𝐯∈𝐕 𝐯≠0b(𝐮, q)/𝐯_𝐕≥βq_Q     for all q ∈ Q.


Therefore, Problem (<ref>) has a unique solution (𝐮, p) ∈𝐕× Q such that

    𝐮_𝐕 + p_Q ≤ C  f_L^2(Ω)

with the constant C depending only on Ω and .







§ VIRTUAL FORMULATION FOR DARCY EQUATIONS







 §.§ Decomposition and the original virtual element spaces


 
We outline the Virtual Element discretization of Problem (<ref>). 
Here and in the rest of the paper the symbol C will indicate a generic positive constant independent of the mesh size that may change at 
each occurrence. Moreover, given any subset ω in ℝ^2 and k ∈ℕ, we will denote by _k(ω) the 
polynomials of total degree at most k defined on ω, with the extended notation _-1(ω)=∅.
Let 𝒯_h_h be a sequence of decompositions of Ω into general polygonal elements K with

    h_K :=  diameter(K) ,   
    h := sup_K ∈𝒯_h h_K .


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


  * (𝐀1) K is star-shaped with respect to a ball of radius ≥ γ  h_K, 

  * (𝐀2) the distance between any two vertexes of K is ≥ c   h_K, 

where γ and c are positive constants. We remark that the hypotheses above, though not too restrictive in many practical cases, 
can be further relaxed, as noted in  <cit.>. 
From now on we assume that  is  piecewise constant with respect to 𝒯_h on Ω. 



Using standard VEM notation, for k ∈, let us define the spaces 



  * _k(K) the set of polynomials on K of degree ≤ k,

  * _k(K) := {v ∈ C^0(∂ K)   s.t   v_|e∈_k(e)   ∀ e ⊂∂ K},

  * 𝒢_k(K):= ∇(_k+1(K)) ⊆ [_k(K)]^2,

  * 𝒢_k(K)^⊥ := 𝐱^⊥[_k-1(K)] ⊆ [_k(K)]^2 with 𝐱^⊥:= (x_2, -x_1).



In <cit.> the authors have introduced a new family of Virtual Elements for the Stokes problem
on polygonal meshes. In particular, by a proper choice of the Virtual space of velocities, the virtual local spaces are associated to a Stokes-like variational problem on each element. The main ideas of the method are


  * the Virtual space contains the space of all the polynomials of the prescribed order plus suitable non polynomial functions,

  * the degrees of freedom are carefully chosen so that the H^1 semi-norm projection onto the space of polynomials can be exactly computed,

  * the choice of the Virtual space of velocities and the associated degrees of freedom guarantee that the final discrete velocity is pointwise divergence-free and more generally the discrete kernel is contained in the continuous one.

In this section we briefly  recall from <cit.> the notations, the main properties of the Virtual spaces and some details of the construction of the  H^1 semi-norm projection.
Let k ≥ 2 the polynomial degree of accuracy of the method, then we define  on each element K ∈𝒯_h the finite dimensional local virtual space

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


where all the operators and equations above are to be interpreted in the distributional sense. 

It is easy to check that
[_k(K)]^2 ⊆𝐖_h^K,  
and that (see <cit.> for the proof) the dimension of 𝐖_h^K is

    ( 𝐖_h^K )    = ([_k(∂ K)]^2) + (𝒢_k-2(K)^⊥) + ( (_k-1(K)) - 1) 
       = 2n_K k + (k-1)(k-2)/2  + (k+1)k/2 - 1.

The corresponding degrees of freedom are chosen prescribing, given a function 𝐯∈𝐖_h^K, the following linear operators 𝐃_𝐕, split into four subsets (see Figure <ref>):




  * 𝐃_𝐕1:  the values of 𝐯 at the vertices of the polygon K,


  * 𝐃_𝐕2: the values of 𝐯 at k-1 distinct points of every edge e ∈∂ K (for example we can take the k-1 internal points of the (k+1)-Gauss-Lobatto quadrature rule  in e, as suggested in <cit.>),


  * 𝐃_𝐕3: the moments of 𝐯

    ∫_K 𝐯·𝐠_k-2^⊥  dK     for all 𝐠_k-2^⊥∈𝒢_k-2(K)^⊥,



  * 𝐃_𝐕4: the moments up to order k-1 and greater than zero of div 𝐯 in K, i.e.


    ∫_K ( div 𝐯)   q_k-1  dK     for all q_k-1∈_k-1(K) /.
 

 


For all K ∈𝒯_h, we introduce the H^1 semi-norm projection Π_k^∇,K𝐖_h^K → [_k(K)]^2, defined by 

    {   ∫_K ∇ 𝐪_k : ∇ (𝐯_h -  Π_k^∇,K𝐯_h)   d K = 0     for all 𝐪_k ∈ [_k(K)]^2,
       Π_0^0,K(𝐯_h -  Π_k^∇,K𝐯_h) = 0  ,
    .
 

where Π_0^0,K is the L^2-projection operator onto the constant functions defined on K. It is immediate to check that the energy projection is well defined and it clearly holds Π_k^∇,K𝐪_k = 𝐪_k for all 𝐪_k ∈_k(K).
Moreover  the operator Π_k^∇,K is computable in terms of the degrees of freedom 𝐃_𝐕 (see equations (27)-(29) in <cit.> and the subsequent discussion).







 §.§ The modified virtual space and the projection Π^0,K_k


Let n a positive integer, then for all K ∈𝒯_h, the L^2-projection Π^0,K_n𝐖^K_h → [_n(K)]^2 is defined by

    ∫_K  𝐪_n· (𝐯_h - Π^0,K_n𝐯_h)    d K = 0     for all 𝐪_n ∈  [_n(K)]^2.

It is possible to check (see Section 3.3 of <cit.> for the proof) that the degrees of freedom 𝐃_𝐕 allow us to compute exactly the L^2-projection Π^0,K_k-2.  On the other hand we can observe that we can not  compute exactly  from the DoFs  the L^2-projection onto the space of polynomials of degree ≤ k.

The goal of the present section is to introduce, taking the inspiration from <cit.>, a new virtual space 𝐕_h^K to be used in place of 𝐖_h^K in such a way that


  * the DoFs 𝐃_𝐕 can still be used for 𝐕_h^K,

  * [_k(K)]^2 ⊆𝐕_h^K,

  * the projection Π_k^0,K𝐕_h^K → [_k(K)]^2 can be exactly computable by the DoFs 𝐃_𝐕.


To construct 𝐕_h^K we proceed as follows: first of all we define an augmented virtual local space 𝐔_h^K by taking

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

Now we define the enhanced Virtual Element space 𝐕_h^K as the restriction  of 𝐔_h^K given by

    𝐕_h^K := {𝐯∈𝐔_h^K   s.t.  (𝐯 - Π^∇,K_k 𝐯,  𝐠_k^⊥)_[L^2(K)]^2 = 0   for all 𝐠_k^⊥∈𝒢_k(K)^⊥/𝒢_k-2(K)^⊥} ,

where the symbol 𝒢_k(K)^⊥/𝒢_k-2(K)^⊥ denotes the polynomials in 𝒢_k(K)^⊥ that are L^2-orthogonal to all polynomials of 𝒢_k-2(K)^⊥.

We proceed by investigating the dimension and by choosing suitable DoFs of the virtual space 𝐕_h. First of all we recall from  <cit.> the following facts 

    ([_k(∂ K)]^2) = 2n_K k,   (_k-1(K)) = k(k+1)/2,   (𝒢_k(K)^⊥) = k(k+1)/2


where n_K is the number of edges of the polygon K. 




The dimension of 𝐔_h^K is

    ( 𝐔_h^K ) 
       = 2n_K k + k(k+1)/2  + (k+1)k/2 - 1.

Moreover as DoFs for 𝐔_h^K we can take the linear operators 𝐃_𝐕 and plus the moments 

    𝐃_𝐔∫_K 𝐯·𝐠_k^⊥  dK     for all 𝐠_k^⊥∈𝒢_k(K)^⊥/𝒢_k-2(K)^⊥.
 



The proof is virtually identical to that given in <cit.> for 𝐖_h^K and it is based (see for instance <cit.>) on the fact that given



  * a polynomial function 𝐠_b ∈ [_k(∂ K)]^2,


  * a polynomial function 𝐡∈𝒢_k(K)^⊥,


  * a polynomial function g ∈_k-1(K)  satisfying the compatibility condition

    ∫_K g   dΩ = ∫_∂ K𝐠_b ·𝐧  ds,



there exists a unique pair (𝐯, s)∈𝐔_h^K×  L^2(K) / such that


    𝐯_|∂ K = 𝐠_b,    div 𝐯 = g,    
     -  Δ𝐯  -  ∇ s = 𝐡.


Moreover, since from <cit.>,  rot𝒢_k(K)^⊥→_k-1(K)

is an isomorphism, 
we can conclude that  the map that associates a given compatible data set   
(𝐠_b,  𝐡,   g) to the velocity field 𝐯 that solves (<ref>) is an injective map. Then

    ( 𝐔_h^K ) = ([_k(∂ K)]^2) + (𝒢_k(K)^⊥) + ( (_k-1(K)) - 1)

and the thesis follows from (<ref>).





The dimension of 𝐕_h^K is equal to that of 𝐖_h^K that is, as in (<ref>)

    ( 𝐕_h^K )    =  2n_K k + (k-1)(k-2)/2  + (k+1)k/2 - 1.

As DoFs in 𝐕_h^K we can take 𝐃_𝐕.



From (<ref>) it is straightforward to check that

    (𝒢_k(K)^⊥/ 𝒢_k-2(K)^⊥) = (𝒢_k(K)^⊥) - (𝒢_k-2(K)^⊥) = 2k -1.

Hence, neglecting the independence of the additional 2k - 1 conditions in (<ref>), it holds that

    ( 𝐕_h^K ) ≥( 𝐔_h^K ) - (2k -1) =  2n_K k + (k-1)(k-2)/2  + (k+1)k/2 - 1 = ( 𝐖_h^K ).

We now observe that a function 𝐯∈𝐕_h^K such that 𝐃_𝐕(𝐯) = 0 is identically zero. Indeed,  from (<ref>), it is immediate to check that in this case the Π_k^∇, K 𝐯 would be zero, implying that all its
moment are zero, in particular, since 𝐯∈𝐕_h^K, all the moments 𝐃_𝐔 of 𝐯 are also zero. Now, from Lemma <ref>, we have that 𝐯 is zero.
Therefore, from (<ref>), we obtain that the dimension of 𝐕_h^K is actually the same of  𝐖_h^K, and that the DoFs 𝐃_𝐕 are unisolvent for 𝐕_h^K.





The degrees of freedom 𝐃_𝐕 allow us to compute exactly the L^2-projection Π^0,K_k𝐕_h → [_k(K)]^2, i.e. the moments

    ∫_K 𝐯·𝐪_k  d K

for all 𝐯∈𝐕_h and for all 𝐪_k∈ [_k(K)]^2. 


Let us set

    𝐪_k = ∇ q_k+1 + 𝐠_k-2^⊥ + 𝐠_k^⊥.

with q_k+1∈_k+1(K) /, 𝐠_k-2^⊥∈𝒢_k-2^⊥(K) and 𝐠_k^⊥∈𝒢_k^⊥(K)/𝒢_k-2^⊥(K). Therefore using the Green formula and since 𝐯∈𝐕_h, we get 

    ∫_K 𝐯·𝐪_k  d K     = ∫_K 𝐯· (∇ q_k+1 +  𝐠_k-2^⊥ + 𝐠_k^⊥)   d K 
        = - ∫_K  div 𝐯  q_k+1  d K  +
    ∫_K 𝐯·𝐠_k-2^⊥  d K +
    ∫_K Π_k^∇, K 𝐯·𝐠_k^⊥  d K  
    + ∫_∂ K q_k+1 𝐯·𝐧  d s.

Now, since div 𝐯 is a polynomial of degree less or equal than k-1 we can reconstruct its value from 𝐃_𝐕4 and compute exactly the first term. The second term is computable from 𝐃_𝐕3. The third term is computable from all the 𝐃_𝐕 using the projection Π_k^∇, K 𝐯. Finally from 𝐃_𝐕1 and 𝐃_𝐕2 we can reconstruct 𝐯 on the boundary and so compute exactly the boundary term. 

For what concerns the pressures we take the standard finite dimensional space

    Q_h^K := _k-1(K)

having dimension

    (Q_h^K) = (_k-1(K))  = (k+1)k/2.

The corresponding degrees of freedom are chosen defining for each q∈ Q_h^K the following linear operators 𝐃_𝐐:




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

    ∫_K q   p_k-1  dK     for all p_k-1∈_k-1(K).




Finally we define the global virtual element spaces as


    𝐕_h := {𝐯∈ [H^1(Ω)]^2    s.t  𝐯·𝐧 = 0   on ∂Ω  and  𝐯_|K∈𝐕_h^K    for all K ∈𝒯_h}
 
and

    Q_h := { q ∈ L_0^2(Ω)   s.t.   q_|K∈  Q_h^K   for all K ∈𝒯_h},

with the obvious associated sets of global degrees of freedom. A simple computation shows that:


    (𝐕_h) = n_P ( (k+1)k/2 -1  +  (k-1)(k-2)/2) 
    + 2(n_V + (k-1) n_E) + (n_V, B + (k-1) n_E, B)


and


    (Q_h) = n_P (k+1)k/2 - 1 ,


where n_P is the number of elements, n_E, n_V (resp., n_E, B, n_V, B) is the number of  internal edges and vertexes (resp., boundary edges and vertexes) in 𝒯_h.
As observed in <cit.>, we remark that


    div 𝐕_h⊆ Q_h .




By definition (<ref>) it is clear that our discrete velocities field is H^1-conforming, in particular we obtain continuous velocities, whereas the natural discretization is only H( div)- conforming. This property, in combination with (<ref>), will make our method suitable for a (robust) extension to the Brinkman problem.





















 §.§ The discrete bilinear forms



The next step in the construction of our method is to define on the virtual spaces 𝐕_h and Q_h a discrete version of the bilinear forms a(·, ·) and b(·, ·) given in (<ref>) and (<ref>).

For simplicity we assume that the tensor  is piecewise constant with respect to the decomposition 𝒯_h, i.e.  is constant on each polygon K ∈𝒯_h.

First of all we decompose into local contributions the bilinear forms a(·,·) and b(·, ·), the norms ·_𝐕 and ·_Q by defining 


    a (𝐮,  𝐯) =: ∑_K ∈𝒯_h a^K (𝐮,  𝐯)     for all 𝐮,  𝐯∈𝐕



    b (𝐯,  q) =: ∑_K ∈𝒯_h b^K (𝐯,  q)     for all 𝐯∈𝐕 and q ∈ Q,


and 


    𝐯_𝐕 =: (∑_K ∈𝒯_h𝐯^2_𝐕, K)^1/2  for all 𝐯∈𝐕,    q_Q =: (∑_K ∈𝒯_hq^2_Q, K)^1/2  for all q ∈ Q.


We now define discrete versions of the bilinear form a(·, ·) (cf. (<ref>)), and of the bilinear form b(·, ·) (cf. (<ref>)). For what concerns b(·, ·), we simply  set 


    b(𝐯, q) = ∑_K ∈𝒯_h b^K(𝐯, q) = ∑_K ∈𝒯_h∫_K  div 𝐯  q   dK     for all 𝐯∈𝐕_h, q ∈ Q_h,


i.e. as noticed in <cit.> we do not introduce any approximation of the bilinear form. We notice that (<ref>)  
is computable from the degrees of freedom 𝐃_𝐕1, 𝐃_𝐕2 and 𝐃_𝐕4, since q is polynomial in each element K ∈𝒯_h. 
On the other hand, the bilinear form a(·, ·) needs to be dealt with in a more careful way.
First of all, by Proposition <ref>, we observe that  for all 𝐪_k ∈ [_k(K)]^2 and for all 𝐯∈𝐕_h^K, the quantity  


    a^K (𝐪_k,  𝐯) = ∫_K^-1 𝐪_k ·𝐯  dK.

is exactly computable by the DoFs.

However, for an arbitrary pair (𝐮,𝐯 )∈𝐕_h^K ×𝐕_h^K, the quantity a_h^K(𝐰, 𝐯) is clearly not computable.  
In the standard procedure of VEM framework, we define a computable discrete local bilinear form


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

approximating the continuous form a^K(·, ·) and satisfying the following properties:




  * 𝐤-consistency: for all 𝐪_k ∈ [_k(K)]^2 and 𝐯_h ∈𝐕_h^K


    a_h^K(𝐪_k, 𝐯_h) = a^K( 𝐪_k, 𝐯_h);



  * stability:  there exist  two positive constants α_* and α^*, independent of h and K, such that, for all 𝐯_h ∈𝐕_h^K, it holds


    α_* a^K(𝐯_h, 𝐯_h) ≤ a_h^K(𝐯_h, 𝐯_h) ≤α^* a^K(𝐯_h, 𝐯_h).




Let ℛ^K 𝐕_h^K ×𝐕_h^K → be a (symmetric) stabilizing bilinear form, satisfying

    c_* a^K(𝐯_h, 𝐯_h) ≤ℛ^K(𝐯_h, 𝐯_h) ≤ c^* a^K(𝐯_h, 𝐯_h)     for all 𝐯_h ∈𝐕_h such that Π_k^0,K𝐯_h= 0

with c_* and c^*  positive  constants independent of h and K.

Then, we can set

    a_h^K(𝐮_h, 𝐯_h) := a^K (Π_k^0,K𝐮_h, Π_k^0,K𝐯_h ) + ℛ^K ((I -Π_k^0,K) 𝐮_h, (I -Π_k^0,K) 𝐯_h )

for all 𝐮_h, 𝐯_h ∈𝐕_h^K.

It is straightforward to check that Definition (<ref>) and properties (<ref>) imply the consistency and the stability of the bilinear form a_h^K(·, ·).


In the construction of the stabilizing form ℛ^K with condition (<ref>)  we essentially require that the stabilizing term ℛ^K(𝐯_h, 𝐯_h) scales as a^K(𝐯_h, 𝐯_h). 
Following the standard VEM technique (cf. <cit.> for more details), denoting with 𝐮̅_h, 𝐯̅_h ∈^N_K
the vectors containing the values of the N_K local degrees of freedom associated to 𝐮_h, 𝐯_h ∈𝐕_h^K, we set

    ℛ^K (𝐮_h, 𝐯_h) =  α^K  𝐮̅_h^T 𝐯̅_h ,
 
where α^K is a suitable positive constant that scales as |K|. For example, in the numerical tests presented in Section <ref>, we have chosen α^K as the mean value of the eigenvalues of the matrix stemming from the  
term a^K (Π_k^0,K𝐮_h, Π_k^0,K𝐯_h ) in  (<ref>).


Finally we define the global approximated bilinear form a_h(·, ·) 𝐕_h ×𝐕_h → by simply summing the local contributions:


    a_h(𝐮_h, 𝐯_h) := ∑_K ∈𝒯_h  a_h^K(𝐮_h, 𝐯_h)     for all 𝐮_h, 𝐯_h ∈𝐕_h.










 §.§ The discrete problem


We are now ready to state the proposed discrete problem. Referring to (<ref>), (<ref>),   (<ref>), and   (<ref>)  we consider the virtual element problem:


    {   find (𝐮_h, p_h) ∈𝐕_h × Q_h, such that
        a_h(𝐮_h, 𝐯_h) + b(𝐯_h, p_h) = 0        for all 𝐯_h ∈𝐕_h,
         b(𝐮_h, q_h) = (f, q_h)        for all q_h ∈ Q_h..



We point out that the symmetry of a_h(·, ·) together with (<ref>)  easily implies that a_h(·, ·)  is  (uniformly) continuous with respect to the  L^2 norm. 

Moreover, as observed in <cit.>, introducing the discrete kernel:


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


it is immediate to check that


    𝐙_h ⊆𝐙 .


Then the bilinear form a_h(·, ·) is also uniformly coercive on the discrete kernel 𝐙_h with respect to the  𝐕 norm.




Moreover as a direct consequence of Proposition 4.3 in <cit.>, we have the following stability result.

 Given the discrete spaces
𝐕_h and Q_h defined in (<ref>) and (<ref>), there exists a positive β̃, independent of h, such that:


    sup_𝐯_h ∈𝐕_h  𝐯_h ≠0b(𝐯_h, q_h)/𝐯_h_𝐕≥β̃q_h_Q     for all q_h ∈ Q_h.




In particular, the  the inf-sup condition of Proposition <ref>, along with property (<ref>), implies that:


    div 𝐕_h = Q_h .


Finally we can state the well-posedness of virtual problem (<ref>).


Problem (<ref>) has a unique solution (𝐮_h, p_h) ∈𝐕_h × Q_h, verifying the estimate

    𝐮_h_𝐕 + p_h_Q≤ C f_0.













§ THEORETICAL RESULTS



We begin by proving an approximation result for the virtual local space 𝐕_h. First of all, let us recall a classical result by Brenner-Scott (see <cit.>).



Let K ∈𝒯_h, then for all 𝐮∈ [H^s+1(K)]^2 with 0 ≤ s ≤ k, there exists a  polynomial function 𝐮_π∈ [_k(K)]^2,  such that

    𝐮 - 𝐮_π_0,K + h_K |𝐮 -𝐮_π |_1,K≤ C h_K^s+1| 𝐮|_s+1,K.



We have the following approximation results (for the  proof see  <cit.>).



Let 𝐮∈𝐕∩ [H^s+1(Ω)]^2 with 0 ≤ s ≤ k. Under the assumption (𝐀1) and (𝐀2) 
on the decomposition 𝒯_h, there exists 𝐮_int∈𝐖_h  such that

    𝐮 - 𝐮_int_0 +  h_K |𝐮 -𝐮_int |_1,K≤ C h_K^s+1| 𝐮|_s+1,K.

where C is a constant independent of h.

For what concerns the pressures, from classic polynomial approximation theory <cit.>, for  q ∈ H^k(Ω) it holds

    inf_𝐪_h∈𝐐_h q - q_h _Q≤ C   h^k   |q|_k.

We are ready to state the following convergence theorem.


Let (𝐮, p) ∈𝐕× Q be the solution of problem (<ref>) and (𝐮_h, p_h) ∈𝐕_h × Q_h be the solution of problem (<ref>). Then it holds

    𝐮 - 𝐮_h_0≤ C   h^k+1  |𝐮|_k+1 ,     and    𝐮 - 𝐮_h_𝐕≤ C   h^k  |𝐮|_k+1 , 
    p - p_h_Q≤ C   h^k ( |𝐮|_k+1 + |p|_k).
  



We begin by remarking that as a consequence of the inf-sup condition  with classical arguments (see for
instance Proposition 2.5 in <cit.>), there exists 𝐮_I ∈𝐕_h such that 

    Π_k-1^0,K ( div 𝐮_I) =  div 𝐮_I = Π_k-1^0,K( div 𝐮)     for all K ∈𝒯_h,
    𝐮 - 𝐮_I_0≤ C  inf_𝐯_h ∈𝐕_h𝐮 - 𝐯_0    and        𝐮 - 𝐮_I_𝐕≤ C  inf_𝐯_h ∈𝐕_h𝐮 - 𝐯_𝐕.

Let us set δ_h = 𝐮_I - 𝐮_h. From (<ref>) and   (<ref>), we have that div δ_h = 0 and thus  δ_h ∈𝐙_h. Now, using (<ref>), (<ref>), (<ref>) and introducing the piecewise polynomial approximation (<ref>) together with (<ref>), we have

    α_*  α δ_h^2_0    ≤α_*   a(δ_h,  δ_h) ≤ a_h(δ_h,  δ_h)  = a_h(𝐮_I,  δ_h) - a_h(𝐮_h,  δ_h) 
        = a_h(𝐮_I,  δ_h) + b(δ_h, p_h) =  a_h(𝐮_I,  δ_h) 
        = ∑_K ∈𝒯_h a_h^K(𝐮_I,  δ_h) =
       ∑_K ∈𝒯_h ( a_h^K(𝐮_I - 𝐮_π,  δ_h) + a^K(𝐮_π , δ_h) ) 
        =  ∑_K ∈𝒯_h ( a_h^K(𝐮_I - 𝐮_π,  δ_h) + a^K(𝐮_π - 𝐮,   δ_h) ) - a( 𝐮,  δ_h)  
        = ∑_K ∈𝒯_h ( a_h^K(𝐮_I - 𝐮_π,  δ_h) + a^K(𝐮_π - 𝐮,   δ_h) ) + b( δ_h   , p)  
        = ∑_K ∈𝒯_h ( a_h^K(𝐮_I - 𝐮_π,  δ_h) + a^K(𝐮_π - 𝐮,   δ_h) ) 
       ≤ C ∑_K ∈𝒯_h ( 𝐮_I - 𝐮_π_0,K + 𝐮 - 𝐮_π_0,K)  δ_h_0, K
       ≤ C   ( 𝐮_I - 𝐮_π_0 + 𝐮 - 𝐮_π_0)  δ_h_0

then

    δ_h_0≤ C  𝐮_I - 𝐮_π_0 + 𝐮 - 𝐮_π_0 .
 
The L^2-estimate follows easily by the triangle inequality.
It is also straightforward to see from (<ref>) and (<ref>) that

    b(𝐮 - 𝐮_h, q_h) = 0     for all q_h ∈ Q_h,

than we get div 𝐮_h = Π_k-1^0,K ( div 𝐮) for all K ∈𝒯_h and therefore

    div (𝐮 - 𝐮_h)_0 = ∑_K ∈𝒯_h div 𝐮 - Π_k-1^0,K ( div 𝐮)_0,K≤ C   h^k   | div 𝐮|_k≤ C  h^k |𝐮|_k+1,

from which the estimate in the 𝐕 norm.

We proceed by analysing the error on the pressure field. 

Let q_h ∈ Q_h, then from the discrete inf-sup condition (<ref>), we infer:


    β̃p_h - q_h_Q ≤sup_𝐯_h ∈𝐕_h  𝐯_h ≠0b(𝐯_h, p_h - q_h)/𝐯_h_V = sup_𝐯_h ∈𝐕_h  𝐯_h ≠0b(𝐯_h, p_h - p) + b(𝐯_h,  p - q_h)/𝐯_h_V.


Since (𝐮,p) and (𝐮_h,p_h) are respectively the solution of (<ref>) and (<ref>), it follows that


    a(𝐮, 𝐯_h) + b(𝐯_h, p) = 0     for all 𝐯_h ∈𝐕_h,
    
    a_h(𝐮_h, 𝐯_h) + b(𝐯_h, p_h) = 0     for all 𝐯_h ∈𝐕_h.


Therefore, we get


    b(𝐯_h, p_h - p)  = a(𝐮, 𝐯_h) - a_h(𝐮_h, 𝐯_h)      for all 𝐯_h ∈𝐕_h.


Using (<ref>), the continuity of a_h(·, ·) and the triangle inequality, we get:


    b(𝐯_h, p_h - p)    = a(𝐮, 𝐯_h)  - a_h(𝐮_h, 𝐯_h)  
    = ∑_K ∈𝒯_h( a^K(𝐮, 𝐯_h)  - a_h^K(𝐮_h, 𝐯_h)) 
       = ∑_K ∈𝒯_h( a^K(𝐮 - 𝐮_π, 𝐯_h)  + 
    a_h^K(𝐮_π  -𝐮_h, 𝐯_h) ) 
       ≤∑_K ∈𝒯_h C ( 𝐮 - 𝐮_π_𝐕,K + 
    (𝐮_π - 𝐮_h)_𝐕,K) 𝐯_h_𝐕,K
       ≤∑_K ∈𝒯_h C  ( 𝐮 - 𝐮_π_𝐕,K + 𝐮 -𝐮_h_𝐕,K) 
    𝐯_h_𝐕,K

where 𝐮_π  is the piecewise polynomial of degree k defined in Lemma <ref>. Then, from estimate (<ref>) and the previous estimate on the velocity error, we obtain

    |b(𝐯_h, p_h - p)|    ≤  C h^k   |𝐮|_k+1 𝐯_h_𝐕.

Moreover, we have

    |b(𝐯_h, p - q_h) | ≤ C p - q_h_Q𝐯_h_𝐕.

Then, using (<ref>) and (<ref>) in (<ref>), we infer


    p_h - q_h_Q ≤ C h^k   |𝐮|_k+1 + C  p - q_h_Q.


Finally, using (<ref>) and the triangular inequality, we get

    p -p_h_Q ≤p -q_h_Q + p_h - q_h_Q ≤ C h^k   |𝐮|_k+1   + C  p - q_h_Q    for all q_h ∈ Q_h.

Passing to the infimum with respect to q_h ∈ Q_h, and using estimate (<ref>), we get the thesis.





We observe that the estimates on the velocity errors in Theorem <ref> do not depend on the continuous pressure, whereas the velocity errors of the classical methods have a pressure contribution. Therefore the proposed scheme belongs to the class of the pressure-robust methods.













§ A STABLE VEM FOR BRINKMAN EQUATIONS





 §.§ The continuous problem

The Brinkman equation describes fluid flow in complex porous media with a viscosity coefficient highly varying so that the flow is dominated by the Darcy equations in some regions of the domain  and by the Stokes equation in others. We consider the Brinkman equation on a polygon Ω⊆^2 with homogeneous Dirichlet boundary
conditions:

    {   μ Δ𝐮  +  ∇ p  +  ^-1𝐮 = 𝐟          in Ω,
        div 𝐮 = 0           in Ω,
       𝐮  = 0          on ∂Ω,.

where 𝐮 and p are the unknown velocity and  pressure fields, μ is the fluid viscosity,  denotes the permeability tensor of the porous media and 𝐟∈ [L^2(Ω)]^2 is the external source term. 

We assume  that   is a symmetric positive definite tensor and that there exist two positive (uniform) constants  λ_1, λ_2 > 0 such that

    λ_1  η^T η≤η^T ^-1η≤λ_2  η^T  η    for all η∈^2.


For what concerns the fluid viscosity we consider 0 < μ≤ C, this include  the case where μ approaches zero and  equation (<ref>) becomes a singular perturbation of the classic Darcy equations.

Let us consider the spaces

    𝐕:= [H_0^1(Ω)]^2,      Q:= L^2_0(Ω)

with the usual norms, and let  A(·, ·) 𝐕×𝐕→ be the bilinear form defined by:

    A (𝐮,  𝐯) :=  a^∇ (𝐮,   𝐯) + a(𝐮,   𝐯) ,     for all 𝐮,  𝐯∈𝐕


where

    a^∇ (𝐮,   𝐯) := ∫_Ωμ ∇𝐮 : ∇𝐯  dx     for all 𝐮,  𝐯∈𝐕

and a(·, ·) is the bilinear form defined in (<ref>).

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

    {   find (𝐮, p) ∈𝐕× Q, such that
        A(𝐮, 𝐯) + b(𝐯, p) = (𝐟, 𝐯)         for all 𝐯∈𝐕,
         b(𝐮, q) = 0        for all q ∈ Q,.

where and b(·, ·) 𝐕× Q →  is the bilinear form in (<ref>) and using standard notation

    (𝐟, 𝐯)  = ∫_Ω𝐟·𝐯  dx.
 
The natural energy norm for the velocities is induced by the symmetric an positive bilinear form A(·, ·) and is defined by (e.g. <cit.>)

    𝐯^2_𝐕, μ := A(𝐯, 𝐯) = μ ∇𝐯_0^2 + ^-1/2𝐯_0^2.

We can observe that the equivalence with the 𝐕 norm  is not uniform, i.e.

    c_1  √(μ) 𝐯_𝐕≤𝐯_𝐕, μ≤ c_2  𝐯_𝐕

where c_1, c_2 here and in the follows denote two positive constant independent of h and μ. 

For what concerns the pressures,  we consider the  norm  (see for instance <cit.>)

    p_Q, μ := sup_𝐯∈𝐕b(𝐯, p)/𝐯_𝐕, μ

Using the inf-sup condition in the usual norm it is possible to check the equivalence between the norms for the pressure but again the equivalence is not uniform, i.e.

    c_1  p_Q ≤p_Q, μ≤c_2/√(μ) p_Q.

Since, considering the modified norm, the bilinear form A(·, ·) is  uniformly continuous and coercive, and inf-sup condition is clearly fulfilled, Problem (<ref>) has a unique solution (𝐮, p) ∈𝐕× Q such that

    𝐮_𝐕, μ + p_Q, μ≤ C  𝐟_𝐕'

where the constant C depends only on Ω.






 §.§ Virtual formulation for Brinkman equations



Mathematically,  Brinkman equations can be viewed as a combination of  the Stokes and the
Darcy equation, that can change from place to place in the computational domain.
Therefore, numerical schemes for  Brinkman equations have to be
carefully designed to accommodate both  Stokes and Darcy simultaneously.
In this section we propose a Virtual Element scheme  that is 
accurate for both Darcy and Stokes flows.

For this goal we combine the ideas  developed in the previous sections with the argument in <cit.>.

Let us consider the virtual spaces 𝐕_h and Q_h (cfr. (<ref>) and (<ref>)). As usual in the VEM framework we need to define a computable approximation of the continuous bilinear forms.  Using obvious notations we split the bilinear form  A(·, ·) as


    A(𝐮,  𝐯) =: ∑_K ∈𝒯_h A^K (𝐮,  𝐯)  = ∑_K ∈𝒯_h( a^∇, K (𝐮,  𝐯) + a^K (𝐮,  𝐯))     for all 𝐮,  𝐯∈𝐕.


We begin by observing that, from <cit.> (in particular c.f. (27)-(29)) and from Section <ref>, A^K(𝐪_k, 𝐯) is computable on the basis of the DoFs 𝐃_𝐕 for all 𝐪_k ∈ [_k(K)]^2 and for all 𝐯∈𝐕_h. 
Starting from this observation we can approximate the continuous form A^K(·, ·) with the bilinear form 
A_h^K(·, ·) 𝐕_h^K ×𝐕_h^K →,
given by

    A_h^K(𝐮, 𝐯) = a_h^∇, K(𝐮, 𝐯) + a_h^K(𝐮, 𝐯)     for all 𝐮, 𝐯∈𝐕_h

where a_h^∇, K(·, ·) is the bilinear form defined in equation (35) in <cit.> and a_h^K(·, ·) is defined in (<ref>).

It is clear that the bilinear form A_h^K(·, ·) satisfies the k-consistency and the stability properties.






















As usual we build the global approximated bilinear form A_h(·, ·) 𝐕_h ×𝐕_h → by simply summing the local contributions.

For what concerns the bilinear form b(·, ·), as observed in Section <ref>, it can be computed exactly.

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 := Π_k^0,K𝐟    for all K ∈𝒯_h,

and consider:

    (𝐟_h, 𝐯_h)  = ∑_K ∈𝒯_h∫_K 𝐟_h ·𝐯_h   dK = ∑_K ∈𝒯_h∫_K Π_k^0,K𝐟·𝐯_h   dK = ∑_K ∈𝒯_h∫_K 𝐟·Π_k^0,K𝐯_h   dK.

We observe that (<ref>) can be exactly computed from 𝐃_𝐕 for all 𝐯_h ∈𝐕_h (see Proposition <ref>). 

Furthermore, the following result concerning a L^2 and  H^1-type norm, can be proved using standard arguments <cit.>.


Let 𝐟_h be defined as in (<ref>),  and let us assume 𝐟∈ H^k+1(Ω). Then, for all 𝐯_h ∈𝐕_h, it holds


    |( 𝐟_h - 𝐟, 𝐯_h ) | ≤ C h^k+1 |𝐟|_k+1𝐯_h_0    and    |( 𝐟_h - 𝐟, 𝐯_h ) | ≤ C h^k+2 |𝐟|_k+1 |𝐯_h|_𝐕.




In the light of the previous definitions, we consider the virtual element approximation of the Brinkman problem:


    {   find (𝐮_h, p_h) ∈𝐕_h × Q_h, such that
        A_h(𝐮_h, 𝐯_h) + b(𝐯_h, p_h) = (𝐟_h,  𝐯_h)      for all 𝐯_h ∈𝐕_h,
         b(𝐮_h, q_h) = 0     for all q_h ∈ Q_h..


Equation (<ref>) is well posed since the discrete bilinear form A_h(·,·) 
is (uniformly) stable with respect to the norm ·_𝐕, μ by construction and the inf-sup condition is fulfilled (the proof follows the guidelines of Proposition 4.2 in <cit.> and the linearity of the Fortin operator). Then we have the following result.


Problem (<ref>) has a unique solution (𝐮_h, p_h) ∈𝐕_h × Q_h, verifying the estimate

    𝐮_h_𝐕, μ + p_h_Q, μ≤ C 𝐟_𝐕'.




We now notice that, if 𝐮∈𝐕 is the velocity solution to Problem (<ref>), then it is the solution
to Problem:

    {   find 𝐮∈𝐙 such that
        A(𝐮, 𝐯) = (𝐟, 𝐯)     for all 𝐯∈𝐙 .

Analogously, if 𝐮_h ∈𝐕_h is the velocity solution to Problem (<ref>), then it is the solution to
Problem:

    {   find 𝐮_h ∈𝐙_h such that
        A_h(𝐮_h, 𝐯_h) = (𝐟_h, 𝐯_h)     for all 𝐯_h ∈𝐙_h .

For what concerns the convergence results  we state the following theorem. The proof can be derived by extending the techniques of the previous section and is therefore omitted.

Let 𝐮∈𝐙 be the solution of problem (<ref>) and 𝐮_h ∈𝐙_h be the solution of
problem (<ref>). Then

    𝐮 - 𝐮_h_𝐕, μ≤  C (√(μ)  h^k + ^-1/2_∞ h^k+1) |u|_k+1 + C  h^k+1  |f|_k+1

Let (𝐮, p) ∈𝐕× Q be the solution of Problem (<ref>) and (𝐮_h, p_h) ∈𝐖_h ∈ Q_h be the solution of Problem (<ref>). Then it holds:

    p - p_h _Q, μ≤ C ( h^k   |u|_k+1 +  h^k/√(μ) |p|_k + h^k+1  |f|_k+1).

The constants C above are independent of h and μ. 









In the last part of this section we present a brief discussion about the construction of a reduced virtual element method for Brinkman equations equivalent to Problem (<ref>) but involving significantly fewer degrees of freedom, especially for large k.
This construction essentially follows the guidelines of Section 5 in <cit.> (where we refer the reader for a deeper presentation).

Let us define the original reduced local virtual spaces, for k≥ 2:

    𝐖_h^K := {𝐯∈ [H^1(K)]^2   s.t  𝐯_|∂ K∈ [_k(∂ K)]^2,   {    -Δ𝐯  -  ∇ s ∈𝒢_k-2(K)^⊥,  
        div 𝐯∈_0(K),
    .     for some s ∈ H^1(K)  }

As before we enlarge the virtual space 𝐖_h^K and we consider

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

Finally we define the enhanced Virtual Element space, the restriction 𝐕_h^K of 𝐔_h^K given by

    𝐕_h^K := {𝐯∈𝐔_h^K   s.t.  (𝐯 - Π^∇,K_k 𝐯,  𝐠_k^⊥)_[L^2(K)]^2 = 0   for all 𝐠_k^⊥∈𝒢_k(K)^⊥/𝒢_k-2(K)^⊥} ,

where as before the symbol 𝒢_k(K)^⊥/𝒢_k-2(K)^⊥ denotes the polynomials i 𝒢_k(K)^⊥ that are L^2-orthogonal to all polynomials of 𝒢_k-2(K)^⊥.

For the pressures we consider the reduced space

    Q_h^K :=  _0(K).


As sets of degrees of freedom for the reduced spaces, combining the argument in Section <ref> and <cit.> we may consider the following. 
For every function 𝐯∈𝐕_h^K we take the following linear operators 𝐃_𝐕, split into three subsets (see Figure <ref>):


  * 𝐃_𝐕1: the values of 𝐯 at each vertex of the polygon K,

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

  * 𝐃_𝐕3: the moments  of 𝐯

    ∫_K 𝐯·𝐠_k-2^⊥  dK     for all 𝐠_k-2^⊥∈𝒢_k-2(K)^⊥.

 



For every q ∈Q_h we consider


  * 𝐃_𝐐: the moment

    ∫_K q    dK.


Therefore we have that:

    ( 𝐕_h^K ) = ([_k(∂ K)]^2) + (𝒢_k-2(K)^⊥) = 2n_K k + (k-1)(k-2)/2,

and

    (Q_h^K) = (_0(K))  = 1,

where n_K is the number of vertexes in K.  


We  define the global reduced virtual element spaces in the standard fashion.

The reduced virtual element discretization of the Brinkman problem (<ref>) is then:


    {   find 𝐮_h ∈𝐕_h and p_h ∈Q_h, such that
        A_h(𝐮_h, 𝐯_h) + b(𝐯_h, p_h) = (𝐟_h, 𝐯_h)        for all 𝐯_h ∈𝐕_h,
         b(𝐮_h, q_h) = 0        for all q_h ∈Q_h..


Above, the bilinear forms A_h(·, ·)  and b(·, ·), and the loading term 𝐟_h  are the same as before. 
The following proposition states the relation between Problem (<ref>) and the reduced Problem (<ref>) (the proof is equivalent to that of Proposition 5.1 in <cit.>).



Let (𝐮_h, p_h) ∈𝐕_h × Q_h be the solution of problem (<ref>) and (𝐮_h, p_h) ∈𝐕_h ×Q_h be the solution of problem (<ref>). Then

    𝐮_h = 𝐮_h       and    p_h|K = Π_0^0, K p_h   for all K ∈𝒯_h.













§ NUMERICAL TESTS




In this section we present two numerical experiments to test the practical performance of the method. The first experiment is focused on the method introduced in Section <ref> for the Darcy Problem, whereas in the second experiment we test the method in Section <ref> for the  Brinkman equations.

Since the VEM velocity solution 𝐮_h is not explicitly known point-wise inside the elements,  we  compute the method error comparing 𝐮 with a suitable polynomial projection of the approximated 𝐮_h. 












In particular we consider the computable error quantities:

    error(𝐮, H^1) := ( ∑_K ∈𝒯_h∇  u -  Π_k-1^0, K (∇  u_h) _0,K^2 )^1/2
     error(𝐮, H( div)) := ( ∑_K ∈𝒯_h div u -   div u_h _0,K^2  + ∑_K ∈𝒯_h u -  Π_k^0, K u_h _0,K^2 )^1/2
     error(𝐮, L^2)  := (∑_K ∈𝒯_h u -  Π_k^0, K u_h _0,K^2 )^1/2
     error(p, L^2)  :=p - p_h_0.


Regarding the computational domain, in our tests we always take the square domain Ω= [0,1] ^2, which is partitioned using the following sequences of polygonal meshes:


  * {𝒱_h}_h: sequence of Voronoi meshes with h=1/4, 1/8, 1/16, 1/32,

  * {𝒯_h}_h: sequence of triangular meshes with h=1/2, 1/4, 1/8, 1/16,

  * {𝒬_h}_h: sequence of square meshes with h=1/4, 1/8, 1/16, 1/32.

  * {𝒲_b}_h: sequence of WEB-like meshes with h= 4/10, 2/10,  1/10, 1/20.

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 non convex WEB-like meshes are composed by hexagons, generated starting from the triangular meshes {𝒯_h}_h and randomly displacing the midpoint of each (non boundary) edge.





In this example we consider the Darcy problem (<ref>) where  we set =I, and we choose  the load term 𝐟 in such a way that the analytical solution is 

    𝐮(x,y) =  -π [ sin(π x) cos(π y); cos(π x) sin(π y) ]     
    p(x,y) = cos(π x) cos(π y).

We analyse the practical performance of the virtual method by studying the errors versus the diameter h of the meshes. In addition we compare the results obtained with the scheme of Section <ref>, labeled as “div-free”, with those obtained with the method in the Appendix, labeled as “non div-free” (in both cases we consider polynomial degrees k=2). 





We notice that the "non div-free" method is a naive extension to the Darcy equation of the inf-sup stable scheme proposed in <cit.>. Since the scheme lacks a uniform ellipticity-on-the-kernel condition, it is not recommended for the problem under consideration. The purpose of the comparison is thus to underline the importance of the property Z_h⊆ Z (cf. Section 3.4) in the present context. 
 


In Figure <ref> and <ref>, we display the results for the sequence of Voronoi    
meshes 𝒱_h. In Figure <ref> and <ref>, we show the results for the sequence of meshes 𝒯_h, while 
in Figure <ref> and <ref> we plot the results for the sequence of 
meshes 𝒬_h, finally in <ref> and <ref> we exhibit the results for the sequence of meshes 𝒲_h.
 


















We notice that the theoretical predictions of Section <ref> and the Appendix are confirmed for both the L^2 norm and the H( div) norm. Note that for the H( div) norm we plot only the error for the “div-free” method since such scheme guarantees by construction, a better approximation of the divergence. Indeed  let u_h (resp. u_h) be the solution obtained with the “div-free” method (“non div-free” method)  then u_h satisfies

    div  u_h =   Π^0,K_k-1 f = Π^0,K_k-1 ( div  u)     for all K ∈𝒯_h

whereas u_h satisfies the same equation only in a projected sense, i.e.

    Π^0,K_k-1  ( div u_h)  =  Π^0,K_k-1  f  = Π^0,K_k-1  ( div u)       for all K ∈𝒯_h.

We can observe that the convergence rate of the L^2 norm for the pressure is optimal also for the “non div-free” method as proved in the Appendix.

Finally, we can observe that using a square mesh decomposition holds a convergence rate that is slightly better than what predicted by the theory. 






In this example we test the Brinkman equation (<ref>)  with different values of the fluid viscosity μ and fixed permeability tensor = I.  We choose the load term 𝐟 and the Dirichlet boundary conditions in such a way that the analytical solution is  

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

The aim of this test is to check the practical performance of the method introduced in Section <ref> in the reduced formulation (c.f. (<ref>). 
In Table <ref> and Table <ref> we display the total amount  of DoFs and the errors for the family of meshes 𝒱_h choosing k=2 respectively for the “div-free” method (cf. Section <ref>) and the “non div-free” method (cf. Reamrk <ref> and the Appendix).

We observe that also in the limit case, when the equation becomes a singular perturbation of the classic Darcy equations (e.g. for “small” μ), the proposed “div-free” method preserves the optimal order of accuracy.

















§ ACKNOWLEDGEMENTS

The author wishes to thank L. Beirão da Veiga and C. Lovadina for several interesting discussions and suggestions on the paper.
The author was 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.










§ APPENDIX: NON DIVERGENCE-FREE VIRTUAL SPACE




We have built a new H^1-conforming (vector valued) virtual space for the velocity vector field different from the more standard one presented in <cit.> for the elasticity problem. 
The topic of the present section is to analyse the extension to the Darcy equation of the scheme of  <cit.>.

Even though the method should not be used for the Darcy problem (cf. Remark 6.1), the numerical experiments have shown an optimal error convergence rate for the pressure variable. In this Section, we theoretically explain such a behaviour, under a convexity assumption on Ω (essentially, a regularity assumption on the problem). To this end, we develop an inverse estimate for the VEM spaces which is interesting on its own, and can be used in other contexts.  

We briefly describe the method by making use of various tools from the Virtual Element technology, and  we  
refer the interested reader to the papers <cit.>) for a deeper presentation.

We consider the local virtual space

    𝐖_h^K := {𝐯∈ [H^1(K)]^2   s.t  𝐯_|∂ K∈ [_k(∂ K)]^2   ,   Δ𝐯∈ [_k-2(K)]^2 }

with local degrees of freedom 𝐃_𝐕:


  * 𝐃_𝐕1: the values of 𝐯 at each vertex of the polygon K,

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

  * 𝐃_𝐕3: the moments  of 𝐯 up to order k-2, i.e.

    ∫_K 𝐯 ·𝐪_k-2  dK     for all 𝐪_k-2∈ [_k-2(K)]^2.
 
 

As observed in <cit.>, the DoFs 𝐃_𝐕 allow us to compute the operator Π_k^∇,K𝐖_h^K → [_k(K)]^2 defined as the analogous of the H^1 semi-norm projection (c.f. (<ref>)).










For all K ∈𝒯_h, the augmented virtual local space 𝐔_h^K is defined by

    𝐔_h^K = {𝐯∈ [H^1(K)]^2    s.t.  𝐯∈ [_k(∂ K)]^2,  Δ𝐯∈ [_k(K)]^2 }.

Now we define the enhanced Virtual Element space, the restriction 𝐕_h^K of 𝐔_h^K given by

    𝐕_h^K := {𝐯∈𝐔_h^K   s.t.  (𝐯 - Π_k^∇,K 𝐯,  𝐪_k )_[L^2(K)]^2 = 0   for all 𝐪∈ [_k(K)/_k-2 (K)]^2} ,

where the symbol _k(K)/_k-2 (K) denotes the polynomials of degree k living on K that are L^2-orthogonal to all polynomials of degree k-2 on K.
The enhanced space 𝐕_h^K has three fundamental properties (see <cit.> for a proof): 


  * [_k(K)]^2 ⊆𝐕_h^K,

  * the set of linear operators 𝐃_𝐕 constitutes a set of DoFs for the space 𝐕_h^K, 

  * the  L^2-projection operator Π^0, K_k𝐕_h^K → [_k(K)]^2 is exactly computable by the DoFs.






Recalling (<ref>) and from <cit.> it holds that dim(  𝐕_h^K ) =  dim(  𝐕_h^K ). 

For the pressures we use the space of the piecewise polynomials Q_h (c.f. (<ref>)).

For what concerns the construction of the approximated bilinear forms, it is straightforward to see that

    b(𝐯, q) = ∑_K ∈𝒯_h b^K(𝐯, q) = ∑_K ∈𝒯_h∫_K  div 𝐯  q   dK = - ∫_K 𝐯·∇ q   dK  + ∫_∂ K q  𝐯·𝐧

is computable from the DoFs for all 𝐯∈𝐕_h, q ∈ Q_h. Moreover using standard arguments <cit.> we can define a computable bilinear form 

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

approximating the continuous form a^K(·, ·), and satisfying the k-consistency (c.f. (<ref>)) and the stability properties (c.f. (<ref>)).  

Finally we define the global approximated bilinear form a_h(·, ·) 𝐕_h ×𝐕_h → by simply summing the local contributions.

By construction (see for instance <cit.>) the discrete bilinear form a_h(·,·) is 
(uniformly) stable with respect to the L^2 norm.

We are now ready to state the proposed discrete virtual element problem:


    {   find (𝐮_h, p_h) ∈𝐕_h × Q_h, such that
       a_h(𝐮_h, 𝐯_h) + b(𝐯_h, p_h) = 0        for all 𝐯_h ∈𝐕_h,
         b(𝐮_h, q_h) = (f, q_h)        for all q_h ∈ Q_h..



We shall first prove an inverse inequality for the virtual element functions in 𝐕_h.



Under the assumption (𝐀1), (𝐀2),  let K ∈𝒯_h and let 𝐯_h ∈𝐕_h^K. Then the following inverse estimate holds

    |𝐯_h|_1,k≤ c_inv  h_K^-1 𝐯_h_0,E

where the constant c_inv is independent of 𝐯_h, h_K and K.



We only sketch the proof, since we follow the guidelines of Lemma 3.1 and 3.3 in <cit.>. Let 𝐯_h ∈𝐕_h^K, then

    |𝐯_h|^2_1,K = ∫_K ∇𝐯_h ·∇𝐯_h = - ∫_K  Δ𝐯_h  𝐯_h + ∫_∂ K𝐯_h  ∇𝐯_h ·𝐧_K.

Under the assumption (𝐀1), (𝐀2) and by Lemma 3.3 in <cit.> we get

    - ∫_K  Δ𝐯_h  𝐯_h  ≤Δ𝐯_h_0,E𝐯_h_0,E≤ C_1   h_K^-1  |𝐯_h|_1,E𝐯_h_0,E

where the constant C_1 is independent of 𝐯_h, h_K and K. For what concerns the second addend in the right side of (<ref>), under the assumptions (𝐀1), (𝐀2), and using Lemma 3.1 in <cit.>, for all 𝐰∈ [H^1/2(∂ K)]^2 the following holds: there exists an extension 𝐰∈ [H^1(K)]^2 of 𝐰 such that

    h_K^-1 𝐰_0, K + |𝐰|_1, K≤ C  𝐰_1/2, ∂ K,

where we consider the scaled norm

    𝐰_1/2, ∂ K := h_K^-1/2 𝐰_0, ∂ K + |𝐰|_1/2, ∂ K.

By definition it holds

    ∫_∂ K𝐯_h  ∇𝐯_h ·𝐧_K ≤𝐯_h_1/2, ∂ K sup_𝐰∈ [H^1/2(∂ K)]^2⟨∇𝐯_h ·𝐧_K , 𝐰⟩/𝐰_1/2, ∂ K.

Now, using the definition (<ref>), an inverse estimate  (𝐯_h is polynomial on ∂ K) and the trace theorem <cit.>, it holds that

    𝐯_h_1/2, ∂ K   = h_K^-1/2 𝐯_h_0, ∂ K + |𝐯_h|_1/2, ∂ K≤ C_2   h_K^-1/2 𝐯_h_0, ∂ K
       ≤ C_2  h_K^-1/2  (  𝐯_h_0, K)^1/2(h_K^-1 𝐯_h_0, K +  |𝐯_h|_1, K)^1/2
       ≤(ϵ + C_2^2/ϵ)  h_K^-1𝐯_h_0, K + ϵ |𝐯_h|_1, K.

for any real ϵ >0. For the last term, using  (<ref>), (<ref>) and (<ref>) we get

    sup_𝐰∈ [H^1/2(∂ K)]^2⟨∇𝐯_h ·𝐧_K , 𝐰⟩/𝐰_h_1/2, ∂ K   ≤ C sup_𝐰∈ [H^1(K)]^2⟨∇𝐯_h ·𝐧_K , 𝐰⟩/h_K^-1 𝐰_0, K + |𝐰|_1, K
       ≤ C ( sup_𝐰∈ [H^1(K)]^2∫_K Δ𝐯_h  𝐰/h_K^-1 𝐰_0, K + |𝐰|_1, K + 
    sup_𝐰∈ [H^1(K)]^2∫_K ∇𝐯_h ·∇𝐰/h_K^-1 𝐰_0, K + |𝐰|_1, K) 
       ≤ C  ( sup_𝐰∈ [H^1(K)]^2∫_K Δ𝐯_h  𝐰/h_K^-1 𝐰_0, K + 
    sup_𝐰∈ [H^1(K)]^2∫_K ∇𝐯_h ·∇𝐰/ |𝐰|_1, K) 
       ≤ C  ( h_K  Δ𝐯_h_0, K + |𝐯_h|_1,K) ≤ C_3   |𝐯_h|_1,K.

From (<ref>) and (<ref>) we can conclude that

    ∫_∂ K𝐯_h  ∇𝐯_h ·𝐧_K ≤( (ϵ + C_2^2/ϵ)  h_K^-1𝐯_h_0, K + ϵ |𝐯_h|_1, K)   C_3|𝐯_h|_1, K

Finally, choosing ϵ = 1/2C_3 and collecting (<ref>) and (<ref>) in (<ref>) we have

    1/2|𝐯_h|_1, K≤( C_1 + 1/2 + 2 C_2^2 C_3^2)   h_K^-1𝐯_h_0, K

from which follows the thesis.



Let us analyse the theoretical properties of the method. We consider the discrete kernel:


    𝐙_h := {𝐯_h ∈𝐕_h   s.t.   b(𝐯_h, q_h) = 0   for all q_h ∈ Q_h} = {𝐯_h ∈𝐕_h   s.t.  Π_k-1^0,K ( div𝐯_h) = 0   for all K ∈𝒯_h},

therefore the divergence-free property is satisfied only in a relaxed (projected) sense. As a consequence the bilinear form a_h(·, ·) is not uniformly coercive on the discrete kernel 𝐙_h; nevertheless it holds the following h-dependent coercivity property


    a_h(𝐯_h,  𝐯_h) ≥α α_*  𝐯_h_0^2 ≥  C    h^2  𝐯_h^2_𝐕

that can be derived by using inverse estimate (<ref>).

Recalling that a_h(·, ·) is continuous with respect the  𝐕 norm and that the discrete inf-sup condition is fulfilled  <cit.>

    sup_𝐯_h ∈𝐕_h  𝐯_h ≠0b(𝐯_h, q_h)/𝐯_h_𝐕≥β̃q_h_Q     for all q_h ∈ Q_h

problem (<ref>) has a unique solution but we expect a worse order of accuracy since the bilinear form  a_h(·, ·) is not uniformly stable. In fact we have the following convergence results that are, perhaps surprisingly, still optimal in the pressure variable.



Let (𝐮, p) ∈𝐕× Q be the solution of problem (<ref>) and (𝐮_h, p_h) ∈𝐕_h × Q_h be the solution of problem (<ref>). Then

    𝐮 -𝐮_h _0≤ C   h^k-1 (|p|_k + h^2 |u|_k+1) ,     and    𝐮 - 𝐮_h _𝐕≤ C  h^k-2 (|p|_k + h^2 |u|_k+1).

Assuming further that Ω is convex, the following estimate holds: 

    p - p_h_Q≤ C   h^k (|p|_k + h^2 |u|_k+1).
  



As observed in the proof of Theorem <ref> the inf-sup condition (<ref>) implies the existence of a function 𝐮_I ∈𝐕_h such that 

    Π_k-1^0,K ( div 𝐮_I) = Π_k-1^0,K( div 𝐮)     for all K ∈𝒯_h,
    𝐮 - 𝐮_I_𝐕≤ C  inf_𝐯_h ∈𝐕_h𝐮 - 𝐯_h_𝐕

Now let us set δ_h = 𝐮_h - 𝐮_I. For what concerns the L^2 norm, using the stability of the bilinear form a_h(·, ·) and (<ref>) together with (<ref>)

    α_*  α δ_h^2_0    ≤α_*   a(δ_h,  δ_h) ≤a_h(δ_h, δ_h) = a_h(𝐮_h, δ_h) - a_h(𝐮_I, δ_h) 
        = - b(δ_h, p_h) - a(𝐮, δ_h) + a(𝐮, δ_h) - a_h(𝐮_I, δ_h)  
        =   b(δ_h, p - p_h) + ( a(𝐮, δ_h) - a_h(𝐮_I, δ_h)) =:
    μ_1(δ_h) + μ_2(δ_h).

By (<ref>) and property (<ref>), it is straightforward to see that

    Π_k-1^0,K ( div 𝐮_h) =  Π_k-1^0,K ( div 𝐮) = Π_k-1^0,K ( div 𝐮_I)     for all K ∈𝒯_h

so that δ_h ∈𝐙_h. Therefore

    μ_1(δ_h) = b(δ_h, p - p_h)=  b(δ_h, p) = b(δ_h, p - q_h)

for all q_h ∈ Q_h. Using the inverse estimate (<ref>) and standard approximation theory we get

    |μ_1(δ_h)| ≤ C δ_h_𝐕inf_q_h ∈ Q_hp - q_h_Q ≤ C  h^-1 δ_h_0  h^k  |p|_k = C  h^k-1 |p|_k δ_h_0.

By standard technique in VEM convergence theory, it holds that 

    |μ_2(δ_h)| ≤ h^k+1  |u|_k+1 δ_h_0.

Collecting  (<ref>) and (<ref>)  in (<ref>) we get the L^2 estimate. Whereas the 𝐕 norm estimate follows from an inverse estimate (<ref>).

For what concerns the estimate on the pressure, let  p_π the piecewise polynomial with respect to 𝒯_h defined by  p_π= Π_k-1^0,K  p for all K ∈𝒯_h. Let us set

    χ_h := 𝐮 - 𝐮_h,      z:= p_π - p     ρ_h:= p_π - p_h
 
From (<ref>) and (<ref>), it is straightforward to see that the couple (χ_h, ρ_h) solves the Darcy problem

    {    a(χ_h, 𝐯_h) + b(𝐯_h, ρ_h) =  (a_h(𝐮_h,  𝐯_h) - a(𝐮_h,  𝐯_h)) + b(𝐯_h, z)        for all 𝐯_h ∈𝐕_h,
         b(χ_h, q_h) = 0        for all q_h ∈ Q_h..

To prove the estimate for the pressure we employ the usual duality argument. Let therefore ϕ be the solution of the auxiliary problem 

    {   Δϕ = ρ_h         on Ω
       ϕ = 0                      on ∂Ω .

that, due to the convexity assumption, satisfies 

    ϕ_2 ≤ C  ρ_h_0

where the constant C depends only on Ω.
For all 𝐯∈𝐕_h  let us denote with 𝐯_I its interpolant defined in (<ref>) and (<ref>). Therefore Green formula together with (<ref>) yields

    ρ_h^2_0      = (ρ_h,  Δϕ) = b(∇ϕ,  ρ_h) = b((∇ϕ)_I,  ρ_h)  
        = (a_h(𝐮_h,   (∇ϕ)_I) - a(𝐮_h,   (∇ϕ)_I)) + b((∇ϕ)_I,   z) -   a(χ_h, (∇ϕ)_I) 
        =: μ_1((∇ϕ)_I) + μ_2((∇ϕ)_I) +  μ_3((∇ϕ)_I).

We analyse separately the three terms. For the first one, using the consistency property of a_h(·, ·), the polynomial approximation of 𝐮 and ∇ϕ, the estimate on the velocity error and (<ref>) we get

    μ_1((∇ϕ)_I)    = a_h(𝐮_h,   (∇ϕ)_I) - a(𝐮_h,   (∇ϕ)_I) 
        = ∑_K ∈𝒯_h ( a_h^K(𝐮_h,   (∇ϕ)_I) - a^K(𝐮_h,   (∇ϕ)_I) ) 
        =  ∑_K ∈𝒯_h(  a^K(𝐮_h - 𝐮_π,   (∇ϕ)_I -  (∇ϕ)_π) - a^K(𝐮_h - 𝐮_π,   (∇ϕ)_I  -  (∇ϕ)_π) ) 
       ≤ C  ∑_K ∈𝒯_h𝐮_h - 𝐮_π_0,K(∇ϕ)_I -  (∇ϕ)_π_0,K
       ≤ C  ∑_K ∈𝒯_h   (𝐮 - 𝐮_h_0,K  + 𝐮 - 𝐮_π_0,K) ((∇ϕ) -  (∇ϕ)_I_0,K + (∇ϕ) -  (∇ϕ)_π_0,K) 
       ≤ C    h^k-1 (|p|_k + h^2|u|_k+1)   h ∇ϕ_1 ≤ C    h^k-1 |p|_k  h ϕ_2 ≤ C   h^k (|p|_k + h^2|u|_k+1)  ρ_h_0

For what concerns the second term we have

    μ_2((∇ϕ)_I)    = b((∇ϕ)_I,   z) = b((∇ϕ)_I - ∇ϕ,   z)  + b(∇ϕ,   z) 
       ≤ C (|∇ϕ - (∇ϕ)_I |_1 + ϕ_2) z_0 ≤ C (|∇ϕ|_1 + ϕ_2) z_0 
       ≤ C   h^k |p|_k ϕ_2 ≤ C   h^k |p|_k ρ_h_0.

Finally, for the third term we begin by observing that from (<ref>) 

    b(χ_h,  ϕ) = b(χ_h,  ϕ - ϕ_π),

for all ϕ_π∈ Q_h, and by the Green formula

    b(χ_h,  ϕ) = - a(χ_h,  ∇ϕ) = - a(χ_h,  ∇ϕ - (∇ϕ)_I)  - a(χ_h,   (∇ϕ)_I).

Therefore, by collecting (<ref>), (<ref>), and using the previous error estimate,  it holds that

    μ_3((∇ϕ)_I)    =   - a(χ_h,   (∇ϕ)_I) = a(χ_h,  ∇ϕ - (∇ϕ)_I)  + b(χ_h,  ϕ - ϕ_π) 
       ≤ C (χ_h_0 ∇ϕ - (∇ϕ)_I_0 + χ_h_𝐕ϕ - ϕ_π_0) 
       ≤ C   h^k-1 (|p|_k + h^2|u|_k+1)   h  ϕ_2 + h^k-2 (|p|_k + h^2|u|_k+1)   h^2  ϕ_2 ) ≤ C h^k (|p|_k + h^2|u|_k+1) ρ_h_0.

Finally by collecting  (<ref>), (<ref>) and (<ref>) in (<ref>) we get the thesis.



tocsection
plain ]