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

1.   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}\in L^{2}(K)$ is defined in each element $K\in \mathcal{T}_{h}$ and the boundary part $v_{b}\in L^{2}(\partial K)$ is defined on the edges/faces of an element $K$, where $\mathcal{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 ($\mathcal{RT}$) element, denoted by $(\mathbb{P}_k, \mathbb{P}_k, \mathcal{RT}_{k})$, and Brezzi-Douglas-Marini ($\mathcal{BDM}$) element, denoted by $(\mathbb{P}_k, \mathbb{P}_{k+1}, [\mathbb{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 $(\mathbb{P}_k, \mathbb{P}_{k}, [\mathbb{P}_{k-1}]^d)$ and $(\mathbb{P}_k, \mathbb{P}_{k-1}, [\mathbb{P}_{k-1}]^d)$ with an integer $k\geq 1$, is most widely used in practical applications. Compared to $\mathcal{RT}$ and $\mathcal{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 $\mathcal{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 $(\mathbb{P}_k, \mathbb{P}_{k}, [\mathbb{P}_{k-1}]^d)$ and $(\mathbb{P}_k, \mathbb{P}_{k-1}, [\mathbb{P}_{k-1}]^d)$, with a stabilizer term is proposed. First, the defect of fast-increasing 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 $(\mathbb{P}_k, \mathbb{P}_{k}, \mathcal{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 $\Omega$ is a polygonal or polyhedral domain in $\mathbb{R}^d\ (d = 2, 3)$, $f\in L^{2}(\Omega)$ and $A$ is a symmetric positive definite matrix-valued function on $\Omega$, i.e., there exist two positive numbers $A_0, \ A_1 > 0$ such that

where $\xi$ is a column vector and $\xi^t$ is the transpose of $\xi$.

Throughout the paper, we will use the standard notations for Sobolev spaces and norms ^[2]. Specifically, for an open bounded domain $D\subset \mathbb{R}^{d}, d = 2, 3$, we use $\|\cdot\|_{W^{m, p}(\Omega)}$ and $|\cdot|_{W^{m, p}(D)}$ to denote the norm and seminorm in $W^{m, p}(D)$ for $m\geq 0$. When $p = 2$, we adopt the notations $\|\cdot\|_{H^{m}(D)}$ and $|\cdot|_{H^{m}(D)}$. In addition, the space $H^{0}(D)$ coincides with $L^{2}(D)$, for which the inner product is denoted by $(\cdot, \cdot)_{D}$ or $(\cdot, \cdot)$ whenever there is no confusion. We introduce the broken Sobolev space for any $m \geq 0$,

equipped with the broken Sobolev norm:

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.

2.   Preliminaries

Let $\mathcal{T}_h$ be a shape regular triangulation of the domain $\Omega$, which is required in classic finite element analysis, see ^[2]. Note that we don't require that the mesh $\mathcal{T}_h$ is quasi-uniform [2,(4.4.13)]. Let $\mathcal{E}_h$ be the set of all edges (or faces), $\mathcal{E}_h^I$ the set of all interior edges and $\mathcal{E}_h^B$ the set of all boundary edges. We denote by $h_K$ the diameter of $K\in\mathcal{T}_h$ and denote by $h = \max_{K\in\mathcal{T}_h}h_{K}$ mesh size for $\mathcal{T}_h$. Let $\mathbb{P}_k(K)$ be the polynomial space of degree no more than $k$ on $K\in \mathcal{T}_h$. Similarly, $\mathbb{P}_k(e)$ denotes the polynomial space of degree no more than $k$ on $e\in \mathcal{E}_h$.

We define the local weak function space on $K\in \mathcal{T}_h$

Definition 1. The weak gradient operator, denoted by $\nabla_{w}$, is a linear operator from $W(K)$ to $[H^{1}(K)]^{d}$, such that for any $v\in W(K)$ and $q\in [H^{1}(K)]^{d}$, the following holds

where $\mathit{\boldsymbol{n}}$ is the outward normal direction to $\partial K$, $(\cdot, \cdot)_{K}$ is the inner product in $L^{2}(K)$, and $\langle \cdot, \cdot\rangle_{\partial K}$ is the inner product in $L^{2}(\partial K)$.

We define a discrete weak gradient operator to approximate $\nabla_{w}$ in a space of vector polynomials.

Definition 2. The discrete weak gradient operator, denoted by $\nabla_{w, k, K}$, is a linear operator from $W(K)$ to $\mathit{\boldsymbol{V}}_k(K)$, such that for any $v\in W(K)$ and $q\in \mathit{\boldsymbol{V}}_k(K)$, the following holds

where $\mathit{\boldsymbol{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 $\nabla_{w}$ the discrete weak gradient operator $\nabla_{w, k, K}$. In what follows, we define the discrete weak gradient space as

where $\mathcal{RT}_k(K)$ is $\mathcal{RT}$ element of order $k$ on $K\in \mathcal{T}_h$, see ^[1,10].

Define the weak Galerkin finite element space associated with $\mathcal{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 $(\mathbb{{P}}_k, \mathbb{{P}}_{k}, \mathcal{{RT}}_k)$.

Remark 2.1. We note that all results developed in this paper are also valid for $(\mathbb{{P}}_k, \mathbb{{P}}_{k+1}, \mathcal{{BDM}}_k)$, i.e., weak function space consisting of $(\mathbb{{P}}_k(K), \mathbb{{P}}_{k+1}(e))$ and the weak gradient space consisting of $\mathcal{BDM}_k(K), k\geq 1$, see [1,(2.3.7)] or [25,Section 5] for more details.

We introduce a set of normal vectors of $\mathcal{E}_h$ as follows:

If $e\in\partial\Omega$, we set $\mathit{\boldsymbol{n}}_{e}$ to be exterior to the domain. Let $K_i^e, i = 1, 2$ be the elements adjacent to $e$, we denote by $[ v_b]_{|e}$ the jump of $v_b$ on the edge $e\in\mathcal{E}_h^I$ whose normal vector $\mathit{\boldsymbol{n}}_{e}$ is oriented from $K_1^e$ to $K_2^e$:

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, i.e., $[v_b]_{|e} = \{v_{b}\}_{|e} = {v_b}_{|e}, \ \forall\ e\in \mathcal{E}_h^B$. In addition, we denote by $|e|$ the length or area of $e$ and by $h_{e}$ the diameter of edge or flat face $e\in \mathcal{E}_h$. For a shape regular mesh, one can see that there exist constants $\kappa$ and $\rho_{e}$ such that $\kappa h_{K}\leq h_{e}$ and $\rho_{e} h_{K}^{d-1} \leq |e| \leq h_{K}^{d-1}, \ e\in \partial K$.

To investigate the approximate properties of WG finite element space $V_h$ and weak gradient space $\mathcal{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 $\Pi_h$, which satisfies the following property [1,Section 2.5.2]: for any $\mathit{\boldsymbol{\tau}}\in H(\mathrm{div}, \Omega)\cap H^{\delta}(\mathcal{T}_h), \ \delta > \frac{1}{2}$, $\Pi_h \mathit{\boldsymbol{\tau}} \in H(\mathrm{div}, \Omega)$; and on each element $K\in \mathcal{T}_h$, one has $\Pi_h (\mathit{\boldsymbol{\tau}}_{|K}) \in \mathcal{RT}_k(K)$ and the following commutative identity holds

which implies

Moreover, there holds for all $m\geq 1$ [1,Proposition 2.5.3]

3.   The interior-penalized weak Galerkin scheme

For any $w, v\in V_h$, we define the following bilinear forms

where $\sigma^{e}$ is the so-called penalty parameter depending on the mesh $\mathcal{T}_h$, and $\beta$ is strictly positive depending on the spatial dimension $d$. In what follows, we take $\epsilon = -1, 0, 1$, and denote $\sigma = \min_{e\in \mathcal{E}_h}\sigma^{e}$.

A class of interior-penalized weak Galerkin approximations of (1.1)-(1.2) can be obtained by finding $u_h = (u_0, u_b) \in V_h$ such that

where

Remark 3.1. In the case of $\epsilon = -1$, $a_{\epsilon}(\cdot, \cdot)$ is symmetry, and we will see that this method converges if $\beta (d-1)\geq 1$ and $\sigma$ large enough; In the case of $\epsilon = 0$, the method converges under the same conditions as the case of $\epsilon = -1$; In the case of $\epsilon = 1$, the method converges for any strictly positive values of $\sigma^{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 $\mathcal{T}_h$ is shape regular, then there exists a constant $C$ such that for any $K\in \mathcal{T}_h$ and edge $e\in \partial K$, for any $\phi\in H^1(K)$, it holds

and for any $\psi\in \mathbb{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.

Proposition 3.1 (Coercivity). For each $v\in V_{h}$, the bilinear form $a_\epsilon(\cdot, \cdot)$ is coercive

with

(1) $\kappa = 1$ if $\epsilon = 1$ and $\sigma^{e} > 0$;

(2) $\kappa = \frac{1}{2}$ if $\epsilon = -1$ or $0$, and $\beta (d-1)\geq 1$ and choosing $\sigma^{e} > 0$ large enough. For instance, one can choose $\sigma^{e}$ as follows

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_{\epsilon}(\cdot, \cdot)$, for each $v\in V_{h}$, there holds

First, in the case of $\epsilon = 1$, by the definition of energy norm (3.2), we have

Next, we focus on the cases of $\epsilon = -1$ and $0$. Without loss of generality, we assume that edge $e\in \mathcal{E}_h$ is shared by two elements $K^{e}_{i}, i = 1, 2$. Trace inequality leads to

Using $|e|\leq h_{K^e_i}^{d-1}$ yields

where we have assumed $\beta (d-1)\geq 1$ and $h\leq 1$. Summing over all edges and using Young's inequality with $\delta > 0$, we obtain

Thus, $a_{\epsilon}(\cdot, \cdot)$ becomes

By taking $\delta = \frac{1}{2}$ for $\epsilon = -1$ and $\delta = 1$ for $\epsilon = 0$, then, for $\sigma^{e}$ large enough (for example, $\sigma^{e}\geq \frac{2C_{t}^{2}A_{1}^{2}n_{0}} {A_{0}}+\frac{1}{2}$ if $\epsilon = -1$ and $\sigma^{e}\geq \frac{C_{t}^{2}A_{1}^{2}n_{0}} {A_{0}}+\frac{1}{2}$ if $\epsilon = 0$), we have the following coercivity result

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.

4.   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 $\Pi_h$, which can be proved by combining [25,Lemma 7.3] with the approximation result in [1,Proposition 2.5.1].

Lemma 4.1. Let $\Pi_h$ be the local projection operator defined in Section 2. For all $u\in H^{s}(\Omega)$ with $s\geq 2$, we have

Next, we state an important result for the projection operator $\Pi_h$.

Lemma 4.2. Let $\mathit{\boldsymbol{\tau}}\in H(\mathrm{div}, \Omega)$ be a smooth vector-valued function and $\Pi_h$ be the local projection operator defined in Section 2. Then, the following identify holds true

Proof. See [10,Lemma 4.2].

In what follows, we define

Testing (1.1) with $v_0$ of $v = (v_0, v_b)\in V_h$ and using (4.3) lead to

which can be rewritten as

By adding

to both sides of (4.5) and noticing that $[Q_bu]_{|e} = 0$ for $e\in \mathcal{E}_h^I$ and $[Q_bu]_{|e} = Q_b(g_{|e})$ for boundary edge $e\in \mathcal{E}_h^B$, we obtain

Subtracting (3.1) from (4.6) leads to

which is the error equation for the IPWG approximations (3.1).

Lemma 4.3. Let $\mathcal{T}_h$ be a shape regular partition. Then for any $u\in H^{s}(K)$ and $v = (v_0, v_b)\in V_h$, we have for any $s\geq 2$

where

Note that the estimate (4.9) holds if $\beta (d-1)\geq 1$.

Proof. The estimate (4.8) is a direct result of Lemma 4.1. We now estimate (4.9). Assume that $\beta (d-1)\geq 1$, which is given in Proposition 3.1. By using the trace inequality and recalling $|e|\leq h_{e}^{d-1}$, we obtain

Let $K^e_1$ and $K^e_2$ be the adjacent elements containing edge $e\in \mathcal{E}_h$. Then, it follows that

Summing over all edges leads to (4.9).

The following theorem follows from Lemma 4.3.

Theorem 4.1. Assume the conditions of Proposition 3.1 hold true and the exact solution $u$ of (1.1)-(1.2) belongs to $H^{s}(\mathcal{T}_h)$ with $s\geq 2$. Let $u_h\in V_h$ be a numerical solution of (3.1). Then, there exits a constant $C > 0$ independent of mesh size $h$ such that

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\in H_0^1(\Omega)$ satisfying

has the usual $H^2$-regularity, i.e., there exists a constant $C > 0$ such that

Theorem 4.2. Under the conditions of Theorem 4.1, and assume that the dual problem (4.11) has the $H^2$-regularity. Then, there exits a constant $C > 0$ independent of the mesh size $h$ such that

Note that for $\epsilon = -1$, the a priori estimate (4.13) is optimal if $k\geq 1$; For $\epsilon = 0, 1$, the estimate is optimal if $k\geq 1$ and $\beta(d-1)-1\geq 2$.

Remark 4.1. One can see that, with sufficient smoothness of $u$ and $f$, the IPWG ($\epsilon = -1$) method possesses a superconvergence with $O(h^{k+2})$. The same superconvergence holds true for IPWG ($\epsilon = 0, 1$) with an over-penalization such that $\beta(d-1)-1\geq 2$.

Proof. Let $w_h = \{w_0, w_b\}\in V_h$ be the IPWG solution of (4.11), then $w_h$ satisfies the weak formulation

and the priori estimate

Thus, we can obtain

We now estimate the second term on the right hand side. By taking

and noticing $[Q_bw] = 0$ on all edges $e\in \mathcal{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 $\mathbb{Q}_h$ is an $L^2$ projection. Next, we turn to the first term. By using $\Pi_h A\nabla u\in H(\mathrm{div}, \Omega)$ and (2.4), an integration by parts leads to

It is easy to see that $T = 0$ when $\epsilon = -1$, hence, one has to estimate $T$ if $\epsilon = 0$ or $1$. The same technique as in (3.6) leads to

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_\epsilon$ and (4.14), we arrive at

which completes the proof.

5.   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 $\Omega = (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\times 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 = \frac{1}{N}$.

Example 1 (A model problem with a smooth solution). Consider the problem (1.1)-(1.2) with the following analytical solution

and the diffusion matrix $A = \begin{bmatrix}1 & 0 \\ 0 & 1 \end{bmatrix}$.

For comparison between the IPWG method and the original WG method, the relative errors in $||| \cdot |||$ and $L^2$-norms, defined as

of the WG method with the element $(\mathbb{P}_k, \mathbb{P}_{k}, \mathcal{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 $\epsilon = -1, 0, 1$ converge optimally with rates $O(h^{k})$ in $|||\cdot |||$-norm. Since the diffusion coefficient $A$ is constant and the domain $\Omega$ is convex, the duality problem possesses the $H^2$-regularity. Thus, from Theorem 4.2, the IPWG solution with $\epsilon = -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 $\epsilon = 0, 1$ without over-penalization.

The relative errors between the exact solution and IPWG ($\epsilon = -1, \beta = 1$) solution in $|||\cdot |||$- and $L^2$-norm are listed in Table 2. Note that, for $\epsilon = 0, 1$, the IPWG methods produce the same numerical results as in the case of $\epsilon = -1$, see Table 2. It can be seen from Tables 2 and 3 that errors of the IPWG ($\epsilon = -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 $\epsilon = 1, \ \sigma^e = 0$ also has optimal rates while there is no theoretical proof.

Example 2 (A model problem with a corner singularity). We consider the Poisson equation defined on the unit square with the exact solution

where $r = \sqrt{x^{2}+y^{2}}$ and $\alpha\in (0, 1]$. As we known

where $\epsilon$ is any small, but positive number. In numerical experiments, we set $\alpha = 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 $\beta$ such that $\beta(d-1)-1 = 0$, i.e., $\beta = 1$ for $d = 2$. From Table 4, we observe all cases of the IPWG method with $\epsilon = -1, 0, 1$ converge optimally at $O(h^{\alpha})$ in $|||\cdot |||$-norm and $O(h^{1+\alpha})$ 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 $\epsilon = -1$, $\beta = 1$, and we set $\sigma = 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 $|||\cdot |||$ and $L^2$-norms are observed.

6.   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 $\sigma^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.

Acknowledgments

The research of the second author was supported in part by the Natural Science Foundation of Gansu Province, China (Grant 18JR3RA290).

Conflict of interest

We declare no conflicts of interest in this paper.