Serendipity Virtual Elements for General Elliptic Equations in Three Dimensions

LourenÇo BEIRÃO DA VEIGAFranco BREZZIFranco DASSI Luisa Donatella MARINIAlessandro RUSSO

(Dedicated to Philippe G.Ciarlet on the occasion of his 80th birthday)

1 Introduction

The virtual element method(VEM for short)was introduced in[5,7]as a generalization of thefinite element method(FEM for short)that allows for very general polygonal and polyhedral meshes,also including non convex and very distorted elements.Differently from standard FEM,the VEM is not based on the explicit construction and evaluation of the basis functions,but rather on a wise choice and a suitable use of the degrees of freedom in order to approximate the operators and the corresponding bilinear forms involved in the problem.The local functions are virtual,in the sense that they are defined,in general,through a partial differential equation(or even a system);they include(but in general are not restricted to)polynomials.However,the non-polynomial functions are never computed in practice,and the accuracy of the method is ensured by the polynomial part of the virtual space.The use of such an approach introduces other potential advantages,such as exact satisfaction of linear constraints as in[12]or[3],and the possibility to build easily discrete spaces of high global regularity[2—3,16].Since its introduction,the VEM has shared a good degree of success and has been applied to a large array of problems.We here mention,in addition to the ones above,a sample of papers[1,4,8,13,17,19,22,27]and refer to[25]for a more complete survey of the existing VEM literature.

Although the construction of the virtual element method for several three dimensional problems is accomplished already in many papers,at the current level of development a detailed presentation of their properties for general second order elliptic operators is still lacking.Moreover very few 3D numerical experiments are available in the literature[19,22—23]and all of them are limited to the lowest order case(k=1)while the only work dealing with the higher order case is[10].

The objective of this work is to extend to the three-dimensional case the work of[8].From another point of view,we could consider it as an extension of[10]to general elliptic second order operators.We have also implemented the serendipity version of the VEM(see[9])applied to faces in three dimensions which allows for a strong reduction of the face degrees of freedom without spoiling the approximation properties.Moreover we also numerically validate the theory with suitable numerical experiments.In doing this we will follow faithfully the construction in[1,5,8,10]and the coding guidelines of[7].

The paper is organized as follows.In Section 2 we introduce the model problem and the virtual element method in three dimensions.The review of the method is complete but brief,and we refer to other contributions in the literature for a more detailed presentation of the scheme.In Section 3 we outline the convergence proof,following straightforwardly the arguments given in the two-dimensional case.Finally,in Section 4 several numerical tests are shown.

2 The Virtual Element Discretization

In the present section we give a very brief overview of the virtual element method in three space dimensions,and in particular to its variant using serendipity elements on faces.More details on several aspects of the method can be found in[1,5,7—8,10].

2.1 The differential problem

Let Ω ⊂ R3be a bounded convex polyhedron and let Γ be its boundary.Let moreover κ and γ be smooth functions Ω → R with κ(x) ≥ κ0＞ 0 for all x ∈ Ω,and let finally b be a smooth vector valued function Ω → R3.

We consider the problem

where f is a given right-hand side in H−1(Ω).

We assume that problem(2.1)is well posed.That is,we assume that the problem is solvable for any f ∈ H−1(Ω),and that the estimate

together with the regularity estimate

holds with a constant C independent of f.

We recall that these assumptions imply that existence and uniqueness hold as well for the(formal)adjoint operator L∗given by

Moreover,they imply that for every g ∈ L2(Ω)there exists a unique ϕ ∈ H2(Ω)∩(Ω)such that L∗ϕ =g,and

for a constant C∗independent of g.Actually,the 2-regularity in(2.3)and in(2.5)is not strictly necessary in order to get the results of the present work,and an s-regularity with s＞1 would be sufficient.Similarly,the convexity assumption on Ω could be by-passed.Here however we are not interested in minimizing the regularity assumptions.

We also point out that the choice of having a scalar diffusion coefficient was done just for simplicity.Having a full diffusion tensor would not change the analysis in a substantial way.Similarly,the use of Dirichlet boundary conditions on Γ is done just for the sake of simplicity,and other boundary conditions can be easily accommodated.

The variational form of our problem reads

where

When convenient,we will also use the notation

and we remark that our assumptions on the coefficients imply that the bilinear form B(·,·)verifies

and hence

It is also easy to check that this,together with(2.2),implies that

for some constant CB＞0 independent of u.On the other hand,it is also well known that(2.9)and(2.10)imply existence and uniqueness of the solution of problem(2.6).

Remark 2.1We point out that,together with(2.9)we also have

that will come at hand later on.

2.2 Some useful notation

In what follows k will denote a positive integer associated with the “polynomial degree” of the virtual element scheme.We will denote by P a polyhedron in R3,while edges,faces,and vertices of P will be indicated by e,F,and ν respectively.We will denote by xP,hPand|P|the centroid,the diameter,and the volume of P,respectively.The set of of polynomials of degree less than or equal to s in P will be indicated by Ps(P).If α =(α1,α2,α3)is a multi-index,we will indicate by mαthe scaled monomial

where|α|= α1+α2+α3.

Remark 2.2Note that,in our computations,we will always use the scaled monomials as a basis for the polynomial spaces.Other choices might be more convenient in some particular cases.See for instance[14,24,26].

The corresponding definitions in the case of a polygon in R2are completely analogous.A face F of a polyhedron is treated as a polygon in two dimensions,using local coordinates(x,y).Edges of polyhedra and polygons are treated in an analogous way as one-dimensional sets.

2.3 Virtual elements on faces

We start by briefly recalling the virtual element spaces on the faces of a generic polyhedron.Given a polygon F(representing a generic face),we define the preliminary virtual space

The standard VEM in three dimensions is actually the one proposed in[1],that carries much less degrees of freedom on faces with respect to the choice(2.14)(see also[10]).Here we avoid the presentation of such space,and recall instead directly its serendipity version,that produces an even larger reduction of the number of degrees of freedom.

2.4 Serendipity version of virtual elements on faces

The Serendipity version of virtual elements was introduced in[9]as a tool to decrease the number of internal degrees of freedom on polygons.We here recall them very briefly and refer to the above paper for a detailed description.

Serendipity virtual elements are based on the introduction of a projection operator

that can be computed using only the boundary degrees of freedom plus a suitable subset of the internal ones:Typically the moments up to a certain degree kF≤k depending on the geometry of the element F.For instance(following[9]),denoting by NSthe number of selected degrees of freedom,one can define the mapping D from~Vk(F)to RNSthat associates to v ∈ ~Vk(F)the vector of its corresponding NSdegrees of freedom.Then one definesv ∈ Pk(F)such that

and the precise formulation for our request on kFis:kFmust be such that(2.18)has a unique solution.

Using such projection operator,one can define the serendipity space

A set of degrees of freedom for the above space are the sets(D1)and(D2)given in(2.15)—(2.16),plus

Remark 2.3It is important to point out that for v∈VkS(F)all the moments

are computable:For polynomials p of degree≤kFwe get them directly form the degrees of freedom(D3′),and the others are obtainable throughv(computable from the degrees of freedom of v)and the definition(2.19).

Remark 2.4It is clear that the bigger is kF,the more degrees of freedom we retain in(F).To minimize their number one would like to choose the smallest possible kFthat ensures the unique solvability of(2.18).On the other hand,the computation of such a kFon every polygon could be a rather heavy burden.We refer to[9]for a deeper discussion of this matter.Following the terminology in[9],denoting by ηFthe minimum number of lines needed to cover the whole boundary of F,we can set kF=k−ηF.Since for strictly convex polygons ηFis equal to the number of edges,and in any case is always≥ 3,a lazy choice could be to take ηF=3 for all the elements,while a stingy choice would compute the exact ηFfor each element.Thefirst choice is robust but leads to a lesser gain in terms of degrees of freedom.The second choice leads to larger computational gain but its performance is mesh dependent(in a sense detailed in[9]).The quest for a cheap way to choose the most convenient kFis not over yet,and we will discuss it more at length in our future works.In Subsection 4.2.1 we will discuss the problem in more details.

Remark 2.5We point out that the serendipity approach depends heavily on the degrees of freedom.The present choice is based on the degrees of freedom of[9].In other circumstances different strategies might be needed,leading to different serendipity spaces(see,e.g.,[6]).

2.5 Virtual elements on polyhedrons

Let Ωhbe a partition of Ω into non-overlapping and conforming polyhedrons.We begin by defining the virtual space locally,on each polyhedron P ∈ Ωh.Note that each face F ∈ ∂P is a two-dimensional polygon.We introduce the following boundary space:

The above space is made of functions that on each face are two-dimensional virtual functions and are continuous across the edges.Once the boundary space is defined,we can construct the local virtual space on P following for instance[1].We recall the procedure briefly and refer to[10]for a more detailed description.We define first the preliminary space

and then we define the local virtual space as

It is easy to check that the following linear operators constitute a set of degrees of freedom for the space Vk(P):

The associated(global)degrees of freedom are the obvious counterpart of the local ones introduced above.

2.6 Discretization of the problem

We start by introducing the discrete counterpart of the bilinear forms(2.7).On every polyhedron P∈Ωhwe follow the same construction introduced in[8]for the two-dimensional case.Wefirst introduce the stabilization form

where dofi(v)is the value of the ith local degree of freedom of the function v,denotes the number of the degrees of freedom(2.24)—(2.27),and the positive constants diare given by the diagonal recipe proposed in[10](which takes into account also the values of κ on P).Then we introduce,for all v,w∈Vk(P),the local approximations of the bilinear forms in(2.7),

Remark 2.6Here we are implicitly assuming that the diffusion coefficient κ is not too small compared to the other coefficients b and γ (or,more precisely,locally to|b|hPand γ).When this is not the case,a suitable stabilizing term inand/orwill also be necessary.

The above bilinear forms are stable in the sense of[5].The respective global bilinear forms ah(·,·),bh(·,·),ch(·,·)are constructed as usual by summing over all elements P of Ωh.

We canfinally state the discrete problem:

where the approximate loading fhis the L2-projection of f on piecewise polynomials of degree k.Note that all the forms and operators appearing above are computable in terms of the degrees of freedom of uhand vh,as observed above.Often we will use

3 Convergence Results

In this section we extend to the three dimensional case the error estimates that were obtained in[9]for the two-dimensional case.Most steps are quite similar,and we just briefly sketch the general path of the proofs.

3.1 Interpolation results

We presentfirst some approximation results concerning the virtual element spaces of the previous section.The results are a simple extension(to the present case of serendipity nodal spaces on faces)of the interpolation estimates shown in[27]for two dimensions and extended in[18]to the three dimensional case(see also[11,15]).

We assume that{Ωh}his a family of meshes,satisfying the following assumptions(typical of the virtual element methods):There exists a positive constant σ such that all elements P of{Ωh}hand all faces of∂P are star-shaped with respect to a sphere(respectively,a disk)of radius bigger than or equal to σhP,hPbeing the diameter of P;moreover,all edges e ∈ ∂P for all P ∈ {Ωh}hhave length bigger than or equal to σhP.

Here and in the sequel C will denote a generic positive constant independent of hP,with different meaning in different occurrences,and generally depending on the coefficients of the operator L.Whenever needed to better follow the steps of the proofs,for a smooth scalar or vector-valued functionℵ,we shall use Cℵto denote a constant depending onℵand possibly on its derivatives up to the needed order.

Proposition 3.1There exists a positive constantC=C(σ,k)such that,for allPinΩhand all smooth enough functionsϕdefined onP,it holds

whereϕI∈ Vkis the interpolant ofϕinVk,i.e.,such thatdofi(ϕ)=dofi(ϕI).

3.2 Continuity results

Concerning the bilinear forms presented in the previous section,we state a continuity result,whose proof is a trivial extension of several classical results in the framework of virtual element methods.

Proposition 3.2The bilinear formBh(·,·)is continuous inVk× Vk,that is,

withCκ,b,γa positive constant depending onκ,b,γbut independent ofh.Moreover,similarly to(2.11),we have

3.3 Approximation of the bilinear forms

The following result will be needed to estimate the difference between continuous and discrete bilinear forms.This is done once and for all in the following preliminary lemma.

Lemma 3.1LetP∈Ωh,letµbe a smooth function onP,and letu,vdenote smooth scalar or vector-valued functions onP.For a genericϕ ∈ L2(P)(or in(L2(P))3)we define,forsinteger≤k,

Then we have the estimate

whereCµis a constant depending onµ.

ProofFor simplifying the notation we will set:=u,:=v.By adding and subtracting terms,and by the definition of projection we have

and the result follows by Cauchy-Schwarz inequality with Cµ= ‖µ‖∞.

The following result follows immediately by a direct application of Lemma 3.1.

Lemma 3.2For allP∈Ωhit holds

From all the above results,proceeding as in[9],we have the following consistency property that is reminiscent of the classical results forfinite element methods(as,e.g.,the crucial Theorem 4.1.4.in[20]).

Proposition 3.3For allusufficiently regular and for allvh∈Vk,it holds

whereuis again theL2-projection ofuontoPk(P).

Remark3.1 We point out that(3.9)holds for a generic vh∈Vk,for which only H1regularity can be used.If for instance vh=vI,that is,vhis the interpolant of a more regular function,(3.9)can be improved.Indeed,proceeding again as in[9]we would have

On the other hand,using(3.7)—(3.8)we can easily show for all u and v in H1(P),

3.4 Stability results

Before going to study the error estimates for our problem,we have to prove the following stability result.

Proposition 3.4The bilinear formBh(·,·)satisfies the following condition(discrete counterpart of(2.10)):There exists anh0＞0and a constantsuch that,for allh＜h0and for alluh∈Vk,

The proof follows the classical path of Schatz[29].Wefirst prove that for every v∗∈there exists a∈Vksuch that

and moreover,there exists a constant C independent of h,such that

This is done following essentially the path of[9]combined with the 3D estimates in[18].Then we recall that for uh∈Vk,using(2.10)we have

Then we take the correspondinggiven by(3.13),and we get

From this,using(3.2)and(3.11)we have

that joined to(3.15)gives

and the result follows easily.

Remark 3.2Clearly,if b=0,and γ=0,(3.12)holds for any h.

3.5 Convergence in H1

We are now ready to prove the following theorem.

Theorem 3.1Forhsufficiently small,problem(2.30)has a unique solutionuh∈Vk,and the following error estimate holds:

withCa constant depending onκ,b,andγbut independent ofh.

ProofThe existence and uniqueness of the solution of problem(2.30),for h small,are a consequence of Proposition 3.4.To prove the estimate(3.17),using(3.12)we have that for h≤h0,there exists a∈Vkverifying

Thefirst term in the right-hand side of(3.19)is bounded by the Cauchy-Schwarz inequality and standard approximation estimates on the load f.The second and fourth terms are bounded using the continuity of Bhand B,respectively,and the third term is bounded by the consistency estimate(3.16).Also using approximation estimates for polynomials and for the virtual element space,we get

and the proof is concluded.

Remark 3.3It is immediate to check that,by the same proof,also the following refined result holds:

3.6 L2estimate

Theorem 3.2Forhsufficiently small,the following error estimate holds:

ProofAs usual,we shall use duality arguments.Let ψ ∈ H2(Ω)∩(Ω)be the solution of the adjoint problem(see(2.4)):

and let ψI∈ Vkbe its interpolant,for which it holds

Then

and the result follows with the usual arguments.Thefirst term is bounded through(3.1),(3.17),and(3.22).For the second term we apply Cauchy-Schwarz and standard approximation estimates.The third and fourth terms are bounded through(3.10),taking k=0 for the third one,and standard approximation estimates.

4 Numerical Results

Figure 1 Three types of discretizations of the domain Ω and cross section of such meshes.

In this section we present some numerical tests.In Subsection 4.1 we focus on the standard VEM approach in 3D(that is the standard construction of[5]but using on faces the advanced space of[1],see for instance[10]),while Subsection 4.2 is devoted to the serendipity VEM approach.

In these examples the domain is the unit cube Ω:=[0,1]3,we take as exact solution the function u(x,y,z):=sin(πx)cos(πy)cos(πz)and we choose

The load term and the boundary data are set in accordance with the above data and solution.

In all the examples we will consider three different discretizations of Ω (see Figure 1):

Structureda structured mesh composed by cubes;

CVTa mesh composed by well-shaped Voronoi cells obtained via a standard Lloyd’s algorithm(see[21]);

Randoma Voronoi mesh composed by distorted cells.

To construct the last two types of meshes we use the c++library voro++(see[28]).Then,for each type of mesh we make a sequence of finer meshes with mesh-size h defined as h:=N being the number of polyhedrons in the mesh and hPthe diameter of the polyhedron P.

We follow the trend of the following errors:

•H1errorcomputed as

where Π∇kuhis the element wise VEM H1-projection on polynomials of degree k defined in(2.23),and|·|1,Ωdenotes the standard H1-seminorm;

•L2errorcomputed as

In accordance with Theorems 3.1—3.2,if we consider a VEM approximation degree k,we expect order k in H1,and k+1 in L2.

4.1 Test case 1:h-analysis with a standard approach

Figure 2 Test case 1:Convergence rates for standard VEM on all meshes.

Figure 2 shows the convergence curves of the errors for each set of meshes,and for various degrees:k=1,2,3,4.From these graphs we can see that both the H1and the L2errors behave as expected.Moreover,the trend of the error is not affected by mesh distortions.Indeed,in all cases the convergence slopes of the Random mesh are close to those obtained via more regular meshes(structured and CVT).

4.2 Test case 2:h-analysis with the serendipity approach

In Subsection 2.4 we proposed the serendipity VEM approach to reduce the computational effort.We consider both the stingy and lazy choice and compare them with the standard VEM approach.

Figure 3 Test case 2:Comparison between VEMOand VEMLfor all meshes.

To make the following discussion clearer,we refer to the standard VEM approach as VEMO,to the stingy choice as VEMS,andfinally to the lazy choice as VEML.

Figure 4 Test case 2:Comparison between VEMOand VEMSfor all meshes.

We recall that,according to Remark 2.4,the two choices correspond to

with ηF=minimum number of straight lines necessary to cover the boundary of F.We focus on two aspects of the serendipity approach.On the one hand,we want to check that the serendipity procedure does not alter the accuracy.On the other hand we estimate the gain by counting the dofs with the standard and the serendipity approach.We count the number of vertex,edge and face dofs(in short,boundary dofs)for the standard VEM,#dof∂,and the serendipity VEM,#.We then define the gain as

We underline that we compute the gain only in terms of the boundary dofs,since the internal(volume)dofs can be removed by static condensation as forfinite elements.

Table 1 Test case 2:Gain of VEMSand VEMLover VEMOfor all meshes.

We show the convergence graphs of the lazy approach(VEML)in Figure 3,and of the stingy approach(VEMS)in Figure 4,together with the standard approach(VEMO).From these graphs we observe that the stingy choice is not so robust with respect to element degeneracies.Indeed,we recover the same convergence rates of the standard case for structured meshes,while the scheme fails to converge for CVT and random meshes,as shown in Figure 4:CVT fails for k=4 and random fails for

and

The lazy approach is definitely more robust,see Figure 3.For all the degrees k we recover the same convergence plots of VEMO(the convergence lines are indistinguishable from their counterpart of a standard VEMO).

In Table 1 we show the gain in terms of boundary dofs.Here,we can appreciate that the gain for the stingy choice is remarkable:For the case k=3 and 4,we save around the 40%of the face dofs.Consequently,if we are dealing with well-shaped meshes,the stingy serendipity approach can tear down the number of dofs.However,we also underline that the gain for the lazy choice is not as large as for the stingy case,but it is still noteworthy:It is at least the 25%for all the cases.

4.2.1 An adaptive stingy choice

Figure 5 Test case 3:Comparison between VEMOand VEMAfor CVT and Random meshes.

In this short paragraph we propose a strategy inspired by[9]to cure the stingy serendipity approach.The idea behind this method is to relax the conditions which determine the value of kFon a face F.Indeed,as explained in[9],the reason for the failures of the stingy choice is due to the presence of very small edges and/or edges laying almost on the same line.The strategy adopted in the code is the following:We fix an angle threshold,θ,and an edge ratio,ρ.Two edges forming an angle bigger than θ are considered as a single edge,and edges having length smaller than ρhFare neglected.If µFis the number of internal angles greater than θ,and ζF

is the number of edges of F with length

the definition of kFis modified as

Wefixed

and

and solved the same problem above.In the following graphs and tables we refer to this approach as VEMA.In Figure 5 we compare the convergence graphs of VEMAwith VEMO,while in Table 2 we collect the gain in terms of boundary degrees of freedom.We do not show the case of structured meshes since we get exactly the same results as the stingy choice.

Table 2 Test case 3:Gain for VEMS,VEMLand VEMAwith CVT and Random meshes.

We observe that this new way to compute kFis robust with respect to element degeneracies.Indeed,all the convergence lines provided by such method are undistinguishable from the standard VEM ones(see Figure 5).Moreover,the gain is now greater than that obtained with the lazy choice and close to the optimal one obtained with the stingy approach,see the highlighted values in Table 2.

[1]Ahmad,B.,Alsaedi,A.,Brezzi,F.,et al.,Equivalent projectors for virtual element methods,Comput.Math.Appl.,66(3),2013,376–391.

[2]Antonietti,P.F.,BeirÃo da Veiga,L.,Scacchi,S.and Verani,M.,A C1virtual element method for the Cahn-Hilliard equation with polygonal meshes,SIAM J.Numer.Anal.,54(1),2016,34–57.

[3]Artioli,E.,de Miranda,C.,Lovadina,C.and Patruno,L.,A stress/displacement virtual element method for plane elasticity problems,Computer Methods in Applied Mechanics and Engineering,325,2017,155–174.

[4]Ayuso,B.,Lipnikov,K.and Manzini,G.,The nonconforming virtual element method,ESAIM:M2AN,50(3),2016,879–904.

[5]BeirÃo da Veiga,L.,Brezzi,F.,Cangiani,A.,et al.,Basic principles of virtual element methods,Math.Models Methods Appl.Sci.,23(1),2013,199–214.

[6]BeirÃo da Veiga,L.,Brezzi,F.,Dassi,F.,et al.,Virtual element approximation of 2d magnetostatic problems,Comput.Methods Appl.Mech.Engrg.,327,2017,173–195.

[7]BeirÃo da Veiga,L.,Brezzi,F.,Marini,L.D.and Russo,A.,The hitchhiker’s guide to the virtual element method,Math.Models Methods Appl.Sci.,24(8),2014,1541–1573.

[8]BeirÃo da Veiga,L.,Brezzi,F.,Marini,L.D.and Russo,A.,Virtual element methods for general second order elliptic problems on polygonal meshes,Math.Models Methods Appl.Sci.,26(4),2016,729–750.

[9]BeirÃo da Veiga,L.,Brezzi,F.,Marini,L.D.and Russo,A.,Serendipity nodal VEM spaces,Comp.Fluids,141,2016,2–12.

[10]BeirÃo da Veiga,L.,Dassi,F.and Russo,A.,High-order virtual element method on polyhedral meshes,Computers&Mathematics with Applications,74(5),2017,1110–1122.

[11]BeirÃo da Veiga,L.,Lovadina,C.and Russo,A.,Stability analysis for the virtual element method,Math.Models Methods Appl.Sci.,27(13),2017,2557–2594.

[12]BeirÃo da Veiga,L.,Lovadina,C.and Vacca,G.,Divergence free Virtual Elements for the Stokes problem on polygonal meshes,ESAIM Math.Model.Numer.Anal.,51,2017,509–535.

[13]Benedetto,M.F.,Berrone,S.,Borio,A.,et al.,A hybrid mortar virtual element method for discrete fracture network simulations,Journal of Computational Physics,306,2016,148–166.

[14]Berrone,S.and Borio,A.,Orthogonal polynomials in badly shaped polygonal elements for the virtual element method,Finite Elem.Anal.Des.,129,2017,14–31.

[15]Brenner,S.,Guan,Q.and Sung,L.,Some estimates for virtual element methods,Computational Methods in Applied Mathematics,17,2017,553–574.

[16]Brezzi,F.and Marini,L.D.,Virtual element methods for plate bending problems,Comput.Methods Appl.Mech.Engrg.,253,2013,455–462.

[17]C´aceres,E.and Gatica,G.N.,A mixed virtual element method for the pseudostress-velocity formulation of the Stokes problem,IMA J.Numer.Anal.,37(1),2017,296–331.

[18]Cangiani,A.,Georgoulis,E.H.,Pryer,T.and Sutton,O.J.,A posteriori error estimates for the virtual element method,Numerische Mathematik,137,2017,857–893.

[19]Chi,H.,BeirÃo da Veiga,L.and Paulino,G.H.,Some basic formulations of the virtual element method(VEM)forfinite deformations,Computer Methods in Applied Mechanics and Engineering,318,2017,148–192.

[20]Ciarlet,P.G.,Thefinite element method for elliptic problems,Studies in Mathematics and Its Applications,4,North-Holland Publishing Co.,Amsterdam,New York,Oxford,1978.

[21]Du,Q.,Faber,V.and Gunzburger,M.,Centroidal voronoi tessellations:Applications and algorithms,SIAM Rev.,41(4),1999,637–676.

[22]Gain,A.L.,Paulino,G.H.,Leonardo,S.D.and Menezes,I.F.M.,Topology optimization using polytopes,Comput.Methods Appl.Mech.Engrg.,293,2015,411–430.

[23]Gain,A.L.,Talischi,C.and Paulino,G.H.,On the Virtual Element Method for three-dimensional linear elasticity problems on arbitrary polyhedral meshes,Comput.Methods Appl.Mech.Engrg.,282,2014,132–160.

[24]Mascotto,L.,A therapy for the ill-conditioning in the virtual element method,ArXiv:1705.10581.

[25]Mascotto,L.,BeirÃo da Veiga,L.,Chernov,A.and Russo,A.,Exponential convergence of the hp virtual element method with corner singularities,Numer.Math.,DOI:10.1007/s00211-017-0921-7.

[26]Mascotto,L.and Dassi,F.,Exploring high-order three dimensional virtual elements:Bases and stabilizations,2017,arXiv:1709.04371.

[27]Mora,D.,Rivera,G.and Rodr´ıguez,R.,A virtual element method for the Steklov eigenvalue problem,Math.Models Methods Appl.Sci.,25(8),2015,1421–1445.

[28]Rycroft,C.H.,Voro++:A three-dimensional voronoi cell library in C++,Chaos,19(4),2009,041111.

[29]Schatz,A.H.,An observation concerning Ritz-Galerkin methods with indefinite bilinear forms,Math.Comp.,28,1974,959–962.