Combining the DPG method with finite elements

TL;DR

This paper introduces a novel discretization scheme combining discontinuous Petrov-Galerkin and finite element methods for complex diffusion-advection-reaction problems, demonstrating its theoretical soundness and practical effectiveness.

Contribution

It presents a new coupled DPG-FEM scheme with a heterogeneous variational formulation, including proof of well-posedness and convergence, for general diffusion-advection-reaction models.

Findings

Proves well-posedness of the coupled scheme

Establishes quasi-optimal convergence rates

Numerical results confirm theoretical predictions

Abstract

We propose and analyze a discretization scheme that combines the discontinuous Petrov-Galerkin and finite element methods. The underlying model problem is of general diffusion-advection-reaction type on bounded domains, with decomposition into two sub-domains. We propose a heterogeneous variational formulation that is of the ultra-weak (Petrov-Galerkin) form with broken test space in one part, and of Bubnov-Galerkin form in the other. A standard discretization with conforming approximation spaces and appropriate test spaces (optimal test functions for the ultra-weak part and standard test functions for the Bubnov-Galerkin part) gives rise to a coupled DPG-FEM scheme. We prove its well-posedness and quasi-optimal convergence. Numerical results confirm expected convergence orders.

Equations

$$\displaystyle Au:={\rm div\,}\Bigl{(}-{\underline{\alpha}}\nabla u+{\boldsymbol{\beta}}u\Bigr{)}+\gamma u$$

$$\bar{\Omega }={\bar{\Omega }}_1\cup {\bar{\Omega }}_2,{\Omega }_1\cap {\Omega }_2=\varnothing ,$$

$${\bar{\Omega }}_1\cap {\bar{\Omega }}_2=\Gamma ,$$

$${\bar{\Omega }}_1\cap \partial \Omega ={\bar{\Gamma }}_1,{\bar{\Omega }}_2\cap \partial \Omega ={\bar{\Gamma }}_2.$$

$$H^{1/2}(\mathcal{S})$$

$$H^{-1/2}(\mathcal{S})$$

$$H_{00}^{1/2}(\mathcal{S})$$

$$H_D^{1/2}(\mathcal{S})$$

$$\Vert \hat{u}{\Vert }_{H^{1/2}(\mathcal{S})}$$

$$\Vert \hat{\sigma }{\Vert }_{H^{-1/2}(\mathcal{S})}$$

$$⟨\hat{u},\mathbf{\tau }\cdot \mathbf{n}{⟩}_{\mathcal{S}}:=\sum _{T\in {\mathcal{T}}_1}⟨\hat{u}{|}_{\partial T},\mathbf{\tau }\cdot {\mathbf{n}}_T{⟩}_{\partial T},⟨\hat{\sigma },v{⟩}_{\mathcal{S}}:=\sum _{T\in {\mathcal{T}}_1}⟨\hat{\sigma }{|}_{\partial T},v{⟩}_{\partial T}.$$

$$\mathbf{\sigma }:=\underline{\alpha }\nabla u-\mathbf{\beta }u\text{on}{\Omega }_1,\hat{u}:={\Pi }_{T\in {\mathcal{T}}_1}u{|}_{\partial T},\hat{\sigma }:={\Pi }_{T\in {\mathcal{T}}_1}(\mathbf{\sigma }\cdot {\mathbf{n}}_T){|}_{\partial T}.$$

$$(u,{\mathrm{div}}_{\mathcal{T}}\mathbf{\tau }+\mathbf{\beta }{\underline{\alpha }}^{-T}\mathbf{\tau }+\gamma v{)}_{{\Omega }_1}+(\mathbf{\sigma },{\nabla }_{\mathcal{T}}v+{\underline{\alpha }}^{-T}\mathbf{\tau }{)}_{{\Omega }_1}-⟨\hat{u},\mathbf{\tau }\cdot \mathbf{n}{⟩}_{\mathcal{S}}-⟨\hat{\sigma },v{⟩}_{\mathcal{S}}=(f,v{)}_{{\Omega }_1}.$$

$$(\underline{\alpha }\nabla u-\mathbf{\beta }u,\nabla w{)}_{{\Omega }_2}+(\gamma u,w{)}_{{\Omega }_2}-⟨{\mathbf{n}}_{{\Omega }_2}\cdot (\underline{\alpha }\nabla u-\mathbf{\beta }u),w{⟩}_{\partial {\Omega }_2}=(f,w{)}_{{\Omega }_2}$$

$$(\mathbf{u},u_2)=(u_1,\mathbf{\sigma },\hat{u},\hat{\sigma },u_2)\in U:=U_1\times H_D^1({\Omega }_2)$$

$$\text{with}{U}_1:=L_2({\Omega }_1)\times {\mathbf{L}}_2({\Omega }_1)\times H_D^{1/2}(\mathcal{S})\times H^{-1/2}(\mathcal{S})$$

$$(u_1,{\mathrm{div}}_{\mathcal{T}}\mathbf{\tau }+\mathbf{\beta }{\underline{\alpha }}^{-T}\mathbf{\tau }+\gamma v{)}_{{\Omega }_1}+(\mathbf{\sigma },{\nabla }_{\mathcal{T}}v+{\underline{\alpha }}^{-T}\mathbf{\tau }{)}_{{\Omega }_1}-⟨\hat{u},\mathbf{\tau }\cdot \mathbf{n}{⟩}_{\mathcal{S}}-⟨\hat{\sigma },v{⟩}_{\mathcal{S}}=(f,v{)}_{{\Omega }_1},$$

$$(\underline{\alpha }\nabla u_2-\mathbf{\beta }{u}_2,\nabla w{)}_{{\Omega }_2}+(\gamma u_2,w{)}_{{\Omega }_2}+⟨\hat{\sigma },w{⟩}_{\Gamma }=(f,w{)}_{{\Omega }_2}$$

$$\text{for any}(\mathbf{v},w)\in V\times H_D^1({\Omega }_2)$$

$$\mathbf{v}=(v,\mathbf{\tau })\text{and}V:=H^1({\mathcal{T}}_1)\times \mathbf{H}(\mathrm{div},{\mathcal{T}}_1).$$

$$b(\mathbf{u},\mathbf{v})$$

$$c_2(u_2,w_2)$$

$$L_1(\mathbf{v})$$

$$d(\mathbf{u},u_2;\mathbf{w},w_2):=⟨\hat{\sigma },w_2{⟩}_{\Gamma }+⟨\hat{\chi },\hat{u}-u_2{⟩}_{\Gamma }+\frac{1}{2}⟨\mathbf{\beta }\cdot {\mathbf{n}}_{{\Omega }_1}(\hat{u}-u_2),\hat{w}+w_2{⟩}_{\Gamma }$$

$$\text{for}(\mathbf{u},u_2),(\mathbf{w},w_2)\in U\text{with}\mathbf{u}=(u_1,\mathbf{\sigma },\hat{u},\hat{\sigma }),\mathbf{w}=(w_1,\mathbf{\chi },\hat{w},\hat{\chi }).$$

$$(\mathbf{u},u_2)=(u_1,\mathbf{\sigma },\hat{u},\hat{\sigma },u_2)\in U:$$

$$b(\mathbf{u},\mathbf{v})$$

$$c_2(u_2,w_2)+d(\mathbf{u},u_2;\mathbf{w},w_2)$$

$$c_1(u_1,w_1)$$

$$\mathbf{u}:=(u_1,\mathbf{\sigma },\hat{u},\hat{\sigma })\in U_1:B\mathbf{u}=L_1.$$

$$⟨\Theta \mathbf{u},\mathbf{v}{⟩}_V=b(\mathbf{u},\mathbf{v})\forall \mathbf{v}\in V.$$

$$a(\mathbf{u},u_2;\mathbf{w},w_2)=L(\mathbf{w},w_2)\forall (\mathbf{w},w_2)\in U$$

$$\text{with}a(\mathbf{u},u_2;\mathbf{w},w_2):=b(\mathbf{u},{\Theta }_{\kappa }\mathbf{w})+c_2(u_2,w_2)+d(\mathbf{u},u_2;\mathbf{w},w_2)$$

$$\text{and}L(\mathbf{w},w_2):=L_1({\Theta }_{\kappa }\mathbf{w})+L_2(w_2).$$

$$a(\mathbf{u},u_2;\mathbf{u},u_2)\gtrsim \Vert (\mathbf{u},u_2){\Vert }_U^2\forall (\mathbf{u},u_2)\in U.$$

$$U^0:=\{(u_1,\mathbf{\sigma },\hat{u},\hat{\sigma };u_2)\in U;\hat{u}{|}_{\Gamma }=u_2{|}_{\Gamma }\}$$

$$d^0(\mathbf{u},w_2):=d(\mathbf{u},u_2;\mathbf{w},w_2)=⟨\hat{\sigma },w_2{⟩}_{\Gamma }\forall (\mathbf{u},u_2)=(u_1,\mathbf{\sigma },\hat{u},\hat{\sigma },u_2),(\mathbf{w},w_2)\in U^0.$$

$$a^0(\mathbf{u},u_2;\mathbf{w},w_2)=L(\mathbf{w},w_2)\forall (\mathbf{w},w_2)\in U^0$$

Full text

Key words: DPG method with optimal test functions, finite element method, domain decomposition, coupling, ultra-weak formulation, diffusion-advection-reaction problem

AMS Subject Classification: 65N30, 35J20

1 Introduction

The discontinuous Petrov-Galerkin method with optimal test functions (DPG method) is an approximation scheme that makes the use of optimal test functions, cf. [1, 5, 7], feasible by considering broken test norms [8]. Optimal test functions are those which maximize discrete inf-sup numbers, and the broken form of test spaces and norms allows for their local calculation or approximation. In this form, the DPG method has been developed by Demkowicz and Gopalakrishnan, see the just cited references [7, 8].

The DPG method has been designed having in mind problems where standard methods suffer from locking phenomena (of small inf-sup numbers) or, otherwise, require specific stabilization techniques. This is particularly the case with singularly perturbed problems where DPG schemes have made some contributions [9, 6, 2, 3, 14]. Nevertheless, in the current form most of the schemes are not cheap to implement. On the one hand, corresponding formulations have several unknowns as is the case with mixed finite elements. On the other hand, the efficient approximation of optimal test functions for singularly perturbed problems is ongoing research. For these reasons, advanced DPG techniques are best used for specific problems whereas finite elements are hard to beat when solving uniformly well-posed problems. Though, it has to be said, that in the latter cases DPG schemes can also be efficient and are competitive in general, cf. the software package developed by Roberts [15].

In this paper we develop a discretization method that combines DPG techniques with standard finite elements. In this way, one can restrict the use of more expensive DPG approximations to regions where they are beneficial. Examples are, e.g., reaction-advection-diffusion problems with small diffusivity in a reduced area, or transmission problems that couple a singularly perturbed problem with an unperturbed problem. In a previous publication [12] we have proposed such a combination with boundary elements to solve transmission problems of the Laplacian in the full space, and studied a singularly perturbed case of reaction diffusion in [11]. In this paper we follow the general framework from [12]. There, the basis is set by a heterogeneous variational formulation consisting of an ultra-weak one in a bounded domain and variational boundary integral equations for the exterior unbounded part. Here, we combine an ultra-weak formulation with a standard variational form. We remark that this approach of combining different variational formulations has been systematically analyzed in [10]. Indeed, it is not essential to use an ultra-weak formulation for the DPG scheme, any well-posed formulation would work. Though, the overall strategy in [10] is to employ DPG techniques throughout whereas we combine different discretization techniques.

Having set our heterogeneous formulation, we proceed to rewrite it by using the so-called trial-to-test operator (which maps the test space to the ansatz space). This is only done for the ultra-weak formulation. The whole system then transforms into one where spaces on the ansatz and test sides are identical. In this way, our heterogeneous variational formulation fits the Lax-Milgram framework just as in [12]. We prove coercivity under the condition that the trial-to-test operator is weighted by a sufficiently large constant. Then, quasi-optimal convergence of a discretized version follows by standard arguments. When proving coercivity we follow steps that are similar to the ones when studying the coupling of DPG with boundary elements. But whereas [12] analyzes only the Laplacian, here we set up the scheme and prove coercivity for a general second-order equation of reaction-advection-diffusion type. Throughout we assume that our problem is uniformly well posed, i.e., we do not study variations for singularly perturbed cases as in [11]. Also note that, since coefficients are variable, transmission problems can be treated the same way by selecting the sub-domains accordingly. One only has to move the possibly non-homogeneous jump data to the right-hand side functional.

The remainder of this paper is as follows. In Section 2 we start by formulating the model problem. A heterogeneous variational formulation is given in §2.1. There, we also state its well-posedness and coercivity (Theorem 1) and briefly mention a simplified case where continuity across the sub-domain interface is incorporated strongly (Corollary 2). The corresponding discrete DPG-FEM scheme is presented in §2.2. Its quasi-optimal convergence is announced in Theorem 3. Most technical details and proofs are given in Section 3. In the last section we report on some numerical experiments.

Furthermore, throughout the paper, suprema are taken over sets excluding the null element, and the notation $A\lesssim B$ is used to say that $A\leq C\cdot B$ with a constant $C>0$ which does not depend on any quantity of interest. Correspondingly, the notation $A\gtrsim B$ is used.

2 Mathematical setting and main results

Let $\Omega\subset\mathbb{R}^{d}$, $d\in\{2,3\}$, be a bounded, simply connected polygonal/polyhedral Lipschitz domain with boundary $\partial\Omega$, and normal vector $\mathbf{n}_{\Omega}$ on $\partial\Omega$ pointing outside of $\Omega$. We consider the following elliptic problem of diffusion-advection-reaction type. Given $f\in L_{2}(\Omega)$ find $u\in H^{1}_{0}(\Omega)$ such that

[TABLE]

Here, $L_{2}(\Omega)$ and $H^{1}_{0}(\Omega)$ denote standard Sobolev spaces, the latter with zero trace on $\partial\Omega$. Furthermore, all coefficients are supposed to be sufficiently regular, with ${\underline{\alpha}}(x)\in\mathbb{R}^{d\times d}$, ${\boldsymbol{\beta}}(x)\in\mathbb{R}^{d}$, $\gamma(x)\in\mathbb{R}$ for $x\in\bar{\Omega}$. We assume that all coefficients are uniformly bounded. Furthermore, we assume that the symmetric part of ${\underline{\alpha}}$ is positive definite and uniformly bounded from below, with minimum eigenvalue larger than or equal to $\alpha_{0}>0$, and that $\frac{1}{2}{\rm div\,}{\boldsymbol{\beta}}+\gamma\geq 0$ in $\Omega$. These conditions imply that the operator $A$ is bounded and coercive on $H^{1}_{0}(\Omega)$.

2.1 Heterogeneous variational formulation

In order to solve (1) by a combination of DPG method and finite elements, we formulate the problem in a heterogeneous way, using different variational forms in different parts of the domain. For ease of illustration, we restrict ourselves to two Lipschitz sub-domains $\Omega_{1}$, $\Omega_{2}$ (again of polygonal/polyhedral form, each with one connected component) with boundaries $\partial\Omega_{1}$, $\partial\Omega_{2}$, as specified in Figure 1. There, also a notation for the boundary pieces is introduced. In particular, $\Gamma$ denotes the interface between the sub-domains. The picture indicates that both sub-domains touch the boundary of $\Omega$ (where the homogeneous Dirichlet condition is imposed), but this is not essential. For instance, one sub-domain, $\Omega_{2}$, can be of annular type so that, in that case, $\partial\Omega\subset\partial\Omega_{2}$ and $\Gamma=\partial\Omega_{1}$. Other combinations can be analyzed without difficulty, also including Neumann conditions. Nevertheless, since our analysis centers around proving coercivity of bilinear forms, we need positivity of the combined advection-reaction term on a sub-domain that does not touch the Dirichlet boundary.

Assumption 1. *For $i=1,2$ there holds:

If $\mathrm{meas}(\Gamma_{i})=0$ then there exists $\beta>0$ such that $\frac{1}{2}{\rm div\,}{\boldsymbol{\beta}}+\gamma\geq\beta$ a.e. in $\Omega_{i}$.*

Standard and broken Sobolev spaces. Essential for the DPG method is to use broken test spaces. Therefore, at this early stage we consider a partitioning $\mathcal{T}_{1}$ of $\Omega_{1}$ into (regular non-intersecting) finite elements $T$ such that $\bar{\Omega}_{1}=\cup\{\bar{}T;\;T\in\mathcal{T}_{1}\}$, and with skeleton $\mathcal{S}:=\{\partial T;\;T\in\mathcal{T}_{1}\}$.

Before describing the variational formulation we introduce the Sobolev spaces we need. For a domain $\omega\subset\Omega$ we use the standard spaces $L_{2}(\omega)$, $H^{1}(\omega)$, $H^{1}_{0}(\omega)$, and $\mathbf{H}(\mathrm{div},\omega)$. The trace operator acting on $H^{1}(\omega)$ will be denoted simply by $(\cdot)|_{\partial\omega}$. Then we define the trace space $H^{1/2}(\partial\omega):=H^{1}(\omega)|_{\partial\omega}$ and its dual space H^{-1/2}(\partial\omega):=\bigl{(}H^{1/2}(\partial\omega)\bigr{)}^{\prime} with canonical norms. The duality pairing on $\partial\omega$ is $\langle\cdot\hskip 1.42262pt,\cdot\rangle_{\partial\omega}$ and extends the $L_{2}(\partial\omega)$ bilinear form. Correspondingly, $(\cdot\hskip 1.42262pt,\cdot)_{\omega}$ is the $L_{2}(\omega)$ bilinear form.

We also need $H^{1}_{D}(\Omega_{i})$ consisting of $H^{1}$-functions with vanishing trace on $\Gamma_{i}$ ($i=1,2$). Vector-valued spaces and functions will be denoted by bold symbols. Connected with $\mathcal{T}_{1}$ we use the product spaces $H^{1}(\mathcal{T}_{1})$ and $\mathbf{H}(\mathrm{div},\mathcal{T}_{1})$ with corresponding product norms.

Now, related with $\mathcal{T}_{1}$ are the skeleton trace spaces

[TABLE]

and

[TABLE]

Here, $\mathbf{n}_{T}$ is the exterior unit normal vector on $\partial T$, and $({\boldsymbol{q}}\cdot\mathbf{n}_{T})|_{\partial T}$ indicates the standard way of defining normal traces of $\mathbf{H}(\mathrm{div},T)$-functions. The notation $\widehat{u}|_{\partial\Omega_{1}}=0$ (resp. $\widehat{u}|_{\Gamma_{1}}=0$) is to be understood in the sense that $\widehat{u}$ is a $\mathcal{T}_{1}$-piecewise trace of an element of $H^{1}_{0}(\Omega_{1})$ (resp. of $H^{1}_{D}(\Omega_{1})$). These trace spaces are equipped with the norms

[TABLE]

and analogously for $H^{1/2}_{00}(\mathcal{S})$ and $H^{1/2}_{D}(\mathcal{S})$. For functions $\widehat{u}\in H^{1/2}(\mathcal{S})$, $\widehat{\sigma}\in H^{-1/2}(\mathcal{S})$ (they are elements of product spaces) and ${\boldsymbol{\tau}}\in\mathbf{H}(\mathrm{div},\mathcal{T}_{1})$, $v\in H^{1}(\mathcal{T}_{1})$ we use the duality pairings

[TABLE]

Heterogeneous formulation in $\Omega_{1}\cup\Omega_{2}$. In $\Omega_{1}$, where the DPG method will be used, we consider an ultra-weak variational formulation. As mentioned before, this is just for illustration as any other formulation of primal, mixed, dual-mixed or strong type can be used and analyzed analogously to our case, cf. [10, Section 2.3].

The ultra-weak formulation requires additional independent unknowns

[TABLE]

Then we test the defining relation of ${\boldsymbol{\sigma}}$ with ${\underline{\alpha}}^{-T}$ and ${\boldsymbol{\tau}}\in\mathbf{H}(\mathrm{div},\mathcal{T}_{1})$, and equation (1) with $v\in H^{1}(\mathcal{T}_{1})$. Integrating by parts element-wise, and substituting the corresponding terms by ${\boldsymbol{\sigma}}$, $\widehat{u}$, and $\widehat{\sigma}$, we obtain

[TABLE]

Here, ${\rm div}_{\mathcal{T}}\,$ and $\nabla_{\mathcal{T}}$ denote the $\mathcal{T}_{1}$-piecewise divergence and gradient operators, respectively.

In $\Omega_{2}$ we use the standard primal formulation

[TABLE]

for $w\in H^{1}_{D}(\Omega_{2})$.

Solving (1) in $\Omega$ is equivalent to solving (in appropriate spaces) (4) and (5) with homogeneous Dirichlet condition on $\partial\Omega$ and transmission conditions on $\Gamma$. These transmission conditions will be imposed in variational form. For the time being, we replace $\mathbf{n}_{\Omega_{2}}\cdot({\underline{\alpha}}\nabla u-{\boldsymbol{\beta}}u)|_{\Gamma}$ by $-\widehat{\sigma}|_{\Gamma}$ in (5). Here, we slightly abuse the notation of $\widehat{\sigma}$ noting that $\langle\widehat{\sigma}\hskip 1.42262pt,v\rangle_{\mathcal{S}}=\langle\widehat{\sigma}\hskip 1.42262pt,v\rangle_{\Gamma}$ for $v\in H^{1}(\Omega_{1})$ with $v|_{\Gamma_{1}}=0$, cf., e.g., [11, Section 2.2].

We formally distinguish between $u_{1}:=u|_{\Omega_{1}}$ and $u_{2}:=u|_{\Omega_{2}}$. Then, our preliminary heterogeneous variational formulation consists in finding

[TABLE]

such that $\widehat{u}|_{\Gamma}=u_{2}|_{\Gamma}$ and

[TABLE]

with

[TABLE]

This formulation can be used to define the combined DPG-FEM discretization, but requires that $\mathcal{T}_{1}$ be compatible across $\Gamma$ with the finite element mesh in $\Omega_{2}$. We therefore replace the continuity constraint $\widehat{u}|_{\Gamma}=u_{2}|_{\Gamma}$ by a variational coupling on $\Gamma$ that is similar to a DG-bilinear form involving jumps and fluxes across element boundaries. To this end we abbreviate

[TABLE]

and define the coupling bilinear form

[TABLE]

The final combined ultra-weak primal formulation of (1) then reads

[TABLE]

We will also need the bilinear form for $\Omega_{1}$ that corresponds to $c_{2}(\cdot,\cdot)$:

[TABLE]

For reference, we explicitly specify the strong form of (8a):

[TABLE]

Following [10] one can show that (8) is equivalent to (1) so that, in particular, (8) has a unique solution. However, since we will use different strategies for solving (8a) and (8b), we need a slightly different representation.

To this end we define the trial-to-test operator $\Theta:\;U_{1}\to V$ by

[TABLE]

Here, $\langle\cdot\hskip 1.42262pt,\cdot\rangle_{V}$ denotes the canonical inner product in $V$. Note that $\Theta=\mathcal{R}^{-1}B$ with Riesz operator $\mathcal{R}:\;V\to V^{\prime}$. Since $B$ is defined on $U_{1}$ without boundary condition along $\Gamma$ it has a non-trivial kernel, and so does $\Theta$. Still, $\Theta:\;U_{1}\to V$ is surjective. Therefore, denoting by $\Theta_{\kappa}:=\kappa\Theta$ the scaled trial-to-test operator (for $\kappa>0$ to be chosen), an equivalent formulation is: For given $\kappa>0$ find $(\boldsymbol{u},u_{2})\in U$ such that

[TABLE]

One of our main results is the following theorem.

Theorem 1**.**

The variational formulation (11) is well posed, and is equivalent to problem (1) in the following sense. If $u\in H^{1}_{0}(\Omega)$ solves (1) then $(\boldsymbol{u},u_{2})=(u_{1},{\boldsymbol{\sigma}},\widehat{u},\widehat{\sigma},u_{2})$, with $u_{i}:=u|_{\Omega_{i}}$ ($i=1,2$) and ${\boldsymbol{\sigma}},\widehat{u},\widehat{\sigma}$ defined by (3), satisfies $(\boldsymbol{u},u_{2})\in U$ and solves (11).

Vice versa, if $(\boldsymbol{u},u_{2})=(u_{1},{\boldsymbol{\sigma}},\widehat{u},\widehat{\sigma},u_{2})\in U$ solves (11) then $u$ defined by $u|_{\Omega_{i}}:=u_{i}$ ($i=1,2$) satisfies $u\in H^{1}_{0}(\Omega)$ and solves (1).

Furthermore, for sufficiently large $\kappa>0$, the bilinear form $a(\cdot,\cdot)$ is $U$-coercive, i.e.,

[TABLE]

Proof.

By the assumptions on $\Omega$, $f$, ${\underline{\alpha}}$, ${\boldsymbol{\beta}}$, and $\gamma$, problem (1) is uniquely solvable. Furthermore, by the derivation of (11), if $u\in H^{1}_{0}(\Omega)$ solves (1) then $(\boldsymbol{u},u_{2})$ as specified in the statement solves (11). This can be seen by integrating by parts and noting that $d(\boldsymbol{u},u_{2};\boldsymbol{w},w_{2})=\langle\widehat{\sigma}\hskip 1.42262pt,w_{2}\rangle_{\Gamma}$ since $\widehat{u}|_{\Gamma}=u_{2}|_{\Gamma}$, cf. (7).

The coercivity of $a(\cdot,\cdot)$ will be shown in Section 3.1 under the assumption that $\kappa>0$ is large enough. It is also straightforward to show that this bilinear form is bounded on $U\times U$, as is the linear functional $L$ on $U$. In that case the Lax-Milgram lemma proves the well-posedness of (11).

Now, since $\kappa$ introduces only a scaling of the test functions $\Theta_{\kappa}\boldsymbol{w}\in V$, the variational formulation (11) is actually independent of $\kappa\not=0$, and so is its well-posedness. ∎

As previously mentioned, the continuity constraint $\widehat{u}|_{\Gamma}=u_{2}|_{\Gamma}$ can also be imposed strongly. In this case the solution space is

[TABLE]

and the coupling bilinear form reduces to

[TABLE]

The variational formulation becomes: For given $\kappa>0$ find $(\boldsymbol{u},u_{2})\in U^{0}$ such that

[TABLE]

Analogously as Theorem 1 one obtains the well-posedness of (13) and coercivity of $a^{0}(\cdot,\cdot)$.

Corollary 2**.**

The variational formulation (13) is well posed, and is equivalent to problem (1) in the following sense. If $u\in H^{1}_{0}(\Omega)$ solves (1) then $(\boldsymbol{u},u_{2})=(u_{1},{\boldsymbol{\sigma}},\widehat{u},\widehat{\sigma},u_{2})$, with $u_{i}:=u|_{\Omega_{i}}$ ($i=1,2$) and ${\boldsymbol{\sigma}},\widehat{u},\widehat{\sigma}$ defined by (3), satisfies $(\boldsymbol{u},u_{2})\in U^{0}$ and solves (13).

Vice versa, if $(\boldsymbol{u},u_{2})=(u_{1},{\boldsymbol{\sigma}},\widehat{u},\widehat{\sigma},u_{2})\in U^{0}$ solves (13) then $u$ defined by $u|_{\Omega_{i}}:=u_{i}$ ($i=1,2$) satisfies $u\in H^{1}_{0}(\Omega)$ and solves (1).

Furthermore, for sufficiently large $\kappa>0$, the bilinear form $a^{0}(\cdot,\cdot)$ is $U^{0}$-coercive, i.e.,

[TABLE]

2.2 Combined DPG-FEM discretization

The coupled DPG-FEM method consists in solving (11) within finite-dimensional subspaces $U_{hp}\subset U$. The indices $h$ and $p$ indicate that this can be piecewise polynomial, conforming spaces of certain polynomial degrees. Specifically, the components of $U_{hp}$ that belong to the unknowns $u_{1},{\boldsymbol{\sigma}},\widehat{u},\widehat{\sigma}$ will be piecewise polynomial with respect to the mesh $\mathcal{T}_{1}$ and its skeleton $\mathcal{S}$. On the other hand, the component of $U_{hp}$ that approximates $u_{2}$ is piecewise polynomial with respect to a mesh $\mathcal{T}_{2}$ in $\Omega_{2}$. In the current form we do not need compatibility of the meshes $\mathcal{T}_{1}$, $\mathcal{T}_{2}$ along $\Gamma$. The discrete scheme then reads: For given $\kappa>0$ find $(\boldsymbol{u}_{hp},u_{2,{hp}})\in U_{hp}$ such that

[TABLE]

Note that this formulation includes the use of optimal test functions for the discretization in $\Omega_{1}$, cf. (8a) and the corresponding terms in (11) with trial-to-test operator $\Theta_{\kappa}$. On the other hand, the part of the problem that belongs to $\Omega_{2}$ is solved by standard finite elements, cf. the corresponding relation (8b).

Our second main result is the following theorem.

Theorem 3**.**

If $\kappa>0$ is sufficiently large then the scheme (15) is uniquely solvable and converges quasi-optimally, i.e.,

[TABLE]

Proof.

The statement is a direct implication of the $U$-coercivity of $a(\cdot,\cdot)$ for large $\kappa$ by Theorem 1, the Lax-Milgram lemma and Cea’s estimate. ∎

Remark 4**.**

We note that also the discrete scheme can be changed to impose strongly the continuity of the approximations of $\widehat{u}$ and $u_{2}$ across $\Gamma$. This only requires compatibility of the meshes $\mathcal{T}_{1}$ and $\mathcal{T}_{2}$ along the interface, conforming subspaces $U_{hp}\subset U^{0}$, and replacing the bilinear form $a(\cdot;\cdot)$ in (15) by the bilinear form $a^{0}(\cdot;\cdot)$, cf. (14). The quasi-optimal error estimate from Theorem 3 then holds analogously.

3 Technical details and proof of coercivity

We start with recalling the $H^{1}_{0}(\Omega)$-coercivity of the full differential operator $A$. This transforms into the following properties of the bilinear forms $c_{2}$, $c_{1}$, cf. (6), (9).

Lemma 5**.**

The bilinear forms $c_{1}(\cdot,\cdot)$ and $c_{2}(\cdot,\cdot)$ satisfy

[TABLE]

for all $u_{i}\in H^{1}_{D}(\Omega_{i})$ ($i=1,2$).

Proof.

Noting that

[TABLE]

there holds for $u_{i}\in H^{1}_{D}(\Omega_{i})$ ($i=1,2$)

[TABLE]

The coercivity property then follows with the positivity of the symmetric part of ${\underline{\alpha}}$ and by using either the Poincaré-Friedrichs inequality and $\frac{1}{2}{\rm div\,}{\boldsymbol{\beta}}+\gamma\geq 0$ in $\Omega_{i}$ (if $\mathrm{meas}(\Gamma_{i})\not=0$) or Assumption 1, i.e., $\frac{1}{2}{\rm div\,}{\boldsymbol{\beta}}+\gamma\geq\beta_{i}>0$ in $\Omega_{i}$ ($i=1,2$). ∎

We continue with some properties of the operator $B$, cf. (10), when restricted to the space incorporating homogeneous Dirichlet boundary conditions on the whole of $\partial\Omega_{1}$, that is,

[TABLE]

Lemma 6**.**

The operator $B:\;U_{1,0}\to V^{\prime}$ is an isomorphism with $\|B\|_{\mathcal{L}(U_{1,0},V^{\prime})}$ and $\|B^{-1}\|_{\mathcal{L}(V^{\prime},U_{1,0})}$ bounded independently of the mesh $\mathcal{T}_{1}$.

Proof.

This is a particular case of the different variational formulations studied in [4, Example 3.7]. More generally, in [4], Carstensen, Demkowicz and Gopalakrishnan proved that “breaking” a continuous variational formulation of a well-posed problem (by introducing broken test spaces) and using canonical trace norms, this does not alter the well-posedness of the formulation. ∎

Let us introduce the trace space $H^{1/2}_{00}(\Gamma):=H^{1}_{0}(\Omega)|_{\Gamma}$ with canonical norm. To simplify the presentation of some technical details we will need the following trace operator,

[TABLE]

The boundedness of this operator is immediate, and is analogous to the case of the Laplacian on a single domain considered in [12, Lemma 4].

Lemma 7**.**

The operator $\mathrm{tr}_{\Gamma}$ is bounded with bound independent of $\mathcal{T}_{1}$.

In the following we identify the kernel of $B$ when acting on the full space $U_{1}$. Let us recall that $A$ is the operator of our problem (1). For given $\varphi\in H^{1/2}_{00}(\Gamma)$ we define its $A$-harmonic extension $\operatorname{\mathcal{E}}(\varphi):=(u_{1},{\boldsymbol{\sigma}},\widehat{u},\widehat{\sigma})\in U_{1}$ by

[TABLE]

Lemma 8**.**

The operator $\operatorname{\mathcal{E}}:\;H^{1/2}_{00}(\Gamma)\to U_{1}$ is bounded with bound independent of $\mathcal{T}_{1}$. Moreover, $\operatorname{\mathcal{E}}$ is a right-inverse of $\mathrm{tr}_{\Gamma}$, and the image of $\operatorname{\mathcal{E}}$ is the kernel of $B$, $\ker B=\operatorname{\mathcal{E}}H^{1/2}_{00}(\Gamma)$.

Proof.

These statements can be proved analogously to the case of the Laplacian, cf. [12, Lemmas 11 and 12]. ∎

Of course, we also need continuity of the bilinear forms $b(\cdot,\cdot)$, $c_{2}(\cdot,\cdot)$ and $d(\cdot,\cdot)$. This is straightforward to show and has already been used in the initial part of the proof of Theorem 1. We only give the statement.

Lemma 9**.**

The bilinear forms $b:\;U_{1}\times V\to\mathbb{R}$, $c_{2}:\;H^{1}(\Omega_{2})\times H^{1}(\Omega_{2})\to\mathbb{R}$, and $d:\;U\times U\to\mathbb{R}$ are all uniformly (in $\mathcal{T}_{1}$) bounded.

3.1 Proof of coercivity, statement (12) in Theorem 1

We are now ready to prove the $U$-coercivity of the bilinear form $a(\cdot,\cdot)$, cf. (11). We adapt the procedure from [12] to our situation.

Let $(\boldsymbol{u},u_{2})=(u_{1},{\boldsymbol{\sigma}},\widehat{u},\widehat{\sigma},u_{2})\in U$ be given. We start with the simple estimate

[TABLE]

By Lemma 8, the $\widehat{u}$-component of $\boldsymbol{u}-\operatorname{\mathcal{E}}\mathrm{tr}_{\Gamma}(\boldsymbol{u})$ has zero trace on $\partial\Omega_{1}$, i.e., $\boldsymbol{u}-\operatorname{\mathcal{E}}\mathrm{tr}_{\Gamma}(\boldsymbol{u})\in U_{1,0},$ cf. (16). Combining Lemmas 6 and 8 this gives

[TABLE]

The last identity is due to the well-known relations of the trial-to-test operator $\Theta$,

[TABLE]

According to Lemma 8, operator $\operatorname{\mathcal{E}}$ is bounded,

[TABLE]

A combination of (18), (19), and (20) then gives

[TABLE]

We continue by considering $\boldsymbol{u}^{e}:=(u_{1}^{e},{\boldsymbol{\sigma}}^{e},\widehat{u}^{e},\widehat{\sigma}^{e}):=\operatorname{\mathcal{E}}\mathrm{tr}_{\Gamma}(\boldsymbol{u})=\operatorname{\mathcal{E}}\widehat{u}|_{\Gamma}$. In particular, there holds $\|\widehat{u}\|_{H^{1/2}_{00}(\Gamma)}\lesssim\|u_{1}^{e}\|_{H^{1}(\Omega_{1})}$. Noting that, cf. (7),

[TABLE]

an application of Lemma 5 gives

[TABLE]

Relation (3.1) can also be written like

[TABLE]

so that the previous bound becomes

[TABLE]

Now, recalling the definitions of $c_{1}(\cdot,\cdot)$ (see (9)) and the extension operator $\operatorname{\mathcal{E}}$ (see (17)), integration by parts yields the expected relation $c_{1}(u_{1}^{e},u_{1}^{e})=\langle\widehat{\sigma}^{e}\hskip 1.42262pt,\widehat{u}\rangle_{\Gamma}$. Therefore, continuing the estimate,

[TABLE]

The last term in (23) can be bounded by duality, the continuity of $H^{-1/2}(\mathcal{S})\ni(\widehat{\sigma}^{e}-\widehat{\sigma})\mapsto(\widehat{\sigma}^{e}-\widehat{\sigma})|_{\Gamma}\in H^{-1/2}(\Gamma)$, and relation (19). This gives

[TABLE]

Combining this bound with (21) and (23), and applying Young’s inequality, we find that

[TABLE]

for a sufficiently large constant $\kappa>0$. This proves the stated coercivity of $a(\cdot,\cdot)$.

4 Numerical experiments

In this section we report on two numerical experiments. In both of them we choose $d=2$ and, starting from a manufactured solution, we compute the right-hand side function $f$. The solution of the second experiment does not satisfy the homogeneous Dirichlet boundary condition. In this case, we use a standard approach and extend the inhomogeneous Dirichlet datum into the domain and then shift the resulting terms to the right-hand side. As discrete trial space we use

[TABLE]

where $\mathcal{T}_{1}$ and $\mathcal{S}$ are a mesh and its skeleton in $\Omega_{1}$ and $\mathcal{T}_{2}$ is a mesh in $\Omega_{2}$. Throughout, we use meshes $\mathcal{T}_{1}$ and $\mathcal{T}_{2}$ which are compatible on the interface $\Gamma$ (although this is not necessary in our analysis). In the definition of $U_{hp}$, $P^{k}(\mathcal{T}_{1})$ denotes the space of $\mathcal{T}_{1}$-piecewise polynomials of degree $k$, $P^{0}(\mathcal{S})$ denotes the space of piecewise constant functions on $\mathcal{S}$, and $S^{1}_{D}(\mathcal{S})\subset H^{1/2}_{D}(\mathcal{S})$ denotes the space of piecewise affine and continuous functions on $\mathcal{S}$ which vanish on $\Gamma_{1}$. The space $S^{1}_{D}(\mathcal{T}_{2})\subset H^{1}_{D}(\Omega_{2})$ is the space of piecewise affine, globally continuous functions on $\mathcal{T}_{2}$ which vanish on $\Gamma_{2}$. The trial-to-test operator $\Theta_{\kappa}=\kappa\mathcal{R}^{-1}B$ with $\mathcal{R}:V\rightarrow V^{\prime}$ being the Riesz operator is approximated using the discrete Riesz operator $\mathcal{R}_{hp}:V_{hp}\rightarrow V_{hp}^{\prime}$ with a finite dimensional space $V_{hp}\subset V$, which we choose to be

[TABLE]

The resulting method is called practical DPG method, and was analyzed in [13]. In the latter work, it was shown that the additional discretization error of using $V_{hp}$ instead of $V$ does not degrade the convergence order. Throughout, we use $\kappa=1$ and do not encounter difficulties with this choice. Note that if $\left(u_{1},{\boldsymbol{\sigma}},\widehat{u},\widehat{\sigma},u_{2}\right)$ denotes the exact solution of (11) and $\left(u_{1,{hp}},{\boldsymbol{\sigma}}_{hp},\widehat{u}_{hp},\widehat{\sigma}_{hp},u_{2,{hp}}\right)\in U_{hp}$ denotes the discrete solution (15), then by definition of the norm $H^{1/2}(\mathcal{S})$ it holds

[TABLE]

where $I_{\mathcal{T}_{1}}\widehat{u}_{hp}\in S^{1}_{D}(\mathcal{T}_{1})$ is the nodal interpolant of $\widehat{u}_{hp}$ with $(I_{\mathcal{T}_{1}}\widehat{u}_{hp})|_{\mathcal{S}}=\widehat{u}_{hp}$. Likewise,

[TABLE]

where $\mathbf{I}_{\mathcal{T}_{1}}\widehat{\sigma}_{hp}\in\mathcal{R}\mathcal{T}_{0}(\mathcal{T}_{1})$ is the lowest-order Raviart-Thomas interpolant of $\widehat{\sigma}_{hp}$, i.e., $(\mathbf{n}_{T}\cdot\mathbf{I}_{\mathcal{T}_{1}}\widehat{\sigma}_{hp})|_{\partial T}=\widehat{\sigma}_{hp}|_{\partial T}$ for any $T\in\mathcal{T}_{1}$. Furthermore, we plot the errors

[TABLE]

as well as the so-called energy error of the DPG part

[TABLE]

cf. (19). In both experiments, we use a sequence of meshes resulting from uniform mesh refinements. The quasi-optimality result of Theorem 3 and well-known approximation results then show that

[TABLE]

Here, $N$ denotes the overall number of degrees of freedom of $U_{hp}$. Hence, $\textrm{err}(\cdot)=\mathcal{O}(N^{-1/2})$ for all of the errors defined above.

4.1 Experiment 1

We choose $\Omega_{1}:=\left(0,1\right)\times\left(0,1\right)$, $\Omega_{2}:=\left(1,2\right)\times\left(0,1\right)$ and use the exact solution

[TABLE]

The remaining parameters in the equation (1) are chosen as ${\underline{\alpha}}=id$, ${\boldsymbol{\beta}}=\left(xy,1\right)^{\top}$, and $\gamma=1-\sin(\pi x)$. In Figure 2 we plot the errors versus the degrees of freedom on a double logarithmic scale. As expected, all the errors behave like $\mathcal{O}(N^{-1/2})$, which is plotted in black without markers. In Figure 3, we plot the error $\widehat{u}-u_{2}$ on the coupling boundary $\Gamma$ for the case with mesh width $1/32$.

4.2 Experiment 2

We choose $\Omega_{1}:=(0.2,0.7)\times(0.2,1.2)$ and $\Omega_{2}:=(0.7,1.2)\times(0.2,1.2)$ and use the exact solution

[TABLE]

The remaining parameters in the equation 1 are chosen as ${\underline{\alpha}}=\varepsilon\cdot id$, ${\boldsymbol{\beta}}=\exp(x)\left(\sin y,\cos y\right)$, $\gamma=0$, and $\varepsilon=0.05$. The exact solution $u$ has a curved layer of moderate width inside $\Omega$, see Figure 4. In Figure 5 we plot the errors versus the degrees of freedom on a double logarithmic scale. Again, as expected, we obtain the convergence order $\mathcal{O}(N^{-1/2})$. In Figure 6, we plot the error $\widehat{u}-u_{2}$ on the coupling boundary $\Gamma$, again for mesh width $1/32$. Note that the layer of the exact solution cuts through $\Gamma$, and this is reflected in Fig. 6.