A Discontinuous Galerkin Method for the Stokes Equation by Divergence-free Patch Reconstruction

TL;DR

This paper introduces a discontinuous Galerkin method using divergence-free patch reconstruction for Stokes flows, transforming the problem into an elliptic system and providing error estimates validated by numerical tests.

Contribution

The paper presents a novel divergence-free patch reconstruction approach within a DG framework that simplifies the Stokes problem into an elliptic system, reducing computational complexity.

Findings

Achieves the same degrees of freedom as the mesh elements.

Provides classical error estimates for the method.

Numerical examples verify the theoretical error bounds.

Abstract

A discontinuous Galerkin method by patch reconstruction is proposed for Stokes flows. A locally divergence-free reconstruction space is employed as the approximation space, and the interior penalty method is adopted which imposes the normal component penalty terms to cancel out the pressure term. Consequently, the Stokes equation can be solved as an elliptic system instead of a saddle-point problem due to such weak form. The number of degree of freedoms of our method is the same as the number of elements in the mesh for different order of accuracy. The error estimations of the proposed method are given in a classical style, which are then verified by some numerical examples.

Figures

Tables

Table 1. Table 4.1. Numerical errors on the quasi-uniform triangle meshes for Example 1.

		h=1.0E-1	5.0E-2		2.5E-2		1.25E-2		6.25E-3
$P^m$	Norm	error	error	order	error	order	error	order	error	order
	$\parallel \cdot {\parallel }_{L^2}$	5.41E-02	2.41E-02	1.16	9.97E-03	1.27	4.20E-03	1.24	1.92E-03	1.13
	$\parallel \cdot {\parallel }_{1,h}$	1.50E+00	8.76E-01	0.77	4.61E-01	0.92	2.55E-01	0.85	1.48E-01	0.77
1	$\parallel \cdot {\parallel }_h$	2.10E+00	1.39E+00	0.59	8.16E-01	0.77	4.97E-01	0.71	3.15E-01	0.65
	$\parallel \cdot {\parallel }_{L^2}$	1.47E-02	2.88E-03	2.34	6.79E-04	2.08	1.72E-04	1.97	4.47E-05	1.94
	$\parallel \cdot {\parallel }_{1,h}$	3.97E-01	1.19E-01	1.73	3.78E-02	1.66	1.28E-02	1.55	4.84E-03	1.40
2	$\parallel \cdot {\parallel }_h$	7.72E-01	3.52E-01	1.13	1.55E-01	1.17	7.53E-01	1.04	3.71E-03	1.01
	$\parallel \cdot {\parallel }_{L^2}$	5.84E-03	8.18E-04	2.83	1.06E-04	2.94	1.88E-05	2.49	2.97E-06	2.66
	$\parallel \cdot {\parallel }_{1,h}$	1.41E-01	2.67E-02	2.40	5.30E-03	2.33	1.48E-03	1.83	3.54E-04	2.06
3	$\parallel \cdot {\parallel }_h$	3.33E-01	9.47E-01	1.81	2.66E-02	1.83	9.20E-03	1.53	3.04E-03	1.59
	$\parallel \cdot {\parallel }_{L^2}$	1.14E-03	1.03E-04	3.46	8.25E-06	3.64	1.08E-06	2.92	1.36E-07	2.99
	$\parallel \cdot {\parallel }_{1,h}$	3.24E-02	4.06E-03	2.99	4.99E-04	3.02	9.96E-05	2.32	1.77E-05	2.48
4	$\parallel \cdot {\parallel }_h$	9.44E-02	1.62E-02	2.54	3.76E-03	2.10	9.15E-04	2.04	2.12E-04	2.10

Table 2. Table 4.2. Numerical errors on the quasi-uniform quadrilateral meshes for Example 1.

		h=1.0E-1	5.0E-2		2.5E-2		1.25E-2		6.25E-3
$P^m$	Norm	error	error	order	error	order	error	order	error	order
	$\parallel \cdot {\parallel }_{L^2}$	1.03E-01	4.41E-02	1.22	1.74E-02	1.33	8.14E-03	1.10	3.67E-03	1.14
	$\parallel \cdot {\parallel }_{1,h}$	1.97E+00	1.08E+00	0.86	5.54E-01	0.96	2.96E-01	0.90	1.62E-01	0.87
1	$\parallel \cdot {\parallel }_h$	3.10E+00	1.92E+00	0.69	1.22E+00	0.64	8.27E-01	0.56	5.50E-01	0.58
	$\parallel \cdot {\parallel }_{L^2}$	3.80E-02	7.71E-03	2.30	1.92E-03	1.99	4.97E-04	1.95	1.21E-04	2.03
	$\parallel \cdot {\parallel }_{1,h}$	7.95E-01	2.06E-01	1.94	7.11E-02	1.53	2.81E-02	1.33	1.10E-02	1.35
2	$\parallel \cdot {\parallel }_h$	1.48E+00	5.82E-01	1.35	2.63E-01	1.14	1.22E-01	1.10	5.58E-03	1.13
	$\parallel \cdot {\parallel }_{L^2}$	1.11E-02	1.88E-03	2.56	2.75E-04	2.77	3.94E-05	2.80	5.96E-06	2.72
	$\parallel \cdot {\parallel }_{1,h}$	2.59E-01	3.31E-02	2.96	6.17E-03	2.42	1.74E-03	1.82	4.92E-04	1.82
3	$\parallel \cdot {\parallel }_h$	4.11E-01	9.71E-02	2.08	3.06E-02	1.66	1.10E-02	1.46	3.82E-03	1.53
	$\parallel \cdot {\parallel }_{L^2}$	8.18E-03	6.85E-04	3.57	5.32E-05	3.68	4.15E-06	3.68	4.04E-07	3.35
	$\parallel \cdot {\parallel }_{1,h}$	1.59E-01	9.69E-03	4.03	1.26E-03	2.94	2.18E-04	2.52	4.06E-05	2.42
4	$\parallel \cdot {\parallel }_h$	2.52E-01	3.55E-02	2.82	8.73E-03	2.02	2.25E-03	1.95	5.62E-04	2.00

Equations

$$\Vert v{\Vert }_{L^2(\partial K)}^2\le C\left(h_K^{-1}\Vert v{\Vert }_{L^2(K)}^2+h_K\Vert \nabla v{\Vert }_{L^2(K)}^2\right),\text{for all}v\in H^1(K).$$

$$\Vert v-\tilde{v}{\Vert }_{L^2(K)}+h_K\Vert \nabla (v-\tilde{v}){\Vert }_{L^2(K)}\le C{h}_K^{r+1}{\left|v\right|}_{H^{r+1}(K)}.$$

$$\Vert \nabla v{\Vert }_{L^2(K)}\le C{h}_K^{-1}\Vert v{\Vert }_{L^2(K)}.$$

$$\Vert v{\Vert }_{L^2(\partial K)}\le C{h}_K^{-1/2}\Vert v{\Vert }_{L^2(K)}.$$

$$\Vert \mathbf{v}{\Vert }_{H^m(D)}^2\triangleq \sum _{i=1}^2\Vert {\mathbf{v}}_i{\Vert }_{H^m(D)}^2,\mathbf{v}\in {\mathbf{H}}^m(D).$$

$${\mathbf{S}}^m(D)=\{\mathbf{v}\in {\mathbf{H}}^m(D):\nabla \cdot \mathbf{v}=0\text{in}D\},$$

$${\mathbb{S}}_m(D)=\{\mathbf{v}\in {\mathbb{P}}_m(D):\nabla \cdot \mathbf{v}=0\text{in}D\}.$$

$$\Vert \mathbf{v}-\tilde{\mathbf{v}}{\Vert }_{L^2(K)}+h_K\Vert \nabla (\mathbf{v}-\tilde{\mathbf{v}}){\Vert }_{L^2(K)}\le C{h}_K^{r+1}{\left|\mathbf{v}\right|}_{H^{r+1}(K)}.$$

$$\mathcal{I}(K)=\left\{{\mathbf{x}}_{K^{\prime }}|K^{\prime }\in S(K)\right\},$$

$$U_h:=\{v\in L^2(\Omega )\mid v{|}_K\in {\mathbb{P}}_0(K)\}.$$

$${\mathcal{S}}_K\mathbf{v}=\arg \underset{\mathbf{p}\in {\mathbb{S}}_m(S(K))}{\min }\sum _{\mathbf{x}\in \mathcal{I}(K)}{\left|\mathbf{v}(\mathbf{x})-\mathbf{p}(\mathbf{x})\right|}^2.$$

$${\mathcal{S}}_K:\mathbf{v}\mapsto {\mathcal{S}}_K\mathbf{v}={\mathcal{S}}_K\tilde{\mathbf{v}},\forall \mathbf{v}\in [C^0(S(K)){]}^2\cap {\mathbf{S}}^m(S(K)),$$

$$\mathbf{p}{|}_{\mathcal{I}(K)}=0\text{implies}\mathbf{p}{|}_{S(K)}\equiv 0.$$

$$\Lambda (m,\mathcal{I}(K))<\infty$$

$$\Lambda (m,\mathcal{I}(K)):=\underset{\mathbf{p}\in {\mathbb{S}}_m(S(K))}{\max }\frac{\Vert \mathbf{p}{\Vert }_{L^{\infty }(S(K))}}{\Vert \mathbf{p}{|}_{\mathcal{I}(K)}{\Vert }_{{\ell }_{\infty }}}.$$

$${\mathcal{S}}_K\mathbf{g}=\mathbf{g}\text{for all}\mathbf{g}\in {\mathbb{S}}_m(S(K)).$$

$$\Vert {\mathcal{S}}_K\mathbf{g}{\Vert }_{L^{\infty }(K)}\le \Lambda (m,\mathcal{I}(K))\sqrt{\#\mathcal{I}(K)}\Vert \mathbf{g}{|}_{\mathcal{I}(K)}{\Vert }_{{\ell }_{\infty }},$$

$$\Vert \mathbf{g}-{\mathcal{S}}_K\mathbf{g}{\Vert }_{L^{\infty }(K)}\le {\Lambda }_m\underset{\mathbf{p}\in {\mathbb{S}}_m(S(K))}{\inf }\Vert \mathbf{g}-\mathbf{p}{\Vert }_{L^{\infty }(S(K))},$$

$$\Vert {\mathcal{S}}_K\mathbf{g}{\Vert }_{L^{\infty }(K)}\le \Vert {\mathcal{S}}_K\mathbf{g}{\Vert }_{L^{\infty }(S(K))}\le \Lambda (m,\mathcal{I}(K))\underset{x\in {\mathcal{I}}_k}{\max }|{\mathcal{S}}_K\mathbf{g}(x)|.$$

$$\sum _{x\in \mathcal{I}(K)}{\left({\mathcal{S}}_K\mathbf{g}(x)\right)}^2\le \sum _{x\in \mathcal{I}}{\left(\mathbf{g}(x)\right)}^2\le \#\mathcal{I}(K)\underset{x\in {\mathcal{I}}_k}{\max }|\mathbf{g}(x){|}^2.$$

$$\Vert {\mathcal{S}}_K\mathbf{g}{\Vert }_{L^{\infty }(K)}\le \Lambda (m,\mathcal{I}(K))\sqrt{\#\mathcal{I}(K)}\Vert \mathbf{g}{|}_{\mathcal{I}(K)}{\Vert }_{{\ell }_{\infty }}.$$

$$\Vert \mathbf{g}-{\mathbf{p}}_0{\Vert }_{L^{\infty }(S(K))}=\underset{\mathbf{p}\in {\mathbb{S}}_m(S(K))}{\inf }\Vert \mathbf{g}-\mathbf{p}{\Vert }_{L^{\infty }(S(K))}.$$

$$\Vert {\mathcal{S}}_K\mathbf{g}-{\mathbf{p}}_0{\Vert }_{L^{\infty }(K)}$$

$$\Vert \mathbf{g}-{\mathcal{S}}_K\mathbf{g}{\Vert }_{L^{\infty }(K)}$$

$$\Vert \mathbf{g}-{\mathcal{S}}_K\mathbf{g}{\Vert }_{L^2(K)}$$

$$\Vert \nabla (\mathbf{g}-{\mathcal{S}}_K\mathbf{g}){\Vert }_{L^2(K)}$$

$$\Vert \mathbf{g}-{\mathcal{S}}_K\mathbf{g}{\Vert }_{L^2(\partial K)}$$

$$\underset{p\in {\mathbb{S}}_m(S(K))}{\inf }\Vert \mathbf{g}-\mathbf{p}{\Vert }_{L^{\infty }(S(K))}\le \Vert \mathbf{g}-\mathbf{\chi }{\Vert }_{L^{\infty }(S(K))}\le C{d}_K^m{\left|\mathbf{g}\right|}_{H^{m+1}(S(K))},$$

$$\Vert \mathbf{g}-{\mathcal{S}}_K\mathbf{g}{\Vert }_{L^2(K)}\le {\left|K\right|}^{1/2}\Vert \mathbf{g}-{\mathcal{S}}_K\mathbf{g}{\Vert }_{L^{\infty }(K)}\le C{\Lambda }_m{h}_K{d}_K^m{\left|\mathbf{g}\right|}_{H^{m+1}(S(K))}.$$

$$\Vert \nabla (\mathbf{g}-{\mathcal{S}}_K\mathbf{g}){\Vert }_{L^2(K)}$$

$$\le C{h}_K^m{\left|\mathbf{g}\right|}_{H^{m+1}(K)}+C{h}_K^{-1}\Vert {\hat{\mathbf{g}}}_m-{\mathcal{S}}_{K}g{\Vert }_{L^2(K)}$$

$$\le C{h}_K^m{\left|\mathbf{g}\right|}_{H^{m+1}(K)}+C{h}_K^{-1}\Vert \mathbf{g}-{\hat{\mathbf{g}}}_m{\Vert }_{L^2(K)}+C{h}_K^{-1}\Vert \mathbf{g}-{\mathcal{S}}_K\mathbf{g}{\Vert }_{L^2(K)}$$

Full text

Keywords: discontinued Galerkin; patch reconstruction; Stokes equation; divergence-free; interior penalty

MSC2010: 65N30;76D07

1. Introduction

The incompressible Stokes equations describe the low Reynolds number flows which are the linearization of the Navier-Stokes equations [18]. In the finite element method for incompressible Stokes problem see, e.g. [8], the major concern is how to impose the incompressibility condition. The conforming mixed finite element methods are usually introduced weakly by the Lagrange multiplier, and solve a saddle points problem. The methods are required to setup the approximation spaces for velocity and pressure to satisfy the in-sup conditions, which is necessary to guarantee the numerical stability. Consequently, the numerical solution obtained is weakly divergence-free, we refer readers to [4] for more details.

The other approach to solve the incompressible Stokes equation is constructing a divergence-free space to automatically fulfill the incompressible condition and eliminate the pressure. There are many efforts on constructing global divergence-free basis function. For example, the Crouzeix-Raviart elements was proposed in [8], which are divergence-free in each element and continuous at the midpoint of element edges. The $H(div)$ finite elements was proposed in [19]. The methods to construct these divergence-free spaces are quite subtle that it is a nontrivial job to extend these methods to a grid with its elements in unusual geometry. An alternative way is abandoning the normal continuity of basis function, and using the locally divergence-free elements which are proposed in [2]. See the recent developments in this fold in [11, 7, 16].

The discontinuous Galerkin (DG) method by patch reconstruction was recently introduced in [13] to solving the elliptic equation, then was developed in [12] and [15] to various model problems. We are motivated to approximate the divergence-free velocity in Stokes problem using the discontinuous space by patch reconstruction. Actually, we propose a locally divergence-free space by imposing the constrain locally in patch reconstruction. With the new space, we adopt the interior penalty method to Stokes problems [10]. The penalty term is introduced on the normal component to control the inconsistent error. Consequently, an elliptic problem is attained to be solved instead of a saddle-point problem. To construct the locally divergence-free space, we use the piecewisely solenoidal polynomial functions only. The space has one degree of freedom in each element. Therefore, the dimension of approximation space does not depend on the polynomial orders. The interpolation error estimate of this space can be obtained naturally, which leads to the approximate error estimate of the Stokes problem following the standard techniques. We note that the new space is a subspace of canonical locally divergence-free DG space. Our method enjoys the advantages of DG method, for example that it works on meshes with polygonal elements. We demonstrate by numerical examples on different meshes.

The rest of this paper is organized as follows. In Section 2, we describe the method to construct the locally divergence-free approximation space and and give the approximation estimate of this space. Then we present the interior penalty method of discontinued Galerkin method for Stokes problems in Section 3. The error estimate of the proposed method is obtained under the DG energy norm. Finally, numerical results for two-dimensional examples are presented in Section 4 to validate our estimates and demonstrate the capacity of our method.

2. Construction of Approximation Space

Let $\Omega$ be a polygonal domain in $\mathbb{R}^{2}$. The mesh $\mathcal{T}_{h}$ is a partition of $\Omega$ with polygons denoted by $K$ that $\bigcup_{K\in\mathcal{T}_{h}}\bar{K}=\Omega$. For any two elements $K_{0}$, $K_{1}\in\mathcal{T}_{h}$, $K_{0}\bigcap K_{1}=\emptyset$ if $K_{0}\neq K_{1}$. Here $h{:}=\max_{K\in\mathcal{T}_{h}}h_{K}$ with $h_{K}$ the diameter of $K$. We denote by $\left\lvert K\right\rvert$ the area of $K$. $\Gamma_{h}$ denotes the union of boundaries of element $K\in\mathcal{T}_{h}$, where $\Gamma_{h}^{0}$ is the union of interior boundaries $\Gamma_{h}^{0}\triangleq\Gamma_{h}\setminus\partial\Omega$. We assume that the partition $\mathcal{T}_{h}$ possesses the following shape regularity conditions  [5, 3]:

1. A1 

There exists an integer $N$ independent of $h$, that any element $K$ admits a sub-decomposition $\widetilde{\mathcal{T}_{h}}|_{K}$ which consists of at most $N$ triangles $T$. 2. A2 

$\widetilde{\mathcal{T}_{h}}$ is a compatible sub-decomposition, that any $T\in\widetilde{\mathcal{T}_{h}}$ is shape-regular in the sense of Ciarlet-Raviart [6]: there exists a real positive number $\sigma$ independent of $h$ such that $h_{T}/\rho_{T}\leq\sigma$, where $\rho_{T}$ is the radius of the greatest ball inscribed in $T$.

Let $D$ to be a subdomain of $\Omega$, which may be an element or an aggregation of the elements belong to $\mathcal{T}_{h}$. Let $H^{m}(D)$ demote the Sobolev space of real valued functions on $D$ with a positive integer $m$. For vector valued functions, we introduce the spaces $\boldsymbol{H}^{m}(D)=[H^{m}(D)]^{2}$. Let $\mathbb{P}_{m}(D)$ be a set of polynomial with total degree not greater than $m$ on domain $D$, and the corresponding vector valued polynomial spaces $[\mathbb{P}_{m}(D)]^{2}$ is denoted by $\boldsymbol{\mathbb{P}}_{m}(D)$.

Assumptions A1 and A2 allow quite general shapes in the partition, such as non-convex or degenerate elements. They also lead to some common used properties and inequalities:

1. M1

$\forall T\in\widetilde{\mathcal{T}_{h}}$, there exists $\rho_{1}\geq 1$ that depends on $N$ and $\sigma$ such that $h_{K}/h_{T}\leq\rho_{1}$. 2. M2

[Agmon inequality] There exists a constant $C$ that depends on $N$ and $\sigma$, but independent of $h_{K}$ such that

[TABLE] 3. M3

[Approximation property] There exists a constant $C$ that depends on $N,r$ and $\sigma$, but independent of $h_{K}$ such that for any $v\in H^{r+1}(K)$, there exists an approximation polynomial $\widetilde{v}\in\mathbb{P}_{r}(K)$ such that

[TABLE] 4. M4

[Inverse inequality] For any $v\in\mathbb{P}_{r}(K)$, there exists a constant $C$ that depends only on $N,r,\sigma$ and $\rho_{1}$ such that

[TABLE] 5. M5

[Discrete trace inequality] For any $v\in\mathbb{P}_{r}(K)$, there exists a constant $C$ that depends only on $N,r,\sigma$ and $\rho_{1}$ such that

[TABLE]

The above four inequalities (2.1), (2.2), (2.3) and (2.4) are hold for vector valued space $\boldsymbol{H}^{m}(D)$ and $\boldsymbol{\mathbb{P}}_{m}(D)$ with corresponding norms,

[TABLE]

We are interested in the spaces of solenoidal vector fields

[TABLE]

and the polynomial solenoidal vector field

[TABLE]

The above inequalities are also hold for the $\boldsymbol{S}^{m}(D)$ and $\mathbb{S}_{m}(D)$, see in [2]. Here we restate the approximation property:

M3 [Approximation property] There exists a constant $C$ that depends on $N,r$ and $\sigma$, but independent of $h_{K}$ such that for any $\boldsymbol{v}\in\boldsymbol{S}^{r+1}(K)$, there exists an approximation polynomial $\widetilde{\boldsymbol{v}}\in\mathbb{S}_{r}(K)$ such that

[TABLE]

In [13], we introduced the reconstruction operator which mapping the piecewise constant space to discontinuous piecewise polynomial space. In this paper, we intend to construct a reconstruction operator $\mathcal{S}$ that embeds the piecewise constant vector valued space to discontinuous piecewise solenoidal vector fields. For each element $K\in\mathcal{T}_{h}$, we assign a sampling node or collocation point $\boldsymbol{x}_{K}\in K$ and element patch $S(K)$. $S(K)$ usually contains $K$ and some elements around $K$. It is quite flexible to assign the sampling nodes and construct the element patch. Usually, we let the barycenter of the element $K$ to be the sampling node, while a perturbation is allowed, cf.  [13]. The element patch is built up by adding Von Neumann neighbors (adjacent edge-neighboring elements) recursively until the size of the element patch reaches a fixed number, for the details we refer to [13, 14].

Let $\mathcal{I}(K)$ denote the set of the sampling nodes belonging to $S(K)$ with $\#\mathcal{I}(K)$ denote its cardinality

[TABLE]

and let $\#S(K)$ be the number of elements belonging to $S(K)$. These two numbers $\#\mathcal{I}(K)$ and $\#S(K)$ are equal to each other. We define $d_{K}{:}=\text{diam}\;S(K)$ and $d=\max_{K\in\mathcal{T}_{h}}d_{K}$.

Denote the piecewise constant space associated with $\mathcal{T}_{h}$ by $U_{h}$, i.e.,

[TABLE]

and let $\boldsymbol{U}_{h}=[U_{h}]^{2}$ to be the piecewise constant vector valued space.

For any $\boldsymbol{v}\in\boldsymbol{U}_{h}$ and for any $K\in\mathcal{T}_{h}$, we reconstruct a high order solenoidal polynomials $\mathcal{S}_{K}\boldsymbol{v}$ of degree $m$ by solving the following discrete least-square problem.

[TABLE]

Though $\mathcal{S}_{K}\boldsymbol{v}$ gives a solenoidal approximation on the entire element patch $S(K)$, but we only use the approximation on the element $K$. Then the reconstruction operator $\mathcal{S}_{K}$ can extended to the function space $[C^{0}(S(K))]^{2}\cap\boldsymbol{S}^{m}(S(K))$, still denote by $\mathcal{S}_{K}$ without ambiguity,

[TABLE]

where $\tilde{\boldsymbol{v}}(\boldsymbol{x}_{K}^{\prime})=\boldsymbol{v}(\boldsymbol{x}_{K}^{\prime}),\forall\boldsymbol{x}_{K}^{\prime}\in\mathcal{I}(K)$.

We make the following assumption on the sampling node set $\mathcal{I}(K)$.

Assumption A  For any $K\in\mathcal{T}_{h}$ and $\boldsymbol{p}\in\mathbb{S}_{m}(S(K))$,

[TABLE]

This assumption can guarantee the uniqueness of the solution of discrete least-square problem (2.6). Obvious, a necessary condition for the solvability of  (2.6) is that $\#\mathcal{I}(K)$ is greater than $(m+1)(m+4)/4$. Assumption A is equal to the following quantitative estimate

[TABLE]

with

[TABLE]

If the mesh is quasi-uniform triangulation and each element patch is convex, the quantitative estimate of the uniform upper bound of $\Lambda(m,\mathcal{I}_{K})$ for the general polynomial is obtained in [14]. The requirements of the uniform upper bound are hard to be satisfied while the polygonal meshes are applied. We refer the readers to [13] which gives the uniform upper bound under milder assumptions, such as polygonal partition and non-convex element patch. Due to the solenoidal polynomial space is the subspace of the general polynomial vector valued space, $\Lambda(m,\mathcal{I}(K))$ in (2.8) is uniformly bounded under the same assumption.

Lemma 2.1**.**

If Assumption A holds, then there exists a unique solution of (2.6), denoted by $\mathcal{S}_{K}\boldsymbol{v}$, for any $K\in\mathcal{T}_{h}$. Moreover $\mathcal{S}_{K}$ satisfies

[TABLE]

The stability property holds true for any $K\in\mathcal{T}_{h}$ and $\boldsymbol{g}\in[C^{0}(S(K))]^{2}$ and $\nabla\cdot\boldsymbol{g}=0$ as

[TABLE]

and the quasi-optimal approximation property is valid in the sense that

[TABLE]

where $\Lambda_{m}{:}=\max_{K\in\mathcal{T}_{h}}\{1+\Lambda(m,\mathcal{I}(K))\sqrt{\#\mathcal{I}(K)}\}$.

Proof.

The reconstruction operator $\mathcal{S}_{K}$ is a projection operator from $[C^{0}(\Omega)]^{2}$ to $\mathbb{S}_{m}(S(K))$ with the discrete $l_{2}$ norm, the identity (2.9) is obvious.

By the Assumption A and definition of $\Lambda(m,\mathcal{I}(k))$ in equation (2.8),

[TABLE]

From the projection property of operator $\mathcal{S}_{K}$, we have

[TABLE]

Combining  (2.12) and (2.13), we get

[TABLE]

Assume that $\boldsymbol{p}_{0}$ is the best approximation of $\boldsymbol{g}$ under $L^{\infty}$ norm, $\boldsymbol{p}_{0}\in\mathbb{S}_{m}(S(K))$, and

[TABLE]

Then, we have following estimate:

[TABLE]

By triangle inequality, the left side of (2.11) can be written as

[TABLE]

Together with the definition of $\Lambda_{m}$, it implies the quasi-optimality (2.11). ∎

Lemma 2.2**.**

If Assumption A holds and $\boldsymbol{g}\in[C^{0}(S(K))]^{2}\cap\boldsymbol{S}^{m+1}(S(K))$, then there exists a constant $C$ that depends on $N,\sigma$, $\gamma$ and $\rho_{1}$ such that

[TABLE]

Proof.

By [2, Theorem 4.3], we take $\boldsymbol{p}=\boldsymbol{\chi}\in\mathbb{S}_{m}$in the right-hand side of (2.11), where $\boldsymbol{\chi}$ is the approximation solenoidal polynomial of order $m$. Then

[TABLE]

where $C$ depends on $N,m,\sigma$ and $\gamma$.

Substituting the above estimate (2.17) into (2.11), we obtain

[TABLE]

which gives (2.14).

Then, assume that $\widehat{\boldsymbol{g}}_{m}$ be the approximation polynomial in (2.2) for function $\boldsymbol{g}$, by the inverse inequality (2.3) and the approximation estimate (2.14), then we have

[TABLE]

This gives (2.15).

Combing the Agmon inequality (2.1),  (2.14) and  (2.15), one has

[TABLE]

Taking the square root of both sides gives (2.16) completes the proof. ∎

A global reconstruction operator $\mathcal{S}$ is defined by $\mathcal{S}|_{K}=\mathcal{S}_{K}$. Given $\mathcal{S}$, we embed $\boldsymbol{U}_{h}$ into a piecewise discontinuous solenoidal polynomial finite element space with its order to be $m$. The approximation space $\boldsymbol{V}_{h}$ is defined by

[TABLE]

Furthermore, the reconstruction operator $\mathcal{S}$ can be extended to function space $[C^{0}(\Omega)]^{2}\cap\boldsymbol{S}^{m}(\Omega)$, and we still denote by $\mathcal{S}$ without ambiguity,

[TABLE]

where $\tilde{\boldsymbol{u}}\in\boldsymbol{U}_{h}$ and $\tilde{\boldsymbol{u}}(\boldsymbol{x}_{K})=\boldsymbol{u}(\boldsymbol{x}_{K})$.

The basis function of $\boldsymbol{V}_{h}$ are given by the following process. Define $e_{K}\in U_{h}$ to be the characteristic function corresponding to $K$,

[TABLE]

Let $\boldsymbol{\lambda}_{K}$ denote the basis function and it is defined by the reconstruction process:

[TABLE]

The reconstruction operator can be wrote explicitly with the given basis functions $\{\boldsymbol{\lambda}_{K}|\forall K\in\mathcal{T}_{h}\}$,

[TABLE]

Next, we will show the implementation of the proposed method by an example on square domain $[0,1]\times[0,1]$. Here a third order solenoidal field reconstruction is considered, the basis functions of the corresponding solenoidal field $\boldsymbol{\xi}_{j}$ are listed as follows,

[TABLE]

[TABLE]

[TABLE]

The total degree of freedoms of the third order solenoidal field is 14, it means the size of element patch at least need to be taken as 7 to guarantee the solvability of the least square problem. Thus we take $\#S(K)$ as 10 and the barycenter as the sampling node. Figure 2.1 shows domain and the corresponding triangulation,

Here we take an element $K_{0}$ as a demonstration element, the element patch of $K_{0}$ is taken as

[TABLE]

and the corresponding set of the sampling nodes is

[TABLE]

The element patch and sampling nodes are shown in Figure 2.1.

For a given function $\boldsymbol{g}=[g_{1},g_{2}]^{T}\in\boldsymbol{S}^{m}(\Omega)$, the least square problem (2.6) is as

[TABLE]

The problem is solved directly by calculating the generalized inverse of a matrix,

[TABLE]

where $A$ and $b$ are defined as follows

[TABLE]

$A$ is a $20\times 14$ matrix, limited by the page space, we only list the first order part, the rest is easy to be complemented. $b$ is a $20\times 1$ vector.

Matrix $(A^{T}A)^{-1}A^{T}$ contains all the information of the basis functions that are defined on element $K_{0}$, the matrix is relevant with $\boldsymbol{\lambda}_{K_{i}},i=0,\cdots,9$. The basis function $\boldsymbol{\lambda}_{K_{0}}$ is determined after the reconstruction are conducted on each element. Figure 2.2 shows the basis function $\boldsymbol{\lambda}_{K_{0}}$. We note that the support of $\boldsymbol{\lambda}_{K_{0}}$ is not equal to the element patch $S(K_{0})$. Insteadly, for any element $K\in\mathcal{T}_{h}$, the support of the basis function $\boldsymbol{\lambda}_{K}$ is related with all the element patches which includes $K$:

[TABLE]

Due to the discontinuity of the reconstructed finite element space, the DG method can be directly employed. We end this section by introducing the average and jump operators which are common used in DG method. Let $e$ be an interior edge shared by two adjacent elements $e=\partial K^{+}\cap\partial K^{-}$. The corresponding unit outward normal vector are denoted by $\boldsymbol{\mathrm{n}}^{+}$, $\boldsymbol{\mathrm{n}}^{-}$. Let $v$ and $\boldsymbol{v}$ be the scalar valued and vector valued functions on $\mathcal{T}_{h}$, respectively. Following the traditional DG notations, we define the $average$ operator $\{\cdot\}$ on element edges as follows:

[TABLE]

with $v^{+}=v|_{K^{+}},\ v^{-}=v|_{K^{-}},\ \boldsymbol{v}^{+}=\boldsymbol{v}|_{K^{+}},\ \boldsymbol{v}^{-}=\boldsymbol{v}|_{K^{-}}$. The $jump$ operator $[[\cdot]]$ on element edges is as

[TABLE]

For $e\in\partial\Omega$, we set

[TABLE]

3. Approximation of Stokes Equation

We consider the Stokes equation with Dirichlet boundary condition:

[TABLE]

where $\boldsymbol{u}$ is the velocity field, $p$ is the pressure, and $\boldsymbol{g}$ is the given boundary value satisfying the compatibility condition $(\boldsymbol{g}\cdot\boldsymbol{n},1)_{\partial\Omega}=0$. We take the local divergence-free reconstruction space $\boldsymbol{V}_{h}$ as the trial and test function spaces. The approximation problem of (3.1) is: Seek $\boldsymbol{u}_{h}\in\boldsymbol{V}_{h}$, such that

[TABLE]

where the bilinear form is defined by

[TABLE]

and

[TABLE]

*where $\eta$ and $\epsilon$ are penalty parameters, corresponding to different penalty term, respectively. *

We note that the weak form (3.2) is inconsistent [10]. Precisely, the solution $\boldsymbol{u}$ of the equation (3.1) does not satisfy the weak form (3.2), saying

[TABLE]

The inconsistency is caused by the discontinuous of normal component $\boldsymbol{v}_{h}\cdot\boldsymbol{n}$ on interior element edges. Instead of satisfying the weak form (3.2), $\boldsymbol{u}$ satisfies

[TABLE]

Therefore, for the interior penalty scheme, the non-consistent penalty parameter $\epsilon$ need to be taken great enough to control the consistency error which is similar to the interior penalty method for elliptic equation [1]. Meanwhile, the penalty parameter $\eta$ is chosen to guarantee the stability of the operator $B_{h}(\cdot,\cdot)$. The magnitude of penalty parameters are taken as follows,

[TABLE]

where $m$ is the degree of polynomials in $\boldsymbol{V}_{h}$.

Let us define some mesh depended semi-norms, for $\forall\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}$,

[TABLE]

and introduce the energy norm,

[TABLE]

By the results in the previous section, we instantly have following approximation estimate about energy norm,

Lemma 3.1**.**

For $\boldsymbol{u}\in[H^{m+1}(\Omega)]^{2}$ to be the solution of equation (3.1), and $\boldsymbol{u}_{I}\in\boldsymbol{V}_{h}$ to be the interpolation of $\boldsymbol{u}$, there exists a constant $C$ that depends on $N$, $\sigma$, $\gamma$ and $m$, such that

[TABLE]

Here we need to clarify that the approximation error is not optimal, which is due to the fact that the BDM space is not a subspace of the approximation space $\boldsymbol{V}_{h}$. The third term of energy norm is dominant in the approximation error. This fact makes the numerical solution not able to attain the optimal accuracy order.

For the consistency error, we have following estimate,

Lemma 3.2**.**

For $\forall\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}$ and $p\in H^{1}(\Omega)$, we have

[TABLE]

Proof.

By Agmon inequality (2.1), we have

[TABLE]

∎

Next, the boundedness and the stability of the operator $B_{h}(\cdot,\cdot)$ can be claimed as

Lemma 3.3**.**

For $\forall\boldsymbol{u}_{h}\in\boldsymbol{V}_{h}$ , $\forall\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}$, and sufficiently large $\eta$ and $\epsilon$ we have

[TABLE]

Now, with the above lemmas, we are ready to give the error estimate of the numerical solution (3.2),

Theorem 3.1**.**

For $\boldsymbol{u}\in[H^{m+1}(\Omega)]^{2}$, $p\in H^{1}(\Omega)$ to be the solution of equation (3.1), and $\boldsymbol{u}_{h}\in\boldsymbol{V}_{h}$ to be the solution of equation (3.2), we then have

[TABLE]

Furthermore, the $L^{2}$ error estimate is as

[TABLE]

Proof.

We split the error with an interpolation function $\boldsymbol{u}_{I}$ in $\boldsymbol{V}_{h}$,

[TABLE]

Then together with the consistency error, we have

[TABLE]

From the stability of bilinear operator $B_{h}(\cdot,\cdot)$, the second term in (3.13) has

[TABLE]

Collecting the above estimates and Lemma 3.2, it gives

[TABLE]

The first term in (3.13) is the interpolation error,

[TABLE]

Collecting estimates (3.13) ,(3.16) and (3.17) together, we can have (3.11).

For the $L^{2}$ error estimate, we define the auxiliary functions $(\boldsymbol{\varphi},q)$ and the adjoint problem

[TABLE]

then apply the test function $\boldsymbol{u}-\boldsymbol{u}_{h}$, we have

[TABLE]

The regularity of Stokes equations implies

[TABLE]

where $\boldsymbol{\varphi}_{BDM}$ is the BDM interpolation function. This gives

[TABLE]

and

[TABLE]

Inserting (3.21) and (3.22) to the equation (3.19),

[TABLE]

Then substituting the energy norm estimate (3.11) into (3.23), the $L^{2}$ error estimate (3.12) is obtained. ∎

4. Numerical Examples

We present some numerical examples to illustrate the effectiveness of the proposed method. One merit of the method is the linear system is maintained with the given mesh while ignore the increase of approximation order. The maximum degree of approximation polynomial is taken as $4$ which is limited by the condition number of linear system. The numerical examples with various polygonal meshes are demonstrated. A direct solver is used to solve the linear algebra systems after discretisation.

4.1. Analytical example

This is an example with a smooth solution to verify the convergence order as the error estimate predicted. The Stokes problem with Dirichlet boundary condition is solved in two dimensional square domain $\Omega=[0,1]\times[0,1]$. The exact velocity field and pressure are given as follows,

[TABLE]

The body force $\boldsymbol{f}$ and boundary condition $\boldsymbol{g}$ are given correspondingly. The quasi-uniform triangular and quadrilateral meshes are used, both of them are generated by the software gmsh[9].

Table 4.1 and 4.2 present the error of the numerical solution in $L^{2}$ norm, $|\cdot|_{1,h}$ seminorm and $\|\cdot\|_{h}$ norm. Figure 4.1 and 4.2 plot the accuracy order under different norm, respectively. The behavior of the proposed method agrees perfectly with the theoretical analysis on the convergence order to be $h^{\frac{m}{2}}$ on $\|\cdot\|_{h}$ norm. And the convergence order of $L^{2}$ norm and $|\cdot|_{1,h}$ seminorm are not optimal as usual situations, which are $m+1$ and $m$, respectively. This is due to the fact that the space $\boldsymbol{V}_{h}$ is not large enough that it can not contain the $BDM$ space as a subspace as the traditional discontinued Galerkin space. The difference between $|\cdot|_{1,h}$ seminorm and $\|\cdot\|_{h}$ norm error implies that the seminorm $|\cdot|_{\diamond}$ is dominate in energy norm. The condition number of resulting matrix is mainly determined by the penalty parameter $\epsilon$. And the exponential growth of $\epsilon$ with the increasing of the degree of the polynomial results in very ill-posed linear system.

4.2. Lid-driven cavity example

The lid-driven cavity is a benchmark test for the incompressible flow that does not have an exact solution. Here we sets the cavity domain as a rectangle $\Omega=[0,1]\times[0,1.5]$. The body force is zero, and a Dirichlet boundary condition is given as

[TABLE]

Figure 4.3 shows the velocity fields and the scaled velocity. The domain is segmented by a quasi-uniform triangle mesh meanwhile the third order polynomial is employed. The numerical results fit the realities. The scaled velocity present the expected phenomenon which the Contra vortices appear in the bottom corner. The subtle vortices will capture with refined mesh.

4.3. Flow around cylinder

Flow around cylinder is another representative benchmark test for the incompressible flow. The numerical setting is slightly different from the traditional setting. The polygonal elements are considered to demonstrate the compatibility of the proposed method. While the computational domain is a rectangle minus a circle $\Omega=[0,1.5]\times[0,1]/(x-0.5)^{2}+(y-0.5)^{2}\leq 0.2^{2}$, which is partitioned into polygonal elements mesh by PolyMesher [17]. The Dirichlet boundary condition is given while with zero body force.

[TABLE]

Figure 4.4 shows the polygonal mesh and the corresponding velocity field. The proposed method provides a smooth solution which is agreed with the expected behavior. and it can handle a wide variety of polygonal elements without additional techniques.

5. Conclusions

The discontinuous Galerkin method by locally divergence-free patch reconstruction for simulation of incompressible Stokes problems is developed. The interior penalty method with solenoidal velocity field allows calculating the velocity field with no presence of pressure. This method can achieve high order accuracy with one degree of freedom per element. In spite of such advantage, our method suffers from the inconsistence penalty term $h^{-(m+1)}$ which leads to an ill-posed linear system for large $m$.