A priori error estimates of regularized elliptic problems

Approximations of the Dirac delta distribution are commonly used to create sequences of smooth functions approximating nonsmooth (generalized) functions, via convolution. In this work, we show a priori rates of convergence of this approximation process in standard Sobolev norms, with minimal regularity assumptions on the approximation of the Dirac delta distribution. The application of these estimates to the numerical solution of elliptic problems with singularly supported forcing terms allows us to provide sharp $H^1$ and $L^2$ error estimates for the corresponding regularized problem. As an application, we show how finite element approximations of a regularized immersed interface method result in the same rates of convergence of its non-regularized counterpart, provided that the support of the Dirac delta approximation is set to a multiple of the mesh size, at a fraction of the implementation complexity. Numerical experiments are provided to support our theories.


Introduction
Singular source terms are often used in partial differential equations (PDE) to model interface problems, phase transitions, or fluid-structure interaction problems. The immersed boundary method (IBM, [29]) is a good example of a model problem where a Dirac delta distribution supported on an immersed fiber or surface is used to capture complex dynamics that are happening on possibly moving interfaces. Similar forcing terms can be used, for example, to model fictitious boundaries in the domain [13], or to couple solids and fluids across non-matching discretizations [17].
When the discretization scheme is based on variational principles, singular sources may be incorporated exactly in the numerical scheme by the definition of their action on test functions [3,4,17]. However, if the discretization scheme is based on finite differences [29] or finite volumes [27], this procedure is cumbersome, and a typical alternative is to regularize the source term by taking its convolution with an approximation of the Dirac delta distribution, or by modifying the differential operators themselves to incorporate the knowledge about the interface [22].
Several works are dedicated to the design of good Dirac approximations [19,34] to use in this regularization process, and their convergence properties are well known in the literature of immersed boundary methods [21,14].
When the source term is a single Dirac distribution, [19] proved convergence in both the weak- * topology of distributions, as well as in a weighted Sobolev norm, provided that certain momentum conditions are satisfied.
In many applications, however, the Dirac distribution is not used directly as a source term, but only in formal convolutions to represent the coupling between overlapping domains. In this context, the resulting source term may or may not be singular, according to the co-dimension of the immersed domain itself. To fix the ideas, consider a source term of the kind The above formalism is used in immersed methods to represent a (possibly singular) coupling between B and Ω. Here δ is a d-dimensional Dirac distribution. If the codimension of the immersed domain B is zero, then the above forcing term reduces to where χ B is the indicator function of B, owing to the distributional definition of δ. However, if the co-dimension of B is greater than zero, the integration over B does not exhaust the singularity of the Dirac distribution: the resulting F is still singular, and it should be interpreted as the distribution whose effect on smooth functions ϕ is given by: f (y)ϕ(y) dy, for all ϕ ∈ C ∞ c (Ω).
In the three dimensional case, assuming that f ∈ L 2 (B), the regularity of F goes from being L 2 (Ω) when the co-dimension of B is zero, to a negative Sobolev space which cannot be smoother than H −1/2 (Ω), H −1 (Ω), and H −3/2 (Ω) when B is a Lipschitz surface, Lipschitz curve, or point, respectively. In all cases, a regularization of F is possible by convolution with an approximate (possibly smooth) Dirac function. In most of the literature that exploits this technique from the numerical point of view, pointwise convergence, truncation, and Taylor expansions are used to argue that high order convergence can be achieved in L p (Ω) norms for the numerical approximation of the regularized problem, provided that some specific conditions are met in the construction of the approximate Dirac function [28,23,24]. Convergence of the regularized solution to the exact solution of the original problem is usually left aside, with the notable exceptions of the work by [19], where only the degenerate case of B being a single point is considered, and the work by [30], where suboptimal estimates in W −1,p (Ω) norm for some p ∈ [1, d d−1 ) are provided for finite element discretizations of the immersed boundary method applied to the Stokes problem.
Even though the mollification (also known as regularization) technique is widely used in the mathematical analysis community, very little is known about the the speed at which regularized functions converge to their non-regularized counterparts in standard Sobolev norms, when non-smooth approximations of the Dirac delta distribution are considered.
In this work we provide convergence results in standard Sobolev norms that mimic closely the a-priori estimates available in finite element analysis, and we investigate the performance of the finite element method using regularized forcing data, in the energy and L 2 (Ω) norms, for a class of singular source terms commonly used in interface problems, fluid structure interaction problems, and fictitious boundary methods.
We start by providing an a priori framework in standard Sobolev spaces for unbounded domains, to study the convergence of regularized functions to their non-regularized counterparts, with minimum requirements on the regularization kernel.
We extend this framework to bounded domains, assuming that the support of the forcing data is away from the physical domain (see the definition in (14)), and we study the convergence speed of elliptic problems with regularized forcing terms to their non-regularized counterparts, with minimum requirements on the support of the forcing term. For compactly supported kernels, in Theorem 14, we provide sharp convergence estimates in the energy norm in terms of powers of the radius of the support. We also provide an L 2 (Ω) error estimate by following a duality argument (or Aubin-Nitsche trick) and using the H 2 interior regularity of a dual problem; we refer to Theorem 15 for more details. We note that although we only consider Dirichlet boundary conditions in this paper, the convergence results that we derive for the regularized problem can be also applied to elliptic problems with other boundary conditions; see Remark 16. As an application, we investigate how the regularization affects the total error of the finite element approximation of an interface problem via immersed methods, and we show that all the estimates we obtain are sharp. Even though a regularization in this case is not necessary when using the finite element method [3], the biggest advantage of using regularization comes from the fact that its numerical implementation is trivial; see Remark 22. Theorem 17 and 18 show that the regularization does not affect the overall performance of the finite element approximation in both energy and L 2 (Ω) norms when choosing the regularization parameter to be a multiple of the maximal size of the quasi-uniform subdivision of Ω.
The rest of the paper is organized as follows. We first introduce some notations in Section 2, and provide general results in both bounded and unbouned domains in Section 3. These results are applied to a general model elliptic problem in Section 4, where we use the results of Section 3 to derive sharp error estimates in energy and L 2 (Ω) norms. An application to immersed methods is provided in Section 5, where we apply the theory to a finite element approximation of an interface problem. Finally we discuss the numerical implementation of the singular forcing data (with and without regularization) and validate our findings via numerical simulations in Section 6.

Notations and preliminaries
In what follows, we set Ω to be a bounded Lipschitz domain in R d with d = 2 or 3. We write A B when A ≤ cB for some constant c independent of regularization and discretization parameters (mentioned in Section 1) in A and B. We indicate A ∼ B when both A B, and B A. For a ∈ R, we denote with a − (or a + ) any real number strictly smaller (or greater) than a.
For x ∈ R d , we use |x| to indicate the euclidean norm of x, and given a normed space X, we denote by X and ·, · X ,X its dual space and the duality pairing, respectively. We also denote by · X and · X the norm of X and the operator norm of X , i.e., 2.1. Sobolev spaces. For s ∈ N and p > 1, we denote by W s,p (Ω), H s (Ω) and L 2 (Ω) the usual Sobolev spaces. For convention we set H 0 (Ω) = L 2 (Ω) and denote with (·, ·) Ω the L 2 (Ω) inner product. We denote with | · | H s (Ω) the usual semi-norm of H s (Ω), and we set | · | H 0 (Ω) := · L 2 (Ω) . Let H 1 0 (Ω) be the set of functions in H 1 (Ω) vanishing on ∂Ω. It is well known that H 1 0 (Ω) is the closure of the space of infinitely differentiable functions with compact support in Ω (denoted by C ∞ c (Ω) or D(Ω)) with respect to the norm of H 1 (Ω). For s ∈ (0, 1), H s (Ω) denotes the space of functions whose Sobolev-Slobodeckij norm where A is either R d or Ω.

Regularization
In the mathematical literature, mollifiers were introduced by Sobolev [33], and formalized by Friedrichs [12], to provide an effective way of regularizing non smooth or singular functions by performing the convolution with a smooth and compactly supported function that integrates to one (the kernel, or mollifier).
The outcome of the procedure is a function that has the same regularity property of the kernel, normally taken to be C ∞ . In this section we relax the requirements of the kernel, and we provide some basic results that can be obtained by regularizing with approximated Dirac distributions that are not necessarily C ∞ . In particular, we consider functions ψ that satisfy the following assumptions: ψ(x) is compactly supported, with support supp(ψ) contained in B r0 (the ball centered in zero with radius r 0 ) for some r 0 > 0; 2. k-th order moments condition: Lemma 2 (Convergence to δ). A function ψ that satisfies Assumption 1 for some k ≥ 0, defines a one-parameter family of Dirac delta approximations, i.e., define where δ(x) is the Dirac delta distribution and the limit must be understood in the space of Schwartz distributions.
Notice that Assumption 1.2 with k = 0 implies that both ψ and δ ε have unit integral. Moreover, we make no additional assumptions on the global regularity of ψ, except from requiring it to be in L ∞ (R d ).
Any function ψ that satisfies Assuption 1 defines a one-parameter family of Dirac delta approximations δ ε , through which it is possible to define a regularization in both R d and Ω: where δ ε is defined as in Eq.
Lemma 4 (L 1 growth control). A Dirac approximation δ ε constructed from a function ψ that satisfies Assumptions 1 (irrespective of k ≥ 0), also satisfies the following polynomial growth condition: where the hidden constant depends on m, d and the choice of ψ.
Proof. By considering the change of variable x = ξε, we observe that We shall consider two common methods for choosing ψ(x). In the first case we choose a one dimensional function ψ ρ ∈ L ∞ (R) so that ψ ρ is supported in [0, 1], and we define the function ψ as where I d is a scaling factor, chosen so that ψ(x) integrates to one. With this choice, ψ satisfies Assumption 1 with k at least equal to one (i.e., it satisfies the first moments conditions). The second construction method is usually referred to as tensor product construction. We start from a L ∞ (R) function ψ 1d that satisfies Assumption 1 for some k in dimension one. Then we define ψ(x) in dimension d by The tensor product approximation ψ(x) satisfies Assumption 1 with r 0 = √ d and the same k of ψ 1d . In particular k = 0 if ψ 1d is not symmetric, and k is equal to at least one if ψ 1d is symmetric. We refer to [19,Section 3] for an in depth discussion on other possible choices of Dirac approximation classes ψ(x) with possibly higher order moment conditions. 3.1. Unbounded domains. We begin by providing some results that follow from an application of Young's inequality for convolutions: Lemma 5 (Young's inequality for convolutions [35]).
, and let v ε be defined by Definition 3. Then there holds Proof. 1 We start by considering v ∈ C ∞ c (R d ), and the case where s is integer, and 1 ≤ s ≤ k + 1. By Taylor expansion, it is possible to expand v(y) around an arbitrary point x in a polynomial part p x and a residual part r x , i.e.: Here α and β are multi-indices. If we regularize v, since the regularization is a linear operator, we obtain: where in the last equality follows from Assumption 1.2 so that p ε x = p x for any polynomial of order up to k.
Applying the definition of r x (y), Fubini's theorem, and by the change of variable ξ = x + t(y − x) for a fixed x ∈ R d , and t ∈ (0, 1), we have: Applying Lemma 5 and 4 we get 2 To show (8) for s = 0, we note that the regularization operator is L 2 -stable, i.e., for any θ ∈ L 2 (R d ), Here we applied again Young's inequality (Lemma 5) together with Lemma 4 for 3 Taking the sup over θ with unit L 2 (R d ) norm, and applying interpolation estimates between s = 0 and s = k +1, the proof is complete by a density argument.
The above lemma immediately implies the following convergence result in R d : Here we identify v in the negative Sobolev space by the duality pairing: Proof. Let us show the desired estimates in three steps.
3.2. Bounded domains. The generalization of the previous results in bounded domains is non-trivial, due to the presence of boundaries. We start by providing some results that work well when restricting v to a region which is strictly contained in Ω. Let this region be defined by a Lipschitz domain ω ⊂ Ω, and assume that there exists a positive constant c 0 such that (14) dist(∂ω, ∂Ω) > c 0 .
The following result relies on the interior regularity of v. More precisely, according to (14), we assume that the regularization parameter satisfies ε ≤ ε 0 ≤ c 0 for some fixed ε 0 , and we set (15) ω Assuming that v ∈ H s (ω ε0 ) ∩ L 1 (Ω) with s ∈ [0, k + 1], we next provide similar error estimates compared with the unbounded case in the previous section. We note that results below are instrumental to our error estimates for the numerical approximation of the model problems in the next section.
, and let v ε be defined as in Definition 3. Then there holds Proof. We denote by · the zero extension from ω, ω ε0 or Ω to R d . Since ε < dist(∂ω, ∂Ω), for any x ∈ ω, B ε (x) ⊂ Ω and For any θ ∈ L 2 (ω) and v ∈ C ∞ 0 (ω ε0 ), we follow the proof of Lemma 6 to get (17) ( Here we apply again Assumption 1.2 for the last equality above. When it comes to the last inequality in (17), we note that for a fixed x ∈ ω and for any y ∈ ω ε0 , the change of variable ξ = t(y − x) + x belongs to ω ε0 for any t ∈ (0, 1). Hence we proceed following the proof of Lemma 6, Step 1 and apply Lemma 5 again to conclude the proof for a positive integer s. Replacing v and θ with v and θ in the proof of Lemma 6, Step 2, we obtain (16) with s = 0. The assertion for any s ∈ [0, k + 1] follows from the interpolation between s = 0 and s = k + 1.
, and let v ε be defined by Definition 3. Then there exists ε 1 > 0 small enough so that for ε < ε 1 , Proof. Let m be a positive integer. Integration by parts yields that for v ∈ C ∞ c (ω ε0 ) with x ∈ ω ε0/2 and for ε < ε 0 /2, Repeating the above argument shows that given a multi-index β so that |β| = m, D β v ε (x) = (D β v) ε (x) for x ∈ ω ε0/2 m when ε < ε 0 /2 m . We then apply (16) in Lemma 8 with s = 0 to ( |β|≤m D β v) ε and a density argument to obtain that Interpolating the above estimate with (16) yields the desired estimate.
Proof. The proof is identical to the proof of Theorem 7, replacing the application of Lemma 6 with that of Lemma 8 and Corollary 9.
The above results are summarized in the following theorem: Theorem 11 (Regularization estimates in Ω).
Define the extension of ω as in (15), i.e., where k is the order of the moments conditions satisfied by δ ε as in Assumption 1.2, we have:

Model problem
We are now in a position to apply the results of Theorem 11 to a model elliptic problem. Let A(x) be a symmetric d × d matrix. We assume that all entries of A(x) are in C 1 (Ω), uniformly bounded, and that A(x) is positive definite, i.e., there exist positive constants a 0 , a 1 satisfying a 0 |ξ| 2 ≤ ξ A(x)ξ ≤ a 1 |ξ| 2 , for all ξ ∈ R d and x ∈ Ω.
We also set c(x) to be a nonnegative function in C 0,1 (Ω). We now define the forcing data for our model problem. We set F ∈ H −1 (Ω). We assume that for some s ∈ (0, 1 2 ], supp s−1 (F ) ⊂ ω. Based on the above definitions, our model problem reads: find the distribution u satisfying To approximate Problem (20) using the finite element method, we shall consider its variational formulation: find u ∈ V := H 1 0 (Ω) satisfying (21) A Regularity. Our error estimates rely on standard regularity results for elliptic problems: given g ∈ V , let T : V → V be the solution operator satisfying We first note that if g ∈ L 2 (Ω), we identify g, · V ,V with (g, ·) Ω and hence T g has the H 2 interior regularity, i.e., given a subset K such that K ⊂ Ω, T g ∈ H 2 (K) and where the hiding constant depends on K and Ω; we refer to [11, Theorem 1 of Section 6.1] and [7, Theorem 5.33] for a standard proof. The following assumption provides the regularity of T g up to the boundary Assumption 12 (elliptic regularity). There exists r ∈ (0, 1] and a positive constant C r satisfying T g H 1+r (Ω) ≤ C r g H −1+r (Ω) .
As an example, consider the case where Ω is a polytope, A(x) is the identity matrix, and c(x) = 0, i.e., A becomes the Dirichlet form Based the regularity results provided by [6], r in Assumption 12 is between 1 2 and 1 and can be decided by the shape of Ω. Assumption 12 also implies that the solution u in (21) belongs to H 1+min{s,r} (Ω) ∩ H 1 0 (Ω). 4.1. Analysis of a regularized problem. Now we are ready to define a regularized problem of (21): find u ε ∈ V satisfying (25) A Remark 13. Notice that we denote with u ε the solution to Problem 25, and in this case the superscript ε does not denote a regularization (hence the different font used for u ε ).
In order to bound the error between u and u ε , we use the boundedness and coercivity of A(·, ·), and obtain that

L 2
(Ω) error estimate. We next show the convergence rate for u ε in L 2 (Ω) norm. To this end, we additionally assume that δ ε satisfies the moments conditions in Assumption 1.2 with k ≥ 1.
Theorem 15 (L 2 (Ω) error estimate). Under the assumptions of Lemma 10, and Assumption 1.2 with k ≥ 1, let u and u ε be the solutions of problems (21) and (25), respectively. For a fixed ε 0 < c 0 and ε < ε 0 /2, there holds Proof. We consider the following dual problem: find z ∈ V such that Hence, we choose v = u − u ε and obtain that Due to the interior regularity of z, u − u ε ∈ H 1 0 (Ω) ⊂ L 2 (Ω) implies that We continue to estimate the right hand side of (26) by where in the second inequality, we invoke Lemma 8 for z. Combing (26) and (27) concludes the proof of the theorem.

Remark 16.
We have to point out that the estimates in Theorem 14 and 15 also hold for Problem (20) with other boundary conditions. These error estimates depend on the smoothness of the test function and of the solution for the dual problem in the neighborhood of ω while ω is away from the boundary. The analysis of the case where ω is attached to the boundary ∂Ω is left aside for future investigation.

Application to immersed methods
The general results presented in the previous section can be applied immediately to immersed interface and immersed boundary methods [15,18,16,2]. In this section, we consider an interface problem whose variational formulation can be written as in the model problem (21), and we shall consider its finite element approximation using the regularization of the forcing data given by the application of Definition 3 to functions in negative Sobolev spaces.

5.
1. An interface problem via immersed methods. Let Γ ⊂ Ω be a closed Lipschitz interface with co-dimension one. We set ω be the domain inside Γ, and we assume that ω satisfies the assumptions introduced in Section 2 (see Figure 1). For x ∈ Γ, ν(x) denotes the normal vector and ∂u ∂ν denotes the corresponding normal derivative. We also use the notation · for the jump across Γ. More precisely speaking, letting u − = u ω and u + = u Ω\ω , the jump across Γ is defined by The above problem can be reformulated in the entire domain Ω using a singular forcing term that induces the correct jump on the gradient of the solution. (cf. [20]). The variational fomulation of the new problem is provided by (21) where A is the Dirichlet form (24) and Here δ denotes the d-dimensional Dirac delta distribution, and Eq. (29) should be interpreted variationally, i.e., We note that for any s in [0, 1 2 ), M is a bounded map from H s−1/2 (Γ) to H s−1 (ω ε1 )∩ H −1 (Ω) and hence the variational formulation makes sense. In fact, for v ∈ H 1 (Ω), Here we apply the trace Theorem (e.g. [26, Theorem 1.5.1.2]) for the second inequality. The above estimate together with (30) shows that F ∈ V . It also implies that F ∈ H s−1 (ω ε1 ) and hence supp(F ) s−1 ⊆ ω. So F satisfies the setting in Section 4 and we can apply the results from Theorem 14 and 15 to the interface problem to get (31) u We note that Fubini's Theorem yields that which is the classical formulation of F ε that can be found in the literature of the Immersed Boundary Method [29], where δ ε is always taken to be even (i.e., k ≥ 1), and F ε is introduced as the regularization of f on Γ, via the d-dimensional Dirac approximation δ ε .

Finite Element Approximation of the regularized problem.
In this section, we consider a finite element approximation of the regularized problem (25). Assume that Ω is polytope and let {T h (Ω)} h>0 be a family of conforming subdivisions of Ω made of simplices with h denoting their maximal size. We assume that T h (Ω) are shape-regular and quasi-uniform in the sense of [10,5]. Let V h be the space of continuous piecewise linear functions subordinate to T h (Ω) that vanish on ∂Ω. Let I h : H 1 0 (Ω) → V h be the Scott-Zhang interpolation [32] which has the following approximation property . The discrete counterpart of the regularized problem (25) reads: find u ε h ∈ V h satisfying For v h ∈ V h , F ε , v h V ,V can be computed by using a quadrature formula on Γ and a quadrature formula on τ ∈ T h (Ω). We refer to the next Section for the details of the implementation. The following theorem shows the error between u and its final approximation u ε h . Theorem 17 (H 1 (Ω) error estimate). Let u and u ε h be the solutions to (21) and (34), respectively. Under Assumption 12 and 1, we have Proof. The coercivity of A(·, ·) implies that A(·, ·) is also V h elliptic. The first Strang's Lemma (see, e.g. [5, Theorem 4.1.1]) yields Setting v h = I h u and invoking (33) together with Assumption 12 and Lemma 10, we conclude that We next show a L 2 (Ω) error estimate between u and u ε h . Theorem 18 (L 2 (Ω) error estimate). Following the settings from Theorem 17, we additionally assume that δ ε satisfies Assumption 1.2 with k ≥ 1. Then we have Proof. We first provide a bound on the error between u ε and u ε h under the regularity assumption of f . In fact, the triangle inequality together with Theorem 15 and 17 implies that Next we bound u ε − u ε h L 2 (Ω) using a duality argument. Let Z ∈ H 1 0 (Ω) satisfy that for all v ∈ H 1 0 (Ω). Hence Z ∈ H 1+r (Ω). Applying Galerkin orthogonality we obtain that This, together with the regularity estimate Z H 1+r (Ω) (Ω) together with the above estimate and the L 2 (Ω) estimate in Theorem 15 concludes the proof of the theorem.
We end this section with some remarks to further explain the error estimates and also for the numerical experiments in the next section.
Remark 19 (knowing the regularity of the solution). If the regularity of the solution is known, e.g., u ∈ H 1+β (Ω) for some β ∈ (0, 1], we can apply the interpolation estimate (33) with s = β in Theorem 17 and 18 to get Remark 20 (choices of ε). Usually we can choose ε = ch q for some q ∈ (0, 1] and for a fixed factor c. Hence, error estimates in the above remark become h min{β+r,r+sq,(1+s)q} .

Numerical Illustrations
We implement the linear system of discrete problem (34) using the deal.II finite element Library [1,25,31]. Before validating our error estimates given by Theorem 17 and 18 via a series of numerical experiments, we want to make some remarks on the computation of the right hand side vector in discrete linear system.
Remark 21 (Approximation of the surface Γ). We shall compute the right hand side data on an approximation of a C 2 interface Γ using the technology of the surface finite element method for the Laplace-Beltrami problem [9,8]. Let Γ h0 be a polytope which consists of simplices with co-dimension one, where h 0 denotes the maximal size of the subdivision. All vertices of these simplices lie on Γ and similar to T h (Ω), we assume that this subdivision of Γ h0 , denoted by T h0 (Γ), is conforming, shape-regular, and quasi-uniform. We compute the forcing data at is the signed distance function for Γ. The error analysis for how the finite element approximation of u in (28) is affected by Γ h0 is out of the scope of this paper. In what follows, we assume that h 0 is small enough compared with h so that the error from the interface approximation does not affect the total error of our test problems.
Remark 22 (Comparing with the usual approach). Based on the previous remark and given a shape function φ h ∈ V h , we can approximate the right hand side data F ε , φ h V ,V using quadrature rules on both τ 1 ∈ T h (Ω) and τ 2 ∈ T h0 (Γ), i.e., Here τ ε 1 follows from the definition (15) and supp(φ h ) denotes the support of φ h . {w j1 , q j1 } Jτ 1 j1=1 and {w j2 , q j2 } Jτ 2 j2=1 are pairs of quadrature weights and quadrature points defined on τ 1 ∈ T h (Ω) and τ 2 ∈ T h0 (Γ), respectively. On the other hand, letting Q be the collection of quadrature points for all τ 2 ∈ T h0 (Γ), the finite element method allows one to approximate directly F, v h V ,V by We note that we usually compute φ h (q j1 ) by using a transformation B τ1 mapping from the reference simplexτ to τ 1 ∈ T h (Ω) and set q j1 = B τ1qj1 and w j1 = w j1 | det B τ1 |, where {ŵ j1 ,q j1 } is the pair of quadrature weights and quadrature points forτ . In the implementation, we store the reference shape functionφ h , the shape valuesφ h (q j1 ) = φ h (q j1 ) and w j1 offline (and similar implementations for τ 2 ) but we have to compute φ h (q j2 ) online byφ h (B −1 τ1 q j2 ). We finally note that searching for q j2 ∈ Q and τ 2 ∈ T h0 (Γ) intersecting with τ 1 can be accelerated by creating R-trees for q j2 and bounding boxes of τ 2 , respectively. Then we use the searching algorithms for intersection with the bounding box of τ 1 ∈ T h (Ω). In our implementation, we use the searching algorithms from the Boost.Geometry.Index library.
6.1. Tests on the unit square domain with different δ ε (x). We test the interface problem (28) on the unit square domain but with a non-homogeneous Dirichlet boundary condition, following the examples provided in [18]. We set Ω = (0, 1) 2 and solve (35) −∆u = 0, in Ω\Γ, 3) , f = 1 0.2 and g = ln(|x − c|). The analytic solution is In view of Assumption 12 and Remark 20, r = 1, s = 1 2 and u ∈ H 3/2 − (Ω). Hence setting ε = h yields where #DoFs stands for the number of degree of freedoms and we used the fact that h ∼ #DoFs −1/d for quasi-uniform meshes. We shall compute u ε h on a sequence of unstructured, quasi-uniform, quadrilateral meshes: we start to compute u ε h1 on a coarse mesh of Ω in Figure 2. Then we compute the next approximated solution in a higher resolution based on the global refinement of the previous mesh. In the meantime, we also take the global refinement of approximated interface according to Remark 21. The right plot in Figure 2 shows the approximated solution on the mesh produced from the coarse mesh with six-time global refinement (744705 degree of freedoms).
In Figure 3 we report L 2 (Ω) and H 1 (Ω) errors against #DoFs using the following types of δ ε (x): • Radially symmetric C 1 : use (5) with ψ ρ (x) = (1 + cos(|πx|))χ B1(0) (x)/2, where χ B1(0) (x) is the characteristic function on B 1 (0); • Tensor product C 1 : use (6) with ψ 1d (x) = (1 + cos(|πx|))χ (−1,1) (x)/2; • Tensor product C ∞ : use (6) with ψ 1d (x) = e 1−1/(1−|x| 2 ) χ (−1,1) (x); • Tensor product L ∞ : use (6) with ψ 1d (x) = 1 2 χ (−1,1) (x). The predicted rates are observed in all four types. Errors in the first three types behave similar while errors for the last type are larger than those in the other cases.   #DoFs − min{1/2,1/3+q/4,3q/4} . Figure 4 provides the coarse mesh of the numerical test and also the solution on the mesh with six-time global refinement (1574913 degrees of freedoms). In Figure 5 we also report L 2 (Ω) and H 1 (Ω) errors against #DoFs with q = 0.2, 0.4, 0.6, 0.8, 1 using the δ ε (x) with the type Tensor product C 1 . We again observed the predicted convergence rates from Figure 5.  Figure 6 shows the unstructured coarse mesh of the unit cube as well as the approximated solution on the mesh after the forth-time global refinement (2324113 degrees of freedoms). In Figure 7 we plot the L 2 (Ω) error decay by setting ε = h and the error decay without using Dirac delta approximations as mentioned in Remark 22. We also plot the CPU time for the computation of the right hand side vector against #DoFs with or without regularization in Figure 7. We note that the the computer we use has 2.2 GHz Intel Core i7 with 16GB memory. Figure 7 shows that both the computation time for the right hand side assembling and error using the regularization approach are comparable to those using the usual method. We have to point out that according to Remark 22, we cannot evaluate B −1 τ1 explicitly when the mesh is unstructured and here we use the Newton iteration instead. So if the geometry of each element is simple such as cube, the computation time without regularization can be reduced significantly.   Figure 7. Plots of L 2 (Ω) error decay (left) and CPU time for right hand side assembling against the number of degree of freedoms using with or without regularization according to Remark 22. We regularize the forcing data using the type Tensor product C1 with ε = h.