High order VEM on curved domains.

# High order VEM on curved domains.

Silvia Bertoluzza IMATI “E. Magenes”, CNR, Pavia (Italy) Micol Pennacchio IMATI “E. Magenes”, CNR, Pavia (Italy)  and  Daniele Prada IMATI “E. Magenes”, CNR, Pavia (Italy)
July 22, 2019
###### Abstract.

We deal with the virtual element method (VEM) for solving the Poisson equation on a domain with curved boundaries. Given a polygonal approximation of the domain , the standard order VEM [hitchVEM], for increasing, leads to a suboptimal convergence rate. We adapt the approach of [BDT] to VEM and we prove that an optimal convergence rate can be achieved by using a suitable correction depending on high order normal derivatives of the discrete solution at the boundary edges of , which, to retain computability, is evaluated after applying the projector onto the space of polynomials. Numerical experiments confirm the theory.

###### Key words and phrases:
Virtual Element method, Nitsche’s method, curved domain
###### 1991 Mathematics Subject Classification:
This paper has been realized in the framework of ERC Project CHANGE, which has received funding from the European Research Council (ERC) under the European Unionâs Horizon 2020 research and innovation programme (grant agreement No 694515), and it was partially supported by INdAM - GNCS

## 1. Introduction

The virtual element method (VEM) is a PDE discretization framework designed to easily handle meshes consisting of very general polygonal or polyhedral elements [basicVEM]. The method can be considered as a generalization of the Finite Element Method (FEM) to polytopal tessellations, in that it looks for the solution in a conforming discretization space with a Galerkin approach. By giving up conformity in the discretization of the bilinear form corresponding to the differential operator, the method manages to avoid the explicit construction of the basis functions (whence the name virtual). Everything is computed directly in terms of the degrees of freedom by resorting to suitable “computable” (in terms of the degrees of freedom) elementwise projectors onto the space of polynomials (see [hitchVEM]), ultimately allowing to define a discrete bilinear form satisfying polynomial exactness and stability properties which allow to prove optimal error estimate for discretization of (arbitrary) order . Different model problems have already been tackled by using VEM ([beirao_parab, Antonietti_VEM_Stokes, Antonietti_VEM_Cahn, beirao_stokes, beirao_Navier_Stokes, perugia_Helmholtz, beirao_elastic, beirao_linear_elasticity, VEM_3D_elasticity, VEM_mixed, VEM_discrete_fracture, VEM_Laplace_Beltrami]), and, while most of the literature deals with the version of the method, the and versions were also discussed and analyzed ([antonietti_p_VEM, beirao_hp, beirao_hp_exponential].

In this paper we consider the problem of extending the method to problems in domains with smooth curved boundaries. As it happens in the Finite Element case, the approximation of the curved domain by straight facets introduces an error that, for higher order methods, can dominate the analysis. Different approaches for the accurate treatment of curved domains in the finite element framework can be found in literature, see e.g. [FEMsurvey_curved]. Among the different possible approaches to such a problem, following the guidelines of [BDT], we choose here to approximate the curved domain with a polygonal domain , while compensating for the discrepancy in the geometry by suitably modifying the bilinear form. Complying with the VEM philosophy, the modified bilinear form will retain the property of being computable in terms of the degrees of freedom. Of course, other approaches are possible. In [VEM_curved_beirao] the authors propose a direct definition of a modified virtual element space that accommodates curved elements, whose boundary matches exactly the boundary of . While loosing exact reproduction of polynomials, the resulting method retains optimality and, contrary to the one we propose here, it is immediately well suited to deal with curved interior interfaces. On the other hand, our approach has the advantage of only requiring, for the implementation, minor modifications with respect to the polygonal case. To the best of our knowledge, other numerical methods that can handle curved polytopal meshes are only [Mimetic_curved] and [HDG_curved].

The paper focuses on a simple elliptic model problem in 2D and it is organized as follows. As the projection method of [BDT] combines Nitsche’s technique for imposing non homogeneous boundary conditions [Nitsche] with the improved accuracy polygonal domain approximation of [Thomee], we start, in Section 2, to adapt Nitsche’s method to the Virtual Element framework and we provide a theoretical analysis of the resulting discretization, proving stability and an error estimate. In Section 3 we introduce and analyze the discretization for the problem on curved domains, proving also in this case stability and optimal error estimate. Finally, in Section 4, we test the method on several test cases, with method of different order. Throughout the paper, we will use the notation (resp. ) to signify that the quantity is bounded from above (resp. from below) by a constant times the quantity , with independent of the mesh size parameter , the diameter the specific shape of the polygon , but possibly depending on the polynomial order of the method, and on the shape regularity constant and appearing in Assumption 2.1.

## 2. The Nitsche’s method in the Virtual Element Context

Before considering the problem of solving a PDE on a domain with a curved boundary, let us discuss how Nitsche’s method for imposing non homogeneous boundary condition [Nitsche] can be applied in the context of the virtual element method. Throughout this section let then denote a bounded polygonal domain. To fix the ideas, we consider the following simple model problem:

 (2.1) −Δu=f, in Ω,u=g, on ∂Ω,

with and . Assume that we are given a family of quasi uniform tessellations of into polygonal elements of diameter . We make the following standard regularity assumptions on the polygons of the tessellation:

###### Assumption 2.1.

There exists constants such that:

1. each element is star-shaped with respect to a ball of radius ;

2. for each element in the distance between any two vertices of is .

Under Assumptions 2.1, several bounds hold uniformly in [Lipnikov]. In particular, in the following we will make use of an inverse inequality on the space of polynomials of order less than or equal to : for all and for all with it holds that

 (2.2) ∥p∥k,K≲hj−kK∥p∥j,K.

Moreover we will make use of the following trace inequality: for all we have

 (2.3) ∥ϕ∥0,∂K≲h−1/2K∥ϕ∥0,K+h1/2K|ϕ|1,K.

We will consider the standard order Virtual Element discretization space ([basicVEM]), whose definition we briefly recall. For each polygon we let the space be defined as

 Bm(∂K)={v∈C0(∂K):v|e∈Pm ∀e∈EK},

where denotes the set of edges of the polygon . We introduce the local VE space as:

 VK,m={v∈H1(K): v|∂K∈Bm(∂K), Δv∈Pm−2(K)}

(with ). The global discrete VE space is then defined as

 (2.4) Vh ={v∈H1(Ω):v|K∈VK,m ∀K∈Th}= ={v∈H1(Ω):∀K∈Th v|∂K∈Bm(∂K), Δv|K∈Pm−2(K)}.

A function in is uniquely determined by the following degrees of freedom

• its values at the vertices of the tessellation;

• (only for ) for each edge , its values at the internal points of the -points Gauss-Lobatto quadrature rule on ;

• (only for ) for each element , its moments in up to order .

For any given function we can then define the unique function such that: a) the values of and at the vertices of the tessellation coincide; b) for each edge , the values of and at the internal points of the -points Gauss-Lobatto quadrature rule on coincide; c) for each element , the moments up to order of and in coincide. The function satisfies the following local approximation bound [basicVEM]: if , with , then

 (2.5) ∥w−wI∥0,K+hK|w−wI|1,K≲hsK|w|s,K.

Let now

 a(ϕ,ψ)=∫Ω∇ϕ⋅∇ψ,aK(ϕ,ψ)=∫K∇ϕ⋅∇ψ,

and let denote the projection defined by

 aK(Π∇Kϕ,p)=aK(ϕ,p),∀p∈Pm,∫KΠ∇Kϕ=∫Kϕ.

We recall that for , can be computed directly from the values of the degrees of freedom ([hitchVEM]) without the need of explicitly constructing it (which would imply somehow solving a partial differential equation), by taking advantage of the identity

 ∫K∇wh⋅∇p=−∫KwhΔp+∫∂Kwh∂p∂νK,

that allows to express the term of the left hand side in terms of the interior moments of (for , is a polynomial of degree less than or equal to ) and of an integral on the boundary (where is a known piecewise polynomial). Letting denote the space of discontinuous piecewise polynomials of order less than or equal to

 P∗m={ϕ∈L2(Ω):ϕ|K∈Pm ∀K∈Th},

we let be defined by assembling, element by element, the ’s:

 Π∇ϕ|K=Π∇K(ϕ|K), ∀K∈Th.

The discretization of (2.1) by the Nitsche’s method would consist in looking for such that for all one has

 a(uh,vh)−∫∂Ωνuhvh−∫∂Ωuhνvh+γh−1∫∂Ωuhvh=∫Ωfvh−∫∂Ωgνvh+γh−1∫∂Ωgvh,

where stands for , denoting the outer normal to . As typical for the Virtual Element method, both at the right hand side and at the left hand side of such an equation we find terms which are not “computable”, that is that can not be computed exactly with only the knowledge of the value of the degrees of freedom of and . Besides the bilinear form , which can be treated by the standard approach, this is the case for all the terms involving and . As usually done in VEM, the bilinear form is then replaced with

 ah(uh,vh)=KaKh(uh,vh),

where

 aKh(ϕ,ψ)=aK(Π∇Kϕ,Π∇Kψ)+SKa(ϕ−Π∇Kϕ,ψ−Π∇Kψ).

We recall that different choices are possible for the bilinear form (see [beirao_stab]), the essential requirement being that it satisfies

 aK(ϕ,ϕ)≲SKa(ϕ,ϕ)≲aK(ϕ,ϕ),∀ϕ∈VK,m  with Π∇Kϕ=0,

so that the local discrete bilinear forms satisfy the following two properties:

• Stability:

 aK(ϕ,ϕ)≲aKh(ϕ,ϕ)≲aK(ϕ,ϕ),∀ϕ∈VK,m
• m-consistency: for any and

 (2.6) aKh(ϕ,p)=aK(ϕ,p).

In the numerical tests performed in Section 4 we made the standard choice of defining in terms of the vectors of local degrees of freedom as the properly scaled euclidean scalar product.

As far as the terms involving the normal derivative are concerned, we treat them by replacing and , boundary edge by boundary edge, respectively with and . Then we can write the Nitsche’s method for the VEM discretization of 3.1 as: find such that for all it holds that

 (2.7) ah(uh,vh)−e∈E∫eνΠ∇(uh)vh−e∈E∫eνΠ∇(vh)uh+γh−1∫∂Ωuhvh=∫Ωfvh−e∈E∫eg(νΠ∇(vh)−γh−1vh),

where is a positive constant and denotes the set of edges of lying on .

We introduce the norm:

 (2.8) \vvvertϕ\vvvert2Ω=|ϕ|21,Ω+h−1∥ϕ∥20,∂Ω

and the space defined as the closure of with respect to the norm . Setting

 (2.9) Bh,γ(ϕ,ψ)=ah(ϕ,ψ)−e∈E∫eνΠ∇(ϕ)ψ−e∈E∫eνΠ∇(ψ)ϕ+γh−1∫∂Ωϕψ,

we start by proving the following lemma:

###### Lemma 2.2.

For all , we have

 (2.10) |Bh,γ(ϕ,ψ)|≲\vvvertϕ\vvvertΩ\vvvertψ\vvvertΩ.

Moreover, there exists such that, for , the bilinear form verifies for all

 (2.11) Bh,γ(ϕ,ϕ)≳\vvvertϕ\vvvert2Ω,

(the implicit constant in the two inequalities depending on ).

###### Proof.

We observe that, since is a polynomial, it is not difficult to verify that the following inverse bound holds:

 (2.12) ∥νΠ∇(ϕ)∥0,e≲h−1/2|Π∇(ϕ)|1,Ke≤h−1/2|ϕ|1,Ke,

where, for , is the unique polygon of the tessellation having as an edge. Then we have

 e∈E∫eνΠ∇(ϕ)ψ≲e∈E∥νΠ∇(ϕ)∥0,e∥ψ∥0,e≲e∈E|ϕ|1,Keh−1/2∥ψ∥0,e.

Obtaining (2.10) is then not difficult. As far as (2.11) is concerned, we have

 (2.13) Bh,γ(ϕ,ϕ)=ah(ϕ,ϕ)+γh−1∥ϕ∥20,∂Ω−2⟨νΠ∇(ϕ),ϕ⟩

where we use the notation

 ⟨ϕ,ψ⟩=e∈E∫eϕψ.

We now have the following bounds

 ⟨νΠ∇(ϕ),ϕ⟩≲e∈E∥νΠ∇(ϕ)∥0,e∥ϕ∥0,e

Thanks to the inverse inequality (2.12), we can write, for arbitrary,

 ⟨νΠ∇(ϕ),ϕ⟩≲e∈Eh−1/2|ϕ|1,Ke∥ϕ∥0,e≲ε2|ϕ|21,Ω+12εh−1∥ϕ∥20,∂Ω.

Substituting into (2.13) we obtain, for a fixed positive constant independent of

 Bh,γ(ϕ,ϕ)≳(1−c1ε)|ϕ|21,Ω+(γ−c1ε)h−1∥ϕ∥20,∂Ω.

We now choose and if with chosen in such a way that , the thesis easily follows. ∎

Existence and uniqueness of the solution of (2.7) easily follow.

We are then able to prove the following result:

###### Theorem 2.3.

If , with , and if we chose , with given by Lemma 2.2, then the following error estimate holds

 \vvvertu−uh\vvvertΩ≲hs−1|u|s,Ω.
###### Proof.

Let denote the VEM interpolant and the projection of onto the space of discontinuous piecewise polynomials. For any , with we have, for with

 (2.14) ∥w−wπ∥j,K≲hk−jK|w|k,K.

For the bound (2.14) can be proven by a standard argument combining the bound for with an inverse inequality. Moreover we have, for ,

 (2.15) ∥ν(w−wπ)∥0,e≲h−1/2|w−wπ|1,K+h1/2|w−wπ|2,K≲hs−3/2|w|s,K.

We set . Using the definition (2.9), summing and subtracting , using the -consistency (2.6) to replace with and then, summing and subtracting , we can write

 (2.16) \vvvertuI−uh\vvvert2Ω≲Bh,γ(uI,dh)−Bh,γ(uh,dh)==K∈ThaKh(uI−uπ,dh)+K∈ThaK(uπ−u,dh)+a(u,dh)−⟨νu,dh⟩+⟨ν(u−Π∇(uI)),dh⟩−⟨u,νΠ∇(dh)⟩+⟨u−uI,νΠ∇(dh)⟩+γh−1⟨u,dh⟩+γh−1⟨uI−u,dh⟩−∫Ωfdh+⟨g,νΠ∇(dh)−γh−1dh⟩=E1+E2+E3+E4+E5.

with

 E1=K∈ThaKh(uI−uπ,dh),E2=K∈ThaK(uπ−u,dh),E3=⟨ν(u−Π∇(uI)),dh⟩ E4=⟨u−uI,νΠ∇(dh)⟩,E5=γh−1⟨uI−u,dh⟩,

where we used that, as is the solution of (2.1),

 a(u,dh)−⟨νu,dh⟩−∫Ωfdh=0,⟨u,νΠ∇(dh)⟩−γh−1⟨u,dh⟩−⟨g,νΠ∇(dh)−γh−1dh⟩=0.

Let us then bound the different components of the error.

Both and are standardly encountered in the analysis of the VEM method, and a bound can be found in the literature (see e.g. [basicVEM]), yielding

 E1≲|dh|1,Ωhs−1|u|s,Ω, and E2≲|dh|1,Ωhs−1|u|s,Ω.

On the other hand we have

 E3≲e∈E∥ν(u−Π∇(uI))∥0,e∥dh∥0,e.

Now, using (2.15), we have

 ∥ν(u−Π∇uI)∥0,e≤∥ν(u−uπ)∥0,e+∥νΠ∇(uπ−uI)∥0,e≲hs−3/2|u|s,K+h−1/2(|uπ−uI|1,K+|uI−u|1,K)

yielding

 E3≲hs−1h−1/2|u|s,Ω∥dh∥0,∂Ω.

As far as is concerned we have

 E4≤e∈E∥u−uI∥0,e∥νΠ∇(dh)∥0,e≲e∈Ehs−1/2|u|s,Keh−1/2|dh|1,Ke≲hs−1|u|s,Ω|dh|1,Ω.

A similar argument yields

 E5≲hs−1|u|s,Ωh−1/2∥dh∥0,∂Ω,

finally giving

 \vvvertuI−uh\vvvert2Ω≲hs−1|u|s,Ω\vvvertuI−uh\vvvertΩ.

Dividing both sides by and using a triangular inequality we get the thesis. ∎

## 3. The Virtual Element Method on domains with curved boundary

Let now consider the solution of the same model problem

 (3.1) −Δu=f, in Ω,u=g, on ∂Ω

with, once again, , , where now is a convex domain with curved boundary assumed, for the sake of convenience, to be of class . In order to solve such a problem by the Virtual Element method, we assume that is approximated by a family of polygonal domains , , each endowed with a quasi uniform shape regular tessellation into polygons with diameter . We assume that all the vertices of lying on also lie on . As is convex this implies that .

### 3.1. The Projection Method of Bramble, Dupont and Thomée

In order to deal with the curved boundary, following [BDT] we apply on the approximating domain a modified version of Nitsche’s method, which takes into account that the boundary data is given on rather than on .

Letting be the outer normal to , for we let denote the non negative scalar such that

 x+δ(x)νh(x)∈∂Ω.

It is known ([BDT]) that, as is smooth and convex, we have that

 (3.2) δh=supx∈∂Ωhδ(x)=o(h2).

We let be the order VEM discretization space relative to the tessellation (defined by (2.4)). Setting , the projection method of [BDT] reads as follows: find such that for all it holds that

 (3.3) Bh,γ(uh,vh)−e∈E∫e(kj=1δjj!jνhΠ∇(uh))(νhΠ∇(vh)−γh−1vh)=∫Ωhfvh−e∈E∫eg(νhΠ∇vh−γh−1vh)

where denotes the -th partial derivative of in the direction and where, for ,

 g(x)=g(x+δ(x)νh(x)).

We then introduce the bilinear form defined by

 (3.4) Ah,γ(ϕ,ψ)=Bh,γ(ϕ,ψ)−e∈E∫e(kj=1δjj!jνhΠ∇(uh))(νhΠ∇(vh)−γh−1vh),

and prove the following lemma:

###### Lemma 3.1.

For all , in it holds that

 (3.5) |Ah,γ(ϕ,ψ)≲\vvvertϕ\vvvertΩh\vvvertψ\vvvertΩh.

Moreover, there exists such that for all the bilinear form verifies for all

 (3.6) Ah,γ(ϕ,ϕ)≳\vvvertϕ\vvvert2Ωh

provided with .

###### Proof.

Let be defined as

 Ck(ϕ,ψ)=e∈E∫e(kj=1δjj!jνhϕ)ψ

so that we can write

 Ah,γ(ϕ,ϕ)=Bh,γ(ϕ,ψ)−Ck(Π∇(ϕ),νhΠ∇(ψ)−γh−1ψ).

As for all , is a piecewise polynomial, using an inverse inequality and the continuity of the operator we get that for and we have

 (3.7) ∥jνhΠ∇K(ϕ)∥0,e≲h1/2−j|Π∇K(ϕ)|1,Ke≲h1/2−j|ϕ|1,Ke,

where, once again, is the unique element of having as an edge. Since , using (3.7) we have

 (3.8) |Ck(Π∇(ϕ),ψ)|≲e∈Ekj=1(δhh)jhj∥jνhΠ∇(ϕ)∥0,e∥ψ∥0,e≲δhh|ϕ|1,Ωhh1/2∥ψ∥0,e

as well as

 (3.9) |Ck(Π∇(ϕ),νhΠ∇(ψ))|≲e∈Ekj=1(δhh)jhj∥jνhΠ∇(ϕ)∥0,e∥νhΠ∇(ψ)∥0,e≲δhh|ϕ|1,Ωh|ψ|1,Ωh.

Combining with (2.10) the bound (3.5) easily follows. Let us now consider (3.6). Combining (2.11) with (3.8) and (3.9) we obtain, for arbitrary and , and fixed positive constants,

 Ah,γ(ϕ,ϕ)=Bh,γ(ϕ,ϕ)−Ck(Π∇(ϕ),νhΠ∇(ϕ))+γh−1Ck(Π∇(ϕ),ϕ)≳(1−c1ε−c2δhh−c3δhhγ)|ϕ|21,Ωh+(γ−c1ε−c3γδhh)h−1∥ϕ∥20,∂Ωh.

We now choose and we fix in such a way that . For , set now . As , we can choose in such a way that for all , and . The thesis easily follows. ∎

Once again, existence and uniqueness of the solution of (3.3) easily follow. Moreover, the error estimate given by the following Theorem hold.

###### Theorem 3.2.

If , with , for and ( and given by Lemma 3.1) the following error estimate holds

 \vvvertu−uh\vvvertΩh≲hs−1|u|s,Ω+h−1/2δk+1∥u∥k+1,∞,Ω.
###### Proof.

Let denote the usual VEM interpolant. Setting and proceeding as in Theorem 2.3, we have

 (3.10) \vvvertuI−uh\vvvert2Ωh≲Ah,γ(uI,dh)−Ah,γ(uh,dh)=K∈ThaKh(uI−uπ,dh)+KaK(uπ−u,dh)+a(u,dh)−⟨νhu,dh⟩+⟨νh(u−Π∇(uI)),dh⟩−⟨u,νhΠ∇(dh)⟩+⟨u−uI,νhΠ∇(dh)⟩+γh−1⟨u,dh⟩+γh−1⟨uI−u,dh⟩−Ck(u,νhΠ∇(dh)−γh−1dh)+Ck(u−Π∇(uI),νhΠ∇(dh)−γh−1dh)−∫Ωhfdh+⟨g,νhΠ∇(dh)−γh−1dh⟩=E1+E2+E3+E4+E5+E6+E7.

with , , , , and as in the proof of Theorem 2.3 and with

 E6=Ck(u−Π∇(uI),νhΠ∇(dh))−γh−1dh),E7=⟨g−kj=0δjj!jνhu,νhΠ∇(dh)−γh−1dh⟩.

The components can be bounded exactly as in Theorem 2.3. Let us then bound the last two terms and . Since , we have

 E6≲e∈Ekj=1hj∥jνh(u−Π∇(uI))∥0,e(∥νhΠ∇(dh)∥0,e+h−1∥dh∥0,e).

Now, using (2.3) as well as (3.7) we have

 ∥jνh(u−Π∇(uI))∥0,e≲∥jνh(u−uπ)∥0,e+∥jνh(Π∇(uπ−uI))∥0,e≲h−1/2|u−uπ|j,Ke+h1/2|u−uπ|j+1,Ke+h1/2−j|uπ−uI|1,Ke≲hs−j−1/2|u|s,Ke,

where is bound by adding and subtracting and using the VEM approximation bounds (2.5) and (2.14), yielding

 E6≲hs−1|u|s,Ω\vvvertdh\vvvertΩh.

Finally, takes into account the approximation of the curved boundary by projection, and, following the paper by Thomée, it can be bound as

 E7≲e∈E∥g−kj=0δjj!jνhu∥0,e(∥νhΠ∇(dh)∥0,e+h−1∥dh∥0,e)≲h−1/2δk+1∥u∥k+1,∞,Ωh\vvvertdh\vvvertΩh.

Assembling the bounds for the seven terms we finally obtain

 \vvvertδh\vvvert2Ωh≲hs−1|u|s,Ωh\vvvertδh\vvvertΩh+h−1/2δk+1∥u∥k+1,∞,Ω\vvvertδh\vvvertΩh.

Dividing both sides by and using a triangular inequality we get the thesis. ∎

As