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

A new stabilizer free weak Galerkin (WG) method is introduced and analyzed for the biharmonic equation. Stabilizing/penalty terms are often necessary in the finite element formulations with discontinuous approximations to ensure the stability of the methods. Removal of stabilizers will simplify finite element formulations and reduce programming complexity. This stabilizer free WG method has an ultra simple formulation and can work on general partitions with polygons/polyhedra. Optimal order error estimates in a discrete H^2 for k> 2 and in L^2 norm for k>2 are established for the corresponding weak Galerkin finite element solutions. Numerical results are provided to confirm the theories.

## Authors

• 17 publications
• 25 publications
06/16/2019

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

A stabilizing/penalty term is often used in finite element methods with ...
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...
06/15/2022

### Neural Control of Discrete Weak Formulations: Galerkin, Least-Squares and Minimal-Residual Methods with Quasi-Optimal Weights

There is tremendous potential in using neural networks to optimize numer...
12/13/2021

### Stabilized finite element methods for a fully-implicit logarithmic reformulation of the Oldroyd-B constitutive law

Logarithmic conformation reformulations for viscoelastic constitutive la...
11/24/2020

### The Weak Galerkin Finite Element Method for the Transport-Reaction Equation

We present and analyze a weak Galerkin finite element method for solving...
11/20/2019

### A Locking-Free P_0 Finite Element Method for Linear Elasticity Equations on Polytopal Partitions

This article presents a P_0 finite element method for boundary value pro...
09/15/2019

### Weak discrete maximum principle of finite element methods in convex polyhedra

We prove that the Galerkin finite element solution u_h of the Laplace eq...
##### 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 the biharmonic equation of the form

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

where is a bounded polytopal domain in .

For the biharmonic problem (1) with Dirichlet and Neumann boundary conditions (2) and (3), the corresponding weak form is given by seeking satisfying and such that

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

where is the subspace of consisting of functions with vanishing value and normal derivative on .

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.

As an alternative approach,  nonconforming  and discontinuous finite element methods have been developed for solving the biharmonic equation over the last several decades. The Morley element [2] is a well-known example of nonconforming element for the biharmonic equation by using piecewise quadratic polynomials. The weak Galerkin finite element methods use discontinuous approximations on general polytopal meshes introduced first in [9]. Many WG finite element methods have been developed for forth order problems [3, 4, 5, 6, 7, 10, 12]. These weak Galerkin finite element methods for (1)-(3) have the following symmetric, positive definite and parameter independent formulation:

 (5) (Δwuh,Δwv)+s(uh,v)=(f,v).

The stabilizer in (5) is necessary to guarantee the well posedness and the convergence of the methods.

The purpose of the work is to further simplify the WG formulation (5) by removing the stabilizer to obtain an ultra simple formulation for the biharmonic equation:

 (6) (Δwuh,Δwv)=(f,v).

We can obtain a stabilizer free WG method (6) by appropriately designing the weak Laplacian . The idea 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.

This new stabilizer free WG method for the forth order problem (2)-(3) has an ultra simple symmetric positive definite formulation (6) and can work on general polytopal meshes. For second order elliptic problems, stabilizer free WG methods have been studied in [1, 9, 11]. However for forth order problems, to the best of our knowledge, this is the first finite element method without any stabilizers for totally discontinuous approximations. Optimal order error estimates in a discrete norm is established for the corresponding WG finite element solutions. Error estimates in the 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.

## 2 Weak Galerkin Finite Element Methods

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 [8] and additional conditions specified in Lemma 3 and Lemma 3. 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 .

First we introduce a set of normal directions on as follows

 (7) Dh={ne: ne is unit and normal to e, e∈Eh}.

Then, we can define a weak Galerkin finite element space for as follows

 (8) Vh={v={v0,vb,vnne}: v0∈Pk(T),vb∈Pk(e),vn∈Pk−1(e),e⊂∂T},

where can be viewed as an approximation of .

Denote by a subspace of with vanishing traces,

 V0h={v={v0,vb,vnne}∈Vh,vb|e=0, vnne⋅n|e=0, e⊂∂T∩∂Ω}.

A weak Laplacian operator, denoted by , is defined as the unique polynomial for that satisfies the following equation

 (9) (Δwv, φ)T=(v0, Δφ)T−⟨vb, ∇φ⋅n⟩∂T+⟨vnne⋅n, φ⟩∂T,∀φ∈Pj(T).

Let , and be the locally defined projections onto , and accordingly on each element and . For the true solution of (1)-(3), we define as

 Qhu={Q0u,Qbu,Qn(∇u⋅ne)ne}∈Vh.
###### Weak Galerkin Algorithm 1

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

 (10) (Δwuh, Δwv)Th=(f,v0)∀ v={v0,vb, vnne}∈V0h.

Let , then on any ,

 (11) Δ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+⟨(∇ϕ⋅ne)ne⋅n,τ⟩∂T−⟨ϕ, ∇τ⋅n⟩∂T = (ϕ,Δτ)T+⟨∇ϕ⋅n, τ⟩∂T−⟨ϕ, ∇τ⋅n⟩∂T = (Δϕ, τ)T=(QhΔϕ, τ)T,

which implies

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

It completes the proof.

## 3 Well Posedness

For any , let

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

We introduce a discrete norm as follows:

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

For any function , the following trace inequality holds true [8],

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

The main goal of this section is to obtain the equivalence of the two norms and . To do so, we need the two following lemmas.

Let be a convex polygon/polyhedron of size with edges/faces , , …. Let , and . Let , , and where is the barycenter of . For any , there is a unique polynomial for some such that

 (16) (q,p)T =0∀p∈Pk−1(T), (17) ⟨∇q⋅n−f,p⟩e1 =0∀p∈Pk(e1), (18) ∥q∥T ≤Ch3/2T∥f∥e1,

where depends on the minimum angle and the smallest ratio , and is defined in (26) below.

We prove is uniquely defined by (16)–(17). Let in (17). As is convex, in the interior of for all . Because of the positive weight, the vanishing weighted inner-product forces on :

 (19) ⟨∇q⋅n,p⟩e1=−1hT⟨λ22⋯λ2nqk,p⟩e1=0∀p∈Pk(e1).

Thus the vanishing weighted inner-product forces on :

 (20) (q,p)T=(λ21⋯λ2nqk−1,p)T=0∀p∈Pk−1(T),

where .

We find some upper bounds and lower bounds of these weight functions .

Let , ( in 2D), be a neighboring edge/face of . Using the distance as its variable, we have

 λi|e1=2he1x,

where is the distance, along , of the point to the intersection of and . Here in 3D, we assume the size of is roughly twice the distance from the barycenter to the intersection edge . To avoid too many constants, we simply assume reasonably . We compute the maximum as

 (21) maxTλi=h⊥ei(T)(he1/2)sinαi≤2hThe1sinαi≤8sinαi≤8sinα0,

where (, and ) is the angle between and , and is the maximal distance of a point on to . For a lower bound, we have

 (22) λi|T0 ≥⎧⎨⎩1516if αi≤π/2,1−√d16sinαi≥12if αi>π/2,

where is a square/cube at middle of with size , cf. Figure 1. We note that other than triangles, for most other polygons. Here in (22), we assumed , where is the space dimension, 2 or 3.

For non-neighboring edges , we have

 λj|e1=⎧⎨⎩1if ej∥e1,2(x+xj)he1+xjotherwise,

where is the arc-length parametrization on toward the extended intersection of and , is the distance on from the an boundary point of to the intersection. Supposing is the only edge/polygonal between and , . Because , it follows that

 (23) maxTλj=h⊥ej(T)(he1/2)sinαi≤2hT(he1+xj)sin(αi+αj)≤8sinα0.

For a lower bound, because and is an edge/polygon in between, we have

 (24) λj|T0 ≥λi|T0≥12.

Together, in (19) and (20), we have, noting ,

 (25) λ22⋯λ2n|T0≥122n−2, and  λ21λ22⋯λ2n|T≤82n−2sin2n−2α0.

Let be the unique solution in (17), i.e., . Letting in (17), cf. (19), we get, by (25),

 1162k122n−2∥~qk∥2e1 ≤122n−2∥~qk∥2e0≤∥λ2⋯λn~qk∥2e0 ≤∥λ2⋯λn~qk∥2e1=⟨λ22⋯λ2n~qk,~qk⟩e1 =−1hT⟨λ22⋯λ2nqk,p⟩e1=⟨f,p⟩e1=⟨f,−hT~qk⟩e1 ≤hT∥f∥e1∥~qk∥e1,

where in the first step we use the fact is a degree polynomial. For the unique solution , we view it as a polynomial on the whole line or whole plane containing . We also extend it to by letting it be constant in the direction orthogonal to . Let be a square/cube of size containing , with one side which contains . It follows that, by (25),

 ∥~qk∥2T ≤∥~qk∥2ST=hT∥~qk∥2Se1≤hT(hThe1)2k∥~qk∥2e1 ≤hT42k∥~qk∥2e1≤hT42k(28k+2n−2hT∥f∥e1)2 =220k+4n−4h3T∥f∥2e1.

We rewrite in terms of this extended polynomial,

 q=λ1λ22⋯λ2nqk=λ1λ22⋯λ2n(λ1qk−1+~qk)for some  qk−1∈Pk−1(T).

Letting in (16), it follows that, by (25),

 ∥qk−1∥2T ≤(hT/he0)2k−2∥qk−1∥2T0≤642k−222n−2(λ22⋯λ2nqk−1,qk−1)T0 ≤22n+12k−14(2hT/he0)2(λ21⋯λ2nqk−1,qk−1)T0,0 ≤22n+12k(2hT/he0)2n+2k−2(λ21⋯λ2nqk−1,qk−1)T =216n+26k−14(λ1λ22⋯λ2n~qk,−qk−1)T ≤216n+26k−1482n−2sin2n−2α0∥~qk∥T∥qk−1∥T,

where is the top half of , and we used the fact and the fact that the integrant on is a degree polynomial. We estimate

 ∥q∥2T =(λ21λ42⋯λ4n(λ1qk−1+~qk),λ1qk−1+~qk)T ≤84n−4sin4n−4α0(λ1qk−1+~qk,λ1qk−1+~qk)T ≤84n−4sin4n−4α02(∥λ1qk−1∥2T+∥~qk∥2T) ≤212n−11sin4n−4α0(∥qk−1∥2T+∥~qk∥2T).

Combining above three bounds, we get

 (26) ∥q∥T≤(212n−11sin4n−4α0)12((222n+26k−20sin2n−2α0)2+1)12∥~qk∥T≤(216n+20k−15sin4n−4α0)12((222n+26k−20sin2n−2α0)2+1)12h3/2T∥f∥e1=:Ch3/2T∥f∥e1.

The proof is completed.

Let the following notations be defined in Lemma 3. For any , there is a unique polynomial for some such that

 (27) (q,p)T =0∀p∈Pk−1(T), (28) ⟨∇q⋅n,p⟩e1 =0∀p∈Pk(e1), (29) ⟨q−g,p⟩e1 =0∀p∈Pk+1(e1), (30) ∥q∥T ≤Ch1/2T∥g∥e1,

where depends on the minimum angle and the smallest ratio , and is defined in (31) below.

For unisovence, letting in (29), we get for some , because the weights . By (28), and thus for some . By (27), and thus .

The upper and lower bounds for are same as that in Lemma 3.

Let be the unique solution in (29), i.e., . Letting in (29), we get, by (25),

 1162k+2122n−2∥~qk+1∥2e1 ≤122n−2∥~qk+1∥2e0≤∥λ2⋯λn~qk+1∥2e0 ≤∥λ2⋯λn~qk+1∥2e1=⟨λ22⋯λ2n~qk+1,~qk+1⟩e1 =⟨g,~qk+1⟩e1≤∥g∥e1∥~qk+1∥e1,

where in the first step we use the fact is a degree polynomial. For the unique solution , we view it as a polynomial on the whole line or whole plane containing . We also extend it to by letting it be constant in the direction orthogonal to . Let be a square/cube of size containing , with one side which contains . It follows that, by (25),

 ∥~qk+1∥2T ≤∥~qk+1∥2ST=hT∥~qk+1∥2Se1≤hT(hThe1)2k+2∥~qk+1∥2e1 ≤hT42k+2∥~qk+1∥2e1≤hT42k+2(28k+2n+6∥g∥e1)2 =220k+4n+16hT∥g∥2e1.

We rewrite in terms of this extended polynomial,

 q=λ22⋯λ2nqk+1=λ22⋯λ2n(λ1qk+~qk+1)for some  qk∈Pk(T).

By (28), we have further, because ,

 q=λ22⋯λ2n(λ21qk−1+~qk+1)for some  qk−1∈Pk−1(T).

Letting in (27), it follows that, by (25),

 ∥qk−1∥2T ≤(hT/he0)2k−2∥qk−1∥2T0≤642k−222n−2(λ22⋯λ2nqk−1,qk−1)T0 ≤22n+12k−14(2hT/he0)2(λ21⋯λ2nqk−1,qk−1)T0,0 ≤22n+12k(2hT/he0)2n+2k−2(λ21⋯λ2nqk−1,qk−1)T =216n+26k−14(λ22⋯λ2n~qk+1,−qk−1)T ≤216n+26k−1482n−2sin2n−2α0∥~qk+1∥T∥qk−1∥T,

where is the top half of , and we used the fact and the fact that the integrant on is a degree polynomial. We estimate

 ∥q∥2T =(λ42⋯λ4n(λ21qk−1+~qk+1),λ21qk−1+~qk+1)T ≤84n−4sin4n−4α0(λ21qk−1+~qk+1,λ21qk−1+~qk+1)T ≤84n−4sin4n−4α02(∥λ21qk−1∥2T+∥~qk+1∥2T) ≤212n−11sin4n−4α0(∥qk−1∥2T+∥~qk+1∥2T).

Combining above three bounds, we get

 (31) ∥q∥T≤(212n−11sin4n−4α0)12((216n+26k−14sin2n−2α0)2+1)12∥~qk∥T≤(216n+20k+5sin4n−4α0)12((216n+26k−14sin2n−2α0)2+1)12hT∥g∥e1=:ChT∥g∥e1.

The proof is completed.

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

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

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

 (33) (Δwv, φ)T = (v0, Δφ)T−⟨vb, ∇φ⋅n⟩∂T+⟨vnne⋅n, φ⟩