The Algebraic Immersed Interface and Boundary Method for Elliptic Equations with Jump Conditions

A new fictitious domain method, the algebraic immersed interface and boundary (AIIB) method, is presented for elliptic equations with immersed interface conditions. This method allows jump conditions on immersed interfaces to be discretized accurately. The main idea is to create auxiliary unknowns at existing grid locations which increases the degrees of freedom of the initial problem. These auxiliary unknowns allow to impose various constraints to the system on interfaces of complex shapes. For instance, the method is able to deal with immersed interfaces for elliptic equations with jump conditions on the solution or discontinuous coefficients with a second order of spatial accuracy. As the AIIB method acts on an algebraic level and only changes the problem matrix, no particular attention to the initial discretization is required. The method can be easily implemented in any structured grid code and can deal with immersed boundary problems too. Several validation problems are presented to demonstrate the interest and accuracy of the method.


Introduction and general motivations
Simulating flows and heat transfer interacting with complex objects on Cartesian structured grids requires an efficient coupling between such grids and the corresponding numerical methods and complex shape interfaces. Such a coupling is often performed thanks to fictitious domain methods, where the computational domain does not match the physical domain. The advantages of this approach are numerous. A second-order accurate discretization of the spatial operators is simple to obtain, grid generation is trivial, and furthermore there is no need to remesh the discretization grid in the case of moving or deformable boundaries. Concerning this last point, fictitious domain methods can be useful even on unstructured grids: Eulerian fixed unstructured grids can fit immobile obstacles, (e.g. a stator of an aircraft motor) while mobile objects (a rotor) are treated with fictitious domain methods. Two particular classes of problems can be drawn: the immersed boundary problems and the immersed interface problems. The firsts deal with complex boundaries, such as flow past objects, where no attention has to be paid to the solution inside the obstacles. The immersed interface problems consider subdomains delimited by interfaces, and the solution is required in both sides of the interface. As particular conditions, such as jump conditions, can be required on the interface, this second class of problems is often more difficult to treat. Let us consider the following model scalar immersed boundary problem with a Dirichlet boundary condition (BC) on the interface Σ (see Fig. 1): A boundary condition is also required on the other part of the boundary ∂Ω 0 so that the whole problem is well-posed. A first approach dealing with immersed boundaries is the distributed Lagrange Multiplier method proposed by Glowinski et al. [9]. Lagrange multipliers are introduced into the weak formulation of the initial elliptic equation to ensure the immersed boundary condition. Cartesian grid [13,22] and Cut-cell [39] methods use a structured grid in the whole domain except near obstacles where unstructured cells are created from structured cells. These methods are hard to implement due to the numerous different space configurations of the intersections between cells and objects. Furthermore, the existence of small cells can induce solver troubles. The immersed boundary method (IBM) was initially presented by Peskin [24,25]. Fictitious boundaries are taken into account through a singular source term defined only near the boundaries. As the source term is weighted with a discrete Dirac function smoothed on a non-zero support, the interface influence is spread over some grid cells. This method is first-order in space and explicit. Another class of IBM, the direct-forcing (DF) method, was initially proposed by Mohd-Yusof [23]. The idea here is to impose a no-slip condition directly on the boundary using a mirrored flow over the boundary. In [5,38], the correct boundary velocity is obtained by interpolating the solution on the boundary and far from the boundary on grid points in the near vicinity of the interface. In [37], Tseng et al. use the same principle but extrapolate the solution in ghost cells outside the domain. This approach can be seen as a generalization of the mirror boundary conditions used in Cartesian staggered grids to impose a velocity Dirichlet condition on pressure nodes. As discussed in [32], this kind of approach seems to be more accurate than [5,38]. The penalty methods for fictitious domains consist in adding specific terms in the conservation equations to play with the order of magnitude of existing physical contributions so as to obtain at the same time and with the same set of equations two different physical properties. The volume penalty method (VPM) [3,2] requires the addition of a penalty term χ ε (u − u D ) in the conservation equations, such that: { −∇ · (a∇u) + χ ε (u − u D ) = f in Ω with χ |Ω0 = 0, χ |Ω1 = 1, for 0 < ε ≪ 1 (1) where ε denotes the penalty parameter which tends to 0. Hence, in Ω 1 the original equation becomes negligible and u = u D is imposed. In ( [14,15]) authors add a Darcy term µ K u to the Naviers-Stokes (NS) equations where µ is the dynamic viscosity and K the permeability. In the fluid medium, K → ∞ so the Darcy term is then negligible and the original set of NS equations is retrieved. In the solid medium, K → 0 and consequently the NS equations tend to u = 0. Classical discretizations of the penalty terms are of first order only since they consider the projected shape of the interface on the Eulerian grid to define the penalty parameters [26]. In [29,31], Sarthou et al. have discretized the volume penalty term with a second order using implicit interpolations as in [37]. This method is called the sub-mesh penalty (SMP) method and has been applied to both elliptic and NS equations. Applied to problem P b , the ghost cell immersed boundary method [37] and SMP method [29] used the first cells in Ω 1 to enhance the accuracy of the solution in Ω 0 . An other approach which considers the extension of the solution is considered in [8,7] by Gibou and Fedkiw. Ghost nodes and simple interpolations are considered, but contrary to the SMP and the IBM-DF methods, only 1D interpolations are used and the operators are rediscretized "by-hand".
Let us now consider a model immersed interface problem with jump interface conditions: A first class of method is the immersed interface methods (IIM) initially introduced by LeVeque and Li [17] and widely described in [19]. This groupe of methods use Taylor series expansion of the solution at discretization points in the vicinity of Σ to modify the discrete operators at these points. Much work has been devoted to the immersed interface method and its numerous applications, such as moving interfaces [11] or Navier-Stokes equations [16]. In [18], Li uses an augmented approach. Additional variables and interface equations are added to the initial linear system. The new variables are the values of jumps at some interface points. This method has been extended to the incompressible Stokes [20] and Navier-Stokes [12].
The Ghost Fluid Method, originally developed by Fedkiw et al. [6,21], introduces ghost nodes where the solution is extended from one side of the interface to the other side. As for IIM, the operator discretization must be modified "by-hand". Zhou et al. overcome this drawback with the matched interface and boundary (MIB) method [42,41,40] by using interface conditions to express the solution at ghost nodes with respect to the solution on physical nodes. Hence, the discretization is automatically performed whatever the discretization scheme. Contrary to [18], the additional equations for these two last methods are not written at "random" points of the interface but at the intersections between the Eulerian grid and the immersed interface. Furthermore, simple Lagrange polynomials are used whereas a more complicated weighted least squares approach is used in [18] to discretize additional equations. In [4], Cisternino and Weynans propose a quite simple method with additional unknowns located at the interface. Interfaces conditions are discretized at these points and are added to the final linear system.
The method presented in this work solves elliptic problems using an augmented method coupled with an auxiliary unknown approach. Contrary to ghost nodes, auxiliary unknowns are present in the linear system. Compact interpolations are used to discretize the additional interface constraints. The method is simple to implement even for interfaces of complex shapes, i.e. not described by analytical equations. Except for the discretization of interface conditions, all operations are automatically performed with algebraic modification or directly by the "black-box" matrix solver. This new method is called the algebraic immersed interface and boundary (AIIB) method. In section 2, the method is presented for immersed boundary problems. Then, the method is extended to immersed interface problems with known solution on the interface. Finally, the method is applied to immersed interfaces with transmission and jump conditions. A special attention is paid to the management of the discretized interface, especially the way to project it onto the Eulerian grid using a fast ray-casting method. In section 3, validation tests and convergence studies are presented. Conclusions and perspectives are finally drawn in section 4.

The algebraic immersed interface and boundary method
The AIIB method is now presented. The method is first formulated for immersed boundary problems when a Dirichlet or a Neumann boundary condition is required. The method is then extended to simple immersed interface problems where the solution is a priori known on the interface. Finally, an extension to jump and transmission conditions is described.

Definitions and notations
Let us consider the original domain of interest denoted by Ω 0 , typically the fluid domain, which is embedded inside a simple computational domain Ω ⊂ R d , d being the spatial dimension of the problem. The auxiliary domain Ω 1 , typically a solid particle or an obstacle, is such that : Ω = Ω 0 ∪ Σ ∪ Ω 1 where Σ is an immersed interface (see Fig. 1). Let n be the unit outward normal vector to Ω 0 on Σ. Our objective is to numerically impose the adequate boundary conditions on the interface Σ. These conditions will be discretized in space on an Eulerian structured mesh covering Ω. As the discretization of the interface or boundary conditions requires interpolations, the following interpolations in 2D: P 2 1 (x, y) = p 1 + p 2 x + p 3 y and Q 2 1 (x, y) = p 1 + p 2 x + p 3 y + p 4 xy are used. In 3D, we use P 3 1 (x, y, z) = p 1 + p 2 x + p 3 y + p 4 z and Q 3 1 (x, y, z) = p 1 + p 2 x + p 3 y + p 4 z + p 5 xy + p 6 yz + p 7 zx + p 8 xyz. An additional interpolation, L 1 1 (x) = p 1 + p 2 x, is also possible to be chosen for 2D and 3D problems. The superscript is the dimension of the interpolation while the subscript is the order of spatial accuracy.   Fig. 2). Our objective is to discretize Dirichlet, Neumann, transmission and jump conditions at these interface points to build a general fictitious domain approach. This method is expected to reach a global second-order spatial accuracy. We shall use the following Eulerian volume fonctions in order to implicitly locate Σ h : • The Heaviside function χ, defined as: This function is built with a point in solid method presented below. The function χ will be used to perform fictitious domain algorithms and to build a level-set function.
• The level-set function ϕ, with: and dist Σ (p) = inf x∈Σ ∥x − p∥. The unsigned distance is computed geometrically. The sign is directly obtained with the discrete Heaviside function χ.
• The colour phase functions C, which is the ratio of a given phase in a control volume. We denote C(x I ) the phase ratio in the control volume centered in x I . This function is approximated from the ϕ function by using the formula proposed by Sussman and Fatemi [36] : New sets of Eulerian points x I are defined near the interface so that each one has a neighbor x J verifying χ J ̸ = χ I (with χ I = χ(x I ) and χ J = χ(x J )), i. e. the segment [x I ; x J ] is cut by Σ h . These Eulerian "interface" points are also sorted according to their location inside Ω 0 or Ω 1 . Two sets {x I , I ∈ N 0 } and {x I , For each x I , I ∈ N 0 or I ∈ N 1 , we associate two unknowns: the physical one denoted as u I and the auxiliary one u * I .

Projection of the Lagrangian shape on the Eulerian grid
The generation of the Lagrangian mesh of the interface is achieved using a computer graphics software. Specific algorithms have been developed to project this Lagrangian grid onto the Eulerian physical grid. In order to obtain the discrete Heaviside function χ, one have to determine which Eulerian points are inside the domain Ω 1 defined by a Lagrangian surface. Such a surface must be closed and not selfintersecting. In [32,15], the authors used a global methodology partly based on [34] where χ is obtained thanks to a PDE. This method suffers from a lack of accuracy and robustness. A Ray-casting method based on the Jordan curve theorem is more adapted and is used in the present work. The principle is to cast a ray from each Eulerian point to infinity and to test the number of intersections between the ray and the Lagrangian mesh. If the number of intersections is odd, the Eulerian point is inside the object, and outside otherwise. The Ray-casting method can be enhanced by classifying elements of the Lagrangian mesh with an octree sub-structure which recursively subdivides the space in boxes. If a ray does not intersect a box, it does not intersect the triangles inside the box. A fast and simple optimization is to test if a given point is in the box bounding the Lagrangian mesh. An improvement of the Ray-casting algorithm, the Thread Ray-casting can be found in [30] and is described by Algorithm 1. Rays are cast from points x I included in a boundary slice S xy of the Eulerian mesh. For each starting point x I , the intersections are stored and sorted according to their z component in a two-entry structure S(I, nsect I ).
For each x I ∈ S xy , the number nsect I of intersections by rows, is not known a priori. If S is an array, a first pass has to be performed to determine the size of S. A better choice is to use chained lists. For

Algorithm 1 Optimized computation of the discrete Heaviside function in 3D
for The Lagrangian points x l of Σ h , l ∈ I are required to couple the Lagrangian surface and the Eulerian grid used to solve the conservation equations. These points can be obtained with two methods. A geometrical computation of the intersections gives the most accurate result. If not optimized the computational cost of this method is not always negligible for some cases.
Using the Level-set function is a faster but less accurate way to obtain the intersection points. Let us consider two Eulerian points x I ∈ Ω 0 and x J ∈ Ω 1 . We denote by d I = d(x I , Σ h ) and d J = d(x J , Σ h ) the unsigned distances between Eulerian points and the interface Σ h . Then, Algorithmic problems can be encountered if the Lagrangian mesh is too complex compared to the Eulerian mesh. For example, two intersecting points x l can be found between x I and x J with the geometric method. In this case, only one intersecting point is considered.
Concerning the use of the Level-set, this function is a projection of the shape on a discrete grid. The local curvature of the projected shape is thus limited by the accuracy of the Eulerian grid. Consequently, no more than one intersecting point can be found between x I and x J with the Level-set.

General principle
Once the shape informations are available on the Eulerian grid, the problem discretization has to be modified to take into account the fictitious domain delimited by an immersed boundary or an immersed interface. The sub-mesh penalty (SMP) method [32,29] was originally designed to treat immersed boundary problems. It could be extended to treat immersed interface problems by symmetrization of the algorithm with introduction of auxiliary unknowns as in the AIIB method. This new method is an enhancement of the SMP method which is also able to solve immersed interface problems. The main idea of the AIIB method is to embed an interface into a given domain by modifying the final matrix only. As no modification of the discretization of the operators is required (contrary to [8,7] and the immersed interface methods [17]), the AIIB method is thus simple to implement.
Let P be a model problem discretized in the whole domain Ω as Au = b where A is a square matrix of order m, u the solution vector and b a source term. The basic idea of the AIIB method is to add new unknowns and equations to the initial linear system so as to take into account additional interface constraints. The new unknowns, so-called the auxiliary or fictitious unknowns and labeled with * , are defined as being the extrapolation of the solution from one side of the interface to the other, and are used to discretize the interface conditions. Hence, the orignal problem A ′ a square matrix of order m + n, with n the number of auxiliary constraints related to the interface conditions. The solution u ′ is decomposed such as u ′ = (u, u * ) T and the source term as The interface constraints are discretized with a (n, m + n) block matrix C and the source term b * .
According to the interface conditions, the regularity of the solution on the interface is often lower than in the rest of the domain. Hence, the discretization of operators with a stencil cutting the interface can induce a great loss of accuracy. The first idea is to consider unknowns u * I , I ∈ N 1 (resp. u * I , I ∈ N 0 ) as the extension of the solution in Ω 0 (resp. Ω 1 ). The initial algebraic link between unknowns from both sides of the interface is cut, and the new link over the interface is obtained thanks to auxiliary unknowns. Practically, matrix coefficients must be modified to take into account the new connectivities. Let α I,J be a coefficient of A at row I, column J and α ′ I,J the new coefficient in A ′ . If I ∈ N 0 and J ∈ N 1 , α ′ I,J = 0 and α ′ I,J * = α I,J , where J * is the index corresponding to u * J . This is exactly the way how we proceed for the practical algorithm. However, this modification can be expressed algebraically with permutation and mask matrices as follows.
We define the two following mask matrices I 1 of dimensions (m, m + n) and I 2 of size (n, m + n) : The matrices A 0 and A 1 are defined such as Finally, the connectivities are changed using the permutation matrices P 0 and P 1 : P 0 is defined to switch row I with row J if I ∈ N 0 , J ∈ N 1 and P 1 to switch row I with row J if I ∈ N 1 , J ∈ N 0 . Hence, the new problem matrix is now defined by: The new problem is A ′ u ′ = b ′ with A ′ written with 4 blocks of various sizes:Ã(m, m), B(m, n), C 1 (n, m), C 2 (n, n). The matrixÃ is thus the modification of the initial matrix A by setting to zero the coefficient , and C 1 and C 2 are the two sub-matrices of the matrix C. The problem can be written as: The entire problem can then be solved to obtain u ′ = (u, u * ) T . However, u * being the auxiliary solution is not required to be computed explicitly . Hence, the Schur complement method can be used to calculate the solution for the physical unknowns only. The final problem is now: The opportunity of such a reduction will be discussed later.

AIIB algorithm for a scalar equation with Dirichlet boundary conditions
For sake of clarity, let us first describe in 2D the AIIB method for the model scalar problem P b with a Dirichlet boundary condition on the interface Σ. For this version of the AIIB algorithm, Ω 0 is the domain of interest and auxiliary unknowns are created in Ω 1 only. Let us consider a point x I , I ∈ N 1 . At location x I , two unknowns coexist: a physical one u I and an auxiliary one u * I . We first describe the case when x I has only one neighbor x J in Ω 0 . The Lagrangian point x l is the intersection between [x I ; x J ] and Σ h (Fig. 1 right). Then, the solution u l = u D (x l ) at the interface is approximated by the P 1 1 interpolation between the Eulerian unknowns u * I and u J : As noticed in [37,7], only a linear interpolation is required to reach a second order of accuracy. If now x I has a second neighbor x K in Ω 0 , the intersection x m between [x I ; x K ] and Σ h is considered with u m = u D (x m ). We choose x p , a new point of Σ h between x l and x m (see Fig. 3 left). The solution u p = u D (x p ) is then imposed using a P 2 1 -interpolation of the values u * I , u J and u K : A Q 2 1 interpolation of u I , u J , u K and u L can be also used by extending the interpolation stencil with the point x L which is the fourth point of the cell of the dual mesh defined by x I , x J and x K (see Fig. 3 left). As a third choice, two independent linear 1D interpolations can be used (one for each direction) for an almost equivalent result. It produces : In this case, two auxiliary unknowns are created. A simple choice for x p is the barycenter between x l and x m where u p = (u l + u m )/2. This particular case enables an easy implementation since we have : A summation of these two constraints gives : what is equivalent to build a constraint imposing u p at x p with a P 2 1 interpolation : Hence, an easy general implementation consists in summing the constraints corresponding to each direction, no matter the number of neighbors of x I . If the elements σ l of Σ h used to define x l and x m are not the same, the barycenter x p of these two points is not necessarily on Σ h , especially for interfaces of strong curvature. However, the distance d(x p , Σ h ) between x p and Σ h varies like O(h 2 ) and so this additional error does not spoil the second-order precision of our discretization. The convergence of this additional error is numerically tested in section (3.3.1). If the curvature of Σ h is small enough relatively to the Eulerian mesh, i.e. if the Eulerian mesh is sufficiently fine, x I almost never has a third or a fourth neighbor in Ω 0 . However, if this case appears, a simple constraint u * I = u B is used with u B being an average of u D at the neighbor intersection points. In any case, by decreasing the Eulerian mesh step h, the number of points x I having more than two neighbors in Ω 0 also decreases. Hence, the present method is suitable to impose a Dirichlet boundary condition on Σ for Ω 0 , when the solution in Ω 1 has no interest. The solution u * I for I ∈ N 1 is an extrapolation of the solution in Ω 0 in order to satisfy the boundary condition on Σ and thus is non-physical. Hence, the solution at the nodes of Ω 1 far from the interface does not impact on the solution in Ω 0 . Nevertheless, the fictitious domain approach computes a non-physical solution in Ω 1 . Correct physical values can be obtained with the initial set of equations together with a volume penalty method such as VPM [14]. The imposed solution can be analytical when possible, or an arbitrary constant value. The computational cost of this approach can be reduced by switching the solving of u I , x I ∈ Ω 1 off, or by totally removing these nodes in the solving matrix.

Symmetric version for Dirichlet interface conditions
The next step is to allow for multiple Dirichlet boundary conditions on both sides of the immersed interface. Thin objects could be treated with this approach. The problem is now: The problem (17) requires for each point x I a physical unknown u I as well as an auxiliary unknown u * I on both sides of the interface. Practically, the AIIB algorithm for a Dirichlet BC is applied a first time with Ω 0 as domain of interest, and auxiliary unknowns are created near Σ h in Ω 1 . As a second step, the Heaviside function is modified as χ := 1 − χ and the algorithm is applied a second time. Now, Ω 1 is the domain of interest and auxiliary unknowns are created near Σ in Ω 0 .

AIIB algorithm for a scalar equation with Neumann boundary conditions
Let us now consider the following model scalar problem with a Neumann BC on the interface Σ: { −∇ · (a∇u) = f in Ω 0 (a · ∇u) · n = g N on Σ The principle is about the same as for Dirichlet BC, and the same interpolations, once derived, can be used to approximate the quantity (a · ∇u) · n. Hence, at any point x l , l ∈ I on Σ h we use (a · ∇u l ) · n ≈ (a · ∇p(x l ) · n). (19) For p ∈ Q 2 1 , we get ∇p(x, y)·n = (p 3 y +p 2 )n x +(p 3 x+p 1 )n y whereas for p ∈ P 2 1 , ∇p(x, y)·n = p 2 n x +p 1 n y is obtained which means that the normal gradient is approximated by a constant over the whole support.
For example, in the configuration of Fig. 3.left, with p ∈ P 2 1 , we have: The diagonal coefficient of the raw related to u * I in C 2 is ( nx hx − ny hy ). The case where nx hx ≈ ny hy leads to numerical instabilities. If we consider the configuration of Fig. 3.left, using the normal vector of the segment [x l , x m ] implies that the signs of n x and n y are always different so the diagonal coefficient is always dominant. The same property occurs for the other cases.
When x I has only one neighbor x J in Ω 0 , the Q 2 1 and P 2 1 interpolations degenerate to L 1 1 interpolations which suit for Dirichlet BC. For Neumann BC, this loss of dimension no longer allows the interface orientation to be accurately taken into account, as one of the components of the normal unit vector disappears from the interfacial constraint. Hence, a third point x K in Ω 0 is caught to build P 2 1 interpolations (see Fig. 3

right). This point is a neighbor of x J and is taken as
In 2D, two choices generally appear, and the point being so that the angle (n, x K − x J ) is in [−π/2; π/2] is taken.

Algebraic elimination using the Schur complement
The Schur complement method allows an algebraic reduction to be performed. For a Dirichlet or Neumann BC, each constraint is written such as only one auxiliary unknown is needed: where u S is the source term. In this case, the matrix C 2 in (8) is diagonal and thus the Schur complement (Ã − BC −1 2 C 1 ) is easy to calculate. Practically, when the algebraic reduction is made,Ã is built directly by the suitable modification of A without considering the extended matrix A ′ . The part −BC −1 2 C 1 is then added toÃ whereas −BC −1 2 b * is added to b. As will be subsequently demonstrated, the algebraic reduction decreases the computational cost of the solver by 10 − 20%.
If only L 1 1 interpolations are used with the algebraic elimination, the matrix obtained with this method is similar to the one obtained in [8] for a Dirichlet problem. However, in this last paper, the auxiliary unknowns are taken into account before the discretization of the operator which requires additional calculations for each discretization scheme.
If P 2 1 interpolations are used, the computed solution in Ω 0 is the same as for the SMP [29] method (when the penalty parameter tends to zero) and the DF-IB method [37]. These methods change the discretization of the initial equation for the nodes x I , I ∈ N 1 . The SMP method uses a penalty term and the DF-IB method uses terms of opposite signs to erase some part of the initial equation. The discretization matrix obtained with both methods is not equivalent to the one obtained with the AIIB method, with or without algebraic reduction. With algebraic reduction, the discretization for the nodes x I , I ∈ N 0 is modified, and without algebraic reduction, both auxiliary and physical unknowns coexist at x I , I ∈ N 1 . The accuracy of these methods will be discussed in the next section.
The present algorithm seems simpler, as the standard discretization of the operators is automatically modified in an algebraic manner. So, various discretization schemes of the spatial operators can be used. However, the discretization of an operator at x I ∈ Ω 0 can only use in Ω 1 the fictitious unknowns and not the physical ones. Hence, the only limitation concerns the stencil of these operators which have to be limited, if centered, to three points by direction.

Application to the Navier-Stokes equations
The SMP method has been applied to the NS equations in [29,31]. For immersed boundary problems, the SMP and the AIIB methods give equivalent results and the AIIB method can be used to immerse obstacles in fluid flows. Both methods can be used for the scalar and the NS equations. In the latter, the procedure is done componentwise for the velocity vector. However, the AIIB method, with L 1 1 interpolations only, cannot be applied to the NS equations on staggered grid (no tests have been performed for a collocated approach). An illustration is given Fig. 4. With such interpolations, two auxiliary unknowns u * I and u * ′ I , I ∈ N 1 , can coexist at the same location x I . Hence, u * I is the natural neighbor of u J and u * ′ I is the natural neighbor of u K . So a problem occurs for the discretization of the inertial term since a node of a given velocity component has to use an auxiliary unknown of an other velocity component. In this case, neither u * I nor u * ′ I are natural neighbors for v l , a velocity unknown in the y direction. No matter which unknown is used, or an average of the two collocated unknowns, the simulation is unstable outside of the Stokes regime. A particular attention has also to be given to the velocity-pressure coupling. If a fractional step method is used, the prediction step is modified by any fictitious domain method to impose an immersed boundary condition for the velocity. Thus, the projection step has to be modified according to the prediction step to remain consistent with the overall problem. Details about this point can be found in [31] where a consistent pressure correction is proposed in the framework of the penalty methods.

The AIIB for immersed interface problems with jump conditions
With the symmetric method described in (2.3.3), the problem can be solved on both sides of the interface when explicit Dirichlet BC are imposed. For many problems, the solution is not a priori known on the interface and some jump transmission conditions on the interface Σ are required. Let us now consider the problem: in Ω + Interface conditions on Σ where the interface conditions are: (a · ∇u) · n Σ = ψ on Σ (23) The notation Σ denotes the jump of a quantity over the interface Σ. In the symmetric version of the AIIB method, a given intersection point x l , l ∈ I, is associated with two auxiliary unknowns on both sides of the interface. Hence, the interface constraints (22) and (23) of (P i ) can be imposed at each intersection point x l by using the two auxiliary unknowns. For example, the I nth row of the matrix A ′ with u * I , I ∈ N 0 , can be used to impose the constraint (22) and the J nth line of the matrix with u * J , J ∈ N 1 , is then used to impose the constraint (23).

The solution constraint
The symmetrized AIIB methods for Dirichlet BC reads : when L 1 1 interpolations are used. With u Σ = u + Σ − u − Σ = φ, we obtain : which is the first constraint to be imposed.

The flux constraint
Following the same idea and using P 2 1 interpolations, for the case presented in Fig. 3 left. Using (23), we obtain: which is the second constraint to be imposed. With such an interpolation, the solution gradient is constant over the whole stencil. As demonstrated later, the second-order accuracy can be reached on Cartesian grids when ψ = 0. Three auxiliary unknowns are thus involved in the discretizations (25) and (27). The auxiliary unknown u * K is also involved in the discretization of (22) and (23) at another intersection point on Σ h . Hence, the whole system A ′ u ′ = b ′ is closed.
On can notice that in [4], the same kind of augmented system is considered. Contrary to the AIIB method, no auxiliary nodes are used andthe spatial discretization at the grid points at the vicinity of the interface has to be modified "by-hand".

Algebraic reduction
Since we need more than one auxiliary unknown to discretize each constraint, the matrix C 2 is not diagonal and a solver has to be used to compute C −1 2 . For the matched interface and boundary (MIB) method, Zhou et al. [42] use a different discretization of the interface conditions which allows an easy algebraic reduction which is directly performed raw by raw.
The algebraic reduction for the immersed interface problems has not been yet implemented. However, the standard discretization of the AIIB method requires a more compact stencil than for the MIB method, and the additional computational time generated by the auxiliary nodes is small. Hence, the lack of algebraic reduction does not seem to a problem.

Numerical results for scalar problems
Elliptic equations are discretized using the standard second-order centered Laplacian. For all problems, similar results have been obtained with a PARDISO direct solver [33], and an iterative BiCGSTAB solver [10], preconditioned under a ILUK method [28]. Unless otherwise mentioned, a numerical domain [−1; 1] × [−1; 1] is used for every simulation. Two discrete errors are used. The discrete relative L 2 error in a domain Ω is defined as: withũ the analytical solution.
The discrete L ∞ error is defined as: Only Ω 0 is taken into account for the immersed boundary problems.

Problem 1
The homogenous 2D Laplace equation is solved. The interface Σ is a centered circle of radius R 1 = 0.5 m with a Dirichlet condition of U 1 = 10 • C. An analytical solution which accounts for the presence of a second circle with a radius R 2 = 2 m and U 2 = 0 • C is imposed on the boundary conditions. The analytical solution is: Accuracy tests are performed with L 1 1 , P 2 1 and Q 2 1 interpolations. Fig. 5 shows the solution and the error map for a 32 × 32 mesh with P 2 1 interpolations. The same results are always obtained with and without algebraic reduction. Fig. 6 shows the convergence of the error for the L 2 and L ∞ norms. For all Figure 5: Solution and error map for problem 1 on a 32 × 32 grid interpolations, the convergence slopes are approximatively 2 for the relative L 2 error. For the L ∞ error, the slopes are about 1.8. The P 2 1 interpolation is the more accurate, followed by the L 1 1 interpolation although it uses more auxiliary points (but a smaller stencil). However, the differences of accuracy between the different interpolations remain small. The performances of the ILUK-BiCG-Stab solver are now benchmarked for the three interpolations with and without algebraic reduction and for the SMP method. Tab. 1 shows the computational times of the matrix inversions (average time in seconds for 25 matrix inversions) and Tab. 2 shows the time ratio between the standard and the reduced matrix. Except for the Q 1 1 interpolation on the 1024×1024 mesh, the differences between the two methods seem to decrease with the size of the matrix. In fact, as interfaces are d − 1 manifolds, the number of intersection points does not increase as fast as the Eulerian points. Hence, the ratio between the size of a reduced and a complete matrix tends to 1. The computational time for the SMP method is quite similar to the one obtained with the AIIB method and algebraic reduction. Figures 7, 8, 9, 10 shows the convergence of the ILUK-BiCG-Stab solver for the seven configurations. The type of interpolation does not significantly impact on solver performances. As expected, the number of solver iterations have to be increased to reach a given residual when the number of computational nodes is also increased.

Problem 2
The 3D equation ∆T = 6 is solved. The solution is T (r) = r 2 . The solution is imposed on an immersed centered sphere of radius 0.2 m. As expected, the second-order code gives the exact solution to almost computer-error accuracy without this inner boundary. Results of the numerical accuracy test with the spherical inner boundary are presented in Fig. 11. The average slope for the L 2 norm is 2.33 and increases for the denser meshes. Even if the method has Figure 11: Curves of errors for problem 2 the same convergence order as the initial discretizeation, computer error accuracy can not be expected, at least because the immersed boundary Σ h is an approximation of the initial sphere Σ.

Problem 3
The 3D equation ∆T = 12r 2 is solved. The solution is T (r) = x 4 + y 4 + z 4 . The results are presented in Fig. 12. For the L ∞ norm, the second order is regularly obtained. For the L 2 norm, the second order is not obtained for the coarsest meshes as the code has not reached its asymptotical convergence domain. As can be noticed by comparing results with and without the AIIB method, this last one does not spoil the convergence order of the code, and the presence of the immersed interface with an analytical solution imposed on Σ h improves the accuracy. For both cases, the numerical solution tends to a second order in space.

Problem 4
The 2D equation ∆T = 4 is solved. The analytical solution is imposed on the boundaries of the domain and a Neuman BC is imposed on a centered circle of radius R = 0.5 m. As can be seen in Fig. 13, the global convergence has an average slope of 1.10 and can be explained by the simplicity of the discretization of the Neumann BC.

Problem 5
The 2D problem P ii with f = −4 and a = 1 is solved. As the equation remains the same in both domains, this problem can be solved without immersed interface method. The analytical solution is u = r 2 . As can be expected with our second order code, computer error is reached for all meshes with or without AIIB method. The difference with problem 2, where the solution is a second-order polynomial too, is that the solution is not explicitly imposed at a given location. In the present case, the interface condition is still correct anywhere in the domain so the approximation of the interface position does not generate errors. Fig. 14 shows that the same result is obtained with a solution jump on a circular interface such as u = r 2 for r > 0.5 and u = r 2 + 1 otherwise.

Problem 6
The same problem as in 3.2.1 is now considered with a discontinuous coefficient a such as a = 10 in Ω 0 and a = 1 in Ω 1 , involving the following analytical solution: Accuracy tests are first performed with the interface almost passing by some grid points (called odd mesh). The interface does not strictly lies on these points, as the shape is shifted by an ϵ = 10 −10 . This configuration is difficult as the interpolations degenerates. Accuracy tests are then performed with a box of length 1.0001 (called even mesh). In this configuration, the interface never passes by a grid point. The results of the numerical accuracy test are presented in Fig. 16. For the odd series of test, the slope is 1.86 for the L 2 and L ∞ errors. For the even series, where no geometrical singularity is present, the slope for both errors is 2.04. Figures 17 shows the solution and the L 2 relative error for a 32 × 32 mesh. As the analytical solution is imposed on the numerical boundary, the error is principally located in the interior subdomain.

Problem 7
The homogenous 2D Laplace equation is considered with the following analytical solution: where Ω 0 and Ω 1 are delimited by Σ a centered circle of radius 0.5 m. Fig. 18 shows that the convergence for both L 2 and L ∞ error are of first order only. The Fig. 19 shows the numerical solution (which is not so different from the analytical solution) and the error map for a 32 × 32 mesh. Hence, the drawback of the compacity of the discretization is a first-order convergence in space only for cases with flux jump.

Convergence
We measure the sensibility of the method with the accuracy of the discretization of the immersed interface. Problem 1 is solved on 32 × 32 and 128 × 128 meshes. Fig. 20 shows the accuracy of the solution with respect to the number of points used to discretized the interface which is here a circle. The reference solutions (Fig. 6) have been computed with an analytical circle. As can be seen, a second order in space is globally obtained. The numerical solutions of reference for the 32 × 32 and 128 × 128 meshes are different but the sensitivity of the error to the number of points in the Lagrangian mesh is almost the same. Figure 18: Convergence of the L 2 relative error and the L ∞ error for problem 7 Figure 19: The solution and the L 2 relative error for problem 7 with a 32 × 32 mesh

The Stanford bunny
This last case demonstrates how a second-order method enhances the representation of the boundary condition compared to a first-order method. The homogenous Laplace problem with a Dirichlet BC T Σ = 10 • C is solved on a 60 × 60 × 50 mesh bounding an obstacle of complex shape (the Stanford bunny). The extension of the solution in Ω 1 is used for the post treatment. Thus, all u J , J ∈ N 1 , are replaced by u * J . Then, the iso-surface T = T Σ gives an idea of the approximation of the boundary condition. Fig. 21 shows the iso-surface for a first order method. As can be seen, the shape of the obstacle endures a rasterization effect as the solution is imposed in the entire control volumes. Fig. 22 shows the iso-surface for the second order AIIB method. The improvement brought by this approach is straightforward. Fig. 23 shows a slice of the solution passing through the bunny. As can be seen, overshoots are present inside the shape which corresponds to the auxiliary values allowing the correct solution at the Lagrangian interface points to be obtained.

Some remarks about the solvers
The kind of interpolation function used and the position of the interface have an impact on the final discretization matrix C ′ , especially on its conditionning. Let us consider an intersection x l of Σ h between Hence, α 1−α tends to 0 when x l tends to x J . As the matrix loses its diagonal dominance, solver problems can be encountered. Tseng et al. [37] proposed changing the interpolation by using a new node which is the image of x I through the interface. In [8,7], authors pointed out this problem and suggest to slightly move the interface to a neighboring point (in our case x J ) if x I is too close to Σ h .
In this case, for the Dirichlet BC, an unknown u * J is created, and the equation at x J is simply u J = u l . For the Neumann BC, the standard interpolation is written at x J with u * J and its neighbor unknowns in Ω 0 .
For the transmission conditions (22)- (23), if ϕ = 0 and ψ = 0, no auxiliary unknown is created and the standard finite-volume centered discretization is used. However, for this case, or for ϕ ̸ = 0 and ψ ̸ = 0, our implementation using ILUK preconditionner or a PARDISO direct solver does not necessarily require such methods, even if α 1−α ≈ 10 −10 . Correct solutions are always obtained.

Discussion and conclusion
A new immersed interface method, the algebraic immersed interface and boundary (AIIB) method, which uses algebraic manipulations and compact stencil discretizations, has been presented. This method is able to treat elliptic equations with discontinuous coefficients and solution jumps over complex interfaces. A second order in space is reached for several configurations with minor modifications of the original code (especially if the reduction of the linear system is not considered).
The aim was to design a method as simple as possible. Contrary to some extensions of the IIM and MIB methods, no particular attention has been paid to treat interfaces with critical curvature or singularities. However, the compact stencil of AIIB method allows it to treat interfaces with quite high curvatures without modification and satisfies the goal of simplicity.
For the immersed boundary problems with a Dirichlet BC, the method has shown a second order of convergence in space for various kinds of interpolations. An algebraic reduction has been applied to accelerate the convergence of the solver. For the Neumann BC, a first order has been obtained.
For the immersed interfaces, a second order of convergence in space is obtained when the jump of the normal flux is zero, even if the equation has discontinuous coefficients.
Future work could be devoted to increase the accuracy of the method when the jump of the normal flux is not zero. It is the main drawback of the method compared to the IIM, MIB method or [4]. A challenge will be to keep a compact stencil. A study of the matrix conditioning would be important too as a strong solver is required for the present method. Another interesting point would be to couple the method with alternative interface conditions such as the Jump Embedded Boundary Conditions proposed in [1] which are : (a · ∇u) · n Σ = β u − g on Σ (36) where u Σ = (u + − u − )/2 denotes the arithmetic mean of the traces of a quantity on both sides of the interface, α, β, h and g are scalar values which can be chosen in order to obtain Dirichlet, Neumann, Fourier of transmission conditions on the interface. At last, a long-term goal is to extend the method to the Navier-Stokes equations with immersed interfaces. To perform fluid/structure coupling, the implicit tensorial penalty method [27,35] can be considered. With this approach, the solid medium is treated as a fluid with a high viscosity. At the fluid/solid interface, average physical quantities are imposed. Such strategy is generally less accurate that methods using polynomial interpolations so a coupling with the AIIB method is desirable.

Appendix : Definition of the interpolation
We explain the calculation of the interpolation coefficients for 2D problems. Let us consider x I , x J , x K and x L which define a dual control cell V ′ I . A p ∈ Q 2 1 interpolation over V ′ I is such that p(x i ) = u i for i = I, J, K, L, and p(x, y) = a 0 + a 1 x + a 2 y + a 3 xy. The following coordinates matrix can be defined : If a is the vector of the interpolation coefficient, p(x, y) = aQ and a = Q −1 p. As each term a i of a is a linear combination of u i , one can write p(x, y) = ∑ i=I,J,K,L α i u i with (x, y) the coordinates of the Lagrangian intersection point x l , l ∈ I. Practically, the four Eulerian points are the four corners of an unit square and x l is easily projected in this new frame by ensuring Fig. 24 illustrates the projection. Then, a unique trivial Q −1 matrix has to be found for each kind of interpolation. The coefficients for the p ∈ Q 2 1 interpolation are the following: In [37], the authors uses three Eulerian points and an interface point to construct the interpolation Figure 24: The projection of the initial square defined by the Eulerian mesh to an unit square providing a more complex linear system which has to be solved for each intersection point.