Error estimates in weighted Sobolev norms for finite element immersed interface methods

TL;DR

This paper investigates error estimates for immersed interface methods in weighted Sobolev norms, demonstrating that suboptimal convergence near interfaces is a local issue and can be mitigated with a distributionally consistent reformulation.

Contribution

It proves that the poor convergence near interfaces is local and shows how weighted norms and reformulation improve global error estimates in immersed boundary methods.

Findings

Error estimates are valid globally when using weighted norms.

Local deterioration of convergence is confined near the interface.

Reformulating the problem improves approximation accuracy.

Abstract

When solving elliptic partial differential equations in a region containing immersed interfaces (possibly evolving in time), it is often desirable to approximate the problem using an independent background discretisation, not aligned with the interface itself. Optimal convergence rates are possible if the discretisation scheme is enriched by allowing the discrete solution to have jumps aligned with the surface, at the cost of a higher complexity in the implementation. A much simpler way to reformulate immersed interface problems consists in replacing the interface by a singular force field that produces the desired interface conditions, as done in immersed boundary methods. These methods are known to have inferior convergence properties, depending on the global regularity of the solution across the interface, when compared to enriched methods. In this work we prove that this detrimental effect on the convergence properties of the approximate solution is only a local phenomenon, restricted to a small neighbourhood of the interface. In particular we show that optimal approximations can be constructed in a natural and inexpensive way, simply by reformulating the problem in a distributionally consistent way, and by resorting to weighted norms when computing the global error of the approximation.

Tables

Table 1. Table 1: Error in the weighted L α 2 subscript superscript 𝐿 2 𝛼 L^{2}_{\alpha} norm ‖ u − u h ‖ 0 , α , Ω subscript norm 𝑢 subscript 𝑢 ℎ 0 𝛼 Ω \|u-u_{h}\|_{0,\alpha,\Omega} for different values of α 𝛼 \alpha in the two-dimensional case.

#dofs		$\alpha =0$		$\alpha =0.1$		$\alpha =0.2$		$\alpha =0.3$		$\alpha =0.4$		$\alpha =0.499$
$25$		$6$	$.6412\cdot {10}^{-2}$	$5$	$.6298\cdot {10}^{-2}$	$4$	$.7772\cdot {10}^{-2}$	$4$	$.0580\cdot {10}^{-2}$	$3$	$.4513\cdot {10}^{-2}$	$2$	$.9441\cdot {10}^{-2}$
$81$		$1$	$.6534\cdot {10}^{-2}$	$1$	$.3347\cdot {10}^{-2}$	$1$	$.0788\cdot {10}^{-2}$	$8$	$.7317\cdot {10}^{-3}$	$7$	$.0789\cdot {10}^{-3}$	$5$	$.7619\cdot {10}^{-3}$
$289$		$7$	$.1702\cdot {10}^{-3}$	$5$	$.3627\cdot {10}^{-3}$	$4$	$.0175\cdot {10}^{-3}$	$3$	$.0157\cdot {10}^{-3}$	$2$	$.2690\cdot {10}^{-3}$	$1$	$.7170\cdot {10}^{-3}$
$1,089$		$2$	$.6533\cdot {10}^{-3}$	$1$	$.8443\cdot {10}^{-3}$	$1$	$.2841\cdot {10}^{-3}$	$8$	$.9595\cdot {10}^{-4}$	$6$	$.2688\cdot {10}^{-4}$	$4$	$.4186\cdot {10}^{-4}$
$4,225$		$9$	$.4960\cdot {10}^{-4}$	$6$	$.1451\cdot {10}^{-4}$	$3$	$.9853\cdot {10}^{-4}$	$2$	$.5918\cdot {10}^{-4}$	$1$	$.6922\cdot {10}^{-4}$	$1$	$.1157\cdot {10}^{-4}$
$16,641$		$3$	$.0996\cdot {10}^{-4}$	$1$	$.8793\cdot {10}^{-4}$	$1$	$.1417\cdot {10}^{-4}$	$6$	$.9562\cdot {10}^{-5}$	$4$	$.2578\cdot {10}^{-5}$	$2$	$.6384\cdot {10}^{-5}$
$66,049$		$1$	$.1688\cdot {10}^{-4}$	$6$	$.5940\cdot {10}^{-5}$	$3$	$.7273\cdot {10}^{-5}$	$2$	$.1131\cdot {10}^{-5}$	$1$	$.2039\cdot {10}^{-5}$	$6$	$.9595\cdot {10}^{-6}$
$263,169$		$4$	$.1721\cdot {10}^{-5}$	$2$	$.1943\cdot {10}^{-5}$	$1$	$.1562\cdot {10}^{-5}$	$6$	$.1104\cdot {10}^{-6}$	$3$	$.2467\cdot {10}^{-6}$	$1$	$.7548\cdot {10}^{-6}$
$1,050,625$		$1$	$.4844\cdot {10}^{-5}$	$7$	$.2710\cdot {10}^{-6}$	$3$	$.5681\cdot {10}^{-6}$	$1$	$.7561\cdot {10}^{-6}$	$8$	$.6942\cdot {10}^{-7}$	$4$	$.3911\cdot {10}^{-7}$

Table 2. Table 2: Error in the weighted H α 1 subscript superscript 𝐻 1 𝛼 H^{1}_{\alpha} norm ‖ u − u h ‖ 1 , α , Ω subscript norm 𝑢 subscript 𝑢 ℎ 1 𝛼 Ω \|u-u_{h}\|_{1,\alpha,\Omega} for different values of α 𝛼 \alpha in the two-dimensional case.

#dofs		$\alpha =0$		$\alpha =0.1$		$\alpha =0.2$		$\alpha =0.3$		$\alpha =0.4$		$\alpha =0.499$
$25$		$1$	$.0430\cdot {10}^0$	$8$	$.8597\cdot {10}^{-1}$	$7$	$.5523\cdot {10}^{-1}$	$6$	$.4486\cdot {10}^{-1}$	$5$	$.5163\cdot {10}^{-1}$	$4$	$.7357\cdot {10}^{-1}$
$81$		$6$	$.5165\cdot {10}^{-1}$	$5$	$.2694\cdot {10}^{-1}$	$4$	$.2683\cdot {10}^{-1}$	$3$	$.4631\cdot {10}^{-1}$	$2$	$.8152\cdot {10}^{-1}$	$2$	$.2984\cdot {10}^{-1}$
$289$		$5$	$.1529\cdot {10}^{-1}$	$3$	$.8646\cdot {10}^{-1}$	$2$	$.9050\cdot {10}^{-1}$	$2$	$.1895\cdot {10}^{-1}$	$1$	$.6556\cdot {10}^{-1}$	$1$	$.2604\cdot {10}^{-1}$
$1,089$		$3$	$.7053\cdot {10}^{-1}$	$2$	$.5813\cdot {10}^{-1}$	$1$	$.8025\cdot {10}^{-1}$	$1$	$.2628\cdot {10}^{-1}$	$8$	$.8867\cdot {10}^{-2}$	$6$	$.3147\cdot {10}^{-2}$
$4,225$		$2$	$.6994\cdot {10}^{-1}$	$1$	$.7514\cdot {10}^{-1}$	$1$	$.1393\cdot {10}^{-1}$	$7$	$.4397\cdot {10}^{-2}$	$4$	$.8858\cdot {10}^{-2}$	$3$	$.2501\cdot {10}^{-2}$
$16,641$		$1$	$.8301\cdot {10}^{-1}$	$1$	$.1107\cdot {10}^{-1}$	$6$	$.7576\cdot {10}^{-2}$	$4$	$.1278\cdot {10}^{-2}$	$2$	$.5387\cdot {10}^{-2}$	$1$	$.5875\cdot {10}^{-2}$
$66,049$		$1$	$.3325\cdot {10}^{-1}$	$7$	$.5323\cdot {10}^{-2}$	$4$	$.2675\cdot {10}^{-2}$	$2$	$.4274\cdot {10}^{-2}$	$1$	$.3911\cdot {10}^{-2}$	$8$	$.1333\cdot {10}^{-3}$
$263,169$		$9$	$.5253\cdot {10}^{-2}$	$5$	$.0205\cdot {10}^{-2}$	$2$	$.6516\cdot {10}^{-2}$	$1$	$.4057\cdot {10}^{-2}$	$7$	$.5118\cdot {10}^{-3}$	$4$	$.1101\cdot {10}^{-3}$
$1,050,625$		$6$	$.7617\cdot {10}^{-2}$	$3$	$.3196\cdot {10}^{-2}$	$1$	$.6328\cdot {10}^{-2}$	$8$	$.0600\cdot {10}^{-3}$	$4$	$.0121\cdot {10}^{-3}$	$2$	$.0531\cdot {10}^{-3}$

Table 3. Table 3: Error in the weighted L α 2 subscript superscript 𝐿 2 𝛼 L^{2}_{\alpha} norm ‖ u − u h ‖ 0 , α , Ω subscript norm 𝑢 subscript 𝑢 ℎ 0 𝛼 Ω \|u-u_{h}\|_{0,\alpha,\Omega} for different values of α 𝛼 \alpha in the three-dimensional case.

#dofs		$\alpha =0$		$\alpha =0.1$		$\alpha =0.2$		$\alpha =0.3$		$\alpha =0.4$		$\alpha =0.499$
$27$		$6$	$.9612\cdot {10}^{-1}$	$6$	$.2832\cdot {10}^{-1}$	$5$	$.6798\cdot {10}^{-1}$	$5$	$.1427\cdot {10}^{-1}$	$4$	$.6649\cdot {10}^{-1}$	$4$	$.2438\cdot {10}^{-1}$
$125$		$1$	$.8304\cdot {10}^{-1}$	$1$	$.5834\cdot {10}^{-1}$	$1$	$.3717\cdot {10}^{-1}$	$1$	$.1901\cdot {10}^{-1}$	$1$	$.0341\cdot {10}^{-1}$	$9$	$.0134\cdot {10}^{-2}$
$729$		$5$	$.7647\cdot {10}^{-2}$	$4$	$.6752\cdot {10}^{-2}$	$3$	$.7965\cdot {10}^{-2}$	$3$	$.0876\cdot {10}^{-2}$	$2$	$.5154\cdot {10}^{-2}$	$2$	$.0573\cdot {10}^{-2}$
$4,913$		$2$	$.0731\cdot {10}^{-2}$	$1$	$.5755\cdot {10}^{-2}$	$1$	$.1994\cdot {10}^{-2}$	$9$	$.1497\cdot {10}^{-3}$	$6$	$.9966\cdot {10}^{-3}$	$5$	$.3802\cdot {10}^{-3}$
$35,937$		$7$	$.6882\cdot {10}^{-3}$	$5$	$.4277\cdot {10}^{-3}$	$3$	$.8402\cdot {10}^{-3}$	$2$	$.7245\cdot {10}^{-3}$	$1$	$.9398\cdot {10}^{-3}$	$1$	$.3922\cdot {10}^{-3}$
$274,625$		$2$	$.6488\cdot {10}^{-3}$	$1$	$.7444\cdot {10}^{-3}$	$1$	$.1517\cdot {10}^{-3}$	$7$	$.6297\cdot {10}^{-4}$	$5$	$.0788\cdot {10}^{-4}$	$3$	$.4177\cdot {10}^{-4}$

Table 4. Table 4: Error in the weighted H α 1 subscript superscript 𝐻 1 𝛼 H^{1}_{\alpha} norm ‖ u − u h ‖ 1 , α , Ω subscript norm 𝑢 subscript 𝑢 ℎ 1 𝛼 Ω \|u-u_{h}\|_{1,\alpha,\Omega} for different values of α 𝛼 \alpha in the three-dimensional case.

#dofs		$\alpha =0$		$\alpha =0.1$		$\alpha =0.2$		$\alpha =0.3$		$\alpha =0.4$		$\alpha =0.499$
$27$		$5$	$.0202\cdot {10}^0$	$4$	$.5494\cdot {10}^0$	$4$	$.1732\cdot {10}^0$	$3$	$.8379\cdot {10}^0$	$3$	$.5388\cdot {10}^0$	$3$	$.2743\cdot {10}^0$
$125$		$3$	$.1626\cdot {10}^0$	$2$	$.7516\cdot {10}^0$	$2$	$.4020\cdot {10}^0$	$2$	$.1003\cdot {10}^0$	$1$	$.8397\cdot {10}^0$	$1$	$.6165\cdot {10}^0$
$729$		$2$	$.1955\cdot {10}^0$	$1$	$.7900\cdot {10}^0$	$1$	$.4623\cdot {10}^0$	$1$	$.1966\cdot {10}^0$	$9$	$.8120\cdot {10}^{-1}$	$8$	$.0792\cdot {10}^{-1}$
$4,913$		$1$	$.6074\cdot {10}^0$	$1$	$.2250\cdot {10}^0$	$9$	$.3566\cdot {10}^{-1}$	$7$	$.1642\cdot {10}^{-1}$	$5$	$.5018\cdot {10}^{-1}$	$4$	$.2513\cdot {10}^{-1}$
$35,937$		$1$	$.1646\cdot {10}^0$	$8$	$.2468\cdot {10}^{-1}$	$5$	$.8560\cdot {10}^{-1}$	$4$	$.1729\cdot {10}^{-1}$	$2$	$.9870\cdot {10}^{-1}$	$2$	$.1579\cdot {10}^{-1}$
$274,625$		$8$	$.1877\cdot {10}^{-1}$	$5$	$.4059\cdot {10}^{-1}$	$3$	$.5802\cdot {10}^{-1}$	$2$	$.3813\cdot {10}^{-1}$	$1$	$.5936\cdot {10}^{-1}$	$1$	$.0801\cdot {10}^{-1}$

Equations

$$-\Delta p=0$$

$$\left[\left[p_{\nu }\right]\right]=f$$

$$\left[\left[p\right]\right]=0$$

$$\left[\left[a\right]\right]:=a^{+}-a^{-},$$

$$p(x)={\int }_{\Gamma }G(x-y)f(y)\mathrm{d}{\Gamma }_y\forall x\in {\mathbb{R}}^n\setminus \Gamma ,$$

$$G(r):=\{\begin{array}{ll}-\frac{1}{2\pi }\log |r|& \text{ when }n=2,\\ \frac{1}{4\pi |r|}& \text{ when }n=3.\end{array}$$

$$-{\Delta }_{x}G(x-y)=\delta (x-y),\forall x\text{ in }{\mathbb{R}}^n\setminus \{y\},$$

$${\int }_{{\mathbb{R}}^n}\delta (x-y)v(x)\mathrm{d}x:=v(y),\forall v\in \mathcal{D}({\mathbb{R}}^n),\forall y\in {\mathbb{R}}^n,$$

$$\begin{array}{rl}-{\Delta }_{x}p(x)=& -{\Delta }_x{\int }_{\Gamma }G(x-y)f(y)\mathrm{d}{\Gamma }_y\\ =& {\int }_{\Gamma }(-{\Delta }_{x}G(x-y))f(y)\mathrm{d}{\Gamma }_y\\ =& {\int }_{\Gamma }\delta (x-y)f(y)\mathrm{d}{\Gamma }_y\forall x\in {\mathbb{R}}^n\setminus \Gamma ,\end{array}$$

$$(\mathcal{M}f)(x):={\int }_{\Gamma }\delta (x-y)f(y)\mathrm{d}{\Gamma }_y.$$

$$⟨\mathcal{M}f,\phi ⟩={\int }_{{\mathbb{R}}^n}\phi (x){\int }_{\Gamma }\delta (x-y)f(y)\mathrm{d}{\Gamma }_y\mathrm{d}x:={\int }_{\Gamma }\phi (y)f(y)\mathrm{d}{\Gamma }_y\forall \phi \in \mathcal{D}({\mathbb{R}}^n).$$

$$-\Delta u=b$$

$$\left[\left[u_{\nu }\right]\right]=f$$

$$\left[\left[u\right]\right]=0$$

$$u=0$$

$$-\Delta z=b$$

$$z=-p$$

$$-\Delta u=b+\mathcal{M}f$$

$$(\mathcal{M}f)(x):={\int }_{\Gamma }\delta (x-y)f(y)\mathrm{d}{\Gamma }_y,$$

$$u=0$$

$$⟨\mathcal{M}f,\phi ⟩:=$$

$$\begin{array}{rl}{\gamma }^{\mathrm{int}}:& H^{s+\frac{1}{2}}({\Omega }^{-})\to H^s(\Gamma ),\\ {\gamma }^{\mathrm{ext}}:& H^{s+\frac{1}{2}}({\Omega }^{+})\to H^s(\Gamma ),\end{array}$$

$$\Vert v{\Vert }_{s,\Gamma }\le C_T\Vert v{\Vert }_{s+\frac{1}{2},\Omega },\forall v\in H^s(\Omega ),0<s<1.$$

$${}_{H^{-s-\frac{1}{2}}(\Omega )}⟨\mathcal{M}f,v{⟩}_{H_0^{s+\frac{1}{2}}(\Omega )}:={\int }_{\Gamma }fv\mathrm{d}\Gamma \forall v\in H_0^{s+\frac{1}{2}}(\Omega ),s\in \left(0,\frac{1}{2}\right],$$

$${}_{H^{-s-\frac{1}{2}}(\Omega )}⟨\mathcal{M}f,v{⟩}_{H_0^{s+\frac{1}{2}}(\Omega )}:={\int }_{\Gamma }fv\mathrm{d}\Gamma \le M\Vert f{\Vert }_{-s,\Gamma }\Vert v{\Vert }_{s+\frac{1}{2},\Omega }.$$

$$\left|{\int }_{\Omega }{\int }_{\Gamma }\delta (x-y)f(y)\phi (x)\mathrm{d}{\Gamma }_y\mathrm{d}{\Omega }_x\right|=\left|{\int }_{\Gamma }f(y)\phi (y)\mathrm{d}{\Gamma }_y\right|\le \Vert f{\Vert }_{-s,\Gamma }\Vert \phi {\Vert }_{s,\Gamma }.$$

$$⟨\mathcal{M}f,\phi ⟩\le C_T\Vert f{\Vert }_{-s,\Gamma }\Vert \phi {\Vert }_{s+\frac{1}{2},\Omega },$$

$$(\nabla u,\nabla v)=(b,v)+⟨\mathcal{M}f,v⟩$$

$$⟨\mathcal{M}f,v⟩:={\int }_{\Gamma }fv\mathrm{d}\Gamma$$

$$|u{|}_{\frac{3}{2}-s,\Omega }\le C\Vert b+\mathcal{M}f{\Vert }_{-\frac{1}{2}-s,\Omega }\le C(\Vert b{\Vert }_{0,\Omega }+\Vert f{\Vert }_{-s,\Gamma }).$$

$$|u{|}_{k+2,\Omega }\le \Vert b+\mathcal{M}f{\Vert }_{k,\Omega },$$

$$W_h^{\ell }$$

$$\begin{array}{rl}|K\cap \Gamma |& \le C_{0}h,\forall K\in {\Omega }_h,\end{array}$$

$$⟨\nabla u_h,\nabla v_h⟩=(b,v_h)+⟨\mathcal{M}f,v_h⟩$$

Full text

A much simpler way to reformulate immersed interface problems consists in replacing the interface by a singular force field that produces the desired interface conditions, as done in immersed boundary methods. These methods are known to have inferior convergence properties, depending on the global regularity of the solution across the interface, when compared to enriched methods.

In this work we prove that this detrimental effect on the convergence properties of the approximate solution is only a local phenomenon, restricted to a small neighbourhood of the interface. In particular we show that optimal approximations can be constructed in a natural and inexpensive way, simply by reformulating the problem in a distributionally consistent way, and by resorting to weighted norms when computing the global error of the approximation.

Error estimates in weighted Sobolev norms for

finite element immersed interface methods

[sissa]Luca Heltai [wias]Nella Rotundo \cortext[cor]Corresponding author. Tel: +49 30 20372-398.

1 Introduction

Interface problems are ubiquitous in nature, and they often involve changes in topology or complex coupling across the interface itself. Such problems are typically governed by elliptic partial differential equations (PDEs) defined on separate domains and coupled together with interface conditions in the form of jumps in the solution and flux across the interface.

It is a general opinion that reliable numerical solutions to interface problems can be obtained using body fitted meshes (possibly evolving in time), as in the Arbitrary Lagrangian Eulerian (ALE) framework Hirt1974 ; DoneaGiulianiHalleux-1982-a . However, in the presence of topological changes, large domain deformations, or freely moving interfaces, these methods may require re-meshing, or even be impractical to use.

Several alternative approaches exist that reformulate the problem using a fixed background mesh, removing the requirement that the position of the interface be aligned with the mesh. These methods originate from the Immersed Boundary Method (IBM), originally introduced by Peskin in Peskin1972 , to study the blood flow around heart valves (see also Peskin2002 , or the review Mittal2005b ), and evolved into a large variety of methods and algorithms.

We distinguish between two main different families of immersed methods. In the first family the interface conditions are incorporated into the finite difference scheme, by modifying the differential operators, or in the finite element space by enriching locally the basis functions to allow for the correct jump conditions in the gradients or in the solution. The second family leaves the discretisation intact, and reformulates the jump conditions using singular source terms.

Important examples of the first family of methods are given by the Immersed Interface Method (IIM) Leveque1994 and its finite element variant Li1998 , or the eXtended Finite Element Method (X-FEM) Sukumar2000 , that exploits partition of unity principles Melenk1996 (see also HansboHansbo-2002-a ; HansboHansbo-2004-a ). For a comparison between IIM and X-FEM see, for example, Vaughan2006 , while for some details on the finite element formulation of the IIM see Li2003a ; Gong2008 ; Hou2013 ; Mu2013 .

The original Immersed Boundary Method Peskin1972 and its variants belong to the second category. Singular source terms are formally written in terms of the Dirac delta distribution, and their discretisation follow two possible routes: i) the Dirac delta distribution is approximated through a smooth function, or ii) the variational definition of the Dirac distribution is used directly in the Finite Element formulation of the problem. For finite differences, the first solution is the only viable option, even though the use of smooth kernels may excessively smear the singularities, leading to large errors in the approximation Hosseini2014 . In finite elements, instead, both solutions are possible. The methods derived from the Immersed Finite Element Method (IFEM) still use approximations of the Dirac delta distribution through the Reproducing Kernel Particle Method (RKPM) ZhangGerstenberger-2004-Immersed-finite-0 .

Variational formulations of the IBM were introduced in BoffiGastaldi-2003-a ; BoffiGastaldi-2007-Numerical-stability-0 ; BoffiGastaldiHeltaiPeskin-2008-a ; Heltai-2008-a , and later generalised in Heltai2012b , where the need to approximate Dirac delta distributions is removed by exploiting directly the weak formulation. Such formulations allow the solution of PDEs with jumps in the gradients without enriching the finite element space, and without introducing approximations of the Dirac delta distribution.

When IBM-like formulations are used to approximate interface problems, a natural deterioration is observed in the convergence of the approximate solution, which is no longer globally smooth, and cannot be expected to converge optimally to the exact solution (see the results section of BoffiGastaldiHeltaiPeskin-2008-a , or Ramiere2008 ). A formal optimal convergence can be observed in special cases Lai2000 , but in general the global convergence properties of these methods is worse when compared to methods where the interface is captured accurately, either by local enrichment of the finite dimensional space, as in the X-FEM or IIM, or by using interface-fitted meshes, as in the ALE method.

In this work we show that this is only partially true, and that optimal approximations can be constructed also when non-body fitted meshes are used, and when no explicit treatment of the jump conditions are imposed in the solution. This can be achieved in a natural and inexpensive way, simply by reformulating the problem in a distributionally consistent way, and by resorting to weighted norms when computing the global error of the approximation. We show here that the deterioration of the error is a purely local phenomena, restricted to a small neighbourhood of the interface itself. In particular we prove that by using suitable powers of the distance function from the interface as weights in weighted Sobolev norms when computing the errors, optimal error estimates can be attained globally.

Weighted Sobolev spaces Kufner1985 ; Turesson2000 , provide a natural framework for the study of the convergence properties of problems with singular sources (see, for example, Agnelli2014 ). These spaces are commonly used in studying problems with singularities in the domain (for example in axisymmetric domains Belhachmi2006 , or in domains with external cusps Duran2009a ) and when the singularities are caused by degenerate or singular behavior of the coefficients of the differential operator Fabes1982a ; Caffarelli2007 ; Cabre2015 . A particularly useful class of weighted Sobolev spaces is given by those spaces whose weights belong to the so-called Muckenhoupt class $A_{p}$ Muckenhoupt1972 . An extensive approximation theory for weighted Sobolev spaces is presented in Nochetto2016 .

The ideas we present in this work are inspired by the works of D’Angelo and Quarteroni DAngelo2008 and D’Angelo Dangelo2012 , where the authors discuss the coupling between one dimensional source terms and three dimensional diffusion-reaction equations DAngelo2008 , and finite element approximations of elliptic problem with Dirac measures Dangelo2012 .

A general setting for the numerical approximation of elliptic problems with singular sources is available in Drelichman2018 and Otarola2018 .

By applying the same principles, we recover optimal error estimates in the approximation of problems with singular sources distributed along co-dimension one surfaces, such as those arising in the variational formulation of immersed methods.

In Section 2 and 5 we outline the problem we wish to solve, and introduce weighted Sobolev spaces. Section 4 is dedicated to the definition of the numerical approximation, and to the proof of the optimal convergence rates in weighted Sobolev norms. Section 6 presents a numerical validation using both two- and three-dimensional examples, while Section 7 provides some conclusions and perspectives.

2 Model interface problem

When approximating problems with interfaces using non-matching grids, one can choose among several possibilities. For example, one could decide to discretise the differential operators using finite differences, or to pose the problem in a finite dimensional space and leave the differential operators untouched, as in the finite element case.

In both cases, if one wants to enforce strongly the interface conditions, it is necessary to modify the differential operators (as in the IIM method Leveque1994 for finite differences) or to enrich the finite dimensional space (as in the X-FEM Sukumar2000 ). A third option, that we will call distributional approach, consists in leaving the space and the differential operators untouched, and to rewrite the jump conditions in terms of singular sources, as in the original Immersed Boundary Method Peskin2002 .

The distributional approach can be tackled numerically either by mollification of Dirac delta distributions, or by applying variational formulations, where the action of the Dirac distributions is applied using its definition to the finite element test functions.

To fix the ideas, consider a Lipschitz, closed, interface $\Gamma\subset\mathbb{R}^{n}$ of co-dimension one (as in Figure 1), and the following model problem:

Problem 1** (Interface only)**

Given a function $f\in H^{-s}(\Gamma)$, $s\in\left(0,\frac{1}{2}\right]$, find a harmonic function $p$ in $\mathbb{R}^{n}\setminus\Gamma$, such that:

[TABLE]

The notation $\left[\!\left[\cdot\right]\!\right]$ is used to indicate the jump across the interface $\Gamma$, and $p_{\nu}$ indicates the normal derivative of $p$, i.e., $\nu\cdot\nabla p$. The direction $\nu$ on $\Gamma$ is used to define precisely the meaning of $\left[\!\left[a\right]\!\right]$ for any quantity $a$, i.e.:

[TABLE]

where $a^{+}$ lies on the same side of $\nu$.

Problem 1 admits a solution that can be constructed explicitly in terms of the boundary integral representation for harmonic functions (see, e.g., HsiaoWendland-2008-a ):

[TABLE]

where $G$ is the fundamental solution of the Poisson problem in $\mathbb{R}^{n}$:

[TABLE]

The function $G$ satisfies, in the distributional sense,

[TABLE]

where $\delta$ is the $n$-dimensional Dirac delta distribution, i.e., the distribution such that

[TABLE]

and $\mathcal{D}(\mathbb{R}^{n})$ is the space of infinitely differentiable functions, with compact support on $\mathbb{R}^{n}$.

If $f$ is at least $H^{-\frac{1}{2}}(\Gamma)$, then $p$ is globally in $H^{1}_{\text{loc}}(\mathbb{R}^{n})$, it is harmonic in the entire $\mathbb{R}^{n}\setminus\Gamma$ (see, e.g., HsiaoWendland-2008-a for a proof), and we can take its Laplacian in the entire $\mathbb{R}^{n}$ in the sense of distributions.

Exploiting the boundary integral representation (3), the distributional laplacian of $p$

[TABLE]

can be formally expressed in terms of a distributional operator $\mathcal{M}$ as

[TABLE]

Notice that in the definition of the operator $\mathcal{M}f$, the Dirac delta distribution $\delta$ is defined through its action on functions in $\mathcal{D}(\mathbb{R}^{n})$, by the distributional definition (6). In $\mathcal{M}f$, the Dirac distribution is convoluted with $f$ on a domain $\Gamma$ of co-dimension one with respect to $\mathbb{R}^{n}$. The resulting distribution is zero everywhere, and singular across $\Gamma$.

This is usual in the Immersed Boundary and Immersed Interface literature, and should be interpreted as the distributional operator whose effect is to take the trace of the test function on $\Gamma$, and apply the duality product on $\Gamma$ with the function $f$:

[TABLE]

Similarly, let us consider a simply connected and convex domain $\Omega$ of $\mathbb{R}^{n}$, $n=2,3$, with Lipschitz boundary $\partial\Omega$, separated into $\Omega^{+}$ and $\Omega^{-}$ by the surface $\Gamma$, as in Figure 2, where we assume that $\Gamma\cap\partial\Omega=\emptyset$.

We use standard notations for Sobolev spaces (see, for example, Adams2003 ), i.e., $H^{s}(A)=W^{s,2}(A)$, for real $s$, where $L^{2}(A)=H^{0}(A)$ and $H^{1}_{0}(A)$ represents square integrable functions on $A$ with square integrable first derivatives and whose trace is zero on the boundary $\partial A$ of the domain of definition $A$.

Problem 2** (Strong)**

Given $b\in L^{2}(\Omega)$ and $f\in H^{-s}(\Gamma)$, $s\in\left(0,\frac{1}{2}\right]$, find a solution $u$ of the problem

[TABLE]

Using the solution $p$ to Problem 1, a solution to Problem 2 can be written as $u=z+p$ where $z$ satisfies

[TABLE]

Since $\Gamma\cap\partial\Omega=\{\emptyset\}$, and away from $\Gamma$ the solution $p$ is harmonic, its restriction on $\partial\Omega$ is at least Lipschitz and continous, and we conclude that problem (11) admits a unique solution which is at least $H^{2}(\Omega)$, thanks to the assumption that $\Omega$ is convex. The solution $u=z+p$ to Problem 2 is globally $H^{1}(\Omega)$, and at least $H^{2}(\Omega\setminus\Gamma)$.

Problem 2 is equivalent to the following distributional formulation, where the jump conditions are incorporated into a singular right hand side.

Problem 3** (Distributional)**

Given $b\in L^{2}(\Omega)$, $f\in H^{-s}(\Gamma)$, $s\in\left(0,\frac{1}{2}\right]$, find the distribution $u$ such that

[TABLE]

The function $u=z+p$ is then a solution to Problem 2 in $\Omega\setminus\Gamma$, and to Problem 3, in the distributional sense, in the entire domain $\Omega$. The distributional definition of $\mathcal{M}f$ is derived in equations (4–7), and it is given by

[TABLE]

where we exploit that $\mathcal{D}(\Omega)\subset\mathcal{D}(\mathbb{R}^{n})$.

This formulation is at the base of the Immersed Boundary Method Peskin2002 . In the IBM, a problem similar to Problem 3 is discretised by Finite Differences, and the Dirac delta distribution is replaced by a regularised delta function, used to interpolate between non-matching sample points on the co-dimension one surface $\Gamma$ and on the domain $\Omega$.

When regularised Dirac distributions are used, a natural deterioration is observed in the convergence of the approximate solution. A formal second order convergence can be observed in special cases Lai2000 , but in general the global convergence properties of these methods is worse when compared to methods where the interface is taken into account explicitly, like in the IIM Leveque1994 . Moreover, the regularisation itself introduces an additional source of approximation, which smears out the singularity, and may deteriorate further the convergence properties of the method Hosseini2014 .

3 Variational formulation

In order to study the regularity and the well posedeness of Problem 3, we begin by showing that $\mathcal{M}f$ belongs to $H^{-1}(\Omega)$, and therefore there exists a solution $u$ to the distributional Problem 3 which is globally in $H^{1}_{0}(\Omega)$.

We begin by recalling standard results for trace operators. For a bounded domain $\Omega^{+}$ (or $\Omega^{-}$) with (part of the) boundary $\Gamma$, given a function $u\in C^{0}(\overline{\Omega^{+}})$ (or $u\in C^{0}(\overline{\Omega^{-}})$), it makes sense to define the restriction of $u$ on $\Gamma$, simply by considering its pointwise restriction. For Sobolev spaces, we recall the following classical result (see, e.g., (mclean2000strongly, , Theorem 3.38), or ciarlet78 ):

Theorem 1** (Trace theorem)**

Let $\Gamma$ be a Lipschitz closed co-dimension one surface, splitting $\Omega$ into $\Omega^{+}$ and $\Omega^{-}$. For $0<s<1$ the interior and the exterior trace operators

[TABLE]

are bounded, linear, and injective mappings, that posses bounded right inverses. If the function $v$ is globally $H^{s+\frac{1}{2}}(\Omega)$, then $\gamma^{\mathrm{int}}v=\gamma^{\mathrm{ext}}v$, and we omit the symbol altogether, simply indicating with $v$ both the function in $H^{s+\frac{1}{2}}(\Omega)$ and its restriction to $H^{s+\frac{1}{2}}(\Gamma)$. There exists a constant $C_{T}>0$ such that

[TABLE]

Theorem 2** (Regularity of $\mathcal{M}$)**

For $\Gamma$ Lipschitz and $s\in\left(0,\frac{1}{2}\right]$, the operator $\mathcal{M}:H^{-s}(\Gamma)\to H^{-s-\frac{1}{2}}(\Omega)$ defined as

[TABLE]

is bounded, i.e., there exists a constant $M$ such that $\forall f\in H^{-s}(\Gamma)$, $\forall v\in H^{s+\frac{1}{2}}_{0}(\Omega)$:

[TABLE]

Proof 3.1**.**

For all $\varphi\in\mathcal{D}(\Omega)$ and $f\in H^{-s}(\Gamma)$, $s\in\left(0,\frac{1}{2}\right]$, we can write

[TABLE]

We apply Theorem 1 to the second argument on the right hand side

[TABLE]

and the first part of the thesis follows with a density argument.

We observe that for any $f\in H^{-s}(\Gamma)$, with $s\leq\frac{1}{2}$, the operator $\mathcal{M}f$ belongs to $H^{-1}(\Omega)$, and can be used as a standard source term in the variational formulation of the Poisson problem. In particular, if $f\in H^{m}(\Gamma)$, with $m\geq 0$, it is also in $H^{-\epsilon}(\Gamma)$ for any $\epsilon>0$.

The final variational problem can be formulated as follows.

Problem 3.2** (Variational).**

Given $b\in L^{2}(\Omega)$ and $f\in H^{-s}(\Gamma)$, $s\in\left(0,\frac{1}{2}\right]$, find $u\in H^{1}_{0}(\Omega)$ such that:

[TABLE]

where

[TABLE]

We indicate with $(\cdot,\cdot)$ the $L^{2}(\Omega)$ inner product, and with $\langle\cdot,\cdot\rangle$ the duality product between $H^{1}_{0}(\Omega)$ and $H^{-1}(\Omega)$.

In particular, we have that the exact solution $u$ satisfies the following regularity result:

Lemma 3.3** (Continuous dependence on data).**

Problem 3.2 is well posed and has a unique solution in $H^{\frac{3}{2}-s}(\Omega)$ that satisfies

[TABLE]

Proof 3.4**.**

For convex domains, with Lipschitz boundary $\partial\Omega$, Problem 3.2 is 2-regular, and admits a unique solution that satisfies the estimate:

[TABLE]

whenever the right hand side of Problem 3.2 is in $H^{k}(\Omega)$, $-1\leq k\leq 0$.

Exploiting Theorem 2 and taking $k=-\frac{1}{2}-s$ in (19), we obtain the thesis.

This formulation does not require any approximation of the Dirac delta, since the regularity of $\mathcal{M}$ is compatible with test functions in $H^{1}_{0}(\Omega)$, allowing a natural approximation by Galerkin methods, using, for example, finite elements BoffiGastaldi-2003-a ; BoffiGastaldi-2007-Numerical-stability-0 ; Heltai-2008-a .

4 Finite element approximation

We consider a decompositions of $\Omega$ into the triangulation $\Omega_{h}$, consisting of cells $K$ (quadrilaterals in 2D, and hexahedra in 3D) such that

$\overline{\Omega}=\cup\{\overline{K}\in\Omega_{h}\}$; 2. 2.

Any two cells $K,K^{\prime}$ only intersect in common faces, edges, or vertices. On $\Omega_{h}$ we define the finite dimensional subspace ${W}^{\ell}_{h}\subset H^{1}_{0}(\Omega)$, such that

[TABLE]

where $\mathcal{Q}^{\ell}(K)$ is a tensor product polynomial space of degree $\ell$ on the cells $K$, and $N_{W}$ is the dimension of the finite dimensional space.

The index $h$ stands for the maximum radius of $K$, and we assume that $\Omega_{h}$ is shape regular, i.e., $\rho_{K}\leq h\leq C\rho_{K}$ where $\rho_{K}$ is the radius of the largest ball contained in $K$, and the inequality is valid for a generic constant $C>0$ independent on $h$. We assume, moreover, that

[TABLE]

where $C_{0}$ is a generic constant independent on $h$, and with $|K\cap\Gamma|$ we indicate the $n-1$ dimensional Hausdorff measure of the portion of $\Gamma$ that lies in $K$. These assumptions are generally satisfied whenever $h$ is sufficiently small to capture all the geometrical features of both $\Gamma$ and $\Omega$.

Using the space $W^{\ell}_{h}$, the finite dimensional version of Problem 3.2 can be written as

Problem 4.1** (Discrete).**

Given $b\in L^{2}(\Omega)$, $f\in H^{-s}(\Gamma)$, $s\in\left(0,\frac{1}{2}\right]$, find $u_{h}\in W^{\ell}_{h}\subset H^{1}_{0}(\Omega)$ such that

[TABLE]

where

[TABLE]

Theorem 3** **(A-priori error estimates)

The finite element solution $u_{h}$ of Problem 4.1 satisfies the following a-priori error estimates for $f\in H^{-s}(\Gamma)$, $s\in\left(0,\frac{1}{2}\right]$

[TABLE]

and

[TABLE]

for $f\in L^{2}(\Gamma)$, where $u$ is the solution to Problem 3.2.

Proof 4.2**.**

For the finite element approximation defined in Problem 4.1, we expect error estimates of the type

[TABLE]

where $u$ is the solution to Problem 3.2 and $u_{h}$ is the solution to Problem 4.1  ciarlet78 . The thesis follows applying Lemma 3.3 with $k=\frac{3}{2}-s$.

The low regularity of $\mathcal{M}f$ affects the numerical approximation of the problem, and produces sub-optimal (w.r.t. to the approximation degree $\ell$) error estimates when standard finite elements are used. This phenomena is known and was observed in the literature of the variational formulation of the Immersed Boundary Method BoffiGastaldiHeltaiPeskin-2008-a ; Heltai2012b , motivating this work.

The solution is expected to be at least $H^{2}(\Omega\setminus\Gamma)$, and it seems reasonable to assume that the numerical solution to such problems is sub-optimal only in a small neighbourhood of $\Gamma$.

In this work we show that the numerical solution obtained by solving Problem 4.1 with the classical finite element method applied directly to Problem 3.2, results in approximate solutions which converge optimally when we measure the error with properly chosen weighted Sobolev norms. Similar ideas are presented in Agnelli2014 , and in Dangelo2012 ; DAngelo2008 . Rigorous proofs of the well posedeness of the approximation of Poisson type problems where the source is given by a singular measure and the domain is a convex polygonal or polyhedral domain are given in Drelichman2018 .

A generalization to non-convex domains and to Stokes problem is available in Otarola2018 , provided that singularities are away from the boundary of the domain $\Omega$.

5 Weighted Sobolev spaces

The set $\Gamma$ has zero measure, so it is reasonable to introduce weighted Sobolev norms, where the weight is chosen to be an appropriate power of the distance from the co-dimension one surface $\Gamma$. Such weight mitigates the jump of the gradient of $u$ across $\Gamma$, allowing a more regular variational formulation of the problem in the entire $\Omega$ in a rigorous and numerically convenient way.

We define the Hilbert space of measurable functions (see, e.g., Kufner1985 ; Turesson2000 )

[TABLE]

equipped with the scalar product

[TABLE]

In (27) and (28), $d$ is the distance between the point $x$ and the surface $\Gamma$, that is

[TABLE]

For any $\alpha$ in $\left(-\tfrac{1}{2},\tfrac{1}{2}\right)$, the weighting function $w:\mathbb{R}^{n}\to\mathbb{R}_{+}$ defined by $w(x):=\text{dist}(x,\Gamma)^{2\alpha}$, is a Muckenhoupt class $A_{2}$-weight, that is

[TABLE]

where $B_{r}(x)$ is the ball centered at $x$ with radius $r$, and $|B|$ is its measure Muckenhoupt1972 ; Kufner1985 ; Turesson2000 .

The identity $\langle u,v\rangle=\langle ud^{\alpha},d^{-\alpha}v\rangle$ implies that $L^{2}_{-\alpha}(\Omega)$ is contained in the dual space of $L^{2}_{\alpha}(\Omega)$. We denote by $\|{u}\|_{0,\alpha,\Omega}$ the norm of a function $u$ in $L^{2}_{\alpha}(\Omega)$, i.e., $\|d^{\alpha}u\|_{0,0,\Omega}$, where $d$ is defined in (29). From now on we will use the notation $\|u\|_{0,0,\Omega}\equiv\|u\|_{0,\Omega}$ for all $u\in L^{2}(\Omega)$. The duality pairing between $L^{2}_{\alpha}(\Omega)$ and $L^{2}_{-\alpha}(\Omega)$ is indicated with $\langle u,v\rangle$. The usual inequality $\langle u,v\rangle=\langle d^{\alpha}u,d^{-\alpha}v\rangle\leq\|u\|_{0,\alpha,\Omega}~{}\|v\|_{0,-\alpha,\Omega}$ follows from Schwartz inequality in $L^{2}(\Omega)$.

Similarly, we define the weighted Sobolev spaces

[TABLE]

where $s\in\mathbb{N}$, $\gamma$ is a multi-index and $D^{\gamma}$ its corresponding distributional derivative. These weighted Sobolev spaces can be equipped with the following seminorms and norms

[TABLE]

and we define the Kondratiev type weighted spaces $V^{s}_{\alpha}(\Omega)$ using the same seminorms, but weighting them differently according to the index of derivation, i.e.,

[TABLE]

equipped with the following norm:

[TABLE]

For example, $|||u|||_{1,\alpha,\Omega}:=|u|_{1,\alpha,\Omega}+||u||_{0,\alpha-1,\Omega}$. The norms in $V^{1}_{\alpha}$ and $H^{1}_{\alpha}$ are equivalent, but not uniformly with respect to $\alpha$ Kufner1985 .

We define the space

[TABLE]

with norm $\|\cdot\|_{1,\alpha,\Omega}$ and we denote with $W^{\prime}_{\alpha}$ its dual space with norm

[TABLE]

Notice that an equivalent definition is obtained using the $|||\cdot|||$ norms.

Lemma 5.1**.**

Given $\alpha\in(-\frac{1}{2},\frac{1}{2})$ and $\epsilon\geq 0$ such that $\epsilon+\alpha\in(-\frac{1}{2},\frac{1}{2})$, the embeddings $H^{m}_{\alpha}\hookrightarrow H^{m}_{\alpha+\epsilon}$ and $W^{\prime}_{\alpha+\epsilon}\hookrightarrow W^{\prime}_{\alpha}$ are continuous.

Proof 5.2**.**

For $\epsilon\geq 0$ the function $d^{2\epsilon}$ is bounded and continuous on $\overline{\Omega}$, by Hölder inequality, we have that

[TABLE]

The thesis follows by applying the inequality in Equation (33) to $D^{\gamma}u$ for all $|\gamma|\leq m$, and we get

[TABLE]

where $C_{\epsilon}=\max_{x\in\Omega}d^{\epsilon}(x)$.

The dual case follows applying the definition of the dual norm.

Theorem 4** **(Isomorphism of $-\Delta$)

The Laplace operator $-\Delta$ is an isomorphism from $W_{\alpha}$ to $W^{\prime}_{-\alpha}$, for any $\alpha$ in $(-\frac{1}{2},\frac{1}{2})$.

Proof 5.3**.**

The proof is an immediate consequence of (Drelichman2018, , Theroem 2.7): taking $k=n-1$ and $p=2$ in that theorem, one obtains that the Laplace operator $-\Delta$ is an isomorphism from $W_{\alpha}$ to $(W_{-\alpha})^{\prime}$ (and we will write this space as $W^{\prime}_{-\alpha}$).

Lemma 5.4** (Weighted space of Dirac terms).**

The operator $\mathcal{M}f$ defined in Problem (12) is in $W^{\prime}_{-\alpha}$ for any $\alpha$ in $[0,\frac{1}{2})$.

Proof 5.5**.**

Since both $z$ and $p$ belong to $H^{1}(\Omega)$, by Lemma 5.1, they also belongs to $H^{1}_{\alpha}(\Omega)$, for $\alpha\in[0,1/2)$. Theorem 4 applied to $z+p$ implies that the function $-\Delta(z+p):=b+\mathcal{M}f$ belongs to $W^{\prime}_{-\alpha}$, independently on the choice of $b\in L^{2}(\Omega)$ and $f\in H^{-s}(\Gamma)$, with $s\in\left(0,\frac{1}{2}\right]$.

The weighted variational problem can be written as follows

Problem 5.6** (Weighted variational).**

Given $b\in L^{2}(\Omega)$ and $f\in H^{-s}(\Gamma)$, for any $\alpha$ in $[0,\frac{1}{2})$, find $u\in W_{\alpha}$ such that

[TABLE]

where

[TABLE]

Lemma 5.7** (Continuous dependence on data – weighted case).**

Problem 5.6 is well posed and has a unique solution that satisfies

[TABLE]

Proof 5.8**.**

The chain of inequalities come directly from Theorem 4, Lemma 5.1, and from Theorem 2. For the second inequality we first consider the term

[TABLE]

In particular,

[TABLE]

since

[TABLE]

For the second term, we exploit the definition of the operator $\mathcal{M}$:

[TABLE]

and use a trace inequality and Lemma 5.1:

[TABLE]

The thesis follows with $C=\max\left\{\sqrt{\|d^{2\alpha}\|_{L^{\infty}(\Omega)}},C_{T}C_{\alpha}\right\}$.

Remark 5.9** (Non-weighted formulation).**

We remark here that the standard variational formulation, presented in Problem 3.2, is a special case of Problem 5.6 when the power of the distance is taken to be zero. In this case $W_{\alpha}=W_{-\alpha}=H^{1}_{0}(\Omega)$ and we recover Problem 3.2. Moreover, the embedding $W_{-\alpha}\hookrightarrow H^{1}_{0}(\Omega)$ implies that the (unique) solution to Problem 3.2 is the same as the solution to Problem 5.6, and it belongs to $W_{\alpha}\cap H^{\frac{3}{2}-s}(\Omega)$ for any $\alpha\in\left[0,\frac{1}{2}\right)$.

Lemma 5.10** (Finite dimensional subspace of $W_{\alpha}$).**

The finite dimensional space ${W}^{\ell}_{h}$ is a subspace of $W_{\alpha}$ for all $\alpha$ in $(-\frac{1}{2},\frac{1}{2})$.

Proof 5.11**.**

The space ${W}^{\ell}_{h}$ is a subspace of $H^{1}_{0}(\Omega)\cap C^{0}(\overline{\Omega})\cap W^{1,\infty}(\Omega)$, i.e., finite element functions are Lipschitz continuous. Following the same lines of Lemma 5.1, we have:

[TABLE]

and similarly for the gradients:

[TABLE]

This allows us to use the same finite dimensional space for both test and trial functions in the numerical implementation, which is identical to the classical Galerkin approximation.

Theorem 5** **(Stability of Problem 5.6)

Let $|\alpha|<\frac{1}{2}$, then there exist $C_{1}$, $C_{2}$, and $C_{3}$ such that

[TABLE]

The Galerkin approximation (Problem 4.1) of Problem 5.6 in $W^{\ell}_{h}$ is stable, and the following error estimate holds

[TABLE]

Proof 5.12**.**

The two inf-sup conditions follow from (Drelichman2018, , Theorem 3.4) setting $k=n-1$, and $p=2$, and observing that $2\alpha=\lambda$. Estimate 40 follows using standard results for Petrov-Galerkin approximations.

From the point of view of the implementation, Problem 4.1 is a standard finite element problem. The only difficulty is given by the integration of the test functions on the surface $\Gamma$, which is not aligned with the grid where $v_{h}$ are defined. This is usually done by some quadrature formulas, where the integration is performed approximately using a fixed number of points on $\Gamma$ (see, e.g., Heltai2012b ). It is possible to choose a quadrature formula such that the error induced by the numerical approximation of the integral over $\Gamma$ is of higher order with respect to the overall order of accuracy, by making sure that the integration on $\Gamma$ is performed by splitting the curve or surface at the boundaries of the elements of the triangulation used for $\Omega$, and by using enough quadrature points.

5.1 Interpolation estimates in weighted Sobolev spaces

The estimate provided in Theorem 5 does not exploit the particular structure of our singular forcing term. We adapt the construction of the interpolation operator from $W_{\alpha}$ to $W_{h}^{\ell}$ presented in Dangelo2012 to the co-dimensional one case, and show that optimal error convergence rates can be achieved when the singular forcing term has the form introduced in this paper.

For each element $K\in\Omega_{h}$ we define the quantities

[TABLE]

We assume that the mesh is quasi-uniform, i.e., there exist two positive constants $c$ and $C$ such that

[TABLE]

We split the mesh in two parts, one containing all elements close to $\Gamma$, that is

[TABLE]

where $\sigma$ is a fixed safety coefficient, and the other, $\Omega_{h}^{\text{out}}:=\Omega_{h}\setminus\Omega_{h}^{\text{in}}$.

It can be shown that

[TABLE]

where the notation $a\lesssim b$ is used to indicate that there exists a constant $C>0$ such that $a\leq Cb$.

We define the folloing discrete norm

[TABLE]

Lemma 5.13**.**

Let $|\alpha|<t$, $t\in[0,\frac{1}{2})$, then the following norms are equivalent

[TABLE]

where the constants of the inequalities depend only on $t$.

Proof 5.14**.**

The proof follows very closely (Dangelo2012, , Lemma 3.2). We consider $\alpha\geq 0$ (the other case follows similarly). Let $K\in\Omega_{h}$ and $x\in K$; we have $d(x)^{2\alpha}\leq(\overline{d}_{K})^{2\alpha}$, so that the first part of the inequality $\|u_{h}\|_{0,\alpha,\Omega}\leq\|u_{h}\|_{h,\alpha}$ follows trivially.

The second part of the inequality follows by distinguishing two cases: If $K\in\Omega^{\text{out}}_{h}$, we use Equation (42), and we have $\overline{d}_{K}^{2\alpha}\|u_{h}\|^{2}_{0,0,K}\lesssim d^{2\alpha}_{K}\|u_{h}\|^{2}_{0,0,K}\leq\|d^{\alpha}u_{h}\|^{2}_{0,0,K}$.

Following (Dangelo2012, , Lemma 3.2), we show that a similar estimate holds true if $K\in\Omega^{\text{in}}_{h}$. Let $\hat{K}$ be the reference element and let $F_{K}:\hat{K}\to K$ be the affine transformation, mapping $\hat{K}$ onto the actual element $K$. Let $\hat{u}_{h}=u_{h}\circ F_{K}$, and let $\hat{\Gamma}_{K}=F_{K}^{-1}(\Gamma)$, such that $\Gamma$ is the $K$ image of $\hat{\Gamma}_{K}$ under $F_{K}$, and let $\hat{d}(\hat{x})=\text{dist}(\hat{x},\hat{\Gamma}_{K})$.

Thanks to shape regularity, the eigenvalues of the Jacobian matrix of $F_{K}$ are uniformly upper and lower bounded by $h_{K}$. Hence, distances are transformed according to $d(F_{K}(\hat{x}))\gtrsim h_{K}\hat{d}(\hat{x})$. As a result,

[TABLE]

Let us introduce the subset $\hat{K}_{\Delta}=\{\hat{x}\in\hat{K}:\text{dist}(\hat{x},\hat{\Gamma}_{K})>\Delta\}$, where $\Delta>0$ is a parameter (see Figure 4); we have

[TABLE]

Note that for $\Delta$ small, $\inf_{K}|\hat{K}_{\Delta}|$ cannot degenerate with respect to $|\hat{K}|\gtrsim 1$.

We can estimate $|\hat{K}|-|\hat{K}_{\Delta}|\lesssim\Delta$, irrespective of the position of $\hat{\Gamma}_{K}$; hence, choosing $\Delta$ sufficiently small, we have $|\hat{K}_{\Delta}|\geq c^{\prime}\Delta|\hat{K}|$, where the constant $c^{\prime}=1-O(\Delta)$ depends on $\Delta$ but not on $K\in\Omega^{\text{in}}_{h}$.

Notice that, when the surface is of co-dimension one, as in our case, the constant $c^{\prime}$ is of order $O(\Delta)$ instead of order $)(\Delta^{2})$ as in the case presented by Dangelo2012 .

Similarly, it can be seen that $\|\hat{u}_{h}\|^{2}_{0,0,\hat{K}_{\Delta}}\geq c_{\Delta}\|\hat{u}_{h}\|^{2}_{0,0,\hat{K}}\quad\forall\hat{u}_{h}\in\mathcal{Q}^{\ell}(\hat{K})$ (it suffices to see that the estimates hold for the local basis functions on $\hat{K}$), where again the constant $c_{\Delta}$ depends only on $\Delta$ and not on the shape of $\hat{K}_{\Delta}$. Hence, using $\alpha\leq t$ and Equation (42), we conclude

[TABLE]

Lemma 5.15**.**

Given $u\in H^{s}(\Omega)$, with $0\leq m\leq s\leq\ell$, $\alpha\in\left[0,\frac{1}{2}\right)$, we have that

[TABLE]

Proof 5.16**.**

By the definition of the weighted norm, we have that

[TABLE]

and similarly for the higher order semi-norms.

A standard scaling argument, on the other hand, implies

[TABLE]

and the thesis follows by a combination of the two inequalities.

5.2 Convergence rate of the finite element approximation

We consider the nodal points $x_{j}$ of the basis functions $\phi_{i}$ that span the space $W^{\ell}_{h}$, i.e., $W^{\ell}_{h}:=\text{span}\{\phi_{i}\}_{i=1}^{N}$, and $\phi_{i}(x_{j})=\delta_{ij}$, where $\delta$ is the Kronecher delta. We define the interpolation operator

[TABLE]

as the operator that coincides with the standard finite element interpolation operator in $\Omega_{h}^{\text{out}}$ (see, for example, ciarlet78 ), and that sets the degrees of freedom whose support belongs to $\Omega^{\text{in}}_{h}$ to zero, i.e.:

[TABLE]

Different interpolation operators could be defined for more general weighted Sobolev spaces, as in Nochetto2016 . In this work we provide a generalisation of a result in (Dangelo2012, , section 3.3), for co-dimension one surfaces.

Theorem 6** **(Properties of $\Pi_{h}$)

Let $\ell\leq k+1$ be such that the embedding $H^{\ell}(\Omega)\hookrightarrow L^{\infty}(\Omega)$ is continuous. For $\alpha\geq\epsilon$, $m\leq\ell$, the operator $\Pi_{h}$ satisfies the following inequalities:

[TABLE]

where $\Delta_{K}$ is the set of all elements $K^{\prime}$ in $\Omega^{\text{out}}_{h}$ that share at least a node with $K$, i.e.,

[TABLE]

Proof 5.17**.**

If $K$ belongs to $\Omega^{\text{out}}_{h}$, the first inequality follows from standard finite element theory (see, for example, ciarlet78 ). Let’s consider then $K$ in $\Omega^{\text{in}}_{h}$. If $K$ does not share at least one node with $\Omega^{\text{out}}_{h}$, $\Pi_{h}u$ is identically zero on $K$, and the second inequality follows trivially. Let us consider then the case in which $K$ shares the node $x_{i}$ with $K^{\prime}\subset\Omega^{\text{out}}_{h}$, and assume that for each element $K^{\prime}$ there exists an affine transformation such that $K^{\prime}=F_{K^{\prime}}(\hat{K})$, and $\hat{u}:=u\circ F_{K^{\prime}}^{-1}$.

In this case we can write

[TABLE]

We start by estimating $|u(x_{i})|$,

[TABLE]

where we used i) standard scaling arguments for Sobolev norms, ii) the fact that since $K^{\prime}$ is in $\Delta_{K}$, and therefore it is sufficiently close to $\Gamma$, we can write $\bar{d}_{K^{\prime}}\lesssim h_{K^{\prime}}$ and $h_{K^{\prime}}\lesssim d_{K^{\prime}}$, iii) and the fact that for a negative power $q^{-}$, we have $|u|_{m,0,K^{\prime}}\leq\bar{d}^{\ \ -q^{-}}_{K^{\prime}}|u|_{m,q^{-},K^{\prime}}$ while for a positive power $q^{+}$ we have $|u|_{m,0,K^{\prime}}\leq d^{\ \ -q^{+}}_{K^{\prime}}|u|_{m,q^{+},K^{\prime}}$.

To estimate the second term in the right-hand side of equation (49), ${|\phi_{i}|_{m,\alpha,K^{\prime}}}$, we use again a scaling argument for Sobolev norms:

[TABLE]

where to obtain the last inequality we have used that ${|\hat{\phi}_{i}|_{m,0,\hat{K}}}{\lesssim 1}$.

Combining (49) and (51), and summing over all $K^{\prime}\in\Delta_{K}$ we get the second inequality of the thesis.

We are now in the position to prove our main result.

Theorem 7

Let $u$ be the exact solution to Problem 5.6, and let $u_{h}$ be the solution to Problem 4.1 in $W_{h}^{\ell}$ with $\ell\geq 1$. In particular $u\in H^{2}(\Omega\setminus\Gamma)\cap H^{\frac{3}{2}-s}(\Omega)\cap W_{\alpha}$, and

[TABLE]

where $s\in\left(0,\frac{1}{2}\right]$, and $\alpha\in\left[0,\frac{1}{2}\right)$.

Proof 5.18**.**

Consider $K\in\Omega^{\text{out}}_{h}$. Using Lemma 5.15 and property (48), we easily obtain

[TABLE]

where we used that $\overline{d}_{K}^{\alpha}\leq|\Omega|^{\alpha}=C$ for $\alpha\geq 0$.

Similarly, for $K$ in $\Omega^{\text{in}}_{h}$, we have

[TABLE]

where we can estimate the first term in the right hand side of equation (54) using lemma 5.15, and the second term using the second property in equation (48). Theorem 5 implies the thesis.

6 Numerical validation

All numerical examples provided in this section were obtained using an open source code based on the deal.II library Bangerth2007 ; Bangerth2016b ; MaierBardelloniHeltai-2016-a ; AlzettaArndtBangerth-2018-a and on the deal2lkit toolkit Sartori2018 .

We construct an artificial problem with a known exact solution, and check the error estimates presented in the previous section. We begin with a simple two-dimensional problem, where we impose the Dirichlet data in order to produce a harmonic solution in the domain $\Omega\setminus\Gamma$, with a jump in the normal gradient along a circular interface, and extend the same test case to the three-dimensional setting.

Recalling our main result (Theorem 7), we have that the exact solution of this problem is in $H^{2}(\Omega\setminus\Gamma)$ (it is in fact analytic everywhere except across $\Gamma$) and globally $H^{\frac{3}{2}}(\Omega)$, therefore we expect a convergence rate in weighted Sobolev spaces of the type:

[TABLE]

Notice that when $\alpha\to 1/2$, then the estimate tends to the optimal case, while when $\alpha=0$, the estimate is classical (suboptimal due to the lack of global regularity in the solution). These situations occur very often in numerical simulations of boundary value problems with interfaces, and the results we present here show that a proper variational formulation of the interface terms results in suboptimal convergence property of the finite element scheme. However such sub-optimality is a local property due to the non-matching nature of the discretisation, and it only influences the solution close to the surface $\Gamma$. If we take this into account when measuring the error, for example using a weighted Sobolev norm as we do in this work, we recover the optimal estimate.

The results we presented in Theorem 7 can be applied also in the case of higher order approximations. However, the theory shows that there would be no improvement in the global convergence rate beyond $\frac{3}{2}-s-m+\alpha$, making the use of polynomial order greater than one redundant. This is not to say that the method would not benefit from higher order polynomiasl away from the co-dimension one surface $\Gamma$, even though the global approximation error would not converge with the expected higher order rate.

6.1 Two-dimensional case

In the two-dimensional case, we define the exact solution to be

[TABLE]

The curve $\Gamma$ is a circle of radius $0.2$ with center in $c$. This function is the solution to the following problem:

[TABLE]

We use a bi-linear finite dimensional space $W^{1}_{h}$, and show a plot of the numerical solution for $h=1/1024$ in Figure 5. We compute the error in the weighted Sobolev norms $\|\cdot\|_{0,\alpha,\Omega}$ and $\|\cdot\|_{1,\alpha,\Omega}$ for values of $h$ varying from $1/4$ to $1/1024$, and values of $\alpha$ varying from zero (standard Sobolev norms in $L^{2}(\Omega)$ and $H^{1}(\Omega)$) to $0.49$.

We report the errors in the weighted $L^{2}_{\alpha}(\Omega)$ norm in Table 1 and in Figure 6, and for the weighted $H^{1}_{\alpha}(\Omega)$ norm in Table 2 and in Figure 7.

From the tables we verify the results of Theorem 7, and we observe rates of convergence in the standard $L^{2}(\Omega)$ and $H^{1}(\Omega)$ Sobolev norms which are coherent with the $H^{3/2}(\Omega)$ global regularity of the solution. In particular we expect a convergence rate of order $3/2$ for the $L^{2}(\Omega)$ norm and $1/2$ for the $H^{1}(\Omega)$ norm. When increasing $\alpha$ to a value close to $1/2$, we observe that the errors in the weighed norms converge to the optimal rates (in this case two and one).

6.2 Three-dimensional case

In the three-dimensional case, we define the exact solution to be

[TABLE]

The surface $\Gamma$ is a sphere of radius $0.2$ with center in $c$. This function is the solution to the following problem:

[TABLE]

We use a bi-linear finite dimensional space $W^{1}_{h}$, and show a plot of the numerical solution for $h=1/128$ in Figure 8. We compute the error in the weighted Sobolev norms $\|\cdot\|_{0,\alpha,\Omega}$ and $\|\cdot\|_{1,\alpha,\Omega}$ for values of $h$ varying from $1/4$ to $1/128$, and values of $\alpha$ varying from zero (standard Sobolev norms in $L^{2}(\Omega)$ and $H^{1}(\Omega)$) to $0.49$.

We report the errors in the weighted $L^{2}_{\alpha}(\Omega)$ norm in Table 3 and in Figure 9, and for the weighted $H^{1}_{\alpha}(\Omega)$ norm in Table 4 and in Figure 10.

From the tables we verify again the results of Theorem 7 in the three dimensional case. We observe rates of convergence in the standard $L^{2}(\Omega)$ and $H^{1}(\Omega)$ Sobolev norms which are coherent with the $H^{3/2}(\Omega)$ global regularity of the solution. In particular we expect a convergence rate of order $3/2$ for the $L^{2}(\Omega)$ norm and $1/2$ for the $H^{1}(\Omega)$ norm. When increasing $\alpha$ to a value close to $1/2$, we observe that the errors in the weighed norms converge to the optimal rates (in this case two and one).

7 Conclusions

One of the major point against the use of immersed boundary methods and their variational counterparts, comes from the unfavourable comparison in convergence rates that can be achieved using matching grid methods (ALE Hirt1974 ; DoneaGiulianiHalleux-1982-a ), or enriching techniques (IIM Leveque1994 , X-FEM Mittal2005b ).

In this work we have shown that this detrimental effect on the convergence properties is only a local phenomena, restricted to a small neighbourhood of the interface. In particular we have proved that optimal approximations can be constructed in a natural and inexpensive way, simply by reformulating the problem in a distributionally consistent way, and by resorting to weighted norms when computing the global error of the approximation, where the weight is an appropriate power of the distance from the interface. Weighted Sobolev spaces Kufner1985 ; Turesson2000 , provide a natural framework for the study of the convergence properties of problems with singular sources Agnelli2014 or problems with singularities in the domain Belhachmi2006 ; Duran2009a .

The method we have presented has the great advantage of not requiring any change in the numerical approximation scheme, which is maintained the same as if no interface were present, requiring only the construction of a special right hand side, incorporating the jump conditions in a simple singularity Peskin2002 ; BoffiGastaldiHeltaiPeskin-2008-a . The analysis is based on results from Drelichman2018 , while the numerical techniques borrow heavily from the formalism and the results presented in Dangelo2012 ; DAngelo2008 .

Applications of this discretisation technique to fluid structure interaction problems is well known and dates back to the early seventies Peskin1972 . Recent developments in finite element variants are available, for example in Heltai2012b ; Roy2015 , and similar constructions are used to impose boundary conditions on ocean circulation simulations Rotundo2016 . Applications of this technique could be used, for example, to allow more general classes of doping profiles in doping optimization problems for semi-conductor devices PeschkaRotundoThomas2016 ; Peschka2018 .

It is still unclear if the same techniques may be used for the treatment of jumps in the solution itself, as in these cases the regularity gain that could be achieved by weighted Sobolev spaces alone may not be sufficient, and we are currently exploring alternative approximation frameworks, following the lines of Nochetto2016 .

Acknowledgments

N.R. acknowledges support by DFG via SFB 787 Semiconductor Nanophotonics: Materials, Models, Devices, project B4 “Multi-dimensional Modeling and Simulation of electrically pumped semiconductor-based Emitters”.