 # Convergence and optimality of an adaptive modified weak Galerkin finite element method

An adaptive modified weak Galerkin method (AmWG) for an elliptic problem is studied in this paper, in addition to its convergence and optimality. The weak Galerkin bilinear form is simplified without the need of the skeletal variable, and the approximation space is chosen as the discontinuous polynomial space as in the discontinuous Galerkin method. Upon a reliable residual-based a posteriori error estimator, an adaptive algorithm is proposed together with its convergence and quasi-optimality proved for the lowest order case. The major tool is to bridge the connection between weak Galerkin method and the Crouzeix-Raviart nonconforming finite element. Unlike the traditional convergence analysis for methods with a discontinuous polynomial approximation space, the convergence of AmWG is penalty parameter free.

## Authors

##### 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

Consider the following model second order elliptic problem

 (1.1) −∇⋅(A∇u) =fin Ω, u =0on ∂Ω,

where is a bounded polygonal or polyhedral domain in . Assume that an initial conforming partition of exists and for all , the coefficient is assumed to be a piecewise constant with respect to this partition.

Weak Galerkin (WG) is a novel numerical method for solving partial differential equations in which classical differential operators (such as gradient, divergence, curl) are approximated in a weak sense. WG method was initially introduced in

[38, 39, 29] for the second order elliptic problem. Since then, the WG method has successfully found its way to many applications, for example, elliptic interface problems , Helmholtz equations [31, 28, 16], biharmonic equations [26, 30], Navier-Stokes equations [23, 20]. In particular, Wang et al.  introduced a modified weak Galerkin (mWG) method for the Poisson equation. The mWG method has been successfully applied to, such as signorini and obstacle problems , Stokes equations . We shall study an mWG in this paper.

The solution to (1.1) may contain singularity. To approximate problem (1.1) efficiently, the general practice is to adopt adaptivity by designing an adaptive finite element cycle through the help of the a posteriori error estimators, a bulk marking strategy, and certain local refinement techniques. Knowing some counterexample of un-convergent algorithms , the convergence analysis of an adaptive algorithm is of fundamental importance as it theoretically guarantees that the correct approximation will be obtained, especially we want to avoid the situation when more computational resources may go wasted after iterations of refinements.

Convergence theory of adaptive finite element method is relatively mature, see  and the references therein. Nevertheless, there are few research results for the a posterior error estimates for WG methods. Chen et al.  presented the a posteriori error estimates for second order elliptic problems; Zhang and Chen  proposed a residual-type error estimator and proved global upper and lower bounds of WG method for second order elliptic problems in a discrete -norm; Li et al.  introduced a simple posteriori error estimator which can be applied to general meshes such as hybrid, polytopal and those with hanging nodes for second order elliptic problems; Mu  presented a posteriori error estimate for the second order elliptic interface problems; Zheng and Xie  discussed a residual-based a posteriori error estimator for the Stokes problem. There are only one research result for a posterior error estimates for mWG methods. Zhang and Lin  proposed a posteriori error estimator for second order elliptic problems. To our best knowledge, there is no literature on the convergence of adaptive mWG methods.

The goal of this paper is to prove the optimal convergence of an AmWG algorithm for a second order elliptic problem. Different from the mWG originally introduced in , we simplify the mWG as follows: for a weak function , the edge/face term is not independent anymore as we choose , i.e., is obtained through averaging the interior discontinuous variable and then projected through to a one degree lower polynomial space. Compared with interior-penalty discontinuous Galerkin(IPDG), the mWG is stable without choosing a sufficiently large penalty parameter. This simplification ( opposed to in ) brings extra difficulty to the analysis of convergence, as a simple port of the workflow presented in , by decomposing the discontinuous approxmation space into a continuous subspace and its orthogonal complement, will introduce a penalty parameter that is not originally in the mWG discretization.

To conquer this difficulty, by introducing an interpolation operator

onto the Crouzeix-Raviart type nonconforming finite element space , we bound the stabilization term and prove an a posteriori error estimate in the discrete -norm. One main ingredient in the convergence analysis of a standard adaptive procedure is the orthogonality of the error to the finite element space, however, such an orthogonality does not hold true for mWG approximations. Instead, a quasi-orthogonality result is established. Hu and Xu  defined a canonical interpolation operator for the lowest Crouzeix-Raviart type nonconforming finite element space and established the quasi-orthogonality property for both the velocity and the pressure in the Stokes problem. The main observation is that the modify weak gradient for a function is equal to the elementwise gradient of the interpolate , namely and we can derive the desired quasi-orthogonality property for the linear mWG.

Another key component to establish the optimality of the adaptive algorithm is the localized discrete upper bound for the a posteriori error estimator. By using a prolongation operator introduced in , we are able to derive the discrete reliability and use it to prove the optimality of the convergent algorithm.

For the a posteriori error analysis of mWG approximations, we mainly follow the Bonito and Nochetto’s work  and Chen, Wang and Ye’s work . For the analysis of the convergence and the optimality of adaptive procedure, we mainly use the Hu and Xu’s work  and Huang and Xu’s work . We do not claim any originality on the proof of convergence and optimality. Instead, the main contribution of this paper is to bound the stabilization term by the element-wise residual and flux jump, as well as to establish a quasi-orthogonality and a discrete upper bound.

The rest of this paper is organized as follows. In Section 2, the definitions of weak gradient and discrete weak gradient are introduced, as well as the modified weak Galerkin finite element spaces and the corresponding bilinear form . In Section 3, a residual-type error estimator is constructed, together with its reliability and efficiency shown. In Section 4, we introduce an adaptive modified weak Galerkin method (AmWG) and prove its convergence and optimality.

## 2. Notations and preliminaries

The goal of this section is to present the modified weak Galerkin (mWG) formulation for (1.1). First, the standard weak Galerkin method is reviewed, then an mWG finite element space and the discretization thereof are introduced.

### 2.1. Weak Galerkin Methods

Given a polygonal/polyhedral element with boundary , the notation defines a weak function on such that and . Subsequently, the weak function space on is defined as

 W(K)={v={v0,vb}:v0∈L2(K),vb∈L2(∂K)}.

Let be the set of polynomials on with degree no more than (). The discrete weak gradient operator is defined for polynomial functions. With slight abuse of notation, when no ambiguity arises, we shall denote as , which should be clear from the context.

###### Definition 2.1 ([27, Definition 2.1]).

The discrete weak gradient operator, denoted by , is defined as the unique polynomial satisfying the following equation

 (2.1) (∇w,K,ℓv,q)K=−(v0,∇⋅q)K+⟨vb,q⋅n⟩∂K,∀q∈(Pℓ−1(K))d.

In the definition above, is the outward normal direction to , is the -inner product of and , and is the -inner produce of and . The differential operators involved are well defined when restricted to one element. In the context of the gradient operator defined across multiple elements in , is the elementwise gradient, i.e., in the element .

In the rest of the paper, we still restrict ourselves to a shape-regular triangulation of . denotes the set of all the edges or faces in , and is the set of all the interior edges or faces. Denote the -dimensional Lebesgue measure. For , its associated patch as . For a set , its associated patch element patch is .

Given a positive integer , the -th order weak Galerkin finite element space on is defined as follows:

 (2.2) VWG(T):={v={v0,vb}:v0|τ∈Pℓ(τ),vb|e∈Pℓ−1(e),e∈E,τ∈T},

and the one with zero boundary condition:

 (2.3) VWG0(T)={v:v∈VWG(T),vb=0 on ∂Ω}.

Note that is single-valued on each .

For , the bilinear form is defined as

 (2.4) aWG(v,w):=∑τ∈T(A∇wv,∇ww)τ+∑τ∈Th−1τ⟨Qbv0−vb,Qbw0−wb⟩∂τ,

where the weak gradient , and is the projection from to on an . Again when it is clear from the context, then we shall omit the degree of the polynomial involved.

The weak Galerkin discretization is: to find a satisfying

 (2.5) aWG(uh,v)=(f,v0),∀v={v0,vb}∈VWG0(T).

The WG discretization (2.5) is well-posed as defines an inner product on the space .

We include a short argument here for completion and the convenience of our readers to show that defines a norm in . Assume that , then on every element and on every , then by definition (2.1)

 0=(∇wv0,∇v0)τ = −(v0,Δv0)τ+⟨vb,∇v0⋅n⟩∂τ = −(v0,Δv0)τ+⟨Qbv0,∇v0⋅n⟩∂τ = −(v0,Δv0)τ+⟨v0,∇v0⋅n⟩∂τ = (∇v0,∇v0)τ,

which implies on . As a result, is constant on , and on as well as on , thus . Therefore, defines a norm and defines an inner product on the space .

The proof above reveals that the boundary part can be set as a polynomial with degree one less than that of the interior one since its presence is only in

. The reduction is beneficial as the interior degrees of freedom can be further eliminated, and a linear algebraic system of smaller size is solved.

### 2.2. Modified weak Galerkin finite element

Let be the common edge/face shared by two elements , and denote by . We assume that globally each

is associated with a fixed unit normal vector

. When , we can get ’s tangential vector by , which is obtained by rotating clockwisely by . Without loss of generality, the is assumed pointing from to . Denote by and the outer unit normal vectors with respect to and , respectively. For a smooth enough scalar function , we define its average and jump on by

 {w}e=(w|τ1+w|τ2)/2,for  e∈Eint; [[w]]\raisebox−2.0pt$e$=w|τ1−w|τ2,for  e∈Eint; [[w]]\raisebox−2.0pt$e$={w}e=w,for  e⊂∂Ω.

Similarly, for an admissible vector function , we have

 {w}e=(w|τ1+w|τ2)/2, for  e∈Eint; [[w⋅ne]]e=w|τ1⋅n1+w|τ2⋅n2, for  e∈Eint; {w}e=w,[[w⋅ne]]\raisebox−2.0pt$e$=w⋅ne, for  e⊂∂Ω.

In the definitions above, and ’s values on are defined for those spaces which yield a well-defined trace, respectively.

Specifically, in discretization (2.5), is chosen to be in the discontinuous polynomial space

 VDG(T):={v:v|τ∈Pℓ(τ),τ∈T},

and the edge/face term . Now a weak function is , which offers a continuous embedding of . The modified weak Galerkin finite element space for is then defined as

 (2.6) V(T):={v={v0,vb}:v0∈VDG(T);vb:=Qb{v0},e∈E},

and

 (2.7) V0(T)={v:v∈V(T),Qb{v0}=0 on ∂Ω}.
###### Remark 2.2.

The original modified weak Galerkin function space was introduced in . The edge/face term is chosen as in , while here to match the reduced-order weak Galerkin scheme in .

The definition of weak gradient is modified accordingly as follows.

###### Definition 2.3.

The modified weak gradient operator acting on any , denoted by on , is defined as the unique polynomial in such that its inner product with any satisfies the following equation:

 (2.8) (∇w,τv,q)τ=−(v,∇⋅q)τ+⟨Qb{v},q⋅n⟩∂τ.

When the triangulation is clear from the context, we shall abbreviate the notation as .

Let be the common edge/face shared by . By using the DG space in the WG bilinear form (2.4), one can simplify the stabilization term on for , henceforth the same argument applies to other edges/faces.

 ⟨Qb(v−{v}),Qb(w−{w})⟩e =⟨Qb(vτ1−vτ1+vτ22),Qb(wτ1−wτ1+wτ22)⟩e =⟨Qb(vτ1−vτ22),Qb(wτ1−wτ22)⟩e =14⟨Qb[[v]],Qb[[w]]⟩e.

Consequently, for , the associated bilinear form is defined as

 (2.9) aT(v,w) := ∑τ∈T(A∇wv,∇ww)τ+4∑τ∈Th−1τ⟨Qb(v−{v}),Qb(w−{v})⟩∂τ = ∑τ∈T(A∇wv,∇ww)τ+∑τ∈Th−1τ⟨Qb[[v]],Qb[[w]]⟩∂τ.

A modified weak Galerkin (mWG) approximation is then to seek satisfying

 (2.10) aT(uT,v)=(f,v),∀v∈V0(T).

It is clear that the modified weak Galerkin finite element scheme (2.10) is also well-posed since the problem is solved in a subspace (the embedding of the DG space) of the original WG space.

###### Remark 2.4.

In the classic IPDG formulation (see e.g., ), the bilinear form is

 (2.11) aIPDGT(vT,wT) = ∑τ∈T(A∇hvT,∇hwT)τ−∑e∈E⟨[[vT]],{A∇hwT⋅n}⟩e −∑e∈E⟨[[wT]],{A∇hvT⋅n}⟩e+∑τ∈Tμh−1τ⟨[[vT]],[[wT]]⟩∂τ.

In comparison, in (2.9), using the definition of the weak gradient, several times of the integration by parts, and the assumption that is a piecewise constant, an equivalent bilinear form of (2.9) reads

 (2.12) aT(vT,wT) = ∑τ∈T(A∇hvT,∇hwT)τ−∑e∈E⟨[[wT]],{A∇wvT⋅n}⟩e −∑e∈E⟨[[vT]],{A∇hwT⋅n}⟩e+∑τ∈Th−1τ⟨Qb[[vT]],Qb[[wT]]⟩∂τ.

One in (2.11) is replaced by in (2.12). The major difference is that the mWG (2.10) is automatically coercive and continuous as the bilinear form induces a norm. While in IPDG (2.11), the penalty parameter being sufficiently large is a necessary condition to achieve the coercivity both theoretically and numerically .

The following quantity is well-defined for and defines a mesh dependent norm on (see )

 (2.13)

With slightly abuse of notation, we are interested in estimating the following error using computable quantities and designing a convergent adaptive algorithm to reduce its magnitude successively:

 (2.14)

## 3. A Posteriori Error Analysis

In this section, we shall prove the reliability and the efficiency of a residual-type error estimator. For , we define the element-wise error estimator as

 (3.2) η2(∇wvT,τ):=h2τA−1τ∥R(∇wvT)∥20,τ

where and for . The element residual is defined as

 R(∇wvT)=f+∇⋅(A∇wvT),

and the normal jump of the weak flux is defined as

 Jn,e(A∇wvT):={[[A∇wvT⋅ne]]e,if  e∈Eint0,otherwise.

For the tangential jumps, when :

 Jt,e(∇wvT):={[[∇wvT⋅te]]\raisebox−2.0pt$e$,if  e∈Eint0,otherwise,

and when

 Jt,e(∇wvT):={[[∇wvT×ne]]\raisebox−2.0pt$e$,if  e∈Eint0,otherwise.

Then the error estimator for the set is defined as

 (3.3) η2(∇wvT,M)=∑τ∈Mη2(∇wvT,τ).

In computing the error estimator, only the information of is used. A further remark on why we use as the notation for the error estimator will be later explained in Remark 4.1 in the context of the convergence analysis.

### 3.1. A space decomposition

In this section, a space decomposition is first introduced to bridge the past result of the full degree jump to with one less degree. Then, similar to , a partial orthogonality for the mWG approximation is introduced to enable the insertion of continuous interpolants to prove the reliability.

For consisting of discontinuous degree :

 (3.4) Vc(T)⊂Vnc(T)⊂VDG(T),

where is the continuous Lagrange finite element space. is a subspace of so that their quotient space is endowed with the induced topology of the seminorm . When , is simply the famous Crouzeix-Raviart finite element space . When , is a generalization of the Crouzeix-Raviart type nonconforming finite element space , which can be viewed a special case of nonconforming virtual element  restricted on triangulations:

 (3.5) Vnc(T):=VDG(T)∩H1,nc(T),

where

 (3.6) H1,nc(T)={v∈∏τ∈TH1(τ):∫e[[v]]qds=0∀q∈Pℓ−1(e),∀e∈E}.

Unlike the virtual element space, where the shape function may not be polynomials, it is shown in [12, 10] that the space can be obtained from by adding locally supported polynomial bases. More specifically, in , the authors constructed a set of nodal basis for a strict subspace of while satisfying the continuity constraint (3.6). However, to serve the purpose of this paper, is to bridge the proofs, it only suffices to know the existence of a set of unisolvent degrees of freedom, while not explicitly construct the dual basis of it. We refer the reader to [2, Section 3.2] for the degrees of freedom for general polytopes which includes the case of triangulations () in this paper.

When restricted a smaller subspace , the weak gradient coincides with the piecewise gradient, and the stabilization is vanished which leads to the following partial orthogonality.

###### Lemma 3.1.

Let and be the solutions of (1.1) and (2.10), respectively, then

 (3.7)
###### Proof.

It follows from that . In addition, for functions when , thus in , which further implies . Lastly since implies that for any , as a result, and the lemma follows. ∎

The following interpolation operator to non-comforming space will play an important role for the analysis.

###### Lemma 3.2.

There exist an interpolation operator , which is locally defined and a projection, as well as a constant depending only on the shape regularity of such that for all the following inequality hold: for ,

 (3.8) ∥∥Da(vT−ITvT)∥∥2τ≲∑e∈Eint(ω(τ))h1−2|a|τ∥∥Qb[[vT]]∥∥2e,∀vT∈VDG(T).
###### Proof.

The proof follows a similar argument as the one in [3, Lemma 6.6], and is presented here for completion. Denote the set of the degrees of freedom functionals in as , and the nodal bases set is corresponding to . Let , where is either a single element or , then we consider the following projection obtained through

 (3.9) (vT−ΠivT,w)ωi=0,∀w∈Vnc(ωi):=Πτ∈ωiPℓ(τ)∩H1,nc(T),

and then the interpolation is defined as:

 (3.10) ITvT(x):=∑i∈Λγi(ΠivT)ϕi(x).

If , we have locally on , thus , and is a projection.

Moreover, implies that , as well as by (3.9), hence by a scaling argument and the equivalence of norms on a finite dimensional space, we have

 (3.11) ∑τ∈ωi∥Da(vT−ΠivT)∥2τ≲∑e∈Eint(ωi)h1−2|a|τ∥Qb[[vT]]∥2e,

where we note that if a nodal basis has support , both sides of the inequality above will be 0, as the projection (3.9) modulo out the element bubbles.

Next, consider on a , and all local degrees of freedom of which the nodal basis having support overlaps with , we have on

 (3.12) vT−IT=(vT−Π1vT)−Nτ∑i=1γi(Π1vT−ΠivT)ϕi.

On , is either defined on (see [2, Section 3.2])

 γi(vT)=1|e|∫evTpds,∀p∈Pℓ−1(e),

or defined on ()

 γi(vT)=1|τ|∫τvTpds,∀p∈Pℓ−2(τ),

and in both cases, we have by a simple scaling check and by (3.11)

 ∥∥γi(Π1vT−ΠivT)∥∥τ≲∥∥Π1vT−ΠivT∥∥τ≲∑e∈Eint(ωi)hτ∥Qb[[vT]]∥2e.

Lastly, combining the estimate above with (3.12), using the fact that , and estimate (3.11) yields the desired estimate. ∎

### 3.2. Reliability

In this section, the global reliability of the estimator in (3.3) is to be shown. The difference between the weak gradient and classical gradient can be controlled by the jump term which becomes a handy tool in our analysis.

###### Lemma 3.3.

For , it holds for any

 ∥∥A1/2(∇wvT−∇hvT)∥∥2M≲∑e∈E(M)Amaxeh−1τ∥Qb[[vT]]∥2e.
###### Proof.

The proof follows from [43, Lemma 2.1] with coefficient added and a more localized version is presented. On , let , applying integration by parts on and the weak gradient definition (2.8) on , we have

 ∫τA(∇wvT−∇hvT)⋅qdx = ∫∂τ(Qb{vT}−vT)Aτq⋅nds = ∫∂τQb({vT}−vT)Aτq⋅nds =

where the plus or minus sign depends on whether the outward normal for coincides with the globally defined normal for this edge/face. Now by a standard trace inequality and an inverse inequality, we have

 ∥∥A1/2(∇wvT−∇hvT)∥∥2τ ≲∑e⊂∂τA1/2τh−1/2τ∥Qb[[vT]]∥eA−1/2τh1/2τ∥Aτq⋅n∥e ≲(∑e⊂∂τAmaxeh−1τ∥Qb[[vT]]∥2e)1/2∥∥A1/2(∇wvT−∇hvT)∥∥τ,

then the desired result follows by canceling a and summing up the element-wise estimate for every . ∎

Next we bound the stabilization term by the element-wise residual and the normal jump of the weak flux. We note that in , though focusing on a different model problem, the -weighted solution jump can be used as the sole error estimator to guarantee the reliability and efficiency up to data oscillation. Nevertheless, the motivation here is to change the dependence of the error indicators on the local mesh size from to , so that a contraction property of the estimator can be proved in Section 4.2 without any saturation assumptions, which is one of the keys to show the convergence of an adaptive algorithm.

###### Lemma 3.4.

Let be the weak solution of (1.1) and be the solution to (2.10), we have

 (3.13) ≲ (∑τ∈Th2τA−1τ∥R(∇wuT)∥20,τ +∑e∈Ehτ(Amaxe)−1∥Jn,e(A∇wuT)∥2e),

where the constant depends on the shape regularity of .

###### Proof.

For simplicity, we denote in the proof. Let that be the nonconforming interpolation defined in Lemma 3.2, using the definition of the modified weak derivative (2.8) element-wisely, we have

 ∑e∈EAeh−1τ∥Qb[[uT]]∥2e = = (f,uT−ITuT)T−(A∇wuT,∇w(uT−ITuT))T = (f,uT−ITuT)T+(∇h⋅(A∇wuT),uT−ITuT)T −∑τ∈Th⟨(A∇wuT)⋅n,{uT−ITuT}⟩∂τ = (f+∇h⋅(A∇wuT),uT−ITuT)T −∑e∈E⟨[[(A∇wuT)⋅n]],{uT−ITuT}⟩e.

By Cauchy-Schwarz inequality, we have

 ∑e∈EAeh−1τ∥Qb[[uT]]∥2e⩽(∑τ∈Th2τA−1τ∥f+∇⋅(A∇wuT)∥2τ)1/2(∑τ∈Th−2τ∥A1/2(uT−ITuT)∥2τ)1/2+(∑e∈EhτA−1e∥∥[[A∇wuT⋅n]]