# A Weak Galerkin Mixed Finite Element Method for second order elliptic
equations on 2D Curved Domains

Yi Liu  School of Mathematical Sciences and Jiangsu Key Laboratory for NSLSCS,
Nanjing Normal University, Nanjing, China<EMAIL_ADDRESS>,  Wenbin
Chen  School of Mathematical Sciences and Shanghai Key Laboratory for
Contemporary Applied Mathematics, Fudan University, Shanghai, China
<EMAIL_ADDRESS>and  Yanqiu Wang  School of Mathematical Sciences and
Jiangsu Key Laboratory for NSLSCS, Nanjing Normal University, Nanjing, China
<EMAIL_ADDRESS>
###### Abstract.

This article concerns the weak Galerkin mixed finite element method (WG-MFEM)
for second order elliptic equations on 2D domains with curved boundary. The
Neumann boundary condition is considered since it becomes the essential
boundary condition in this case. It is well-known that the discrepancy between
the curved physical domain and the polygonal approximation domain leads to a
loss of accuracy for discretization with polynomial order $\alpha>1$. The
purpose of this paper is two-fold. First, we present a detailed error analysis
of the original WG-MFEM for solving problems on curved domains, which exhibits
an $O(h^{1/2})$ convergence for all $\alpha\geq 1$. It is a little surprising
to see that even the lowest-order WG-MFEM ($\alpha=1$) experiences a loss of
accuracy. This is different from known results for the finite element method
(FEM) or the mixed FEM, and appears to be a combined effect of the WG-MFEM
design and the fact that the outward normal vector on the polygonal
approximation domain is different from the one on the curved domain. Second,
we propose a remedy to bring the approximation rate back to optimal by
employing two techniques. One is a specially designed boundary correction
technique. The other is to take full advantage of the nice feature that weak
Galerkin discretization can be defined on polygonal meshes, which allows the
curved boundary to be better approximated by multiple short edges without
increasing the total number of mesh elements. Rigorous analysis shows that a
combination of the above two techniques renders optimal convergence for all
$\alpha$. Numerical results further confirm this conclusion.

###### Key words and phrases:

weak Galerkin method, polygonal mesh, curved domain, mixed formulation

###### 2020 Mathematics Subject Classification:

65N15, 65N30

Liu and Wang are supported by the NSFC grant 12171244. Chen is supported by
NSFC grant 12071090.

##  1\. Introduction

Many practical problems arising in science and engineering are posed on
domains with curved boundaries. When such problems are approximated on
polygonal or polyhedral computational domains, the geometric difference
between the two leads to a loss of approximation accuracy [34, 35] for high-
order elements. To resolve this issue, a straight-forward idea is to reduce
the geometric error down to the same level of the approximation error. Popular
methods following this track include the isoparametric finite element method
[21, 25] and the isogeometric analysis [20, 24]. However, due to their
specialized design, neither of them can be applied to meshes consisting of
polygons or polyhedra.

In the past two decades, discretizations on polygonal and polyhedral meshes
have gained considerable attention in the scientific computing community.
Various numerical schemes have been proposed, including the virtual element
method (VEM) (see [2] and references therein), the discontinuous Galerkin
method (DG) [18, 22, 31], and the weak Galerkin method (WG) [38, 39, 30], to
name a few. Very recently, researchers start to apply these discretizations to
curved domains, which requires innovative techniques, with its reason
explained above. One emerging method is the boundary correction technique,
which may have its root date back to a 1972 paper by Bramble, Dupont and
Thomée [10]. The idea is to use normal-directional Taylor expansion, in most
cases just a linear approximation, to correct function values on the boundary.
Burman et. al. [13, 14] proposed the technique for a CutFEM discretization in
2018. It was soon applied to VEM by Bertoluzza et. al. [6]. In 2019, Cheung
et. al. [16] proposed a polynomial extension FEM which is based on an averaged
Taylor expansion. In a series of papers starting from 2018, Main and Scovazzi
[26, 27, 1] designed a shifted boundary FEM on non-fitted meshes, i.e.,
boundary nodes of the mesh may not lie on the curved physical boundary.
However, their method uses 1st-order (linear) Taylor expansion and hence only
works for linear elements. Finally, we mention an earlier but closely related
work [19], where instead of Taylor expansion the authors used a path
integration to achieve a similar ‘boundary correction’ effect.

A totally different track, first proposed for VEM by Beirão da Veiga et. al.
[4] in 2019, is to define the discretization directly on curved mesh elements.
This is possible because of the ‘skeletal’ style design of VEM, where the
degrees of freedom (dofs) in the interior and on the boundary of each mesh
element are defined separately. In [4], dofs as well as related shape
functions on curved boundary edges are defined using the parameter in the
parametric equation of the curved boundary. Hence the shape functions on
curved edges are no longer polynomials in the physical space. Later a
modification was proposed [3] which uses the restriction of physical
polynomials to define shape functions on curved boundary edges. In 2021, Mu
applied the idea to the primal WG discretization [28]. Because of the direct
use of the curved boundary, the implementation of these methods requires
numerical integration formulae on curved mesh elements, as well as a mapping
between each curved boundary segment and its flat counterpart in the
parametric space. The theoretical analysis is also more complicated as it has
to deal with the parametric mapping.

In this paper, we consider the weak Galerkin mixed finite element method (WG-
MFEM) on 2D curved domains. Since its first debut, the WG method has been
quickly applied to various situations [15, 29, 36, 37, 39, 40, 41], among
which [39] focuses on the WG-MFEM for second order elliptic equations and [15]
uses WG-MFEM to solve coupled Darcy-Stokes equations. However, there seems to
be no rigorous analysis of the original WG-MFEM on curved domains. The first
objective of this paper is to fill this gap. Although analysis of
discretizations in the primal formulation on curved domains has long been
well-known, study of the mixed formulation is relatively rare. We shall first
clarify that in either cases, the main difficulty lies in how to impose the
essential boundary condition. For the mixed formulation of second order
elliptic equations, it is the Neumann boundary condition that becomes
essential. A subtlety arises as the Neumann boundary condition involves the
outward normal vector on the boundary, which is different on the curved
boundary and its polygonal approximation. To our knowledge, such a problem, as
well as subsequent mixed-FEM error analysis on curved domains, was first
studied by Bertrand et. al. [7, 8] for the Raviart-Thomas element in 2014.
Both the analysis and numerical results show that the loss of accuracy only
occurs for high-order elements. Therefore it was a little surprising when we
found that the original WG-MFEM yields only an $O(h^{1/2})$ convergence for
all $\alpha$, which means an accuracy loss even for the lowest-order WG-MFEM
discretization. This appears to be a combined effect of the WG-MFEM design and
the outward normal vector issue mentioned above.

We then propose a remedy to bring the approximation rate back to optimal,
which is the second objective of this paper. The remedy employs two
techniques. One is a specially designed boundary correction technique that
treats the difference between outward normal vectors on the curved domain and
the polygonal computational domain. Unlike the boundary correction [13, 14,
16, 26, 27, 1] designed for Dirichlet boundary conditions, ours mainly deals
with the normal component of the flux. The implementation is easy and does not
involve integration on curved regions. But one has to be careful about the
discrete compatibility condition, noticing that the solution to the pure
Neumann boundary problem is not unique. The other technique is a simple
strategy proposed in [33], which uses multiple short straight edges to obtain
a better geometric approximation to the curved boundary. Since the WG
discretizations can be defined on polygonal meshes, the curved boundary is
then better approximated by polygonal elements with multiple short edges.
Mathematically, this multiple short edge approach is still a ‘linear’
approximation and hence the number of short edges required will definitely
increase with $\alpha$. The main advantage of this strategy lies in its
simplicity in the implementation. The algorithm itself remains untouched. One
only needs to provide a new polygonal mesh consisting of multiple short
boundary edges which better approximates the curved domain. We show that
combining the above two techniques and using suitable meshes, the modified WG-
FEM reaches optimal approximation rates for all $\alpha$.

The paper is organized as follows. In Section 2, we introduce the model
problem, the notation and mesh assumptions. In Section 3, we present a
rigorous analysis of the original WG-MFEM on curved domains. The modified WG-
FEM and its theoretical analysis are given in Section 4. Finally, numerical
results are presented in Section 5.

##  2\. Model problem and mesh assumptions

Let $\Omega\subset\mathbb{R}^{2}$ be a general domain with Lipschitz-
continuous and possibly curved boundary $\partial\Omega$. Consider the
Poisson’s equation in its mixed form: Given $g$ in $L^{2}(\Omega)$, find
functions ${\bm{u}}$ and $p$ such that

(2.1) $\displaystyle{\bm{u}}+\nabla p$ $\displaystyle=\bm{0}\quad$
$\displaystyle\text{in }\Omega,$ $\displaystyle\nabla\cdot{\bm{u}}$
$\displaystyle=g\quad$ $\displaystyle\text{in }\Omega.$

Here we conveniently use bold face characters to denote vectors or vector-
valued functions. Equip system (2.1) with the homogeneous Neumann boundary
condition as follows:

(2.2) ${\bm{u}}\cdot\widetilde{\bf n}=0\qquad\text{on }\partial\Omega,$

where $\widetilde{\bf n}$ is the unit outward normal vector on
$\partial\Omega$. System (2.1)-(2.2) is well-posed as long as the following
compatibility condition holds:

(2.3) $\int_{\Omega}g\,\mathrm{d}x=0.$

Moreover, the solution is unique assuming $\int_{\Omega}p\,\mathrm{d}x=0$.

Denote by $H^{m}(D)$, $m\geq 0$, the usual Sobolev space defined on an open
bounded domain $D\subset\Omega$, and endow it with the inner-product
$(\cdot,\cdot)_{m,D}$, the norm $\|\cdot\|_{m,D}$ and the seminorm
$|\cdot|_{m,D}$. Let $H_{0}^{m}(D)=\overline{C_{0}^{\infty}}^{H^{m}(D)}$. When
$m=0$, the space $H^{0}(D)$ is identical to $L^{2}(D)$. In this case, we use
$(\cdot,\cdot)_{D}$ to denote the $L^{2}$ inner-product on $D$. The above
notation also extends to a curve or edge $e$ in $\overline{D}$. Moreover,
denote by $\langle\cdot,\cdot\rangle_{e}$ the duality pair on $e$. Denote by
$H(\text{div},D)$ the space of vector-valued functions with all its components
and divergence in $L^{2}(D)$.

Define spaces

$\displaystyle{\bm{V}}$ $\displaystyle=\\{{\bm{v}}\in
H(\text{div},\Omega)\,|\,{\bm{v}}\cdot\widetilde{\bf n}=0\text{ on
}\partial\Omega\\},$ $\displaystyle\Psi$
$\displaystyle=L^{2}_{0}(\Omega)\triangleq\left\\{q\in
L^{2}(\Omega)\,\bigg{|}\,\int_{\Omega}q\,\mathrm{d}x=0\right\\}.$

Then, the mixed variational formulation of the Poisson’s equation with
homogeneous Neumann boundary condition can be written as: Find
$({\bm{u}},p)\in{\bm{V}}\times\Psi$ such that

(2.4)
$\begin{cases}a({\bm{u}},{\bm{v}})+b({\bm{v}},p)=0\quad&\forall\,{\bm{v}}\in{\bm{V}},\\\
b({\bm{u}},q)=-(g,q)_{\Omega}\quad&\forall\,q\in\Psi,\end{cases}$

where

$a({\bm{u}},{\bm{v}})=({\bm{u}},{\bm{v}})_{\Omega},\qquad
b({\bm{v}},q)=-(\nabla\cdot{\bm{v}},q)_{\Omega}.$

The existence and uniqueness of a weak solution to the mixed problem (2.4) can
be found in [9].

Let $\mathcal{T}_{h}$ be a body-fitted partition of $\Omega$ consisting of
polygons. By ‘body-fitted’ we mean that each boundary edge of
$\mathcal{T}_{h}$ has its two end points lying on $\partial\Omega$. Let
$\overline{\Omega}_{h}=\cup_{K\in\mathcal{T}_{h}}K$. If $\Omega$ is a
polygonal domain in $\mathbb{R}^{2}$, the domains $\Omega_{h}$ and $\Omega$
are identical in the case of a fitted mesh, i.e., $\Omega_{h}=\Omega$. But
when $\Omega$ has curved boundary, $\Omega_{h}$ differs from $\Omega$. It’s
easy to see that when $\Omega$ is convex, we have $\Omega_{h}\subset\Omega$.
In general, $\Omega_{h}$ can be viewed as an approximation to $\Omega$. For
the sake of simplicity, we only present theoretical analysis for the case when
$\Omega$ is convex. Generalization to non-convex domains can be done but
requires a non-trivial use of extension operators. In Section 5, a numerical
example in non-convex domain is presented, which exhibits the same behavior as
experiments in convex domains.

Denote by $h_{K}$ the diameter of each element $K\in\mathcal{T}_{h}$, and let
$h=\max_{K\in\mathcal{T}_{h}}h_{K}$. Denote by $K_{0}$ and $\partial K$ the
interior and the boundary of $K$, respectively. Let $\mathcal{E}_{h}$ be the
set of all edges in $\mathcal{T}_{h}$, which are straight edges. For each edge
$e\in\mathcal{E}_{h}$, denote by $h_{e}$ its length. Let $\mathcal{E}^{I}_{h}$
and $\mathcal{E}^{B}_{h}$ be the set of all interior and boundary edges,
respectively, in $\mathcal{T}_{h}$. Denote
$s=\max_{e\in\mathcal{E}_{h}^{B}}h_{e}$. It is obvious that $s\leq h$. Later
we shall allow $s$ to be much smaller than $h$ in order to get a more accurate
approximation to the curved boundary of $\Omega$.

Note that each $e\in\mathcal{E}_{h}^{B}$ has its two end points lying on
$\partial\Omega$. Denote by $\tilde{e}\subset\partial\Omega$ the (short)
section intersected by these two end points, and by $M_{e}$ the crescent-
shaped region surrounded by $e$ and $\tilde{e}$, as shown in Figure 1. It is
possible that part of $\partial\Omega$ is indeed flat. In this case
$\tilde{e}=e$ and $M_{e}$ is just empty. We use $\partial\Omega_{h}$ to denote
the boundary of $\Omega_{h}$, which consists of all $e\in\mathcal{E}_{h}^{B}$.
Recall that $\widetilde{\bf n}$ is the unit outward normal vector on
$\partial\Omega$, i.e., on curved edge $\tilde{e}$. Denote by ${\bf n}$ the
unit outward normal vector on $\partial\Omega_{h}$, i.e., on straight edge
$e\in\mathcal{E}_{h}^{B}$. Later we also need unit outward normal vectors on
$\partial K$ for each $K\in\mathcal{T}_{h}$. Since $\partial K$ consists of
straight edges, they are still denoted by ${\bf n}$ and should not bring any
ambiguity.

Denote by $\mathcal{T}_{h}^{B}$ all mesh elements containing at least one edge
in $\mathcal{E}_{h}^{B}$. It is possible that an element
$K\in\mathcal{T}_{h}^{B}$ contains multiple edges in $\mathcal{E}_{h}^{B}$.

Following [31, 33], we introduce a set of regularity assumptions on the
partition $\mathcal{T}_{h}$.

A1. There exists a positive constant $C_{1}$ such that each element
$K\in\mathcal{T}_{h}$ is star-shaped with respect to a ball
$\mathcal{B}_{K}\subset K$ with radius $\rho_{K}$ satisfying

$\rho_{K}\geq C_{1}h_{K}.$

A2. There exists a positive constant $C_{2}$ such that each element
$K\in\mathcal{T}_{h}$ has at least one edge $e$ with length

$h_{e}\geq C_{2}h_{K}.$

A3. The mesh is quasi-uniform, that is, there exists a positive constant
$C_{3}$ such that for every element $K\in\mathcal{T}_{h}$ one has

$1\leq\dfrac{h}{h_{K}}\leq C_{3}.$

A4. There exist positive constants $C_{4}$ and $C_{4}^{\prime}$ such that: for
each edge $e\in\mathcal{E}_{h}^{B}$ and corresponding $\tilde{e}$, there is a
one-to-one map $\bm{\Gamma}:e\rightarrow\tilde{e}$ defined in a local
coordinate system $\hat{x}$-$\hat{y}$ (see Figure 1) by

$\bm{\Gamma}(\hat{x},0)=(\hat{x},\gamma(\hat{x}))\in\tilde{e}\qquad\forall\,\hat{x}\in[0,h_{e}],$

where the function $\gamma$ satisfies

$\sup_{\hat{x}\in[0,h_{e}]}|\gamma(\hat{x})-\hat{x}|\leq C_{4}s^{2},$

and

$\sup_{\hat{x}\in[0,h_{e}]}|{\bf n}(\hat{x},0)-\widetilde{\bf
n}(\hat{x},\gamma(\hat{x}))|\leq C_{4}^{\prime}s.$

Figure 1. The map $\bm{\Gamma}$.

A5. The boundary edges in $\mathcal{E}_{h}^{B}$ are quasi-uniform, i.e., there
exists a positive constant $C_{5}$ such that

$s\leq C_{5}\min_{e\in\mathcal{E}_{h}^{B}}h_{e}.$

A6. There exists a positive constant $C_{6}$ such that: for every
$e\in\mathcal{E}_{h}^{B}$ and the unique polygon $K\in\mathcal{T}_{h}^{B}$
having $e$ as an edge, one can draw a triangle $P(e)$ with base $e$ and the
center of ball $\mathcal{B}_{K}$ (defined in assumption A1) as apex; the
height $h_{K,e}$ of $P(e)$ obviously satisfies $\rho_{K}\leq h_{K,e}\leq
h_{K}$; denote by $\hat{{\bm{x}}}$ the local coordinate system with $e$ and
the height of $P(e)$ as abscissa axis and ordinate axis, and define a linear
transformation $F$ by

$F(\hat{{\bm{x}}})=\left(\begin{aligned} h_{e}^{-1}&&0\\\ 0&&h^{-1}_{K,e}\\\
\end{aligned}\right)\hat{{\bm{x}}},$

then the triangle $\mathcal{O}=F(P(e))$ satisfies

  * •

$\text{diam}(\mathcal{O})=O(1)$;

  * •

The radius ratio of the circumscribed circle and the inscribed circle of
$\mathcal{O}$ is less than or equal to $C_{6}$.

We further assume that triangles $P(e)$ from all $e\in\partial
K\cap\partial\Omega_{h}$ form a finite overlapping, in the sense that each
point in $K$ can only be covered by no more than $M$ such triangles.

######  Remark 2.1.

A1-A2 are polygonal shape regularity conditions proposed in [39, 31]. A3 is
the quasi-uniform assumption. These three are standard mesh assumptions used
in WG discretizations. It is well-known that A4 always holds as long as
$\partial\Omega$ is piecewise $C^{2}$ continuous [35, 7]. Conditions A5-A6 are
extra requirements for the case of curved domains, which were first proposed
in [33]. We shall see later how these conditions are used in the analysis.

######  Remark 2.2.

Note that conditions A1-A6 allow the existence of small edges, and hence allow
the approximation of the curved boundary by multiple straight short edges.

######  Remark 2.3.

A less obvious but important consequence of A5-A6 is that, each
$K\in\mathcal{T}_{h}^{B}$ contains at most $O(\frac{h}{s})$ edges in
$\mathcal{E}_{h}^{B}$, i.e., boundary edges in $\mathcal{E}_{h}^{B}$ should
never be in a zigzag formation.

For the sake of brevity, throughout this article, we write $x\lesssim y$ and
$x\gtrsim y$ in place of $x\leq Cy$ and $x\geq Cy$, respectively, for a
positive constant $C$ independent of the discretization parameters. Moreover,
$x\approx y$ means that there exist positive constants $c,C$ such that $cy\leq
x\leq Cy$. When required, the constants will be written explicitly.

In the entire paper, we always assume that the mesh $\mathcal{T}_{h}$
satisfies assumptions A1-A6. This guarantees the trace inequality, the inverse
inequality and a few other important lemmas on the curved domains.

######  Lemma 2.4.

(Trace inequality [31, 33, 39]) For $K\in\mathcal{T}_{h}$ and $v\in H^{1}(K)$,
one has

$\displaystyle\|v\|_{0,\partial K}^{2}\lesssim h^{-1}\|v\|_{0,K}^{2}+h\|\nabla
v\|_{0,K}^{2}.$

Given a non-negative integer $j$, we denote by $P_{j}(K)$ the space of
polynomials with degree less than or equal to $j$ on $K\in\mathcal{T}_{h}$.

######  Lemma 2.5.

(Inverse inequality [31, 33, 39]) For $K\in\mathcal{T}_{h}$ and $v\in
P_{j}(K)$, one has

$\displaystyle\|\nabla v\|_{0,K}\lesssim h^{-1}\|v\|_{0,K},$

where the hidden constant in $\lesssim$ may depend on $j$ but not on $h$ or
the shape of $K$.

Note that if $K\in\mathcal{T}_{h}^{B}$, polynomials in $P_{j}(K)$ can be
naturally extend to any adjacent crescent-shaped region $M_{e}$, or vice
versa.

######  Lemma 2.6.

For $K\in\mathcal{T}_{h}^{B}$ and $v\in P_{j}(K)$, one has

$\sum_{e\in\partial K\cap\partial\Omega_{h}}\|v\|_{0,M_{e}}^{2}\lesssim
h^{-1}s^{2}\|v\|_{0,K}^{2}.$

###### Proof.

Take the center of $\mathcal{B}_{K}$, as defined in Assumption A1, and draw a
ball with radius $h_{K}+C_{4}s^{2}$. Denote the larger ball by
$\mathcal{B}_{K}^{\prime}$. Then, using Assumption A4 and the fact that each
$K\in\mathcal{T}_{h}^{B}$ contains at most $O(\frac{h}{s})$ edges in
$\mathcal{E}_{h}^{B}$ (see Remark 2.3), we have

$\sum_{e\in\partial
K\cap\partial\Omega_{h}}\|v\|_{0,M_{e}}^{2}\lesssim\sum_{e\in\partial
K\cap\partial\Omega_{h}}s^{3}\|v\|_{L^{\infty}(\mathcal{B}_{K}^{\prime})}^{2}\\\
\lesssim\frac{h}{s}s^{3}\|v\|_{L^{\infty}(\mathcal{B}_{K}^{\prime})}^{2}.$

Finally, using Assumption A1, a scaling argument and the fact that $s\leq
h\leq O(1)$, we have

$\|v\|_{L^{\infty}(\mathcal{B}_{K}^{\prime})}^{2}\lesssim\|v\|_{L^{\infty}(\mathcal{B}_{K})}^{2}\lesssim
h^{-2}\|v\|_{0,\mathcal{B}_{K}}^{2}\leq h^{-2}\|v\|_{0,K}^{2}.$

This completes the proof of the lemma. ∎

Following [11], we have

######  Lemma 2.7.

For $v\in H^{1}(\Omega)$, one has

(2.5) $\sum_{e\in\mathcal{E}_{h}^{B}}\|v\|_{0,M_{e}}^{2}\lesssim
s^{2}\|v\|_{1,\Omega}^{2}.$

Moreover, when $v\in H_{0}^{1}(\Omega)$, for each $e\in\mathcal{E}_{h}^{B}$
one has

(2.6) $\|v\|_{0,M_{e}}\lesssim
s^{2}\|v\|_{1,M_{e}},\qquad\qquad\|v\|_{0,e}\lesssim s\|v\|_{1,M_{e}}.$

###### Proof.

Inequality (2.9) in [11] states that

$\displaystyle\|v\|_{0,M_{e}}^{2}\lesssim
s^{2}\|v\|_{0,\tilde{e}}^{2}+s^{4}\|v\|_{1,M_{e}}^{2}.$

Taking the summation and applying the trace inequality on $\Omega$, we get
(2.5). For $v\in H_{0}^{1}(\Omega)$, the first inequality in (2.6) follows
immediately from $\|v\|_{0,\partial\Omega}=0$. The second inequality in (2.6)
is exactly inequality (2.11) in [11]. ∎

The following version of the Bramble-Hilbert lemma has been proved in Chapter
4 of [12].

######  Lemma 2.8.

(Bramble-Hilbert) For $K\in\mathcal{T}_{h}$ and $v\in H^{m}(K)$, there exists
an averaged Taylor polynomial $v^{m}_{K}$ of degree less than or equal to
$m-1$ satisfying

$\displaystyle|v-v^{m}_{K}|_{s,K}\lesssim
h_{K}^{m-s}|v|_{m,K},\quad\text{for}\ s\ =\ 0,\ldots,m.$

where the hidden constant in $\lesssim$ may depend on $m$ but not on $h$ or
the shape of $K$.

##  3\. Weak Galerkin discretization

Now we introduce the WG-MFEM discretization for system (2.1)-(2.2) in a form
presented in [15]. On each $K\in\mathcal{T}_{h}$, denote by $P_{j}(K_{0})$ or
$P_{j}(K)$ the set of polynomials with degree less than or equal to $j$.
Likewise, on each $e\in\mathcal{E}_{h}$, let $P_{j}(e)$ be the set of
polynomials of degree no more than $j$. We define the weak Galerkin spaces

$\displaystyle\bm{W}_{h}=\big{\\{}$
$\displaystyle{\bm{v}}=\\{{\bm{v}}_{0},{\bm{v}}_{b}\\}\in[L^{2}(\Omega_{h})]^{2}\times[L^{2}(\mathcal{E}_{h})]^{2}\text{
such that}$
$\displaystyle{\bm{v}}_{0}|_{K_{0}}\in[P_{\alpha}(K_{0})]^{2}\text{ for
}K\in\mathcal{T}_{h},$ $\displaystyle{\bm{v}}_{b}|_{e}=v_{b}{\bf n}_{e}\text{,
where }v_{b}\in P_{\beta}(e)\text{, for }e\in\mathcal{E}_{h}\big{\\}},$
$\displaystyle{\bm{V}}_{h}=\big{\\{}$
$\displaystyle{\bm{v}}\in\bm{W}_{h}\textrm{ satisfying
}{\bm{v}}_{b}|_{e}=0\text{ for }e\in\mathcal{E}_{h}^{B}\big{\\}},$

where $\alpha$ and $\beta$ are given non-negative integers and ${\bf n}_{e}$
is a prescribed normal direction on each edge $e\in\mathcal{E}_{h}$. A key
feature of the weak Galerkin discretization is that a function
${\bm{v}}\in{\bm{V}}_{h}$ takes separate values ${\bm{v}}_{0}$ on the interior
of each $K$ and ${\bm{v}}_{b}$ on edges. Define

$\Psi_{h}=\big{\\{}q\in L^{2}_{0}(\Omega_{h})\text{ such that }q|_{K}\in
P_{\sigma}(K)\ \text{for}\ K\in\mathcal{T}_{h}\big{\\}},$

where $\sigma$ is a given non-negative integer. Moreover, assume that

(3.1) $\beta-1\leq\sigma\leq\beta=\alpha.$

Condition (3.1) is imposed to make sure that the weak Galerkin discretization
has desired stability and approximation properties, as will become clear in
the analysis to be given later.

On each $K\in\mathcal{T}_{h}$, define the weak divergence
$\nabla_{w}\cdot{\bm{v}}\in P_{\beta}(K)$ for ${\bm{v}}\in{\bm{V}}_{h}$ by

$(\nabla_{w}\cdot{\bm{v}},q)_{K}=-({\bm{v}}_{0},\nabla
q)_{K}+\langle{\bm{v}}_{b}\cdot{\bf n},q\rangle_{\partial
K}\qquad\forall\,q\in P_{\beta}(K).$

Define the discrete bilinear forms
$a_{h}:\>{\bm{V}}_{h}\times{\bm{V}}_{h}\to\mathbb{R}$ and
$b_{h}:\>{\bm{V}}_{h}\times\Psi_{h}\to\mathbb{R}$ by

$\displaystyle a_{h}({\bm{u}},{\bm{v}})$
$\displaystyle=\left({\bm{u}}_{0},{\bm{v}}_{0}\right)_{\Omega_{h}}+\rho\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\left\langle\left({\bm{u}}_{0}-{\bm{u}}_{b}\right)\cdot{\bf
n},\left({\bm{v}}_{0}-{\bm{v}}_{b}\right)\cdot{\bf n}\right\rangle_{\partial
K},$ $\displaystyle b_{h}({\bm{v}},q)$
$\displaystyle=-(\nabla_{w}\cdot{\bm{v}},q)_{\Omega_{h}},$

in which $\rho$ is a positive constant.

######  Remark 3.1.

The second part in $a_{h}(\cdot,\cdot)$ is a stabilization term. However,
different from the discontinuous Galerkin method, the stabilization parameter
$\rho$ in the weak Galerkin discretization can be chosen arbitrarily without
affecting the approximation results. In practice, one can simply set $\rho=1$.
For convenience, we denote the stabilization part by

(3.2)
$c({\bm{u}},{\bm{v}})=\rho\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\left\langle\left({\bm{u}}_{0}-{\bm{u}}_{b}\right)\cdot{\bf
n},\left({\bm{v}}_{0}-{\bm{v}}_{b}\right)\cdot{\bf n}\right\rangle_{\partial
K}.$

The weak Galerkin formulation for system (2.1)-(2.2) can now be written as
follows: Find ${\bm{u}}_{h}\in{\bm{V}}_{h}$ and $p_{h}\in\Psi_{h}$ such that

(3.3)
$\begin{cases}a_{h}({\bm{u}}_{h},{\bm{v}})+b_{h}({\bm{v}},p_{h})=0\quad&\forall\,{\bm{v}}\in{\bm{V}}_{h},\\\
b_{h}({\bm{u}}_{h},q)=-(g,q)_{\Omega_{h}}\quad&\forall\,q\in\Psi_{h}.\end{cases}$

######  Remark 3.2.

Recall that for the continuous problem (2.1)-(2.2) to be well-posed, a
compatibility condition (2.3) is necessary which requires $g$ to be mean-value
free on $\Omega$. One may wonder whether the discrete problem requires a
similar compatibility condition or not. We point out that the compatibility
mechanism works differently for finite-dimensional problems and the key is to
have

(3.4) $b_{h}({\bm{v}},1)=0\qquad\forall\,{\bm{v}}\in{\bm{V}}_{h},$

which is obviously true in our case according to the definitions of
$b_{h}(\cdot,\cdot)$ and $\nabla_{w}\cdot$. Note that $(1,q)_{\Omega_{h}}=0$
for all $q\in\Psi_{h}$. Hence the second equation in (3.3) is the same as

(3.5) $b_{h}({\bm{u}}_{h},q)=-(\tilde{g},q)_{\Omega_{h}}\quad\forall
q\in\Psi_{h},$

where $\tilde{g}=g-\frac{1}{|\Omega_{h}|}\int_{\Omega_{h}}g\,\mathrm{d}x$ is
mean-value free on $\Omega_{h}$. Equation (3.5) is now ‘compatible’ in the
traditional sense as it also holds for $q\equiv const$. From the theoretical
point of view, there is no difference should one choose to use (3.3) or to
replace its second equation by (3.5).

For finite-dimensional problems, what the ‘compatibility’ condition may
actually affect is the implementation procedure. In practice, it is not
convenient to compute a basis for $\Psi_{h}$, which needs to be mean-value
free on $\Omega_{h}$. One usually drops the mean-value free condition while at
the same time expecting the resulting stiffness matrix $M$ to have a rank $1$
deficiency (under the discrete inf-sup condition to be proved later). Now, if
the ‘compatible’ equation (3.5) is used in the implementation, it guarantees
that the right-hand side vector of the linear system is orthogonal to
$ker(M^{T})$. Hence the system is solvable, i.e., compatible. Using elementary
linear algebra, one immediately sees that even if the ‘non-compatible’ second
equation of (3.3) is used in the implementation, it just yields a right-hand
side vector not orthogonal to $ker(M^{T})$. In this case, we know that the
kernel consists of exactly $q\equiv const$. Therefore a pure algebraic post-
process after assembling the entire linear system, i.e., making the right-hand
vector orthogonal to $ker(M^{T})$, can easily resolve this issue and render
the linear system ‘compatible’.

Due to the above explanation, we do not need to worry about the
‘compatibility’ of the discrete system either theoretically or in the
implementation, as long as (3.4) holds.

To analyze the well-posedness and approximation properties of the weak
Galerkin discretization (3.3), we first define the norm on ${\bm{V}}_{h}$ as
follows:

$\displaystyle\|{\bm{v}}\|_{{\bm{V}}_{h}}=$
$\displaystyle\left(\|{\bm{v}}_{0}\|_{\Omega_{h}}^{2}+\rho\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\left\|\left({\bm{v}}_{0}-{\bm{v}}_{b}\right)\cdot{\bf
n}\right\|_{\partial
K}^{2}\right)^{1/2}\quad\forall\,{\bm{v}}\in{\bm{V}}_{h}.$

It is obvious that $a_{h}({\bm{v}},{\bm{v}})=\|{\bm{v}}\|_{{\bm{V}}_{h}}^{2}$
for all ${\bm{v}}\in{\bm{V}}_{h}$. We shall show that $\|\cdot\|_{\bm{V_{h}}}$
is a well-defined norm. Indeed, if $\|{\bm{v}}\|_{{\bm{V}}_{h}}=0$ for some
${\bm{v}}\in{\bm{V}}_{h}$; i.e.,

$\left({\bm{v}}_{0},{\bm{v}}_{0}\right)_{\Omega_{h}}+\rho\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\left\langle\left({\bm{v}}_{0}-{\bm{v}}_{b}\right)\cdot{\bf
n},\left({\bm{v}}_{0}-{\bm{v}}_{b}\right)\cdot{\bf n}\right\rangle_{\partial
K}=0,$

one has ${\bm{v}}_{0}\equiv\bm{0}$ on each element $K$ and
$({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf n}=0$ on each edge
$e\in\mathcal{E}_{h}$. This leads to $0={\bm{v}}_{b}\cdot{\bf n}=v_{b}{\bf
n}_{e}\cdot{\bf n}$ and consequently ${\bm{v}}_{b}\equiv\bm{0}$ on each edge
$e\in\mathcal{E}_{h}$. Therefore $\|\cdot\|_{\bm{V_{h}}}$ is a norm on
${\bm{V}}_{h}$.

On each $K\in\mathcal{T}_{h}$, denote by ${\bm{Q}}_{0},\ \mathbb{Q}_{h}$ and
$\pi_{h}$ the $L^{2}$ projections onto $\left[P_{\alpha}(K)\right]^{2}$,
$P_{\sigma}(K)$ and $P_{\beta}(K)$, respectively. On each
$e\in\mathcal{E}_{h}$, denote by $Q_{b}$ the $L^{2}$ projection onto
$P_{\beta}(e)$. On the entire $\Omega_{h}$, we use the same notation to denote
the combination of the above local projections. For
${\bm{v}}\in[H^{1}(\Omega_{h})]^{2}$, define
${\bm{Q}}_{b}{\bm{v}}=(Q_{b}({\bm{v}}\cdot{\bf n}_{e})){\bf n}_{e}$ for each
$e\in\mathcal{E}_{h}$. Then we define a projection
${\bm{Q}}_{h}:\>[H^{1}(\Omega_{h})]^{2}\to\bm{W}_{h}$ by

${\bm{Q}}_{h}{\bm{v}}=\\{{\bm{Q}}_{0}{\bm{v}},{\bm{Q}}_{b}{\bm{v}}\\}.$

One has

######  Lemma 3.3.

(Lemma 3.12 in [15]) Assume that $\alpha$ and $\beta$ satisfy (3.1), the
following commutative property holds:

(3.6) $\displaystyle\nabla_{w}\cdot\left({\bm{Q}}_{h}{\bm{v}}\right)$
$\displaystyle=\pi_{h}(\nabla\cdot{\bm{v}})\qquad\forall\,{\bm{v}}\in[H^{1}(\Omega_{h})]^{2}.$

The above commutative property is one of the key features of the weak Galerkin
discretization and has played an important role in its theoretical analysis
[15, 39]. However, the situation is quite different when $\Omega$ is a curved
domain, because for ${\bm{v}}\in[H^{1}(\Omega)]^{2}$ satisfying
${\bm{v}}\cdot\widetilde{\bf n}|_{\partial\Omega}=0$ one does not get
${\bm{Q}}_{h}{\bm{v}}\cdot{\bf n}|_{\partial\Omega_{h}}=0$. The projection
${\bm{Q}}_{h}$ only maps ${\bm{v}}$ into $\bm{W}_{h}$, but not ${\bm{V}}_{h}$.

To partly remedy this, we introduce

$\widetilde{\bm{Q}}_{b}{\bm{v}}=\begin{cases}\bm{0}\quad&\textrm{on
}e\in\mathcal{E}_{h}^{B},\\\ {\bm{Q}}_{b}{\bm{v}}&\textrm{on all other
edges},\end{cases}$

and define a modified projection
$\widetilde{\bm{Q}}_{h}:\>[H^{1}(\Omega_{h})]^{2}\to{\bm{V}}_{h}$ by

$\widetilde{\bm{Q}}_{h}{\bm{v}}=\\{{\bm{Q}}_{0}{\bm{v}},\widetilde{\bm{Q}}_{b}{\bm{v}}\\}.$

Obviously, the price to pay is that $\widetilde{\bm{Q}}_{h}$ no longer
satisfies the highly desired commutative property (3.6). This brings a lot
trouble to the theoretical analysis of the weak Galerkin formulation (3.3), as
we will see later.

Using the Bramble-Hilbert lemma 2.8, we easily get the following result:

######  Lemma 3.4.

Let ${\bm{u}}\in[H^{r+1}(\Omega)]^{2}$ and $p\in H^{t+1}(\Omega)$, where
$0\leq r\leq\alpha$ and $0\leq t\leq\sigma$. Then

$\displaystyle\sum_{K\in\mathcal{T}_{h}}\|{\bm{u}}-{\bm{Q}}_{0}{\bm{u}}\|^{2}_{K}+\sum_{K\in\mathcal{T}_{h}}h_{K}^{2}\|\nabla({\bm{u}}-{\bm{Q}}_{0}{\bm{u}})\|^{2}_{K}\lesssim
h^{2(r+1)}|{\bm{u}}|^{2}_{r+1,\Omega_{h}},$
$\displaystyle\sum_{K\in\mathcal{T}_{h}}\|p-\mathbb{Q}_{h}p\|^{2}_{K}+\sum_{K\in\mathcal{T}_{h}}h_{K}^{2}\|\nabla(p-\mathbb{Q}_{h}p)\|^{2}_{K}\lesssim
h^{2(t+1)}|p|^{2}_{t+1,\Omega_{h}}.$

###  3.1.  Existence and uniqueness of the discrete solution

Since $a_{h}(\cdot,\cdot)$ is indeed an inner-product on ${\bm{V}}_{h}$,
according to the standard theory of mixed finite element methods [9], System
(3.3) admits a unique solution as long as the following discrete inf-sup
condition holds:

(3.7)
$\sup_{{\bm{v}}\in{\bm{V}}_{h}}\frac{|b_{h}({\bm{v}},q)|}{\|{\bm{v}}\|_{{\bm{V}}_{h}}}\gtrsim\|q\|_{0,\Omega_{h}}\qquad\forall\,q\in\Psi_{h}.$

The standard way to prove the discrete inf-sup condition is to use the
continuous inf-sup condition together with a stable projection onto the
discrete space, as is the case in [39]. However, following this strategy is
not easy on domains with curved boundary. Below we shall explain why. Firstly,
the continuous inf-sup condition must be defined on $\Omega$ instead of
$\Omega_{h}$ so that the constant in it does not depend on $h$ or $s$.
However, a function $q\in\Psi_{h}$ is mean value free on $\Omega_{h}$ but not
necessarily on $\Omega$, and hence can not be used directly in the continuous
inf-sup condition. Secondly, as pointed out earlier, the projection
${\bm{Q}}_{h}$ satisfies the commutative property (3.6) but only maps
functions into $\bm{W}_{h}$, while $\widetilde{\bm{Q}}_{h}$ maps functions
into ${\bm{V}}_{h}$ at the price of violating the commutative property.
Neither of them can serve directly in a traditional proof of the discrete inf-
sup condition. In fact, the proof of the inf-sup condition (3.7) turns out to
be non-trivial:

######  Lemma 3.5.

The discrete $\inf$-$\sup$ condition (3.7) holds when $h$ is sufficiently
small.

###### Proof.

Each $q\in\Psi_{h}$ is a function defined on $\Omega_{h}$. Note that it can be
naturally extended to $\Omega$ by filling the gap $\Omega\backslash\Omega_{h}$
with the same polynomial values on neighboring mesh elements. For simplicity,
we still denote this extension by $q$. Now $q$ is mean-value free on
$\Omega_{h}$. Define
$\bar{q}=q-\frac{1}{|\Omega|}\int_{\Omega\setminus\Omega_{h}}q\,\mathrm{d}x$.
Obviously $\bar{q}$ is mean-value free on $\Omega$. Moreover, by the triangle
inequality, Lemma 2.6, and the fact that $s\leq h\leq O(1)$, one has

$\displaystyle\|\bar{q}\|_{0,\Omega}^{2}$
$\displaystyle=\bigg{(}\|q\|_{0,\Omega_{h}}^{2}+\sum_{e\in\mathcal{E}_{h}^{B}}\|q\|_{0,M_{e}}^{2}\bigg{)}+\frac{1}{|\Omega|}\left(\int_{\Omega\setminus\Omega_{h}}q\,\mathrm{d}x\right)^{2}$
$\displaystyle\lesssim\|q\|_{0,\Omega_{h}}^{2}+h^{-1}s^{2}\|q\|_{0,\Omega_{h}}^{2}+|\Omega\backslash\Omega_{h}|\int_{\Omega\backslash\Omega_{h}}|q|^{2}\,\mathrm{d}x$
$\displaystyle=\|q\|_{0,\Omega_{h}}^{2}+h^{-1}s^{2}\|q\|_{0,\Omega_{h}}^{2}+|\Omega\backslash\Omega_{h}|\sum_{e\in\mathcal{E}_{h}^{B}}\|q\|_{0,M_{e}}^{2}$
$\displaystyle\lesssim\|q\|_{0,\Omega_{h}}^{2}+h^{-1}s^{2}\|q\|_{0,\Omega_{h}}^{2}\lesssim\|q\|_{0,\Omega_{h}}^{2}.$

It can be shown (see, e.g.,[23]) that there exists a
$\bm{w}\in\left[H_{0}^{1}(\Omega)\right]^{2}$ such that
$\nabla\cdot\bm{w}=\bar{q}$ and
$\|\bm{w}\|_{1,\Omega}\lesssim\|\bar{q}\|_{0,\Omega}\lesssim\|q\|_{0,\Omega_{h}}$.
By Lemma 3.3 and the fact that $\sigma\leq\beta,$ we have

$\left(\nabla_{w}\cdot{\bm{Q}}_{h}\bm{w},q\right)_{\Omega_{h}}=\left(\pi_{h}(\nabla\cdot\bm{w}),q\right)_{\Omega_{h}}=(\nabla\cdot\bm{w},q)_{\Omega_{h}}=(\bar{q},q)_{\Omega_{h}}=\|q\|_{0,\Omega_{h}}^{2}.$

Combine the above equation with the definitions of ${\bm{Q}}_{h}$,
$\widetilde{\bm{Q}}_{h}$ and $\nabla_{w}\cdot$, one gets

$\displaystyle\left|\left(\nabla_{w}\cdot\widetilde{\bm{Q}}_{h}\bm{w},q\right)_{\Omega_{h}}\right|$
$\displaystyle=\left|\left(\nabla_{w}\cdot{\bm{Q}}_{h}\bm{w},q\right)_{\Omega_{h}}-\sum_{K\in\mathcal{T}_{h}^{B}}\langle{\bm{Q}}_{b}\bm{w}\cdot{\bf
n},q\rangle_{\partial K\cap\partial\Omega_{h}}\right|$
$\displaystyle\geq\|q\|_{0,\Omega_{h}}^{2}-\left(\sum_{K\in\mathcal{T}_{h}^{B}}h_{K}^{-1}\|{\bm{Q}}_{b}\bm{w}\cdot{\bf
n}\|_{0,\partial
K\cap\partial\Omega_{h}}^{2}\right)^{1/2}\left(\sum_{K\in\mathcal{T}_{h}^{B}}h_{K}\|q\|_{0,\partial
K\cap\partial\Omega_{h}}^{2}\right)^{1/2}.$

Check the right-hand side of the above inequality. First, by the trace
inequality and the inverse inequality (lemmas 2.4-2.5), it is clear that

(3.8) $\left(\sum\limits_{K\in\mathcal{T}_{h}^{B}}h_{K}\|q\|_{0,\partial
K\cap\partial\Omega_{h}}^{2}\right)^{1/2}\lesssim\|q\|_{0,\Omega_{h}}.$

Next, using $\bm{w}|_{\partial\Omega}=\bm{0}$ and Lemma 2.7, we have

(3.9)
$\displaystyle\left(\sum_{K\in\mathcal{T}_{h}^{B}}h_{K}^{-1}\|{\bm{Q}}_{b}\bm{w}\cdot{\bf
n}\|_{0,\partial K\cap\partial\Omega_{h}}^{2}\right)^{1/2}$
$\displaystyle\leq\left(\sum_{K\in\mathcal{T}_{h}^{B}}h_{K}^{-1}\|\bm{w}\cdot{\bf
n}\|_{0,\partial K\cap\partial\Omega_{h}}^{2}\right)^{1/2}$
$\displaystyle\lesssim h^{-\frac{1}{2}}s\|\bm{w}\|_{1,\Omega}\lesssim
h^{-\frac{1}{2}}s\|q\|_{0,\Omega_{h}}.$

Hence when $h$ is sufficiently small,

$\displaystyle\left|(\nabla_{w}\cdot\widetilde{\bm{Q}}_{h}\bm{w},q)\right|\geq\left(1-O(h^{-\frac{1}{2}}s)\right)\|q\|^{2}_{0,\Omega_{h}}\geq\left(1-O(h^{\frac{1}{2}})\right)\|q\|^{2}_{0,\Omega_{h}}\gtrsim\|q\|^{2}_{0,\Omega_{h}}.$

Using the fact that ${\bm{Q}}_{0}{\bm{w}}={\bm{Q}}_{b}{\bm{Q}}_{0}{\bm{w}}$
and the approximation property of ${\bm{Q}}_{0}$, one gets

$\displaystyle
h_{K}^{-1}\left\|\left({\bm{Q}}_{0}\bm{w}-{\bm{Q}}_{b}\bm{w}\right)\cdot{\bf
n}\right\|_{0,\partial K}^{2}$ $\displaystyle\leq
h_{K}^{-1}\left\|{\bm{Q}}_{0}\bm{w}-{\bm{Q}}_{b}\bm{w}\right\|_{0,\partial
K}^{2}$ $\displaystyle\leq
h_{K}^{-1}\left\|{\bm{Q}}_{0}\bm{w}-\bm{w}\right\|_{0,\partial K}^{2}$
$\displaystyle\lesssim\|\nabla\bm{w}\|_{0,K}^{2}.$

We then have

$\displaystyle\|\widetilde{\bm{Q}}_{h}\bm{w}\|_{{\bm{V}}_{h}}=$
$\displaystyle\left(\|{\bm{Q}}_{0}\bm{w}\|_{0,\Omega_{h}}^{2}+\rho\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\|\left({\bm{Q}}_{0}\bm{w}-\widetilde{\bm{Q}}_{b}\bm{w}\right)\cdot{\bf
n}\|_{0,\partial K}^{2}\right)^{1/2}$ $\displaystyle\leq$
$\displaystyle\left(\|{\bm{Q}}_{0}\bm{w}\|_{0,\Omega_{h}}^{2}+\rho\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\|\left({\bm{Q}}_{0}\bm{w}-{\bm{Q}}_{b}\bm{w}\right)\cdot{\bf
n}\|_{0,\partial K}^{2}\right)^{1/2}$
$\displaystyle+\left(\rho\sum_{K\in\mathcal{T}_{h}^{B}}h_{K}^{-1}\|{\bm{Q}}_{b}\bm{w}\cdot{\bf
n}\|_{0,\partial K\cap\partial\Omega_{h}}^{2}\right)^{1/2}$
$\displaystyle\lesssim\|{\bm{w}}\|_{1,\Omega}+h^{-\frac{1}{2}}s\|{\bm{w}}\|_{1,\Omega}$
$\displaystyle\lesssim\|q\|_{0,\Omega_{h}}.$

Combining the above gives

$\sup_{{\bm{v}}\in{\bm{V}}_{h}}\frac{|b_{h}({\bm{v}},q)|}{\|{\bm{v}}\|_{{\bm{V}}_{h}}}\geq\frac{|b_{h}(\widetilde{\bm{Q}}_{h}{\bm{w}},q)|}{\|\widetilde{\bm{Q}}_{h}{\bm{w}}\|_{{\bm{V}}_{h}}}\gtrsim\|q\|_{0,\Omega_{h}}.$

This completes the proof of the lemma. ∎

By the standard theory of the mixed finite elements [9], we know that system
(3.3) admits a unique solution. Moreover, the discrete inf-sup condition
ensures that the discrete operator in (3.3) is stable in the sense that that
unique solution to

$\begin{cases}a_{h}({\bm{u}}_{h},{\bm{v}})+b_{h}({\bm{v}},p_{h})=F({\bm{v}})\quad&\forall{\bm{v}}\in{\bm{V}}_{h},\\\
b_{h}({\bm{u}}_{h},q)=G(q)\quad&\forall q\in\Psi_{h},\end{cases}$

satisfies

(3.10)
$\|{\bm{u}}_{h}\|_{{\bm{V}}_{h}}+\|p_{h}\|_{\Psi_{h}}\lesssim\|F\|_{{\bm{V}}_{h}^{\prime}}+\|G\|_{\Psi_{h}^{\prime}}.$

The stability result (3.10) is essential to the error analysis to be given
next.

###  3.2. Error analysis

To analyze the approximation error of the weak Galerkin discretization (3.3),
we first give two lemmas that help to simplify the derivation of error
equations. The proof of these two lemmas use quite standard techniques and
will be given in Appendix A.

######  Lemma 3.6.

The solution ${\bm{u}}$ and $p$ to system (2.1)-(2.2) satisfy

$a_{h}\left(\widetilde{\bm{Q}}_{h}{\bm{u}},{\bm{v}}\right)+b_{h}\left({\bm{v}},\mathbb{Q}_{h}p\right)=c({\bm{Q}}_{h}{\bm{u}},{\bm{v}})+l_{s}({\bm{v}})-l_{\mathrm{div}}({\bm{v}})\qquad\forall\,{\bm{v}}\in{\bm{v}}_{h},$

where $c(\cdot,\cdot)$ is the stabilization term defined in (3.2) and the
linear functionals $l_{s}(\cdot)$ and $l_{\mathrm{div}}(\cdot)$ are defined by

$\displaystyle l_{s}({\bm{v}})$
$\displaystyle=\rho\sum_{K\in\mathcal{T}_{h}^{B}}h_{K}^{-1}\left\langle{\bm{u}}\cdot{\bf
n},\left({\bm{v}}_{0}-{\bm{v}}_{b}\right)\cdot{\bf n}\right\rangle_{\partial
K\cap\partial\Omega_{h}},$ $\displaystyle l_{\mathrm{div}}({\bm{v}})$
$\displaystyle=\sum_{K\in\mathcal{T}_{h}}\langle({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n},p-\mathbb{Q}_{h}p\rangle_{\partial K}.$

######  Lemma 3.7.

The solution ${\bm{u}}$ to problem (2.1)-(2.2) satisfies

$b_{h}\left(\widetilde{\bm{Q}}_{h}{\bm{u}},q\right)=-(g,q)_{\Omega_{h}}+l_{b}(q)\qquad\forall\,q\in\Psi_{h},$

where the linear functional $l_{b}(\cdot)$ is defined by

$l_{b}(q)=\sum_{K\in\mathcal{T}_{h}^{B}}\langle{{\bm{u}}}\cdot{{\bf
n}},q\rangle_{\partial K\cap\partial\Omega_{h}}.$

With the aid of lemmas 3.6-3.7, we can easily derive the error equations. Let
$({\bm{u}},p)$ be the solution to problem (2.1)-(2.2), and
$({\bm{u}}_{h},p_{h})$ be the solution to the weak Galerkin discretization
(3.3). Define

(3.11)
${\bm{e}}_{{{\bm{u}}}}=\widetilde{\bm{Q}}_{h}{{\bm{u}}}-{\bm{u}}_{h}=\left\\{{\bm{Q}}_{0}{{\bm{u}}}-{{\bm{u}}}_{0},\widetilde{\bm{Q}}_{b}{{\bm{u}}}-{{\bm{u}}}_{b}\right\\},\quad
e_{p}=\mathbb{Q}_{h}p-p_{h}.$

Using lemmas 3.6-3.7, we clearly have

(3.12)
$\begin{cases}a_{h}\left({\bm{e}}_{{{\bm{u}}}},{\bm{v}}\right)+b_{h}\left({\bm{v}},e_{p}\right)=c({\bm{Q}}_{h}{\bm{u}},{\bm{v}})+l_{s}({\bm{v}})-l_{\text{div}}({\bm{v}})\qquad&\forall\,{\bm{v}}\in{\bm{V}}_{h},\\\
b_{h}\left({\bm{e}}_{{{\bm{u}}}},q\right)=l_{b}(q)\quad&\forall\,q\in\Psi_{h}.\end{cases}$

We shall first derive the upper bounds for the right-hand side of (3.12).

######  Lemma 3.8.

Let $u\in{\bm{V}}\cap[H^{r+1}(\Omega)]^{2}$ and $p\in\Psi\cap H^{t+1}(\Omega)$
be the solution to problem (2.1)-(2.2), where $0\leq r\leq\alpha$ and $0\leq
t\leq\sigma$. Then

$\displaystyle|c({{\bm{Q}}}_{h}\bm{{u}},{\bm{v}})|$ $\displaystyle\lesssim
h^{r}\|{\bm{u}}\|_{r+1,\Omega}\|{\bm{v}}\|_{{\bm{V}}_{h}}\quad$
$\displaystyle\forall\,{\bm{v}}\in{\bm{V}}_{h},$
$\displaystyle|l_{\operatorname{div}}({\bm{v}})|$ $\displaystyle\lesssim
h^{t+1}\|p\|_{t+1,\Omega}\|{\bm{v}}\|_{{\bm{V}}_{h}}\quad$
$\displaystyle\forall\,q\in\Psi_{h}.$

###### Proof.

Note that ${\bm{Q}}_{0}{\bm{u}}={\bm{Q}}_{b}({\bm{Q}}_{0}{\bm{u}})$ on each
edge in $\mathcal{T}_{h}$. Hence by the Cauchy-Schwarz inequality, the trace
inequality (Lemma 2.4) and the approximation property of ${\bm{Q}}_{0}$, we
have

(3.13) $\displaystyle|c\left({\bm{Q}}_{h}{\bm{u}},{\bm{v}}\right)|$
$\displaystyle=\left|\rho\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\left\langle\left({\bm{Q}}_{0}{\bm{u}}-{\bm{Q}}_{b}{\bm{u}}\right)\cdot{\bf
n},\left({\bm{v}}_{0}-{\bm{v}}_{b}\right)\cdot{\bf n}\right\rangle_{\partial
K}\right|$
$\displaystyle\leq\rho\left(\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\left\|{\bm{Q}}_{0}{\bm{u}}-{\bm{Q}}_{b}{\bm{u}}\right\|_{0,\partial
K}^{2}\right)^{1/2}\left(\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\left\|({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n}\right\|_{0,\partial K}^{2}\right)^{1/2}$
$\displaystyle\leq\rho\left(\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\left\|{\bm{Q}}_{0}{\bm{u}}-{\bm{u}}\right\|_{0,\partial
K}^{2}\right)^{1/2}\left(\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\left\|\left({\bm{v}}_{0}-{\bm{v}}_{b}\right)\cdot{\bf
n}\right\|_{0,\partial K}^{2}\right)^{1/2}$ $\displaystyle\lesssim
h^{r}\|{\bm{u}}\|_{r+1,\Omega}\|{\bm{v}}\|_{{\bm{V}}_{h}}.$

Similarly, the upper bound for $l_{\mathrm{div}}({\bm{v}})$ follows directly
from Lemma 2.4 and the approximation property of $\mathbb{Q}_{h}$, that is,

$\displaystyle|l_{\mathrm{div}}({\bm{v}})|$
$\displaystyle=\left|\sum_{K\in\mathcal{T}_{h}}\langle({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n},p-\mathbb{Q}_{h}p\rangle_{\partial K}\right|$
$\displaystyle\lesssim\left(\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\left\|({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n}\right\|_{0,\partial K}^{2}\right)^{1/2}h^{t+1}\|p\|_{t+1,\Omega_{h}}$
$\displaystyle\lesssim h^{t+1}\|p\|_{t+1,\Omega}\|{\bm{v}}\|_{{\bm{V}}_{h}}.$

This completes the proof of the lemma. ∎

######  Lemma 3.9.

Assume that ${\bm{u}}\in{\bm{V}}\cap[H^{2}(\Omega)]^{2}$ is the solution to
problem (2.1)-(2.2), then

(3.14) $\displaystyle|l_{b}(q)|$ $\displaystyle\lesssim
h^{-\frac{3}{2}}s^{2}\|{\bm{u}}\|_{2,\Omega}\|q\|_{0,\Omega_{h}}$
$\displaystyle\forall\,q\in\Psi_{h},$ $\displaystyle|l_{s}({\bm{v}})|$
$\displaystyle\lesssim
h^{-\frac{1}{2}}s\|{\bm{u}}\|_{2,\Omega}\|{\bm{v}}\|_{{\bm{V}}_{h}}$
$\displaystyle\forall\,{\bm{v}}\in{\bm{V}}_{h}.$

Figure 2. $M_{e}$ in the local coordinate system.

###### Proof.

For each $e\in\mathcal{E}_{h}^{B}$, the crescent-shaped region $M_{e}$ is
surrounded by the straight edge $e$ and the curved edge $\tilde{e}$. Let ${\bf
n}_{e}$ and ${\bf n}_{\tilde{e}}$ be unit normal vectors on $e$ and
$\tilde{e}$, respectively, with directions shown in Figure 2. It is just for
convenience that we draw both ${\bf n}_{e}$ and ${\bf n}_{\tilde{e}}$ as
downward pointing, and the direction has no affect to the final result.

As mentioned in the beginning of the proof of Lemma 3.5, a function
$q\in\Psi_{h}$ can be naturally extended to a function defined on $\Omega$,
which is still denoted by $q$ for simplicity. Using Lemma 2.6 and the inverse
inequality we get

$\sum_{e\in\mathcal{E}_{h}^{B}}\|q\|_{0,M_{e}}^{2}\lesssim
h^{-1}s^{2}\|q\|_{0,\Omega_{h}}^{2},\qquad\sum_{e\in\mathcal{E}_{h}^{B}}\|\nabla
q\|_{0,M_{e}}^{2}\lesssim h^{-3}s^{2}\|q\|_{0,\Omega_{h}}^{2}.$

Combining the above and noticing that the solution ${\bm{u}}$ to problem (2.1)
satisfies ${\bm{u}}\cdot{\bf n}_{\tilde{e}}|_{\tilde{e}}=0$, we have

$\displaystyle|l_{b}(q)|=\left|-\sum_{e\in\mathcal{E}_{h}^{B}}\langle{\bm{u}}\cdot{\bf
n}_{e},q\rangle_{e}\right|$
$\displaystyle=\left|-\sum_{e\in\mathcal{E}_{h}^{B}}\bigg{(}\langle{\bm{u}}\cdot{\bf
n}_{e},q\rangle_{e}-\langle{\bm{u}}\cdot{\bf
n}_{\tilde{e}},q\rangle_{\tilde{e}}\bigg{)}\right|$
$\displaystyle=\left|-\sum_{e\in\mathcal{E}_{h}^{B}}\bigg{(}(\nabla\cdot{\bm{u}},q)_{M_{e}}+({\bm{u}},\nabla
q)_{M_{e}}\bigg{)}\right|$
$\displaystyle\leq\|\nabla\cdot{\bm{u}}\|_{0,\Omega\setminus\Omega_{h}}\|q\|_{0,\Omega\setminus\Omega_{h}}+\|{\bm{u}}\|_{0,\Omega\setminus\Omega_{h}}\|\nabla
q\|_{0,\Omega\setminus\Omega_{h}}$ $\displaystyle\lesssim
s\|{\bm{u}}\|_{2,\Omega}\,h^{-\frac{1}{2}}s\|q\|_{0,\Omega_{h}}+s\|{\bm{u}}\|_{1,\Omega}\,h^{-\frac{3}{2}}s\|q\|_{0,\Omega_{h}}$
$\displaystyle\lesssim
h^{-\frac{3}{2}}s^{2}\|{\bm{u}}\|_{2,\Omega}\|q\|_{0,\Omega_{h}},$

where in the second last step we have used Lemma 2.7. This completes the proof
of the first inequality in (3.14).

For any ${\bm{v}}=\\{{\bm{v}}_{0},{\bm{v}}_{b}\\}\in{\bm{V}}_{h}$, note that
$(({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf n}_{e})|_{e}$ is a polynomial of
$\hat{x}$ in the local coordinate system $\hat{x}$-$\hat{y}$ as shown in
Figure 2. This polynomial form can be viewed as the restriction of a two-
variable polynomial $(({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf n}_{e})^{*}$, which
depends only on $\hat{x}$, on edge $e$. Then

$\displaystyle\sum_{e\in\mathcal{E}_{h}^{B}}\|(({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n}_{e})^{*}\|_{0,M_{e}}^{2}$
$\displaystyle=\sum_{e\in\mathcal{E}_{h}^{B}}\int_{0}^{h_{e}}\int_{0}^{\gamma(\hat{x})}|(({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n}_{e})^{*}|^{2}\,\mathrm{d}\hat{y}\,\mathrm{d}\hat{x}$
$\displaystyle\lesssim\sum_{e\in\mathcal{E}_{h}^{B}}\int_{0}^{h_{e}}s^{2}|({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n}_{e}|^{2}\,\mathrm{d}\hat{x}$ $\displaystyle\lesssim
hs^{2}\|{\bm{v}}\|_{{\bm{V}}_{h}}^{2}.$

Similarly,

$\displaystyle\sum_{e\in\mathcal{E}_{h}^{B}}\|\nabla(({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n}_{e})^{*}\|_{0,M_{e}}^{2}$
$\displaystyle=\sum_{e\in\mathcal{E}_{h}^{B}}\int_{0}^{h_{e}}\int_{0}^{\gamma(\hat{x})}|\frac{\partial(({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n}_{e})^{*}}{\partial\hat{x}}|^{2}\,\mathrm{d}\hat{y}\,\mathrm{d}\hat{x}$
$\displaystyle\lesssim\sum_{e\in\mathcal{E}_{h}^{B}}s^{2}|({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n}_{e}|^{2}_{1,e}$
$\displaystyle\lesssim\sum_{e\in\mathcal{E}_{h}^{B}}\|({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n}_{e}\|^{2}_{0,e}$ $\displaystyle\lesssim h\|{\bm{v}}\|_{{\bm{V}}_{h}}^{2}.$

Combining the above and using the fact that ${\bm{u}}\cdot{{\bf
n}}_{\tilde{e}}|_{\tilde{e}}=0$ together with Lemma 2.7, we have

$\displaystyle|l_{s}({\bm{v}})|$
$\displaystyle=\left|\rho\sum_{K\in\mathcal{T}_{h}^{B}}h_{K}^{-1}\langle{\bm{u}}\cdot{\bf
n},({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf n}\rangle_{\partial
K\cap\partial\Omega_{h}}\right|$ $\displaystyle\lesssim
h^{-1}{\sum_{e\in\mathcal{E}_{h}^{B}}\bigg{|}\langle{\bm{u}}\cdot{{\bf
n}}_{e},({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n}_{e}\rangle_{e}-\langle{\bm{u}}\cdot{{\bf
n}}_{\tilde{e}},(({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n}_{e})^{*}\rangle_{\tilde{e}}}\bigg{|}$
$\displaystyle=h^{-1}{\sum_{e\in\mathcal{E}_{h}^{B}}\bigg{|}(\nabla\cdot{{\bm{u}}},(({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n}_{e})^{*})_{M_{e}}+({{\bm{u}}},\nabla(({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n}_{e})^{*})_{M_{e}}}\bigg{|}$ $\displaystyle\leq
h^{-1}\bigg{(}s\|{\bm{u}}\|_{2,\Omega}\,h^{\frac{1}{2}}s\|{\bm{v}}\|_{{\bm{V}}_{h}}+s\|{\bm{u}}\|_{1,\Omega}\,h^{\frac{1}{2}}\|{\bm{v}}\|_{{\bm{V}}_{h}}\bigg{)}$
$\displaystyle\lesssim
h^{-\frac{1}{2}}s\|{\bm{u}}\|_{2,\Omega}\|{\bm{v}}\|_{{\bm{V}}_{h}}.$

This completes the proof of the lemma. ∎

Now, we are able to present the error estimate:

######  Theorem 3.10.

Let ${\bm{u}}$ and $p$ satisfy the conditions in Lemma 3.8 and Lemma 3.9. The
error ${\bm{e}}_{{\bm{u}}}$ and $e_{p}$ satisfy

$\left\|{\bm{e}}_{{{\bm{u}}}}\right\|_{{\bm{V}}_{h}}+\left\|e_{p}\right\|_{0,\Omega_{h}}\lesssim
h^{r}\|{\bm{u}}\|_{r+1,\Omega}+h^{t+1}\|p\|_{t+1,\Omega}+h^{-\frac{1}{2}}s\|{\bm{u}}\|_{2,\Omega},$

where $0\leq r\leq\alpha$ and $0\leq t\leq\sigma$.

###### Proof.

The theorem follows immediately from the mixed finite element theory,
Inequality (3.10), the error equation (3.12), lemmas 3.8-3.9, and the fact
that $s\leq h$. ∎

The error bound in Theorem 3.10 consists of two parts: an approximation error
$h^{r}\|{\bm{u}}\|_{r+1,\Omega}+h^{t+1}\|p\|_{t+1,\Omega}$ and a consistency
error $h^{-\frac{1}{2}}s\|{\bm{u}}\|_{2,\Omega}$. The consistency error comes
purely from the discrepancy between $\Omega$ an $\Omega_{h}$. It vanishes if
$\Omega=\Omega_{h}$. Here we are only interested in the case when $\Omega$ has
curved boundary and hence $\Omega\neq\Omega_{h}$. To fully examine the effect
of the consistency error, we set $\alpha=\beta\geq 1$, $\sigma=\alpha-1$ and
assume that the exact solution satisfies
${\bm{u}}\in[H^{\alpha+1}(\Omega)]^{2}$, $p\in H^{\alpha}(\Omega)$, so that
the approximation error reaches its optimal order $O(h^{\alpha})$. Then by
adjusting $s$, we get the following results:

######  Corollary 3.11.

Assuming the polynomial approximation error reaches the optimal
$O(h^{\alpha})$ with $\alpha\geq 1$. On domains with curved boundary, when
$s=O(h)$, the errors $\bm{e}_{{\bm{u}}}$ and $e_{p}$ satisfy

$\left\|{\bm{e}}_{{{\bm{u}}}}\right\|_{{\bm{V}}_{h}}+\left\|e_{p}\right\|_{0,\Omega_{h}}\lesssim
O(h^{\frac{1}{2}}).$

Numerical results in Section 5 will show that the result in Corollary 3.11 is
indeed sharp.

######  Corollary 3.12.

Assuming that the polynomial approximation error reaches the optimal
$O(h^{\alpha})$ with $\alpha\geq 1$. On domains with curved boundary, when
$s=O(h^{\alpha+\frac{1}{2}})$, the errors $\bm{e}_{{\bm{u}}}$ and $e_{p}$
reach an optimal

$\left\|{\bm{e}}_{{{\bm{u}}}}\right\|_{{\bm{V}}_{h}}+\left\|e_{p}\right\|_{0,\Omega_{h}}\lesssim
O(h^{\alpha}).$

When $h$ is small, setting $s=O(h^{\alpha+\frac{1}{2}})$ means to approximate
the curved boundary $\partial\Omega$ with multiple short edges. Numerical
results in Section 5 will show that for moderate $\alpha$ (our experiments
take $\alpha=1,\,2,\,3$), this approach seem to work quite well.

##  4\. A modified weak Galerkin method

In this section, we construct a modified weak Galerkin discretization for
problem (2.1)-(2.2), for which the consistency error improves from
$O(h^{-\frac{1}{2}}s)$ to $O(h^{-\frac{1}{2}}s^{2})$.

Previously we have defined ${\bf n}$ and $\widetilde{\bf n}$ as the outward
normal on $e\in\mathcal{E}_{h}^{B}$ and $\tilde{e}$, respectively. By using
the map $\bm{\Gamma}$ as shown in Figure 1, the definition of $\widetilde{\bf
n}$ can be pulled to $e$ through $\widetilde{\bf n}\circ\bm{\Gamma}$. We still
denote this pullback by $\widetilde{\bf n}$, i.e., $\widetilde{\bf n}$ is now
also well-defined on $e\in\mathcal{E}_{h}^{B}$. For simplicity of notation, on
interior edges $e\in\mathcal{E}_{h}^{I}$ we just set $\widetilde{\bf n}={\bf
n}$. Now we introduce the following modified bilinear forms

(4.1) $\displaystyle c_{1}({\bm{u}},{\bm{v}})$
$\displaystyle=\rho\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\left\langle\left({\bm{u}}_{0}-{\bm{u}}_{b}\right)\cdot\widetilde{\bf
n},\left({\bm{v}}_{0}-{\bm{v}}_{b}\right)\cdot\widetilde{\bf
n}\right\rangle_{\partial K},$ $\displaystyle a_{h,1}({\bm{u}},{\bm{v}})$
$\displaystyle=\left({\bm{u}}_{0},{\bm{v}}_{0}\right)_{\Omega_{h}}+c_{1}({\bm{u}},{\bm{v}}),$

while the definition $b_{h}({\bm{v}},q)=-(\nabla_{w}\cdot{\bm{v}},q)$ remains
unchanged.

######  Remark 4.1.

The definition of $b_{h}(\cdot,\cdot)$ depends on $\nabla_{w}\cdot$ which in
turn depends on the choice of outward normal vectors on each
$K\in\mathcal{T}_{h}$. We emphasize that the definition of $\nabla_{w}\cdot$
still uses ${\bf n}$ instead of $\widetilde{\bf n}$. Hence
$b_{h}(\cdot,\cdot)$ remains unchanged.

The modified weak Galerkin formulation for problem (2.1)-(2.2) reads as
follows: Find ${\bm{u}}_{h}\in{\bm{V}}_{h}$ and $p_{h}\in\Psi_{h}$ such that

(4.2)
$\begin{cases}a_{h,1}({\bm{u}}_{h},{\bm{v}})+b_{h}({\bm{v}},p_{h})=0\quad&\forall\,{\bm{v}}\in{\bm{V}}_{h},\\\
b_{h}({\bm{u}}_{h},q)-\langle{\bm{u}}_{0}\cdot{\bf
n}-\overline{{\bm{u}}_{0}\cdot{\bf
n}},q\rangle_{\partial\Omega_{h}}=-(g,q)_{\Omega_{h}}\quad&\forall\,q\in\Psi_{h},\end{cases}$

where $\langle{\bm{u}}_{0}\cdot{\bf n}-\overline{{\bm{u}}_{0}\cdot{\bf
n}},q\rangle_{\partial\Omega_{h}}$ in the second equation is a boundary
correction term with

$\overline{{\bm{u}}_{0}\cdot{\bf
n}}|_{e}=\frac{1}{|e|}\int_{e}{\bm{u}}_{0}\cdot{\bf
n}\,\mathrm{d}s,\quad\forall e\in\mathcal{E}_{h}^{B}.$

######  Remark 4.2.

The term $\overline{{\bm{u}}_{0}\cdot{\bf n}}$ is added for two purposes.
First, it ensures that the left-hand side of the second equation of (4.2)
vanishes for $q=const$. Then, under the discrete inf-sup condition to be
proved later, the kernels of both the stiffness matrix for (4.2) and its
transpose are $q\equiv const$. As discussed earlier in Remark 3.2, one does
not need to worry about the ‘compatibility’ issue for the modified problem
(4.2). Second, it gives an $O(s)$ asymptotic rate to the boundary correction
term, which will play an important role in the current proof of the discrete
inf-sup condition.

Similar to the analysis in Section 3, its not hard to see that

$\|{\bm{v}}\|_{\bm{V}_{h,1}}=\sqrt{a_{h,1}({\bm{v}},{\bm{v}})}$

is a well-defined norm on ${\bm{V}}_{h}$. Moreover, we have

######  Lemma 4.3.

One has $\|{\bm{v}}\|_{\bm{V}_{h,1}}\approx\|{\bm{v}}\|_{{\bm{V}}_{h}}$ for
all ${\bm{v}}\in{\bm{V}}_{h}$.

###### Proof.

By the mesh assumptions A4 and A6, the trace and the inverse inequality
(lemmas 2.4-2.5), the facts that ${\bm{v}}_{b}=\bm{0}$ on $\partial\Omega_{h}$
and each $K\in\mathcal{T}_{h}^{B}$ contains at most $O(\frac{h}{s})$ edges in
$\mathcal{E}_{h}^{B}$, one gets

$\rho\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\left\|\left({\bm{v}}_{0}-{\bm{v}}_{b}\right)\cdot(\widetilde{\bf
n}-{\bf n})\right\|_{0,\partial
K\cap\partial\Omega_{h}}^{2}\lesssim\|{\bm{v}}_{0}\|_{0,\Omega_{h}}^{2}.$

The lemma then follows from the above inequality, the definitions of
$\|\cdot\|_{{\bm{V}}_{h}}$ and $\|\cdot\|_{\bm{V}_{h,1}}$, and the triangle
inequality. ∎

Denote the left-hand side of the second equation in (4.2) by

$b_{h,1}({\bm{u}}_{h},q)=b_{h}({\bm{u}}_{h},q)-\langle{\bm{u}}_{0}\cdot{\bf
n}-\overline{{\bm{u}}_{0}\cdot{\bf n}},q\rangle_{\partial\Omega_{h}}.$

System (4.2) is non-symmetric. Non-symmetric mixed systems have been studied
in [5, 17, 32]. We follow the theoretical framework in [32]. Since
$a_{h,1}(\cdot,\cdot)$ is coercive on the entire ${\bm{V}}_{h}$, stability of
(4.2) only requires that both $b_{h}(\cdot,\cdot)$ and $b_{h,1}(\cdot,\cdot)$
satisfy the discrete inf-sup condition (see [32] and also Remark 4.2.7 in
[9]), which we prove in the following lemma:

######  Lemma 4.4.

When $h$ and $\frac{s}{h}$ are sufficiently small, both $b_{h}(\cdot,\cdot)$
and $b_{h,1}(\cdot,\cdot)$ satisfy the discrete $\inf$-$\sup$ condition under
the modified norm $\|\cdot\|_{\bm{V}_{h,1}}$, i.e.,

(4.3)
$\displaystyle\sup_{{\bm{v}}\in{\bm{V}}_{h}}\frac{|b_{h}({\bm{v}},q)|}{\|{\bm{v}}\|_{\bm{V}_{h,1}}}$
$\displaystyle\gtrsim\|q\|_{0,\Omega_{h}}\qquad\forall\,q\in\Psi_{h},$
$\displaystyle\sup_{{\bm{v}}\in{\bm{V}}_{h}}\frac{|b_{h,1}({\bm{v}},q)|}{\|{\bm{v}}\|_{\bm{V}_{h,1}}}$
$\displaystyle\gtrsim\|q\|_{0,\Omega_{h}}\qquad\forall\,q\in\Psi_{h}.$

###### Proof.

Using lemmas 3.5 and 4.3, one immediately gets the discrete inf-sup condition
for $b_{h}(\cdot,\cdot)$ when $h$ is sufficiently small. As for
$b_{h,1}(\cdot,\cdot)$, through a careful examination of the proof of Lemma
3.5, it is not hard to see that one only needs to prove (using the same
${\bm{w}}$ as in the proof of Lemma 3.5)

(4.4)
$\displaystyle\left(\sum_{K\in\mathcal{T}_{h}^{B}}h_{K}^{-1}\|{\bm{Q}}_{0}\bm{w}\cdot{\bf
n}-\overline{{\bm{Q}}_{0}\bm{w}\cdot{\bf n}}\|_{0,\partial
K\cap\partial\Omega_{h}}^{2}\right)^{1/2}\left(\sum_{K\in\mathcal{T}_{h}^{B}}h_{K}\|q\|_{0,\partial
K\cap\partial\Omega_{h}}^{2}\right)^{1/2}$ $\displaystyle\leq
C(h,s)\|q\|_{0,\Omega_{h}}^{2}$ $\displaystyle<\|q\|_{0,\Omega_{h}}^{2},$

with $C(h,s)$ strictly bounded below $1$. By (3.8), this reduces to proving a
bound for
$\sum_{K\in\mathcal{T}_{h}^{B}}h_{K}^{-1}\|{\bm{Q}}_{0}\bm{w}\cdot{\bf
n}-\overline{{\bm{Q}}_{0}\bm{w}\cdot{\bf n}}\|_{0,\partial
K\cap\partial\Omega_{h}}^{2}$, which can be achieved by using the trace
inequality and the inverse inequality

$\displaystyle\left(\sum_{K\in\mathcal{T}_{h}^{B}}h_{K}^{-1}\|{\bm{Q}}_{0}\bm{w}\cdot{\bf
n}-\overline{{\bm{Q}}_{0}\bm{w}\cdot{\bf n}}\|_{0,\partial
K\cap\partial\Omega_{h}}^{2}\right)^{\frac{1}{2}}$
$\displaystyle\lesssim\left(\sum_{K\in\mathcal{T}_{h}^{B}}h_{K}^{-1}\left(\sum_{e\subset\partial
K\cap\partial\Omega_{h}}s^{2}\|\nabla({\bm{Q}}_{0}\bm{w})\|_{0,e}^{2}\right)\right)^{\frac{1}{2}}$
$\displaystyle\lesssim\left(s^{2}\sum_{K\in\mathcal{T}_{h}^{B}}h_{K}^{-1}\|\nabla({\bm{Q}}_{0}\bm{w})\|_{0,\partial
K\cap\partial\Omega_{h}}^{2}\right)^{\frac{1}{2}}$
$\displaystyle\lesssim\left(s^{2}\sum_{K\in\mathcal{T}_{h}^{B}}h^{-2}\|\nabla({\bm{Q}}_{0}\bm{w})\|_{0,K}^{2}\right)^{\frac{1}{2}}$
$\displaystyle\lesssim\frac{s}{h}\|\bm{w}\|_{1,\Omega},$

Thus the coefficient $C(h,s)$ in (4.4) is of order $O(\frac{s}{h})$. It can be
strictly bounded below $1$ if $\frac{s}{h}$ is sufficiently small.  The rest
of the proof follows directly from the proof of Lemma 3.5. ∎

######  Remark 4.5.

When the curved boundary of $\partial\Omega$ is approximated by multiple short
edges, we usually have $\frac{s}{h}$ being sufficiently small and hence Lemma
4.4 holds. The proof does not work for the case of $s=h$, i.e., when no
refinement is imposed on boundary edges. However, several numerical
experiments to be presented in Section 5 show that the modified discrete
scheme appears to be stable even when $s=h$. We suspect that a better discrete
inf-sup condition, with less restrictions on $s$ and $h$, can be proved using
more advanced skills.

######  Remark 4.6.

The proof of Lemma 4.4 depends on the $O(s)$ approximation provided by
$\overline{{\bm{Q}}_{0}\bm{w}\cdot{\bf n}}$, which is one of the reasons for
introducing this term, as pointed out in Remark 4.2. An interesting question
is whether Lemma 4.4 holds or not with the piecewise average
$\overline{{\bm{Q}}_{0}\bm{w}\cdot{\bf n}}$ replaced by a global average on
$\partial\Omega_{h}$, which also satisfies the compatibility condition
discussed in Remark 4.2. Our numerical experiments (not reported in this
article) suggest that this is possible. But it remains to be proved and the
proof appears to be non-trivial.

According to Theorem 3.1 in [32], the inf-sup conditions (4.3) ensure that the
discrete problem (4.2) admits a unique solution when $h$ and $\frac{s}{h}$ are
sufficiently small. Moreover, the discrete operator in (4.2) is stable in the
sense that the unique solution to

$\begin{cases}a_{h,1}({\bm{u}}_{h},{\bm{v}})+b_{h}({\bm{v}},p_{h})=F({\bm{v}})\quad&\forall\,{\bm{v}}\in{\bm{V}}_{h},\\\
b_{h,1}({\bm{u}}_{h},q)=G(q)\quad&\forall\,q\in\Psi_{h},\end{cases}$

satisfies

(4.5)
$\|{\bm{u}}_{h}\|_{\bm{V}_{h,1}}+\|q\|_{\Psi_{h}}\lesssim\|F\|_{{\bm{V}}_{h}^{\prime}}+\|G\|_{\Psi_{h}^{\prime}}.$

###  4.1. Error analysis

Next we analyze the error of the modified discretization (4.2). Again we first
derive the error equations. An argument similar to the one used in the proof
of Lemma 3.6 shows that

######  Lemma 4.7.

The solution ${\bm{u}}$ and $p$ to problem (2.1) satisfy

$a_{h,1}\left(\widetilde{\bm{Q}}_{h}{\bm{u}},{\bm{v}}\right)+b_{h}\left({\bm{v}},\mathbb{Q}_{h}p\right)=c_{1}({\bm{Q}}_{h}{\bm{u}},{\bm{v}})+l_{s,1}({\bm{v}})-l_{\mathrm{div}}({\bm{v}})\qquad\forall\,{\bm{v}}\in{\bm{V}}_{h},$

where $l_{\mathrm{div}}(\cdot)$ is defined as in Lemma 3.6 and
$l_{s,1}(\cdot)$ is defined by

$l_{s,1}({\bm{v}})=\rho\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\left\langle{\bm{Q}}_{b}{\bm{u}}\cdot\widetilde{\bf
n},\left({\bm{v}}_{0}-{\bm{v}}_{b}\right)\cdot\widetilde{\bf
n}\right\rangle_{\partial K\cap\partial\Omega_{h}}.$

###### Proof.

The proof is similar to the one for Lemma 3.6 given in Appendix A. One simply
replaces (A.3) by

$\displaystyle c_{1}(\widetilde{\bm{Q}}_{h}{\bm{u}},{\bm{v}})$
$\displaystyle=\rho\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\langle({\bm{Q}}_{0}{\bm{u}}-\widetilde{\bm{Q}}_{b}{\bm{u}})\cdot\widetilde{\bf
n},({\bm{v}}_{0}-{\bm{v}}_{b})\cdot\widetilde{\bf n}\rangle_{\partial K},$
$\displaystyle=c_{1}({\bm{Q}}_{h}{\bm{u}},{\bm{v}})+\rho\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\langle{{\bm{Q}}}_{b}{\bm{u}}\cdot\widetilde{\bf
n},({\bm{v}}_{0}-{\bm{v}}_{b})\cdot\widetilde{\bf n}\rangle_{\partial
K\cap\partial\Omega_{h}},$
$\displaystyle=c_{1}({\bm{Q}}_{h}{\bm{u}},{\bm{v}})+l_{s,1}({\bm{v}}).$

∎

Let ${\bm{u}}$, $p$ be the solution to problem (2.1), and ${\bm{u}}_{h}$,
$p_{h}$ be the solution to the modified weak Galerkin formulation (4.2).
Consider the error ${\bm{e}}_{{\bm{u}}}$ and $e_{p}$ defined in (3.11). By
lemmas 4.7 and 3.7, the error to problem (4.2) satisfies

(4.6)
$\begin{cases}a_{h,1}\left({\bm{e}}_{{\bm{u}}},{\bm{v}}\right)+b_{h}\left({\bm{v}},e_{p}\right)=c_{1}({\bm{Q}}_{h}{\bm{u}},{\bm{v}})+l_{s,1}({\bm{v}})-l_{\text{div}}({\bm{v}})\quad&\forall{\bm{v}}\in{\bm{V}}_{h},\\\
b_{h,1}\left({\bm{e}}_{{\bm{u}}},q\right)=l_{b,1}(q)+l_{b,2}(q)\quad&\forall
q\in\Psi_{h},\end{cases}$

where

$\displaystyle l_{b,1}(q)$
$\displaystyle=\langle({\bm{u}}-{\bm{Q}}_{0}{\bm{u}})\cdot{\bf
n}-\overline{({\bm{u}}-{\bm{Q}}_{0}{\bm{u}})\cdot{\bf
n}},q\rangle_{\partial\Omega_{h}},$ $\displaystyle l_{b,2}(q)$
$\displaystyle=\langle\overline{{\bm{u}}\cdot{\bf
n}},q\rangle_{\partial\Omega_{h}}.$

We first derive an upper bound for $c_{1}(\cdot,\cdot)$, $l_{s,1}(\cdot)$ and
$l_{b,1}(\cdot)$.

######  Lemma 4.8.

For any $0\leq r\leq\alpha$ and $0\leq t\leq\sigma$, assume that
${\bm{u}}\in{\bm{V}}\cap[H^{\max\\{r+1,3\\}}(\Omega)]^{2},p\in\Psi\cap
H^{t+1}(\Omega)$ are the solutions to problem (2.1)-(2.2), and
$({\bm{v}},q)\in({\bm{V}}_{h},\Psi_{h})$; then we have

$\displaystyle c_{1}({\bm{Q}}_{h}{\bm{u}},{\bm{v}})$ $\displaystyle\lesssim
h^{r}\|{\bm{u}}\|_{r+1,\Omega}\|{\bm{v}}\|_{\bm{V}_{h,1}},$ $\displaystyle
l_{s,1}({\bm{v}})$
$\displaystyle\lesssim\left({h^{-\frac{1}{2}}s^{2}\|{\bm{u}}\|_{3,\Omega}}+h^{r}\|{\bm{u}}\|_{r+1,\Omega}\right)\|{\bm{v}}\|_{\bm{V}_{h,1}},$
$\displaystyle l_{b,1}(q)$ $\displaystyle\lesssim
h^{r}\|{\bm{u}}\|_{r+1,\Omega}\|q\|_{0,\Omega_{h}}.$

###### Proof.

The bound for $c_{1}(\cdot,\cdot)$ follows from a trivial modification of
(3.13). The bound for $l_{b,1}(\cdot)$ follows from the Cauchy-Schwarz
inequality, the trace inequality, the approximation property of ${\bm{Q}}_{0}$
and Inequality (3.8). We are left with estimating $l_{s,1}(\cdot)$. Note that

$\displaystyle l_{s,1}({\bm{v}})$
$\displaystyle=\rho\sum_{K\in\mathcal{T}_{h}^{B}}h_{K}^{-1}\left\langle{\bm{Q}}_{b}{\bm{u}}\cdot\widetilde{\bf
n},\left({\bm{v}}_{0}-{\bm{v}}_{b}\right)\cdot\widetilde{\bf
n}\right\rangle_{\partial K\cap\partial\Omega_{h}}$
$\displaystyle\lesssim\left(\sum_{K\in\mathcal{T}_{h}^{B}}h^{-1}_{K}\|{\bm{Q}}_{b}{\bm{u}}\cdot\widetilde{\bf
n}\|_{0,\partial
K\cap\partial\Omega_{h}}^{2}\right)^{1/2}\|{\bm{v}}\|_{\bm{V}_{h,1}},$
$\displaystyle\lesssim\left(\sum_{K\in\mathcal{T}_{h}^{B}}h^{-1}_{K}\|{\bm{u}}\cdot\widetilde{\bf
n}\|_{0,\partial
K\cap\partial\Omega_{h}}^{2}+\sum_{K\in\mathcal{T}_{h}^{B}}h^{-1}_{K}\|({\bm{u}}-{\bm{Q}}_{b}{\bm{u}})\cdot\widetilde{\bf
n}\|_{0,\partial
K\cap\partial\Omega_{h}}^{2}\right)^{1/2}\|{\bm{v}}\|_{\bm{V}_{h,1}}$
$\displaystyle\lesssim\left(\sum_{K\in\mathcal{T}_{h}^{B}}h^{-1}_{K}\|{\bm{u}}\cdot\widetilde{\bf
n}\|_{0,\partial
K\cap\partial\Omega_{h}}^{2}+h^{2r}\|{\bm{u}}\|_{r+1}^{2}\right)^{1/2}\|{\bm{v}}\|_{\bm{V}_{h,1}},$

where in the last step we used the same argument as in (3.13), together with
the fact that
$\|{\bm{u}}-{\bm{Q}}_{b}{\bm{u}}\|_{e}\leq\|{\bm{u}}-{\bm{Q}}_{0}{\bm{u}}\|_{e}$
on each edge.

Now we only need to estimate
$\sum_{K\in\mathcal{T}_{h}^{B}}h^{-1}_{K}\|{\bm{u}}\cdot\widetilde{\bf
n}\|_{0,\partial K\cap\partial\Omega_{h}}^{2}$. We use the same argument as in
the proof of Lemma 2.4.4 in [33]. Reference [33] is written in Chinese. For
reader’s convenience, we present the full detail of the estimate below. The
key is to use ${\bm{u}}\cdot\widetilde{\bf n}|_{\partial\Omega}=0$ and the map
$\bm{\Gamma}$ defined in Assumption A4. In the local coordinate system
$\hat{x}$-$\hat{y}$ as shown in Figure 1, the value of
${\bm{u}}\cdot\widetilde{\bf n}$ at a point $(\hat{x},0)\in e$ satisfies

$\displaystyle\left|({\bm{u}}\cdot\widetilde{\bf n})(\hat{x},0)\right|$
$\displaystyle=\left|\left({\bm{u}}(\hat{x},\gamma(\hat{x}))-\int_{0}^{\gamma(\hat{x})}\frac{\partial{\bm{u}}}{\partial\hat{y}}\,\mathrm{d}\hat{y}\right)\cdot\widetilde{\bf
n}(\hat{x},\gamma(\hat{x}))\right|$
$\displaystyle=\left|\left(\int_{0}^{\gamma(\hat{x})}\frac{\partial{\bm{u}}}{\partial\hat{y}}\,\mathrm{d}\hat{y}\right)\cdot\widetilde{\bf
n}(\hat{x},\gamma(\hat{x}))\right|$ $\displaystyle\lesssim
s^{2}\|\nabla{\bm{u}}\|_{L^{\infty}(\Omega)}.$

Therefore, on each $e\in\mathcal{E}_{h}^{B}$, one has

$\|{\bm{u}}\cdot\widetilde{\bf
n}\|_{e}^{2}=\int_{0}^{h_{e}}\left|({\bm{u}}\cdot\widetilde{\bf
n})(\hat{x},0)\right|^{2}\,\mathrm{d}\hat{x}\lesssim
s^{5}\|\nabla{\bm{u}}\|_{L^{\infty}(\Omega)}^{2}.$

Hence

$\sum_{K\in\mathcal{T}_{h}^{B}}h^{-1}_{K}\|{\bm{u}}\cdot\widetilde{\bf
n}\|_{0,\partial K\cap\partial\Omega_{h}}^{2}\lesssim
h^{-1}s^{4}\|\nabla{\bm{u}}\|_{L^{\infty}(\Omega)}^{2}\lesssim
h^{-1}s^{4}\|{\bm{u}}\|^{2}_{3,\Omega}.$

Combining the above, this completes the proof of the lemma. ∎

Then we derive an upper bound for $l_{b,2}(\cdot)$.

######  Lemma 4.9.

For any $0\leq r\leq\alpha$ and $0\leq t\leq\sigma$, assume that ${\bm{u}}$ is
the solution to (2.1)-(2.2), and $({\bm{v}},q)\in({\bm{V}}_{h},\Psi_{h})$;
then we have

$\displaystyle l_{b,2}(q)$ $\displaystyle\lesssim
h^{-\frac{1}{2}}s^{2}\|{\bm{u}}\|_{3,\Omega}\|q\|_{0,\Omega_{h}}.$

###### Proof.

We first derive an upper bound for $\overline{{\bm{u}}\cdot{\bf n}}|_{e}$ on
each $e\in\mathcal{E}_{h}^{B}$ using integration by parts

$\displaystyle\big{|}\,\overline{{\bm{u}}\cdot{\bf n}}|_{e}\,\big{|}$
$\displaystyle=\left|\frac{1}{|e|}\int_{e}{\bm{u}}\cdot{\bf
n}\,\mathrm{d}s\right|=\frac{1}{|e|}\left|\int_{e}{\bm{u}}\cdot{\bf
n}\,\mathrm{d}s-\int_{\tilde{e}}{\bm{u}}\cdot\widetilde{\bf
n}\,\mathrm{d}s\right|$
$\displaystyle=\frac{1}{|e|}\left|\int_{M_{e}}\nabla\cdot{\bm{u}}\,\mathrm{d}x\right|$
$\displaystyle\lesssim s^{-1}|M_{e}|\|{\bm{u}}\|_{1,\infty,\Omega}$
$\displaystyle\lesssim s^{2}\|{\bm{u}}\|_{3,\Omega}.$

Combining the above with the trace and the inverse inequalities, we have the
bound for $l_{b,2}(q)$

$\displaystyle l_{b,2}(q)$ $\displaystyle=\langle\overline{{\bm{u}}\cdot{\bf
n}},q\rangle_{\partial\Omega_{h}}\leq\sup_{e\in\mathcal{E}_{h}^{B}}\big{|}\,\overline{{\bm{u}}\cdot{\bf
n}}|_{e}\,\big{|}\int_{\partial\Omega_{h}}|q|ds$ $\displaystyle\lesssim
s^{2}\|{\bm{u}}\|_{3,\Omega}\|q\|_{0,\partial\Omega_{h}}$
$\displaystyle\lesssim
h^{-\frac{1}{2}}s^{2}\|{\bm{u}}\|_{3,\Omega}\left(\sum_{K\in\mathcal{T}_{h}^{B}}h_{K}\|q\|_{0,\partial
K\cap\partial\Omega_{h}}^{2}\right)^{\frac{1}{2}}$ $\displaystyle\lesssim
h^{-\frac{1}{2}}s^{2}\|{\bm{u}}\|_{3,\Omega}\|q\|_{\Omega_{h}}.$

∎

Now, we are able to derive the error estimate:

######  Theorem 4.10.

Let ${\bm{u}}$ and $p$ satisfy the condition in Lemma 4.8. The error
$\bm{e}_{{\bm{u}}}$ and $e_{p}$ satisfy

$\left\|{\bm{e}}_{{\bm{u}}}\right\|_{\bm{V}_{h,1}}+\left\|e_{p}\right\|_{0,\Omega_{h}}\lesssim
h^{r}\|{\bm{u}}\|_{r+1,\Omega}+h^{t+1}\|p\|_{t+1,\Omega}+h^{-\frac{1}{2}}s^{2}\|{\bm{u}}\|_{3,\Omega}.$

where $0\leq r\leq\alpha$ and $0\leq t\leq\sigma$.

###### Proof.

The estimate follows from the theoretical framework in [32], the stability
result (4.5), the error equation (4.6), lemmas 3.8 and 4.8-4.9. ∎

Similar to Corollary 3.12, one can adjust $s$ to improve the error estimate. A
simple calculation shows that an optimal $O(h^{\alpha})$ convergence requires
$s=O(\sqrt{h^{\alpha+\frac{1}{2}}})$, which is less demanding comparing to the
$s=O(h^{\alpha+\frac{1}{2}})$ in Corollary 3.12. However, we still need to
check the requirements for the discrete inf-sup condition stated in Lemma 4.4.
When $\alpha\geq 2$, by taking $s=O(h^{\alpha+\frac{1}{2}})$ one gets
$\frac{s}{h}\leq O(h^{\frac{1}{4}})$ which goes to $0$ as $h$ goes to $0$,
i.e., $\frac{s}{h}$ is sufficiently small whe $h$ is small. It is different
for the case $\alpha=1$, where
$O(\sqrt{h^{\alpha+\frac{1}{2}}})=O(h^{\frac{3}{4}})>O(h)$. Since we always
have $s\leq h$, optimal convergence holds as long as $\frac{s}{h}$ is
sufficiently small.

######  Corollary 4.11.

Assuming the polynomial approximation error reaches the optimal
$O(h^{\alpha})$ with $\alpha\geq 1$. For $\alpha=1$, when $h$ and
$\frac{s}{h}$ are sufficiently small, the errors $\bm{e}_{{\bm{u}}}$ and
$e_{p}$ reach an optimal

$\left\|{\bm{e}}_{{{\bm{u}}}}\right\|_{\bm{V}_{h,1}}+\left\|e_{p}\right\|_{0,\Omega_{h}}\lesssim
O(h).$

For $\alpha\geq 2$, when $s=O(\sqrt{h^{\alpha+\frac{1}{2}}})$ and $h$ is
sufficiently small, the errors $\bm{e}_{{\bm{u}}}$ and $e_{p}$ reach an
optimal

$\displaystyle\left\|{\bm{e}}_{{{\bm{u}}}}\right\|_{\bm{V}_{h,1}}+\left\|e_{p}\right\|_{0,\Omega_{h}}\lesssim
O(h^{\alpha}).$

######  Remark 4.12.

As mentioned previously in Remark 4.5, we suspect that the condition of
$\frac{s}{h}$ being sufficiently small is not necessary, but are not yet able
to prove it. Indeed, numerical results in Section 5 will show that the
estimate in Corollary 4.11 is sharp in the case $\alpha=1$ with just $s=O(h)$.

##  5\. Numerical examples

In this section, we present a sequence of numerical examples to validate the
accuracy of the original weak Galerkin method (3.3) and the modified weak
Galerkin method (4.2). All experiments are done with polynomial degrees
$\alpha=\beta=j$ and $\sigma=j-1$, for $j\geq 1$, which, according to theorems
3.10 and 4.10, are the best choice to minimize the approximation error.
Throughout this section, we shall address such scheme using
$P_{j}$-$P_{j}$-$P_{j-1}$, with various $j$s.

Example 5.1. Rectangular domain

We first quickly present numerical results on the square domain
$\Omega=(0,1)\times(0,1)$, which shall serve as a comparison group. Uniform
triangular meshes are used. The exact solution is,

${\bm{u}}(x,y)=\left(\begin{aligned} \pi\sin(\pi x)\cos(\pi y)\\\ \pi\cos(\pi
x)\sin(\pi y)\end{aligned}\right),\,\qquad p(x,y)=\cos(\pi x)\cos(\pi y),$

which satisfies a homogeneous Neumann boundary condition ${\bm{u}}\cdot{\bf
n}=0\text{ on }\partial\Omega$ and $\int_{\Omega}p\,\mathrm{d}x=0$. In Figure
3, it can be seen that the $P_{j}$-$P_{j}$-$P_{j-1}$ original WG schemes
exhibits an optimal $(j+1)$th order convergence, as predicted.

Figure 3. Convergence rates (from left to right; top to bottom) for the
$P_{1}$-$P_{1}$-$P_{0}$, $P_{2}$-$P_{2}$-$P_{1}$, $P_{3}$-$P_{3}$-$P_{2}$ and
$P_{4}$-$P_{4}$-$P_{3}$ original WG scheme (3.3) on a rectangular domain.

Example 5.2. Unit disk

The domain is a unit disk $\Omega=\\{(x,y)\,|\,x^{2}+y^{2}<1\\}$, and the
exact solution is

${\bm{u}}(x,y)=\left(\begin{aligned} 3x^{2}&+y^{2}-3\\\
&2xy\end{aligned}\right),\,\qquad p(x,y)=3x-x(x^{2}+y^{2}),$

which satisfies a homogeneous Neumann boundary condition and
$\int_{\Omega}p\,\mathrm{d}x=0$.

Three meshes, as shown in Figure 4, are considered. The leftmost is a plain
triangular mesh with $s=O(h)$. According to corollaries 3.11, 4.11 and Remark
4.12, on this mesh we expect to have $O(h^{\frac{1}{2}})$ convergence for all
original WG schemes, optimal $O(h)$ convergence for the lowest-order modified
WG scheme, and $O(h^{\frac{3}{2}})$ convergence for all high-order modified WG
schemes. The mesh in the middle is a polygonal mesh with each curved side
further divided into $\lceil h^{{\frac{1}{2}-j}}\rceil$ short straight edges,
where $\lceil\cdot\rceil$ means the ceiling. In other words, we have
$\frac{h}{s}=O(h^{{\frac{1}{2}-j}})$ which is equivalent to
$s=O(h^{j+\frac{1}{2}})$. In Figure 4 the illustration is given for $j=2$.
According to Corollary 3.12, this is the requirement for the original WG to
exhibit optimal convergence. The rightmost is a polygonal mesh with each
curved side divided into $\lceil h^{\frac{3-2j}{4}}\rceil$ short straight
edges, or in other words, $s=O(\sqrt{h^{j+\frac{1}{2}}})$. Again the
illustration is given for $j=2$. One can immediately see that it uses less
short edges than the mesh in the middle. According to Corollary 4.11, this
guarantees that the modified WG scheme has optimal convergence.

Figure 4. Triangular and polygonal meshes on the unit disk.

In figures 5-6, convergence rates of the original and modified WG schemes on
the left mesh of Figure 4 are reported. The results agree with the conclusion
in corollaries 3.11 and 4.11. Higher-order schemes with $j>2$ (we tested $j$
up to $4$) behave the same as the case of $j=2$, and hence are omitted to save
space. We notice an interesting fact that $\|e_{p}\|_{0,\Omega_{h}}$ appears
to have a better convergence rate than
$\|{\bm{e}}_{{{\bm{u}}}}\|_{\bm{V}_{h}}$, which awaits further study in the
future.

Figure 5. Convergence rates for the $P_{1}$-$P_{1}$-$P_{0}$ and
$P_{2}$-$P_{2}$-$P_{1}$ original WG scheme (3.3) on the left mesh in Figure 4.

Figure 6. Convergence rates for the $P_{1}$-$P_{1}$-$P_{0}$ and
$P_{2}$-$P_{2}$-$P_{1}$ modified WG scheme (4.2) on the left mesh in Figure 4.

We then test the original and the modified WG scheme, respectively, on the
middle and the right meshes in Figure 4. The results are reported in figures
7-8. Both are optimal, as predicted in corollaries 3.12 and 4.11.

Figure 7. Convergence rates for the $P_{1}$-$P_{1}$-$P_{0}$ and
$P_{2}$-$P_{2}$-$P_{1}$ original WG scheme (3.3) on the middle mesh in Figure
4.

Figure 8. Convergence rates for the $P_{2}$-$P_{2}$-$P_{1}$,
$P_{3}$-$P_{3}$-$P_{2}$ and $P_{4}$-$P_{4}$-$P_{3}$ modified WG scheme (4.2)
on the right mesh in Figure 4.

Example 5.3. Ring

Finally, we consider a ring domain
$\Omega=\\{(r,\theta)\,|\,r\in(\frac{1}{2},1)\text{, }\theta\in[0,2\pi]\\}$.
This domain is non-convex and thus is not covered by the theoretical analysis
in the paper. However, numerical results to be presented next show that the
convergence rates are the same as results in convex domains. We set the exact
solution to be

${\bm{u}}(r,\theta)=\left(\begin{aligned} &-\frac{r\sin 2\theta\,(8r-9)}{2}\\\
9r+9r&(\sin\theta)^{2}-8r^{2}(\sin\theta)^{2}-4r^{2}-6\end{aligned}\right),\qquad
p(r,\theta)=(4r^{3}+6r-9r^{2})\sin\theta,$

which satisfies a homogeneous Neumann boundary condition and
$\int_{\Omega}p\,\mathrm{d}x=0$.

Again, three meshes are considered, as shown in Figure 9. The leftmost is a
triangular mesh with $s=O(h)$. The middle and the right are polygonal meshes
with $s=O(h^{j+\frac{1}{2}})$ and $s=O(\sqrt{h^{j+\frac{1}{2}}})$,
respectively. Numerical results of the original and the modified WG schemes on
these meshes are reported in figures 10-13, which exhibit the same convergence
rates as in convex domains.

Figure 9. Triangular and polygonal meshes on the ring.

Figure 10. Convergence rates for the $P_{1}$-$P_{1}$-$P_{0}$ and
$P_{2}$-$P_{2}$-$P_{1}$ original WG scheme (3.3) on the left mesh in Figure 9.

Figure 11. Convergence rates for the $P_{1}$-$P_{1}$-$P_{0}$ and
$P_{2}$-$P_{2}$-$P_{1}$ modified WG scheme (4.2) on the left mesh in Figure 9.

Figure 12. Convergence rates for the $P_{1}$-$P_{1}$-$P_{0}$ and
$P_{2}$-$P_{2}$-$P_{1}$ original WG scheme (3.3) on the middle mesh in Figure
9.

Figure 13. Convergence rates for the $P_{2}$-$P_{2}$-$P_{1}$,
$P_{3}$-$P_{3}$-$P_{2}$ and $P_{4}$-$P_{4}$-$P_{3}$ modified WG scheme (4.2)
on the right mesh in Figure 9.

##  Appendix A Proof of lemmas 3.6 and 3.7

Proof of Lemma 3.6. Testing the first equation in (2.1) with ${\bm{v}}_{0}$ in
${\bm{v}}\in{\bm{V}}_{h}$ and using integration by parts, one gets

(A.1) $\displaystyle 0=$
$\displaystyle({\bm{u}},{\bm{v}}_{0})_{\Omega_{h}}+(\nabla
p,{\bm{v}}_{0})_{\Omega_{h}}$ $\displaystyle=$
$\displaystyle({\bm{u}},{\bm{v}}_{0})_{\Omega_{h}}+\sum_{K\in\mathcal{T}_{h}}(-(\nabla\cdot{\bm{v}}_{0},p)_{K}+\langle{\bm{v}}_{0}\cdot{\bf
n},p\rangle_{\partial K})$ $\displaystyle=$
$\displaystyle({\bm{u}},{\bm{v}}_{0})_{\Omega_{h}}+\sum_{K\in\mathcal{T}_{h}}(-(\nabla\cdot{\bm{v}}_{0},p)_{K}+\langle({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n},p\rangle_{\partial K}),$

where, in the last step we have used ${\bm{v}}_{b}=\bm{0}$ on
$\partial\Omega_{h},$ and the continuity of $p$ across edges in
$\mathcal{E}_{h}^{I}$. Since $\mathbb{Q}_{h}$ and ${\bm{Q}}_{h}$ are
projections, one gets

(A.2) $\displaystyle b_{h}({\bm{v}},\mathbb{Q}_{h}p)$
$\displaystyle=-(\nabla_{w}\cdot{\bm{v}},\mathbb{Q}_{h}p)$
$\displaystyle=\sum_{K\in\mathcal{T}_{h}}\big{(}({\bm{v}}_{0},\nabla(\mathbb{Q}_{h}p))_{K}-\langle{\bm{v}}_{b}\cdot{\bf
n},\mathbb{Q}_{h}p\rangle_{\partial K}\big{)}$
$\displaystyle=\sum_{K\in\mathcal{T}_{h}}\big{(}-(\nabla\cdot{\bm{v}}_{0},\mathbb{Q}_{h}p)_{K}+\langle({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n},\mathbb{Q}_{h}p\rangle_{\partial K}\big{)}$
$\displaystyle=\sum_{K\in\mathcal{T}_{h}}\big{(}-(\nabla\cdot{\bm{v}}_{0},p)_{K}+\langle({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n},\mathbb{Q}_{h}p\rangle_{\partial K}\big{)}.$

Using the definitions of $\widetilde{\bm{Q}}_{h}$ and ${\bm{Q}}_{h}$, we have

(A.3) $\displaystyle c(\widetilde{\bm{Q}}_{h}{\bm{u}},{\bm{v}})$
$\displaystyle=\rho\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\langle({\bm{Q}}_{0}{\bm{u}}-\widetilde{\bm{Q}}_{b}{\bm{u}})\cdot{\bf
n},({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf n}\rangle_{\partial K},$
$\displaystyle=c({\bm{Q}}_{h}{\bm{u}},{\bm{v}})+\rho\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\langle{{\bm{Q}}}_{b}{\bm{u}}\cdot{\bf
n},({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf n}\rangle_{\partial
K\cap\partial\Omega_{h}},$
$\displaystyle=c({\bm{Q}}_{h}{\bm{u}},{\bm{v}})+\rho\sum_{K\in\mathcal{T}_{h}}h_{K}^{-1}\langle{\bm{u}}\cdot{\bf
n},({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf n}\rangle_{\partial
K\cap\partial\Omega_{h}}$
$\displaystyle=c({\bm{Q}}_{h}{\bm{u}},{\bm{v}})+l_{s}({\bm{v}}).$

Combining (A.1)-(A.3), one gets

(A.4) $\displaystyle
a_{h}(\widetilde{\bm{Q}}_{h}{{\bm{u}}},{\bm{v}})+b_{h}({\bm{v}},\mathbb{Q}_{h}p)$
$\displaystyle=({\bm{Q}}_{0}{{\bm{u}}},{\bm{v}}_{0})_{\Omega_{h}}+c(\widetilde{\bm{Q}}_{h}{{\bm{u}}},{\bm{v}})+b_{h}({\bm{v}},\mathbb{Q}_{h}p)$
$\displaystyle=({{\bm{u}}},{\bm{v}}_{0})_{\Omega_{h}}+c({\bm{Q}}_{h}{{\bm{u}}},{\bm{v}})+l_{s}({\bm{v}})$
$\displaystyle\qquad+\sum_{K\in\mathcal{T}_{h}}(-(\nabla\cdot{\bm{v}}_{0},p)_{K}+\langle({\bm{v}}_{0}-{\bm{v}}_{b})\cdot{\bf
n},\mathbb{Q}_{h}p\rangle_{\partial K})$
$\displaystyle=c({\bm{Q}}_{h}{\bm{u}},{\bm{v}})+l_{s}({\bm{v}})-l_{\mathrm{div}}({\bm{v}}).$

This completes the proof of the lemma. ∎

Proof of Lemma 3.7. By Lemma 3.3, we have

$b_{h}({\bm{Q}}_{h}{{\bm{u}}},q)=-(\nabla_{w}\cdot({\bm{Q}}_{h}{{\bm{u}}}),q)_{\Omega_{h}}=-(\pi_{h}(\nabla\cdot{{\bm{u}}}),q)_{\Omega_{h}}=-(\nabla\cdot{{\bm{u}}},q)_{\Omega_{h}}=-(g,q)_{\Omega_{h}}.$

Using the definitions of $\widetilde{\bm{Q}}_{h}$, ${\bm{Q}}_{h}$ and
$\nabla_{w}\cdot$, one gets

$\displaystyle b_{h}(\widetilde{\bm{Q}}_{h}{{\bm{u}}},q)$
$\displaystyle=b_{h}({\bm{Q}}_{h}{{\bm{u}}},q)-(\nabla_{w}\cdot(\widetilde{\bm{Q}}_{h}{{\bm{u}}}-{\bm{Q}}_{h}{{\bm{u}}}),q)_{\Omega_{h}}$
$\displaystyle=-(g,q)_{\Omega_{h}}+\sum_{K\in\mathcal{T}_{h}^{B}}\langle{\bm{Q}}_{b}{{\bm{u}}}\cdot{{\bf
n}},q\rangle_{\partial K\cap\partial\Omega_{h}}$
$\displaystyle=-(g,q)_{\Omega_{h}}+l_{b}(q),$

where the last step follows from the property of the projection
${\bm{Q}}_{b}$. This completes the proof of the lemma. ∎

## References

  * [1] N. Atallah, C. Canuto and G. Scovazzi.  The second-generation shifted boundary method and its numerical analysis. Comput. Methods Appl. Mech. Engrg., 372 (2020), Article 113341.
  * [2] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. Marini, and A. Russo.  Basic principles of virtual element methods.  Math. Models Meth. Appl. Sci., 23(2013), 199-214.
  * [3] L. Beirão da Veiga, F. Brezzi, L. Marini, and A. Russo.  Polynomial preserving virtual elements with curved edges. Math. Models Meth. Appl. Sci., 30 (2020), 1555-1590.
  * [4] L. Beirão da Veiga, A. Russo, and G. Vacca.  The virtual element method with curved edges. ESAIM: M2AN, 53 (2019), pp.375-404.
  * [5] C. Bernardi, C. Canuto, and Y. Maday.  Generalized Inf-Sup conditions for Chebyshev spectral approximation of the stokes problem.  SIAM J. Numer. Anal., 25 (1988), 1237-1265.
  * [6] S. Bertoluzza, M. Pennacchio, and D. Prada.  High order VEM on curved domains. Atti Accad. Naz. Lincei Rend. Lincei Mat. Appl., 30 (2019), 391-412.
  * [7] F. Bertrand, S. Munzenmaier, and G. Starke.  First-order system least squares on curved boundaries: Lowest-order Raviart-Thomas elements.  SIAM J. Numer. Anal., 52 (2014), 880-894.
  * [8] F. Bertrand, S. Munzenmaier, and G. Starke.  First-order system least squares on curved boundaries: Higher-order Raviart-Thomas elements.  SIAM J. Numer. Anal., 52 (2014), 3165-3180.
  * [9] D. Boffi, F. Brezzi, and M. Fortin.  Mixed finite element methods and applications.  Springer, 2013.
  * [10] J. Bramble, T. Dupont, and V. Thomée.  Projection methods for Dirichlet’s problem in approximating polygonal domains with boundary-value corrections.  Math. Comp., 26 (1972), 869-879.
  * [11] J. Bramble and J. King.  A robust finite element method for nonhomogeneous Dirichlet problems in domains with curved boundaries.  Math. Comp., 63 (1994), 1-17.
  * [12] S. Brenner and L. Scott.  The mathematical theory of finite element methods.  Springer, 2008.
  * [13] E. Burman, P. Hansbo and M. Larson.  A cut finite element method with boundary value correction. Math. Comp., 87 (2018), 633-657.
  * [14] E. Burman, P. Hansbo and M. Larson.  A cut finite element method with boundary value correction for the incompressible Stokes’ equations. In: Radu, F.A., Kumar, K., Berre, I., Nordbotten, J.M., Pop, I.S. (eds.), Numerical Mathematics and Advanced Applications ENUMATH 2017, pp. 183-192. Springer, 2019.
  * [15] W. Chen, F. Wang, and Y. Wang.  Weak Galerkin method for the coupled Darcy-Stokes flow.  IMA J. Numer. Anal., 36 (2016), 897-921.
  * [16] J. Cheung, M. Perego, P. Bochev, and M. Gunzburger.  Optimally accurate higher-order finite element methods for polytopial approximations of domains with smooth boundaries.  Math. Comp., 88 (2019), 2187-2219.
  * [17] P. Ciarlet, J. Huang, and J. Zou.  Some observations on generalized saddle-point problems.  SIAM J. Matrix Anal. Appl., 25 (2003), 224-236.
  * [18] B. Cockburn, J. Gopalakrishnan, and R. Lazarov.  Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems.  SIAM J. Numer. Anal., 47 2009), 1319-1365.
  * [19] B. Cockburn and M. Solano.  Solving Dirichlet boundary-value problems on curved domains by extensions from subdomains.  SIAM J. Sci. Comp., 34 (2012), A497-A519.
  * [20] J. Cottrell, T. Hughes, and Y. Bazilevs.  Isogeometric analysis: toward integration of CAD and FEA.  John Wiley & Sons, 2009.
  * [21] I. Ergatoudis, B. Irons, and O. Zienkiewicz.  Curved, isoparametric,“quadrilateral” elements for finite element analysis.  Int. J. Solids Struct., 4 (1968), 31-42.
  * [22] G. Gassnera, F. Lörchera, C.-D. Munza, and J. Hesthaven.  Polymorphic nodal elements and their application in discontinuous Galerkin methods.  J. Comp. Phys., 228 (2009), 1573-1590.
  * [23] V. Girault and P. Raviart.  Finite element methods for Navier-Stokes equations: theory and algorithms.  Springer, 1986.
  * [24] T. Hughes, J. Cottrell, and Y. Bazilevs.  Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement.  Comput. Methods Appl. Mech. Engrg., 194 (2005), 4135-4195.
  * [25] M. Lenoir.  Optimal isoparametric finite elements and error estimates for domains involving curved boundaries.  SIAM J. Numer. Anal., 23 (1986), 562-580.
  * [26] A. Main and G. Scovazzi.  The shifted boundary method for embedded domain computations. Part I: Poisson and Stokes problems. J. Comp. Phys., 372 (2018), 972-995.
  * [27] A. Main and G. Scovazzi.  The shifted boundary method for embedded domain computations. Part II: Linear advection-diffusion and incompressible Navier-Stokes equations. J. Comp. Phys., 372 (2018), 996-1026.
  * [28] L. Mu.  Weak Galerkin finite element with curved edges. J. Comp. Appl. Math., 381 (2021), Article 113038.
  * [29] L. Mu, J. Wang, and X. Ye.  Weak Galerkin finite element methods for the biharmonic equation on polytopal meshes.  Numer. Methods Part. Diff. Eqs., 30 (2014), 1003-1029.
  * [30] L. Mu, J. Wang, and X. Ye.  Weak Galerkin finite element methods on polytopal meshes.  Inter. J. Numer. Anal. Model., 12 (2015), 31-53.
  * [31] L. Mu, X. Wang, and Y. Wang.  Shape regularity conditions for polygonal/polyhedral meshes, exemplified in a discontinuous Galerkin discretization.  Numer. Methods Part. Diff. Eqs., 31 (2015), 308-325.
  * [32] R. Nicolaides.  Existence, uniqueness and approximation for generalized saddle point problems.  SIAM J. Numer. Anal., 19 (1982), 349-357.
  * [33] D. Pang.  Weak Garlerkin finite element for curved regions.  Master’s Thesis, School of Science, Nanjing Normal University, 2020\.
  * [34] G. Strang and A. Berger.  The change in solution due to change in domain.  In: Partial differential equations, pp. 199-205, 1973.
  * [35] V. Thomée.  Polygonal domain approximation in Dirichlet’s problem.  J. Appl. Math., 11 (1973), 33-44.
  * [36] T. Tian, Q. Zhai, R. Zhang.  A new modified weak Galerkin finite element scheme for solving the stationary Stokes equations. J. Comp. Appl. Math., 329 (2018), 268-279.
  * [37] C. Wang, J. Wang, R. Wang, R. Zhang.  A Locking-Free Weak Galerkin Finite Element Method for Elasticity Problems in the Primal Formulation. J. Comp. Appl. Math., 307 (2016), 346-366.
  * [38] J. Wang and X. Ye.  A weak Galerkin finite element method for second-order elliptic problems.  J. Comp. Appl. Math., 241 (2013), 103-115.
  * [39] J. Wang and X. Ye.  A weak Galerkin mixed finite element method for second order elliptic problems.  Math. Comp., 83 (2014), 2101-2126.
  * [40] J. Wang and X. Ye.  A weak Galerkin finite element method for the Stokes equations.  Adv. Comp. Math., 42 (2016), 155-174.
  * [41] J. Zhang, K. Zhang, J. Li, and X. Wang.  A weak Galerkin finite element method for the Navier-Stokes equations.  Commun. Comput. Phys, 23 (2018), 706-746.