A sharp interface method using enriched finite elements for elliptic interface problems

We present an immersed boundary method for the solution of elliptic interface problems with discontinuous coefficients which provides a second-order approximation of the solution. The proposed method can be categorised as an extended or enriched finite element method. In contrast to other extended FEM approaches, the new shape functions get projected in order to satisfy the Kronecker-delta property with respect to the interface. The resulting combination of projection and restriction was already derived in Höllbacher and Wittum (TBA, 2019a) for application to particulate flows. The crucial benefits are the preservation of the symmetry and positive definiteness of the continuous bilinear operator. Besides, no additional stabilisation terms are necessary. Furthermore, since our enrichment can be interpreted as adaptive mesh refinement, the standard integration schemes can be applied on the cut elements. Finally, small cut elements do not impair the condition of the scheme and we propose a simple procedure to ensure good conditioning independent of the location of the interface. The stability and convergence of the solution will be proven and the numerical tests demonstrate optimal order of convergence.

tion, diffusion problems and elasticity. The application of our interest is the motion of incompressible immiscible fluids.
Interface problems can be discretised using fitted or unfitted finite element spaces. In a fitted approach, the variational formulation is based on the continuous bilinear form and therefore inherits all its properties. Especially stability and optimal order of convergence can be guaranteed. However, in case of moving interfaces, the mesh generation process becomes expensive.
In order to avoid expensive re-meshing the immersed boundary method (IBM) has become a popular method to ease the representation of the geometry by allowing the interface to cut the elements. The original IBM was formulated by Peskin [28] for fluidstructure interactions. New Lagrangian points (LP) are defined along the interface. Since the singular forces arising at the interface can be understood as delta functions, Peskin introduces suitable discrete delta functions (DDF) for the Lagrangian degrees of freedom which spread the forces into the surrounding domain. The discontinuity will be smoothed out and therefore one drawback is a minor accuracy. In [30,31] a generalized IBM with higher order approximation is developed. The discrete delta function acts as a link between the moving interface and fixed Eulerian grid. Since it employs explicit expressions for the body force, it is categorized into the class of direct forcing schemes and can be interpreted as fitted method. A good overview is given in [25].
An alternative approach are unfitted finite elements. As for the IBM of Peskin, the Eulerian mesh does not resolve the interface. By introducing local modifications, enrichments or extensions of the finite element spaces, the features of the interface can be approximated properly. The most prominent example is the extended finite element method (XFEM) first introduced by Moës et al. [26] to model the discontinuities without adapting the grid to the interface. In [13,29] XFEM is applied to flow problems in order to model the discontinuity of the pressure. In the recent work of Kirchhart et al. [22] theoretical analysis of the XFEM is applied to the interface Stokes problem. Other developements include the generalized finite element method [21] and the unfitted Nitsche method by Hansbo and Hansbo [15], also called cut finite element method (CutFEM). In contrast to the early XFEM approaches slightly different enrichment functions are utilised for the description of the discontinuity at the interface and Nitsche's method, cf. [27], is applied to impose the interface conditions weakly. Further penalty terms need to be introduced to stabilize the system. CutFEM was first derived for the elliptic interface problem by Hansbo and Hansbo [15], later Becker et al. [3] and Hansbo et al. [17] developed the method for a Stokes interface problem and in [6,16] the weak coupling approach was applied to fluid-structure interaction. In [6] the problem of pressure oscillation at the interface is treated by introducing penalty terms. A drawback of these methods is that the structure of the operators becomes more complicated due to additional enrichment. Furthermore, some enrichments produce ill-conditioning due to "small" cut elements for which the volume on one of the parts of the cut element gets very small and equally the support of the associated shape function. Additional penalisation terms yields operators which differ from the original elliptic operator. Since in the literature the terms XFEM and CutFEM are not perfectly distinguishable we will use the term CutFEM because we mainly refer to the extended finite element spaces as defined in [15]. One focus of this article is the comparison of differently modified local finite element spaces.
The immersed boundary method derived in this work lies in-between fitted and unfitted methods. After enriching and modifying the local finite element spaces in a proper way the proposed spaces on cut elements are conforming to the interface. This resembles an enrichment by discrete delta-like functions. But in contrast to Peskins approach they can be embedded into the function space of the Eulerian mesh. The result is a fitted method and the essential benefit of the enrichment is that desireable porperties of the continuous bilinear form like symmetry and positive definiteness are preserved. In particular, stability results can be derived easily without the necessity to introduce additional penalisation terms. Furthermore, our enrichment can be interpreted as a standard Galerkin scheme on an adaptively refined mesh. This perspective enables to apply the standard integration schemes which simplifies the assembling procedure on cut elements. Finally, a good condition of the discrete system can be provided by slight changes dependent on the location of the interface. But in contrast to other extended methods, where the conditioning is impaired on small cut elements, our enriched spaces features a natural stabilisation on cut elements with small support of the shape function. Our enrichment can be summarized as projection onto a space which inherits the interface under the condition to form a partition of unity (PU) which yields a reduction of the shape functions surrounding the interface. The projection and the reduction are therefore the key ingredients of the proposed immersed boundary method.
In an earlier work, c.f. [19], the construction of the proposed enrichment spaces were derived via a vertex-centered finite volume method (FVM). For the formulation of a vertex-centered FV (finite volume element method, box method) usually a second, dual mesh is introduced. Using first-order trial functions on the primal mesh (given triangulation) and piecewise constant test functions on the corresponding dual mesh the FV scheme can be formulated in variational form (Petrov-Galerkin method). Therefore we will further refer to it as a Petrov-Galerkin finite volume method (PG-FVM). Expoliting the Petrov-Galerkin interpretation of this method enables the development of the associated enriched finite element scheme. As a consequence, our method is applicable to FE schemes and to vertex-centered FV schemes. It should be emphasized that another common FV discretisation, the so-called cell-centered FV scheme, follows another approach. The technique which will be derived in this paper refers to the vertex-centered scheme. Since one focus of this article is the comparison of differently modified local finite element spaces, we emphasize that also the DDF can be interpreted as an extension of the function spaces used for the approximation of the solution.
The paper is organized as follows: In Sect. 2 we introduce the enriched spaces for a general elliptic interface problem of diffusion with discontinuous coefficients, since it likewise reflects many interface problem of interest. A suitable finite element (projFEM) and finite volume scheme (projFVM) will be derived. In Sect. 3 we compare the spaces with those of CutFEM and DDF-IBM. In Sect. 4 the properties of the enriched spaces are described in more detail. In Sect. 5 consistency, symmetry, stability and convergence of the scheme are proven. Finally, we present numerical results in Sect. 6 for the projFVM.

Mathematical formulation
For the formulation of our immersed boundary method we consider the elliptic equation with discontinuous coefficients [u] = 0 on Γ , on a polygonal and convex domain Ω : (Fig. 1). The jump condition on Γ is given by [α∇u · n] := (α 1 ∇u 1 − α 2 ∇u 2 ) · n with normal vector n on Γ pointing from Ω 1 to Ω 2 . We assume α j > 0 to be constant in each subdomain Ω j , j = 1, 2. We want to emphasize that, with regard to the applications of interest, we focus on moderate jumps of the coefficients. Additional treatment might be necessary for α 1 /α 2 1 or 1. Further, we assume that the interface is sufficiently resolved by the mesh. In [7,15,29] different assumptions are stated to similarly assure this property.
In [19] we derived an enrichment of the finite element spaces on cut elements which yields conforming finite elements w.r.t. the interface. Therein the scheme was applied to particulate flows. For this application the correct discretisation of the forces as a vectorial measure acting at the particle interface is crucial. The main motivation therefore was the consistency of the gradients of the shape functions instead of the their values. In this paper, we want to apply this gradient-consistent enrichment to the simpler Laplace equation with embedded interface. For this application, the second feature of the projected spaces being a fitted enrichment will be more important. We also compare the proposed enrichement with the extended spaces of the CutFEM and the DDF-IBM of Peskin.

Locally enriched finite element spaces
We will derive a conforming finite element formulation of the interface problem already defined for the application to particulate flows in [19]. Therefore, we define enriched spaces on the elements cut by the interface.
We first introduce some notations for quantities related to the mesh and the immersed interface. Let T h be a shape-regular, simplicial triangulation of Ω and Γ ⊂ Ω the interface which usually does not coincide with the boundary of the elements of T h . Let X h,Γ := T ∈T h ∂ T ∩ Γ be the set of all intersecting points of Γ with the edges of elements. The (d-1)-dimensional convex hull of all points of X h,Γ yields a piecewise planar approximation Γ h of Γ . Let further Ω 1,h and Ω 2,h be the corresponding subdomains of Ω satisfying ∂Ω 1,h ∩ ∂Ω 2,h = Γ h . We introduce the set of all cut elements T h,Cut := {T ∈ T h : T ∩ Γ = ∅}. By means of Ω 1,h and Ω 2,h all cut elements T ∈ T h,Cut can be decomposed into the parts lying on either side of Γ h , i.e. T 1 h,Cut := {T ∩ Ω 1,h : T ∈ T h,Cut } and T 2 h,Cut := {T ∩ Ω 2,h : T ∈ T h,Cut }, respectively (see Fig. 2). Based on the sub-elements of T 1 h,Cut and T 2 h,Cut we define the new grid T * h can be interpreted as a special non-regular refinement of T h with respect to Γ . In order to construct appropriate function spaces we further introduce some notations for the different set of vertices. Let X h and X * h be the set of all vertices of T h and T * h , respectively, and X h,Cut the vertices of all cut elements. The set of vertices X h,Cut can be decomposed into the set X h,Γ of vertices on Γ h and those being part of X h and belonging to either side Ω j , j = 1, 2, of the interface. This yields Following standard finite element theory we define a basis {ϕ i } i∈I h of nodal func- The resulting enriched finite element space will be denoted by V * h . It should be noted that T 1 h,Cut and T 2 h,Cut will contain not only triangles or tetrahedrons, but also quadrilaterals or hexahedrons, prisms and pyramids, respectively (see Sect. 4.3 for a detailed description of considered elements). As a consequence, the enriched space V * h on T * h will contain not only linear but also multi-linear shape functions.
Remark 1 (Non-extended enrichment of finite element spaces) We emphasize that the described enrichment does not lead to an extension of the domain across the boundary, since it introduces the degrees of freedom on the interface.

projFEM: a conforming finite element formulation of the interface conditions
We define a slightly adapted coefficient We can directly write down the projected FEM (projFEM) as follows: with In contrast to α h the interface condition (g, v h ) Γ will be evaluated along the original interface Γ and not be replaced by Γ h . The resulting finite element scheme is a Galerkin formulation on the entire domain Ω with respect to the enriched space V * h and additional right hand side (g, v) Γ . We emphasize that testing with ϕ * k ∈ V * h with k ∈ I h,Γ yields the following kind-of-weak form of the interface condition (2):

projFVM: a conforming finite volume formulation of the interface condition
We now apply the enriched space V * h in order to solve the interface problem in the PG-FV formulation. In [19] we proceeded in inverse direction and started with the vertex-centered PG-FV formulation in order to derive the enriched finite element spaces. This illustrates nicely, that because of their close relation either direction can be appropriate, depending on the purpose and application.
Choosing different trial and test spaces yields so-called Petrov-Galerkin schemes. The choice of piecewise constant test functions introduces additional boundary integrals in the weak formulation (as in the context of DG schemes). As a consequence, the resulting Petrov-Galerkin scheme is a finite volume formulation balancing along the according boundaries. Application of the enriched space V * h finally yields the following projected finite volume formulation (projFVM): with bilinear form The control volumes w.r.t T h will be denoted by B i accordingly. We emphasize that for a given simplicial mesh T h the set B h is required to form a partition of the domain Ω, but not necessarily is a simplicial decomposition or even identical with T h . In the case of the vertex-centered PG-FVM used in this work each B i is constructed around the vertex x i ∈ T h by connecting the barycenters of all neighbouring edges, faces and volumes to a convex hull enclosing x i (see Fig. 3). The indicator functions  Fig. 4a and those with respect to T * h in Fig. 4b. We emphasize that for x i ∈ X j h,Cut ⊂ X h , j = 1, 2, the support of the according shape functions ϕ * i gets reduced to the cut parts T j h,Cut , j = 1, 2 of an element and therefore the related B * i gets reduced accordingly (compare the colored control volumes in Fig. 4a and b). Consequently, the partition B * h contains the additional interface-enclosing control volumes B * k associated to the additional vertices x * k ∈ X h,Γ satisfying Γ ∩ B i = ∅, as depicted shaded in grey in Consequently, the enrichment by the grey control volumes is compatible with the enrichment V * h of the finite element test space. We can interpret the space of piecewise constant functions Π(V * h ) as the enriched test space for the PG-FVM with respect to the interface Γ .
In analogy to (7) we can state: Therefore, the pointwise interface condition (2) is replaced by equilibrated fluxes across the boundary of the control volume B * k , enclosing the part Γ ∩ B * k of the interface, as depicted in Fig. 4b

Relating the Galerkin FEM and the PG-FVM
For the special choice of control volumes as drawn in Fig. 3 the resulting Petrov-Galerkin scheme is known to satisfy on all simplicial elements T ∈ T * h with accordingly piecewise linear functions u, ϕ i ∈ P 1 (T * h ) and normalized test function ϕ i (x i ) = 1. This identity was first proven by Bank and Rose [1] for the two-dimensional case and later also for arbitrary dimension by Chen [8], Xu and Zou [35] and Hackbusch [14]. Therefore, for piecewise linear spaces on simplices it reproduces the stiffness matrix for the Laplace (assuming piecewise constant coefficients on each triangle). We write down the following, generalized identification between a FE : with Rest(h) = 0 due to (10) in the case of the original, simplicial mesh T h . In [36] the estimate Rest(h) = O(h 2 ) was derived for rectangular meshes. The only requirement is an approximation of α(x) by a piecewise constant representation. In the works of Ye [36] and Chou and Kwak [9] that strong relation serves for the analysis of finite volume schemes for the Stokes equations. In [19,20] we similarly exploited this identification for the formulation of the FEM for particulate flow. The identification (11) then yields a FE (u, ϕ i ) ≈ a FV (u, χ i ) for all i ∈ I * h and with regard to (7) and (9) we can particularly state for all k ∈ I h,Γ that with equality on simplicial elements. In other words, testing on a volume by integrating on the support ω k := supp(ϕ * k ), k ∈ I h,Γ can be replaced by balancing along a surface by integrating along the boundary ∂ B * k of that special interface-enclosing control volumes of the PG-FVM in (8). We emphasize the relation B * k ⊂ ω k . We want to mention the slight difference in the right hand sides in (7) compared to (9). For simplicial meshes we get identity of the right hand sides for f ∈ P 0 (T h ), see e.g. [4], but in general both may differ.
Due to the construction of x * k ∈ X h,Γ the associated B * k in particular satisfy the relation independent of the original grid T h , see also Fig. 4b. Finally, condition (3) is satisfied as well, since the shape functions ϕ * k , k ∈ I h,Γ are continuous across the interface. We emphasize that their gradients have a jump along the discrete interface Γ h ∩ B * k , not along Γ and together with the adaption of α h the solution inherits the discontinuous gradient across Γ h instead of Γ . In Sect. 5 we prove that these adaptions are still consistent.

Comparison of projFEM/projFVM with CutFEM and DDF-IBM
Since the immersed boundary method derived in this work shares ideas of the IBM by Peskin [28] as well as the extended finite element techniques of CutFEM we will compare the according function spaces. First, we emphasize the following: A DDF-IBM introduces new nodes into the system. But since explicit expressions of the boundary conditions are formulated as force terms, these nodes are not treated as independent degrees of freedom. In contrast, XFEM techniques as CutFEM introduce additional shape functions on the cut elements and thereby new degrees of freedom. Since the shape functions are commonly defined w.r.t. the underlying Eulerian grid, the new degrees of freedom are not located directly on the immersed interface but the domain gets extended across the interface. Main motivation of both strategies is to maintain the standard discretisation scheme on the background mesh.
For a proper comparison with CutFEM and DDF-IBM we construct the enriched space V * h anew and apply three distinct steps starting with the original Eulerian space V h .
Step 1-Projection Define the set of new nodes as Lagrangian points (LP) of T h with respect to Γ (see the blue vertices on Γ in Fig.  5a). The resulting discrete interface Γ h defines a suitable projection of Γ onto the grid T h .
Step 2-Embedding For all k ∈ I h,Γ define functions δ k associated to each point x k ∈ X h,Γ forming the Lagrangian space L Γ := span{δ k } k∈I h,Γ . For the sake of stability and consistency it is reasonable to embed the Lagrangian space L Γ into an appropriate Eulerian (finite element) space with respect to T h . Since X h,Γ is the set of intersections of Γ with the edges of triangles T ∈ T h , the definition yields a natural embedding L Γ ⊂ V * h into the enriched V * h which differs from the original space V h only on cut elements.
Step 3-Reduction We request that the enriched space forms a partition of unity (PU) (see also Sect. 4, Remark 2). Therefore, the basis functions ϕ i ∈ V h in the nearinterface nodes x i ∈ X j h,Cut , j = 1, 2, which share support with the functions δ k , need to be projected onto corresponding shape functions with respect to the adapted mesh T * h . We end up with reduced shape functions satisfying ϕ * i (x)| Ω j,h ≡ 0 for all i / ∈ I j h,Cut , j = 1, 2, see Fig. 7a. We denote by   Eulerian space with respect to T h and Γ . By means of the Lagrangian shape functions δ k ∈ L Γ and the reduced shape functions ϕ * i ∈ V − h we finally define the enriched finite element space as the direct sum Remark that for the reduced space we in general obtain V − h ⊂ V h due to the bilinear cut elements.
Based on these three sub-steps we can compare our approach with two common immersed interface methods.

XFEM: embedding without reduction and projection
In the CutFEM approach of Hansbo and Hansbo [15] the discontinuity across Γ gets introduced into the local space of each cut element by doubling its degrees of freedom, see Fig. 5b. This is equivalent to extending each domain Ω j across Γ -presuming the smoothness of the inner boundary. Subsequently, the two function spaces are restricted onto Ω 1 and Ω 2 which reconstructs the discontinuity. The doubling enables the retainment of the shape functions on the original mesh T h and we obtain the locally extended space by V ext with the according restriction operators R j : L 2 (Ω) → L 2 (Ω j ), j = 1, 2, see Fig. 6b. By construction, the space V ext h forms a PU, which is important to assure mass conservation, see Remark 2. However, the drawback of that construction is an asymmetry on the restricted elements R j (T ). Re-symmetrisation and stabilisation become necessary. Besides that, the gradient near the interface will remain unchanged independent of the location of the interface, see also Fig. 7b. For the detailed discussion of that observation we refer to Sect. 4.1.

DDF-IBM: Projection without embedding and reduction
In common DDF-IBM the unadapted discrete scheme is employed on the background Eulerian grid. In addition, LPs {y i } N i=1 are introduced, usually distributed equidistantly on the immersed boundary, see Fig. 5c, together with the associated discrete delta Fig. 6c. Defining the Lagrangian space L Γ := span{δ i } N i=1 of all discrete delta functions we get a representation of the whole space as We want to emphasize that the LPs commonly do not serve as distinct degrees of freedom but as interpolation functions for prescribed forces which get distributed from the interface onto the Eulerian grid points. This technically yields V DDF h = V h . Moreover, the discrete delta functions are not necessarily chosen in accordance with the shape functions on the Eulerian grid, i.e. δ i / ∈ V h . Most DDF-IBM approaches even do not define shape functions on the Eulerian grid since DDF-IBM is often applied for finite difference schemes on cartesian grids. The embedding into a common space is not intended. In addition, the set of functions V h × L Γ does mostly not form a PU. The DDF are defined to provide a good approximation of the forces rather than of the solution itself.

Properties of the enriched space
The projection and reduction step are the building blocks for the construction of T * h when deriving it from the originally given mesh T h . The direct consequence is the Kronecker-delta property of the enriching shape funcitons w.r.t. the interface. Its implications for the properties of the discrete scheme will be explained in the following.
Furthermore, the enriched spaces by construction form a partition of unity. In particular, projection and reduction are mutually dependent when claiming the partition of unity property, i.e. projection entails reduction and vice versa. The following remark emphasizes the importance of requesting partition of unity.

Remark 2 (Partition of Unity) The PU-property
is a necessary and sufficient condition to guarantee mass conservation, see e.g. [23]. The standard finite element shape functions defined on the original mesh T h satisfy the PU property. After the introduction of the shape functions in the projected LPs the reduction step becomes necessary in order to preserve the PU property.

The Kronecker-delta property: consistency of the gradients
For application to particulate flows in [19] the construction of the enriched spaces by means of projection and reduction was driven by the requirement of consistent gradients of the shape functions near the interface. In the context of fluid problems directions of forces at the interface play a crucial role and these were shown to be related to the gradient of the test space. For more details we refer to [19].

The Kronecker-delta property: natural stabilisation by steeper gradients
With regard to scalar elliptic problems we want to point out the impact of the projection and restriction on the stability. Since the enriched spaces can be interpreted as a special non-regular, interface-adapted refinement, the resulting finite element spaces are conforming. As a consequence, the formulation is based on the continuous bilinear form and inherits its nice properties: First, the symmetry is maintained and not impaired by additional penalisation terms. Second and more important, the stability of the bilinear operator can be easily guaranteed. In fact, for good approximation properties and conditioning the shape regularity of T * h is the only criterion that needs to be ensured, see Sect. 5.1. And in particular, small cut elements are not an issue. If the area of one part of the cut element is small, the small support of the shape function results in an almost-dependency between the degrees of freedom if the enrichment is based on shape functions w.r.t. to the original mesh. In contrast, V * h provides a natural stabilisation on small cut elements because of its gradients: For linear finite elements it is |∇ϕ i | T = 1/h i (with h i being the height of T w.r.t. x i and the opposite edge) and therefore they scale inversely proportional to the area of the element. For illustration, we consider the two dimensional case of a triangle T ∈ T h of the original mesh. The local coupling matrix A ∈ R 3×3 , with A i j := a(ϕ i , ϕ j ), is given as follows: Considering the simplified case that the interface cuts the two edges of T in the same relation, i.e. Γ h || e i for one i ∈ {1, 2, 3}. For the corresponding objects on the cut part T this yieldsẽ i = λe i andh j = λh j for some 0 < λ < 1. Hence, it isẽ i /h i = e i /h i andÃ i j = A i j . Therefore, the condition of the matrix is not impaired because of small support. For the general case we still getÃ i j ≈ A i j . Therefore, the shape regularity of a mesh is the only condition for the discrete system to be well-conditioned. A criterion for a suitable adaption of the scheme will be derived in Sect. 5.1.

Existence of a local finite element space on cut elements
Appropriate nodal shape functions for x * k ∈ X h,Γ with its support on the elements T ∈ T j h,Cut , j = 1, 2, need to be provided. For piecewise linear finite element spaces V h := P 1 (T h ) on simplicial triangulations suitable local spaces can be defined quite easily. In two dimensions the simplex is a triangle which gets cut into two triangles or a triangle and a quadrilateral (see Fig. 5a) on which correspondingly a linear and bi-linear ansatz exists, see Fig. 6a. In three dimensions the simplex is a tetrahedron and gets cut into combinations of tetrahedra, prisms and pyramids, see Fig. 8. For all these elements local shape functions of second order exist. The details on the local spaces used for our computations can be found in [32]. Consequently, for T j ∈ T j h,Cut , j = 1, 2, we get the optimal approximation error [5] with constant c > 0 and h T j being half the diameter of the triangle T j . Our approach certainly will afford more effort for the definition of suitable local spaces on the cut elements, if the mesh contains also quadrilaterals or octahedrons in two or three dimensions, respectively. For all considered applications the simplicial grid was appropriate.

Cut elements as new reference elements
Because of the reduction step the discretisation w.r.t. the original mesh T h can not be applied for the near-interface nodes. However, by the use of reference elements for the assembling process the usual assembling can be applied: Algorithms for finite element schemes on unstructured meshes usually exploit the possibility to map each (unstructured) element of the mesh onto a reference element. The assembling of the discrete the interface; 3 edges (top) and 4 edges (bottom) cut by the interface.
interface; 2 opposite edges cut by the interface.  With regard to Sect. 4.3 we can assume that for the sub-elements the discretisation on the reference element is well defined. As a consequence, the same standard routines for unstructured elements can be applied to the assembling on each cut element T ∈ T h,Cut . This simply requires calling these routines twice and with respect to the corresponding new coordinates. Moreover, we maintain local stencils and optimal convergence properties as will be shown in Sect. 5.3. We emphasize, that in particular no special integration rules on the cut elements (due to e.g. restriction and change of support) need to be defined.

Numerical analysis
In this section we will derive theoretical results for the projFEM in (6). Unfortunately, due to R(h) = 0 for the space V * h , the theoretical results can not directly be transferred to the projFVM in (8). However, since the enriched mesh T * h contains only a few multilinear elements along the interface compared to the global number of elements, it is reasonable to assume that Rest(h) = O(h 2 ) also holds for the projFVM, provided that T * h again is shape regular. We therefore already state for the projFVM in (8) that the results derived in the following subsections, based on finite element theory, together with the identification (11) at least allow us to expect a similar convergence behaviour for the projFVM. The numerical results presented in Sect. 6 apply the projFVM to elliptic problems and confirm our reasoning.

Conditioning
In the case of a moving interface the elements of T h can be cut in arbitrary ways. Consequently, the shape-regularity of the mesh T * h is not guaranteed (i.e. the existence of a constant κ > 0 independent of the location of the interface, s.t. h T /ρ T ≤ κ, with h T and ρ T being the circumcenter and incircle-radius of an arbitrary T ∈ T * h ). That is a potential source of ill-conditioning of our discrete scheme. In order to preserve shape-regularity we simply re-define those points x * k ∈ X h,Γ whose distance to a point on the Eulerian mesh falls below a certain threshold. This way, also the approximation property will not be impaired. To formulate a suitable criterion for a re-definition we choose a lower bound D > 0 satisfying D < h min , with h min := min{h T : T ∈ T h } being the minimal circumcenter in the original grid. The re-definition then reads: The adapted triangulation T * h satisfies the condition h T /κ ≤ ρ T . The proof can be found in [18]. The resulting discrete interface Γ h approximates Γ less accurately. However, for the choice of D := h 2 < h min , with h being the mean circumcenter of the given grid, we can conclude that the resulting piecewise planar interface In Sect. 5.3 we will show that the corresponding discrete scheme maintains optimal approximation order. A higher condition number in favor of a stable solution process for the linear system appears to be a valid strategy confirmed by the numerical simulations in Sect. 6. In common XFEM approaches the ill-conditioning of the discrete system is similarly an issue. But in contrast to the shape-regularity above, already small but shape-regular cut elements are critical, see Sect. 4.2. There we also explained that the stabilising effect of the steeper gradients of the projected and reduced shape functions (which in particular are not defined w.r.t. the original mesh). The Hansbo-averageing is a general strategy to handle the issue of ill-conditioning due to small support for the CutFEM spaces. We want to mention it in order to contrast it to the proposed re-definition in (13): Remark 3 (Hansbo-averageing) Formulations based on Nitsche's method [27] introduce an averaging operator {·} to re-weight the shape functions on the cut elements by with weights κ 1 and κ 2 satisfying κ 1 + κ 2 = 1 in order to preserve consistency of the discrete scheme. For T j := supp(ϕ j ) ∩ T , j = 1, 2, the so called Hansbo-averaging, cf. [15], defines κ j := T j /T . Hence, the weights depend on the size of the relative area of each sub-element which resembles the parameter d i, j in (13). Other weighting strategies are defined in [17,29,33] likewise depending on the properties of the subelements.
Remark 4 (Neglection vs. Redefinition) A weighting rule which defines κ 1 := 0 and κ 2 := 1 in (15) in the case of very small support, i.e. T 1 /T ≈ 0, corresponds to neglecting the presence of the interface on the cut element T . Optimal approximation of the solution can be preserved if a criterion for neglecting is defined suitably, cf. [29]. It should be emphasized that the described re-definition in (13) resembles the strategy of neglecting but does not correspond to it: Consider the case of a two-dimensional triangle cut by the interface in such a way that after application of criterion (13) the whole triangle lies on one side of the interface. Therefore, the triangle will not be part of the enrichment. In contrast to neglecting the interface, the re-definition of the LP leads to a slight displacement of the interface. In particular, also on the neighbouring elements the shape of Γ h gets displaced. The essential consequence is the continuity of the interface across neighbouring cut elements.

Coercivity and symmetry
The coercivity of a FE (·, ·) on V * h directly follows from the coercivity of the continuous bilinear form. That directly implies the important advantage of our approach that no additional stabilisation is needed which can cause numerical fluxes across the interface. With regard to the coercivity of a FV (·, ·) on V * h , we state the following: Remark 5 (Coercivity of a FV (·, ·)) Due to (11) the coercivity of the finite volume scheme is potentially impaired only on the multi-linear elements. Based on our assumption that Rest(h) = O(h 2 ) there exists a h FV > 0 and an according constant C(h FV ) such that the coercivity is satisfied with constant C(h FV ) and for all h ≤ h FV .
A further advantage of the proposed projFEM and projFVM is the avoidance of the weak boundary condition term ([α∇u · n], v) Γ so that the symmetry of the scheme is preserved. Since the original problem is self-adjoint, the lack of symmetry and therefore adjoint-consistency can lead to impaired numerical results. For that reason Nitsche's method retrieves symmetry by introducing suitable additional terms.

A priori error analysis
In order to derive an a priori error estimate for problem (5) we in particular need to bound the approximation error. The shape-regularity is guaranteed by means of 13. This enables to use the standard nodal interpolation operator I h : For all elements T / ∈ T h,Cut being not cut by the interface, we get the optimal approximation error [5] with constant C > 0 and h T being half the diameter of the triangle T , see also (12).
For a complete a priori analysis two difficulties need to be addressed: the solution u is not smooth across the interface Γ and due to Γ h ⊂ Γ the piecewise linear, discrete interface is not aligned with the real interface. It is therefore necessary to derive appropriate interpolation estimates on the according elements T ∈ T h,Cut .

Lemma 1 (Interpolation Error on
Proof Let T c ∈ T h,Cut . By means of 16 it is sufficient to derive an estimate for ∇(u − I h u) 0,T c . The proof is based on the reasoning in the work of Frei and Richter [11], Frei [10] and Basting and Prignitz [2]. Their proofs rely on a continuous extensionũ i ∈ H 2 (Ω) of u ∈ H 2 (Ω i ), i = 1, 2. Such an extension exists, see [34] and satisfies Insertion of ±ũ i and ± I hũi yields Since all Lagrange nodes lie on the interface, it is I hũi = I h u and therefore the according term of the insertion vanishes. To derive an estimate for ∇(u −ũ i ) 0,T c we can apply the same argumentation as in [10,11]. Herein, they use a Poincare-like estimate to bound ∇u 0,T c . Combined with the continuity of the extension (17) this yields [10,11] To derive an estimate for the interpolation error ∇(ũ i − I hũi ) 0,T c we can use the standard estimate stated in (12), since it is T c = T 1 ∪ T 2 , with T 1 ∈ T 1 h,Cut , T 2 ∈ T 2 h,Cut . Together with an enlargement onto Ω and again the continuity of the extension (17) we get Combining the last estimate with (18) and (19) completes the proof.
The final discretisation error bound can be derived by combining the interpolation error stated in Lemma 1 with the consistency error. That is the common procedure in Strang's lemmata. Let u ∈ H 2 (Ω 1,2 ) and u h ∈ V * h be solutions of problem (5) and (6), respectively. For arbitrary, but fixed coefficients α 1 and α 2 let further c α be the according coercivity constant, i.e.   Second order convergence for the interface problem (24) In the work of Frei and Richter [11], Frei [10] and Basting and Prignitz [2] slightly different discrete spaces are suggested, where they similarly introduce Lagrange nodes on the interface and derive the same error estimates. Furthermore, for the CutFEM spaces, comparable approximation error bounds have been derived in Reusken [29] and Hansbo and Hansbo [15].

Numerical experiments
In this section we apply the projFVM introduced in Sect. 2 to the problem (1)-(4) for different test cases.

Experiment 1
This test case was already investigated by Liu et al. [24]. The computational domain is the unit square Ω = (−1, 1) 2 and the interface Γ is defined as the circle with radius R = 0.5 and center x m = (0.0, 0.0). As diffusion coefficients we choose α 1 = α 2 = 1. The analytical solution for Δu = 0 is given by The Dirichlet data is defined according to the exact solution u(x) yielding the non-zero jump condition g(x) = 2. The L 2 -errors obtained under global mesh refinement are depicted in Fig. 9a. We get a convergence of order O(h 2 ) in accordance with Theorem 1.

Experiment 2
This test case was already investigated with variations by Frei and Richter [11] and Gangl and Langer [12]. The computational domain is the unit square Ω = (−1, 1) 2 and the interface Γ is defined as the circle with radius R = 0.4 and center x m = (0.0, 0.0). As diffusion coefficients we choose α 1 = 10 and α 2 = 1. The analytical solution is given by The right hand side and the Dirichlet data are defined according to the exact solution u(x). Remark that the weak discontinuity in the solution along Γ combined with the jumping coefficients yield a jump condition equal to zero, i.e. g(x) = 0. The L 2 -errors obtained under global mesh refinement are depicted in Fig. 10a. We get a convergence of order O(h 2 ) in accordance with Theorem 1.
The same experiment was performed with shifted center x m = (−0.08, 0.3) and same radius R = 0.4. A comparable L 2 -error was obtained and is depicted in Fig. 11a.

Conclusions
We proposed enriched finite element spaces which yield an interface fitted scheme for elliptic interface problems. Therefore, the derived method provides the advantage of an immersed interface method avoiding costy remeshing. In [19] the enrichment was derived for particulate flows. The motivation was a correct discretisation of the forces acting on the boundary. Therein, the focus was the gradient of the enriching shape functions. In this paper we could show that applied to a general elliptic interface problem we get the same order of approximation as comparable enriched methods. But in contrast to those, the discrete operator even herits the symmetry and positive definiteness of the continuous operator since it is a fitted method.
We further compared our method with the enrichment of a CutFEM and the idea of Peskins original IBM by looking at the local finite element spaces used in the according approach. Instead of extending and restricting the domain at the interface for the CutFEM, we apply as we call it projection and reduction. By opposing the shape functions on cut elements we could argue, that as already for the application in [19] the gradient of the enriching shape functions has a crucial role. Its slight change depicted in Figs. 6 and 7 inherits a natural stabilisation on small cut elements.
In summary, a deeper insight into the properties of the enriching shape functions could provide a physically meaningful formulation of the interface problem with good numerical properties and without the need for additional, usually non-physical penalty terms.
In the case of deformable interfaces a good approximation of an immersed interface and the forces due to gradients becomes even more relevant. Therefore, the next step will be the application of the derived approach to multiphase flows.