A modified $P_1$ - immersed finite element method

In recent years, the immersed finite element methods (IFEM) introduced in \cite{Li2003}, \cite{Li2004} to solve elliptic problems having an interface in the domain due to the discontinuity of coefficients are getting more attentions of researchers because of their simplicity and efficiency. Unlike the conventional finite element methods, the IFEM allows the interface cut through the interior of the element, yet after the basis functions are altered so that they satisfy the flux jump conditions, it seems to show a reasonable order of convergence. In this paper, we propose an improved version of the $P_1$ based IFEM by adding the line integral of flux terms on each element. This technique resembles the discontinuous Galerkin (DG) method, however, our method has much less degrees of freedom than the DG methods since we use the same number of unknowns as the conventional $P_1$ finite element method. We prove $H^1$ and $L^2$ error estimates which are optimal both in order and regularity. Numerical experiments were carried out for several examples, which show the robustness of our scheme.


Introduction
In recent years, there have been some developments of immersed finite element methods for elliptic problems having an interface. These methods use meshes which do not necessarily align with the discontinuities of the coefficients [19], [20], thus violate a basic principle of triangulations in the conventional finite element methods [4], [10]. However, when the basis functions are modified so that they satisfy the interface conditions, they seem to work well [9], [19], [20]. These methods were extended to the case of Crouzeix-Raviart P 1 nonconforming finite element method [11] by Kwak et al. [17], and to the problems with nonzero jumps in [7]. Some related works on interface problems can be found in [5], [15], [16], [18], [21], [22], [25].
On the other hand, the discontinuous Galerkin methods (DG) where one uses completely discontinuous basis functions were developed and have been studied extensively, see [1], [2], [12], [23] and references therein. The DG methods work quite well for problems with discontinuous coefficient in the sense that they capture the sharp changes of the solutions well, yet they require large number of unknowns and the meshes have to be aligned with the discontinuity.
The purpose of this paper is to combine the advantages of the two methods. We use a DG type idea of adding the consistency terms to the IFEM, thus proposing a modified version of IFEM based on piecewise linear conforming basis functions on triangular grids. In spirit, it resembles [14] in the sense that the standard linear basis functions are used for noninterface elements and line integrals are added, but in our method the line integrals along the edges, not along the interface, are added. Furthermore, our method incorporate the flux jump conditions to the basis functions hence requires no extra unknowns along the interface as in [14].
We prove error estimates in the mesh dependent H 1 -norm and L 2 -norm which are optimal both in the order and the regularity. We carry out various numerical Figure 1. A domain Ω with interface tests to confirm our theory and compare the performance with the unmodified scheme.

Preliminaries
Let Ω be a connected, convex polygonal domain in R 2 which is divided into two subdomains Ω + and Ω − by a C 2 interface Γ = ∂Ω + ∩ ∂Ω − , see Figure 1. We assume that β(x) is a positive function bounded below and above by two positive constants. Although our theory applies to the case of nonconstant β(x), we assume β(x) is piecewise constant for the simplicity of presentation: there are two positive constants β + , β − such that β(x) = β + on Ω + and β(x) = β − on Ω − . Consider the following elliptic interface problem with the jump conditions along the interface where f ∈ L 2 (Ω) and u ∈ H 1 0 (Ω) and the bracket [·] Γ means the jump across the interface: [u] Γ := u| Ω + − u| Ω − .
For any domain D, we let W m p (D) be the usual Sobolev space with (semi)-norms and denoted by | · | m,p,D and · m,p,D . When p = 2, we write H m (D) := W m 2 (D) with the (semi)-norms | · | m,D and · m,D . Let H 1 0 (Ω) be the subspace of H 1 (Ω) with zero trace on the boundary.
Let {T h } be the usual shape regular triangulations of the domain Ω by the triangles of maximum diameter h which may not be aligned with the interface Γ. . We also need some subspaces of H 2 (T ) and H 2 (Ω) satisfying the jump conditions: Throughout the paper, the constants C, C 0 , C 1 , etc., are generic constants independent of the mesh size h and functions u, v but may depend on the problem data β, f and Ω, and are not necessarily the same on each occurrence.
Theorem 2.1. Assume that f ∈ L 2 (Ω). Then the variational problem (2.4) has a unique solution u ∈ H 2 (Ω) which satisfies We briefly review the immersed finite element space based on the P 1 -Lagrange basis functions ( [19], [20]). We call an element T ∈ T h an interface element if the interface Γ passes through the interior of T , otherwise we call it a noninterface element. Let T I h be the collection of all interface elements. We assume that the interface meets the edges of an interface element at no more than two points.
We construct the local basis functions on each element T of the partition T h . For a noninterface element T ∈ T h , we simply use the standard linear shape functions on T whose degrees of freedom are functional values on the vertices of T , and use S h (T ) to denote the linear spaces spanned by the three nodal basis functions on T : S h (T ) = span{ φ i : φ i is the standard linear shape function } We let S h (Ω) denote the space of usual continuous, piecewise linear polynomials with vanishing boundary values. Now we consider a typical interface element T ∈ T I h whose geometric configuration is given as in Fig. 2. Here the curve between the two points D and E is a part of the interface and DE is the line segment connecting the intersections of the interface and the edges.
These are continuous, piecewise linear functions on T satisfying the flux jump condition along DE, whose uniqueness and existence are known [9], [19]. We denote by S h (T ) the space of functions generated byφ i , i = 1, 2, 3 constructed above. Next we define the global immersed finite element space S h (Ω) to be the set of all where A i , i = 1, 2, 3 are the vertices of T and we callÎ h u the local interpolant of u in S h (T ). We naturally extend it to H 2 Γ (Ω) by (Î h u)| T =Î h (u| T ) for each T . Then we have the following approximation property [17], [20].
for all u ∈ H 2 Γ (Ω). Now the IFEM based on the conforming P 1 -finite element method introduced in [19], [20] reads: The error estimate for this scheme is shown in [9].

Modified P 1 -IFEM
In this section, we modify the P 1 -IFEM above by adding the line integrals for jumps of fluxes and functional values. The method resembles the discontinuous Galerkin methods (see [2], [13], [23] and references therein) which use completely discontinuous basis functions, but the degrees of freedom in our method are much smaller than the DG methods since our method has the same number of basis functions as the conventional P 1 -FEM.
In order to describe the new method, we need some additional notations. Let the collection of all the edges of T ∈ T h be denoted by E h and we split E h into two disjoint sets; is the set of edges lying in the interior of Ω, and E b h is the set of edges on the boundary of Ω. In particular, we denote the set of edges cut by the interface Γ by E I h . For every e ∈ E o h , there are two element T 1 and T 2 sharing e as a common edge. Let n Ti , i = 1, 2 be the unit outward normal vector to the boundary of T i , but for the edge e, we choose a direction of the normal vector, say n e = n T1 and fix it once and for all. For functions v defined on T 1 ∪ T 2 , we let [·] e and {·} e denote the jump and average across e respectively, i.e. [ We also need the mesh dependent norm ||| · ||| on the space H h (Ω), Multiplying both sides of the equation (2.1) by v ∈ H 1 (T ), applying Green's formula and adding, we get By using the preassigned normal vectors n e and adding the unharmful term ǫ e {β∇v· n e } e [u] e for any ǫ, we see the above equation becomes Now, for each ǫ = 0, ǫ = −1 and ǫ = 1, we define the modified P 1 -IFEM for the problem (2.1)-(2.3): . This is similar to a class of DG methods, corresponding to IP, SIPG, NIPG and OBB ( [1], [13], [12], [3]), if ǫ = 0, ǫ = −1, ǫ = 1, and ǫ = 0, σ = 0, respectively. Remark 4.2. We take any positive σ if ǫ = 1, and σ > 0 large enough if ǫ = 0 or −1. In our scheme, a small positive σ or even σ = 0 works for all the cases we have tested. This is in contrast to the usual DG schemes, where sufficiently large σ is necessary. Thus our scheme is robust in choosing σ. The reason seems to be that, unlike the usual DG, the a h (u, v) form is coercive and the stability term j σ (u, v) is necessary only to control the line integral term b ǫ (u, v).

Error analysis
In this section, we prove an optimal order of error estimates in H 1 -and L 2 -norm of our schemes. For simplicity, we present the case with ǫ = −1 only. All other cases are similar.
Lemma 5.1. There exist positive constants C 0 , C 1 , C 2 independent of the function v such that for all v ∈ P k (T ), Now we show the following interpolation error estimate for the mesh dependent norm ||| · |||.
Proof. By (5.2) and Proposition 3.1, we see the first term satisfies For the next term, we fix an edge e and let T 1 and T 2 be two elements which share e as the common edge. Apply (5.2) to the restriction u −Î h u| e∩∂Ti , i = 1, 2, and use Proposition 3.1 to get Multiplying h 1/2 and adding with respect to all T ∈ T h we get Thus the proof is complete.

H 1 -error analysis.
First we check that the modified P 1 -IFEM is consistent.
Lemma 5.5. Let u be the solution of (2.1)-(2.3) and let u m h be the solution of (4.3). For any v h ∈ S h (Ω), we have In other words, a ǫ (u − u m h , v h ) = 0. Proof. By (4.2), the definition of the a ǫ form and the homogeneous jump condition of u, we have Now we can prove the H 1 -error estimate which is optimal both in order and the regularity.

Numerical Experiments
For numerical tests, we solve the problem (2.1)-(2.3) on the rectangular domain Ω = [−1, 1] × [−1, 1] partitioned into unform right triangles having mesh size h. Three types of interface problems are considered with various values of parameter β. We measure u − u h 0 , u − u h 1,h and e u − u h 0,e which are very close to the theoretical orders of convergence, 2, 1 and 1.5 respectively. Although not reported, we also measured e ∂(u − u h )/∂n 0,e , the orders of which agree with the theoretical value 0.5.  Tables 1 and 2, we see the modified method performs better than the original P 1 -IFEM when β + /β − = 10. Example 6.2 (Sharp corner). In this example, we consider an interface with a sharp corner having interior angle 2θ. Let Γ be the zero set of L(x, y) = −y 2 + ((x − 0.6) tan θ) 2 (x + 0.4) for x ≤ 0.6. We test the cases with θ = 10 • and θ = 45 • when β + /β − = 10 and 1000. The exact solution is u = L(x,y) β . We see that both methods work well when θ = 10 • , but modified method works better when θ = 45 • and β + /β − = 10; See Tables 3, 4, 5 and 6.   Table 3. Example 5.2 (Sharp corner) : θ = 10 • , β − = 1, β + = 10 Example 6.3 (Variable coefficient). Finally, we consider the case with variable coefficient. We take the level set of L(x, y) = x 2 /(0.9) 2 + y 2 /(0.5) 2 − 1.0 as an interface. The exact solution is u = L(x, y)/β(x, y) where In this case, both methods show an optimal order of convergence. The modified method performs slightly better in energy norm and edge norm; See Table 7.