# On stabilizer-free weak Galerkin finite element methods on polytopal meshes

A stabilizing/penalty term is often used in finite element methods with discontinuous approximations to enforce connection of discontinuous functions across element boundaries. Removing stabilizers from discontinuous finite element methods will simplify the formulations and reduce programming complexity significantly. The goal of this paper is to introduce a stabilizer free weak Galerkin finite element method for second order elliptic equations on polytopal meshes in 2D or 3D. This new WG method keeps a simple symmetric positive definite form and can work on polygonal/polyheral meshes. Optimal order error estimates are established for the corresponding WG approximations in both a discrete H^1 norm and the L^2 norm. Numerical results are presented verifying the theorem.

## Authors

• 17 publications
• 25 publications
07/19/2019

### A stabilizer free weak Galerkin method for the Biharmonic Equation on Polytopal Meshes

A new stabilizer free weak Galerkin (WG) method is introduced and analyz...
07/01/2019

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

A conforming discontinuous Galerkin (DG) finite element method has been ...
08/28/2020

### A stabilizer free weak Galerkin finite element method on polytopal mesh: Part II

A stabilizer free weak Galerkin (WG) finite element method on polytopal ...
10/22/2019

### Finite Element Methods for Maxwell's Equations

We survey finite element methods for approximating the time harmonic Max...
03/30/2021

### Hybrid high-order and weak Galerkin methods for the biharmonic problem

We devise and analyze two hybrid high-order (HHO) methods for the numeri...
01/08/2021

### Projection in negative norms and the regularization of rough linear functionals

In order to construct regularizations of continuous linear functionals a...
06/06/2021

### Leveraging spectral analysis to elucidate membrane locking and unlocking in isogeometric finite element formulations of the curved Euler-Bernoulli beam

In this paper, we take a fresh look at using spectral analysis for asses...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

We consider Poisson equation with a homogeneous Dirichlet boundary condition in

dimension as our model problem for the sake of clear presentation. This stabilizer free weak Galerkin method can be extended to solve other partial differential equations. The Poisson problem seeks an unknown function

satisfying

 (1) −Δu = f,inΩ, (2) u = 0,on∂Ω,

where is a polytopal domain in .

The weak form of the problem (1)-(2) is to find such that

 (3) (∇u,∇v)=(f,v),∀v∈H10(Ω).

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

 (4) (∇uh,∇v)=(f,v),∀v∈Vh,

where is a finite dimensional subspace of . The functions in are required to be continuous that makes the classic conforming finite element formulation (4) less flexible in construction of high order elements and in mesh generation. In contrast, finite element methods using discontinuous approximations have two significant advantages: 1. Easy construction of high order elements and avoiding constructing some special elements such as conforming elements; 2. Easy working on general meshes. Therefore, discontinuous finite element methods are the most active research area in the context of finite element methods for the past two decades. Discontinuous approximation was first used in finite element procedure as early as in 1970s [2, 7, 14, 20]. Local discontinuous Galerkin methods were introduced in [6]. Then a paper [1] in 2002 provides a unified analysis of discontinuous Galerkin (DG) finite element methods for Poisson equation. More discontinuous finite element methods have been developed such as hybridizable discontinuous Galerkin (HDG) method [5], mimetic finite differences method [11], hybrid high-order (HHO) method [13], weak Galerkin (WG) method [17] and references therein.

One obvious disadvantage of discontinuous finite element methods is their rather complex formulations which are often necessary to ensure continuity of discontinuous solutions across element boundaries. Most of discontinuous finite element methods have one or more stabilizing terms to guarantee stability and convergence of the methods. Existing of a stabilizing terms further complicates formulations. Complexity of discontinuous finite element methods makes them difficult to be implemented and to be analyzed. The purpose of this paper is to obtain a finite element formulation close to its original PDE weak form (3) for discontinuous polynomials . We believe that finite element formulations for discontinuous approximations can be as simple as follows:

 (5) (∇wuh,∇wv)=(f,v),

if , an approximation of gradient, is appropriately defined. The formulation (5) can be viewed as a revised version of (3) for discontinuous approximations. In fact such an ultra simple formulation (5) has been achieved for one kind of WG method in [17], and for a conforming DG method in [21], on triangular/tetrahedral meshes. The lowest order WG method developed in [17] has been improved in [12] for convex polygonal meshes, in which non-polynomial functions are used for computing weak gradient.

In this paper, we develop a discontinuous finite element method that has an ultra simple weak formulation (5) on polytopal meshes for any polynomial degree . It is achieved by starting with the WG method in [9, 16]. It is interesting to see that stabilizer can be removed by raising the degree of polynomials in the definition of the weak gradient. The degree of polynomial of the weak gradient depends on the space dimension and the shape of element. An higher polynomial degree for the weak gradient will not change the size, neither the sparsity of the stiffness matrix. On the other side, without a stabilizer, this new WG method reduces programming complexity by half in two dimensional space and a lot more than half in three dimensional space by avoiding lower-dimensional integration completely. Optimal order error estimates are established for the corresponding WG approximations in both a discrete norm and the norm. Numerical results are presented to verify the theorem.

## 2 Weak Galerkin Finite Element Schemes

Let be a partition of the domain consisting of polygons in two dimension or polyhedra in three dimension satisfying a set of conditions specified in [18]. Denote by the set of all edges or flat faces in , and let be the set of all interior edges or flat faces. For every element , we denote by its diameter and mesh size for .

We start by introducing weak function on element such that

 v={v0,inT,vb,on∂T.

If is continuous on , then .

For a given integer , let be the weak Galerkin finite element space associated with defined as follows

 (6) Vh={v={v0,vb}:v0|T∈Pk(T), vb|e∈Pk(e), e⊂∂T,T∈Th}

and its subspace is defined as

 (7) V0h={v: v∈Vh, vb=0 on ∂Ω}.

We would like to emphasize that any function has a single value on each edge .

For given and , a weak gradient () is defined as the unique polynomial satisfying

 (8) (∇wv,q)T=−(v0,∇⋅q)T+⟨vb,q⋅n⟩∂T,∀q∈[Pj(T)]d.

Let and be the element-wise defined projections onto and with on respectively. Define . Let be the element-wise defined projection onto on each element .

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.
###### Weak Galerkin Algorithm 1

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

 (9) (∇wuh,∇wv)Th=(f,v0),∀v={v0,vb}∈V0h.

For , then on any ,

 (10) ∇wϕ=Qh∇ϕ.

Using (8) and integration by parts, we have that for any

 (∇wϕ,q)T = −(ϕ,∇⋅q)T+⟨ϕ,q⋅n⟩∂T = (∇ϕ,q)T=(Qh∇ϕ,q)T,

which implies the desired identity (10).

## 3 Well Posedness

For any , let

 (11) |||v|||2=(∇wv,∇wv)Th.

We introduce a discrete semi-norm as follows:

 (12) ∥v∥1,h=⎛⎝∑T∈Th(∥∇v0∥2T+h−1T∥v0−vb∥2∂T)⎞⎠12.

The following lemma indicates that is equivalent to the norm in (11).

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

 (13) C1∥v∥1,h≤|||v|||≤C2∥v∥1,h.

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

 (14) (∇wv,q)T=(∇v0,q)T+⟨vb−v0,q⋅n⟩∂T,∀q∈[Pj(T)]d.

By letting in (14) we arrive at

 (∇wv,∇wv)T=(∇v0,∇wv)T+⟨vb−v0,∇wv⋅n⟩∂T.

From the trace inequality (33) and the inverse inequality we have

 ∥∇wv∥2T ≤ ∥∇v0∥T∥∇wv∥T+∥v0−vb∥∂T∥∇wv∥∂T ≤ ∥∇v0∥T∥∇wv∥T+Ch−1/2T∥v0−vb∥∂T∥∇wv∥T,

which implies

 ∥∇wv∥T≤C(∥∇v0∥T+h−1/2T∥v0−vb∥∂T),

and consequently

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

Next we will prove . For and , by (8) and integration by parts, we have

 (15) (∇wv,q)T=(∇v0,q)T+⟨vb−v0,q⋅n⟩∂T.

We like to find a such that,

 (16) (∇v0,q)T=0,⟨vb−v0,q⋅n⟩∂T=⟨vb−v0,q⋅n⟩e.

One can construct such satisfying (16) by choosing a large enough polynomial degree . Let us give a detailed construction. For a 2D -side polygon , we list its boundary edges as , , …, . Let

. Using finite element’s moments dof (degrees of freedom), we define a

by

 (17) (q,p)T =0 ∀p∈[Pk−1(T)]2, (18) ⟨q⋅n,p⟩e =⟨vb−v0,p⟩e ∀p∈Pn+k−1(e), (19) ⟨q⋅t,p⟩e =0 ∀p∈Pn+k−1(e), (20) ⟨q,p⟩ei =0 ∀p∈[Pn+k−1−i(ei)]2,i=1,...,n−1,

where

is a tangent vector to edge

. We note that the above linear system of equations is square. We show its uniqueness of solution when the right hand side vector vanishes. By (18)–(19), , where and such that . Sequentially, by (20), , where and such that . (17) implies .

For a -side polyhedron in 3D, we list its face polygons as , , …, . Let and be two orthogonal, tangent vectors of . Let . is uniquely defined by

 (21) (q,p)T =0 ∀p∈[Pk−1(T)]3, (22) ⟨q⋅n,p⟩e =⟨vb−v0,p⟩e ∀p∈Pn+k−1(e), (23) ⟨q⋅tj,p⟩e =0 ∀p∈Pn+k−1(e), j=1,2, (24) ⟨q,p⟩ei =0 ∀p∈[Pn+k−1−i(ei)]3,i=1,...,n−1.

The proof of uniqueness is identical to the proof above for the 2D case.

Substituting the above solution into (15), we get

 (25) (∇wv,q)T=⟨vb−v0,q⋅n⟩e.

It follows from Cauchy-Schwarz inequality that

 |⟨vb−v0,q⋅n⟩e|≤C∥∇wv∥T∥q∥T.

On the reference element , as all norms are equivalent on a finite dimensional vector space

 V={^q∣the unique solution in (???)--% (???) for a given ^q⋅^n|^e∈Pj(^e) },

we get

 ∥^q∥^T≤C∥^q⋅^n∥^e∀^q∈V,

with an absolute constant independent of but depending on the worst shape of , cf. Stenberg [15] and Wang-Ye [18]. After an affine mapping, it follows

 ∥q∥T≤Ch1/2T∥q⋅n∥e.

Then

 (26) |⟨vb−v0,q⋅n⟩e|≤Ch1/2T∥∇wv∥T∥q⋅n∥e.

We would get

 ∥q⋅n∥2e =∥ΠPk(e)(q⋅n)∥2e+∥ΠPk(e)⊥(q⋅n)∥2e =∥vb−v0∥2e,

where is the orthogonal projection to on the edge or polygon . Thus, we rewrite (26) as

 (27) ∥vb−v0∥2e≤Ch1/2T∥∇wv∥T∥q⋅n∥e≤Ch1/2T∥∇wv∥T∥vb−v0∥e.

Repeating (27) on all boundary edges/polygons, we have

 (28) h−1/2T∥v0−vb∥∂T≤C∥∇wv∥T.

It follows from the trace inequality, the inverse inequality and (28),

 ∥∇v0∥2T≤∥∇wv∥T∥∇v0∥T+Ch−1/2T∥v0−vb∥∂T∥∇v0∥T≤C∥∇wv∥T∥∇v0∥T.

Combining the above estimate and (28), by the definition (12), we prove the lower bound of (13) and complete the proof of the lemma.

The weak Galerkin finite element scheme (9) has a unique solution.

If and are two solutions of (9), then would satisfy the following equation

 (∇wεh,∇wv)=0,∀v∈V0h.

Note that . Then by letting in the above equation we arrive at

It follows from (13) that implies , or equivalently, and . Thus , i.e., . This completes the proof of the lemma.

## 4 Error Estimates in Energy Norm

Let and . Next we derive an error equation that satisfies.

For any , then,

 (29) (∇weh,∇wv)Th=ℓ(u,v),

where

 ℓ(u,v) = ⟨(∇u−Qh∇u)⋅n,v0−vb⟩∂Th.

For , testing (1) by and using the fact that , we arrive at

 (30) (∇u,∇v0)Th−⟨∇u⋅n,v0−vb⟩∂Th=(f,v0).

It follows from integration by parts, (8) and (10) that

 (31) (∇u,∇v0)Th = (Qh∇u,∇v0)Th = = (Qh∇u,∇wv)Th+⟨v0−vb,Qh∇u⋅n⟩∂Th = (∇wu,∇wv)Th+⟨v0−vb,Qh∇u⋅n⟩∂Th.

Combining (30) and (31) gives

 (32) (∇wu,∇wv)Th = (f,v0)+ℓ(u,v).

The error equation follows from subtracting (9) from (32),

 (∇weh,∇wv)Th=ℓ(u,v),∀v∈V0h.

This completes the proof of the lemma.

For any function , the following trace inequality holds true (see [18] for details):

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

Assume that is shape regular. Then for any and , we have

 (34) |ℓ(w,v)| ≤

Using the Cauchy-Schwarz inequality, the trace inequality (33) and (13), we have

 |ℓ(w,v)| = ∣∣ ∣∣∑T∈Th⟨(∇w−Qh∇w)⋅n,v0−vb⟩∂T∣∣ ∣∣ ≤ C∑T∈Th∥(∇w−Qh∇w)∥∂T∥v0−vb∥∂T ≤ C⎛⎝∑T∈ThhT∥(∇w−Qh∇w)∥2∂T⎞⎠12⎛⎝∑T∈Thh−1T∥v0−vb∥2∂T⎞⎠12 ≤

which proves the lemma.

Let be the weak Galerkin finite element solution of (9). Assume the exact solution . Then, there exists a constant such that

 (35) |||u−uh|||≤Chk|u|k+1.

It is straightforward to obtain

 = (∇weh,∇weh)Th = (∇wu−∇wuh,∇weh)Th = (∇wQhu−∇wuh,∇weh)Th+(∇wu−∇wQhu,∇weh)Th = (∇weh,∇wϵh)Th+(∇wu−∇wQhu,∇weh)Th.

We will bound each terms in (4). It follows from (8), integration by parts and the trace inequality,

 (∇w(u−Qhu),q)T = −(u−Q0u,∇⋅q)T+⟨u−Qbu,q⋅n⟩∂T = (∇(u−Q0u),q)T+⟨Q0u−Qbu,q⋅n⟩∂T ≤ ∥∇(u−Q0u)∥T∥q∥T+Ch−1/2∥u−Q0u∥∂T∥q∥T ≤ Chk|u|k+1,T∥q∥T.

Letting in the above equation and taking summation over , we have

 (37) |||u−Qhu|||≤Chk|u|k+1.

Letting in (29) and using (34) and (37), we have

 (38) |(∇weh,∇wϵh)Th| = |ℓ(u,ϵh)| ≤ Chk|u|k+1|||ϵh||| ≤ Chk|u|k+1|||Qhu−uh||| ≤ Chk|u|k+1(|||Qhu−u|||+|||u−uh|||) ≤

The estimate (37) implies

 (39) |(∇wu−∇wQhu,∇weh)Th| ≤ C|||u−Qhu||||||eh||| ≤

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

 |||eh|||≤Chk|u|k+1,

which completes the proof.

## 5 Error Estimates in L2 Norm

The standard duality argument is used to obtain error estimate. Recall . The considered dual problem seeks satisfying

 (40) −ΔΦ = e0,inΩ.

Assume that the following -regularity holds

 (41) ∥Φ∥2≤C∥e0∥.

Let be the weak Galerkin finite element solution of (9). Assume that the exact solution and (41) holds true. Then, there exists a constant such that

 (42) ∥u−u0∥≤Chk+1|u|k+1.

Testing (40) by and using the fact that give

 (43) ∥e0∥2 = −(ΔΦ,e0) = (∇Φ, ∇e0)Th−⟨∇Φ⋅n, e0−eb⟩∂Th.

Setting and in (31) yields

 (44) (∇Φ,∇e0)Th=(∇wΦ,∇weh)Th+⟨Qh∇Φ⋅n, e0−eb⟩∂Th.

Substituting (44) into (43) gives

 (45) ∥e0∥2 = (∇weh, ∇wΦ)Th−⟨(∇Φ−Qh∇Φ)⋅n, e0−eb⟩∂Th = (∇weh, ∇wQhΦ)Th+(∇weh, ∇w(Φ−QhΦ))Th+ℓ(Φ,eh) = ℓ(u,QhΦ)+(∇weh, ∇w(Φ−QhΦ))Th+ℓ(Φ,eh).

Next we will estimate all the terms on the right hand side of (45). Using the Cauchy-Schwarz inequality, the trace inequality (33) and the definitions of and we obtain

 |ℓ(u,QhΦ)| ≤ ∣∣⟨(∇u−Qh∇u)⋅n,Q0Φ−QbΦ⟩∂Th∣∣ ≤ ⎛⎝∑T∈Th∥(∇u−Qh∇u)∥2∂T⎞⎠1/2⎛⎝∑T∈Th∥Q0Φ−QbΦ∥2∂T⎞⎠1/2 ≤ C⎛⎝∑T∈Thh∥(∇u−Qh∇u)∥2∂T⎞⎠1/2⎛⎝∑T∈Thh−1∥Q0Φ−Φ∥2∂T⎞⎠1/2 ≤ Chk+1|u|k+1|Φ|2.

It follows from (35) and (37) that

 |(∇weh, ∇w(Φ−QhΦ))Th| ≤ ≤ Chk+1|u|k+1|Φ|2.

Similarly, it follows from (13) and (35) that

 |ℓ(Φ,eh)|=∣∣ ∣∣∑T∈Th⟨(Qh∇Φ−∇Φ)⋅n, e0−eb⟩∂T∣∣ ∣∣ ≤ Chk+1|u|k+1∥Φ∥2.

Combining all the estimates above with (45) yields

 ∥e0∥2≤Chk+1|u|k+1∥Φ∥2.

The estimate (42) follows from the above inequality and the regularity assumption (41). We have completed the proof.

## 6 Numerical Experiments

We solve the following Poisson equation on the unit square:

 (46) −Δu=2π2sinπxsinπy,(x,y)∈Ω=(0,1)2,

with the boundary condition on .

In the first computation, the level one grid consists of two unit right triangles cutting from the unit square by a forward slash. The high level grids are the half-size refinements of the previous grid. The first levels of grids are plotted in Figure 1. We apply finite element methods and list the error and the order of convergence in the following table, Table 1. The numerical results confirm the convergence theory.

In the next computation, we use a family of polygonal grids (with 12-side polygons) shown in Figure 2. We can see, from Table 2 that the polynomial degree for the weak gradient must be a little bigger, as indicated in the proof of (16). The convergence history confirms the theory.