Residual-based a posteriori error estimation for immersed finite element methods

In this paper we introduce and analyze the residual-based a posteriori error estimation of the partially penalized immersed finite element method for solving elliptic interface problems. The immersed finite element method can be naturally utilized on interface-unfitted meshes. Our a posteriori error estimate is proved to be both reliable and efficient with reliability constant independent of the location of the interface. Numerical results indicate that the efficiency constant is independent of the interface location and that the error estimation is robust with respect to the coefficient contrast.

dimensional elliptic interface problems [16,18,21,26,32,34]. Recently, the partially penalized immersed finite element method was introduced in [25]. Compared to classical IFEM, the partially penalized IFEM contains normal flux jump terms on interface edges to ensure the consistency of the scheme. In addition, the new IFEM includes a stabilization term on interface edges to guarantee the stability of the scheme. The partially penalized IFEM significantly improves the numerical approximation, especially the accuracy around the interface. The optimal a priori error estimate is theoretically proved for partially penalized IFEM in the energy norm [25].
We note that the partially penalized IFEM can solve elliptic interface problems accurately on uniform Cartesian meshes provided that the exact solution of interface problems is piecewise smooth, and the contrast of the coefficient is "moderate". However, for interface problems that also involve singularity or steep gradient, the partially penalized IFEM alone may not be efficient to obtain an accurate approximation on uniform meshes. In such cases, it is necessary to apply certain adaptive mesh refinement (AMR) strategy to IFEM. The goal for AMR is to obtain an approximate solution within prescribed error tolerance with minimum computational cost which is particularly rewarding for interface problems with non-smooth solutions.
The key success for AMR is the a posteriori error estimation which provides both global and local information on the approximation error. Moreover, in many applications even without the intention of performing adaptive mesh refinement, the a posteriori error estimation is also important in assessing the quality of the simulation by providing an effective error control. We note that the a priori error estimation of IFEM is getting mature in the past decade, but the a posteriori error estimation of IFEM is still in the infancy. In this paper, we will develop and analyze a residualbased a posteriori error estimate for the partially penalized IFEM for second-order elliptic interface problems in the two dimensional space.
Comparing to the residual-based error estimation for the classical finite element solution [6], the newly introduced estimator in this work additionally includes the jump of tangential derivative of numerical solutions on interface edges besides the standard terms of element residual and jump of normal flux on all edges. This is necessary since the numerical solution is in general discontinuous across the interface edges due to the construction of the IFEM basis functions. Moreover, the new error estimator also includes the geometrical fitting error due to the polygonal approximation of the curved interface. Theoretically, we prove that the error estimation is both globally reliable and locally efficient. To prove reliability, we use the Helmholtz decomposition and a L 2 representation technique recently introduced in [8,9]; more-over, we introduce a new type of Clément-type interpolation in the IFEM space that allows us to take advantage of the error equation of IFEM. In the efficiency analysis for IFEM, the technique using the standard bubble functions [33] is invalid because the edge jumps of the normal flux and the tangential derivative become piecewise constant on the interface edges. Instead, we prove the efficiency in two different approaches that aims to provide an optimal efficiency constant for both regular and irregular interface edges.
The rest of the article is organized as follows. In Section 2, we recall the partially penalized IFEM for elliptic interface problems. In Section 3, we introduce our residual-based error estimator specially designed for IFEM. Section 4 and Section 5 are dedicated to the analysis of global reliability and local efficiency, respectively. Finally, in Section 6, we present several numerical experiments to test the performance of our a posteriori error estimators.

Interface Problems and Partially
Assume that meas(Γ D ) > 0. We consider the elliptic interface problem: Here, f ∈ L 2 (Ω), g N ∈ L 2 (Γ N ), and n is the unit vector outward normal to ∂Ω. The notations ∇ and ∇· are the gradient and divergence operators, respectively. Furthermore, assume that Ω is separated by a closed smooth interface curve Γ into Ω + and Ω − such that Ω = Ω + ∪ Γ ∪ Ω − . The diffusion coefficient α is assumed to be a positive piecewise constant function as follows Denote by ρ = α + α − the ratio of the coefficient jump. The solution is assumed to satisfy the following interface jump conditions: where the jump of a function v across the interface Γ is defined by We use the standard notations for the Sobolev spaces. Let Then the variational problem for (2.1) is to find u ∈ H 1 D (Ω) such that where (·, ·) ω is the L 2 inner product on ω. The subscript ω is omitted when ω = Ω.

Triangulation.
In this paper, we only consider the triangular meshes in two dimensions. Let T = {K} be a triangulation of Ω that is regular but not necessarily body-fitted. Denote the set of all vertices of the triangulation T by where N I is the set of all interior vertices, and N D and N N are the sets of vertices on Γ D and Γ N , respectively. Denote the set of all edges of the triangulation T by where E I is the set of all interior edges and E D and E N are the sets of boundary edges on Γ D and Γ N , respectively. For each element K ∈ T , denote by h K the diameter of K, and by N K and E K the sets of all vertices and edges on K, respectively.
For simplicity, we assume that the interface cuts the partition with the following properties: (I) If Γ meets an edge at more than one point, then this edge is part of Γ.
(II) If the case (I) does not occur, then Γ must intersect a triangle at two points, and these two points must be on different edges of this triangle. Based on the above assumption, all triangular elements in the partition can be categorized into two classes: non-interface elements that either has no intersection with Γ or Γ ∩ K ⊂ ∂K, and interface elements whose interior is cut through by Γ.
Denote the set of all interface elements by T int . For each interface triangle K we let Γ K = Γ ∩ K andΓ K be the line segment approximating Γ K by connecting two endpoints of Γ K . Let K + = K ∩ Ω + and K − = K ∩ Ω − . Also we letK + andK − be the two sub-elements of K separated byΓ K . From the setting above, it is easy to see thatK ± is either a triangle or a quadrangle. Also we define which is the region enclosed by Γ K andΓ K . Under the assumption that Γ K is C 2continuous, the area of the S K is of at least O(h 3 K ) [21]. Figure 2.1 provides an illustration of a typical triangular interface element. For an edge F ∈ E, if F is cut through by Γ, i.e., F ∩ Γ = ∅ and F ⊂ Γ, then F is called an interface edge and denote by E int the set of all such interface edges. For each F ∈ E, denote by h F the length of F . Denote by n F = (n 1 , n 2 ) and t F = (−n 2 , n 1 ) the unit vectors normal and tangential to F , respectively. Let K F,1 and K F,2 be the two elements sharing the common edge F ∈ E I such that the unit vector out normal to K F,1 coincides with n F . When F ⊂ ∂Ω, n F is the unit outward vector normal to ∂Ω, and denote by K F,1 the boundary element having the edge F . For a function v that is defined on K F,1 ∪ K F,2 , denote its traces on F by v| 1 F and v| 2 F restricted on K F,1 and K F,2 , respectively. Define the jump of a function v on the edge F by and the average of a function v on the edge F by It is easy to verify that

IFEM approximation.
For simplicity, we assume that the interface does not intersect with the boundary, i.e., E int ⊂ E I . Letα be an approximation of α such thatα For each interface element K ∈ T int , define the local IFE space bỹ where P 1 (w) is the space of all polynomial functions in w of degree no more than 1. The global IFE space S(T ) is then defined to include all functions such that 1. v| K ∈P 1 (K) for all K ∈ T int and v| K ∈ P 1 (K) for all K ∈ T /T int , and 2. v is continuous at every vertex z ∈ N . Note that for each z ∈ N , there exists a unique IFE nodal basis function [21,24], denoted byλ z ∈ S(T ), such that where δ is the Kronecker delta function.
The partially penalized IFEM solution for the interface problem is to find where the bilinear form a h (w, v) is defined by Here may take the values −1, 0, and 1, corresponding to symmetric, incomplete, and non-symmetric IFEM. The constant γ is the stability parameter and needs to be chosen large enough for symmetric and incomplete IFEMs guarantee the coercivity. For non-symmetric IFEM, the constant γ is only required to be positive. For more details on the partially penalized IFEM, we refer readers to [25]. Remark 2.1. By the definition ofα andΓ K it is easy to see thatα = α on all F ∈ E int .
2.3. Inconsistency error. Due to the geometrical approximation of the interface curve Γ by a polygonal interfaceΓ = K∈T intΓ K , the following geometrical inconsistent error exists. By (2.6) and integration by parts we have for any v ∈ S D (T ) Remark 2.2. If the interface Γ is a polygon such thatα = α on each interface element, then the term (2.7) vanishes. In this case, the partially penalized IFEM scheme is consistent. In general, for a curved interface Γ, the global convergence of IFEM will not be affected by such linear approximation of the interface, since the partially penalized IFEM uses piecewise linear approximation [7].
3. Residual-Based A Posteriori Error Estimation. In this section, we introduce the residual-based error estimator for the partially penalized IFEM. We note that the classical residual-based a posteriori error estimation for conforming finite element methods on fitted meshes consists of element residual and the jump of the normal flux on edges. For the IFEM, it is also necessary to include the jump of the tangential derivative on interface edges since the IFEM solution may not be continuous across interface edges.
Define the normal flux jump of u T on each edge by and the tangential derivative jump of u T on each edge by Note that on each interface edge F ∈ E int both j n,F and j t,F are piecewise constant.
For all K ∈ T we define the local error indicator η K by The global error estimator η is then defined by Remark 3.1. If K is a non-interface element, i.e., E K ∩ E int = ∅, the first and second terms in (3.1) vanish. In this case, the local error indicator η K is identical to the residual-based error indicator for the classical body-fitting conforming finite element method [6].

Global Reliability.
In this section, we establish the reliability bound of the global estimator η given in (3.2). For each z ∈ N , let ω z be the union of all elements sharing z as a common vertex. To this end, let N Γ be the set of all vertices z such that meas d−1 (ω z ∩ Γ) > 0, and define where α z = α(x)| ωz and f z is the average value of f on ω z . Remark 4.1. The first term in H f (T ) is a higher order term for f ∈ L 2 (Ω) [12]. It is also well known that for linear finite element methods the edge residuals are dominant. In our adaptive algorithm the element residual is also omitted.

Helmholtz decomposition. Let
: Proof. The proof can be referred to [3,15]. Here we also sketch a proof for the convenience of readers. Letẽ σ := α∇u −α∇ h u T and e T = u − u T . Assume that φ ∈ H 1 (Ω) solves the following equation: Then there exists ψ ∈ H 1 (Ω) such that . By integration by parts and the boundary conditions it is easy to check that This completes the proof of the lemma.
Lemma 4.2. Let φ and ψ be given in (4.1). Then we have the following error representations in the weighted semi-H 1 norm: for any v ∈ S D (T ) and Proof.
Let v ∈ S D (T ) be arbitrary. Applying (4.2), (4.1), and integration by parts gives Applying integration by parts again gives The last equality used (2.5), φ − v = 0 on Γ D , and the facts that By integration by parts and (2.6) we also have which, together with (4.5) and (4.6), gives (4.3).
To prove (4.4), by (4.2), (4.1), integration by parts, and the facts that Hence, we obtain (4.4). This completes the proof of the lemma.

Modified Clément-type interpolation. Define a modified Clément-type interpolation operator
where π z is defined by where λ z andλ z are the classical barycentric hat function and the linear IFEM nodal basis function of S(T ) associated to z, respectively. Note that The following lemma provides the approximation and stability properties of the modified Clément-type interpolation operator.
Then there exists a constant C > 0 that is independent of the mesh size and the location of the interface such that where ω K is the union of all elements sharing at least one vertex with K.
Proof. Without loss of generality, assume that K ∈ T is an interior element. To prove (4.11), by the partition of unity, the triangle inequality, the boundedness of IFE basis functions λ z ∞,K ≤ C(ρ) (Theorem 2.4, [21]), and (4.10), we have To prove (4.12), by the partition of unity, the triangle inequality, the fact that ∇λ z ∞,K ≤ Ch −1 K (Theorem 2.4, [21]), and (4.10), we have Finally, (4.13) follows from the partition of unity, the triangle inequality, the trace inequality, and (4.10): This completes the proof of the lemma. Lemma 4.4. There exists a constant C > 0, independent of the mesh size and the location of the interface, such that and Proof. We first prove the results on the reference element K formed by vertices ((0, 0), (1, 0), (0, 1)) and let F be the edge of K on the x-axis. Without loss of generality, let (a, 0), 0 < a < 1, be the interface point on F and takes the value 0 at both endpoints (0, 0) and (1, 0). By direct calculations, we have Regarding the H 1/2 -norm, we have where K 1 is the subtriangle fromed by ((0, 0), (a, 0), (0, 1)) and It is easy to verify that which, together with the scaling argument, gives (4.14) and (4.15). This completes the proof of the lemma. Lemma 4.5. Let φ be given in (4.1). There exists a constant C independent of the mesh size and the location of the interface such that Proof. By Lemma 4.2, Let v = φ I ∈ S D (T ) be the modified Clément-type interpolation defined in (4.8) of φ. By the partition of unity, the fact thatλ z = λ z for z ∈ N \ N Γ , the Cauchy-Schwarz inequality, (4.9), and (4.10), we have For any F ∈ (E I ∪ E N ) \ E int , it follows from the continuity of φ − φ I on F , Cauchy-Schwarz inequality, and (4.13) that For any F ∈ E int , it then follows from the Cauchy-Schwarz and Young's inequality and (4.13) that where ω F = ω K F,1 ∪ ω K F,2 , which, together with (4.18) and the Cauchy-Schwarz inequality, yields To bound I 3 , we apply the following trace inequality for functions in the IFEM space (see Lemma 3.2 in [25]). Let v ∈ S D (T ) be arbitrary, then It now follows from the Cauchy-Schwarz inequality, (4.20), and (4.12) that Together with (4.14) we obtain To bound I 4 , we use the fact that the Cauchy-Schwarz and triangle inequalities, and (4.13) to obtain which, combining with (4.14), yields Finally, (4.16) is a direct consequence of (4.17), (4.19), (4.21), (4.22) and the Young's inequality. This completes the proof of the lemma.
Lemma 4.6. Let ψ be given in (4.1). There exists a constant C independent of the mesh size and the location of the interface such that Proof. For the first term in (4.4) it follows from the duality inequality, (4.15), the trace and Young's inequalities that For the second term in (4.4) it follows from the Cauchy-Schwarz inequality that (4.26) Proof. First by adding and subtracting proper terms we have (4.27) By (4.1) and (4.2) we have (4.28) Applying the Cauchy-Schwarz and Young's inequalities also gives

Local Efficiency.
In this section, we establish the efficiency bound for the error indicators η K defined in (3.1) for every element K ∈ T . For non-interface elements, the efficiency bound for η K is well known (see [6]) and the key technique is using the local edge and element bubble functions. However, the same technique without modification becomes invalid for interface elements because jumps of the normal flux and the tangential derivative on the interface edges become piecewise constant. To overcome this difficulty it is natural to design more localized bubble functions that allow us to derive the efficiency bounds on F + and F − separately for each F ∈ E int .
For each F ∈ E int we first design auxiliary elements and edge bubble functions associated to F − . ChooseK F,1 ⊂ K F,1 andK F,2 ⊂ K F,2 to be two regular triangular sub-elements such that F − is a common edge of bothK F,1 andK F,2 , to be more precise, Note thatK F,1 andK F,2 are not necessarily in K − F,1 and K − F,2 , respectively. The key requirement here is to make sub-elementsK F,1 ,K F,2 regular while K − F,1 and K − F,2 may not be in general. For each sub-elementK ∈ {K F,1 ,K F,2 }, define the auxiliary element bubble function υK such that (i) υK ∈ P 3 (K), (ii) υK| ∂K ≡ 0 and (iii) υK takes the value 1 at the barycentric center ofK. We also define the auxiliary edge bubble function It is easy to verify that υK forK ∈ {K F,1 ,K F,2 } and υ F − uniquely exist. Let w F − = (j n,F | F − )v F − and wK = vKf K with f K being the average of f on K.
Applying the continuity of α∇u · n F , the divergence theorem, and the Cauchy-Schwarz inequality gives which, combining with the properties of w F − that We now establish the efficiency bound for the element residual f 0,K forK ∈ {K F,1 ,K F,2 }. By the property of wK, the triangle inequality, the divergence theorem, and the Cauchy-Schwarz inequality, we have which, combining with similar properties for wK: By (5.1) and (5.2), we have

Adding proper weights gives
and the constant C is independent of the mesh size and the location of the interface but may depend on ρ. Similarly, one can prove that

and, after adding proper weights, that
Similarly, by defining auxiliary bubble functions on F + , we also have the following local efficiency results on F + : and For the first case we assume that for each F ∈ E int , there exist positive constants λ and Λ such that Lemma 5.1. Let u and u T be the solution in (2.3) and (2.6), respectively. Then there exists a constant C 1 that is independent of the mesh size but may depend on λ, Λ, and ρ such that Proof. In the case that K is not an interface element, (5.8) is a well known result (see, e.g., [8]). In the case that K is an interface element, by (5.3) and (5.5), we have where C 1 is independent of the mesh size but might depend on ρ, λ, and Λ. Similarly, by (5.4) and (5.6), we have where C 1 is independent of the mesh size but might depend on ρ, λ, and Λ. This completes the proof of the lemma. In the above lemma the efficient constant C 1 blows up when the λ and Λ become extreme. However, we note that in the extreme circumstances the related basis functions for IFEM are very "close" to classical FE nodal basis functions thanks to its consistency with standard FE basis functions (see [21] for more detail). The partially penalized IFEM solution then should be also "close" to the classical finite element solution on fitted meshes. In the following lemma we prove the efficiency in a different approach that yields a bounded coefficient constant when the λ and Λ are extreme.
Lemma 5.2. There exists a constant C 2 that is independent of the mesh size and the location of the interface such that for each interface edge F ∈ E int the following efficiency bound holds: Proof. Without loss of generality assume that h + F ≥ h − F . Let K ∈ {K F,1 , K F,2 }. DefineK ⊂ K be the triangle such that (i) F + is a complete edge ofK and (ii) the vertex in K that is opposite to F , denote by z F,K , is also a vertex ofK. It is obvious thatK is regular. Also defineû K be a piecewise linear function such thatû K ≡ u T on K + (orK), ifK ⊂ K + (or K + ⊂K) and thatû K | F is linear. It is obvious that u K uniquely exists. Finally defineû T such thatû T | K =û K , K ∈ {K F,1 , K F,2 }. By the triangle inequality, we have Then by the definition ofû T and a direct computation, we have for each K ∈ {K F,1 , K F,2 }: where ρ K − denotes that radius of the ball inscribed in K − . This, combining with the fact that [21,24] for the consistence of FE and IFE functions when one piece of the interface element is small.) yields Now applying the fact that [[D nûT ]]| F is a constant and the standard efficiency result onK F,1 andK F,2 gives (5.14) Combing (5.11)-(5.14) we have where the constant C 2 is independent of the mesh size and the location of the interface.
Similarly we can prove that where the constant C 2 is also independent of the mesh size and the location of the interface. Finally (5.10) is a direct result of (5.15) and (5.16). Remark 5.1. (5.15) indicates that for the case when λ and Λ are extreme the term h 1/2 F h − F becomes negligible, and, therefore, C 2 can be used as an effective efficiency constant.
Theorem 5.3. Let u and u T be the solution in (2.3) and (2.6), respectively. The following efficiency bound holds for any K ∈ T int : χ is the characteristic function and C 1 and C 2 are the efficiency constants in Lemma 5.1 and Lemma 5.2, respectively.
6. Numerical Results. In this section, we report some numerical results to demonstrate the performance of the residual-based error estimator for partially penalized IFEM.
For the first three examples, we consider a diffusion interface problem with a smooth elliptical interface curve which has been reported in [25,27]. Let Ω = [−1, 1] 2 , and the interface Γ is an ellipse centered at (x 0 , y 0 ) = (0, 0) with horizontal semiaxis a = π 6.28 and the vertical semi-axis b = 3 2 a. The interface separates Ω into two sub-domains, denoted by Ω − and Ω + such that The exact solution to this interface problem is Here β ± > 0 are the diffusion coefficients, and p > 0 is the regularity parameter. In the following, we use ρ = β + β − to denote the ratio of the coefficient jump. Our adaptive mesh refinement follows the standard procedure: We solve the interface problem using IFEM in (2.6), then we compute the residualbased error indicator η K on each element by (3.1). We adopt the equilibration marking strategy, i.e., construct a minimal subsetT of T such that where the threshold θ = 0.5. Finally we refine the marked triangles by newest vertex bisection [28]. The initial mesh is formed by first partitioning the domain into a 4 × 4 congruent rectangles, then cutting each rectangle into two congruent triangles by connecting its diagonal with positive slope.
Example 6.1 (Piecewise smooth solution with moderate jump). In this example, we let ρ = 100 and p = 5 which represents a moderate coefficient contrast and piecewise smooth solution. In Figure 6.1, we list, from left to right, some typical meshes of similar number of elements and degrees of freedom (DOF) generated by the uniform IFEM, the adaptive IFEM, and the adaptive FEM on unfitted mesh [13]. We observe that there is not much local mesh refinement around the interface for the adaptive IFEM in the middle of Figure 6.1, comparing with the uniform mesh in the left plot. This indicates that the errors of the partially penalized IFEM on uniform mesh are almost equally distributed, and IFEM itself can resolve the interface accurately for moderate coefficient jump. Whereas using the finite element method on unfitted meshes requires much more local mesh refinement around the interface (see the plot on the right of Figure 6.1) to resolve the non-smooth feature of interface problems.  In Figure 6.2, we report the convergence of these three methods and the residual-based error estimator for the adaptive IFEM. The slopes of the log(DOF)-log( α 1/2 (∇u− ∇ h u T ) 0,Ω ) and the log(DOF)-log(η) for the adaptive IFEM are both very close to −1/2, which indicates the optimal-order decay of errors with respect to the number of unknowns and, hence, the efficiency of our local error indicators. We use the following efficiency index, to test the efficiency of our residual-based error estimator. The eff-index is very stable at every mesh level and the value is around 3. We note that the errors of uniform IFEM are very close to the errors of adaptive IFEM. This again indicates that the IFEM itself can resolve the interface accurately for moderate coefficient jumps and piecewise smooth solutions. However, using the standard FEM, the magnitudes of errors are much larger than those of IFEM, with similar degrees of freedom. Example 6.2 (Piecewise smooth solution with large jump). In this example, we test the large jump case by choosing ρ = 10 6 . In this case, the true solution possesses a very steep gradient at the interface.
The left plot of Figure 6.3 shows a typical mesh for the adaptive IFEM, which, compared with Figure 6.1, has much denser refinement around the interface. In the right plot of Figure 6.3, we observe the optimal-rate decay of the errors and the estimators. Nevertheless, even if the convergence rate is optimal for the uniform IFEM, the magnitudes of errors are significantly larger than the errors of adaptive IFEM. Hence, applying adaptive mesh refinement is computationally more efficient for interface problems with large coefficient jump even for IFEM. The efficiency indices are between 2.5 -3.5 on all meshes except the first few coarse ones, and they become more stable (close to 3) as the computations reach the asymptotical convergence region. This phenomenon indicates the robustness of the error estimation with respect to the ratio of coefficient jump.
Furthermore, the numerical solutions and the error surfaces on the uniform and adaptive meshes with similar DOFs are depicted in Figure 6.4 and Figure 6.5, respectively. We can observe that the error is significantly diminished for the adaptive solution. The left plot of Figure 6.6 shows a typical mesh of the adaptive IFEM, and we observe that the mesh is densely refined around the interface as well as the point of singularity. The right plot of Figure 6.6 shows the optimal-rate decay of the errors and the estimators of our adaptive IFEM. Again, the averaging effectivity index is close to 3, which is similar to that in previous examples. This indicates the uniform effectivity of the error estimate with respect to the type of elements, i.e., interface and non-interface elements. The comparison of error in the right plot of Figure 6.6 shows a stronger superiority of the adaptive mesh refinement than the uniform mesh refinement. In fact, the convergence of IFEM with uniform mesh refinement is not optimal, due to the singularity of the solution.
The numerical solutions and error surfaces for the adaptive IFEM and uniform IFEM are depicted in Figure 6.7 and Figure 6.8, respectively. It is easy to see that the Example 6.4 (Solution with complicated interface shape). In this example, we consider an interface problem with a more complicated interfacial shape. The exact solution has a petal-shaped interface and it is defined through the following level set function: Due to the complexity of interface shape, we start the AMR procedure with a finer initial mesh, a 16 × 16 Cartesian triangular mesh. A typical mesh is depicted on the left plot of Figure 6.9. Comparing with Example 6.1, in which the interface is an ellipse (see Figure 6.1), the refinement around interface is denser. This is because the larger curvature of the interface causes the larger value of the inconsistency term in the error indicator in (3.1).
The convergence plot depicted in the right plot of Figure 6.9 indicates the optimalrate decay for both the errors and estimators. The errors of adaptive solution are a little smaller than the errors of the uniform solution with similar degrees of freedom, although the latter also converge in optimal rate. Moreover, the efficient index for this example is close to 3 which is similar to those in the previous examples. Numerical solutions and error surfaces of uniform and adaptive IFEMs are reported in Figure  6.10 and Figure 6.11, respectively. We can again see that the errors of adaptive IFEM are smaller than the errors of uniform IFEM given similar degrees of freedom. Example 6.5 (Effect of the inconsistent error terms). In this example, we will explore the effect of inconsistent error term for different interface geometries by revisiting the Example 6.1 and Example 6.4 which has a simple ellipse interface and a more complicated petal interface. We consider the following error indicator First we compare the convergences using the error estimator ξ and η for Example 6.1. In the left plot of Figure 6.12, the estimators ξ and η are very similar at all meshes, and the errors of the corresponding IFEM solutions guided by these two estimators are also close. We believe that due to the simple and smooth interface shape of this example (an ellipse interface), the inconsistent error term is negligible. Next, we use the new error estimator ξ in Example 6.4 In which the geometry of interface is more complicated. As we can see in the right plot of Figure 6.12, the error estimators η and ξ and the corresponding IFEM solutions show notable differences, especially on the first few coarse meshes. In this sense, including the geometrical correction (inconsistent error) term leads to better error indication particularly on coarse meshes. We also note that as the mesh is adaptively refined, the error estimators η and ξ become closer, as well as the IFEM solutions leading by these two estimators. Example 6.6 (Additional Comments on Large Jump Scenarios). In this test, we revisit the large jump scenario using test problem from the Example 6.2. To show that necessity of performing adaptive mesh strategy for IFEM, we use the true error in energy norm as our error indicator, i.e., η * K = α 1/2 (∇u − ∇ h u T ) 0,K , and the global error estimator is the true error in the energy norm, i.e., An adaptive mesh using the error estimator η * with similar number of triangles are shown in the left plot of Figure 6.13. In both cases the refinement is concentrated around the interface, which is similar to the adaptive mesh in Figure 6.3. This again shows that the partially penalized IFEM itself may not be sufficient to obtain accurate solution for extremely large jumps of the coefficient. In this case, the adaptive mesh refinement is more advantageous. The right plot of Figure 6.13 shows the convergence of the errors governed by the estimator η and the true error η * . We can see that both converge in an optimal rate, although using the true error as estimator gives slightly more accurate solutions.