A Simplified Weak Galerkin Finite Element Method: Algorithm and Error Estimates
Abstract
In this article a simplified weak Galerkin finite element method is developed for the Dirichlet boundary value problem of convectiondiffusionreaction equations. The simplified weak Galerkin method utilizes only the degrees of freedom on the boundary of each element and, hence, has significantly reduced computational complexity over the regular weak Galerkin finite element method. A stability and some optimal order error estimates in the and norms are established for the corresponding numerical solutions. Numerical results are presented to verify the theory error estimates and a superconvergence phenomena on rectangular partitions.
onvectiondiffusionreaction equations, simplified weak Galerkin, finite element methods, error estimates.
Primary, 65N30, 65N15; Secondary, 35J50
1 Introduction
This paper is concerned with the development of a simplified formulation for the weak Galerkin finite element method for second order elliptic equations. For simplicity, consider the model problem that seeks an unknown function satisfying
(1)  
(2) 
where is a bounded polytopal domain in with boundary , is the diffusion coefficient, is the convection, and is the reaction coefficient in relevant applications. We assume that is sufficient smooth, , and is piecewise smooth with respect to a partition of the domain. For wellposedness of the problem (1)(2), we assume , , and
(3) 
for a constant .
The model problem (1)(2) arises from many scientific applications such as fluid flow in porous media. Mostly importantly, this model problem has served, and still serves, the scientific computing community as a testbed in the search and design of new and efficient computational algorithms for partial differential equations. The classical Galerkin finite element method (see, e.g., [10, 28, 16]) is particularly a numerical technique originated from the study of elliptic problems closed related to (1)(2) or its variations. In the last three decades, various finite element methods using discontinuous trial and test functions, including discontinuous Galerkin (DG) methods and weak Galerkin (WG) methods, have been developed for numerical solutions of partial differential equations. These developments were often tested over testbed problems such as (1)(2) before they were generalized or applied to more complex problems in science and engineering. The DG method, also known as the interior penalty method in different contexts, was originated in early 70s of the last century for a numerical study of model problems such as (1)(2); see, e.g., [3, 14, 25, 38] for early incubations and [1, 13, 17, 27] for a detailed discussion and recent developments.
The weak Galerkin finite element method is a recently developed discretization framework for partial differential equations [36, 37, 24, 34]. With new concepts referred to as weak differential operators (e.g., weak gradient, weak curl, weak Laplacian etc.) and weak continuity through the use of various stabilizers, the method allows the use of totally discontinuous functions and provides stable numerical schemes that are parameterindependent or free of locking [33]. For the convectiondiffusionreaction equation (1)(2), the recent work in the context of weak Galerkin includes the algorithm developed and analyzed in [9], the one in [20] for singularly perturbed problems, and an earlier one in [39]. The WG finite element method has been rapidly developed and applied to several different types of problems, including second order elliptic problems, the Stokes and NavierStokes equations, the biharmonic and elasticity equations, divcurl systems and the Maxwell’s equations, etc. The latest development of the WG methods is the primedual formulation for problems that are either nonsymmetric or do not have variational forms friendly for numerical use. Details on the new developments can be found in [30] for second order elliptic equations in nondivergence form, [31] for the FokkerPlanck equation, and [32] for elliptic Cauchy problems.
The typical WG method for the model problem (1)(2) seeks weak finite element approximations satisfying and
(4) 
for all test functions satisfying , where is an interpolation of the Dirichlet boundary data, is the discrete weak gradient operator, and is a properly selected stabilizer that gives weak continuities for the numerical solutions. The numerical solution consists of two components: the approximation on each element and the approximation on the boundary of each element. To reduce the computational complexity, some hybridized formulations have been introduced in [22, 29] for the method when applied to the diffusion equation and the biharmonic equation through the elimination of the degrees of freedom associated with the unknown function locally on each element. In the superconvergence study for WG [18] on rectangular elements, this hybridized formulation was further simplified in the description of the numerical algorithm, yielding a simplified weak Galerkin (SWG) finite element scheme for the diffusion equation. In our further investigation of the SWG to the convectiondiffusionreaction equation (1), we came to the conclusion that SWG represents a new discretization scheme that is different from the usual WG through a simple elimination of the unknown . As a result, we believe that a systematic study of the SWG for the convectiondiffusionreaction problem (1)(2) should be conducted for its stability and convergence. This paper is in response to this observation and shall provide a mathematical theory for the stability and the convergence of the simplified weak Galerkin finite element method for the model problem (1)(2). We believe that the result of this paper can be extended to other types of modeling equations.
The paper is organized as follows: In Section 2, we shall describe the simplified weak Galerkin finite element method for (1)(2) on general polygonal partitions. In Section 3, we shall present a computational formula for the element stiffness matrices and the element load vectors from SWG. In Section 4, we provide a mathematical theory for the stability and wellposedness of the SWG scheme. Sections 5 and 6 are devoted to a discussion of the error estimates in a discrete and the norm for the numerical solutions. Finally, in Section 7, we present some numerical results to demonstrate the efficiency and accuracy of the SWG method.
Throughout the rest of the paper, we assume and shall use the standard notations for Sobolev spaces and norms [10, 16]. For any open set , and denote the norm and innerproduct in the Sobolev space consisting of square integrable partial derivatives up to order . When or , we shall drop the corresponding subscripts in the norm and innerproduct notation.
2 Algorithm on Polymesh
Assume that the domain is of polygonal type and is partitioned into nonoverlap polygons that are shape regular. For each , denote by its diameter and by the number of edges. For each edge , denote by the midpoints and the outward normal direction of (see Fig. 1 for an illuatration). The meshsize of is defined as .
Let be a piecewise constant function defined on the boundary of , i.e.,
with being a constant. We define the weak gradient of on by:
(5) 
where is the length of the edge and is the area of the element . It is not hard to see that the weak gradient satisfies the following equation:
(6) 
for all constant vector . Here and in what follows of the paper, stands for the usual inner product in .
Denote by the space of piecewise constant functions on . The global finite element space is constructed by patching together all the local elements through single values on interior edges. The subspace of consisting of functions with vanishing boundary value is denoted as .
We use the conventional notation of for the space of polynomials of degree on . For each , we associate it with a linear extension in , denoted as , satisfying
(7) 
It is easy to see that is well defined by (7), and its computation is local and straightforward. In fact, can be viewed as an extension of from to through a leastsquares fitting.
On each element , we introduce the following bilinear forms:
(8)  
(9)  
(10) 
For simplicity, we set
(11) 
for . We further introduce the stabilizer
(12) 
where is the projection operator onto ; namely is the average of on each edge. In particular, is welldefined and takes the average of the Dirichlet data on each boundary edge.
3 Element Stiffness Matrices
The simplified weak Galerkin finite element method (13) is userfriendly in computer implementation. In this section, we present a formula for the computation of the element stiffness matrices and the element load vector on general polygonal elements.
Let be a polygonal element of sides. Denote by the vector representation of given by . Then, the element stiffness matrix and the element load vector for the SWG scheme (13) are given in a block matrix form as follows:
(14) 
where the block components in (14) are given by:

,

, with ,

, with ,

, with ,

, with ,

and ,

and are given by
Here is any point on the plane (e.g., the center of as a specific case), is the midpoint of , is the length of edge , is the unit outward normal vector on , and is the area of the element .
From (13), the element stiffness matrix on consists of two submatrices corresponding to the following forms:
The bilinear form is composed of three bilinear forms given by (11). The rest of this section is devoted to a computation of the element stiffness matrices for each of the bilinear forms involved.
3.1 The stiffness matrix for
For the element stiffness matrix corresponding to , the key is to compute and which can be accomplished through its definition (7); readers are referred to [21] for a detailed derivation. Specifically, let be the center of T (or any point on the plane), the extension can be represented as follows:
where
(15) 
From and (15), we have
(16) 
Let be the basis function corresponding to the edge of :
Then the coefficient for is given by
It follows that
(17) 
where is the identity matrix of size .
3.2 The stiffness matrix for
For a computation of the element stiffness matrix corresponding to the bilinear form , we have from the weak gradient formula (5) that
which leads to the block matrix in the element stiffness matrix.
3.3 The stiffness matrix for
Recall that the bilinear form is given by
Note that the extension has the following representation:
where
(18) 
Thus, with , we have from the weak gradient formula (5) that
(19) 
For simplicity, we introduce the following functions:
(20) 
Then, the equation (19) indicates that the element stiffness matrix corresponding to the bilinear form is given by
3.4 The stiffness matrix for
Recall that the bilinear form is given by
Thus, the element stiffness matrix corresponding to has the following formula:
where is the function defined in (20).
3.5 The element load vector
Finally, the element load vector can be obtained from
for .
4 Stability and WellPosedness
The SWG scheme (13) can be derived from the classical weak Galerkin finite element method [36, 24, 37] by eliminating the degrees of freedom associated with the interior of each element when and . But for the general case of and , the SWG finite element method (13) is different from the weak Galerkin schemes in existing literature. It is thus necessary to provide a mathematical theory for the stability and wellposedness of the numerical scheme (13).
Let be a shaperegular polygonal partition of the domain . There exists a constant such that
(21)  
(22) 
Moreover, the following Poincarétype estimate holds true:
(23) 
From the formula (6) for the weak gradient, we have for any constant vector that
which gives
Hence, by letting we arrive at
which verifies (21).
Next, from the usual error estimate for the projection operator and the estimate (21), we have
It follows that
(24) 
which verifies the estimate (22).
To derive the inequality (23), we note the following discrete Poincaré inequality:
Combining the above estimate with (21) and (21) gives rise to the desired inequality (23). This completes the proof of the lemma.
On each element , the following identity holds true:
(25) 
where is the average of on the element .
From the formula (5), we have
(26) 
Note that
Substituting the above identity into (26) yields
(27) 
which leads to the identify (25).
In the finite element space , we introduce the following seminorm:
(28) 
We claim that defines a norm in the closed subspace . It suffices to show that for any satisfying . In fact, if , then from (28) we have
It follows that on each element
(29) 
for . Thus,
so that has constant value on each element . By using (29) we see that on each edge, which, together with the fact that on , leads to in .
For the model problem (1), assume that and the condition (3) is satisfied. Then, the bilinear form is bounded and coercive in the finite element space ; i.e., there exist constants and such that
(30)  
(31) 
provided that the meshsize of is sufficiently small.