# A conforming DG method for the biharmonic equation on polytopal meshes

A conforming discontinuous Galerkin finite element method is introduced for solving the biharmonic equation. This method, by its name, uses discontinuous approximations and keeps simple formulation of the conforming finite element method at the same time. The ultra simple formulation of the method will reduce programming complexity in practice. Optimal order error estimates in a discrete H^2 norm is established for the corresponding finite element solutions. Error estimates in the L^2 norm are also derived with a sub-optimal order of convergence for the lowest order element and an optimal order of convergence for all high order of elements. Numerical results are presented to confirm the theory of convergence.

• 17 publications
• 25 publications
07/01/2019

### A conforming discontinuous Galerkin finite element method: Part II

A conforming discontinuous Galerkin (DG) finite element method has been ...
01/20/2022

### A stabilizer-free C^0 weak Galerkin method for the biharmonic equations

In this article, we present and analyze a stabilizer-free C^0 weak Galer...
12/06/2021

### A high order explicit time finite element method for the acoustic wave equation with discontinuous coefficients

In this paper, we propose a novel high order explicit time discretizatio...
11/10/2019

12/11/2020

### A mass-lumping finite element method for radially symmetric solution of a multidimensional semilinear heat equation with blow-up

This study presents a new mass-lumping finite element method for computi...
08/21/2018

### A Hybridized Discontinuous Galerkin Method for A Linear Degenerate Elliptic Equation Arising from Two-Phase Mixtures

We develop a high-order hybridized discontinuous Galerkin (HDG) method f...
12/16/2021

### A high-order residual-based viscosity finite element method for the ideal MHD equations

We present a high order, robust, and stable shock-capturing technique fo...

## 1 Introduction

We consider the biharmonic equation of the form

 (1) Δ2u = finΩ, (2) u = 0on∂Ω, (3) ∂u∂n = 0on∂Ω,

where is a bounded polytopal domain in .

The weak formulation of the boundary value problem (2) and (3) is seeking satisfying

 (4) (Δu,Δv)=(f,v)∀v∈H20(Ω).

The conforming finite element method for the problem (1)-(3) keeps the same simple form as in (4): find such that

 (5) (Δuh,Δv)=(f,v)∀v∈Vh.

However, it is known that -conforming methods require -continuous piecewise polynomials on a simplicial meshes, which imposes difficulty in practical computation. Due to the complexity in the construction of -continuous elements, -conforming finite element methods are rarely used in practice for solving the biharmonic equation.

An approach of avoiding construction of -conforming elements is to use discontinuous approximations. Due to the flexibility of discontinuous Galerkin (DG) finite element methods in element constructions and in mesh generations, many finite element methods have been developed using totally discontinuous polynomials. Here we are only interested in interior penalty discontinuous Galerkin (IPDG) methods since the proposed the method shares the same finite element spaces with IPDG method. For the biharmonic equation, interior penalty discontinuous Galerkin finite element methods have been studied in [1, 2, 3, 4, 5, 8]. One obvious disadvantage of discontinuous finite element methods is their rather complicated formulations which are often necessary to guarantee well posedness and convergence of the methods. For example, the symmetric IPDG method for the biharmonic equation with homogenous boundary conditions [1, 2] has the following formulation:

 (6) (Δuh,Δv)Th + ∑∫e({∇Δuh}⋅[v]+{∇Δv}⋅[uh])ds + ∑∫e({Δuh}⋅[∇v]+{Δv}⋅[∇uh])ds + ∑∫e(σ[uh]⋅[v]+τ[∇uh][∇v])ds=(f,v),

where and are two parameters that need to be tuned.

The purpose of this work is to introduce a conforming DG finite element method for the biharmonic equation which has the following ultra simple formulation without any stabilizing/penalty terms and other mixed terms of lower dimension integrations in (6):

 (7) (Δwuh, Δwv)=(f,v),

where is called weak Laplacian, an approximation of . The formulation (7) can be viewed as a counterpart of (5) for discontinuous approximations. The conforming DG method was first introduced in [12, 13] for second order elliptic equations, which, by name, means the method using the finite element spaces of DG methods and the simple formulations of conforming methods. This new finite element method shares the same finite element space with the IPDG methods but having much simpler formulation. This simple formulation can be obtained by defining weak Laplacian appropriately. The idea here is to raise the degree of polynomials used to compute weak Laplacian . Using higher degree polynomials in computation of weak Laplacian will not change the size, neither the global sparsity of the stiffness matrix. Optimal order error estimates in a discrete for and in norm for are established for the corresponding finite element solutions. Numerical results are provided to confirm the theories.

## 2 A Conforming DG Finite Element Method

Let be a partition of the domain consisting of polygons in two dimension or polyhedra in three dimension satisfying a set of conditions defined in [9] and additional conditions specified in [14]. Denote by the set of all edges or flat faces in , and let be the set of all interior edges or flat faces.

For simplicity, we adopt the following notations,

 (v,w)Th = ∑T∈Th(v,w)T=∑T∈Th∫Tvwdx, ⟨v,w⟩∂Th = ∑T∈Th⟨v,w⟩∂T=∑T∈Th∫∂Tvwds.

Let consist all the polynomials degree less or equal to defined on .

We define a finite element space for as follows

 (8) Vh={v∈L2(Ω): v|T∈Pk(T)T∈Th}.

Let and be two polygons/polyhedrons sharing if . Let and

be scalar and vector valued functions, the jumps

and are defined as

 (9) [v]=v|T1n1+v|T2n2,[q]=q|T1⋅n1+q|T2⋅n2,

and the averages and are defined as

 (10) {v}=12(v|T1+v|T2){q}=12(q|T1+q|T2).

If is on , then

 (11) {v}=0,{q}=0,[v]=vn,[q]=q⋅n.

The new conforming DG finite element method for the biharmonic equation (1)-(3) is defined as follows.

###### Weak Galerkin Algorithm 1

A numerical approximation for (1)-(3) can be obtained by seeking satisfying the following equation:

 (12) (Δwuh, Δwv)Th=(f,v)∀v∈Vh.

Next we will discuss how to compute the weak Laplacian and in (12). The concept of weak derivative was first introduced in [10, 9] for weak functions in weak Galerkin methods and was modified in [7, 11]. A weak Laplacian operator, denoted by , is defined as the unique polynomial for that satisfies the following equation

 (13) (Δwv, φ)T=(v, Δφ)T−⟨{v}, ∇φ⋅n⟩∂T+⟨{∇v}⋅n, φ⟩∂T,∀φ∈Pj(T).

Let , then on any ,

 (14) Δwϕ=Qh(Δϕ),

where is a locally defined projections onto on each element . It is not hard to see that for any we have

 (Δwϕ, τ)T = (ϕ, Δτ)T+⟨{∇ϕ}⋅n,τ⟩∂T−⟨{ϕ}, ∇τ⋅n⟩∂T = (ϕ,Δτ)T+⟨∇ϕ⋅n, τ⟩∂T−⟨ϕ, ∇τ⋅n⟩∂T = (Δϕ, τ)T=(QhΔϕ, τ)T,

which implies

 (15) Δwϕ=Qh(Δϕ).

We complete the proof.

## 3 Well Posedness

First we define a semi-norm as

 (16) |||v|||2=(Δwv,Δwv)Th.

Then we introduce a discrete norm as follows:

 (17) ∥v∥2,h=⎛⎝∑T∈Th(∥Δv∥2T+h−3T∥[v]∥2∂T+h−1T∥[∇v]∥2∂T)⎞⎠12.

The following lemma indicates that the two norms and are equivalent.

First we need the following trace inequality. For any function , the trace inequality holds true (see [9] for details):

 (18) ∥φ∥2e≤C(h−1T∥φ∥2T+hT∥∇φ∥2T).

There exist two positive constants and such that for any , we have

 (19) C1∥v∥2,h≤|||v|||≤C2∥v∥2,h.

For any , it follows from the definition of weak Laplacian (13) and integration by parts that

 (20) (Δwv, φ)T = (v, Δφ)T−⟨{v}, ∇φ⋅n⟩∂T+⟨{∇v}⋅n, φ⟩∂T = −(∇v, ∇φ)T+⟨v−{v}, ∇φ⋅n⟩∂T+⟨{∇v}⋅n, φ⟩∂T = (Δv, φ)T+⟨v−{v}, ∇φ⋅n⟩∂T+⟨({∇v}−∇v)⋅n, φ⟩∂T.

By letting in (20) we arrive at

 ∥Δwv∥2T = (Δv, Δwv)T+⟨v−{v}, ∇(Δwv)⋅n⟩∂T+⟨({∇v}−∇v)⋅n, Δwv⟩∂T

It is easy to see that the following equations hold true for on with ,

 (21) ∥v−{v}∥e=∥[v]∥eife⊂∂Ω,∥v−{v}∥e=12∥[v]∥eife∈E0h.

and

 (22) ∥(∇v−{∇v})⋅n∥e=∥[∇v]∥eife⊂∂Ω,|(∇v−{∇v})⋅n∥e=12∥[∇v]∥eife∈E0h.

From the trace inequality (18), (21)-(22) and the inverse inequality we have

 ∥Δwv∥2T ≤ ∥Δv∥T∥Δwv∥T+∥v−{v}∥∂T∥∇(Δwv)∥∂T+∥({∇v}−∇v)⋅n∥∂T∥Δwv∥∂T ≤ C(∥Δv∥T+h−3/2T∥[v]∥∂T+h−1/2T∥[∇v]∥∂T)∥Δwv∥T,

which implies

 ∥Δwv∥T≤C(∥Δv∥T+h−3/2T∥[v]∥∂T+h−1/2T∥[∇v]∥∂T),

and consequently

 |||v|||≤C2∥v∥2,h.

Next we will prove

 ∑h−3T∥[v]∥2∂T≤C|||v|||.

It follows from (20) that for any ,

 (23) (Δwv, φ)T = (Δv, φ)T+⟨v−{v}, ∇φ⋅n⟩∂T + ⟨({∇v}−∇v)⋅n, φ⟩∂T.

By Lemma 3.1 in [14], there exist a such that for ,

 (Δv,φ0)T=0,⟨({∇v}−∇v)⋅n, φ0⟩∂T=0, ⟨v−{v},∇φ0⋅n⟩∂T∖e=0,⟨v−{v},∇φ0⋅n⟩∂T=∥v−{v}∥2e.

and

 ∥φ0∥T≤Ch3/2T∥v−{v}∥e.

Letting in (23) yields

 ∥v−{v}∥2e = (Δwv, φ0)T≤∥Δwv∥T∥φ0∥T≤Ch3/2T∥Δwv∥T∥v−{v}∥e.

which implies

 ∥v−{v}∥e≤Ch3/2T∥Δwv∥T.

Taking the summation of the above equation over and using (21), one has

 (24)

Similarly, by Lemma 3.2 in [14], we can have

 (25) ∑T∈Thh−1T∥[∇v]∥2∂T≤C|||v|||2.

Finally, by letting in (23) we arrive at

 ∥Δv∥2T = (Δv, Δwv)T−⟨v−{v}, ∇(Δwv)⋅n⟩∂T −⟨({∇v}−∇v)⋅n, Δwv⟩∂T.

Using the trace inequality (18), the inverse inequality and (24)-(25), one has

 ∥Δv∥2T ≤ C∥Δwv∥T∥Δv∥T,

which gives

 ∑T∈Th∥Δv∥2T ≤

We complete the proof.

The finite element scheme (12) has a unique solution.

It suffices to show that the solution of (12) is trivial if . It follows that

 (Δwuh,Δwuh)Th=0.

Then the norm equivalence (19) implies , i.e.

Therefore, is a smooth harmonic function on and on . Thus we have , which completes the proof.

## 4 An Error Equation

Let . Next we derive an error equation that satisfies.

For any , we have

 (26) (Δweh,Δwv)Th=ℓ1(u,v)+ℓ2(u,v),

where

 ℓ1(u,v) = ⟨∇(QhΔu−Δu)⋅n,v−{v}⟩∂Th, ℓ2(u,v) = ⟨Δu−QhΔu,(∇v−{∇v})⋅n⟩∂T.

Testing (1) by and using the fact that and and integration by parts, we arrive at

 (f,v) = (Δ2u,v)Th = (Δu,Δv)Th−⟨Δu,∇v⋅n⟩∂Th+⟨∇(Δu)⋅n,v⟩∂Th = (Δu,Δv)Th−⟨Δu,(∇v−{∇v})⋅n⟩∂Th + ⟨∇(Δu)⋅n,v−{v}⟩∂Th.

Next we investigate the term in the above equation. Using (14), integration by parts and the definition of weak Laplacian (13), we have

 (Δu,Δv)Th=(QhΔu,Δv)Th = (v,Δ(QhΔu))T+⟨∇v⋅n, QhΔu⟩∂T−⟨v,∇(QhΔu)⋅n⟩∂T = (Δwv, QhΔu)T−⟨v−{v},∇(QhΔu)⋅n⟩∂T+⟨(∇v−{∇v})⋅n,QhΔu⟩∂T = (Δwu, Δwv)T−⟨v−{v},∇(QhΔu)⋅n⟩∂T+⟨(∇v−{∇v})⋅n,QhΔu⟩∂T.

Combining the above two equations gives

 (28) (f,v) = (Δ2u,v)Th = (Δwu, Δwv)Th−⟨v−{v},∇(QhΔu−Δu)⋅n⟩∂Th − ⟨(∇v−{∇v})⋅n,Δu−QhΔu⟩∂T,

which implies that

 (Δwu, Δwv)Th=(f,v)+ℓ1(u,v)+ℓ2(u,v).

The error equation follows from subtracting (12) from the above equation,

 (Δweh, Δwv)Th=ℓ1(u,v)+ℓ2(u,v).

We have proved the lemma.

## 5 An Error Estimate in H2

We start this section by defining some approximation operator. Let be the element-wise defined projection onto on each element .

Let and . There exists a constant such that the following estimates hold true:

 (29) ⎛⎝∑T∈ThhT∥Δw−QhΔw∥2∂T⎞⎠12≤Chk−1∥w∥k+1, (30) ⎛⎝∑T∈Thh3T∥∇(Δw−QhΔw)∥2∂T⎞⎠12≤Chk−1(∥w∥k+1+hδk,2∥w∥4).

Here is the usual Kronecker’s delta with value when and value otherwise.

The above lemma can be proved by using the trace inequality (18) and the definition of . The proof can also be found in [6].

Let , and . There exists a constant such that

 (31) |ℓ1(w,v)| ≤ Chk−1(∥w∥k+1+hδk,2∥w∥4)|||v|||. (32) |ℓ2(w,v)| ≤

Using the Cauchy-Schwartz inequality, (29), (30), (21), (22) and (19), we have

 (33) ℓ1(w,v) = ∣∣ ∣∣∑T∈Th⟨∇(Δw−QhΔw)⋅n,v−{v}⟩∂T∣∣ ∣∣ ≤ ⎛⎝∑T∈Thh3T∥∇(Δw−QhΔw)∥2∂T⎞⎠12⎛⎝∑T∈Thh−3T∥v−{v}∥2∂T⎞⎠12 ≤ ⎛⎝∑T∈Thh3T∥∇(Δw−QhΔw)∥2∂T⎞⎠12⎛⎝∑T∈Thh−3T∥[v]∥2∂T⎞⎠12 ≤ Chk−1(∥w∥k+1+hδk,2∥w∥4)|||v|||,

and

 ℓ2(w,v) = ∣∣ ∣∣∑T∈Th⟨Δw−QhΔw,(∇v−{∇v})⋅n⟩∂T∣∣ ∣∣ ≤ ⎛⎝∑T∈ThhT∥Δw−QhΔw∥2∂T⎞⎠12⎛⎝∑T∈Thh−1T∥[∇v]∥2∂T⎞⎠12 (34) ≤

We complete the proof.

Let , then

 (35) |||w−Qhw|||≤Chk−1|w|k+1.

For any , it follows from (13), integration by parts, (18) and inverse inequality that for ,

 ∥Δw(w−Qhw)∥2T = (Δw(w−Qhw),Δw(w−Qhw))T = (w−Qhw,Δ(Δw(w−Qhw)))T−⟨{w−Qhw},∇(Δw(w−Qhw))⋅n⟩∂T + ⟨{∇w−∇Qhw}⋅n,Δw(w−Qhw)⟩∂T ≤ C(h−2T∥w−Qhw∥T+h−3/2T∥w−Qhw∥∂T + h−1/2T∥∇w−∇Qhw∥∂T)∥Δw(w−Qhw)∥T ≤ Chk−1|w|k+1,T∥Δw(w−Qhw)∥T,

which implies

 ∥Δw(w−Qhw)∥T ≤ Chk−1|w|k+1,T.

Taking the summation over , we have proved the lemma.

Let be the finite element solution arising from (12). Assume that the exact solution . Then, there exists a constant such that

 (36) |||u−uh|||≤Chk−1(∥u∥k+1+hδk,2∥u∥4).

Let . Then it is straightforward to obtain

 = (Δweh,Δweh)Th = (Δweh,Δw(u−uh))Th = (Δweh,Δw(Qhu−uh))Th+(Δweh,Δw(u−Qhu))Th = (Δweh,Δwϵh)Th+(Δweh,Δw(u−Qhu))Th.

We will bound the term on right hand side of (5) first. Letting in (26) and using (31)-(32) and (35), we have

 (38) |(Δweh,Δwϵh)Th| ≤ |ℓ1(u,ϵh)|+|ℓ2(u,ϵh)| ≤ Chk−1(∥u∥k+1+hδk,2∥u∥4)|||ϵh||| ≤ Chk−1(∥u∥k+1+hδk,2∥u∥4)(|||u−Qhu|||+|||u−uh|||) ≤

To bound the second term on right hand side of (5), we have by (35),

 (39) |(Δweh,Δw(u−Qhu))Th| ≤ C|||u−Qhu||||||eh||| ≤

Combining the estimates (38) and (39) with (5), we arrive

 |||eh|||≤Chk−1(∥u∥k+1+hδk,2∥u∥4),

which completes the proof.

## 6 Error Estimates in L2 Norm

In this section, we will obtain an error bound for the finite element solution in norm.

The dual problem considered has the following form,

 (40) Δ2w = ehinΩ, (41) w = 0on∂Ω, (42) ∇w⋅n = 0on∂Ω.

Assume that the regularity holds,

 (43) ∥w∥4≤C∥eh∥.

Let be the finite element solution arising from (12). Assume that the exact solution and (43) holds true. Then, there exists a constant such that

 (44) ∥u−uh∥≤Chk+1−δk,2(∥u∥k+1+hδk,2∥u∥4).

Testing (40) by and using the fact that and and integration by parts, we arrive at

 ∥eh∥2 = (Δ2w,eh)Th = (Δw,Δeh)Th−⟨Δw,∇eh⋅n⟩∂Th+⟨∇(Δw)⋅n,eh⟩∂Th = (Δw,Δeh)Th−⟨Δw,(∇eh−{∇eh})⋅n⟩∂Th + ⟨∇(Δw)⋅n,eh−{eh}⟩∂Th. = (QhΔw,Δeh)Th+(Δw−QhΔw,Δeh)Th − ⟨Δw,(∇eh−{∇eh})⋅n⟩∂Th+⟨∇(Δw)⋅n,eh−{eh}⟩∂Th.

It follows from integration by parts, the definition of weak Laplacian (13) and (14),

 (QhΔw,Δeh)Th = (eh,Δ(QhΔw))T+⟨∇eh⋅n, QhΔw⟩∂T−⟨eh,∇(QhΔw)⋅n⟩∂T = (Δweh, QhΔw)T−⟨eh−{eh},∇(QhΔw)⋅n⟩∂T + ⟨(∇eh−{∇eh})⋅n,QhΔw⟩∂T = (Δww, Δweh)T−⟨eh−{eh},∇(QhΔw)⋅n⟩∂T + ⟨(∇eh−{∇eh