An arbitrary order mixed finite element method with boundary value correction for the Darcy flow on curved domains

Many practical problems arising in science and engineering are posed on the curved domain $\Omega$. Classical finite element methods defined on a polygonal approximation domain $\Omega_{h}$ often suffer from an additional geometric error due to the difference between $\Omega$ and $\Omega_{h}$. This leads to a loss of accuracy for higher-order discretizations [29, 30]. Moreover, when the problem is equipped with an essential boundary condition, effective ways to transfer the boundary condition from $\partial\Omega$ to $\partial\Omega_{h}$ must be designed. Different strategies have been proposed to solve this problem, for example, the isoparametric finite element method [20, 22] and the isogeometric analysis [21, 17]. Another strategy is the boundary value correction method [8], which shifts the essential boundary condition from $\partial\Omega$ to $\partial\Omega_{h}$, and results in a modified variational formulation. One advantage of the boundary value correction method is that it avoids using curved mesh elements and thus reduces the complexity of implementation. Recently, there are many works utilizing the boundary correction strategy, including the discontinuous Galerkin method via extensions from subdomains [15, 16], transferring technique based on the path integral [26], the shifted boundary method [24, 1], the cutFEM [12], the boundary-corrected virtual element method [2] and the boundary-corrected weak Galerkin method [23], etc.

For the mixed finite element method (MFEM), the Neumann boundary condition becomes essential. Again, higher-order MFEM suffers from a loss of accuracy on curved domains. In [3, 4, 5], the authors studied a parametric Raviart-Thomas MFEM on curved domains, which is a generalization of the isoparametric method to the MFEM. In [27], a cutFEM method is proposed in the curved domain, which is based on a primal mixed formulation [7]. The authors of [27] show that the cutFEM method reaches suboptimal convergence.

In this paper, we propose a new boundary corrected MFEM based on the primal mixed formulation [7], and prove that it has optimal convergence rate. Similarly to [18], we weakly impose the Neumann boundary condition. Then, a boundary value correction technique is designed to pull the Neumann boundary data from $\partial\Omega$ to $\partial\Omega_{h}$. However, unlike the boundary correction in [8, 1, 24], which considers the Dirichlet boundary condition, the Neumann boundary condition involves the outward normal vector on the boundary, which is different in $\partial\Omega$ and $\partial\Omega_{h}$. This poses extra difficulty in the design and analysis of our scheme. Finally, following [14], we added a term $(\mathrm{div}\,\cdot,\mathrm{div}\,\cdot)$ to the discrete formulation to increase stability.

The paper is organized as follows. In Section 2, we introduce some notation and settings. In Section 3, we describe the model problem and introduce the boundary value correction method. In Section 4, the discrete space and the discrete form are defined and analyzed. In Section 5, we prove the optimal error estimate. In Section 6, numerical results are presented. Finally, we draw our conclusions in Section 7.

Let $\Omega$ be a bounded open set in $\mathbb{R}^{2}$ with Lipschitz continuous and piecewise $C^{2}$ boundary $\Gamma$. Let $\mathcal{T}_{h}$ be a body-fitted triangulation of $\Omega$, i.e., all boundary vertices of $\mathcal{T}_{h}$ lie on $\Gamma$. We assume that $\mathcal{T}_{h}$ is geometrically conforming, shape regular, and quasi-uniform. Denote by $\Omega_{h}$ the polygonal region occupied by the triangular mesh $\mathcal{T}_{h}$, and by $\Gamma_{h}$ the boundary of $\Omega_{h}$. When $\Omega$ has a curved boundary, the polygonal region $\Omega_{h}$ does not coincide with $\Omega$, as shown in Fig. 1. Denote by $\mathcal{E}_{h}^{b}$ the set of all boundary edges in $\mathcal{T}_{h}$ and by $\mathcal{T}_{h}^{b}$ all mesh elements that contain at least one edge in $\mathcal{E}_{h}^{b}$. For each $e\in\mathcal{E}_{h}^{b}$, denote by $\tilde{e}\subset\Gamma$ the curved edge cut by two end points of $e$ (see Fig. 1). We further assume that each $\tilde{e}$ is $C^{2}$ continuous, i.e. $\tilde{e}$ cannot cross the connection point of two $C^{2}$ continuous pieces of $\Gamma$. Note that $\tilde{e}$ may coincide with $e$ when part of $\Gamma$ is flat.

We use $\Omega_{h}^{e}$ to denote a crescent-shaped region surrounded by $e$ and $\tilde{e}$, as shown in Fig. 1. There are three possibilities: (a) $\Omega_{h}^{e}\subset\Omega\backslash\Omega_{h}$; (b) $\Omega_{h}^{e}\subset\Omega_{h}\backslash\Omega$; (c) $e=\tilde{e}$ and $\Omega_{h}^{e}=\emptyset$. If $\Omega$ is convex, we have $\Omega_{h}\subset\Omega$. But when $\Omega$ is not convex, there is no inclusion relationship between $\Omega$ and $\Omega_{h}$.

Throughout the paper, we assume that the mesh $\mathcal{T}_{h}$ is geometrically conforming, shape regular, and quasi-uniform [10]. For each triangle $K\in{\mathcal{T}}_{h}$, denote by $h_{K}$ the diameter of $K$. Let $h=\max_{K\in{\mathcal{T}}_{h}}h_{K}$. Similarly, for each edge $e\in\mathcal{E}_{h}^{b}$, denote by $h_{e}$ the length of $e$. Denote by $P_{k}(K)$ the space of polynomials on $K\in\mathcal{T}_{h}$ with degree less than or equal to $k$. We use the standard notation for Sobolev spaces [10]. Let $H^{m}(S)$, for $m\in\mathbb{R}$ and $S\subset\mathbb{R}^{2}$ be the usual Sobolev space equipped with the norm $\|\cdot\|_{m,S}$ and the seminorm $|\cdot|_{m,S}$. When $m=0$, the space $H^{0}(S)$ coincides with the square integrable space $L^{2}(S)$. Denote by $L_{0}^{2}(S)$ the subspace of $L^{2}(S)$ consisting of functions with mean value zero. The above notation extends to a portion $s\subset\Gamma$ or $s\subset\Gamma_{h}$. For example, $\|\cdot\|_{m,s}$ is the Sobolev norm on $H^{m}(s)$. We also use the notation $(\cdot,\cdot)_{S}$ to indicate the $L^{2}$ inner-product on $S\subset\mathbb{R}^{2}$, and $\langle\cdot,\cdot\rangle_{s}$ to indicate the duality pair on $s\subset\Gamma$ or $s\subset\Gamma_{h}$.

All of the above notations extend to vector-valued functions, following the usual rule of product spaces. We use bold face letters to distinguish the vector-valued function spaces from the scalar-valued ones. For example, $\bm{P}_{k}(K)=[P_{k}(K)]^{2}$ and $\bm{H}^{m}(S)=[H^{m}(S)]^{2}$ are the spaces of vector-valued functions with each component in $P_{k}(K)$ and $H^{m}(S)$, respectively. Define

The norm in $H(\mathrm{div},S)$ is defined by

Throughout the paper, we use $\lesssim$ to denote less than or equal up to a constant, and the analogous notation $\gtrsim$ to denote greater than or equal up to a constant.

Finally, we introduce some well-known inequalities.

In this section, we briefly describe the model problem and the boundary value correction method [8]. Consider the Darcy’s equation with the Neumann boundary condition

with $f\in L^{2}(\Omega)$ and $g_{N}\in H^{-1/2}(\Gamma)$. Here, $\bm{\mathop{}\mathrm{n}}$ denotes the unit outward normal on $\Gamma$. Problem (3.1) admits a unique weak solution $(\bm{u},p)\in H(\mathrm{div},\Omega)\times L_{0}^{2}(\Omega)$ as long as the following compatibility condition holds:

The weak formulation of Problem (3.1) can be written as: Find $(\bm{u},p)\in H(\mathrm{div},\Omega)\times L_{0}^{2}(\Omega)$ satisfying $\bm{u}\cdot\bm{\mathop{}\mathrm{n}}=g_{N}$ on $\Gamma$, such that

where $H_{0}(\mathrm{div},\Omega)=\{\bm{v}\in H(\mathrm{div},\Omega):\bm{v}\cdot\bm{n}=0\textrm{ on }\Gamma\}$. The existence and uniqueness of a weak solution to the mixed problem (3.3) can be found in [10].

We use the Brezzi-Douglas-Marini (BDM) finite element method [11] on triangulation $\mathcal{T}_{h}$ to discretize Equation (3.3). The discretization is defined on $\Omega_{h}$ rather than on $\Omega$. Note that in (3.3), the Neumann boundary condition $\bm{u}\cdot\bm{\mathop{}\mathrm{n}}=g_{N}$ on $\Gamma$ becomes essential. This causes the main difficulty in its finite element discretization when $\Gamma$ is curved and hence not equal to $\Gamma_{h}$. Another difficulty is that the compatibility condition (3.2) holds only on $\Omega$ but not on $\Omega_{h}$. We shall take care of these issues one by one. First, we use a boundary correction technique [8] to transform the boundary data from $\Gamma$ to $\Gamma_{h}$.

Assume that there exists a map $\rho_{h}:\Gamma_{h}\rightarrow\Gamma$ such that

as shown in Fig. 2, where $\bm{\nu}_{h}$ is the unit vector pointing from $\bm{x}_{h}\in\Gamma_{h}$ to $\rho_{h}(\bm{x}_{h})\in\Gamma$, and $\delta_{h}(\bm{x}_{h})=|\rho_{h}(\bm{x}_{h})-\bm{x}_{h}|$. For convenience, we denote $\bm{x}:=\rho_{h}(\bm{x}_{h})$. Define $\widetilde{\bm{\mathop{}\mathrm{n}}}(\bm{x}_{h}):=\bm{\mathop{}\mathrm{n}}\circ\rho_{h}(\bm{x}_{h})$, i.e. the unit outward normal vector on $\bm{x}\in\Gamma$ is pulled to the corresponding $\bm{x}_{h}$. Denote by $\bm{\mathop{}\mathrm{n}}_{h}$ the unit outward normal vector on $\Gamma_{h}$. Since $\Gamma$ is piecewise $C^{2}$ continuous, there exists a map $\rho_{h}$ satisfying [12]

For a sufficiently smooth vector-valued function $\bm{v}$ defined between $\Gamma$ and $\Gamma_{h}$, consider its $m$-th order Taylor expansion along $\bm{\nu_{h}}$

where $\partial_{\bm{\nu}_{h}}^{j}$ is the $j$-th order directional derivative and the remainder term satisfies

For simplicity of notation, denote

Both $T^{m}\bm{v}$ and $T_{1}^{m}\bm{v}$ are functions defined on $\Gamma_{h}$. Denote by $\widetilde{\bm{v}}(\bm{x}_{h}):=\bm{v}\circ\rho_{h}(\bm{x}_{h})$ the pullback of $\bm{v}$ from $\Gamma$ to $\Gamma_{h}$. Then clearly

Assume that (3.5) holds, we have the following two lemmas.

Lemmas 3.1-3.2 can be easily extended to vector-valued functions.

For a given integer $k\geq 1$, define the BDM finite element space on $\mathcal{T}_{h}$ by

and denote $Q_{0h}:=Q_{h}\cap L_{0}^{2}(\Omega_{h})$.

For $K\in\mathcal{T}_{h}$, define the local interpolation $I_{K}:\bm{H}^{s}(K)\rightarrow BDM_{k}(K),s>1/2$, by

where $b_{K}=\lambda_{1}\lambda_{2}\lambda_{3}$ is a bubble function, and $\lambda_{i},\,i=1,2,3$, are the barycentric coordinates associated with $K$. Denote the global interpolation by $I_{h}=\prod_{K\in\mathcal{T}_{h}}I_{K}:H^{s}(\Omega_{h})\to V_{h}$.

Denote by $\Pi_{k-1}^{0,K}$ and $\Pi_{k}^{0,e}$ the $L^{2}$ orthogonal projections onto $P_{k-1}(K)$ and $P_{k}(e)$, respectively. Define the global projection $\Pi_{h}=\left(\prod_{K\in\mathcal{T}_{h}}\Pi_{k-1}^{0,K}\right):L^{2}(\Omega_{h})\to Q_{h}$.

By the lemmas 2.1,2.2,2.4, it is not hard to see that

Finally, we recall the extension property of Sobolev spaces:

In Lemma 4.4, the hidden constant in $\lesssim$ may depend on the shape of $\Omega$. Throughout the paper, we shall only use Lemma 4.4 on $\Omega$ but not on $\Omega_{h}$, in order to ensure that the hidden constant does not depend on $h$. For brevity, we also denote the extension by $v^{E}=Ev$. The extension and its notation can be extended to vector-valued functions.