A family of interior-penalized weak Galerkin methods for second-order elliptic equations

Abstract: Interior-penalized weak Galerkin (IPWG) finite element methods are proposed and analyzed for solving second order elliptic equations. The new methods employ the element (Pk,Pk,RT k), with dimensions of space d = 2, 3, and the optimal a priori error estimates in discrete H1-norm and L2-norm are established. Moreover, provided enough smoothness of the exact solution, superconvergence in H1 and L2 norms can be derived. Some numerical experiments are presented to demonstrate flexibility, effectiveness and reliability of the IPWG methods. In the experiments, the convergence rates of the IPWG methods are optimal in L2-norm, while they are suboptimal for NIPG and IIPG if the polynomial degree is even.


Introduction
The weak Galerkin finite element method (WG-FEM) is a class of non-conforming FEMs, in which the differential operators, such as gradient, divergence, etc., are approximated by their weak forms. The initial WG method based on a simplicial mesh was introduced by Wang and Ye [25] for solving general second-order elliptic equations. A new WG method with a stabilizer based on polygonal/polyhedral elements was proposed by Mu et al. [13,14] , which brings great convenience and flexibility in mesh generation and assembly of the stiffness matrix. This great progress allowed WG-FEMs to flourish in solving many important PDEs, such as interface problems [11,17,21], Helmholtz equation [12,28], Maxwell equation [16,19,22], biharmonic equation [23], linear elasticity problem [24], Darcy-Stokes and Darcy flow [4,8,27], stochastic equations [33,34], parabolic equations [3,7], etc.
Due to the great flexibility in an arbitrary finite element shape and variety of weak finite element spaces, many interesting variations of WG methods arose, such as mixed weak Galerkin [26], least-squares-based weak Galerkin (LSWG) [15], stabilizer-free weak Galerkin [31], modified weak Galerkin (MWG) [30] and over-penalized weak Galerkin (OPWG) [9,10], etc. A WG method with a modified stabilizer, named after over-relaxed WG method [20,21], exhibits great ability in handling elliptic (interface) problems with low regularity, and is proven to have super-closeness properties in solving second-order elliptic problems [29]. By using a Schur complement formulation, the degrees of freedom of weak Galerkin method with boundary continuity [32] can be reduced to the level of the continuous Galerkin (CG) methods. For second-order elliptic equations, each weak function denoted by v = (v 0 , v b ) consists of two parts: the interior part v 0 ∈ L 2 (K) is defined in each element K ∈ T h and the boundary part v b ∈ L 2 (∂K) is defined on the edges/faces of an element K, where T h is a finite element partition. The gradient operator is approximated by a discrete weak gradient consisting of piecewise vector-valued polynomials defined on each element of a finite element partition. Two kinds of elements are commonly employed. The first kind consists of, see [25] for example, Raviart-Thomas (RT ) element, denoted by (P k , P k , RT k ), and Brezzi-Douglas-Marini (BDM) element, denoted by (P k , P k+1 , [P k+1 ] d ), where the three components represent polynomial spaces for two parts of a weak function and its weak gradient, respectively. The second kind, including (P k , P k , [P k−1 ] d ) and (P k , P k−1 , [P k−1 ] d ) with an integer k ≥ 1, is most widely used in practical applications. Compared to RT and BDM elements, a stabilizer is required to control the discontinuity of interior part and edge part of weak function within each element of finite element partition consisting of an arbitrary polygon or polyhedron shapes.
By taking completely discontinuous approximation functions as shape functions, the WG method inherits many advantages of the DG method, like complex geometries and various boundary conditions can be treated. Traditionally, WG-FEMs are based on the weak functions with the single-valued boundary part on each interior edge of a partition. We proposed and analyzed the OPWG method based on RT element for second-order elliptic problems [10], where the shape function along the interior edges/faces are double-valued. The main idea there is to integrate WG-FEM and DG-FEM, expecting to inherit more advantages of the DG methods. In [9], an OPWG method, based on element (P k , P k , [P k−1 ] d ) and (P k , P k−1 , [P k−1 ] d ), with a stabilizer term is proposed. First, the defect of fastincreasing condition numbers is perfectly settled by using a simple block-diagonal preconditioner. The reduction of the degree of freedoms by a Schur complement algorithm makes the OPWG method attractive and comparable to the HDG and the primal DG methods (cf [9]). Besides, it makes an adaptive approximation or hp-WG possible and the OPWG method is promising in parallel computing and approximating discontinuous solutions of PDEs.
In this paper, we propose a family of interior-penalized weak Galerkin (IPWG) methods based on element (P k , P k , RT k ). First, the ill-conditioned system due to the over-penalization in [9,10] is avoided by introducing new terms in the variational formulation. Second, numerical experiments in Section 5 confirm that the converge rates in L 2 -norm of all IPWG methods are optimal with respect to the mesh size h, while for non-symmetric interior penalty Galerkin (NIPG) or incomplete interior penalty Galerkin (IIPG), they are suboptimal if the polynomial degree is even [18]. Here, we first introduce the second-order elliptic problem with Dirichlet boundary condition: where Ω is a polygonal or polyhedral domain in R d (d = 2, 3), f ∈ L 2 (Ω) and A is a symmetric positive definite matrix-valued function on Ω, i.e., there exist two positive numbers A 0 , A 1 > 0 such that where ξ is a column vector and ξ t is the transpose of ξ.
Throughout the paper, we will use the standard notations for Sobolev spaces and norms [2]. Specifically, for an open bounded domain D ⊂ R d , d = 2, 3, we use · W m,p (Ω) and | · | W m,p (D) to denote the norm and seminorm in W m,p (D) for m ≥ 0. When p = 2, we adopt the notations · H m (D) and | · | H m (D) . In addition, the space H 0 (D) coincides with L 2 (D), for which the inner product is denoted by (·, ·) D or (·, ·) whenever there is no confusion. We introduce the broken Sobolev space for any m ≥ 0, The rest of this paper is organized as follows. In Section 2, we recall the definition of the weak gradient and define jumps of weak functions, as well as some projection operators. In Section 3, we introduce the IPWG formulations. In Section 4, the a priori error estimates in energy norm and L 2 -norm are derived. In the last section, some numerical experiments are presented.

Preliminaries
Let T h be a shape regular triangulation of the domain Ω, which is required in classic finite element analysis, see [2]. Note that we don't require that the mesh T h is quasi-uniform [2, (4.4.13)]. Let E h be the set of all edges (or faces), E I h the set of all interior edges and E B h the set of all boundary edges. We denote by h K the diameter of K ∈ T h and denote by h = max K∈T h h K mesh size for T h . Let P k (K) be the polynomial space of degree no more than k on K ∈ T h . Similarly, P k (e) denotes the polynomial space of degree no more than k on e ∈ E h .
We define the local weak function space on Definition 1. The weak gradient operator, denoted by ∇ w , is a linear operator from W(K) to [H 1 (K)] d , such that for any v ∈ W(K) and q ∈ [H 1 (K)] d , the following holds where n is the outward normal direction to ∂K, (·, ·) K is the inner product in L 2 (K), and ·, · ∂K is the inner product in L 2 (∂K).
We define a discrete weak gradient operator to approximate ∇ w in a space of vector polynomials.
Definition 2. The discrete weak gradient operator, denoted by ∇ w,k,K , is a linear operator from W(K) to V k (K), such that for any v ∈ W(K) and q ∈ V k (K), the following holds where V k (K) is a subspace of vector-valued polynomials of degree no more than k in element K.
For simplicity of notation, we always denote by ∇ w the discrete weak gradient operator ∇ w,k,K . In what follows, we define the discrete weak gradient space as [1,10].
Define the weak Galerkin finite element space associated with T h as where v b is a double-valued function on each interior edge/face. For convenience, We refer the element we use as (P k , P k , RT k ).
Remark 2.1. We note that all results developed in this paper are also valid for (P k , P k+1 , BDM k ), i.e., weak function space consisting of (P k (K), P k+1 (e)) and the weak gradient space consisting of ] or [25, Section 5] for more details.
We introduce a set of normal vectors of E h as follows: D = {n e : n e is unit and normal to e, e ∈ E h }.
If e ∈ ∂Ω, we set n e to be exterior to the domain. Let K e i , i = 1, 2 be the elements adjacent to e, we denote by [v b ] |e the jump of v b on the edge e ∈ E I h whose normal vector n e is oriented from K e 1 to K e 2 : and denote by {v b } the average of v b : We also extend the definition of the jump and average to the edges that belong to the domain boundary, h . In addition, we denote by |e| the length or area of e and by h e the diameter of edge or flat face e ∈ E h . For a shape regular mesh, one can see that there exist constants κ and ρ e such that κh K ≤ h e and ρ e h d−1 To investigate the approximate properties of WG finite element space V h and weak gradient space RT k , we introduce the following three L 2 projections: We combine Q 0 and Q b by writing Q h = (Q 0 , Q b ). It is easy to see that (see [25]): To complete error analysis, we need a divergence conforming projection operator Π h , which satisfies the following property [1, Section 2.5.2]: for any τ ∈ H(div, Ω); and on each element K ∈ T h , one has Π h (τ |K ) ∈ RT k (K) and the following commutative identity holds

The interior-penalized weak Galerkin scheme
For any w, v ∈ V h , we define the following bilinear forms where σ e is the so-called penalty parameter depending on the mesh T h , and β is strictly positive depending on the spatial dimension d. In what follows, we take = −1, 0, 1, and denote σ = min e∈E h σ e . A class of interior-penalized weak Galerkin approximations of (1.1)-(1.2) can be obtained by Remark 3.1. In the case of = −1, a (·, ·) is symmetry, and we will see that this method converges if β(d − 1) ≥ 1 and σ large enough; In the case of = 0, the method converges under the same conditions as the case of = −1; In the case of = 1, the method converges for any strictly positive values of σ e .
To justify the well-posedness of (3.1), we equip V h with the norm [10]: The following classic trace inequality can be derived from a trace theorem and a scaling argument, see, e.g., [2,5].
Lemma 3.1 (Trace inequality). Suppose that the triangulation T h is shape regular, then there exists a constant C such that for any K ∈ T h and edge e ∈ ∂K, for any φ ∈ H 1 (K), it holds

3)
and for any ψ ∈ P k (K), it holds We mimic the procedure of [18, Section 2.7.1] for interior-penalized DG method to prove the following coercivity.
Here, C t is from the trace inequality and n 0 is the maximum number of neighbors an element could have (n 0 = 3 if d = 2 and n 0 = 4 if d = 3).
Proof. By the definition of a (·, ·), for each v ∈ V h , there holds First, in the case of = 1, by the definition of energy norm (3.2), we have Next, we focus on the cases of = −1 and 0. Without loss of generality, we assume that edge e ∈ E h is shared by two elements K e i , i = 1, 2. Trace inequality leads to Using |e| ≤ h d−1 where we have assumed β(d − 1) ≥ 1 and h ≤ 1. Summing over all edges and using Young's inequality with δ > 0, we obtain (3.6) Thus, a (·, ·) becomes By taking δ = 1 2 for = −1 and δ = 1 for = 0, then, for σ e large enough (for example, σ e ≥ The following uniqueness is a direct consequence of Proposition 3.1.

Proposition 3.2 (Uniqueness).
Under the conditions given in Proposition 3.1, the interior-penalized weak Galerkin methods (3.1) have one and only one solution.

Error estimates
The goal of this section is to establish error estimates for the IPWG methods (3.1). The error is measured in two norms: the triple-bar norm as defined in (3.2) and the standard L 2 norm.
We first state two approximation properties of the operator Π h , which can be proved by combining [ Next, we state an important result for the projection operator Π h .

Lemma 4.2.
Let τ ∈ H(div, Ω) be a smooth vector-valued function and Π h be the local projection operator defined in Section 2. Then, the following identify holds true In what follows, we define e h = (e 0 , e b ) := Q h u − u h .
which can be rewritten as (4.5)

By adding
to both sides of (4.5) and noticing that [Q b u] |e = 0 for e ∈ E I h and [Q b u] |e = Q b (g |e ) for boundary edge e ∈ E B h , we obtain (4.6) Subtracting (3.1) from (4.6) leads to which is the error equation for the IPWG approximations (3.1).

Lemma 4.3.
Let T h be a shape regular partition. Then for any u ∈ H s (K) and v = (v 0 , v b ) ∈ V h , we have for any s ≥ 2 Note that the estimate (4.9) holds if β(d − 1) ≥ 1.
Proof. The estimate (4.8) is a direct result of Lemma 4.1. We now estimate (4.9). Assume that β(d − 1) ≥ 1, which is given in Proposition 3.1. By using the trace inequality and recalling |e| ≤ h d−1 e , we obtain Let K e 1 and K e 2 be the adjacent elements containing edge e ∈ E h . Then, it follows that Summing over all edges leads to (4.9).
The following theorem follows from Lemma 4.3.
In the rest of the section, we shall derive an error estimate in L 2 -norm for the IPWG methods (3.1) by using a duality argument. Suppose the dual problem: Find a solution w ∈ H 1 0 (Ω) satisfying −∇ · (A∇w) = e 0 , in Ω, (4.11) has the usual H 2 -regularity, i.e., there exists a constant C > 0 such that w H 2 (Ω) ≤ C e 0 L 2 (Ω) .  Note that for = −1, the a priori estimate (4.13) is optimal if k ≥ 1; For = 0, 1, the estimate is optimal if k ≥ 1 and β(d − 1) − 1 ≥ 2. We now estimate the second term on the right hand side. By taking and noticing [Q b w] = 0 on all edges e ∈ E h , we obtain from the error equation (4.7) The second term on the right hand side clearly is zero by noticing that the coefficient A is piecewise constant and Q h is an L 2 projection. Next, we turn to the first term. By using Π h A∇u ∈ H(div, Ω) and (2.4), an integration by parts leads to

(4.16)
It is easy to see that T = 0 when = −1, hence, one has to estimate T if = 0 or 1. The same technique as in (3.6) leads to (4.17) Hence, we can obtain from (4.16), (4.17) and (4.1) Substituting the above estimate into (4.15) and using the continuity of a and (4.14), we arrive at which completes the proof.

Numerical experiments
In this section, we present a series of computational examples to numerically investigate the asymptotic convergence behaviour of the proposed IPWG methods.
We take Ω = (0, 1) 2 , and multilevel uniform triangular meshes are employed. The meshes are generated in the following way. First, we partition the square domain into N × N subsquares uniformly; then we divide each subsquare into two triangles by the diagonal line with a negative slope, completing the construction of uniformly refined triangular meshes, see Figure 1. Set the mesh size h = 1 N . For comparison between the IPWG method and the original WG method, the relative errors in · and L 2 -norms, defined as of the WG method with the element (P k , P k , RT k ) in Table 1, from which one can see that the WG method converges optimally. It follows from Theorem 4.1, all IPWG solutions with = −1, 0, 1 converge optimally with rates O(h k ) in · -norm. Since the diffusion coefficient A is constant and the domain Ω is convex, the duality problem possesses the H 2 -regularity. Thus, from Theorem 4.2, the IPWG solution with = −1 should converge optimally given by O(h k+1 ) in L 2 -norm; and converge sub-optimally in L 2 -norm given by O(h k ) for = 0, 1 without over-penalization.
The relative errors between the exact solution and IPWG ( = −1, β = 1) solution in · -and L 2 -norm are listed in Table 2. Note that, for = 0, 1, the IPWG methods produce the same numerical results as in the case of = −1, see Table 2. It can be seen from Tables 2 and 3 that errors of the IPWG ( = −1, 0, 1) solutions are numerically the same as those of the WG method, and all the convergence rates are optimal without any over-penalization. The numerical results are fairly in agreement with the theoretical estimates given by Theorems 4.1-4.2. In addition, one can see that from Table 3, the IPWG method with = 1, σ e = 0 also has optimal rates while there is no theoretical proof.
where r = x 2 + y 2 and α ∈ (0, 1]. As we known where is any small, but positive number. In numerical experiments, we set α = 0.5. The problem is tested on the uniform and locally refined meshes to investigate the convergence behaviour of the IPWG method, respectively. We set β such that β(d − 1) − 1 = 0, i.e., β = 1 for d = 2. From Table 4, we observe all cases of the IPWG method with = −1, 0, 1 converge optimally at O(h α ) in · -norm and O(h 1+α ) in L 2 -norm, or equivalently where the Dof denotes the degrees of freedom. Next, we employ a series of locally-refined meshes generated by Gmsh [6], see Figure 2, to illustrate the convergence behaviour of the IPWG methods for Example 2. For simplicity, we set = −1, β = 1, and we set σ = 1, 8, 16 for k = 0, 1, 2, respectively. The convergence rates of the relative error respect to Dof 1/2 are plotted in Figure 3. One can see that the IPWG method on the locally-refined mesh performs better than the uniform mesh in the sense that higher convergence rates in · and L 2 -norms are observed.

Conclusions
We have presented a family of IPWG methods for the second-order elliptic problems and established optimal a priori error estimates. The superconvergence of the new method is revealed. The new method has many in common with the DG method: they share a similar numerical formulation; they have the same conditions for the uniqueness, see Proposition 3.2; and the penalty parameter σ e which, in practice, needs to be adjusted according to the mesh used because there is still no explicit formula for it.
Yet, it is worthy to note that the new method possesses its own merits. First, the IPWG method exhibits more stability on the polynomial order than IPDG. In detail, numerical experiments for Example 1 state the convergence rates in H 1 and L 2 norms of all IPWG methods are optimal (in fact they are superconvergent) without the over-penalization, no matter the polynomial degree is even or odd. In contrast, for the NIPG and IIPG methods, the converge rates are suboptimal if the polynomial degree is even [18]. Second, Theorem 4.2 reveals superconvergence of the IPWG method if the exact solution is smooth enough. To the end, the degrees of freedom defined on the interior of each element can be reduced by a Schur complement formulation, see [9], and only the unknowns on each edge/face of the mesh left.