Table Of ContentVirtual Element Method for fourth order problems:
L2−estimates
Claudia Chinosia, L. Donatella Marinib
6
1
aDipartimento di Scienze e Innovazione Tecnologica, Universit`a del Piemonte Orientale,
0
Viale Teresa Michel 11, 15100, Alessandria, Italy email: [email protected]
2
n bDipartimento di Matematica, Universita` di Pavia, and IMATI-CNR, Via Ferrata 1,
a 27100, Pavia, Italy email: [email protected]
J
7
2
] Abstract
A
N
WeanalysethefamilyofC1-VirtualElementsintroducedin[7]forfourth-order
.
h problems and prove optimal estimates in L2 and in H1 via classical duality
t
a arguments.
m
Keywords: Fourth-order problems problems; C1-Virtual Elements
[
1
v
4
1. Introduction
8
4
7
TheVirtualElementMethod(VEM),introducedin[2]andfurtherdeveloped
0
in [1], can be seen as the extension of the Finite Element Method (FEM) to
.
1
decompositions into almost arbitrary polygons and/or polyhedra. Since the
0
first paper in 2013 ([2]) the Virtual Element approach has been applied to
6
1 a number of applications: linear elasticity in two and three dimensions ([3]
: and [9], respectively), general advection-diffusion-reaction problems, both in
v
primal [5] and in mixed form [4], Helmholtz problem [12], and plate bending
i
X problems in the Kirchhoff-Love formulation [7]. In [7] a family of elements was
r constructed and analysed, showing the ductility of the approach to design C1-
a
elements. Optimal convergence rates were proved in the energy norm, i.e., in
H2. Namely, order k−1, with k ≥ 2 whenever the discrete space V contains
h
locallypolynomialsofdegreek. Inthepresentpaperweproveoptimalestimates
also in H1 and in L2, obtained via classical duality arguments, and we provide
numerical results confirming the theoretical estimates.
We point out that the use of C1-approximations is of interest not only for
plate bending problems, although the family of elements we are dealing with
was originally introduced having in mind plates. In many other applications
Preprint submitted to Elsevier January 28, 2016
the presence of fourth order operators calls for higher continuity. For example,
Cahn-Hilliard equation for phase separation, or Navier-Stokes equations in the
stream-vorticityformulationcontainthebiharmonicoperator,exactlyasinplate
bending problems, which we will refer to throughout the paper.
An outline of the paper is as follows. In Section 2 we state the continuous
problem and fix some notation. In Section 3 we recall the VEM-approximation
and the convergence result given in [7]. In particular, in Subsection 3.3 we
propose a different approximation of the loading term, more suited for deriving
optimalestimatesinL2 andH1. InSections4and5weproveerrorestimatesin
H1 and in L2, respectively. Numerical results are presented in Section 6, and a
comparisonwiththeclassicalClough-TocherandReduced-Clough-Tocherfinite
elements is carried out.
Throughout the paper we shall use the common notation for the Sobolev
spaces Hm(D) for m a non-negative integer and D an open bounded domain.
In particular (see e.g. [11], [8]) the L2(D) scalar product and norm will be
indicated by (·,·) or (·,·) and (cid:107)·(cid:107) or (cid:107)·(cid:107) , respectively. When D ≡Ω
0,D D 0,D D
the subscript D will often be omitted. Finally, P will denote the space of
k
polynomials of degree ≤ k, with the convention that P = {0}, and C will
−1
denote a positive constant independent of the mesh size h.
2. The Continuous Problem
Let Ω ⊂ R2 be a convex polygonal domain occupied by the plate, let Γ be
its boundary, and let f ∈L2(Ω) be a transversal load acting on the plate. The
Kirchoff-Love model for thin plates (see e.g. [8]) corresponds to look for the
transversal displacement w, the solution of
D∆2w = f in Ω, (2.1)
where D =Et3/12(1−ν2) is the bending rigidity, t the thickness, E the Young
modulus, and ν the Poisson’s ratio. Assuming for instance the plate to be
clamped all over the boundary, equation (2.1) is supplemented with the bound-
ary conditions
∂w
w = =0 on Γ. (2.2)
∂n
The variational formulation of (2.1)-(2.2) is:
(cid:40)
Find w ∈V :=H2(Ω) solution of
0
(2.3)
a(w,v)=(f,v) ∀v ∈H2(Ω),
0 0
where the energy bilinear form a(·,·) is given by
(cid:104) (cid:90) (cid:90) (cid:105)
a(w,v)=D (1−ν) w v dx+ν ∆w∆vdx . (2.4)
/ij /ij
Ω Ω
2
In (2.4) v = ∂v/∂x , i = 1,2, and we used the summation convention of
/i i
repeated indices. Setting (cid:107)v(cid:107) := |v| , it is easy to see that, thanks to the
V 2,Ω
boundary conditions in V and to the Poincar´e inequality, this is indeed a norm
on V. Moreover
∃M >0 such that a(u,v)≤M(cid:107)u(cid:107) (cid:107)v(cid:107) u,v ∈V, (2.5)
V V
∃α>0 such that a(v,v)≥α(cid:107)v(cid:107)2 v ∈V. (2.6)
V
Hence, (2.3) has a unique solution, and (see, e.g. [11])
(cid:107)w(cid:107) ≤C(cid:107)f(cid:107) . (2.7)
V L2(Ω)
3. Virtual Element discretization
We recall the construction of the family of elements given in [7], and the
estimates there obtained. The family of elements depends on three integer
indices (r, s, m), related to the degree of accuracy k ≥2 by:
r =max{3,k}, s=k−1, m=k−4. (3.1)
Let T be a decomposition of Ω into polygons K, and let E be the set of edges
h h
inT . Wedenotebyh thediameterofK,i.e.,themaximumdistancebetween
h K
any two vertices of K. On T we make the following assumptions (see e.g. [2]):
h
H1 there exists a fixed number ρ > 0, independent of T , such that for every
0 h
element K (with diameter h ) it holds
K
i) K is star-shaped with respect to all the points of a ball of radius ρ h , and
0 K
ii) every edge e of K has length |e|≥ρ h .
0 K
3.1. Definition of the discrete space V
h
On a generic polygon K we define the local virtual element space as
∂v
V(K):={v ∈H2(K): v ∈P (e), ∈P (e) ∀e∈∂K, ∆2v ∈P (K)}.
|e r ∂n|e s m
Then the global space V is given by
h
V ={v ∈V :v ∈V(K), ∀K ∈T }. (3.2)
h |K h
A function in V is uniquely identified by the following degrees of freedom:
h
• The value of v(ξ) ∀ vertex ξ (3.3)
3
• The values of v (ξ) and v (ξ) ∀ vertex ξ (3.4)
/1 /2
(cid:90)
• For r >3, the moments q(ξ)v(ξ)dξ ∀q ∈P (e), ∀e∈E (3.5)
k−4 h
e
(cid:90)
• For s>1, the moments q(ξ)v (ξ)dξ ∀q ∈P (e), ∀e∈E (3.6)
/n s−2 h
e
(cid:90)
• For m≥0, the moments q(x)v(x)dx ∀q ∈P (K) ∀K. (3.7)
m
K
Proposition 3.1. In each element K the d.o.f. (3.3)–(3.7) are unisolvent.
Moreover, (3.3), (3.4), and (3.5) uniquely determine a polynomial of degree ≤r
on each edge of K, the degrees of freedom (3.4) and (3.6) uniquely determine a
polynomial of degree ≤s on each edge of K, and the d.o.f. (3.7) are equivalent
to prescribe Π0 v in K, where, for m a nonnegative integer,
m
Π0 v is the L2(K)−projector operator onto the space P (K). (3.8)
m m
Remark 3.1. We recall that our assumptions on T allow to define, for every
h
smooth enough w, an “interpolant” in V with the right interpolation properties.
h
More precisely, if g (w), i = 1,2,...G are the global d.o.f. in V , there exists a
i h
unique element w ∈V such that
I h
g (w−w )=0 ∀i=1,2,....G. (3.9)
i I
Moreover, by the usual Bramble-Hilbert technique (see e.g. [8]) and scaling
arguments (see e.g. [6]) we can prove that
(cid:107)w−w (cid:107) ≤ Chβ−s|w| s=0,1,2, 3≤β ≤k+1 (3.10)
I s,Ω β,Ω
(with C >0 independent of h) as in the usual Finite Element framework.
3.2. Construction of a
h
We need to define a symmetric discrete bilinear form which is stable and
consistent. More precisely, denoting by aK(·,·) the restriction of a (·,·) to a
h h
generic element K, the following properties must be satisfied (see [2]). For all
h, and for all K in T ,
h
• k-Consistency: ∀p ∈P , ∀v ∈V(K)
k k
aK(p ,v)=aK(p ,v). (3.11)
h k k
• Stability: ∃ two positive constants α and α∗, independent of h and of K,
∗
such that
∀v ∈V(K) α aK(v,v)≤aK(v,v)≤α∗aK(v,v). (3.12)
∗ h
4
The symmetry of a , (3.12) and the continuity of aK imply the continuity of
h
aK:
h
(cid:16) (cid:17)1/2(cid:16) (cid:17)1/2
aK(u,v)≤ aK(u,u) aK(v,v)
h h h (3.13)
≤ α∗M(cid:107)u(cid:107) (cid:107)v(cid:107) for all u and v in V(K).
2,K 2,K
In turn, (3.12) and (3.13) easily imply
∀v ∈V α a(v,v)≤a (v,v)≤α∗a(v,v), (3.14)
h ∗ h
and
a (u,v)≤ α∗M(cid:107)u(cid:107) (cid:107)v(cid:107) for all u and v in V . (3.15)
h V V h
Next,letΠK :V(K)−→P (K)⊂V(K)betheoperatordefinedasthesolution
k k
of
aK(ΠKψ,q)=aK(ψ,q) ∀ψ ∈V(K), ∀q ∈P (K)
k k
(cid:90) (cid:90) (3.16)
(ΠKk ψ−ψ)=0, ∇(ΠKk ψ−ψ)=0.
∂K ∂K
We note that for v ∈P (K) the first equation in (3.16) implies (ΠKv) =v
k k /ij /ij
for i,j =1,2, that joined with the second equation gives easily
ΠKv =v ∀v ∈P (K). (3.17)
k k
Hence, ΠK is a projector operator onto P (K). Let then SK(u,v) be a sym-
k k
metric positive definite bilinear form, verifying
c aK(v,v)≤SK(v,v)≤c aK(v,v), ∀v ∈V(K) with ΠKv =0, (3.18)
0 1 k
for some positiveconstants c , c independent of K and h . Werefer to [7]for
0 1 K
a precise choice of SK(u,v). We just recall that SK(u,v) can simply be taken
as the euclidean scalar product associated to the degrees of freedom, properly
scaled to satisfy (3.18). Then set
aK(u,v):=aK(ΠKu,ΠKv)+SK(u−ΠKu,v−ΠKv). (3.19)
h k k k k
Clearly the bilinear form (3.19) verifies both the consistency property (3.11)
and the stability property (3.12).
We postpone the construction of the right-hand side, and recall the conver-
gence result of [7].
Theorem 3.1. Under assumptions H1 on the decomposition the discrete prob-
lem:
(cid:40)
Find w ∈V solution of
h h
(3.20)
a (w ,v )=<f ,v > ∀v ∈V
h h h h h h h
has a unique solution w . Moreover, for every approximation w of w in V
h I h
and for every approximation w of w that is piecewise in P , we have
π k
(cid:16) (cid:17)
(cid:107)w−w (cid:107) ≤C (cid:107)w−w (cid:107) +(cid:107)w−w (cid:107) +(cid:107)f −f (cid:107)
h V I V π h,V h V(cid:48)
h
5
where C is a constant depending only on α, α , α∗, M and, with the usual
∗
notation, the norm in V(cid:48) is defined as
h
<f −f ,v >
(cid:107)f −f (cid:107) := sup h h . (3.21)
h Vh(cid:48) vh∈Vh (cid:107)vh(cid:107)V
3.3. Construction of the right-hand side
In order to build the loading term < f ,v > for v ∈ V in a simple and
h h h h
easy way it is convenient to have internal degrees of freedom in V , and this
h
means,accordingto(3.1)and(3.7),thatk ≥4isneeded. In[7]suitablechoices
were made for different values of k, enough to guarantee the proper order of
convergence in H2. Namely,
(cid:107)w−w (cid:107) ≤Chk−1(cid:107)w(cid:107) . (3.22)
h V k+1
In order to derive optimal estimates in L2 and H1 we need to make different
choices. To this end, following [1], we modify the definition (3.2) of V . For
h
k ≥ 2, and r and s related to k by (3.1), let Wk(K) be the new local space,
given by
∂v
Wk(K):={v ∈H2(K): v ∈P (e), ∈P (e)∀e∈∂K, ∆2v ∈P (K)}.
|e r ∂n|e s k−2
For k =2 we define the new global space as
(cid:90) (cid:90)
W2 ={v ∈V :v ∈W2(K),and v dx= ΠKv dx ∀K ∈T }, (3.23)
h |K k h
K K
and for k ≥3
Wk ={v ∈V :v ∈Wk(K), and
h |K
(3.24)
(cid:82) vp dx=(cid:82) ΠKvp dx, α=k−3,k−2 ∀K ∈T }.
K α K k α h
In (3.24) p are homogeneous polynomials of degree α. It can be checked that
α
the d.o.f. (3.3)–(3.7) are the same, but the added conditions on the moments
allownowtocomputetheL2−projectionofanyv ∈Wk ontoP (K)∀K,and
h k−2
not only onto P (K) as before. Taking then f = L2−projection of f onto
k−4 h
the space of piecewise polynomials of degree k−2, that is,
f =Π0 f on each K ∈T ,
h k−2 h
the right-hand side in (3.20) can be exactly computed:
(cid:90) (cid:90)
(cid:88) (cid:88)
<f ,v >= f v dx≡ (Π0 f)v dx
h h h h k−2 h
K∈Th K K∈Th K
(cid:90)
(cid:88)
= f(Π0 v )dx.
k−2 h
K∈Th K
6
Moreover, standard L2 orthogonality and approximation estimates yield
(cid:90)
(cid:88)
<f ,v >−(f,v ) = (Π0 f −f)(v −Π0 v )dx
h h h k−2 h k−2 h
K∈T(cid:88)h K (3.25)
≤C hk−1|f| (cid:107)v −Π0 v (cid:107) .
K k−1,K h k−2 h 0,K
K∈Th
4. Estimate in H1
We shall use duality arguments, both for deriving estimates in H1 and in
L2. In view of this, let us recall some regularity results for the problem
D∆2ψ =g in Ω, ψ =ψ =0 on ∂Ω. (4.1)
/n
Since Ω is a convex polygon, it holds (see [10])
g ∈H−1(Ω)=⇒ ψ ∈H3(Ω), (cid:107)ψ(cid:107) ≤C(cid:107)g(cid:107) , (4.2)
3 −1
and
∃p with 0<p≤1 such that
(4.3)
g ∈L2(Ω)=⇒ ψ ∈H3+p(Ω), (cid:107)ψ(cid:107) ≤C(cid:107)g(cid:107) .
3+p 0
The value of p depends on the maximum angle in Ω. Moreover, there exists a
θ <π such that, for all θ ≤θ it holds p=1, thus giving ψ ∈H4(Ω).
0 0
We shall prove the following result.
Theorem 4.1. Let w be the solution of (2.3), and let w be the solution of
h
(3.20). Then
(cid:88)
|w−w | ≤Chk(|w| +( |f|2 )1/2), (4.4)
h 1 k+1 k−1,K
K∈Th
with C a positive constant independent of h.
Proof. Let ψ ∈H2(Ω) be the solution of (4.1) with g =−∆(w−w ):
0 h
D∆2ψ =−∆(w−w ) in Ω. (4.5)
h
By (4.2) we have
(cid:107)ψ(cid:107) ≤C(cid:107)∆(w−w )(cid:107) ≤C|w−w | . (4.6)
3 h −1 h 1
Let ψ be the interpolant of ψ in W2, for which it holds (see (3.10))
I h
(cid:107)ψ−ψ (cid:107) ≤Ch3−m(cid:107)ψ(cid:107) , m=0,1,2. (4.7)
I m 3
7
Integrating by parts, using (4.5), adding and subtracting ψ , and using (2.3)
I
and (3.20) we have:
|w−w |2 =−(∆(w−w ),w−w ) =(D∆2ψ,w−w )
h 1 h h 0 h 0
=a(w−w ,ψ−ψ )+a(w−w ,ψ )
h I h I
(4.8)
=a(w−w ,ψ−ψ )+[(f,ψ )−<f ,ψ >]
h I I h I
+[a (w ,ψ )−a(w ,ψ )]=:T +T +T .
h h I h I 1 2 3
Thefirsttermiseasilyboundedthrough(2.5),andthen(3.22),(4.7),and(4.6):
T ≤Chk−1(cid:107)w(cid:107) h(cid:107)ψ(cid:107) ≤Chk(cid:107)w(cid:107) |w−w | . (4.9)
1 k+1 3 k+1 h 1
For T we use (3.25) with v =ψ . Standard interpolation estimates give
2 h I
(cid:107)ψ −Π0 ψ (cid:107) ≤(cid:107)ψ −Π0ψ (cid:107)
I k−2 I 0,K I 0 I 0,K
≤(cid:107)ψ −ψ(cid:107) +(cid:107)ψ−Π0ψ(cid:107) +(cid:107)Π0(ψ−ψ )(cid:107)
I 0,K 0 0,K 0 I 0,K
≤Ch |ψ|
K 1,K
which inserted in (3.25) gives
(cid:88)
T ≤Chk( |f|2 )1/2|w−w | . (4.10)
2 k−1,K h 1
K∈Th
ItremainstoestimateT . Addingandsubtractingw (=Π0w)andusing(3.11),
3 π k
then adding and subtracting ψ =Π0ψ and using again (3.11) we have
π 2
(cid:88)
T = (aK(w ,ψ )−aK(w ,ψ ))
3 h h I h I
K
(cid:88)
= (aK(w −w ,ψ )+aK(w −w ,ψ ))
h h π I π h I
K
(cid:88)
= (aK(w −w ,ψ −ψ )+aK(w −w ,ψ −ψ )).
h h π I π π h I π
K
From (3.15), (2.5), standard approximation estimates, (3.22) and (4.7) we de-
duce
T ≤Chk|w| |w−w | . (4.11)
3 k+1 h 1
Inserting (4.9), (4.10), and (4.11) in (4.8) we have the result (4.4)
5. Estimate in L2
We shall prove the following result.
8
Theorem 5.1. Let w be the solution of (2.3), and let w be the solution of
h
(3.20). Then
(cid:16) (cid:88) (cid:17)
h2 |w| +( |f|2 )1/2 for k =2
3 1,K
(cid:107)w−wh(cid:107)0 ≤C (cid:16) K∈Th (cid:88) (cid:17) (5.1)
hk+p |w|k+1+( |f|2k−1,K)1/2 for k ≥3,
K∈Th
with C >0 independent of h, and p the regularity index given in (4.3).
Proof. Weshalltreatonlythecasesk ≥3,thereasonbeingthat,ifβistheorder
ofconvergenceinH2,theexpectedorderinL2isgivenbymin{2β,β+2}. Hence,
for k = 2 we can expect not more than order 2, which is a direct consequence
of the H1−estimate (4.4):
(cid:16) (cid:88) (cid:17)
for k =2, (cid:107)w−w (cid:107) ≤C|w−w | ≤Ch2 |w| +( |f|2 )1/2 . (5.2)
h 0 h 1 3 1,K
K∈Th
Let then k ≥3, and let ψ ∈H2(Ω) be the solution of (4.1) with g =w−w :
0 h
D∆2ψ =w−w in Ω. (5.3)
h
By the regularity assumption (4.3) we have
(cid:107)ψ(cid:107) ≤C(cid:107)w−w (cid:107) . (5.4)
3+p h 0
Let ψ be the interpolant of ψ in W3, for which it holds
I h
(cid:107)ψ−ψ (cid:107) ≤Ch3+p−m(cid:107)ψ(cid:107) , m=0,1,2. (5.5)
I m 3+p
Then, from (5.3) and proceeding as we did in (4.8) we have
(cid:107)w−w (cid:107)2 =(D∆2ψ,w−w ) =a(ψ,w−w )
h 0 h 0 h
=a(w−w ,ψ−ψ )+[(f,ψ )−<f ,ψ >] (5.6)
h I I h I
+[a (w ,ψ )−a(w ,ψ )]=:T +T +T .
h h I h I 1 2 3
The rest of the proof follows exactly the steps used for proving Theorem 4.1.
Thus, from (3.22), (5.5) and (5.4),
T ≤Chk−1(cid:107)w(cid:107) h1+p(cid:107)ψ(cid:107) ≤Chk+p(cid:107)w(cid:107) (cid:107)w−w (cid:107) . (5.7)
1 k+1 3+p k+1 h 0
For the term T we use again (3.25) with v =ψ , which now gives
2 h I
(cid:107)ψ −Π0 ψ (cid:107) ≤(cid:107)ψ −Π0ψ (cid:107) ≤Ch2 |ψ| ,
I k−2 I 0,K I 1 I 0,K K 2,K
so that
(cid:88)
T ≤Chk+1( |f|2 )1/2(cid:107)w−w (cid:107) . (5.8)
2 k−1,K h 0
K∈Th
9
Finally, proceeding exactly as for (4.11),
T ≤Chk−1|w| h1+p(cid:107)ψ(cid:107) ≤Chk+p|w| (cid:107)w−w (cid:107) . (5.9)
3 k+1 3+p k+1 h 0
Collecting (5.7)–(5.9) in (5.6) gives
(cid:88)
(cid:107)w−w (cid:107) ≤Chk+p(|w| +( |f|2 )1/2)
h 0 k+1 k−1,K
K∈Th
and the proof is concluded.
6. Numerical results
Inordertoassessaccuracyandperformanceofvirtualelementsforplates,we
presentnumericaltestsusingthefirsttwoelementsofthefamilyheredescribed.
The corresponding polynomial degree indices, defined in (3.1), are r =3, s=1,
m = −2 and r = 3, s = 2, m = −1. Thus, the elements are named VEM31
and VEM32, respectively. The degrees of freedom, chosen according to the
definitions(3.3)-(3.7),arethevaluesofthedisplacementanditsfirstderivatives
at the vertices ((3.3) and (3.4)) for VEM31, and the same degrees of freedom
(3.3)-(3.4)plusthemomentoforderzeroofthenormalderivative(see(3.7))for
VEM32. The two elements are presented in Figure 1. They are the extensions
topolygonalelementsoftwowell-knownfiniteelementsforplates: theReduced
Hsieh-Clough-Tocher triangle (labelled CLTR), and the Hsieh-Clough-Tocher
triangle (labelled CLT) (see e.g. [8]), respectively. As a test problem we solve
(2.1)-(2.2) with Ω = unit square and f chosen to have as exact solution the
function w =x2(x−1)2y2(y−1)2. As a first test, we compare the behaviour
ex
of virtual and finite elements; for this we take a sequence of uniform meshes of
N ×N ×2 equal right triangles (N =4,8,16,32), and we plot the convergence
curves of the error in L2, H1 and H2 produced by the virtual elements VEM31
and VEM32, and the finite elements CLTR and CLT respectively. Figure 2
shows the relative errors in L2, H1 and H2 norm against the mesh diameter
(h = 0.3536, h = 0.1768, h = 0.0884, h = 0.0442). The convergence rates are
obviouslythesame,althoughinallcasesthevirtualelementsseemtoperforma
littlebetter. Next,wetestthebehaviourofthevirtualelementsonasequenceof
random Voronoi polygonal tessellations of the unit square in 25, 100, 400, 1600
polygonswithmeandiametersh=0.3071, h=0.1552, h=0.0774, h=0.0394
respectively. (Figure 3 shows the 100 and 1600-polygon meshes). In Figure 4
wecomparetheconvergencecurvesinL2, H1 andH2 normobtainedusingthe
virtual elements VEM31 and VEM32 on the Voronoi meshes and on uniform
triangular meshes.
10
