The Stokes complex for Virtual Elements in three dimensions

The present paper has two objectives. On one side, we develop and test numerically divergence free Virtual Elements in three dimensions, for variable ``polynomial'' order. These are the natural extension of the two-dimensional divergence free VEM elements, with some modification that allows for a better computational efficiency. We test the element's performance both for the Stokes and (diffusion dominated) Navier-Stokes equation. The second, and perhaps main, motivation is to show that our scheme, also in three dimensions, enjoys an underlying discrete Stokes complex structure. We build a pair of virtual discrete spaces based on general polytopal partitions, the first one being scalar and the second one being vector valued, such that when coupled with our velocity and pressure spaces, yield a discrete Stokes complex.


Introduction
The Virtual Element Method (VEM) was introduced in [16,17] as a generalization of the Finite Element Method (FEM) allowing for general polytopal meshes. Nowadays the VEM technology has reached a good level of success; among the many papers we here limit ourselves in citing a few sample works [41,22,21,34,10,51,6,3]. It was soon recognized that the flexibility of VEM allows to build elements that hold peculiar advantages also on more standard grids. One main example is that of "divergence-free" Virtual Elements for Stokes-type problems, initiated in [14,5] and further developed in [15,57]. An advantage of the proposed family of Virtual Elements is that, without the need of a high minimal polynomial degree as it happens in conforming FEM, it is able to yield a discrete divergence-free (conforming) velocity solution, which can be an interesting asset as explored for Finite Elements in [44,46,53,48,45]. For a wider look in the literature, other VEM for Stokes-type problems can be found in [50,29,28,42,31,35] while different polygonal methods for the same problem in [49,33,24,38].
The present paper has two objectives. On one side, we develop and test numerically for the first time the divergence free Virtual Elements in three dimensions (for variable "polynomial" order k). These are the natural extension of the two-dimensional VEM elements of [14,15], with some modification that allows for a better computational efficiency. We first test the element's performance for the Stokes and (diffusion dominated) Navier-Stokes equation for different kind of meshes (such as Voronoi, but also cubes and tetrahedra) and then show a specific test that underlines the divergence free property (in the spirit of [48,15]).
The second, and perhaps main, motivation is to show that our scheme, also in three dimensions, enjoys an underlying discrete Stokes complex structure. That is, a discrete structure of the kind where the image of each operator exactly corresponds to the kernel of the following one, thus mimicking the continuous complex with Σ(Ω) denoting functions of L 2 (Ω) with curl in H 1 (Ω) [43,4]. Discrete Stokes complexes has been extensively studied in the literature of Finite Elements since the presence of an underlying complex implies a series of interesting advantages (such as the divergence free property), in addition to guaranteeing that the discrete scheme is able to correcly mimic the structure of the problem under study [37,36,8,7,9,27,40,39,52]. This motivation is therefore mainly theoretical in nature, but it serves the important purpose of giving a deeper foundation to our method. We therefore build a pair of virtual discrete spaces based on general polytopal partitions of Ω, the first one W h which is conforming in H 1 (Ω) and the second one Σ h conforming in Σ(Ω), such that, when coupled with our velocity and pressure spaces, yield a discrete Stokes complex. We also build a set of carefully chosen associated degrees of freedom. This construction was already developed in two dimensions in [19], but here things are more involved due to the much more complex nature of the curl operator in 3D when compared to 2D. In this respect we must underline that, to the best of the authors knowledge, no Stokes exact complex of the type above exists for conforming Finite Elements in three dimensions. There exist FEM for different (more regular) Stokes complexes, but at the price of developing cumbersome elements with a large minimal polynomial degree (we refer to [46] for an overview) or using a subdivision of the element [32]. We finally note that our construction holds for a general "polynomial" order k ≥ 2. The paper is organized as follows. After introducing some notation and preliminaries in Section 2, the Virtual Element spaces and the associated degrees of freedom are deployed in Section 3. In Section 4 we prove that the introduced spaces constitute an exact complex. In Section 5 we describe the discrete problem, together with the associated projectors and bilinear forms. In Section 6 we provide the numerical tests. Finally, in the appendix we prove a useful lemma.

Notations and preliminaries
In the present section we introduce some basic tools and notations useful in the construction and theoretical analysis of Virtual Element Methods.
Throughout the paper, we will follow the usual notation for Sobolev spaces and norms [1]. Hence, for an open bounded domain ω, the norms in the spaces W s p (ω) and L p (ω) are denoted by · W s p (ω) and · L p (ω) respectively. Norm and seminorm in H s (ω) are denoted respectively by · s,ω and |·| s,ω , while (·, ·) ω and · ω denote the L 2 -inner product and the L 2 -norm (the subscript ω may be omitted when ω is the whole computational domain Ω).

Basic notations and mesh assumptions
From now on, we will denote with P a general polyhedron having V vertexes V , e edges e and f faces f . For each polyhedron P , each face f of P and each edge e of f we denote with: n f P (resp. n P ) the unit outward normal vector to f (resp. to ∂P ), n e f (resp. n f ) the unit vector in the plane of f that is normal to the edge e (resp. to ∂f ) and outward with respect to f , t e f (resp. t f ) the unit vector in the plane of f tangent to e (resp. to ∂f ) counterclockwise with respect to n f P , τ f 1 and τ f 2 two orthogonal unit vectors lying on f and such that τ f t e a unit vector tangent to the edge e.
Notice that the vectors t e f , t f , τ f 1 and τ f 2 actually depend on P (we do not write such dependence explicitly for lightening the notations).
In the following O will denote a general geometrical entity (element, face, edge) having diameter h O .
Let Ω be the computational domain that we assume to be a contractible polyhedron (i.e. simply connected polyhedron with boundary ∂Ω which consists of one connected component), with Lipschitz boundary. Let { Ω h } h be a sequence of decompositions of Ω into general polyhedral elements P where h := sup P ∈Ω h h P .
We suppose that for all h, each element P in Ω h is a contractible polyhedron that fulfils the following assumptions: (A1) P is star-shaped with respect to a ball B P of radius ≥ h P , (A2) every face f of P is star-shaped with respect to a disk B f of radius ≥ h P , where is a uniform positive constant. We remark that the hypotheses (A1), (A2) and (A3), though not too restrictive in many practical cases, can be further relaxed, as investigated in [18,25,26,30].
The total number of vertexes, edges, faces and elements in the decomposition Ω h are denoted respectively with L V , L e , L f , L P .
For any mesh object O and for n ∈ N we introduce the spaces: for m ≤ n, denotes the polynomials in P n (O) with monomials of degree strictly greater than m.
Moreover for any mesh object O of dimension d we define and thus dim( P n\m (O)) = π n,d − π m,d . In the following the symbol will denote a bound up to a generic positive constant, independent of the mesh size h, but which may depend on Ω, on the "polynomial" order k and on the shape constant in assumptions (A1), (A2) and (A3).

Vector calculus & de Rham complexes
Here below we fix some additional notation of the multivariable calculus. Three dimensional operators. In three dimensions we denote with x = (x 1 , x 2 , x 3 ) the independent variable. With a usual notation the symbols ∇ and ∆ denote the gradient and Laplacian for scalar functions, while ∆, ∇, ε, div and curl denote the vector Laplacian, the gradient and the symmetric gradient operator, the divergence and the curl operator for vector fields. Note that on each polyhedron P the following useful polynomial decompositions hold [P n (P )] 3 = curl(P n+1 (P )) ⊕ x P n−1 (P ) .
Tangential operators. Let f be a face of a polyhedron P , we denote with x f := (x f 1 , x f 2 ) the independent variable (i.e. a local coordinate system on f associated with the axes τ 1 f and τ 2 f ). The tangential differential operators are denoted by a subscript f . Therefore the symbols ∇ f and ∆ f denote the gradient and Laplacian for scalar functions, while ∆ f , ∇ f , and div f denote the vector Laplacian, the gradient operator and the divergence for vector fields on f (with respect to the coordinate x f ). Furthermore for a scalar function ϕ and a vector field v := (v 1 , v 2 ) we set The following 2-d polynomial decompositions hold Noticing that v f is a 3-d vector field tangent to f , with a slight abuse of notations we define the 2-d vector field v τ on ∂P , such that on each face f its restriction to the face f satisfies The 3-d function v and its 2-d tangential restriction v τ are related by the curl-rot compatibility condition Moreover the Gauss theorem ensures the following rot-tangent component relation Finally, for any scalar function v defined in P , we denote with v τ the scalar function defined in On a generic mesh object O with geometrical dimension d, on a face f and on a polyhedron P we define following the functional spaces: with the "homogeneous counterparts" Remark 2.1. Notice that for each face f ∈ ∂P , the vector fields v f and v ∧ n f P are different. In fact both lie in the plane of the face f , but v f is π/2-rotation in f (with respect to the axes τ 1 and τ 2 ) of v ∧ n f P . However v f = 0 if and only if v ∧ n f P = 0. For that reason in the definition of Ψ 0 (P ), we consider a slightly different, but sustantially equivalent, set of homogeneous boundary conditions to that considered in literature [43,20,4].
Recalling that a sequence is exact if the image of each operator coincides with the kernel of the following one, and that P is contractible, from (2) and (3) it is easy to check that the following sequence is exact [11]: where i denotes the mapping that to every real number r associates the constant function identically equal to r and 0 is the mapping that to every function associates the number 0. The three dimensional de Rham complex with minimal regularity (in a contractible domain Ω) is provided by [8,37] In this paper we consider the de Rham sub-complex with enhanced smoothness [39] that is suitable for the Stokes (Navier-Stokes) problem. Therefore our goal is to construct conforming (with respect to the decomposition Ω h ) virtual element spaces that mimic the complex (7), i.e. are such that is an exact sub-complex of (7). To the best of our knowledge, no conforming finite elements sub-complex of (7) exists (see for instance [46]).

The virtual element spaces
The present section is devoted to the construction of conforming virtual element spaces (8) that compose the virtual sub-complex (9). As we will see, the space W h consists of the lowest degree three dimensional nodal VEM space [2,13], whereas the spaces V h and Q h (that are the spaces actually used in the discretization of the problem) are the three dimensional counterparts of the inf-sup stable couple of spaces introduced in [57,15]. Therefore the main novelty of the present section is in the construction of the Σ-conforming space Σ h .
In order to facilitate the reading, we present the spaces in the reverse order, from right to left in the sequence (9). In particular, in accordance with (9), the space Σ h will be careful designed to fit curl Σ h ⊆ V h . We stress that the readers mainly interested on the virtual elements approximation of the three dimensional Navier-Stokes equation (and not on the virtual de Rham sequence) can skip Subsection 3.3, Subsection 3.4 and Section 4.
One essential idea in the VEM construction is to define suitable (computable) polynomial projections. For any n ∈ N and each polyhedron/face O we introduce the following polynomial projections: with obvious extension for vector functions Π ∇,O Let k ≥ 2 be the polynomial degree of accuracy of the method. We recall that, in standard finite element fashion, the virtual element spaces are defined element-wise and then are assembled in such a way the global regularity requirements are satisfied.

Scalar L 2 -conforming space
We start our construction with the rightmost discrete space Q h in (9). Since we are not requiring any smoothness on Q h , the local space Q h (P ) is simply defined by having dimension (cf. (1)) dim(Q h (P )) = π k−1,3 . The corresponding DoFs are chosen, defining for each q ∈ Q h (P ) the following linear operators • D Q : the moments up to order k − 1 of q, i.e., P q p k−1 dP for any p k−1 ∈ P k−1 (P ).
The global space is given by It is straightforward to see that the dimension of Q h is

Vector H 1 -conforming VEM space
The subsequent space in the de Rham complex (9) is the vector-valued H 1 -conforming virtual element space V h . The construction of V h has to combine two main ingredients: • to define a 3-d version of the space [15] that fits the conformity requirement (that definition follows the guidelines of Appendix of reference [14]); • "to play" with the enhanced technique [2] in order to achieve the computability of the polynomial projections stated in Proposition 5.1.
We first consider on each face f of the element P , the face space and the boundary space that is a modification of the standard boundary nodal VEM [11]. Indeed the "super-enhanced" constraints (the last line in the definition (14)) are needed to exactly compute the polynomial projection Π 0,f k+1 (see Proposition 5.1). On the polyhedron P we define the virtual element space V h (P ) (15) The definition above is the 3-d counterpart of the virtual elements [15], in particular we remark that the enhancing constraints (the last line in (15)) are necessary to achieve the computability of the L 2 -projection Π 0,P k (see Proposition 5.1). Moreover, notice that the space V h (P ) contains [P k (P )] 3 and this will guarantee the good approximation property of the space (cf. Theorem 5.1).

Moreover, the following linear operators D V , split into five subsets constitute a set of DoFs for
Proof. We only sketch the proof since it follows the guidelines of Proposition 3.1 in [57] for the analogous 2-d space. First of all, recalling (1) and polynomial decomposition (2), simple computations yield and therefore number( Now employing Proposition 2 and Remark 5 in [2], it can be shown that the DoFs which in turn implies (recalling (15)) Now the result follows by proving that D V (v) = 0 implies that v is identically zero, that can be shown first works on ∂P and then inside P . As a consequence the linear operators D V are unisolvent for V h and in particular dim(V h (P )) = number(D V ).
The global space V h is defined by gluing the local spaces with the obvious associated sets of global DoFs: The dimension of V h is given by We also consider the discrete local kernel and the corresponding global version A crucial observation is that, extending to the 3-d case the result in [14], the proposed discrete spaces (12) and (18) are such that div V h ⊆ Q h . As a consequence the considerable kernel inclusion holds The inclusion here above and explicit computations (cf. (16)) yield that The notable property (21) leads to a series of important advantages, as explored in [48,46,15].
Remark 3.1. In the third line of Definition (15) the H 1 -seminorm projection Π ∇,P k can be actually replaced by any polynomial projection Π P k that is computable on the basis of the DoFs D V (in the sense of Proposition 5.1). This change clearly propagates throughout the rest of the analysis (see Definitions (26) and (37)). An analogous observation holds also for the operator Π ∇,f k in the third line of Definition (14). The present remark allows to make use of computationally cheaper projections, as done in the numerical tests of Section 6.

Vector Σ-conforming VEM space
In the present subsection we consider the construction of the Σ-conforming virtual space Σ h in (9). As mentioned before, this brick constitutes the main novelty in the foundation of the virtual de Rham sequence (9). The core ideas in building such space are the following: • the space is careful designed to satisfy curl Σ h = Z h ; • the DoFs are conveniently chosen in order to have a direct correspondence between the curl of the Lagrange-type basis functions of Σ h and the Lagrange basis functions of V h ; • the boundary space and the boundary DoFs are picked in accordance with the global conformity requirements ensuing from the regularity of the space Σ.
We start by introducing on each face f ∈ ∂P the face space and the boundary space (23) Concerning the differential problem in definition (22), we recall that on simply connected polygon f , given two sufficiently regular functions g and h defined on f and a sufficiently regular function ω defined on ∂f , the problem is well posed if and only if, in accordance with (5), the following holds On the polyhedron P we define the virtual space: We stress that the variational problem stated in (26) is coupled with the non homogeneous version of the boundary conditions in [43,4]. In fact, in order to force Σ-conforming regularity, for any function ϕ ∈ Σ h (P ), we need to prescribe curl ϕ and ϕ τ on ∂P . We address the well-posedness of the biharmonic problem in definition (26) in the Appendix. Note that, in accordance with the target curl Σ h ⊆ V h , the second and the last line in definition (26) are the curl version of the first and last line in definition (15). Whereas we will see that curl of the solutions of the biharmonic problem in (26) are solutions to the Stokes problem in (15) (see Proposition 4.2).
Moreover, the following linear operators D Σ , split into five subsets constitute a set of DoFs for Σ h (P ): Proof. We start the proof counting the number of the linear operators D Σ . Using similar computations as in (16) we have: For sake of simplicity, we prove that D Σ constitutes a set of DoFs for the non-enhanced space associated with Σ h (P ), i.e. the space Σ h (P ) obtained by dropping the last line in (26) (the enhanced constraints) and by taking in the biharmonic system Once the proof for Σ h (P ) is given, the extension to the original space Σ h (P ) easily follows by employing standard techniques for VEM enhanced spaces (see [2] and Proposition 5.1 in [19]). Employing Theorem 6.1, given • h ∈ S k (∂P ) satisfying the compatibility condition (cf. (4)) there exists a unique function ϕ ∈ Ψ(P ) such that where the last term (− dim( B k (f )) f ) ensues from the compatibility condition (28). We calculate the addenda in the right hand side of (29). Regarding the first term in (29), we preliminary note that the following characterization ensues from the exact sequence (6) and polynomial decomposition (2) Employing again the exact sequence (6), curl restricted to x ∧ [P k−3 (P )] 3 is actually an isomorphism, therefore from (30) and (16) follows that From definitions (22) and (23) and since problem (24) is well-posed, direct computations yield where the −1 in the formula above is due to the compatibility condition (25).
Having proved that number(D Σ ) is equal to dim( Σ h (P )), in order to validate that the linear operators D Σ constitute a set of DoFs for Σ h (P ) we have to check that they are unisolvent. Let ϕ ∈ Σ h (P ) such that D Σ (ϕ) = 0, we need to show that ϕ is identically zero. It is straightforward that D 1 Recalling the well known results for nodal boundary spaces [11], it is quite obvious to check that (33) D 3 Σ (ϕ) = 0 implies (curl ϕ) f = 0 for any f ∈ ∂P .
Remark 3.2. A careful inspection of Theorem 6.1 (see also Remark 5.1 in [43] and [20]) reveals that the space (26) admits the equivalent formulation The global space Σ h is defined by collecting the local spaces Σ h (P ), i.e.
The global set of DoFs is the global counterpart of D Σ , in particular the choice of DoFs D Σ establishes the conforming property curl Σ h ⊆ [H 1 (Ω)] 3 . The dimension of Σ h is given by

Scalar H 1 -conforming VEM space
In the present section we briefly define the H 1 -conforming space W h in the virtual complex (9). The space W h consists of low order nodal VEM [11].
We first introduce the low order boundary space and then we consider the VEM space on the polyhedron P with the associated set of DoFs: • D W : the values of v at the vertexes of the polyhedron P .
It is straightforward to see that the dimension of W h (P ) is given by dim(W h (P )) = V . The global space is obtained by collecting the local spaces with the obvious associated DoFs. The dimension of W h thus is given by

The virtual elements de Rham sequence
The aim of the present section is to show that the set of virtual spaces introduced in Section 3 realizes the exact sequence (9).   (41) and (38) respectively. Then Proof. Essentially we need to prove that For what concerns the inclusion (i1), every w ∈ W h clearly satisfies Therefore we need to verify that (∇ w) |P ∈ Σ h (P ) for any P ∈ Ω h . Notic0e that the tangential component of ∇ w verifies From definition (39) and (22), for any f ∈ ∂P we infer that, recalling (42), implies (∇ w) τ ∈ S k (f ). Moreover w ∈ C 0 (∂P ) entails and thus (cf. definition (23) Furthermore definition (40) implies Collecting (43) and (44) in definition (26), we easily obtain (i1).
We prove now the property (i2). Consider ϕ ∈ Σ h such that curl ϕ = 0. Since (7) is an exact sequence, there exists unique (up to constant) w ∈ H 1 (Ω) such that ∇ w = ϕ. Therefore for any face f in the decomposition Ω h , the tangential component of ∇ w satisfies (cf. definition (22) Hence on each face f the function w fulfils (∇ w · t e ) |e = (ϕ · t e ) |e ∈ P 0 (e) on any e ∈ ∂f , From (45) it follows that w restricted to the mesh skeleton is continuous and piecewise linear. Thus the function w is well defined (single valued) on the vertexes of the decomposition Ω h and D V ( w) makes sense. Let now w ∈ W h be the interpolant function of w in the sense of DoFs, i.e. the function uniquely determined by Inclusion (i1) guarantees that ∇ w ∈ Σ h . Hence, by Proposition 3.2, w realizes (i2) if and only if D Σ (∇ w) = D Σ (ϕ). Being curl(∇ w) = curl ϕ = 0, this reduce to verify that For any edge e in the decomposition Ω h , we denote with ν 2 and ν 1 the two endpoints of e, with t e pointing from ν 1 to ν 2 . Therefore, from (46) and (45), we infer that concludes the proof. (38) and (20) respectively. Then

Proposition 4.2. Let Σ h and Z h be the spaces defined in
Proof. The proof follows by showing the following points: Let us analyse the inclusion (i1). Let ϕ ∈ Σ h , clearly curl ϕ ∈ Z(Ω). Therefore we need to verify that (curl ϕ) |P ∈ V h (P ) for any P ∈ Ω h . It is evident that the second and the last line in definition (26) correspond respectively to the curl version of the first and last line in definition (15). Hence it remains to show that curl ϕ is the velocity solution of the Stokes problem associated with definition (15) on each element P . A careful inspection of the biharmonic problem in definition (26), imply that the following are equivalent P ∆ ϕ · ∆ ψ dP = P p k−1 · ψ dP ∀ψ ∈ Ψ 0 (P ), (by definition (26)) In particular the last equation is still valid considering all ψ ∈ Ψ 0 (P ) ∩ Z(P ). Therefore, using the identity ∆ = −curl curl + ∇ div and an integration by parts (coupled with the homogeneous boundary condition curl ψ = 0 on ∂P ), it can be proved that the following are equivalent P ∆ ϕ · ∆ ψ dP = P (x ∧ p k−1 ) · curl ψ dP ∀ψ ∈ Ψ 0 (P ) ∩ Z(P ), P ∆ ϕ · (−curl curl ψ) dP = P (x ∧ p k−1 ) · curl ψ dP ∀ψ ∈ Ψ 0 (P ) ∩ Z(P ), P −∆(curl ϕ) · curl ψ dP = P (x ∧ p k−1 ) · curl ψ dP ∀ψ ∈ Ψ 0 (P ) ∩ Z(P ).
Exploiting Lemma 5.1 in [43], for every z ∈ Z 0 (P ) there exists ψ ∈ Ψ 0 (P ) ∩ Z(P ) such that z = curl ψ. Therefore the last equation in (47) is equivalent to P ∇(curl ϕ) : ∇z dP = P (x ∧ p k−1 ) · z dP for all z ∈ Z 0 (P ), and therefore v = curl ϕ is the velocity solution of a Stokes problem as in Definition (15). That concludes the proof for (i1).
Let us consider the interpolant ϕ ∈ Σ h of ϕ in the sense of DoFs, i.e. the function uniquely determined by (cf. Proposition 3.2) Property (i1) ensures curl ϕ ∈ Z h . Therefore, employing Proposition 3.1, ϕ realizes (i2) if and only if D V (curl ϕ) = D V (v). Is it straightforward to check that for any face f .
In order to show that the two quantities above are equal we exploit the same computations in (34) and (48) This ends the proof. (18) and (12) respectively. Then

Proposition 4.3. Let V h and Q h be the spaces defined in
Proof. We follow same strategy adopted in the previous propositions and show that The inclusion (i1) is trivial. Regarding the point (i2), since (7) is an exact sequence, for any q ∈ Q h there exists v ∈ [H 1 (Ω)] 3 such that div v = q. Now let v ∈ V h the function uniquely determined by (cf. Proposition 3.1) Notice that being v ∈ [H 1 (Ω)] 3 the face moments in (49) and D 5 V ( v) are actually well defined. Therefore for any P ∈ Ω h we infer P (div v) p k−1 dP = P (div v) p k−1 dP = P q p k−1 dP for all p k−1 ∈ P k−1\0 (P ). (50) Moreover employing the divergence theorem, (49) implies Notice that (50) and (51) coincide with D Q (div v) = D Q (q) that coupled with div v ∈ Q h (from (i1)) concludes the proof.

Virtual Elements for the 3-d Navier-Stokes equation
We consider the steady Navier-Stokes equation on a polyhedral domain Ω ⊆ R 3 with homogeneous Dirichlet boundary conditions: where ν > 0 represents the viscosity, f ∈ [L 2 (Ω)] 3 is the external force and For sake of simplicity we here consider Dirichlet homogeneous boundary conditions, different boundary conditions can be treated as well.
It is well known [54] that in the diffusion dominated regime the Navier-Stokes equation (52) has a unique solution (u, p) with Moreover Problem (52) can be formulated in the equivalent kernel form:

Discrete forms and load term approximation
In this subsection we briefly describe the construction of a discrete version of the bilinear form a(·, ·) given in (53) and trilinear form c(·; ·, ·) given in (54). We can follow in a rather slavish way the procedure initially introduced in [16] for the laplace problem and further developed in [15] for flow problems. First, we decompose into local contributions the bilinear form a(·, ·) and the trilinear form c(·; ·, ·) by considering a(u, v) =: for all w, u, v ∈ [H 1 (Ω)] 3 . As usual in VEM framework the discrete counterpart of the continuous forms above are defined starting from the polynomial projections defined in (10) and (11). The following proposition extends to the 3-d case the result for the bi-dimensional spaces [14,57]. (14) and (15) respectively. The DoFs D V allow us to compute exactly the face projections

Proposition 5.1. Let B k (f ) and V h (P ) be the spaces defined in
for any f ∈ ∂P , and the element projections in the sense that, given any v h ∈ V h (P ), we are able to compute the polynomials using only, as unique information, the DoFs values Proof. The computability of the face projections is a direct application of Remark 5 in [2].
Concerning the element projections we here limit to prove the last item, the first two follow analogous techniques. By definition of L 2 -projection (10), in order the determine, for any v ∈ V h (P ), the polynomial Π 0,P k v we need to compute From polynomial decomposition (2) we can write The first and the last integrals are computable being Π ∇,P k v and Π 0,f k+1 computable. The second addend corresponds to the DoFs D 4 V . For the third addend we observe that, since div v is a polynomial of degree less or equal than k − 1 we can exactly reconstruct its value from the DoFs D 5 V and the normal face moments in D 3 V .
On the basis of the projections above, following a standard procedure in the VEM framework, we define the computable (in the sense of Proposition 5.1) discrete local forms and the approximated right hand side for all w, u, v ∈ V h (P ), where clearly and the symmetric stabilizing form S P : The condition above essentially requires the stabilizing term S P (v, v) to scale as |v| 2 1,P . For instance, a standard choice for the stabilization is the D-recipe stabilization introduced in [13].
Remark 5.1. The H 1 -seminorm projection Π ∇,P k in the stabilization term of Definition (56) can be replaced by any polynomial projection Π P k that is computable on the basis of the DoFs D V (in the sense of Proposition 5.1). A possible choice will be explored in Section 6.
The global virtual forms and the global approximated right-hand side are defined by simply summing the local contributions: for all w, u, v ∈ V h .

The discrete problem
Referring to the discrete spaces (18), (12), the discrete forms and the approximated load term (59) and the div form (55), the virtual element approximation of the Navier-Stokes equation is given by where V h,0 := V h ∩ [H 1 0 (Ω)] 3 and Q h,0 := Q h ∩ L 2 0 (Ω). Recalling the kernel inclusion (21), Problem (60) can be also formulated in the equivalent kernel form with the obvious notation Z h,0 := Z h ∩ [H 1 0 (Ω)] 3 . Combining the arguments in [14,15,26] it is possible to show that the virtual space V h has an optimal interpolation order of accuracy with respect to the degree k, and that the couple of spaces (V h , Q h ) is inf-sup stable [23]. The following existence and convergence theorem extends the analogous result for the bi-dimensional case [15].
for suitable functions F, H, K independent of h.
Note that, as a consequence of the important property (21), there is no direct dependence of the velocity error on the pressure solution.
Remark 5.2. Since Proposition 4.2 yields an explicit characterization of Z h as curl Σ h , one could follow (61) and build an equivalent curl (discrete) formulation (see for instance Problem (77) in [19]). Such approach is less appealing in 3-d since the curl operator has a non trivial kernel and thus some stabilization or additional Lagrange multiplier would be needed in the formulation. Moreover this approach does not seem to be competitive in terms of number of DoFs with the reduced version of the method (see Subsection 5.3). As a consequence, we do not explore any scheme resulting from the curl formulation.

Reduced spaces and reduced scheme
In the present section we briefly show that Problem (60) is somehow equivalent to a suitable reduced problem entangling relevant fewer DoFs, especially for large k. This reduction is analogous to its two-dimensional counterpart in Section 5 in [14] and Section 5.2 in [57].
The core idea is that D 5 V (u h ) = 0, where u h denotes the solution of (60), and therefore such degrees of freedom (and also the associated pressures) can be trivially eliminated from the system. Hence on each polygon P , let us define the reduced local spaces: and Q h (P ) := P 0 (P ) .
Exploiting the same tools of Proposition 3.1 it can be proved that the linear operators D V , split into four subsets, defined by constitute a set of DoFs for V h (P ). Concerning the space Q h (P ), it is straightforward to see that dim( Q h (P )) = 1 with unique DoF D Q defined by D Q (q) := P q dP . The global spaces V h and Q h are obtained in the standard fashion by gluing the local spaces: We remark that by construction Z h ⊆ V h , therefore employing Proposition 4.1 and Proposition 4.2 we can state the following result.
The following proposition is easy to check and states the relation between Problem (60) and the reduced Problem (66).

Numerical validation
In this section we numerically verify the proposed discretization scheme. Before dealing with such examples, we briefly describe an alternative (computationally cheaper) projection adopted in the implementation of the method. Then we outline the polyhedral meshes and the error norms used in the analysis.

An alternative DoFs-based projection
In the light of Remark 3.1 and Remark 5.1, the aim of the present subsection is to exhibit an alternative projection to be used in the place of the standard H 1 -seminorm projection Π ∇,P k in (15) and (56) that will turn out to be very easy to implement. An analogous alternative projection could also be used to substitute Π ∇,f k in (14). For any element P ∈ Ω h , let NDoF := dim(V h (P )). Then referring to Proposition 3.1 we set i.e. D v is the vector containing the degree of freedom values D V associated to v. We consider: • the DoFs-projection Π D,P Notice that Π D,P n is a special case of the serendipity projection introduced in [12]. Although the projection Π D,P n may seem awkward on paper, it is quite simple and cheap to implement on the computer (since it is nothing but an euclidean projection with respect to the degree of freedom vectors). Indeed, it can be checked that the matrix formulation Π D,P n of the operator Π D,P n acting from V h (P ) to V h (P ) (containing [P n (P )] 3 ) with respect to the basis V (cf. [17], formula (3.18)) is where D ∈ R NDoF×3πn,3 is the matrix defined by (cf. [17], formula (3.17)) D i,α := D V ,i (m α ) for i = 1, . . . , NDoF and α = 1, . . . , 3π n,3 , where using standard VEM notation, m α denotes the scaled monomial with x B barycenter of the polyhedron P , and α 1 , α 2 and α 3 suitable multi-indexes.

Meshes and error norms
We consider the standard [0, 1] 3 cube as domain Ω and we make four different discretizations of such domain: a) Structured refers to meshes composed by structured cubes inside the domain, Figure 1 (a). b) Tetra is a constrained Delaunay tetrahedralization of Ω, Figure 1  We would like to underline that the last type of mesh will severely test the robustness of the proposed method. Indeed, Random meshes are characterized by elements whose faces can be very small and distorted, see the detail in Figure 1 (d).
The tetrahedral meshes are generated via tetgen [56], while the last two are obtained by exploiting the c++ library voro++ [55]. To analyze the error convergence rate, we make, for each family, a sequence of four meshes with decreasing size. For each mesh we define the mesh-size as h : Let (u, p) and (u h , p h ) be the continuous and discrete VEM solution of the Stokes (or Navier-Stokes problem) under study. To evaluate how this discrete solution is close to the exact one, we use the following error measures, that make use of the local projection described in Proposition 5.1: • H 1 -velocity error: the theoretical expected convergence rate is h k (cf. (62)); • L 2 -pressure error: the expected rate is h k (cf. (63)).

Numerical tests
In this subsection we consider three different tests. In the first two examples, we numerically verify the theoretical trend of all the errors for a Stokes and Navier-Stokes problem. Finally, we propose two benchmark examples for the Stokes equation with the property of having the velocity solution in the discrete space V h . It is well known that classical mixed finite element methods lead in this situations 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 5.1, estimate (62)).

Example 1 (Stokes problem).
In this section we solve the Stokes problem on the unit cube [0 , 1] 3 , the discreted version being as in (60)  We consider the Structured, CVT and Random meshes. In Figure 2 we show the behaviour of the errors e u H 1 and e p L 2 . The slope of such errors are the expected ones, O(h k ) see Theorem 5.1. Moreover, for each approximation degree k the convergence lines associated with different meshes are close to each other and this represents a numerical evidence that the proposed method is robust with respect to the adopted meshes.  Example 2 (Navier-Stokes problem). In this subsection we consider the Navier-Stokes problem described in Equation (60) with Dirichlet boundary conditions. We consider the same discretization of the unit cube of the previuos example, i.e. the set of meshes Structured, CVT and Random. We define the right hand side and the boundary conditions in such a way that the exact solution is u(x, y, z) :=   sin(πx) cos(πy) cos(πz) cos(πx) sin(πy) cos(πz) −2 cos(πx) cos(πy) sin(πz)   and p(x, y, z) := sin(2πx) sin(2πy) sin(2πz) .
We solve the nonlinear problem by using standard Newton-Rapson iterations with a stopping criterion based on the displacement convergence test error with a tolerance tol=1e-10, i.e. until ||x n − x n+1 || < tol ||x n || where x n refers to the solution at the n-step. In Figure 3 we show the convergence lines of the H 1 error on the velocity and the L 2 error on the pressure, respectively. In all these cases we have the predicted trend: h k for the velocity and h k for the pressure, see Theorem 5.1. Moreover, also in this case the lines are close to each other varying the mesh discretization, expecially for the velocity solution. Note that, for the pressure solution and random meshes, higher order case k = 3, there seems to be a loss of accuracy at the second step. We believe this is due to difficulties related to the Newton convergence iterates (the associated linear system getting quite badly conditioned) since random meshes have a very bad geometry and we are reaching near the memory limit of our platform. We where unable to run a further step due to memory limits. Improving this aspect, possibly by exploring more advanced solvers or changing the adopted virtual element basis [34], is beyond the scope of the present paper.

Example 3 (Benchmark Problems).
In this paragraph, inspired by [47], we consider a particular example to numerically show the an advantage of the proposed method. It is well known that the error on the velocity field of standard inf-sup stable elements for the Stokes equation is pressure dependent [23]. Consequently, the accuracy of the discrete solution u h is affected by the discrete pressure error. As already shown for the two-dimensional case in [14], also in the three-dimensional case we do not have such dependency on the error, i.e. the error on the discrete velocity field u h does not depend on the pressure, but only on the velocity u and on the load term f (see Theorem 5.1, estimate (62)). Note that the present method, although div-free, is not pressure-robust in the sense of [47] since the error on the velocities is indirectly affected by the pressure through the loading approximation term [14]. Nevertheless it is still much better then the inf-sup stable element in this respect, as the accuracy of the load approximation (being a known quantity) can be easily improved.
To numerically verify such property we consider two Stokes problems where the exact velocity field is contained in V h u(x, y, z) := where k is the VEM approximation degree, but we vary the solution on the pressure. More specifically we will consider these two pressure solutions: a polynomial pressure p 1 (x, y, z) := x k y + y k z + z k x − 3 2(k + 1) , and an analytic pressure p 2 (x, y, z) := sin(2πx) sin(2πy) sin(2πz) .
Note that in both cases, since p i ∈ Q h for i = 1, 2, a standard inf-sup stable element of analogous polynomial degree would obtain O(h k ) error for the velocities in the H 1 norm even if u ∈ V h . In the first case, the velocity is a polynomial vector field of degree k, while the pressure is a polynomial of degree k and the load term f is a polynomial of degree k. In such configuration the presented VEM scheme yields the exact solution up to machine precision for the velocity field. Indeed, the velocity virtual element space contains polynomials of degree k and, since the load term is a polynomial of degree k, the term H(f , ν) in Equation (62) is close to the machine precision, i.e. we approximate exactly the load term f (cf. Definition (58)), so the error on u is close to the machine precision.
In table of Figure 4 left, we collect the errors e u H 1 only for the coarsest meshes composed by 27 and 68 elements for the Structured and Tetra meshes, respectively. In the second case the velocity is still a polynomial of degree k, but, since the pressure is a sinusoidal function, now the right hand side f is not a polynomial. Even if the velocity virtual element space contains the polynomial of degree k, the error is affected by the term H(f , ν) in Equation (62), i.e. it is affected by the polynomial approximation we make of the load term so the expected error for e u H 1 is h k+2 , which is still much better than O(h k ). In Figure 5 we show the convergence lines for both e u H 1 and e p L 2 . The error trends are the expected ones: we get O(h 4 ), O(h 5 ) and O(h 6 ) for degrees k = 2, 3 and 4, respectively, while we get O(h k ) for the pressure. In the last step of the H 1 norm error, the error is higher than expected (this behaviour is due to machine algebra effects since we are in a range of very small errors). The biharmonic problem coupled with the non homogeneous boundary conditions find ϕ ∈ Ψ(P ), such that P ∆ ϕ · ∆ ψ dP = P f · ψ dP for all ψ ∈ Ψ 0 (P ), ∂P ϕ · n P df = 0 , has a unique solution ϕ.
Proof. Let us consider the following auxiliary problem find ϕ ∂ ∈ Ψ(P ), such that ∂P ϕ ∂ · n P df = 0 , We construct by hand a suitable ϕ ∂ that satisfies (71). Let us consider the Stokes-type problem defined on P then by Theorem 3.4 [43], there exists a vector potential ϕ g (possibly not unique) satisfying      ϕ g ∈ Ψ(P ), such that curl ϕ g = u in P , div ϕ g = 0 in P .
Therefore it can be shown that there exists ζ ∈ H 1 (∂P ) such that Now we consider the elliptic problem ∆ ω = 0 in P , ω = ζ on ∂P .