A cut finite element method for a model of pressure in fractured media

Abstract

We develop a robust cut finite element method for a model of diffusion in fractured media consisting of a bulk domain with embedded cracks. The crack has its own pressure field and can cut through the bulk mesh in a very general fashion. Starting from a common background bulk mesh, that covers the domain, finite element spaces are constructed for the interface and bulk subdomains leading to efficient computations of the coupling terms. The crack pressure field also uses the bulk mesh for its representation. The interface conditions are a generalized form of conditions of Robin type previously considered in the literature which allows the modeling of a range of flow regimes across the fracture. The method is robust in the following way: (1) Stability of the formulation in the full range of parameter choices; and (2) Not sensitive to the location of the interface in the background mesh. We derive an optimal order a priori error estimate and present illustrating numerical examples.

1 Introduction

The numerical modelling of flow in fractured porous media is important both in environmental science and in industrial applications. It is therefore not surprising that it is currently receiving increasing attention from the scientific computing community. Here we are interested in models where the fractures are modelled as embedded surfaces of dimension $d-1$ in a d dimensional bulk domain. Models on this type of geometries of mixed dimension are typically obtained by averaging the flow equations across the width of the fracture and introducing suitable coupling conditions for the modelling of the interaction with the bulk flow. Such reduced models have been derived for instance in [1, 25, 29]. The coupling conditions in these models typically take the form of a Robin type condition. The physical properties of the coupling enters as parameters in this interface condition. The size of these parameters can vary with several orders of magnitude depending on the physical properties of the crack and of the material in the porous matrix. This makes it challenging to derive methods that both are flexible with respect to mesh geometries and robust with respect to coupling conditions. A wide variety of different strategies for the discretisation of fractured porous media flow has been proposed in the literature. One approach is to use a method that allows for nonconforming coupling between the bulk mesh and the fracture mesh [3], or even arbitrary polyhedral elements in the bulk mesh in order to be able to mesh the fractures easily. This latter approach has been developed using discontinuous Galerkin methods [2], virtual element methods [18] and high order hybridised methods [14].

Herein we will consider an unfitted approach, drawing on previous work [4, 11, 13] where flow in fractured porous media was modelled in the situation where the pressure is a globally continuous function. When using unfitted finite element methods, the bulk mesh can be created completely independently of the fractures. Instead the finite element space is modified locally to allow for discontinuities across fractures and interface conditions are typically imposed weakly, or using methods similar to Nitsche’s method. For other recent work using unfitted methods we refer to [27], where a stabilized Lagrange multiplier method is considered for the interface coupling and [15] where a mixed method is considered for the Darcy’s equations both in the bulk and on the surface.

The upshot here, compared to [11] is that the pressure in the crack has its own approximation field, allowing accurate approximation of problems where there is a pressure jump between the bulk and the fracture, and that the interface conditions are imposed in a way allowing for the full range of parameter values in the Robin condition, without loss of stability or order of approximation. We use the variant of the interface modelling considered in [29], that was also recently applied for the numerical modelling in [2]. In these models we may obtain a wide range of parameter values in the interface condition and we therefore develop a method which handle the full range of values and produces approximations with optimal order convergence. The approach is inspired by the work of Stenberg [26] and may be viewed as a version of the Nitsche method that can handle Robin type conditions and which converges to the standard Nitsche method when the Robin parameter tends to infinity. Previous applications of this approach in the context of fitted finite elements include [22, 37].

The finite element spaces are constructed starting from a standard mesh equipped with a finite element space. For each geometric domain (subdomains and interface) we mark all the elements intersected by the domain and then we restrict the finite element space to that set to form a finite element space for each domain. This procedure leads to cut finite elements and we use stabilization to ensure that the resulting form associated with the method is coercive and that the stiffness matrix is well conditioned. The stabilization is of face or ghost penalty type [6, 7, 28], and is added both to the bulk and interface spaces. Previous related work work on cut finite element methods include the interface problem [21]; overlapping meshes [23]; coupled bulk-surface problems [8, 11, 12, 20]; mixed dimensional problems [9], and surface partial differential equations [7, 32]. For a general introduction to cut finite element methods we refer to [4].

The outline of the paper is as follows: In Sect. 2 we introduce the model problem, show an elliptic regularity result which is robust with respect to the critical parameters in the interface condition, and discuss the relation between our formulation of the interface conditions and previous work; in Sect. 3 we formulate the cut finite element method; in Sect. 4 we prove the basic properties of the formulation and in particular an optimal order a priori error estimate which is uniform in the full range of interface parameters; and in Sect. 5 we present numerical results.

2 The model problem

2.1 Governing equations

Let $\varOmega$ be a convex polygonal domain in $\mathbb {R}^d$, $d=2$ or 3, with boundary $\partial \varOmega$ and exterior unit normal n. Let $\varGamma$ be a smooth embedded interface in $\varOmega$ and let $n_\varGamma$ be a unit normal to $\varGamma$, which partitions $\varOmega$ into two subdomains $\varOmega _1$ and $\varOmega _2$ with exterior unit normals $n_1$ and $n_2$. We assume that $\varGamma$ is a closed surface without boundary residing in the interior of $\varOmega$, more precisely we assume that there is $\delta _0>0$ such that the distance between $\varGamma$ and $\partial \varOmega$ is larger than $\delta$. We consider for simplicity the case with homogeneous Dirichlet conditions on $\partial \varOmega$.

The problem takes the form: find $u_i:\varOmega _i \rightarrow \mathbb {R}$ and $u_\varGamma :\varGamma \rightarrow \mathbb {R}$ such that

$$\begin{aligned} -\nabla \cdot A_i \nabla u_i= & {} f_i \qquad \text {in}\; \varOmega _i \end{aligned}$$

(1)

$$\begin{aligned} -\nabla _\varGamma \cdot A_\varGamma \nabla _\varGamma u_\varGamma= & {} f_\varGamma - \llbracket n \cdot A \nabla u \rrbracket \qquad \text {on}\; \varGamma \end{aligned}$$

(2)

$$\begin{aligned} n \cdot A \nabla u + B (u - u_\varGamma )= & {} 0 \qquad \text {on}\;\varGamma \end{aligned}$$

(3)

$$\begin{aligned} u= & {} 0 \qquad \text {on}\;\partial \varOmega \end{aligned}$$

(4)

Here the jump (or sum) of the normal fluxes is defined by

$$\begin{aligned} \llbracket n \cdot A \nabla v \rrbracket = \sum _{i=1}^2 n_i \cdot A_i \nabla v_i \end{aligned}$$

(5)

In the interface condition (3), B is a $2\times 2$ symmetric matrix valued function with eigenvalues $\lambda _i$ such that $\lambda _i \in [c_\lambda ,\infty )$, with $c_\lambda$ a positive constant, which means that B is uniformly positive definite on $\varGamma$,

$$\begin{aligned} c_\lambda \Vert x \Vert ^2_{\mathbb {R}^2} \le x\cdot B \cdot x\qquad \forall x \in \mathbb {R}^2 \end{aligned}$$

(6)

We also used the notation

$$\begin{aligned} n \cdot A \nabla v = \left[ \begin{matrix} n_1 \cdot A_1 \nabla v_1\\ n_2 \cdot A_2 \nabla v_2 \end{matrix} \right] ,\qquad v - v_\varGamma = \left[ \begin{matrix} v_1 - v_\varGamma \\ v_2 - v_\varGamma \end{matrix} \right] \end{aligned}$$

(7)

and thus in component form (3) reads

$$\begin{aligned} \left[ \begin{matrix} n_1 \cdot A_1 \nabla u_1 \\ n_2 \cdot A_2 \nabla u_2 \end{matrix} \right] + B \left[ \begin{matrix} u_1 - u_\varGamma \\ u_2 - u_\varGamma \end{matrix} \right] = \left[ \begin{matrix} 0 \\ 0 \end{matrix} \right] \end{aligned}$$

(8)

The coefficients $A_1$, $A_2$, are smooth uniformly positive definite symmetric $d \times d$ matrices, $A_\varGamma$ is smooth tangential to $\varGamma$ and uniformly positive definite on the tangent space of $\varGamma$, so that

$$\begin{aligned} \sum _{i=1}^2 \Vert \nabla v_i \Vert ^2_{\varOmega _i} + \Vert \nabla _\varGamma v_\varGamma \Vert ^2_\varGamma \lesssim \sum _{i=1}^2 (A_i \nabla v_i, \nabla v_i )_{\varOmega _i} + (A_\varGamma \nabla _\varGamma v_\varGamma , \nabla _\varGamma v_\varGamma )_\varGamma \end{aligned}$$

(9)

where $\lesssim$ denotes less or equal up to a constant. Finally, we assume $f_i\in L_2(\varOmega _i)$ and $f_\varGamma \in L_2(\varGamma )$.

Remark 1

Several generalizations are possible on the external boundary. For instance, we may let the interface intersect the boundary of $\varOmega$. In this case we let $\nu$ denote the unit exterior conormal to $\varGamma \cap \partial \varOmega$, i.e. $\nu$ is tangent to $\varGamma$ and normal to $\partial \varOmega \cap \varGamma$, and we assume that $\nu \cdot n \ge c > 0$ for some constant c so that the interface is transversal to $\partial \varOmega$. We may then enforce the Dirichlet condition $u_\varGamma = g_\varGamma$ on $\partial \varOmega \cap \varGamma$ (see [10]) or some other standard boundary condition.

Remark 2

In practical modeling we may want to take the thickness of the interface inte account. Assuming that the permeability matrix in an interface of thickness t takes the form

$$\begin{aligned} A|_{U_{t/2}(\varGamma )} = A^e_\varGamma + a^e_{\varGamma } n_\varGamma \otimes n_\varGamma \end{aligned}$$

(10)

where $n_\varGamma$ is a unit normal vector field to $\varGamma$, $U_{t/2}(\varGamma )$ is the set of points with distance less than t/2 to $\varGamma$, $v^e$ denotes the extension of a function v on $\varGamma$ that is constant in the normal direction, $A_\varGamma$ is the tangential tangential permeability tensor, and finally $a_{\varGamma ,n}$ is the permeability across the interface. Also assuming that $f = f_\varGamma ^e$ and $u = u^e$ in $U_{t/2}(\varGamma )$, the equation on the interface (2) may be modelled as follows

$$\begin{aligned} -\nabla _\varGamma \cdot t A_\varGamma \nabla _\varGamma u_\varGamma= & {} t f_\varGamma - \llbracket n \cdot A \nabla u \rrbracket \qquad \text {on}~\varGamma \end{aligned}$$

(11)

Note that the last term on the right hand side does not scale with t since it accounts for flow into the crack from the bulk domains.

Remark 3

We comment on how our interface condition (3) relates to the condition in [29] and later reformulated, see [2], in terms of averages and jumps of the bulk fields across the interface. The interface conditions in [29], Eqs. (3.18) and (3.19), take the form

$$\begin{aligned} \xi n_1 \cdot A_1 \nabla v_1 - (1-\xi ) n_2 \cdot A_2 \nabla v_2= & {} \alpha (v_\varGamma - v_1) \end{aligned}$$

(12)

$$\begin{aligned} \xi n_2 \cdot A_2 \nabla v_2 - (1-\xi ) n_1 \cdot A_1 \nabla v_1= & {} \alpha (v_\varGamma - v_2) \end{aligned}$$

(13)

where $\xi$ and $\alpha$ are parameters. The parameter $\alpha$ is related to physical properties of the interface as follows

$$\begin{aligned} \alpha = \frac{2 a_{\varGamma ,n}}{t} \end{aligned}$$

(14)

where $a_{\varGamma ,n}$ is the permeability coefficient across the interface $\varGamma$ and t is the thickness of the interface, see (3.8) in [29]. In matrix form we obtain

$$\begin{aligned} \left[ \begin{matrix} \xi &{} \xi - 1 \\ \xi - 1 &{} \xi \end{matrix} \right] \left[ \begin{matrix} n_1 \cdot A_1 \nabla v_1 \\ n_2 \cdot A_2 \nabla v_2 \end{matrix} \right] + \left[ \begin{matrix} \alpha &{} 0 \\ 0 &{} \alpha \end{matrix} \right] \left[ \begin{matrix} v_1 - v_\varGamma \\ v_2 - v_\varGamma \end{matrix} \right] = 0 \end{aligned}$$

(15)

which leads to

$$\begin{aligned} B = \frac{1}{2\xi - 1} \left[ \begin{matrix} \xi &{} 1 - \xi \\ 1 - \xi &{} \xi \end{matrix} \right] \left[ \begin{matrix} \alpha &{} 0 \\ 0 &{} \alpha \end{matrix} \right] = \frac{\alpha }{2\xi - 1} \left[ \begin{matrix} \xi &{} 1 -\xi \\ 1 - \xi &{} \xi \end{matrix} \right] \end{aligned}$$

(16)

We note that we have the eigen pairs

$$\begin{aligned} B e_1 = \underbrace{\frac{\alpha }{2\xi - 1}}_{\lambda _1} e_1, \qquad B e_2 =\underbrace{\alpha }_{\lambda _2} e_2 \end{aligned}$$

(17)

with the corresponding eigen vectors defined by

$$\begin{aligned} e_1 = \frac{1}{\sqrt{2}} \left[ \begin{matrix} 1\\ 1 \end{matrix} \right] \quad \text{ and } \quad e_2 = \frac{1}{\sqrt{2}} \left[ \begin{matrix} 1\\ -1 \end{matrix} \right] \end{aligned}$$

(18)

and thus B is positive definite for $\xi >1/2$, singular for $\xi =1/2$, and indefinite for $\xi <1/2$. It is therefore natural to consider the case when $\alpha >0$ and $\xi >1/2$. We remark that when $\alpha$ tends to infinity (zero) both eigenvalues tend to infinity (zero) and when $\xi$ tends to 1/2 from above one eigenvalue tends to infinity. It is therefore important to construct a method which is robust in the full range $\lambda _i \in (0,\infty )$ of possible values for the two eigenvalues.

To see the relation to the formulation of the interface conditions in [2] we first note that we have the expansions

$$\begin{aligned}&\left[ \begin{matrix}n_1 \cdot A_1 \nabla v_1\\ n_2 \cdot A_2 \nabla v_2 \end{matrix} \right] = 2^{-1/2} \llbracket n\cdot A \nabla \rrbracket ~ e_1 + 2^{1/2} \langle \langle n\cdot A \nabla v \rangle \rangle ~ e_2 \end{aligned}$$

(19)

$$\begin{aligned}&\left[ \begin{matrix} v_1 - v_\varGamma \\ v_2 - v_\varGamma \end{matrix} \right] = 2^{1/2} (\langle \langle v \rangle \rangle - v_\varGamma ) ~ e_1 + 2^{-1/2} \llbracket v \rrbracket ~ e_2 \end{aligned}$$

(20)

where the jumps and averages of the bulk fields across the the interface are defined by

$$\begin{aligned}&\displaystyle \llbracket n \cdot A \nabla v \rrbracket = \sum _{i=1}^2 n_i \cdot A_i \nabla v_i, \qquad \llbracket v \rrbracket = v_1 - v_2 \end{aligned}$$

(21)

$$\begin{aligned}&\displaystyle \langle \langle n \cdot A \nabla v \rangle \rangle = \frac{1}{2}(n_1 \cdot A_1 \nabla v_1 - n_2 \cdot A_2 \nabla v_2), \qquad \langle \langle v \rangle \rangle = \frac{1}{2}( v_1 + v_2 ) \end{aligned}$$

(22)

Using the expansions (19) and (20) together with (17) and matching the coefficients associated with each eigenvector we obtain the interface conditions

$$\begin{aligned}&\displaystyle \llbracket n \cdot A\nabla v \rrbracket + \frac{2 \alpha }{2\xi - 1 } (\langle \langle v \rangle \rangle - v_\varGamma ) =0 \end{aligned}$$

(23)

$$\begin{aligned}&\displaystyle \langle \langle n \cdot A \nabla v \rangle \rangle + \frac{\alpha }{2}\llbracket v \rrbracket = 0 \end{aligned}$$

(24)

which are precisely the conditions used in [2].

Remark 4

The geometry of the interface may be generalized in several ways, which is needed in practical modeling of systems of cracks. For instance, we may consider bifurcating cracks where a so called Kirchhoff condition holds along the intersection, cracks that meet the boundary, and cracks which are piecewise smooth. We refer to [5, 8, 11, 24], for details on how to construct CutFEM for bifurcating cracks and how to handle the Kirchhoff condition weakly in a systematic manner. The regularity of the exact solution may be locally lower due to nonconvex corners and edges and therefore there may be a need for adaptive mesh refinement. More difficult to handle are cracks with boundary in the interior of the domain since they may lead to singularities in the bulk field, see [16, 17], and furthermore the geometry of the crack tip plays an important since it determines the boundary conditions at the crack tip see [30] for details. The properties of the solutions to problems with crack boundaries are interesting future research topics. We do however remark that the finite element method may be directly extended to cracks with boundaries if we assume a homogeneous Neumann condition at the crack tip, see [31] for numerical studies using a simpler but related method.

2.2 Weak form

Define the function spaces

$$\begin{aligned} V= & {} V_1 \oplus V_2 \oplus V_\varGamma \end{aligned}$$

(25)

$$\begin{aligned} V_i= & {} \{ v_i \in H^1(\varOmega _i) : v =0 \text { on}~~\partial \varOmega \cap \partial \varOmega _i \}\qquad i =1,2 \end{aligned}$$

(26)

$$\begin{aligned} V_\varGamma= & {} H^1(\varGamma ) \end{aligned}$$

(27)

and let $v\in V$ denote the vector $v=(v_1,v_2,v_\varGamma )$. We will also use the notation $\tilde{V}$ for functions $v \in V$ such that $v_i \in H^{\frac{3}{2}+\epsilon }(\varOmega _i)$, $i=1,2$, and $v_\varGamma \in H^{\frac{3}{2}+\epsilon }(\varGamma )$, with $\epsilon >0$. Using partial integration on $\varOmega _i$ we obtain

$$\begin{aligned}&\sum _{i=1}^2 (f_i,v_i)_{\varOmega _i} \end{aligned}$$

(28)

$$\begin{aligned}&\quad = \sum _{i=1}^2 (-\nabla \cdot A_i \nabla u_i, v_i)_{\varOmega _i} \end{aligned}$$

(29)

$$\begin{aligned}&\quad = \sum _{i=1}^2 ( A_i \nabla u_i, \nabla v_i)_{\varOmega _i} - (n_i \cdot A_i \nabla u_i,v_i)_{\partial \varOmega _i} \end{aligned}$$

(30)

$$\begin{aligned}&\quad = \sum _{i=1}^2 ( A_i \nabla u_i, \nabla v_i)_{\varOmega _i} \end{aligned}$$

(31)

$$\begin{aligned}&\qquad - (n_i \cdot A_i \nabla u_i,v_i -v_\varGamma )_{\partial \varOmega _i\cap \varGamma } - (n_i \cdot A_i \nabla u_i, v_\varGamma )_{\partial \varOmega _i\cap \varGamma } \end{aligned}$$

(32)

$$\begin{aligned}&\quad = \sum _{i=1}^2 ( A_i \nabla u_i, \nabla v_i)_{\varOmega _i} - (n\cdot A \nabla u, v - v_\varGamma )_{\varGamma } - (\llbracket n \cdot A \nabla u \rrbracket , v_\varGamma )_{\varGamma } \end{aligned}$$

(33)

$$\begin{aligned}&\quad = \sum _{i=1}^2 ( A_i \nabla u_i, \nabla v_i)_{\varOmega _i} + (B(u-u_\varGamma ), v - v_\varGamma )_{\varGamma } \end{aligned}$$

(34)

$$\begin{aligned}&\qquad + (A_\varGamma \nabla _\varGamma u_\varGamma , \nabla _\varGamma v_\varGamma )_{\varGamma } - (f_\varGamma ,v_\varGamma )_\varGamma \end{aligned}$$

(35)

Thus we arrive at the weak problem: find $u=(u_1,u_2,u_\varGamma ) \in V$ such that

$$\begin{aligned} \boxed {\mathcal {A}(u,v) = L(v)\quad \forall v \in V} \end{aligned}$$

(36)

where the forms are defined by

$$\begin{aligned} \mathcal {A}(u,v)= & {} \sum _{i=1}^2 ( A_i \nabla u_i, \nabla v_i)_{\varOmega _i} \end{aligned}$$

(37)

$$\begin{aligned}&+ (A_\varGamma \nabla _\varGamma u_\varGamma , \nabla _\varGamma v_\varGamma )_{\varGamma } + (B(u-u_\varGamma ), v - v_\varGamma )_{\varGamma } \end{aligned}$$

(38)

$$\begin{aligned} L(v)= & {} \sum _{i=1}^2 (f_i,v_i)_{\varOmega _i} + (f_\varGamma ,v_\varGamma )_\varGamma \end{aligned}$$

(39)

2.3 Existence, uniqueness, and regularity

Introducing the energy norm

$$\begin{aligned} |||v |||^2 = \sum _{i=1}^2 \Vert v \Vert ^2_{H^1(\varOmega _i)} + \Vert v_\varGamma \Vert ^2_{H^1(\varGamma )} + \Vert v - v_\varGamma \Vert ^2_{B,\varGamma } \end{aligned}$$

(40)

on V, we directly find using a Poincaré inequality and the Cauchy–Schwarz inequality that the form A is coercive and continuous

$$\begin{aligned} |||v |||^2 \lesssim \mathcal {A}(v,v), \qquad \mathcal {A}(v,w) \lesssim |||v \!|||\, |||w |||\end{aligned}$$

(41)

Furthermore, L is a continuous functional on V and it follows from the Lax–Milgram Lemma that there is a unique solution in V to (36).

Lemma 1

In the case considered here where $\varGamma$ is a smooth, closed interface the model problem (36) satisfies the elliptic regularity estimate

$$\begin{aligned} \boxed { \Vert u_1 \Vert _{H^2(\varOmega _1)} + \Vert u_2 \Vert _{H^2(\varOmega _2)} + \Vert u_\varGamma \Vert _{H^2(\varGamma )} \lesssim \Vert f_1 \Vert _{\varOmega _1} + \Vert f_2 \Vert _{\varOmega _2} + \Vert f_\varGamma \Vert _\varGamma } \end{aligned}$$

(42)

Under the additional assumption that B is a constant matrix and $A_i = a_i \mathbb {I}_{[d\times d]}$ $i=1,2$ and $A_\varGamma = a_\varGamma \mathbb {I}_{[(d-1)\times (d-1)]}$ with $a_i \in \mathbb {R}^+$ and $\mathbb {I}_{[d\times d]}$ the identity matrix, then the constant in (42) is independent of the coefficients of B.

Remark 5

The assumptions that $A_i = a_i \mathbb {I}_{[d\times d]}$ and B is constant along $\varGamma$ can be relaxed to smoothly varying coefficients, with additional technical work. See “Appendix C” for the case of variable B with uniformly bounded derivatives. We have not included the full analysis of the general case in the manuscript to keep the presentation at a reasonable level of complexity.

To prove (42) we first recall a partial integration formula from [19], see Eq. 3.1.1.1 on p. 134. For completeness we include a proof based on tangential calculus, which is in line with the notation used in this paper, in “Appendix A”.

Lemma 2

Let $\omega \subset \mathbb {R}^d$ be a domain with $C^2$ boundary $\partial \omega$ and exterior unit normal n. For $w \in [H^2(\omega )]^d$ it holds

$$\begin{aligned} (\nabla \cdot w, \nabla \cdot w)_\omega= & {} (w\otimes \nabla ,\nabla \otimes w)_\omega + 2(w_n, \text {div}_T w_T)_{\partial \omega } \end{aligned}$$

(43)

$$\begin{aligned}&+(w_T,w_T)_{\kappa ,\partial \omega } + (w_n,w_n)_{\text {tr}(\kappa ),\partial \omega } \end{aligned}$$

(44)

Here $w = w_T + w_n n$ is the decomposition of w into the tangential and normal components of the vector field w in an open neighborhood of the boundary; $\kappa = \nabla \otimes n = \nabla ^2 \zeta$ is the tangential curvature tensor of $\partial \omega$ and $\zeta$ is the signed distance function of $\partial \omega$ such that $n = \nabla \zeta$; and $\text {div}_T(w) = \text {tr}(w \otimes \nabla _T)$ is the surface divergence on $\partial \omega$.

Proof of Lemma 1

The proof consists of three steps: (1) Use the Lax–Milgram lemma to show existence and B independent stability in $H^1$. (2) Show that the $H^1$ solution is in fact in $H^2$. (3) Apply the partial integration identity (43) to derive a B independent estimate for the $H^2$ norm.

We will neglect the exterior boundary and focus our attention on the interface condition. The extension to the convex polygonal exterior boundary can be handled using standard techniques, see [19, Theorem 4.3.1.4] for a proof in the two dimensional case. For brevity we will also employ the notation $\nabla _n v = n \cdot \nabla v$ in the proof.

Step 1. Using the Lax–Milgram lemma there is a unique solution $(u_1 ,u_2,u_\varGamma ) \in H^1(\varOmega _2) \oplus H^1(\varOmega _1) \oplus H^1(\varGamma )$, which satisfies the energy bound

$$\begin{aligned} \sum _{i=1}^2 \Vert u_i \Vert ^2_{H^1(\varOmega _i)} + \Vert u_\varGamma \Vert ^2_{H^1(\varGamma )} + \Vert u - u_\varGamma \Vert ^2_{B,\varGamma } \lesssim \sum _{i=1}^2 \Vert f_i \Vert ^2_{\varOmega _i} + \Vert f_\varGamma \Vert ^2_\varGamma \end{aligned}$$

(45)

with hidden constant independent of B.

Step 2. Since $u_i \in H^1(\varOmega _i)$, $i=1,2,$ and $u_\varGamma \in H^1(\varGamma )$ we have $B(u - u_\varGamma )\vert _{\varGamma } \in [H^{\frac{1}{2}}(\varGamma )]^2$ and using (3), $\nabla _n u \in [H^{\frac{1}{2}}(\varGamma )]^2$. This means that the right hand side of (2) is in $L^2$ and hence $u_\varGamma \in H^2(\varGamma )$ by elliptic regularity. Considering once again (3) we see that in each subdomain the solution coincides with a single domain solution with a Robin condition with data in $H^{\frac{1}{2}}(\varGamma )$ on $\varGamma$. By the elliptic regularity of the Robin problem we can then conclude that $u_i \in H^2(\varOmega _i)$ and therefore $(u_1,u_2,u_\varGamma ) \in H^2(\varOmega _1) \otimes H^2(\varOmega _2) \otimes H^2(\varGamma )$.

Step 3. Setting $\omega = \varOmega _i$ and $w_i = a_i^{\frac{1}{2}}\nabla u_i$ in the partial integration identity (43) and summing the contributions from $\varOmega _1$ and $\varOmega _2$ we obtain

$$\begin{aligned} \sum _{i=1}^2 a_i \Vert \varDelta u_i\Vert ^2_{\varOmega _i}= & {} \sum _{i=1}^2 a_i \Vert \nabla ^2 u_i \Vert ^2_{\varOmega _i} + 2\sum _{i=1}^2( a_i \nabla _n u_i, \varDelta _\varGamma u_i)_{\partial \varOmega _i \cap \varGamma } \end{aligned}$$

(46)

$$\begin{aligned}&+\sum _{i=1}^2 (a_i \nabla _\varGamma u_i,\nabla _\varGamma u_i)_{\kappa ,\varGamma } + \sum _{i=1}^2(a_i \nabla _n u_i,\nabla _n u_i)_{\text {tr}(\kappa ),\varGamma } \end{aligned}$$

(47)

Recalling that $a_i \varDelta u_i = f_i$ and rearranging the terms we

$$\begin{aligned}&\sum _{i=1}^2 a_i \Vert \nabla ^2 u_i \Vert ^2_{\varOmega _i} + 2\underbrace{\sum \nolimits _{i=1}^2( a_i \nabla _n u_i, \varDelta _\varGamma u_i)_{\partial \varOmega _i \cap \varGamma }}_{I} \end{aligned}$$

(48)

$$\begin{aligned}&\;\; \le \sum _{i=1}^2 a_i^{-1} \Vert f_i \Vert ^2_{\varOmega _i} +\sum _{i=1}^2 \underbrace{|(a_i \nabla _\varGamma u_i, \nabla _\varGamma u_i)_{\kappa ,\varGamma }|}_{\textit{II}_i} + \sum _{i=1}^2 \underbrace{|(a_i \nabla _n u_i,\nabla _n u_i)_{\text {tr}(\kappa ),\varGamma }|}_{\textit{III}_i}\nonumber \\ \end{aligned}$$

(49)

We now proceed with estimates from below of term I, and estimates of $\textit{II}_i$ and $\textit{III}_i$. We let C denote a generic constant non necessarily the same in all occurrences.

Term I. Adding and subtracting $u_\varGamma$ we get

$$\begin{aligned} I= & {} \sum _{i=1}^2 ( a_i\nabla _{n_i} u_i, \varDelta _\varGamma (u_i - u_\varGamma ) )_{\partial \varOmega _i} + \sum _{i=1}^2 (a_i \nabla _{n_i} u_i, \varDelta _\varGamma u_\varGamma )_{\partial \varOmega _i} \end{aligned}$$

(50)

$$\begin{aligned}= & {} - ( B[u - u_\varGamma ],\varDelta _\varGamma [u - u_\varGamma ] )_{\varGamma } + \sum _{i=1}^2(\llbracket a \nabla _n u \rrbracket , \varDelta _\varGamma u_\varGamma )_{\varGamma } \end{aligned}$$

(51)

$$\begin{aligned}= & {} ( \nabla _\varGamma B[u - u_\varGamma ], \nabla _\varGamma [u - u_\varGamma ] )_{\varGamma } + \sum _{i=1}^2 (a_\varGamma \varDelta _\varGamma u_\varGamma + f_\varGamma , \varDelta _\varGamma u_\varGamma )_{\varGamma } \end{aligned}$$

(52)

$$\begin{aligned}\ge & {} \Vert \nabla _\varGamma [u - u_\varGamma ] \Vert ^2_{B,\varGamma } + a_\varGamma \frac{1}{2} \Vert \varDelta _\varGamma u_\varGamma \Vert ^2_\varGamma - \frac{1}{2} a_\varGamma ^{-1}\Vert f_\varGamma \Vert ^2_\varGamma \end{aligned}$$

(53)

Here we used the assumption that B is constant, which gives

$$\begin{aligned} -( B[u - u_\varGamma ], \varDelta _\varGamma [u - u_\varGamma ] )_{\varGamma }= & {} (\nabla _\varGamma (B[u - u_\varGamma ]), \nabla _\varGamma [u - u_\varGamma ] )_{\varGamma } \end{aligned}$$

(54)

$$\begin{aligned}= & {} ( B\nabla _\varGamma [u - u_\varGamma ], \nabla _\varGamma [u - u_\varGamma ] )_{\varGamma } \end{aligned}$$

(55)

$$\begin{aligned}= & {} \Vert \nabla _\varGamma [ u - u_\varGamma ] \Vert ^2_{B,\varGamma } \end{aligned}$$

(56)

For the second term

$$\begin{aligned} (a_\varGamma \varDelta _\varGamma u_\varGamma + f_\varGamma , \varDelta _\varGamma u_\varGamma )_{\varGamma }&\ge a_\varGamma \Vert \varDelta _\varGamma u_\varGamma \Vert ^2_\varGamma - \Vert f_\varGamma \Vert _\varGamma \Vert \varDelta _\varGamma u_\varGamma \Vert _{\varGamma } \end{aligned}$$

(57)

$$\begin{aligned}&\ge \frac{1}{2} a_\varGamma \Vert \varDelta _\varGamma u_\varGamma \Vert ^2_\varGamma - \frac{1}{2} a_\varGamma ^{-1} \Vert f_\varGamma \Vert ^2_\varGamma \end{aligned}$$

(58)

Term $\textit{II}_i$. Using interpolation [35, Proposition 3.1] followed by the trace inequality

$$\begin{aligned} \Vert w \Vert _{H^s(\varGamma )} \le C \Vert w \Vert _{H^{s+1/2}(\varOmega _i)}, \quad s > 0 \end{aligned}$$

(59)

we get

$$\begin{aligned} \textit{II}_i= & {} |(a_i \nabla _\varGamma u_i,\nabla _\varGamma u_i)_{\kappa ,\varGamma }| \end{aligned}$$

(60)

$$\begin{aligned}\le & {} \Vert \kappa \Vert _{L^\infty (\varGamma )} (a_i \nabla _\varGamma u_i,\nabla _\varGamma u_i)_{\varGamma } \end{aligned}$$

(61)

$$\begin{aligned}\le & {} C a_i \Vert u_i \Vert _{H^{3/2}(\varGamma )} \Vert u_i \Vert _{H^{1/2}(\varGamma )} \end{aligned}$$

(62)

$$\begin{aligned}\le & {} C a_i \Vert u_i \Vert _{H^2(\varOmega _i)} \Vert u_i \Vert _{H^1(\varOmega _i)} \end{aligned}$$

(63)

$$\begin{aligned}\le & {} \frac{1}{4} a_i \Vert u_i \Vert ^2_{H^2(\varOmega _i)} + C a_i \Vert u_i \Vert ^2_{H^1(\varOmega _i)} \end{aligned}$$

(64)

Term $\textit{III}_i$. Again using interpolation and then by applying the trace inequalities (59) and

$$\begin{aligned} \Vert \nabla _{n_i} u_i\Vert _{H^{-1/2}(\varGamma )} \le C(\Vert \nabla u_i\Vert ^2_{\varOmega _i} + \Vert \varDelta u_i\Vert ^2_{\varOmega _i})^{1/2} \end{aligned}$$

(65)

we get

$$\begin{aligned} \textit{III}_i= & {} |a_i (\nabla _n u_i,\nabla _n u_i)_{\text {tr}(\kappa ),\varGamma }| \end{aligned}$$

(66)

$$\begin{aligned}\le & {} \Vert \text {tr}(\kappa )\Vert _{L^\infty (\varGamma )} a_i (\nabla _n u_i,\nabla _n u_i)_\varGamma \end{aligned}$$

(67)

$$\begin{aligned}\le & {} C a_i \Vert \nabla _n u_i\Vert _{H^{1/2}(\varGamma )} \Vert \nabla _n u_i\Vert _{H^{-1/2}(\varGamma )} \end{aligned}$$

(68)

$$\begin{aligned}\le & {} C a_i \Vert \nabla u_i\Vert _{H^{1}(\varOmega _i)} (\Vert \nabla u_i \Vert ^2_{\varOmega _i} + \Vert \varDelta u_i \Vert ^2_{\varOmega _i} )^{1/2} \end{aligned}$$

(69)

$$\begin{aligned}\le & {} \frac{1}{4} a_i \Vert u_i\Vert ^2_{H^{2}(\varOmega _i)} + C a_i (\Vert \nabla u_i \Vert ^2_{\varOmega _i} + \Vert \varDelta u_i \Vert ^2_{\varOmega _i} ) \end{aligned}$$

(70)

Here we used the estimate $\Vert \nabla _n u_i \Vert _\varGamma \le C \Vert \nabla u_i \Vert _{\varOmega _i}$, which we prove in “Appendix B”.

Conclusion. Starting from (49) and using the estimate (53) to bound term I from below and estimates (64) and (70) to bound terms $\textit{II}_i$ and $\textit{III}_i$ from above we obtain

$$\begin{aligned}&\sum _{i=1}^2 a_i \Vert \nabla ^2 u_i \Vert ^2_{\varOmega _i} + \Vert \nabla _\varGamma [u - u_\varGamma ] \Vert ^2_{B,\varGamma } + \frac{1}{2} a_\varGamma \Vert \varDelta _\varGamma u_\varGamma \Vert ^2_\varGamma - \frac{1}{2} a_\varGamma ^{-1} \Vert f_\varGamma \Vert ^2_\varGamma \end{aligned}$$

(71)

$$\begin{aligned}&\quad \le \sum _{i=1}^2 a_i^{-1} \Vert f_i \Vert ^2_{\varOmega _i} +\sum _{i=1}^2 \frac{1}{4} a_i \Vert u_i \Vert ^2_{H^2(\varOmega _i)} + C a_i \Vert u_i \Vert ^2_{H^1(\varOmega _i)} \end{aligned}$$

(72)

$$\begin{aligned}&\qquad + \sum _{i=1}^2 \frac{1}{4} a_i \Vert u_i\Vert ^2_{H^2(\varOmega _i)} + C a_i (\Vert \nabla u_i \Vert ^2_{\varOmega _i} + \Vert \varDelta u_i \Vert ^2_{\varOmega _i} ) \end{aligned}$$

(73)

Reorganizing the terms and using the identity $a_i \varDelta u_i = f_i$, and writing $\Vert u_i \Vert ^2_{H^2(\varOmega _i)} = \Vert \nabla ^2 u_i \Vert ^2_{\varOmega _i} + \Vert u_i \Vert ^2_{H^1(\varOmega _i)}$ we get

$$\begin{aligned}&\sum _{i=1}^2 a_i \Vert \nabla ^2 u_i \Vert ^2_{\varOmega _i} + \Vert \nabla _\varGamma [u - u_\varGamma ] \Vert ^2_{B,\varGamma } + a_\varGamma \Vert \varDelta _\varGamma u_\varGamma \Vert ^2_\varGamma \end{aligned}$$

(74)

$$\begin{aligned}&\qquad \le \frac{1}{2} a_\varGamma ^{-1} \Vert f_\varGamma \Vert ^2_\varGamma + \sum _{i=1}^2 C a_i^{-1}\Vert f_i \Vert ^2_{\varOmega _i} + \sum _{i=1}^2 C a_i \Vert u_i \Vert ^2_{H^1(\varOmega _i)} \end{aligned}$$

(75)

$$\begin{aligned}&\qquad \lesssim a_\varGamma ^{-1} \Vert f_\varGamma \Vert ^2_\varGamma + \sum _{i=1}^2(1 + a_i^{-1}) \Vert f_i \Vert ^2_{\varOmega _i} \end{aligned}$$

(76)

where we used the energy stability (45) in the final step, and the hidden constant depends on the trace constants and the curvature constant $\kappa$.

3 A robust finite element method

3.1 The mesh and finite element spaces

To formulate the finite element method we introduce the following notation:

• Let $\mathcal {T}_{h,0}$ be a quasiuniform mesh on $\varOmega$ with mesh parameter $h \in (0,h_0]$. Define the active meshes

  $$\begin{aligned} \mathcal {T}_{h,i}= & {} \{ T \in \mathcal {T}_{h,0} : T \cap \varOmega _i \ne \emptyset \} \quad i=1,2 \end{aligned}$$

  (77)
  $$\begin{aligned} \mathcal {T}_{h,\varGamma }= & {} \{ T \in \mathcal {T}_{h,0} : T \cap \varGamma \ne \emptyset \} \end{aligned}$$

  (78)

  associated with the bulk domains $\varOmega _i$, $i=1,2,$ and the interface $\varGamma$, and the domains covered by the meshes

  $$\begin{aligned} O_{h,i} = \cup _{T\in \mathcal {T}_{h,i}} \quad i=1,2, \qquad O_{h,\varGamma } = \cup _{T\in \mathcal {T}_{h,\varGamma }} \end{aligned}$$

  (79)

• Let

  $$\begin{aligned} \mathcal {T}_{h,i}(\varGamma ) = \{ T \in \mathcal {T}_{h,i} : T \cap \varGamma \ne \emptyset \} \end{aligned}$$

  (80)

• Let $\mathcal {F}_{h,i}$ be the set of all interior faces in $\mathcal {T}_{h,i}$ associated with an element in $\mathcal {T}_{h,i}(\partial \varOmega _i) = \{ T \in \mathcal {T}_{h,i} : T \cap \partial \varOmega _i \ne \emptyset \}$.

• Let $\mathcal {F}_{h,\varGamma }$ be the set of all interior faces in $\mathcal {T}_{h,\varGamma }$ and $\mathcal {K}_{h,\varGamma } = \{ K = T \cap \varGamma : T \in \mathcal {T}_{h,\varGamma } \}$.

• Let $V_{h,0}$ be the space of continuous piecewise linear functions on $\mathcal {T}_{h,0}$ and define

  $$\begin{aligned} V_{h,i} = V_{h,0}|_{\mathcal {T}_{h,i}}\quad i=1,2, \quad V_{h,\varGamma } = V_{h,0}|_{\mathcal {T}_{h,\varGamma }} \end{aligned}$$

  (81)

  and

  $$\begin{aligned} V_h = V_{h,1} \oplus V_{h,2} \oplus V_{h,\varGamma } \end{aligned}$$

  (82)

  For $V_{h,1}$ we also impose the homogeneous boundary condition on $\partial \varOmega$ strongly, i.e., we require $v = 0$ on $\partial \varOmega$.

3.2 Standard formulation

The standard finite element method takes the form: find $u_h = (u_{h,1},u_{h,2},u_{h,\varGamma }) \in V_h = V_{h,1} \oplus V_{h,2} \oplus V_{h,\varGamma }$ such that

$$\begin{aligned} \mathcal {A}^S_h(u_h,v) = L(v) \quad \forall v \in V_h \end{aligned}$$

(83)

Here the form $\mathcal {A}^S_h$ is defined by

$$\begin{aligned} \mathcal {A}^S_h = \mathcal {A} + s_h \end{aligned}$$

(84)

where $s_h$ is a stabilization term of the form

$$\begin{aligned} s_h = s_{h,1} + s_{h,2} + s_{h,\varGamma } \end{aligned}$$

(85)

with

$$\begin{aligned} s_{h,i}(v,w) = \sum _{F \in \mathcal {F}_{h,i} } h_F \Vert \zeta (A_i)\Vert _{\infty ,F} ([ n\cdot \nabla v ], [ n\nabla w ])_{F} \qquad i=1,2 \end{aligned}$$

(86)

where $\zeta (X)$ denotes the maximum eigenvalue of the matrix X,

$$\begin{aligned} s_{h,\varGamma }(v,w)= & {} \sum _{F \in \mathcal {F}_{h,\varGamma } } h_F \Vert \zeta (A_\varGamma )\Vert _{\infty ,F\cap \varGamma } ([ n\cdot \nabla v ], [ n \cdot \nabla w ])_{\mathcal {F}_{h,\varGamma }} \end{aligned}$$

(87)

$$\begin{aligned}&+ \sum _{T \in \mathcal {T}_{h,\varGamma } }h_K^2\Vert \zeta (A_\varGamma )\Vert _{\infty ,K\cap \varGamma } (n_\varGamma \cdot \nabla v, n_\varGamma \cdot \nabla w)_{T \cap \varGamma } \end{aligned}$$

(88)

where for a face sharing two elements $T_1$ and $T_2$ the jump $[n\cdot \nabla v]$ is defined by

$$\begin{aligned}{}[n \cdot \nabla v ] = n_1 \cdot \nabla v_1 + n_2 \cdot \nabla v_2 \end{aligned}$$

(89)

where $n_i$ is the exterior normal to $\partial T_i$ and $v_i = v|_{T_i}$.

3.2.1 Properties of the stabilization terms

The rationale for the design of the stabilizing terms is that they improve the stability, while remaining consistent for sufficiently smooth solutions.

Accuracy relies on the following consistency property that is immediate from the definitions above. For any function $v \in H^{\frac{3}{2}+\epsilon }(O_{h,i})$ there holds $s_{h,i}(v,w) = 0$ for all $w \in V_{h,i}+H^{\frac{3}{2}+\epsilon }(O_{h,i})$, $i=1,2$. For any function $v \in H^{\frac{3}{2}+\epsilon }(O_{h,\varGamma })$, such that $n_\varGamma \cdot \nabla v =0$ on $\varGamma$ there holds $s_{h,\varGamma }(v,w) = 0$ for all $w \in V_{h,\varGamma }+H^{\frac{3}{2}+\epsilon }(O_{h,\varGamma })$.

The stability properties are well known and we collect them in the following Lemma.

Lemma 3

There are constants such that

$$\begin{aligned} \Vert \nabla v \Vert ^2_{A_i, O_{h,i} } \lesssim \Vert \nabla v \Vert ^2_{A_i,\varOmega _i} + \Vert v \Vert ^2_{s_{h,i}} \quad i =1,2 \end{aligned}$$

(90)

and

$$\begin{aligned} \Vert \nabla _\varGamma v \Vert ^2_{A_\varGamma , O_{h,\varGamma }} \lesssim \Vert \nabla _\varGamma v \Vert ^2_{A_\varGamma ,\varGamma } + \Vert v \Vert ^2_{s_{h,\varGamma }} \end{aligned}$$

(91)

where we introduced the (semi) norm $\Vert v \Vert ^2_{s_h} = s_h(v,v)$.

Proof

See [6, 7, 28], with minor modifications to account for the varying coefficients.

Remark 6

Observe that the hidden constants in Lemma 3 depend on the variation of the $A_i$ and $A_\varGamma$.

Remark 7

The finite element methods we develop in this paper directly extends to higher continuous piecewise polynomial spaces with the modification that the stabilizing terms controls jumps in higher derivatives across faces and for the surface stabilization higher order normal derivatives must be added, see [28] for full details.

3.3 Robust formulation

The stabilizing terms ensure robustness irrespective of the intersection of the fracture and the mesh, but they do not counter instabilities due to degenerate B. Our aim is to design a formulation which is robust in the case when the eigenvalues of B degenerate. Indeed as we saw above as $\xi$ approaches 1/2, $\lambda _1$ blows up. For clarity we recall the abstract boundary condition

$$\begin{aligned} n\cdot A\nabla v + B (v - v_\varGamma ) = 0 \end{aligned}$$

(92)

where we now assume that the matrix B is a positive definite symmetric $2\times 2$ matrix with eigenvalues $\lambda _i$ and eigenvectors $e_i$, such that $\lambda _i \in (0,\infty )$ and thus one or both eigenvalues may become very large or small. To handle this situation we instead enforce

$$\begin{aligned} B^{-1} n\cdot A\nabla v + (v - v_\varGamma ) = 0 \end{aligned}$$

(93)

weakly using a modified Nitsche method. This approach was originally developed in [26] where fitted finite element approximation of Robin conditions were considered.

Derivation of an Alternative Weak Form. As before we have the identity

$$\begin{aligned} L(v)= & {} \underbrace{\sum _{i=1}^2 ( A_i \nabla u_i, \nabla v_i)_{\varOmega _i} + (A_\varGamma \nabla _\varGamma u_\varGamma , \nabla _\varGamma v_\varGamma )_{\varGamma }}_{=:\mathcal {A}_1(u,v)} - (n\cdot A \nabla u, v - v_\varGamma )_{\varGamma } \end{aligned}$$

(94)

$$\begin{aligned}= & {} \mathcal {A}_1(u,v) - (n\cdot A \nabla u, v - v_\varGamma )_{\varGamma } \end{aligned}$$

(95)

where we introduced the bilinear form $\mathcal {A}_1$ for brevity. To enforce the interface conditions we proceed as follows

$$\begin{aligned} L(v)= & {} \mathcal {A}_1(u,v) - (n\cdot A \nabla u, v - v_\varGamma )_{\varGamma } \end{aligned}$$

(96)

$$\begin{aligned}= & {} \mathcal {A}_1(u,v) + (n\cdot A \nabla u, B^{-1} (n \cdot A \nabla v) )_{\varGamma } \end{aligned}$$

(97)

$$\begin{aligned}&-\, \,(n\cdot A \nabla u, B^{-1} (n \cdot A \nabla v) + (v - v_\varGamma ) )_{\varGamma } \end{aligned}$$

(98)

$$\begin{aligned}= & {} \mathcal {A}_1(u,v) + (n\cdot A \nabla u, B^{-1} (n \cdot A \nabla v) )_{\varGamma } \end{aligned}$$

(99)

$$\begin{aligned}&-\,\, (n\cdot A \nabla u, B^{-1} (n \cdot A \nabla v) + (v - v_\varGamma ) )_{\varGamma } \end{aligned}$$

(100)

$$\begin{aligned}&-\, \,(B^{-1} (n \cdot A \nabla u) + (u - u_\varGamma ), n\cdot A \nabla v )_{\varGamma } \end{aligned}$$

(101)

$$\begin{aligned}&+\, (B^{-1} (n \cdot A \nabla u) + (u - u_\varGamma ) , \tau (B^{-1} (n \cdot A \nabla v) + (v - v_\varGamma )) )_{\varGamma } \end{aligned}$$

(102)

where the last two terms are zero due to the interface condition and the resulting form on the right hand side is symmetric. Furthermore, $\tau$ is a stabilization parameter (a $2\times 2$ matrix) of the form

$$\begin{aligned} \tau = \sum _{i=1}^2 \tau _i e_i \otimes e_i , \quad \tau _i = \frac{\lambda _i \beta }{ \lambda _i h + \beta } \quad i=1,2 \end{aligned}$$

(103)

where $\beta$ is a positive parameter and we recall that $\lambda _i$ and $e_i$ are the eigenvalues and eigenvectors of B. The parameter $\beta$ is chosen so that

$$\begin{aligned} \Vert n \Vert ^2_{A,\infty ,\varGamma } = \sum _{i=1}^2 \Vert n_i \Vert ^2_{A_i,\infty ,\varGamma } \lesssim \beta \end{aligned}$$

(104)

where $\Vert w \Vert _{A_i,\infty ,\varGamma } = \Vert (w \cdot A_i \cdot w)^{1/2} \Vert _{\infty , \varGamma }$. We next show some technical estimates for the stabilization parameter, in particular, we show that $\tau$ is uniformly bounded as a function of $\lambda _i\in (0,\infty )$.

Lemma 4

For all eigenvalues $\lambda _i \in (0,\infty )$, $i=1,2,$ of B, the following estimates related to the stabilization parameter $\tau$ hold

$$\begin{aligned}&\Vert B^{-1} \tau B^{-1} - B^{-1}\Vert _{L^\infty (\varGamma )} \le \frac{h}{\beta } \end{aligned}$$

(105)

$$\begin{aligned}&\Vert (B^{-1} \tau - I) \tau ^{-1/2} \Vert _{L^\infty (\varGamma )} \le \left( \frac{h}{\beta }\right) ^{1/2} \end{aligned}$$

(106)

$$\begin{aligned}&\Vert \tau \Vert _{L^\infty (\varGamma )} \le \frac{\beta }{h} \end{aligned}$$

(107)

Proof

First we recall that for any symmetric matrix D it holds

$$\begin{aligned} \Vert D \Vert _{\mathbb {R}^d} \lesssim \max _i |\gamma _i| \end{aligned}$$

(108)

where $\gamma _{i}$ are the eigenvalues of D. To prove (105) we write B in terms of its eigenvalues $\lambda _i$ and eigenvectors $e_i$,

$$\begin{aligned} B = \sum _{i=1}^2 \lambda _i e_i \otimes e_i \end{aligned}$$

(109)

and using the definition (103) of $\tau$ we obtain the identity

$$\begin{aligned} B^{-1} \tau B^{-1} - B^{-1} = \sum _{i=1}^2 \left( \frac{\tau _i}{\lambda _i} - 1 \right) \frac{1}{\lambda _i} e_i \otimes e_i \end{aligned}$$

(110)

Here we have the following estimate of the eigenvalues

$$\begin{aligned} \left| \left( \frac{\tau _i}{\lambda _i} - 1 \right) \frac{1}{\lambda _i} \right| = \left| \left( \frac{\beta }{\lambda _i h + \beta } - 1 \right) \frac{1}{\lambda _i} \right| = \frac{h}{\lambda _i h + \beta } \le \frac{h}{\beta } \end{aligned}$$

(111)

which in view of (108) completes the verification of (105). Next, for (106) we have

$$\begin{aligned} (B^{-1} \tau - I) \tau ^{-1/2} = \sum _{i=1}^2 \left( \frac{\tau _i}{\lambda _i} - 1 \right) \frac{1}{\tau _i^{1/2}} e_i \otimes e_i \end{aligned}$$

(112)

and

$$\begin{aligned}&\left| \left( \frac{\tau _i}{\lambda _i} - 1 \right) \frac{1}{\tau _i^{1/2}} \right| = \left| \left( \frac{\beta }{\lambda _i h + \beta } - 1 \right) \left( \frac{\lambda _i h + \beta }{\lambda _i \beta } \right) ^{1/2} \right| \end{aligned}$$

(113)

$$\begin{aligned}&\qquad = \frac{\lambda _i h}{\lambda _i h + \beta }\left( \frac{\lambda _i h + \beta }{\lambda _i \beta } \right) ^{1/2} = \left( \frac{\lambda _i h}{\lambda _i h + \beta }\right) ^{1/2} \left( \frac{h}{\beta } \right) ^{1/2} \le \left( \frac{h}{\beta } \right) ^{1/2} \end{aligned}$$

(114)

which proves (106). The final bound (107) is a direct consequence of the definition of $\tau$ and the estimate

$$\begin{aligned} \frac{\lambda _i \beta }{\lambda _i h + \beta } \le \frac{\lambda _i \beta }{\lambda _i h} \le \frac{\beta }{h} \end{aligned}$$

(115)

Remark 8

The choice of $\tau _i$ can be further refined as follows

$$\begin{aligned} \tau _i = \frac{\lambda _i \beta _i}{ \lambda _i h + \beta _i} \quad i=1,2 \end{aligned}$$

(116)

with

$$\begin{aligned} \sum _{j=1}^2 \Vert n_j \Vert ^2_{A_j,\infty ,\varGamma } |e_{ij}|^2 \lesssim \beta _i \end{aligned}$$

(117)

where $e_i = [e_{i1}\; e_{i2}]^T$. This approach is beneficial in situations where the components of $e_i$ are very different and there is a large difference between the $\Vert n_j \Vert ^2_{A_j,\infty ,\varGamma }$ with $j=1$ and $j=2$.

The robust finite element method. Find $u_h \in V_h$ such that

$$\begin{aligned} \boxed {\mathcal {A}^R_h(u_h,v)=\mathcal {A}^R(u_h,v) + s_h(u_h,v) = L(v) \qquad \forall v \in V_h} \end{aligned}$$

(118)

where

$$\begin{aligned}&\mathcal {A}^R(v,w) \end{aligned}$$

(119)

$$\begin{aligned}&\quad = \mathcal {A}_1(v,w) + (n\cdot A \nabla v, B^{-1} (n \cdot A \nabla w) )_{\varGamma } \end{aligned}$$

(120)

$$\begin{aligned}&\qquad - (n\cdot A \nabla v, B^{-1} (n \cdot A \nabla w) + (w - w_\varGamma ) )_{\varGamma } \end{aligned}$$

(121)

$$\begin{aligned}&\qquad - (B^{-1} (n \cdot A \nabla v) + (v - v_\varGamma ), n\cdot A \nabla w )_{\varGamma } \end{aligned}$$

(122)

$$\begin{aligned}&\qquad + (B^{-1} (n \cdot A \nabla v) + (v - v_\varGamma ) , \tau (B^{-1} (n \cdot A \nabla w) + (w - w_\varGamma )) )_{\varGamma } \end{aligned}$$

(123)

It follows by the design of $\mathcal {A}^R$ that for a sufficiently smooth exact solution $u \in \tilde{V}$ of the problem (36) there holds

$$\begin{aligned} \mathcal {A}(u,v) = \mathcal {A}^R(u,v) = L(v), \quad \forall v \in (V \cap H^2(\varOmega _1\cup \varOmega _2 \cup \varGamma ) + V_h \end{aligned}$$

(124)

As a consequence we immediately get the Galerkin orthogonality

Lemma 5

Let $u \in \tilde{V}$ be the solution of (36) and $u_h \in V_h$ the solution of (118) then there holds

$$\begin{aligned} \mathcal {A}^R(u - u_h,v) = s_h(u_h,v)\quad \forall v \in V_h \end{aligned}$$

(125)

Proof

The proof follows by combining (124) and (118).

4 Error estimates

4.1 The energy norm

We introduce the energy norm

$$\begin{aligned} |||v |||_h^2= & {} \sum _{i=1}^2 \Vert \nabla v_i \Vert ^2_{A_i,\varOmega _i} + h \Vert \nabla v_i \Vert ^2_{A_i,\varGamma } + \Vert v \Vert ^2_{s_{h}} \end{aligned}$$

(126)

$$\begin{aligned}&+ \,\Vert \nabla _\varGamma v_\varGamma \Vert ^2_{A_\varGamma ,\varGamma } + \Vert v - v_\varGamma \Vert ^2_{\tau ,\varGamma } \end{aligned}$$

(127)

where $\Vert w \Vert ^2_{\psi ,\omega } = \int _\omega \psi w^2$ is the $\psi$ weighted $L^2$ norm over the set $\omega$.

4.2 Interpolation error estimates

We begin by introducing interpolation operators and derive the basic approximation error estimates. Then collecting the estimates we show an interpolation error estimate in the energy norm (126127). The basic idea in the construction of the interpolation operators is to use an extension of the solution outside of the domain and then employ a syandard weak type interpolation operator. The extension operator is stable with respect to Sobolev norms.

The Scott–Zhang interpolant. Given a mesh $\mathcal {T}_h$ covering a domain $O_h$ and the space of piecewise linear continuous finite elements $W_h$, a Scott–Zhang interpolation operator $\pi _{h,SZ} : H^1(\varOmega _h) \rightarrow W_h$ satisfies the element wise estimate

$$\begin{aligned} \Vert v - \pi _{h,SZ} v \Vert _{H^m(T)} \lesssim h^{2-m} \Vert v \Vert _{H^2(\mathcal {N}(T))}, \quad m=0,1 \end{aligned}$$

(128)

where $\mathcal {N}(T)$ is the set of all elements in $\mathcal {T}_{h,i}$ that share a node with T. More precisely, we employ a Scott–Zhang interpolation operator that averages on all elements sharing a node for interior nodes, and for nodes at the boundary $\partial \varOmega$ the average is taken over all faces on the boundary sharing the node leading to exact preservation of homogeneous boundary conditions. See [33] for further details.

Bulk domain fields. It is shown in [34, Section 2.3, Theorem 5] that there is an extension operator $E_i: H^s(\varOmega _i)\rightarrow H^s(\mathbb {R}^d)$, not dependent on $s\ge 0$, which is stable in the sense that

$$\begin{aligned} \Vert E_i v_i \Vert _{H^s(\mathbb {R}^d)} \lesssim \Vert v_i \Vert _{H^s(\varOmega _i)} \end{aligned}$$

(129)

We define the interpolation operator $\pi _{h,i}:H^1(\mathcal {N}(\mathcal {T}_{h,i})) \rightarrow V_{h,i}$ by

$$\begin{aligned} \pi _{h,i} v_i = \pi _{h,i,SZ} E_i v \end{aligned}$$

(130)

where $\pi _{h,i,SZ}:H^1(\mathcal {N}(\mathcal {T}_{h,i})) \rightarrow V_{h,i}$ is the Scott–Zhang interpolant and $\mathcal {N}(\mathcal {T}_{h,i})$ is the set of elements sharing a node with an element in $\mathcal {T}_{h,i}$. We then have the error estimate

$$\begin{aligned} \boxed { \Vert v_i - \pi _{h,i} v \Vert _{H^m(\varOmega _i)} \lesssim h^{2-m} \Vert v_i \Vert _{H^2(\varOmega _i)}\quad m=0,1} \end{aligned}$$

(131)

Proof

Using the notation $\rho _i = v_i - \pi _{h,i} v_i$ we obtain

$$\begin{aligned} \Vert \rho _i \Vert _{H^m(\varOmega _i)} \lesssim \Vert \rho _i \Vert _{H^m(\mathcal {T}_{h,i})} \lesssim h^{2-m} \Vert E_i u_i \Vert _{H^2(\mathcal {N}(\mathcal {T}_{h,i}))} \lesssim h^{2-m} \Vert u_i \Vert ^2_{H^2(\varOmega _i)} \end{aligned}$$

(132)

where we used the interpolation error estimate (128) and finally the stability (129) of the extension operator $E_i$.

Interface field. Let $p_\varGamma : U_\delta (\varGamma ) \rightarrow \varGamma$ be the closest point mapping from the tubular neighborhood $U_\delta (\varGamma ) =\{x: \text{ dist }(x,\varGamma ) < \delta \}$ to $\varGamma$, which is well defined for all $\delta \in (0,\delta _0]$ for some $\delta _0>0$. Define the extension operator $E_\varGamma : L^2(\varGamma ) \rightarrow L^2(U_{\delta }(\varGamma ))$ by $E_\varGamma v = v \circ p_\varGamma$. Since $\varGamma$ is smooth we have the stability estimate

$$\begin{aligned} \Vert E_\varGamma v_\varGamma \ \Vert _{H^s(U_\delta (\varGamma ))} \lesssim \delta ^{1/2} \Vert v_\varGamma \Vert _{H^s(\varGamma )} \end{aligned}$$

(133)

Observe also that since $n_\varGamma \cdot \nabla E_\varGamma v_\varGamma = 0$ by construction, and then assuming $s>3/2$ in (133), we see that

$$\begin{aligned} s_{h,\varGamma }(E_\varGamma v_\varGamma ,w) = 0, \qquad \forall w \in V_{h,\varGamma }+H^{\frac{3}{2}+\epsilon }(O_{h,\varGamma }) \end{aligned}$$

(134)

We define the interpolation operator $\pi _{h,\varGamma }:H^1(\varGamma ) \rightarrow V_{h,\varGamma }$ by

$$\begin{aligned} \pi _{h,\varGamma } v_\varGamma =( \pi _{h,\varGamma ,SZ} E_\varGamma v_\varGamma )|_{O_{h,\varGamma }} \end{aligned}$$

(135)

Here $\pi _{h,\varGamma ,SZ}: H^1(\mathcal {N}(\mathcal {T}_{h,\varGamma })) \rightarrow V_{h,\varGamma }$ is the Scott–Zhang interpolation operator defined on the set $\mathcal {N}(\mathcal {T}_{h,\varGamma })$ of all elements that share a node with an element in $\mathcal {T}_{h,\varGamma }$. We also note that there is $\delta \sim h$ such that

$$\begin{aligned} \mathcal {N}(\mathcal {T}_{h,\varGamma }) \subset U_\delta (\varGamma ) \end{aligned}$$

(136)

We have the error estimate

$$\begin{aligned} \boxed { \Vert v - \pi _{h,\varGamma } v \Vert _{H^m(\varGamma )} \lesssim h^{2-m} \Vert v \Vert _{H^2(\varGamma )}\quad m=0,1} \end{aligned}$$

(137)

Proof

Using the element wise trace inequality

$$\begin{aligned} \Vert w \Vert ^2_{\varGamma \cap T} \lesssim h^{-1} \Vert w \Vert ^2_{T} + h \Vert \nabla w \Vert ^2_{T} \quad w \in H^1(T) \end{aligned}$$

(138)

see [23, 36], we obtain after summation over all elements $T \in \mathcal {T}_{h,\varGamma }$,

$$\begin{aligned} \Vert w \Vert ^2_{\varGamma } \lesssim h^{-1} \Vert w \Vert ^2_{\mathcal {T}_{h,\varGamma }} + h \Vert \nabla w \Vert ^2_{\mathcal {T}_{h,\varGamma }} \qquad w \in H^1(\mathcal {T}_{h,\varGamma }) \end{aligned}$$

(139)

Applying (139) with $w = \nabla _\varGamma ^m \rho$ gives

$$\begin{aligned} \Vert \nabla _\varGamma ^m \rho \Vert ^2_{\varGamma }\lesssim & {} h^{-1} \Vert \nabla ^m \rho \Vert ^2_{\mathcal {T}_{h,\varGamma }} +\delta \Vert \nabla ^{m+1} \rho \Vert ^2_{\mathcal {T}_{h,\varGamma }} \end{aligned}$$

(140)

$$\begin{aligned}\lesssim & {} (h^{-1} h^{2(2-m)} +h h^{2(1-m)}) \Vert E_\varGamma v\Vert ^2_{H^2(\mathcal {N}(\mathcal {T}_{h,\varGamma }))} \end{aligned}$$

(141)

$$\begin{aligned}\lesssim & {} (h^{-1} h^{2(2-m)} +h h^{2(1-m)}) \Vert E_\varGamma v\Vert ^2_{H^2(U_\delta (\varGamma ))} \end{aligned}$$

(142)

$$\begin{aligned}\lesssim & {} \delta (h^{-1} h^{2(2-m)} +h h^{2(1-m)}) \Vert E_\varGamma v\Vert ^2_\varGamma \end{aligned}$$

(143)

$$\begin{aligned}\lesssim & {} h^{2(2-m)} \Vert v \Vert ^2_{H^2(\varGamma )} \end{aligned}$$

(144)

where we used, the interpolation error estimate (128), the inclusion (136), and the stability (133) of the extension operator $E_\varGamma$.

We finally define the interpolation operator $\pi _h: V \rightarrow V_h$ as follows

$$\begin{aligned} \pi _h v = (\pi _{h,1} E_1 v_1, \pi _{h,2} E_2 v_2, \pi _{h,\varGamma } E_\varGamma v_\varGamma ) \end{aligned}$$

(145)

Lemma 6

There is a constant not dependent on the matrix B, in the interface condition (3), such that

$$\begin{aligned} \boxed { |||v - \pi _h v |||_h \lesssim h \left( \sum _{i=1}^2 \Vert v_i \Vert _{H^2(\varOmega _i)} + \Vert v_\varGamma \Vert _{H^2(\varGamma )} \right) } \end{aligned}$$

(146)

Proof

Let $v - \pi _h v = \rho$ be the interpolation error. Using the triangle inequality, the estimate $\tau \lesssim h^{-1}$, see (107), and the fact that the coefficients $A_i$ are uniformly bounded, we obtain

$$\begin{aligned} |||\rho |||_h^2= & {} \sum _{i=1}^2 \Vert \nabla \rho _i \Vert ^2_{A_i,\varOmega _i} + h \Vert \nabla \rho _i \Vert ^2_{A_i,\varGamma } \end{aligned}$$

(147)

$$\begin{aligned}&+\, \Vert \nabla _\varGamma \rho _\varGamma \Vert ^2_{A_\varGamma ,\varGamma } +\Vert \rho _i - \rho _\varGamma \Vert ^2_{\tau ,\varGamma } + \Vert \rho \Vert ^2_{s_h} \end{aligned}$$

(148)

$$\begin{aligned}\lesssim & {} \sum _{i=1}^2 \underbrace{\Vert \nabla \rho _i \Vert ^2_{\varOmega _i}}_{I_i} + \underbrace{h \Vert \nabla \rho _i \Vert ^2_{\varGamma } + h^{-1} \Vert \rho _i \Vert ^2_{\varGamma }}_{II_i} + \underbrace{\Vert \rho _i \Vert ^2_{s_{h,i}}}_{III_i} \end{aligned}$$

(149)

$$\begin{aligned}&+\, \underbrace{\Vert \nabla _\varGamma \rho _\varGamma \Vert ^2_{\varGamma } +h^{-1} \Vert \rho _\varGamma \Vert ^2_{\varGamma }}_{IV} + \underbrace{\Vert \rho _\varGamma \Vert ^2_{s_{h,\varGamma }}}_{V} \end{aligned}$$

(150)

We proceed with estimates of the terms on the right hand side. Using (131) we directly have

$$\begin{aligned} I_i = \Vert \nabla \rho _i \Vert ^2_{\varOmega _i} \lesssim h^2 \Vert v_i \Vert _{H^2(\varOmega _i)} \end{aligned}$$

(151)

Next using the trace inequality

$$\begin{aligned} \Vert w \Vert ^2_\varGamma \lesssim \delta ^{-1} \Vert w \Vert ^2_{U_\delta (\varGamma )\cap \varOmega _i} + \delta \Vert \nabla w \Vert ^2_{U_\delta (\varGamma )\cap \varOmega _i } \qquad w \in H^1(U_\delta (\varGamma ) \cap \varOmega _i) \end{aligned}$$

(152)

with $\delta \sim h$ we get

$$\begin{aligned} II_i= & {} h \Vert \nabla \rho _i \Vert ^2_{\varGamma } + h^{-1} \Vert \rho _i \Vert ^2_{\varGamma } \end{aligned}$$

(153)

$$\begin{aligned}\lesssim & {} \delta ^{-1} h \Vert \nabla \rho _i \Vert ^2_{U_\delta (\varGamma )\cap \varOmega _i} + \delta h \Vert \nabla ^2 \rho _i \Vert ^2_{U_\delta (\varGamma )\cap \varOmega _i} \end{aligned}$$

(154)

$$\begin{aligned}&+\, \delta ^{-1} h^{-1} \Vert \rho _i \Vert ^2_{U_\delta (\varGamma )\cap \varOmega _i} + \delta h^{-1} \Vert \nabla \rho _i \Vert ^2_{U_\delta (\varGamma )\cap \varOmega _i} \end{aligned}$$

(155)

$$\begin{aligned}\lesssim & {} \delta ^{-1} h^{-1}\Vert \rho _i \Vert ^2_{\varOmega _i} + (\delta h^{-1} + \delta ^{-1} h ) \Vert \nabla \rho _i \Vert ^2_{\varOmega _i} + \delta h \Vert \nabla ^2 \rho _i \Vert ^2_{\varOmega _i} \end{aligned}$$

(156)

$$\begin{aligned}\lesssim & {} h^2 \Vert v_i \Vert ^2_{H^2(\varOmega _i)} \end{aligned}$$

(157)

where we at last used (131). For Term IV, applying (137) yields,

$$\begin{aligned} IV = h \Vert \nabla \rho _i \Vert ^2_{\varGamma } + h^{-1} \Vert \rho _i \Vert ^2_{\varGamma }&\lesssim \sum _{m=0}^2 h^{2(m-1)} \Vert \nabla ^m \rho _i \Vert ^2_{\varOmega _i} \end{aligned}$$

(158)

Term $\textit{III}_i$ is estimate using the element wise trace inequality

$$\begin{aligned} \Vert w \Vert _F^2 \lesssim h^{-1} \Vert w \Vert ^2_T + h \Vert \nabla w \Vert ^2_T \end{aligned}$$

(159)

which gives

$$\begin{aligned} \textit{III}_i= & {} \Vert \rho _i \Vert ^2_{s_{h,i}} \lesssim \sum _{F\in \mathcal {F}_{h,i}} h \Vert [ n\cdot \nabla \rho _i]\Vert ^2_{F} \end{aligned}$$

(160)

$$\begin{aligned}\lesssim & {} \sum _{F\in \mathcal {F}_{h,i}} \Vert \nabla \rho _i \Vert ^2_{\mathcal {T}_{h,i}(F)} + h \Vert \nabla ^2 \rho _i \Vert ^2_{\mathcal {T}_{h,i}(F)} \lesssim \Vert v \Vert ^2_{H^m(O_{h,i})} \end{aligned}$$

(161)

where $\mathcal {T}_{h,i}(F)$ is the set of two elements that share the face F. In a similar way show that

$$\begin{aligned} V = \Vert \rho _\varGamma \Vert ^2_{s_{h,\varGamma }} \lesssim \Vert v \Vert ^2_{H^m(O_{h,i})} \end{aligned}$$

(162)

see [7, 28] for details.

4.3 Continuity and coercivity

We utilize the bounds on the stabilization parameter $\tau$ provided by Lemma 4 to prove continuity and coercivity of the form $A_h$.

Lemma 7

There is a constant independent of the eigenvalues of B, such that for all $v,w\in \tilde{V} + V_h$,

$$\begin{aligned} \boxed { \mathcal {A}^R_h (v, w ) \lesssim |||v \!|||_h |||w |||_h } \end{aligned}$$

(163)

There is a constant independent of the eigenvalues of B, such that for all $v \in V_h$,

$$\begin{aligned} \boxed { |||v |||_h^2 \lesssim \mathcal {A}^R_h (v, v) } \end{aligned}$$

(164)

Proof

(163). Starting from the definition (119), expanding the terms in $\mathcal {A}^R$, and using Cauchy–Schwarz we obtain

$$\begin{aligned} \mathcal {A}^R(v,w)= & {} \sum _{i=1}^2 (A_{i}\nabla v_i , \nabla w_i)_{\varOmega _i} +(A_\varGamma \nabla _\varGamma v_\varGamma , \nabla _\varGamma w_\varGamma )_\varGamma \end{aligned}$$

(165)

$$\begin{aligned}&+\, ( (n\cdot A \nabla v), (B^{-1} \tau B^{-1} - B^{-1}) ( n\cdot A \nabla w ) )_\varGamma \end{aligned}$$

(166)

$$\begin{aligned}&+\, ( (n\cdot A \nabla v), (B^{-1} \tau - I) ( w - w_\varGamma ) )_\varGamma \end{aligned}$$

(167)

$$\begin{aligned}&+\, ( (n\cdot A \nabla w), (B^{-1} \tau - I) ( v - v_\varGamma ) )_\varGamma \end{aligned}$$

(168)

$$\begin{aligned}&+\, (( v - v_\varGamma ), \tau ( w - w_\varGamma ))_\varGamma \end{aligned}$$

(169)

$$\begin{aligned}\le & {} \sum _{i=1}^2\Vert \nabla v_i \Vert _{A_i,\varOmega _i} \Vert \nabla w_i\Vert _{A_i,\varOmega _i} +\Vert \nabla _\varGamma v_\varGamma \Vert _{A_\varGamma ,\varGamma } \Vert \nabla _\varGamma w_\varGamma \Vert _{A_i,\varGamma } \end{aligned}$$

(170)

$$\begin{aligned}&+\, \Vert n\cdot A \nabla v\Vert _\varGamma \Vert B^{-1} \tau B^{-1} - B^{-1}\Vert _{L^\infty (\varGamma )} \Vert n\cdot A \nabla w \Vert _{\varGamma } \end{aligned}$$

(171)

$$\begin{aligned}&+\, \Vert n\cdot A \nabla v\Vert _\varGamma \Vert ( B^{-1} \tau - I) \tau ^{-1/2}\Vert _{L^\infty (\varGamma )} \Vert w - w_\varGamma \Vert _{\tau ,\varGamma } \end{aligned}$$

(172)

$$\begin{aligned}&+\, \Vert n\cdot A \nabla w\Vert _\varGamma \Vert (B^{-1} \tau - I)\tau ^{-1/2} \Vert _{L^\infty (\varGamma )} \Vert v - v_\varGamma \Vert _{\tau ,\varGamma } \end{aligned}$$

(173)

$$\begin{aligned}&+\, \Vert v - v_\varGamma \Vert _{\tau ,\varGamma } \Vert w - w\Vert _{\tau ,\varGamma } \end{aligned}$$

(174)

$$\begin{aligned}= & {} \bigstar \end{aligned}$$

(175)

Using the estimates (105)–(106) we obtain

$$\begin{aligned} \bigstar\le & {} \sum _{i=1}^2 \Vert \nabla v_i\Vert _{A_i,\varOmega _i} \Vert \nabla w_i\Vert _{A_i,\varOmega _i} +\Vert \nabla _\varGamma v_\varGamma \Vert _{A_\varGamma ,\varGamma } \Vert \nabla _\varGamma w_\varGamma \Vert _{A_\varGamma ,\varGamma } \end{aligned}$$

(176)

$$\begin{aligned}&+\, \beta ^{-1} h \Vert n\cdot A \nabla v\Vert _\varGamma \Vert n\cdot A \nabla w \Vert _{\varGamma } \end{aligned}$$

(177)

$$\begin{aligned}&+\, \beta ^{-1/2} h^{1/2} \Vert n\cdot A \nabla v\Vert _\varGamma \Vert w - w_\varGamma \Vert _{\tau ,\varGamma } \end{aligned}$$

(178)

$$\begin{aligned}&+\, \beta ^{-1/2} h^{1/2} \Vert n\cdot A \nabla w\Vert _\varGamma \Vert v - v_\varGamma \Vert _{\tau ,\varGamma } \end{aligned}$$

(179)

$$\begin{aligned}&+\,\Vert v - v_\varGamma \Vert _{\tau ,\varGamma } \Vert w - w\Vert _{\tau ,\varGamma } \end{aligned}$$

(180)

$$\begin{aligned}\le & {} \sum _{i=1}^2 \Vert \nabla v_i\Vert _{A_i,\varOmega _i} \Vert \nabla w_i\Vert _{A_i,\varOmega _i} +\Vert \nabla _\varGamma v_\varGamma \Vert _{A_\varGamma ,\varGamma } \Vert \nabla _\varGamma w_\varGamma \Vert _{A_\varGamma ,\varGamma } \end{aligned}$$

(181)

$$\begin{aligned}&+\, (\beta ^{-1} \Vert n\Vert ^2_{A,\infty ,\varGamma }) h^{1/2} \Vert \nabla v\Vert _{A,\varGamma } h^{1/2} \Vert \nabla w \Vert _{\varGamma } \end{aligned}$$

(182)

$$\begin{aligned}&+\, (\beta ^{-1} \Vert n\Vert ^2_{A,\infty ,\varGamma } )^{1/2} h^{1/2} \Vert \nabla v\Vert _{A,\varGamma } \Vert w - w_\varGamma \Vert _{\tau ,\varGamma } \end{aligned}$$

(183)

$$\begin{aligned}&+\, (\beta ^{-1} \Vert n\Vert ^2_{A,\infty ,\varGamma } )^{1/2} h^{1/2} \Vert \nabla w\Vert _{A,\varGamma } \Vert v - v_\varGamma \Vert _{\tau ,\varGamma } \end{aligned}$$

(184)

$$\begin{aligned}&+\,\Vert v - v_\varGamma \Vert _{\tau ,\varGamma } \Vert w - w\Vert _{\tau ,\varGamma } \end{aligned}$$

(185)

$$\begin{aligned}\lesssim & {} |||v \!|||_h |||w |||_h \end{aligned}$$

(186)

where we used the bound $\beta ^{-1} \Vert n\Vert ^2_{A,\infty ,\varGamma } \lesssim 1$, see (104). By the Cauchy–Schwarz inequality we have $s_h(v,w) \lesssim |||v \!|||_h |||w |||_h$.

(164). To prove the coercivity we have the identity

$$\begin{aligned} \mathcal {A}^R_h(v,v)= & {} \sum _{i=1}^2 (A_{i}\nabla v_i , \nabla v_i)_{\varOmega _i} +(A_\varGamma \nabla _\varGamma v_\varGamma , \nabla _\varGamma v_\varGamma )_\varGamma + s_h(v,v) \end{aligned}$$

(187)

$$\begin{aligned}&+\, ( (n\cdot A \nabla v), (B^{-1} \tau B^{-1} - B^{-1}) ( n\cdot A \nabla v ) )_\varGamma \end{aligned}$$

(188)

$$\begin{aligned}&+\, 2( (n\cdot A \nabla v), (B^{-1} \tau - I) ( v - v_\varGamma ) )_\varGamma \end{aligned}$$

(189)

$$\begin{aligned}&+\, (( v - v_\varGamma ), \tau ( v - v_\varGamma ))_\varGamma \end{aligned}$$

(190)

$$\begin{aligned}\ge & {} \sum _{i=1}^2 \Vert \nabla v_i\Vert _{A_i,\varOmega _i}^2 +\Vert \nabla _\varGamma v_\varGamma \Vert _{A_\varGamma ,\varGamma }^2 + \Vert v\Vert ^2_{s_h} \end{aligned}$$

(191)

$$\begin{aligned}&-\, \beta ^{-1} \Vert n\Vert ^2_{A,\infty ,\varGamma } h \Vert \nabla v\Vert ^2_{A,\varGamma } \end{aligned}$$

(192)

$$\begin{aligned}&-\, 2 (\beta ^{-1} \Vert n\Vert _{A,\infty ,\varGamma }^2)^{1/2} h^{1/2} \Vert \nabla v\Vert _{A,\varGamma } \Vert v - v_\varGamma \Vert _{\tau ,\varGamma } \end{aligned}$$

(193)

$$\begin{aligned}&+\, \Vert v - v_\varGamma \Vert ^2_{\tau ,\varGamma } \end{aligned}$$

(194)

We conclude the argument as usual by estimating the negative terms as follows

$$\begin{aligned}&\beta ^{-1} \Vert n\Vert ^2_{A,\infty ,\varGamma } h \Vert \nabla v\Vert ^2_{A,\varGamma } + 2 (\beta ^{-1} \Vert n\Vert _{A,\infty ,\varGamma }^2)^{1/2} h^{1/2} \Vert \nabla v\Vert _{A,\varGamma } \Vert v - v_\varGamma \Vert _{\tau ,\varGamma } \end{aligned}$$

(195)

$$\begin{aligned}&\qquad \le 3 \beta ^{-1} \Vert n\Vert ^2_{A,\infty ,\varGamma } h \Vert \nabla v\Vert ^2_{A,\varGamma } + \frac{1}{2} \Vert v - v_\varGamma \Vert ^2_{\tau ,\varGamma } \end{aligned}$$

(196)

$$\begin{aligned}&\qquad \le 3 \beta ^{-1} \Vert n\Vert ^2_{A,\infty ,\varGamma } C_I\left( \sum _{i=1}^2 \Vert \nabla v_i\Vert ^2_{A_i, \varOmega _i} + \Vert v \Vert ^2_{s_{h,i}} \right) + \frac{1}{2} \Vert v - v_\varGamma \Vert ^2_{\tau ,\varGamma } \end{aligned}$$

(197)

$$\begin{aligned}&\qquad \le \frac{1}{2} \left( \sum _{i=1}^2 \Vert \nabla v_i\Vert ^2_{A_i, \varOmega _i} + \Vert v \Vert ^2_{s_{h,i}} \right) + \frac{1}{2} \Vert v - v_\varGamma \Vert ^2_{\tau ,\varGamma } \end{aligned}$$

(198)

Here we used the inverse estimate

$$\begin{aligned} h \Vert \nabla v_i \Vert ^2_{A_i,\varGamma } \le C_I( \Vert \nabla v_i\Vert ^2_{A_i, \varOmega _i} + \Vert v \Vert ^2_{s_{h,i}} ) \end{aligned}$$

(199)

which follows from the inverse bound

$$\begin{aligned} h \Vert \nabla v_i \Vert ^2_{A_i,\varGamma \cap T} \lesssim h \Vert \nabla v_i \Vert ^2_{\varGamma \cap T} \lesssim \Vert \nabla v_i \Vert ^2_{T} \lesssim \Vert \nabla v_i \Vert ^2_{A_i,T} \end{aligned}$$

(200)

together with (90), and finally, we chose $\beta$ large enough to guarantee that

$$\begin{aligned} 3 \beta ^{-1} \Vert n\Vert ^2_{A,\infty ,\varGamma } C_I \le \frac{1}{2} \end{aligned}$$

(201)

We conclude that

$$\begin{aligned} \mathcal {A}^R_h(v,v) \ge \frac{1}{2} |||v |||_h^2 \end{aligned}$$

(202)

which completes the proof.

4.4 A priori error estimates

In this section we prove error estimates for the approximate solution $u_h$.

Theorem 1

Let $u \in \tilde{V}$ be the solution of (36) and $u_h \in V_h$ be the solution of (118). Then there is a constant not dependent on the matrix B in the interface condition (3) such that

$$\begin{aligned} \boxed {|||u - u_h |||_h \lesssim h \left( \sum _{i=1}^2 \Vert u_i \Vert _{H^2(\varOmega _i)} + \Vert u_\varGamma \Vert _{H^2(\varGamma )} \right) } \end{aligned}$$

(203)

Under the assumption that B is a constant positive definite matrix and $A_i = a_i \mathbb {I}_{[d \times d]}$ for $i=1,2,$ and $A_\varGamma = a_i \mathbb {I}_{[(d-1) \times (d-1)]}$ we also have the estimate

$$\begin{aligned} \boxed {|||u - u_h |||_h \lesssim h \left( \sum _{i=1}^2 \Vert f_i \Vert _{L^2(\varOmega _i)} + \Vert f_\varGamma \Vert _{L^2(\varGamma )} \right) } \end{aligned}$$

(204)

with constant independent of B.

Proof

First we decompose the error in the approximation error and the discrete error $u-u_h = u - \pi _h u + \pi _h u - u_h$ and note that by the triangle inequality

$$\begin{aligned} |||u - u_h |||_h \lesssim |||u - \pi _h u |||_h+ |||\pi _h u - u_h |||_h. \end{aligned}$$

(205)

The first term on the right hand side is bounded by (146). For the second term on the right hand side, using coercivity (164), Galerkin orthogonality (125), and continuity (163) we obtain

$$\begin{aligned} |||\pi _h u - u_h |||_h^2\lesssim & {} \mathcal {A}^R_h(\pi _h u - u_h, \pi _h u - u_h ) \end{aligned}$$

(206)

$$\begin{aligned}= & {} \mathcal {A}^R(\pi _h u-u, \pi _h u - u_h ) + s_h(\pi _h u, \pi _h u - u_h ) \end{aligned}$$

(207)

$$\begin{aligned}\lesssim & {} |||\pi _h u-u |||_h |||\pi _h u - u_h |||_h. \end{aligned}$$

(208)

In the last inequality we used that if $u^e = (E u_1, E u_2, E_\varGamma u_\varGamma ) \in \tilde{V}$ then

$$\begin{aligned} s_h(\pi _h u, \pi _h u - u_h )=s_h(\pi _h u - u^e, \pi _h u - u_h ) \lesssim |||\pi _h u-u |||_h |||\pi _h u - u_h |||_h. \end{aligned}$$

(209)

Thus

$$\begin{aligned} |||u -u_h |||_h \lesssim |||u - \pi _h u |||_h \lesssim h \left( \sum _{i=1}^2 \Vert u_i \Vert _{H^2(\varOmega _i)} + \Vert u_\varGamma \Vert _{H^2(\varGamma )} \right) \end{aligned}$$

(210)

where we used the interpolation error estimate (146). To conclude we apply the regularity estimate (42).

Corollary 1 Under the same assumptions as for Theorem 1 there holds

$$\begin{aligned} s_h(u_h,u_h) \lesssim h \left( \sum _{i=1}^2 \Vert u_i \Vert _{H^2(\varOmega _i)} + \Vert u_\varGamma \Vert _{H^2(\varGamma )} \right) . \end{aligned}$$

(211)

Proof

Using the triangle inequality we see that

$$\begin{aligned} \Vert u_h\Vert _{s_h} \le \Vert \pi _h u\Vert _{s_h}+\Vert \pi _h u - u_h\Vert _{s_h} \end{aligned}$$

(212)

The second term on the right hand side is bounded by the arguments of Theorem 1. For the first term on the right hand side recall that by the consistency properties of $s_h$ and the construction of $\pi _h u$ there holds

$$\begin{aligned} s_h(\pi _h u,\pi _h u) =s_h(u^e-\pi _h u, u^e - \pi _h u). \end{aligned}$$

(213)

We conclude the proof by applying (161) and (162).

The following error estimate in the $L^2$-norm also holds.

Theorem 2

Let $u \in \tilde{V}$ be the solution of (36) and $u_h \in V_h$ be the solution of (118). Assume that B is a constant positive definite matrix and $A_i = a_i \mathbb {I}_{[d \times d]}$ for $i=1,2,$ and $A_\varGamma = a_i \mathbb {I}_{[(d-1) \times (d-1)]}$. Then there holds

$$\begin{aligned} \boxed { \Vert u_h - u\Vert _{\varOmega } + \Vert u_{h,\varGamma } - u_\varGamma \Vert _\varGamma \lesssim h^2 \left( \sum _{i=1}^2 \Vert f_i\Vert _{L^2(\varOmega _i)} + \Vert f_\varGamma \Vert _{L^2(\varGamma )}\right) } \end{aligned}$$

(214)

Proof

For $\psi _{i} \in L^2(\varOmega _i)$ and $\psi _\varGamma \in L^2(\varGamma )$, let $\varphi =(\varphi _1,\varphi _2, \varphi _\varGamma ) \in \tilde{V}$ be the weak solution to

$$\begin{aligned} \mathcal {A}(v,\varphi ) = \sum _{i=1}^2 (\psi _{\varOmega _i},v_i)_{\varOmega _i} + (\psi _\varGamma ,v_\varGamma )_\varGamma \qquad \forall v \in V \end{aligned}$$

(215)

Then using the regularity result (42) we have

$$\begin{aligned} \Vert \varphi _1 \Vert _{H^2(\varOmega _1)} + \Vert \varphi _2 \Vert _{H^2(\varOmega _2)} + \Vert \varphi _\varGamma \Vert _{H^2(\varGamma )} \lesssim \Vert \psi _\varOmega \Vert _{\varOmega } + \Vert \psi _\varGamma \Vert _\varGamma \end{aligned}$$

(216)

with constant independent of B. Let $e = (u_1-u_{1,h},u_2-u_{2,h},u_\varGamma -u_{\varGamma ,h})$ be the error and observe that using (124),

$$\begin{aligned}&(\psi _\varOmega ,u_1-u_{1,h})_{\varOmega _1} +(\psi _\varOmega ,u_2-u_{2,h})_{\varOmega _2} + (\psi _\varGamma ,u_\varGamma -u_{h,\varGamma })_\varGamma \end{aligned}$$

(217)

$$\begin{aligned}&\qquad = \mathcal {A}(e,\varphi ) = \mathcal {A}^R(e,\varphi ) \end{aligned}$$

(218)

Applying now the Galerkin orthogonality (125) we see that

$$\begin{aligned} \mathcal {A}^R(e,\varphi )= & {} \mathcal {A}^R(e,\varphi -\pi _h \varphi ) + \mathcal {A}^R(e,\pi _h \varphi ) \end{aligned}$$

(219)

$$\begin{aligned}= & {} \mathcal {A}^R(e,\varphi -\pi _h \varphi ) + s_h(u_h,\pi _h \phi ) \end{aligned}$$

(220)

$$\begin{aligned}= & {} \mathcal {A}^R(e,\varphi -\pi _h \varphi ) + s_h(e,\phi - \pi _h \phi ) \end{aligned}$$

(221)

$$\begin{aligned}= & {} \mathcal {A}^R_h(e,\varphi -\pi _h \varphi ) \end{aligned}$$

(222)

By the continuity (163) we can bound the right hand side,

$$\begin{aligned} \mathcal {A}^R_h(e,\varphi -\pi _h \varphi ) \lesssim |||e |||_h |||\varphi -\pi _h \varphi |||_h \end{aligned}$$

(223)

Then applying the approximation (146) and the regularity (216) we have

$$\begin{aligned} (\psi _\varOmega ,u_1-u_{1,h})_{\varOmega _1} +(\psi _\varOmega ,u_2-u_{2,h})_{\varOmega _2} + (\psi _\varGamma ,u_\varGamma -u_{h,\varGamma })_\varGamma \lesssim h |||e |||_h \end{aligned}$$

(224)

We conclude by applying Theorem 1 in the right hand side and taking the supremum over the functions $(\psi _1,\psi _2,\psi _\varGamma )$ such that $\sum _{i=1}^2 \Vert \psi _i \Vert ^2_{\varOmega _i} + \Vert \psi _\varGamma \Vert ^2_\varGamma =1$.

5 Numerical examples

In this section we illustrate the properties of the model and method by presenting some numerical results. In all examples we used $\beta =10$ as a stabilization parameter.

5.1 Convergence and robustness with respect to conditioning

We consider a simple example with known exact solution: the domain $(0,1)\times (0,1)$ is cut in half along a vertical line at $x=1/2$. We take $A_1 = A_2 = A_{\varGamma } = I$ and choose a problem with exact solution $u=x(1-x)y(1-y)$. This solution corresponds (without coupling) to the source terms

$$\begin{aligned} f_i=2x(1-x)+2y(1-y)\quad \quad i=1,2 \end{aligned}$$

(225)

Since the normal derivative of the exact solution is zero at $x=1/2$, it does not contribute to the source term on the interface. We choose $f_{\varGamma } = 1/2$ corresponding to $u_{\varGamma } =y(1-y)/4$, and thus $u_{\varGamma }=u$ at $x=1/2$. We apply zero Dirichlet boundary conditions on u and on $u_{\varGamma }$ (imposed on the boundary of the band of elements intersected by (1/2, y)). This is now the solution of (1)–(4) independent of B. A sample discrete solution is shown in Fig. 1 with $u_{\varGamma }$ shown as a red line. We did not impose gradient jumps on the band (second term in $s_{h,\varGamma }$), normal stabilization proved sufficient in this case.

Fig. 1

Elevation of the computed solution on a particular mesh (for $\alpha =1$, $\xi =1$)

Fig. 2

Convergence in $L_2(\varOmega )$ and in $L_2(\varGamma )$ for varying $\alpha$ with $\xi =1$. Dashed line has inclination 1:2

Fig. 3

Convergence in $L_2(\varOmega )$ and in $L_2(\varGamma )$ for varying $\xi$ with $\alpha =1$. Dashed line has inclination 1:2

Fig. 4

Convergence in $H^1(\varOmega )$ and in $H^1(\varGamma )$ for varying $\alpha$ with $\xi =1$. Dashed line has inclination 1:1

Fig. 5

Convergence in $H^1(\varOmega )$ and in $H^1(\varGamma )$ for varying $\xi$ with $\alpha =1$. Dashed line has inclination 1:1

Fig. 6

Left: Condition number as a function of meshsize. Dashed line has inclination 1:2. Right: condition numbers on a fixed mesh with varying $\alpha$ using the robust method (118) and the non-robust method (83)

Fig. 7

Elevation of the solution on $\varOmega$ and the band containing $\varGamma$ for $\gamma = 0$

Fig. 8

Elevation of the solution on $\varOmega$ and the band containing $\varGamma$ for $\gamma = 10^{-2}$

Fig. 9

Elevation of the solution on $\varOmega$ and the band containing $\varGamma$ for $\gamma = 1$

In Figs. 2, 3, 4 and 5 we show convergence for different choices of parameters in different norms. The method is completely robust with optimal convergence for all choices. In Fig. 6 we show the variation of the condition number (left) with respect to mesh refinement and choice of $\alpha$. The condition number is $O(h^{-2})$ as expected and does not grow with $\alpha$. We also show (right) the effect of using the non-robust method (83) which shows a linear dependece on $\alpha$ on a fixed mesh, while no such effect is present in the robust method. This robustness is important since $\alpha$ physically depends on the crack width [29] which is expected to be small.

5.2 Effect of gradient jump stabilization

This example is taken from [29] with domain is $(0,2)\times (0,1)$ with Dirichlet data $u=1$ at $x=2$ and $u=0$ at $x=0$. Homogeneous Neumann data were applied at $y=0$ and $y=1$. Data were $f_i=f_\varGamma =0$, $A_1 = A_2= I$ and $A_{\varGamma } = a_{\varGamma }d\, I$ with $a_{\varGamma }=2\times 10^{-3}$ for $1/4< y< 3/4$, $a_{\varGamma }=1$ elsewhere, and with $t=0.01$ (the thickness of the crack). Following [29] we then set $\alpha = 2a_{\varGamma }/d$.

Fig. 10

Computational mesh with interface indicated

Fig. 11

Elevation for $d=10^{-2}$

Fig. 12

Elevation for $d=10^{-3}$

Fig. 13

Elevation for $d=10^{-4}$

To show the effect of stabilization, we chose to scale $s_{h,i}$ and $s_{h,\varGamma }$ by a parameter $\gamma$. We retained $\beta = 10$ and normal stabilization on the band. In Figs. 7, 8 and 9 we show the effect of the parameter $\gamma$. When $\gamma = 0$ the jump in diffusion on the interface leads to slight instabilities at $y=1/4$ and $y=3/4$ which are visible to the eye. These are less pronounced for $\gamma =10^{-2}$ and not significant for $\gamma =1$. The overall solution agrees with that of [29].

5.3 Physical effect of crack width

Finally, we show the effect of the crack width with respect to the solution. We used a domain $(0,1)\times (0,1)$ with a quarter circle crack, shown on the computational mesh in Fig. 10. The data were $A_1 = 5\, I$ (inside the circle) $A_2= I$ (outside the circle) and $a_{\varGamma }=0.1$ with definitions as in Example 5.2. Dirichlet data $u=1$ at $x=1$ and $u=0$ at $x=0$ were used (also on the band) and homogeneous Neumann data on the remaining boundaries. In Figs. 11, 12 and 13 we see the effect of decreasing the interface width by one order of magnitude between figures. The solution rapidly tends to a continuous state.

Appendices

A Proof of Lemma 2

We recall that $\omega \subset \mathbb {R}^d$ is a domain with $C^2$ boundary $\partial \omega$ and we shall show that for $w \in [H^2(\omega )]^d$ it holds

$$\begin{aligned} (\nabla \cdot w, \nabla \cdot w)_\omega= & {} (w\otimes \nabla ,\nabla \otimes w)_\omega + 2(w_n, \text {div}_T w_T)_{\partial \omega } \end{aligned}$$

(226)

$$\begin{aligned}&+(w_T,w_T)_{\kappa ,\partial \omega } + (w_n,w_n)_{\text {tr}(\kappa ),\partial \omega } \end{aligned}$$

(227)

Here $w = w_T + w_n n$ is the decomposition of the vector field w into the tangential and normal components in an open neighborhood $U_\delta (\partial \omega )$ consisting of points in $\omega$ with distance less than some $\delta >0$ of the boundary and the normal field on $\partial \omega$ is extended to $U_\delta (\partial \omega )$ by composition with the closest point mapping associated with $\partial \omega$. Furthermore, $\kappa = \nabla \otimes n = \nabla ^2 \gamma$ is the tangential curvature tensor of $\partial \omega$, $\gamma$ is the signed distance function associated with $\partial \omega$ such that $n = \nabla \gamma$. The surface divergence $\text {div}_T(w)$ on $\partial \omega$ is defined by $\text {div}_T(w) = \text {tr}(w \otimes \nabla _T)$.

Proof

Using partial integration we obtain

$$\begin{aligned}&(\nabla \cdot w, \nabla \cdot w)_\omega = -(w, \nabla (\nabla \cdot v) )_\omega + (w_n, \nabla \cdot w)_{\partial \omega } \end{aligned}$$

(228)

$$\begin{aligned}&\quad = -(w,(\nabla \otimes w)\cdot \nabla )_\omega + ( w_n, \nabla \cdot w)_{\partial \omega } \end{aligned}$$

(229)

$$\begin{aligned}&\quad = (w\otimes \nabla ,\nabla \otimes w)_\omega \underbrace{-(w,(\nabla \otimes w)\cdot n )_{\partial \omega } + (w_n, \nabla \cdot w)_{\partial \omega }}_{\bigstar } \end{aligned}$$

(230)

where we employed the identity

$$\begin{aligned} \nabla (\nabla \cdot w ) = (\nabla \otimes w )\cdot \nabla \end{aligned}$$

(231)

To evaluate the boundary contribution $\bigstar$ we first note that using the product rule

$$\begin{aligned} \nabla w_n = \nabla (w \cdot n ) = (\nabla \otimes w)\cdot n + (\nabla \otimes n) \cdot w = (\nabla \otimes w)\cdot n + \kappa \cdot w \end{aligned}$$

(232)

and thus

$$\begin{aligned} (\nabla \otimes w)\cdot n = \nabla w_n - \kappa \cdot w \end{aligned}$$

(233)

We also note that if $w=w_T$ is tangential we have $w_n = 0$ and

$$\begin{aligned} (\nabla \otimes w_T)\cdot n = - \kappa \cdot w_T \end{aligned}$$

(234)

which also imply

$$\begin{aligned} n \cdot (\nabla \otimes w_T)\cdot n = (\nabla _n w_T)\cdot n= 0 \end{aligned}$$

(235)

Using (233) the boundary contribution takes the form

$$\begin{aligned} \bigstar= & {} -(w,(\nabla \otimes w)\cdot n )_{\partial \omega } + (w_n, \nabla \cdot w)_{\partial \omega } \end{aligned}$$

(236)

$$\begin{aligned}= & {} (w, \kappa \cdot w )_{\partial \omega } \underbrace{-( w, \nabla w_n )_{\partial \omega } + (w_n, \nabla \cdot w)_{\partial \omega }}_{\bigstar \bigstar } \end{aligned}$$

(237)

Continuing $\bigstar \bigstar$ and writing $w = w_T + w_n n$ we get

$$\begin{aligned} \bigstar \bigstar&-( w, \nabla w_n )_{\partial \omega } + (w_n, \nabla \cdot w)_{\partial \omega } \end{aligned}$$

(238)

$$\begin{aligned}&= -( w_n n, \nabla w_n )_{\partial \omega } + (w_n, \nabla \cdot (w_n n) )_{\partial \omega } \end{aligned}$$

(239)

$$\begin{aligned}&\quad -( w_T, \nabla w_n )_{\partial \omega } + (w_n, \nabla \cdot w_T)_{\partial \omega } \end{aligned}$$

(240)

$$\begin{aligned}&= -( w_n n, \nabla w_n )_{\partial \omega } + (w_n, (\nabla \cdot w_n)n + w_n \nabla \cdot n )_{\partial \omega } \end{aligned}$$

(241)

$$\begin{aligned}&\quad -( w_T, \nabla w_n )_{\partial \omega } + (w_n, \nabla \cdot w_T)_{\partial \omega } \end{aligned}$$

(242)

$$\begin{aligned}&= (w_n,\text {tr}(\kappa ) w_n)_{\partial \omega } \underbrace{ -( w_T, \nabla w_n )_{\partial \omega } + (w_n, \nabla \cdot w_T)_{\partial \omega } }_{\bigstar \bigstar \bigstar } \end{aligned}$$

(243)

where we used the product rule $\nabla \cdot (w_n n ) = (\nabla w_n) n + w_n \nabla \cdot n$, and the identity $\nabla \cdot n = \text {tr}(n \otimes \nabla ) = \text {tr}(\kappa )$. Finally, for the last term $\bigstar \bigstar \bigstar$, we split the gradient $\nabla w$ in its tangential and normal components $\nabla w = \nabla _T w + n \nabla _n w$ and starting from the definition of the divergence as the trace of the derivative of w we get

$$\begin{aligned} \nabla \cdot w_T= & {} \text {tr}(w_T \otimes \nabla ) \end{aligned}$$

(244)

$$\begin{aligned}= & {} \text {tr}(w_T \otimes \nabla _T) + \text {tr}(w_T \otimes n \nabla _n) \end{aligned}$$

(245)

$$\begin{aligned}= & {} \text {div}_T w_T + \underbrace{n \cdot (w_T \otimes \nabla )\cdot n}_{=0} \end{aligned}$$

(246)

$$\begin{aligned}= & {} \text {div}_T w_T \end{aligned}$$

(247)

where $\text {div}_T$ is the tangential divergence and we used (235). Finally, using Green’s formula

$$\begin{aligned} -(\text {div}_T w_T, v)_{\partial \omega } = (w_T, \nabla _T v)_{\partial \omega } \end{aligned}$$

(248)

on the closed boundary $\partial \omega$ we obtain

$$\begin{aligned} \bigstar \bigstar \bigstar= & {} -( w_T, \nabla _T w_n )_{\partial \omega } + (w_n, \nabla \cdot w_T)_{\partial \omega } \end{aligned}$$

(249)

$$\begin{aligned}= & {} 2( \text {div}_T ( w_T ) , w_n )_{\partial \omega } \end{aligned}$$

(250)

which concludes the proof.

B Proof of A trace inequality

We shall prove the trace inequality

$$\begin{aligned} \Vert \nabla _n v \Vert _{H^{1/2}(\varGamma )} \lesssim \Vert \nabla v \Vert _{H^1(\varOmega _i)} \end{aligned}$$

(251)

We begin by introducing some notation. Let $\gamma$ be the signed distance function associated with $\varGamma$, let $U_\delta (\varGamma )$ be the tubular neighborhood consisting of points x with distance $|\gamma (x)|$ less than $\delta$ to $\varGamma$, let p be the closest point mapping $U_\delta (\varGamma ) \rightarrow \varGamma$, and recall that since $\varGamma$ is smooth there is $\delta _0>0$ such that p is a well defined function on $U_{\delta _0}(\varGamma )$. Let $\psi _{\delta }:\varOmega _i \rightarrow [0,1]$ be a non negative smooth cut off function with support in $U_{\delta }(\varGamma ) \cap \varOmega _i$, such that $\psi = 1$ on $\varGamma$ and $\Vert \psi \Vert _{L^\infty (\varOmega _i)} \le C \delta ^{-1}$. Let $n^e = n \circ p$ be the extension of n to $U_{\delta _0}(\varGamma )$ and using the definition $\Vert v \Vert _{H^{1/2}(\varGamma )} = \inf _{w \in H^1(\varOmega _i), w = v {\mathrm{on}} \varGamma } \Vert w \Vert _{H^1(\varOmega _i)}$ of the $H^{1/2}(\varGamma )$ norm we conclude that taking $w =\psi _\delta n^e \cdot \nabla u_i$ we have the inequality

$$\begin{aligned} \Vert n \cdot \nabla u_i \Vert ^2_{H^1/2(\varGamma )} \le \Vert \psi _\delta n^e \cdot \nabla u_i \Vert ^2_{\varOmega _i} + \Vert \nabla (\psi _\delta n^e \cdot \nabla u_i ) \Vert ^2_{\varOmega _i} \end{aligned}$$

(252)

Since $\Vert \psi _\delta \Vert _{L^\infty (\varOmega _i)} \le 1$ and $\Vert n^e \Vert _{L^\infty (U_{\delta _0}(\varGamma )) \cap \varOmega _i} = 1$ we have

$$\begin{aligned} \Vert \psi _\delta n^e \cdot \nabla u_i \Vert ^2_{\varOmega _i} \le \Vert \nabla u_i \Vert ^2_{\varOmega _i} \end{aligned}$$

(253)

and using the product rule followed by the triangle inequality and standard estimates

$$\begin{aligned} \Vert \nabla ( \psi _\delta n^e \cdot \nabla u_i ) \Vert ^2_{\varOmega _i}\le & {} \Vert \nabla \psi _\delta \Vert _{L^\infty (U_\delta (\varGamma )\cap \varOmega _i)} \Vert \nabla u_i \Vert ^2_{\varOmega _i} \end{aligned}$$

(254)

$$\begin{aligned}&+\, \Vert \nabla n^e \Vert _{L^\infty (U_\delta (\varGamma )\cap \varOmega _i)} \Vert \nabla u_i \Vert ^2_{\varOmega _i} + \Vert \nabla ^2 u_i \Vert ^2_{\varOmega _i} \end{aligned}$$

(255)

$$\begin{aligned}\le & {} C^2 \delta ^{-2} \Vert \nabla u_i \Vert ^2_{\varOmega _i} + C^2_\kappa \Vert \nabla u_i \Vert ^2_{\varOmega _i} + \Vert \nabla ^2 u_i \Vert ^2_{\varOmega _i} \end{aligned}$$

(256)

$$\begin{aligned}\lesssim & {} \Vert \nabla u_i \Vert ^2_{H^1(\varOmega _i)} \end{aligned}$$

(257)

where C denotes a generic constant and we used the properties of the cut off function. Finally, to compute $D n^e$ we recall that the closest point mapping takes the form $p(x) = x - \gamma (x) n^e(x)$, and using the chain rule together with the fact that $\nabla \gamma = n^e$ we get

$$\begin{aligned} D n^e = D n D p = D n ( I - n^e \otimes n^e - \gamma D n^e ) \end{aligned}$$

(258)

Rearranging terms and using the fact that $Dn = \kappa$ we get

$$\begin{aligned} ( I + \gamma \kappa ) D n^e = \kappa ( I - n^e \otimes n^e) = \kappa \end{aligned}$$

(259)

since $\kappa$ is tangential and therefore

$$\begin{aligned} D n^e = ( I + \gamma \kappa )^{-1} \kappa \end{aligned}$$

(260)

For $\gamma \le \delta$, with $\delta$ small enough we have $\Vert D n^e(x) \Vert \lesssim \Vert \kappa (p(x) \Vert \lesssim 1$, which completes the proof.

C Variable interface coefficient in the regularity estimate

We show that the regularity estimate (42) also holds for a variable coeffcient Bin the interface condition such that $B \in [W^1_\infty (\varGamma )]^{2\times 2}$. We replace identities (54)–(56) with the following estimate

$$\begin{aligned}&( B[u - u_\varGamma ], \varDelta _\varGamma [u - u_\varGamma ] )_{\varGamma } = ((B[u - u_\varGamma ])\otimes \nabla _\varGamma , [u - u_\varGamma ]\otimes \nabla _\varGamma )_{\varGamma } \end{aligned}$$

(261)

$$\begin{aligned}&\quad = ( B ([u - u_\varGamma ] \otimes \nabla _\varGamma ) ,[u - u_\varGamma ]\otimes \nabla _\varGamma )_{\varGamma } \end{aligned}$$

(262)

$$\begin{aligned}&\qquad + ( (\nabla _\varGamma \otimes B)[u - u_\varGamma ], \nabla _\varGamma \otimes [u - u_\varGamma ] )_{\varGamma } \end{aligned}$$

(263)

$$\begin{aligned}&\quad \ge \Vert [ u - u_\varGamma ] \otimes \nabla _\varGamma \Vert ^2_{B,\varGamma } \end{aligned}$$

(264)

$$\begin{aligned}&\qquad - \Vert (\nabla _\varGamma \otimes B) B^{-1}\Vert _{L^\infty (\varGamma )} \Vert [u - u_\varGamma ]\Vert _{B,\varGamma } \Vert [u - u_\varGamma ] \otimes \nabla _\varGamma \Vert _{B,\varGamma } \end{aligned}$$

(265)

$$\begin{aligned}&\quad \ge \frac{1}{2} \Vert [ u - u_\varGamma ] \otimes \nabla _\varGamma \Vert ^2_{B,\varGamma } -\frac{1}{2} \Vert (\nabla _\varGamma \otimes B) B^{-1}\Vert ^2_{L^\infty (\varGamma )} \Vert [u - u_\varGamma ]\Vert _{B,\varGamma }\nonumber \\ \end{aligned}$$

(266)

Here we note that the negative term on the right hand side can be estimated using the energy norm and using the fact that B is positive definite (6) we also have $\Vert (\nabla _\varGamma \otimes B) B^{-1}\Vert _{L^\infty (\varGamma )} \lesssim c_\lambda ^{-1}\Vert (\nabla _\varGamma \otimes B) \Vert _{L^\infty (\varGamma )}$.