A new over-penalized weak Galerkin finite element method. II: Elliptic interface problems (Q2033589)

The goal of this article is to introduce and analyse a mesh-based numerical method for computing the solution to discontinuous elliptic interface problems. Interface problems involving non-homogenous boundary conditions and jumping coefficients such as the one considered here can be slightly delicate to solve and analyse -- one reason being the lack of global solution regularity -- although there is now a considerable literature on this topic (see, e.g., [\textit{J. Guzmán} et al., Math. Comput. 85, No. 301, 2071--2098 (2016; Zbl 1342.65218)] or [\textit{X. Claeys} et al., Adv. Comput. Math. 43, No. 5, 1075--1101 (2017; Zbl 1381.65094)] for an approach based on boundary integral equations and the references therein). The specific problem considered by the authors in this paper is the following: Given a polygonal domain \(\Omega \subset \mathbb{R}^2\) that can be partitioned as \(\Omega= \Omega_1 \cup \Omega_2\) with \(\Omega_1, \Omega_2 \subset \mathbb{R}^2\) being Lipschitz and having common interface \(\Gamma:= \overline{\Omega_1} \cap \overline{\Omega_2}\), we seek to the solution to the interface problem: \[ \begin{split} -\nabla \cdot (A \nabla u)&=f \qquad \text{in } \Omega,\\ u&=g \qquad \text{on } \partial \Omega \setminus \Gamma,\\ \llbracket u \rrbracket_{\Gamma} &=\psi \qquad \text{on } \Gamma,\\ \llbracket A \nabla u \cdot \boldsymbol{\eta}\rrbracket_{\Gamma}&=\phi \qquad \text{on } \Gamma. \end{split}\tag{1} \] Here, \(A=A(x) \colon \Omega \rightarrow \mathbb{R}^{2 \times 2}\) is a symmetric matrix-valued function that is assumed positive definite and piecewise constant (w.r.t. \(\Omega_1\) and \(\Omega_2\)). Moreover \(\llbracket u \rrbracket_{\Gamma}\) and \(\llbracket A\nabla u \cdot \eta\rrbracket_{\Gamma}\) denote the usual Dirichlet jump and co-normal derivative jumps of \(u \in H^1(\Omega)\) respectively, i.e., \begin{align*} \llbracket u\rrbracket_{\Gamma}:=& u\vert_{\partial \Omega_1} + u \vert_{\partial \Omega_2},\\ \llbracket A\nabla u \rrbracket_{\Gamma}:=& (A \nabla u)\vert_{\partial \Omega_1}\cdot \boldsymbol{\eta}_{1} + (A \nabla u)\vert_{\partial \Omega_2}\cdot \boldsymbol{\eta}_{2}, \end{align*} where \(\eta_{i}\) denotes the unit outward pointing normal vector on \(\Omega_i\). The numerical method proposed by the authors to solve the transmission problem (1) can be viewed as part of the family of so-called `weak Galerkin' methods. The original `weak Galerkin' method was introduced by [\textit{J. Wang} and \textit{X. Ye}, J. Comput. Appl. Math. 241, 103--115 (2013; Zbl 1261.65121)] as a numerical procedure for solving source problems involving second order elliptic PDEs. The core idea of this method is to replace classical gradients that appear in the discrete formulation of the underlying PDE with so-called `weak discrete gradients'. Indeed, given a polygon \(K \subset \mathbb{R}^2\) and \(v:=(v_0, v_b) \in \big(L^2(K), H^{\frac{1}{2}}(\partial K)\big)\), Ye and Wang define the weak gradient of \(\nabla_{\mathrm{w}} v\) as a linear functional on \(H(\operatorname{div}, K)\) such that \begin{align*} \langle \nabla_{\mathrm{w}} v, q \rangle = -\int_{K} v_0 \nabla \cdot q \, dK + \int_{\partial K} v_b q\cdot \boldsymbol{\eta}\, ds \qquad \forall q \in H(\operatorname{div}, K). \end{align*} Clearly, discrete weak gradients can now be defined in a natural manner by considering finite-dimensional subspaces of \(H(\operatorname{div}, K)\). Equipped with the notion of discrete weak gradients, Ye and Wang proposed a classical mesh-based Galerkin method where the test and trial space are (piecewise polynomial-based) finite-dimensional subspaces of \(\big(L^2(K), H^{\frac{1}{2}}(\partial K)\big)\), and the classical gradients in the underlying bilinear form are replaced with discrete weak gradients. In the current paper, the authors introduce a mesh-based, over-penalised and stabilised weak Galerkin method in the spirit of Ye and Wang, i.e., by using the notion of weak discrete gradients -- adapted of course to the setting of the transmission problem (1). The remainder of the paper is classical numerical analysis. The authors prove that their proposed over-penalised weak Galerkin method is well-posed, i.e., for any meshwidth \(h >0\) there exists a discrete solution \(u_h\) to the over-penalised weak Galerkin equations (the existence argument is due to [\textit{L. Mu} et al., IMA J. Numer. Anal. 35, No. 3, 1228--1255 (2015; Zbl 1323.65116)]). Next, the authors derive convergence rates for the discrete solution \(u_h\) that depend on the mesh-width and maximal polynomial degree of the approximation space. In particular, the authors show that for a shape-regular and body-fitted partition \(\mathcal{T}_h\) of \(\Omega\) (see, [\textit{J. Wang} and \textit{X. Ye}, Math. Comput. 83, No. 289, 2101--2126 (2014; Zbl 1308.65202)]), if the exact solution \(u\) to Equation (1) has local Sobolev regularity of order \(k+1\) on \(\Omega_1\) and \(\Omega_2\), i.e., \(u\vert_{\Omega_1} \in H^{k+1}(\Omega_1)\) and \(u\vert_{ \Omega_2} \in H^{k+1}(\Omega_2)\), then there exists \(C>0\) independent of \(u\) such that for all mesh-widths \(h> 0\) \begin{align*} ||| u -u_h||| \leq C \left(h^k + h^{(\beta_0-1)/2}\right) \left(\Vert u \Vert_{H^{k+1}(\Omega_1)} +\Vert u \Vert_{H^{k+1}(\Omega_2)}\right), \end{align*} where \(u_h\) denotes the discrete solution constructed using piecewise polynomial of maximal order \(K\) in the over-penalised and stabilised weak Galerkin method, \(\beta \geq 1\) is a user-specified parameter that appears in the penalty term, and \(||| \cdot |||\) is a mesh-dependent energy norm. In particular, the authors recover optimal convergence rates if \(\beta \geq 2k+1\). As a final step in the analysis, the authors also derive similar \(L^2\) error estimates for the discrete solution \(u_h\), which are likewise optimal if \(\beta \geq 2k+1\). The paper ends with a series of numerical experiments that support the theoretical results presented in the paper.