A Simplified Weak Galerkin Finite Element Method: Algorithm and Error Estimates
In this article a simplified weak Galerkin finite element method is developed for the Dirichlet boundary value problem of convection-diffusion-reaction 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.
onvection-diffusion-reaction equations, simplified weak Galerkin, finite element methods, error estimates.
Primary, 65N30, 65N15; Secondary, 35J50
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
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 well-posedness of the problem (1)-(2), we assume , , and
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 parameter-independent or free of locking . For the convection-diffusion-reaction equation (1)-(2), the recent work in the context of weak Galerkin includes the algorithm developed and analyzed in , the one in  for singularly perturbed problems, and an earlier one in . 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 Navier-Stokes equations, the biharmonic and elasticity equations, div-curl systems and the Maxwell’s equations, etc. The latest development of the WG methods is the prime-dual 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  for second order elliptic equations in nondivergence form,  for the Fokker-Planck equation, and  for elliptic Cauchy problems.
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  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 convection-diffusion-reaction 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 convection-diffusion-reaction 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 well-posedness 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 inner-product 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 inner-product notation.
2 Algorithm on Polymesh
Assume that the domain is of polygonal type and is partitioned into non-overlap 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:
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:
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
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 least-squares fitting.
On each element , we introduce the following bilinear forms:
For simplicity, we set
for . We further introduce the stabilizer
where is the projection operator onto ; namely is the average of on each edge. In particular, is well-defined 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 user-friendly 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:
where the block components in (14) are given by:
, with ,
, with ,
, with ,
, with ,
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 sub-matrices 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  for a detailed derivation. Specifically, let be the center of T (or any point on the plane), the extension can be represented as follows:
From and (15), we have
Let be the basis function corresponding to the edge of :
Then the coefficient for is given by
It follows that
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:
Thus, with , we have from the weak gradient formula (5) that
For simplicity, we introduce the following functions:
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
4 Stability and Well-Posedness
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 well-posedness of the numerical scheme (13).
Let be a shape-regular polygonal partition of the domain . There exists a constant such that
Moreover, the following Poincaré-type estimate holds true:
From the formula (6) for the weak gradient, we have for any constant vector that
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
which verifies the estimate (22).
To derive the inequality (23), we note the following discrete Poincaré inequality:
On each element , the following identity holds true:
where is the average of on the element .
From the formula (5), we have
Substituting the above identity into (26) yields
which leads to the identify (25).
In the finite element space , we introduce the following semi-norm:
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
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 .
provided that the meshsize of is sufficiently small.
Recall that for any we have
The boundedness estimate (30) is then straightforward from the usual Cauchy-Schwarz and the inequality (23). We shall focus on the derivation of the coercivity inequality (31) in the rest of the proof.