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

Arthur Sarthou^1*, Stéphane Vincent², Philippe Angot³, Jean-Paul Caltagirone⁴

¹Buf Compagnie, Paris, France.
²Laboratoire MSME, Université Gustave Eiffel, Champs-sur-Marne, France.
³Institut de Mathématiques de Marseille, Centrale Marseille, Marseille, France.
⁴Institut Polytechnique de Bordeaux, Université de Bordeaux, Talence, France.

DOI: 10.4236/ojfd.2020.103015                

Abstract

A new simple 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 with a good accuracy on a compact stencil. Auxiliary unknowns are created at existing grid locations to increase the degrees of freedom of the initial problem. These auxiliary unknowns allow imposing 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.

Keywords

Fictitious Domain, Immersed Interface Method, Immersed Boundary Method, Penalty Methods, Finite Volumes, Elliptic Equations, Jump Embedded Conditions

1. 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 first deals 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 a domain of interest ${\Omega }_0$, typically a fluid domain, which is embedded inside a computational domain of simple shape $\Omega \subset {ℝ}^d$, d being the spatial dimension of the problem. The auxiliary domain ${\Omega }_1$, typically a solid particle or an obstacle, is such that as shown in Figure 1: $\Omega ={\Omega }_0\cup \Sigma \cup {\Omega }_1$ where $\Sigma$ is an immersed interface, and $n$ is the unit outward normal vector to ${\Omega }_0$ on $\Sigma$. Let us consider the following model scalar immersed boundary problem in ${\Omega }_0$ with a Dirichlet boundary condition (BC) on the interface $\Sigma$ :

${\mathcal{P}}_b\{\begin{array}{l}-\nabla \cdot \left(a\nabla u\right)=f\text{in}{\Omega }_0\\ u_{\mathrm{<mo>\; |\; </mo>}\Sigma }=u_D\text{on}\Sigma \end{array}$

A boundary condition is also required on the other part of the boundary $\partial {\Omega }_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. [1]. Lagrange multipliers are introduced into the weak formulation of the initial elliptic equation to ensure the immersed boundary condition.

Figure 1. Definition of the subdomains and the interface.

Cartesian grid [2] [3] and Cut-cell [4] 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 [5] [6]. 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 [7]. The idea here is to impose a no-slip condition directly on the boundary using a mirrored flow over the boundary. In [8] [9], 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 [10], 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 [11], this kind of approach seems to be more accurate than [8] [9].

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) [12] requires the addition of a penalty term $\frac{\chi }{\epsilon }\left(u-u_D\right)$ in the conservation equations, such that:

$\{\begin{array}{l}-\nabla \cdot \left(a\nabla u\right)+\frac{\chi }{\epsilon }\left(u-u_D\right)=f\text{in}\Omega \\ \text{with}{\chi }_{|{\Omega }_0}=0,{\chi }_{|{\Omega }_1}=1,\text{for}0<\epsilon \ll 1\end{array}$ (1)

where $\epsilon$ denotes the penalty parameter which tends to 0. Hence, in ${\Omega }_1$ the original equation becomes negligible and $u=u_D$ is imposed. In [13] [14] authors add a Darcy term $\frac{\mu }{K}u$ to the Navier-Stokes (NS) equations where $\mu$ is the dynamic viscosity and K the permeability. In the fluid medium, $K\to \infty$ so the Darcy term is then negligible and the original set of NS equations is retrieved. In the solid medium, $K\to 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 [15]. In [16] [17], Sarthou et al. have discretized the volume penalty term with a second order using implicit interpolations as in [10]. This method is called the sub-mesh penalty (SMP) method and has been applied to both elliptic and NS equations.

Applied to problem ${\mathcal{P}}_b$, the ghost cell immersed boundary method [10] and SMP method [16] used the cells in ${\Omega }_1$ in the vicinity of the interface to enhance the accuracy of the solution in ${\Omega }_0$.

Another approach which considers the extension of the solution is considered in [18] [19] 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:

$\left({\mathcal{P}}_i\right)\{\begin{array}{ll}-\nabla \cdot \left(a\nabla u\right)=f\hfill & \text{in}\Omega \hfill \\ {〚u〛}_{\Sigma }=\phi \hfill & \text{on}\Sigma \hfill \\ {〚\left(a\cdot \nabla u\right)\cdot n〛}_{\Sigma }=\psi \hfill & \text{on}\Sigma \hfill \end{array}$

A first class of method is the immersed interface methods (IIM) initially introduced by LeVeque and Li [20] and widely described in [21]. This group of methods uses Taylor series expansion of the solution at discretization points in the vicinity of $\Sigma$ 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 [22] or Navier-Stokes equations [23]. In [24], Li uses an augmented variables 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 [25] and Navier-Stokes [26].

The Ghost Fluid Method, originally developed by Fedkiw et al. [27] [28], 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 [29] [30] [31] 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 [24], 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 [24] to discretize additional equations. In [32], 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. Its main contribution is its simplicity that eases its implementation, and its generality that allows the method to be applied to various boundary and interface conditions. Contrary to ghost nodes, auxiliary unknowns are present in the linear system and the discretization of the interface conditions is independant from the discretization of the initial equation.

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.

2. 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.

2.1. Definitions and Notations

Our objective is to numerically impose the adequate conditions on the interface $\Sigma$ for problems similar to ${\mathcal{P}}_b$ and ${\mathcal{P}}_i$. These conditions will be discretized in space on an Eulerian structured mesh covering $\Omega$. As the discretization of the interface or boundary conditions requires interpolations, the following interpolations in 2D: ${\mathbb{P}}_1^2\left(x\mathrm{,}y\right)=p_1+p_{2}x+p_{3}y$ and ${\mathbb{Q}}_1^2\left(x\mathrm{,}y\right)=p_1+p_{2}x+p_{3}y+p_{4}xy$ are used. In 3D, we use ${\mathbb{P}}_1^3\left(x\mathrm{,}y\mathrm{,}z\right)=p_1+p_{2}x+p_{3}y+p_{4}z$ and ${\mathbb{Q}}_1^3\left(x\mathrm{,}y\mathrm{,}z\right)=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, ${\mathbb{L}}_1^1\left(x\right)=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.

The computational domain $\Omega$ is approximated with a curvilinear mesh $T_h$ composed of $N\times M$ ( $\times L$ in 3D) cell-centered finite volumes ( ${\mathcal{V}}_I$ ) for $I\in \mathcal{E}$, $\mathcal{E}$ being the set of indexes of the Eulerian orthogonal curvilinear structured mesh. Let $x_I$ be the vector coordinates of the center of each volume ${\mathcal{V}}_I$. In 2D, the horizontal and vertical mesh steps are respectively $h_x$ and $h_y$. This grid is used to discretize the conservation equations. A dual grid is introduced for the management of the AIIB method. The grid lines of this dual cell-vertex mesh are defined by the network of the cell centers $x_I$. The volumes of the dual mesh are denoted by ( ${{\mathcal{V}}^{\prime }}_I$ ). The Eulerian unknowns are noted $u_I$ which are the approximated values of $u\left(x_I\right)$, i.e. the solution at the cell centers $x_I$.

The discrete interface ${\Sigma }_h$, hereafter called the Lagrangian mesh, is given by a discretization of the original interface $\Sigma$. It is described by a piecewise linear approximation of $\Sigma$ : ${\Sigma }_h=\left\{{\sigma }_l\subset {ℝ}^{d-1}\mathrm{,}l\in {\mathcal{L}}_f\right\}$, ${\mathcal{L}}_f$ being the set of indexes of the Lagrangian mesh and K being the cardinal of ${\mathcal{L}}_f$. Typically, ${\sigma }_l$ are segments in 2D and triangles in 3D. The vertices of each face ${\sigma }_l$ are denoted by $x_{l\mathrm{,}i}$ for $i=\mathrm{1,}\cdots \mathrm{,}d$ and the set of all vertices is $\left\{x_l\mathrm{,}l\in {\mathcal{L}}_v\right\}$, ${\mathcal{L}}_v$ being the set of indexes of the vertices. The intersection points between the grid lines of the Eulerian dual mesh and the faces ${\sigma }_l$ of the Lagrangian mesh are denoted by $\left\{x_i\mathrm{,}i\in \mathcal{I}\right\}$ (see Figure 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 functions in order to implicitly locate ${\Sigma }_h$ :

• The Heaviside function $\chi$, defined as:

$\chi \left(x\right)=\{\begin{array}{l}1\text{if}x\in {\Omega }_1\hfill \\ 0\text{otherwise}\hfill \end{array}$ (2)

This function is built with a point in solid method presented below. The function $\chi$ will be used to perform fictitious domain algorithms and to build a level-set function.

• The level-set function $\varphi$, with:

$\varphi \left(x\right)=\{\begin{array}{l}-{\text{dist}}_{\Sigma }\left(x\right)\text{if}x\in {\Omega }_1\hfill \\ {\text{dist}}_{\Sigma }\left(x\right)\text{otherwise}\hfill \end{array}$ (3)

and ${\text{dist}}_{\Sigma }\left(p\right)=in{f}_{x\in \Sigma }\Vert x-p\Vert$. The unsigned distance is computed geometrically. The sign is directly obtained with the discrete Heaviside function $\chi$.

• The colour phase functions C, which is the ratio of a given phase in a control volume. We denote $C\left(x_I\right)$ the phase ratio in the control volume centered in $x_I$. This function is approximated from the $\varphi$ function by using the formula proposed by Sussman and Fatemi [33]:

Figure 2. Definition of the discretization kernels for the AIIB method.

$C\left(x\right)\approx \{\begin{array}{l}1\text{if}\varphi \left(x\right)>h\hfill \\ 0\text{if}\varphi \left(x\right)<-h\hfill \\ \frac{1}{2}\left(1+\frac{\varphi \left(x\right)}{h}+\frac{1}{\pi }\sin \left(\pi \varphi \left(x\right)/h\right)\right)\text{otherwise}\hfill \end{array}$ (4)

New sets of Eulerian points $x_I$ are defined near the interface so that each one has a neighbor $x_J$ verifying ${\chi }_J\ne {\chi }_I$ (with ${\chi }_I=\chi \left(x_I\right)$ and ${\chi }_J=\chi \left(x_J\right)$ ), i.e. the segment $\left[x_I\mathrm{<mo>\; ;\; </mo>}{x}_J\right]$ is cut by ${\Sigma }_h$. These Eulerian “interface” points are also sorted according to their location inside ${\Omega }_0$ or ${\Omega }_1$. Two sets $\left\{x_I\mathrm{,}I\in {\mathcal{N}}_0\right\}$ and $\left\{x_I\mathrm{,}I\in {\mathcal{N}}_1\right\}$ are thus obtained, where ${\mathcal{N}}_0=\left\{I\mathrm{,}{x}_I\in {\Omega }_0\mathrm{,}{x}_I\in Neighb\left(x_J\right)\mathrm{,}{\chi }_I\ne {\chi }_J\mathrm{,}{x}_J\in {\Omega }_1\right\}$ and ${\mathcal{N}}_1=\left\{I\mathrm{,}{x}_I\in {\Omega }_1\mathrm{,}{x}_J\in Neighb\left(x_I\right)\mathrm{,}{\chi }_I\ne {\chi }_J\mathrm{,}{x}_J\in {\Omega }_0\right\}$.

For each $x_I$, $I\in {\mathcal{N}}_0$ or $I\in {\mathcal{N}}_1$, we associate two unknowns: the physical one denoted as $u_I$ and the auxiliary one $u_I^{\mathrm{<mo>\; *\; </mo>}}$.

2.2. 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 $\chi$, one has to determine which Eulerian points are inside the domain ${\Omega }_1$ defined by a Lagrangian surface. Such a surface must be closed and not self-intersecting. In [11] [14], the authors used a global methodology partly based on [34] where $\chi$ 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 [35] [36]. Instead of casting rays in any direction, rays are all aligned on a principal direction of the grid such as each ray passes through a line of points $x_I$. Hence, only one ray is cast for each line of points allowing the complexity of the algorithm to be reduced by one order.

The Lagrangian points $x_l$ of ${\Sigma }_h$, $l\in \mathcal{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\in {\Omega }_0$ and $x_J\in {\Omega }_1$. We denote by $d_I=d\left(x_I\mathrm{,}{\Sigma }_h\right)$ and $d_J=d\left(x_J\mathrm{,}{\Sigma }_h\right)$ the unsigned distances between Eulerian points and the interface ${\Sigma }_h$. Then, $x_l=\left(x_I{d}_J+x_J{d}_I\right)/\left(d_I+d_J\right)$.

Algorithmic problems can be encountered if the Lagrangian mesh is too complex compared to the Eulerian mesh. In some cases, more than one segment of the Lagrangian mesh can intersect a same grid segment $\left[x_I\mathrm{<mo>\; ;\; </mo>}{x}_J\right]$, and more than one intersecting points can then 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.

2.3. The AIIB for Dirichlet and Neumann Conditions

2.3.1. 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 [11] [16] 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 [18] [19] and the immersed interface methods [20]), the AIIB method is thus simple to implement.

Let $\mathcal{P}$ be a model problem discretized in the whole domain $\Omega$ 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 original problem $Au=b$ becomes $A^{\prime }{u}^{\prime }=b^{\prime }$, with $A^{\prime }$ a square matrix of order $m+n$, with n the number of auxiliary constraints related to the interface conditions. The solution $u^{\prime }$ is decomposed such as $u^{\prime }={\left(u\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\right)}^{\text{T}}$ and the source term as $b^{\prime }={\left(b\mathrm{,}{b}^{\mathrm{<mo>\; *\; </mo>}}\right)}^{\text{T}}$. The interface constraints are discretized with a $\left(n\mathrm{,}m+n\right)$ block matrix C and the source term $b^{\mathrm{<mo>\; *\; </mo>}}$. The discretizations of the constraints for various interface and boundary conditions are given in Sections 2.3.2 and 2.4.

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^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}I\in {\mathcal{N}}_1$ (resp. $u_I^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}I\in {\mathcal{N}}_0$ ) as the extension of the solution in ${\Omega }_0$ (resp. ${\Omega }_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 ${\alpha }_{I\mathrm{,}J}$ be a coefficient of A at row I, column J and ${{\alpha }^{\prime }}_{I\mathrm{,}J}$ the new coefficient in $A^{\prime }$. If $I\in {\mathcal{N}}_0$ and $J\in {\mathcal{N}}_1$, ${{\alpha }^{\prime }}_{I\mathrm{,}J}=0$.

Hence, there is no more direct link between the initial unknowns in ${\Omega }_0$ and ${\Omega }_1$. However, the initial value of ${\alpha }_{I\mathrm{,}J}$ is reused to express the new relation between the original unknown $u_I$ and its new neighbor $u_J^{\mathrm{<mo>\; *\; </mo>}}$. Consequently, ${{\alpha }^{\prime }}_{I,J^{*}}={\alpha }_{I,J}$, where $J^{\mathrm{<mo>\; *\; </mo>}}$ is the index corresponding to $u_J^{\mathrm{<mo>\; *\; </mo>}}$.

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 $\left(m\mathrm{,}m+n\right)$ and $I_2$ of size $\left(n\mathrm{,}m+n\right)$ :

$I_1=\left(\begin{array}{ccccccc}1& 0& \cdots & 0& \cdots & \cdots & 0\\ 0& \ddots & & ⋮& \ddots & & ⋮\\ ⋮& & \ddots & 0& & \ddots & ⋮\\ 0& \dots & 0& 1& 0& \cdots & 0\end{array}\right)$ (5)

$I_2=\left(\begin{array}{ccccccc}0& \cdots & 0& 1& 0& \cdots & 0\\ 0& \ddots & & 0& \ddots & & ⋮\\ ⋮& & \ddots & ⋮& & \ddots & 0\\ 0& \dots & 0& 0& \cdots & 0& 1\end{array}\right)$ (6)

The matrices $A_0$ and $A_1$ are defined such as $A_0+A_1=A$, $A_0\left(I,J\right)=A\left(I,J\right)$ if $I\in {\mathcal{N}}_0$, else $A_0\left(I,J\right)=0$. Similarly $A_1\left(I,J\right)=A\left(I,J\right)$ if $I\in {\mathcal{N}}_1$ else $A_1\left(I,J\right)=0$. 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\in {\mathcal{N}}_0$, $J\in {\mathcal{N}}_1$ and $P_1$ to switch row I with row J if $I\in {\mathcal{N}}_1$, $J\in {\mathcal{N}}_0$. Hence, the new problem matrix is now defined by:

$A^{\prime }=I_1^{\text{T}}\left(P_0\left(A_0{I}_1\right)+P_1\left(A_1{I}_1\right)\right)+I_2^{\text{T}}C$ (7)

The new problem is $A^{\prime }{u}^{\prime }=b^{\prime }$ with $A^{\prime }$ written with 4 blocks of various sizes: $\tilde{A}\left(m\mathrm{,}m\right)$, $B\left(m\mathrm{,}n\right)$, $C_1\left(n\mathrm{,}m\right)$, $C_2\left(n\mathrm{,}n\right)$. The matrix $\tilde{A}$ is thus the modification of the initial matrix A by setting to zero the coefficient ${\alpha }_{I\mathrm{,}J}$ if $\chi \left(x_I\right)\ne \chi \left(x_J\right)$, and $C_1$ and $C_2$ are the two sub-matrices of the matrix C. The problem can be written as:

$\left(\begin{array}{cc}\tilde{A}& B\\ C_1& C_2\end{array}\right)\left(\begin{array}{c}u\\ u^{\mathrm{<mo>\; *\; </mo>}}\end{array}\right)=\left(\begin{array}{c}b\\ b^{\mathrm{<mo>\; *\; </mo>}}\end{array}\right)$ (8)

The entire problem can then be solved to obtain $u^{\prime }={\left(u\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\right)}^{\text{T}}$. However, $u^{\mathrm{<mo>\; *\; </mo>}}$ 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:

$\left(\tilde{A}-B{C}_2^{-1}{C}_1\right)u=b-B{C}_2^{-1}{b}^{\mathrm{<mo>\; *\; </mo>}}$ (9)

The opportunity of such a reduction will be discussed later.

2.3.2. AIIB Algorithm for a Scalar Equation with Dirichlet Boundary Conditions

For the sake of clarity, let us first describe in 2D the AIIB method for the model scalar problem ${\mathcal{P}}_b$ with a Dirichlet boundary condition on the interface $\Sigma$. For this version of the AIIB algorithm, ${\Omega }_0$ is the domain of interest and auxiliary unknowns are created in ${\Omega }_1$ only. Let us consider a point $x_I\mathrm{,}I\in {\mathcal{N}}_1$. At location $x_I$, two unknowns coexist: a physical one $u_I$ and an auxiliary one $u_I^{\mathrm{<mo>\; *\; </mo>}}$. We first describe the case when $x_I$ has only one neighbor $x_J$ in ${\Omega }_0$. The Lagrangian point $x_l$ is the intersection between $\left[x_I\mathrm{<mo>\; ;\; </mo>}{x}_J\right]$ and ${\Sigma }_h$ (Figure 3 right). Then, the solution $u_l=u_D\left(x_l\right)$ at the interface is approximated by the ${\mathbb{P}}_1^1$ interpolation between the Eulerian unknowns $u_I^{\mathrm{<mo>\; *\; </mo>}}$ and $u_J$ :

$u_l={\alpha }_I{u}_I^{*}+{\alpha }_J{u}_J\text{with}0<{\alpha }_I,{\alpha }_J<1\text{and}{\alpha }_I+{\alpha }_J=1$ (10)

To our knowledge, it has been verified that for schemes close to the one presented here (including ghost-like schemes [19] and IBM-Direct Forcing schemes [10]) only a linear interpolation is required to reach a second order of accuracy as a Dirichlet BC is imposed here. It has been theoretically demonstrated in [37] that under the assumption of pointwise stability of the scheme the discretization at the boundary can be two orders less accurate than in the rest of the domain without reducing the global convergence order. It is important to notice this does not apply for immersed interface conditions [32].

If now $x_I$ has a second neighbor $x_K$ in ${\Omega }_0$, the intersection $x_m$ between $\left[x_I\mathrm{<mo>\; ;\; </mo>}{x}_K\right]$ and ${\Sigma }_h$ is considered with $u_m=u_D\left(x_m\right)$. We choose $x_p$, a new point of ${\Sigma }_h$ between $x_l$ and $x_m$ (see Figure 3 left). The solution $u_p=u_D\left(x_p\right)$ is then imposed using a ${\mathbb{P}}_1^2$ -interpolation of the values $u_I^{\mathrm{<mo>\; *\; </mo>}}$, $u_J$ and $u_K$ :

$u_p={\alpha }_I{u}_I^{*}+{\alpha }_J{u}_J+{\alpha }_K{u}_K,0<{\alpha }_I,{\alpha }_J,{\alpha }_K<1,{\alpha }_I+{\alpha }_J+{\alpha }_K=1$ (11)

A ${\mathbb{Q}}_1^2$ 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 Figure 3 left). As a third

(a) (b)

Figure 3. Example of selection of points for Dirichlet (a) and Neumann (b) constraints.

choice, two independent ${\mathbb{L}}_1^2$ linear 1D interpolations can be used (one for each direction) for an almost equivalent result. It produces:

$\{\begin{array}{l}{u}_l={\alpha }_I{u}_I^{*}+{\alpha }_J{u}_J\text{with}0<{\alpha }_I,{\alpha }_J<1\text{and}{\alpha }_I+{\alpha }_J=1\\ u_m={{\alpha }^{\prime }}_I{u}_I^{*}{}^{\prime }+{\alpha }_K{u}_K\text{with}0<{{\alpha }^{\prime }}_I,{\alpha }_K<1\text{and}{{\alpha }^{\prime }}_I+{\alpha }_J=1\end{array}$ (12)

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=\left(u_l+u_m\right)/2$. A summation of the two constraints from (12) gives:

${\alpha }_I{u}_I^{*}+{\alpha }_J{u}_J+{\alpha }_{I^{\prime }}{u}_I^{*}+{\alpha }_K{u}_K=u_l+u_m$ (13)

which is equivalent to build a constraint imposing $u_p$ at $x_p$ with a ${\mathbb{P}}_1^2$ interpolation:

$\begin{array}{l}\frac{\left({\alpha }_I+{{\alpha }^{\prime }}_I\right)u_I^{*}+{\alpha }_J{u}_J+{\alpha }_K{u}_K}{2}=u_p,\\ \text{with}0<\frac{{\alpha }_I+{\alpha }_{I^{\prime }}}{2},\frac{{\alpha }_J}{2},\frac{{\alpha }_K}{2}<1,\frac{{\alpha }_I+{\alpha }_{I^{\prime }}}{2}+\frac{{\alpha }_J}{2}+\frac{{\alpha }_K}{2}=1\end{array}$ (14)

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 ${\sigma }_l$ of ${\Sigma }_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 ${\Sigma }_h$, especially for interfaces of strong curvature. However, the distance $d\left(x_p\mathrm{,}{\Sigma }_h\right)$ between $x_p$ and ${\Sigma }_h$ varies like $\mathcal{O}\left(h^2\right)$ 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 ${\Sigma }_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 ${\Omega }_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 ${\Omega }_0$ also decreases.

Concerning the C submatrix, it is filled with the coefficients $\alpha$. One line corresponds to the discretization of the boundary constraint related to one auxiliary node. Each value of the second member $b^{\prime }$ is the boundary value $u_l$ related on auxiliary point.

Hence, the present method is suitable to impose a Dirichlet boundary condition on $\Sigma$ for ${\Omega }_0$, when the solution in ${\Omega }_1$ is of no interest. The solution $u_I^{\mathrm{<mo>\; *\; </mo>}}$ for $I\in {\mathcal{N}}_1$ is an extrapolation of the solution in ${\Omega }_0$ in order to satisfy the boundary condition on $\Sigma$ and thus is non-physical. Hence, the solution at the nodes of ${\Omega }_1$ far from the interface does not impact on the solution in ${\Omega }_0$ as the unknowns $u_I$ for $I\in {\mathcal{N}}_1$ have no more algebraic relations with the nodes $u_I$ for $I\in {\mathcal{N}}_0$. Nevertheless, if nothing specific is performed, a solution with no interest is computed by default in ${\Omega }_1$. A prescribed solution inside ${\Omega }_1$ can be obtained easily with a volume penalty method such as VPM [13] where initial equations are penalized to impose a desired analytical solution. However, a desired solution in ${\Omega }_1$ can also be obtained by modifying the solution obtained by solver. Then, the computational cost of the approach can be reduced by switching the solving of $u_I\mathrm{,}{x}_I\in {\Omega }_1$ off, or by totally removing these nodes in the solving matrix.

2.3.3. 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:

$\{\begin{array}{ll}-\nabla \cdot \left(a\nabla u\right)=f\hfill & \text{in}\Omega \hfill \\ u_{|\Sigma }^{-}=u_D\hfill & \text{on}\Sigma \hfill \\ u_{|\Sigma }^{+}=u_G\hfill & \text{on}\Sigma \hfill \end{array}$ (15)

The problem (15) requires for each point $x_I$ a physical unknown $u_I$ as well as an auxiliary unknown $u_I^{\mathrm{<mo>\; *\; </mo>}}$ on both sides of the interface.

Practically, the algebraic modification part of the AIIB algorithm for a Dirichlet BC is applied a first time with ${\Omega }_0$ as domain of interest, and auxiliary unknowns are created near ${\Sigma }_h$ in ${\Omega }_1$. As a second step, the Heaviside function is modified as $\chi :=1-\chi$ and the algorithm is applied a second time. Now, ${\Omega }_1$ is the domain of interest and auxiliary unknowns are created near $\Sigma$ in ${\Omega }_0$.

Then, both problems are solved at the same time with the resolution of the linear system.

2.3.4. 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 $\Sigma$ :

$\{\begin{array}{ll}-\nabla \cdot \left(a\nabla u\right)=f\hfill & \text{in}{\Omega }_0\hfill \\ \left(a\nabla u\right)\cdot n=g_N\hfill & \text{on}\Sigma \hfill \end{array}$ (16)

The principle is about the same as for Dirichlet BC, and the same interpolations, once derived, can be used to approximate the quantity $\left(a\nabla u\right)\cdot n$. Hence, at any point $x_l\mathrm{,}l\in \mathcal{I}$ on ${\Sigma }_h$ we use

$\left(a\nabla u_l\right)\cdot n\approx \left(a\nabla p\left(x_l\right)\cdot n\right)\mathrm{.}$ (17)

For $p\in {\mathbb{Q}}_1^2$, we get $\nabla p\left(x\mathrm{,}y\right)\cdot n=\left(p_{3}y+p_2\right)n_x+\left(p_{3}x+p_1\right)n_y$ whereas for $p\in {\mathbb{P}}_1^2$, $\nabla p\left(x\mathrm{,}y\right)\cdot 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 Figure 3(a), with $p\in {\mathbb{P}}_1^2$, we have:

$\nabla p\left(x,y\right)\cdot n=\frac{u_I^{*}-u_J}{h_x}{n}_x+\frac{u_K-u_I^{*}}{h_y}{n}_y=u_I^{*}\left(\frac{n_x}{h_x}-\frac{n_y}{h_y}\right)+u_J\frac{n_x}{h_x}+u_K\frac{n_y}{h_y}$ (18)

The diagonal coefficient of the row related to $u_I^{\mathrm{<mo>\; *\; </mo>}}$ in $C_2$ is $\left(\frac{n_x}{h_x}-\frac{n_y}{h_y}\right)$. The case where $\frac{n_x}{h_x}\approx \frac{n_y}{h_y}$ leads to numerical instabilities. However, it can happens for this configuration case only if $n_x$ and $n_y$ have the same sign, which will not happen if the interface is sufficiently smooth. For a smooth enough interface, $n_x>0$ and $n_y<0$. A smoothed interface can be simply obtained using the normal vector of the segment $\left[x_l\mathrm{,}{x}_m\right]$ 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 stencil configurations.

When $x_I$ has only one neighbor $x_J$ in ${\Omega }_0$, the ${\mathbb{Q}}_1^2$ and ${\mathbb{P}}_1^2$ interpolations degenerate to ${\mathbb{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 ${\Omega }_0$ is caught to build ${\mathbb{P}}_1^2$ interpolations (see Figure 3 right). This point is a neighbor of $x_J$ and is taken as $\left[x_I\mathrm{,}{x}_J\right]\perp \left[x_J\mathrm{,}{x}_K\right]$. In 2D, two choices generally appear, and the point being so that the angle $\left(n\mathrm{,}{x}_K-x_J\right)$ is in $\left[-\pi /2\mathrm{<mo>\; ;\; </mo>}\pi /2\right]$ is taken.

2.3.5. 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 in such a way that only one auxiliary unknown is needed:

$u_I^{*}={\displaystyle \underset{J\in \mathcal{N}}{\sum }}{\alpha }_J{u}_J+u_S$ (19)

where $u_S$ is the source term. In this case, the matrix $C_2$ in (8) is diagonal and thus the Schur complement $\left(\tilde{A}-B{C}_2^{-1}{C}_1\right)$ is easy to calculate. Practically, when the algebraic reduction is made, $\tilde{A}$ is built directly by the suitable modification of A without considering the extended matrix $A^{\prime }$. The part $-B{C}_2^{-1}{C}_1$ is then added to $\tilde{A}$ whereas $-B{C}_2^{-1}{b}^{\mathrm{<mo>\; *\; </mo>}}$ is added to b. As it will be subsequently demonstrated, the algebraic reduction decreases the computational cost of the solver by 10% - 20%.

If only ${\mathbb{L}}_1^1$ interpolations are used with the algebraic elimination, the matrix obtained with this method is similar to the one obtained in [18] 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 ${\mathbb{P}}_1^2$ interpolations are used, the computed solution in ${\Omega }_0$ are very similar (within machine error if the solver perfectly converges) to those that can be obtained with the SMP [16] method (when the penalty parameter is sufficiently small) and the DF-IB method [10]. These methods discretize equivalent boundary constraints, but the way the constraint is applied to the initial discretization of the equation for the nodes $x_I\mathrm{,}I\in {\mathcal{N}}_1$ differs. 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\mathrm{,}I\in {\mathcal{N}}_0$ is modified, and without algebraic reduction, both auxiliary and physical unknowns coexist at $x_I\mathrm{,}I\in {\mathcal{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\in {\Omega }_0$ can only use in ${\Omega }_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.

2.4. 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 $\Sigma$ are required. Let us now consider the problem:

$\left({\mathcal{P}}_i\right)\{\begin{array}{ll}-\nabla \cdot \left(a\nabla u\right)=f\hfill & \text{in}\Omega \hfill \\ +\text{Interface}\text{conditions}\hfill & \text{on}\Sigma \hfill \end{array}$

where the interface conditions are:

$〚u_{\Sigma }〛=\phi \text{on}\Sigma$ (20)

${〚\left(a\nabla u\right)\cdot n〛}_{\Sigma }=\psi \text{on}\Sigma$ (21)

The notation ${〚〛}_{\Sigma }$ denotes the jump of a quantity over the interface $\Sigma$. In the symmetric version of the AIIB method, a given intersection point $x_l\mathrm{,}l\in \mathcal{I}$, is associated with two auxiliary unknowns, one on each side of the interface. Hence, the interface constraints (20) and (21) of $\left({\mathcal{P}}_i\right)$ can be imposed at each intersection point $x_l$ by using the two auxiliary unknowns. For example, the row I of the matrix $A^{\prime }$ with $u_I^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}I\in {\mathcal{N}}_0$, can be used to impose the constraint (20) and the line J of the matrix with $u_J^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}J\in {\mathcal{N}}_1$, is then used to impose the constraint (21).

2.4.1. The Solution Constraint

The symmetrized AIIB methods for Dirichlet BC reads:

$\{\begin{array}{l}{u}_{\Sigma }^{+}={\alpha }_1{u}_I+{\alpha }_2{u}_J^{*}\\ u_{\Sigma }^{-}={\alpha }_1{u}_I^{*}+{\alpha }_2{u}_J\end{array}$ (22)

when ${\mathbb{L}}_1^1$ interpolations are used. With $u_{\Sigma }=u_{\Sigma }^{+}-u_{\Sigma }^{-}=\phi$, we obtain:

${\alpha }_1{u}_I+{\alpha }_2{u}_J^{\mathrm{<mo>\; *\; </mo>}}-{\alpha }_1{u}_I^{\mathrm{<mo>\; *\; </mo>}}-{\alpha }_2{u}_J=\phi$ (23)

which is the first constraint to be imposed.

2.4.2. The Flux Constraint

Following the same idea and using ${\mathbb{P}}_1^2$ interpolations,

$\{\begin{array}{l}\left(a\cdot \nabla u_{\Sigma }^{+}\right)\cdot n=a^{+}\left(\frac{u_I-u_J^{*}}{h_x}{n}_x+\frac{u_K^{*}-u_I}{h_y}{n}_y\right)\\ \left(a\cdot \nabla u_{\Sigma }^{-}\right)\cdot n=a^{-}\left(\frac{u_I^{*}-u_J}{h_x}{n}_x+\frac{u_K-u_I^{*}}{h_y}{n}_y\right)\end{array}$ (24)

for the case presented in Figure 3 left. The normal is pointing from the + to the − sides. Using (21), we obtain:

$a^{+}\left(\frac{u_I-u_J^{*}}{h_x}{n}_x+\frac{u_K^{*}-u_I}{h_y}{n}_y\right)-a^{-}\left(\frac{u_I^{*}-u_J}{h_x}{n}_x-\frac{u_K-u_I^{*}}{h_y}{n}_y\right)=\psi$ (25)

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 $\psi =0$.

Three auxiliary unknowns are thus involved in the discretizations (23) and (25). The auxiliary unknown $u_K^{\mathrm{<mo>\; *\; </mo>}}$ is also involved in the discretization of (20) and (21) at another intersection point on ${\Sigma }_h$. Hence, the whole system $A^{\prime }{u}^{\prime }=b^{\prime }$ is closed.

One can notice that in [32], the same kind of augmented system is considered. Contrary to the AIIB method, no auxiliary nodes are used and the spatial discretization at the grid points at the vicinity of the interface has to be modified “by-hand”.

2.4.3. 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_2^{-1}$. This matrix is not diagonally dominant. For instance, if the discretized solution constraint (23) is inserted in $C_2$ for the row related to $u_I^{\mathrm{<mo>\; *\; </mo>}}$, column $I_{\mathrm{<mo>\; *\; </mo>}}$ has a coefficient ${\alpha }_1$ while the coefficient of the column $J^{\mathrm{<mo>\; *\; </mo>}}$ has a coefficient ${\alpha }_2$. In many situations, ${\alpha }_2$ is greater that ${\alpha }_1$, and the situation can be critical if ${\alpha }_2\gg {\alpha }_1$. The equivalent behaviour can happen with the discretized flux constraint (25). This point is discussed later in Section 3.4.

For the matched interface and boundary (MIB) method, Zhou et al. [29] use a different discretization of the interface conditions which allows an easy algebraic reduction which is directly performed row by row.

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 then the need for a reduction seems to be less critical. More generally, as auxiliary points are in the vicinity of a surface immersed in a volume, their count is one order below the total count of grid points. It suggests that the number of unknowns added to the initial matrix is low compared to the total number of unknowns.

3. 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 [38], and an iterative BiCGSTAB solver [39], preconditioned under a ILUK method [40]. Unless otherwise mentioned, a numerical domain $\left[-\mathrm{1\; <mo>\; ;\; </mo>1}\right]\times \left[-\mathrm{1\; <mo>\; ;\; </mo>1}\right]$ is used for every simulation. Two discrete errors are used.

The discrete relative $L^2$ error in a domain $\Omega$ is defined as:

${\Vert u\Vert }_{L_{rel}^2\left(\Omega \right)}=\frac{{\Vert u-\tilde{u}\Vert }_{L^2\left(\Omega \right)}}{{\Vert \tilde{u}\Vert }_{L^2\left(\Omega \right)}}$ (26)

$={\left(\frac{{\displaystyle \underset{x_I\in \Omega }{\sum }}meas\left({\mathcal{V}}_I\right){\left|u_I-\tilde{u}\left(x_I\right)\right|}^2}{{\displaystyle \underset{x_I\in \Omega }{\sum }}meas\left({\mathcal{V}}_I\right){\left|\tilde{u}\left(x_I\right)\right|}^2}\right)}^{\frac{1}{2}}$ (27)

with $\tilde{u}$ the analytical solution.

The discrete $L^{\infty }$ error is defined as:

${\Vert u\Vert }_{L^{\infty }\left(\Omega \right)}={\max }_{x_I\in \Omega }\left|u_I-\tilde{u}\left(x_I\right)\right|$ (28)

without AIIB method, the present code only reach a first-order spatial error convergence rate for cases with irregular interfaces. With the AIIB method, the code is expected to reach a second-order convergence rate as linear interpolations are used to approximate the solution at the vicinity of the interface.

Only ${\Omega }_0$ is taken into account for the immersed boundary problems. For all cases, $\Sigma$ is the boundary of ${\Omega }_1$ as illustrated in Figure 1. If not explicitly written, $\phi$ and $\psi$ have a default value of 0. For all cases, the analytical solution is imposed for the cells at the boundary of the numerical domain.

3.1. Immersed Boundary Problems

3.1.1. Problem 1

The homogenous 2D Laplace equation is solved. The interface $\Sigma$ is a centered circle of radius $R_1=0.5$ with a Dirichlet condition of $U_1=10$. An analytical solution accounts for the presence of a second circle with a radius $R_2=2$ and $U_2=0$ is imposed on the boundary conditions. The analytical solution is:

$u\left(r\right)=\left(U_2-U_1\right)\frac{1-\ln \left(R_1\right)}{\ln \left(R_2\right)-\ln \left(R_1\right)}\ln \left(r\right)+U_1$ (29)

Accuracy tests are performed with ${\mathbb{L}}_1^1$, ${\mathbb{P}}_1^2$ and ${\mathbb{Q}}_1^2$ interpolations. Figure 4 shows the solution and the error map for a 32 × 32 mesh with ${\mathbb{P}}_1^2$ interpolations. The same results are always obtained with and without algebraic reduction. Figure 5 shows the convergence of the error for the $L^2$ and $L^{\infty }$ norms. For all interpolations, the convergence slopes are approximately 2 for the relative $L^2$ error. For the $L^{\infty }$ error, the slopes are about 1.8. The ${\mathbb{P}}_1^2$ interpolation is the more accurate, followed by the ${\mathbb{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

Figure 4. Solution and error map for problem 1 on a 32 × 32 grid.

Figure 5. Curves of errors for Section 3.1.1.

solver are now benchmarked for the three interpolations with and without algebraic reduction and for the SMP method. Table 1 shows the computational times of the matrix inversions (average time in seconds for 25 matrix inversions) and Table 2 shows the time ratio between the standard and the reduced matrix. Except for the ${\mathbb{Q}}_1^2$ 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. Figure 6 shows the convergence of the ILUK-BiCG-Stab solver for

Table 1. Computational times in seconds for problem 1. Tests are performed with three different interpolations with (red) and without (std) algebraic reduction, and compared to the SMP method.

Figure 6. Residual against iterations of ILUK solver for problem 1 with a 256 × 256 mesh.

the seven interpolations with a 256 × 256 mesh. The same behavior has been observed for other grid resolutions. The type of interpolation does not significantly impact on solver performances. As expected, the number of solver iterations has to be increased to reach a given residual when the number of computational nodes is also increased.

As a conclusion on the interpolation type, ${\mathbb{L}}_1^1$ interpolation seems to be the most desirable as it is generally slightly easier to implement, and requires less memory than a ${\mathbb{Q}}_1^1$ interpolation. To explain that, one can suppose that the boundary constraint forces the points of the stencil to follow a polynomial solution. As the real solution is not polynomial, adding points to the stencil increases the zone where the analytical solution is more difficult to match. For the next problems, ${\mathbb{P}}_1^1$ will be used for constraints on the solution.

3.1.2. Problem 2

The 3D equation $\Delta T=6$ is solved. The solution is $T\left(r\right)=r^2$. The solution is imposed on an immersed centered sphere of radius 0.2. Results of the numerical accuracy test with the spherical inner boundary are presented in Figure 7. The average slope for the $L^2$ norm is 2.33 and increases for the denser meshes. Even if the method has a convergence order similar to the order of the polynomial solution, computer error accuracy cannot be expected because the immersed boundary ${\Sigma }_h$ is an approximation of the initial interface $\Sigma$. The analytical solution

Table 2. Ratio of computational times for reduced and standard matrices for Section 3.1.1.

Figure 7. Curves of errors for problem 2.

(the value at the boundary) is imposed on the boundary, but the boundary location is approximated so an error which has the order of the discretization error of the interface is introduced.

3.1.3. Problem 3

The 3D equation $\Delta T=12r^2$ is solved. The solution is $T\left(r\right)=x^4+y^4+z^4$. Results with and without an immersed interface (a circle of radius $r=1$ ) with the analytical solution imposed on it. The results are presented in Figure 8. For the $L^{\infty }$ norm, the second order is regularly obtained. For the $L^2$ norm, the convergence order tends to a second order with the size of the mesh. 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 ${\Sigma }_h$ improves the accuracy. For both cases, the numerical solution tends to a second order in space.

3.1.4. Problem 4

The 2D equation $\Delta 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$. As can be seen in Figure 9, the global convergence has an average slope of 1.10 and can be explained by the simplicity of the discretization of the Neumann BC.

Figure 8. Curves of errors for problem 3.

Figure 9. Curves of errors for problem 4.

3.2. Immersed Interface Problems

3.2.1. Problem 5

The 2D problem ${\mathcal{P}}_i$ 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 instead of imposing the solution on an approximated shape, we impose a solution jump. If the numerical nodes stay in the correct side of the interface (where they should be if the interface were not approximated), the location of the jump does not influence the numerical solution.

In Figure 10 the same result is obtained with a solution jump $\phi =-1$ on a circular interface such as $u=r^2$ for $r>0.5$ ( ${\Omega }_1$ ) and $u=r^2+1$ otherwise ( ${\Omega }_0$ ).

An equivalent quality of result is obtained (see Figure 11) with $\Sigma$ such as:

Figure 10. The solution and the error for problem 5 with a 32 × 32 mesh.

Figure 11. The solution and the error for problem 5 with a 64 × 64 mesh.

$\{\begin{array}{l}x\left(\alpha \right)=\left(0.5+0.2\sin \left(5\alpha \right)\right)\cos \left(\alpha \right)\\ y\left(\alpha \right)=\left(0.5+0.2\sin \left(5\alpha \right)\right)\sin \left(\alpha \right)\end{array}$ (30)

with $\alpha \in \left[\mathrm{0,2}\pi \right]$. The small stencil of the method allows interfaces with relatively strong curvatures to be used.

3.2.2. Problem 6

For this case, $\left({\mathcal{P}}_i\right)$ is now considered with a discontinuous coefficient a such as $a=10$ in ${\Omega }_0$ and $a=1$ in ${\Omega }_1$, involving the following analytical solution:

$u\left(r\right)=\{\begin{array}{l}{r}^2\text{in}{\Omega }_0\\ \frac{r^2}{10}+\frac{0.9}{4}\text{in}{\Omega }_1\end{array}$ (31)

Accuracy tests are first performed with the interface almost passing by some grid points (the mesh is referred as centered). The interface does not strictly lies on these points, as the shape is shifted by a value of 10⁻¹⁰. This configuration is difficult as the interpolations degenerates. Accuracy tests are then performed with a numerical domain of size 1.0001 (called shifted mesh as the center of $\Sigma$ does not match the center of the numerical domain). In this configuration, the interface never passes by a grid point. The numerical errors on the solution are presented in Figure 12. For the centered mesh, the slope is 1.86 for the $L^2$ and $L^{\infty }$ errors. For the even series, where no geometrical singularity is present, the slope for both errors is 2.04.

Figure 13 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.

3.2.3. Problem 7

The homogenous 2D Laplace equation is considered with the following analytical solution:

$u\left(x\mathrm{,}y\right)=\{\begin{array}{l}0\text{in}{\Omega }_0\\ {\text{e}}^{x}cos\left(y\right)\text{in}{\Omega }_1\end{array}$ (32)

Figure 12. Curves of errors for the centered and shifted meshes for problem 6.

Figure 13. The solution and the error for problem 6 with a 33 × 33 mesh.

where ${\Omega }_0$ and ${\Omega }_1$ are delimited by $\Sigma$ a centered circle of radius 0.5. The jumps coefficients are $\phi =-{\text{e}}^{x}cos\left(y\right)$, $\psi ={\text{e}}^{x}cos\left(y\right)n_x-{\text{e}}^{x}sin\left(y\right)n_y$ and $a=1$.

Figure 14 shows that the convergence for both $L^2$ and $L^{\infty }$ error are of first order only. Figure 15 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. As for the Neumann BC, the simplicity of the discretization of the flux condition seems to be the cause of this low order. Contrary to the discretization of the solution jump, the precise position of the interface is not present which induces a first order of convergence. Generally, other approaches [20] [29] [32] use the precise position of the interface (at the expanse of a larger stencil) and can reach higher orders with all flux conditions.

3.3. Shape Management

3.3.1. 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. Figure 16 shows the accuracy of the solution with respect to the number of points used to discretize the interface which is here a circle. The reference solutions (Figure 5) 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.

3.3.2. The Stanford Bunny

This last case demonstrates how a method with a second-order convergence rate

Figure 14. Convergence of the $L^2$ relative error and the $L^{\infty }$ error for problem 7.

Figure 15. The solution and the $L^2$ relative error for problem 7 with a 32 × 32 mesh.

on spatial accuracy enhances the representation of the boundary condition compared to a first-order method such as the VPM which imposes a prescribed solution at the center of the control volume as shown in Equation (1).

The homogenous Laplace problem with a Dirichlet BC $T_{\Sigma }=10$ is solved on a $60\times 60\times 50$ mesh bounding an obstacle of complex shape (the Stanford bunny). The extension of the solution in ${\Omega }_1$ is used for the post treatment. Thus, all $u_J\mathrm{,}J\in {\mathcal{N}}_1$, are replaced by $u_J^{\mathrm{<mo>\; *\; </mo>}}$. Then, the iso-surface $T=T_{\Sigma }$ gives an idea of the approximation of the boundary condition. Figure 17 shows the iso-surface for a first order method. As can be seen, the shape of the obstacle endures a

Figure 16. Convergence of the error with respect to the number of elements forming the Lagrangian shape for problem 1.

Figure 17. Iso-surface $T=10$ for the Stanford bunny with a first order method.

rasterization effect as the solution is imposed in the entire control volumes. Figure 18 shows the iso-surface for the second order AIIB method. The improvement brought by this approach is straightforward. The graphic processing has been applied to both methods. Figure 19 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.

3.4. 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 $A^{\prime }$, especially on its conditioning. Let us consider an intersection $x_l$ of ${\Sigma }_h$ between two points $x_J\mathrm{,}J\in {\mathcal{N}}_0$ and $x_I\mathrm{,}I\in {\mathcal{N}}_1$. A Dirichlet BC $u_l$ is imposed on it. The constraint constructed with a ${\mathbb{L}}_1^1$ interpolation is $\left(1-\alpha \right)u_J+\alpha u_I^{\mathrm{<mo>\; *\; </mo>}}$, with $\alpha =\frac{x_l-x_J}{x_I-x_J}$. Hence, $\frac{\alpha }{1-\alpha }$ tends to 0 when $x_l$ tends to $x_J$. As the matrix loses its diagonal dominance, solver

Figure 18. Iso-surface $T=10$ for the Stanford bunny with a second order method.

Figure 19. Iso-surface $T=10$ and a slice of the solution.

problems can be encountered. Tseng et al. [10] proposed changing the interpolation by using a new node which is the image of $x_I$ through the interface. In [18] [19], 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 ${\Sigma }_h$.

In our approach, a simple solution is to shift the interface such as the interface matches $x_I$. The loss of accuracy depends on the threshold of distance between $x_I$ and ${\Sigma }_h$ under which the solver is unstable. In our case, test 3.2.2 suggests that a threshold below the error of the discretization of $\Sigma$ into the piecewise linear interface ${\Sigma }_h$ can be taken.

In the case of interface shift, for the Dirichlet BC, an unknown $u_J^{\mathrm{<mo>\; *\; </mo>}}$ 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^{\mathrm{<mo>\; *\; </mo>}}$ and its neighbor unknowns in ${\Omega }_0$.

For the transmission conditions (20)-(21), if $\varphi =0$ and $\psi =0$, no auxiliary unknown is created and the standard finite-volume centered discretization is used. However, for this case, or for $\varphi \ne 0$ and $\psi \ne 0$, our implementation using ILUK preconditioner or a PARDISO direct solver does not necessarily require such methods, even if $\frac{\alpha }{1-\alpha }\approx {10}^{-10}$. Correct solutions are always obtained.

4. Discussion and Conclusions

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.

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. For a non-zero flux jump, the method reaches a first order.

The main advantage of the method is its simplicity and its ease of implementation. First, the method does not require to modify the original discretization of the equation and works only at the level of the discretization matrix. Second, the method has a compact stencil and interfaces with high curvature can be treated with no particular attention.

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 [32]. 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 [41] which are:

${〚\left(a\cdot \nabla u\right)\cdot n〛}_{\Sigma }=\alpha {\bar{u}}_{\mathrm{<mo>\; |\; </mo>}\Sigma }-h\text{on}\Sigma$ (33)

${\overline{\left(a\cdot \nabla u\right)\cdot n}}_{\Sigma }=\beta u-g\text{on}\Sigma$ (34)

where ${\bar{u}}_{\Sigma }=\left(u^{+}+u^{-}\right)/2$ denotes the arithmetic mean of the traces of a quantity on both sides of the interface; $\alpha$, $\beta$, h and g are scalar values which can be chosen in order to obtain Dirichlet, Neumann, Fourier transmission conditions on the interface.

As the method acts on an algebraic level, it can be extended to other implicit approaches where the solution is obtained through the resolution of a linear system. It would be interesting to apply this method to finite element approaches. The interface conditions currently written for finite volumes could be discretized in a finite element manner.

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 [42] [43] 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 than methods using polynomial interpolations so a coupling with the AIIB method is desirable.

Acknowledgements

The authors greatly thank the reviewers for their accurate and pertinent questions and remarks that helped to improve this article. The authors also thank Lisl Weynans for the helpful discussions.

Definition of the Parameters

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Glowinski, R., Pan, T.W., Hesla, T.I. and Joseph, D.D. (1999) A Distributed Lagrange Multiplier/Fictitious Domain Method for Particulate Flows. International Journal of Multiphase Flow, 25, 755-794. https://doi.org/10.1016/S0301-9322(98)00048-2
[2]	Johansen, H. and Colella, P. (1998) A Cartesian Grid Embedded Boundary Method for Poisson’s Equation on Irregular Domains. Journal of Computational Physics, 147, 60-85. https://doi.org/10.1006/jcph.1998.5965
[3]	McCorquodale, P., Colella, P. and Johansen, H. (2001) A Cartesian Grid Embedded Boundary Method for the Heat Equation on Irregular Domains. Journal of Computational Physics, 173, 620-635. https://doi.org/10.1006/jcph.2001.6900
[4]	Ye, T., Mittal, R., Udaykumar, H.S. and Shyy, W. (1999) An Accurate Cartesian Grid Method for Viscous Incompressible Flows with Complex Immersed Boundaries. Journal of Computational Physics, 156, 209-240. https://doi.org/10.1006/jcph.1999.6356
[5]	Peskin, C.S. (1972) Flow Patterns around Heart Valves: A Numerical Method. Journal of Computational Physics, 10, 252-271. https://doi.org/10.1016/0021-9991(72)90065-4
[6]	Peskin, C.S. (2002) The Immersed Boundary Method. Acta Numerica, 11, 479-517. https://doi.org/10.1017/S0962492902000077
[7]	Mohd-Yusof, J. (1997) Combined Immersed Boundary/B-Spline Methods for Simulations of Flows in Complex Geometries. Technical Report, NASA ARS/Stanford University CTR, Stanford.
[8]	Fadlun, E.A., Verzicco, R., Orlandi, P. and Mohd-Yusof, J. (2000) Combined Immersed-Boundary Finite-Difference Methods for Three-Dimensional Complex Flow Simulations. Journal of Computational Physics, 161, 35-60. https://doi.org/10.1006/jcph.2000.6484
[9]	Verzicco, R., Iaccarino, G., Fatica, M. and Orlandi, P. (2001) Flow in an Impeller-Stirred Tank Using an Immersed Boundary Method. AIChE Journal, 50, 1109-1118.
[10]	Tseng, Y.-H. and Ferziger, J.H. (2003) A Ghost-Cell Immersed Boundary Method for Flow in Complex Geometry. Journal of Computational Physics, 192, 593-623. https://doi.org/10.1016/j.jcp.2003.07.024
[11]	Sarthou, A., Vincent, S., Caltagirone, J.-P. and Angot, P. (2008) Eulerian-Lagrangian Grid Coupling and Penalty Methods for the Simulation of Multiphase Flows Interacting with Complex Objects. International Journal for Numerical Methods in Fluids, 56, 1093-1099. https://doi.org/10.1002/fld.1661
[12]	Angot, P., Bruneau, C.-H. and Fabrie, P. (1999) A Penalization Method to Take into Account Obstacles in Incompressible Viscous Flows. Numerische Mathematik, 81, 497-520. https://doi.org/10.1007/s002110050401
[13]	Khadra, K., Angot, P., Parneix, S. and Caltagirone, J.-P. (2000) Fictitious Domain Approach for Numerical Modelling of Navier-Stokes Equations. International Journal for Numerical Methods in Fluids, 34, 651-684. https://doi.org/10.1002/1097-0363(20001230)34:8%3C651::AID-FLD61%3E3.0.CO;2-D
[14]	Lacanette, D., Vincent, S., Sarthou, A., Malaurent, P. and Caltagirone, J.P. (2009) An Eulerian/Lagrangian Method for the Numerical Simulation of Incompressible Convection Flows Interacting with Complex Obstacles: Application to the Natural Convection in the Lascaux Cave. International Journal of Heat and Mass Transfer, 52, 2528-2542. https://doi.org/10.1016/j.ijheatmasstransfer.2008.12.028
[15]	Ramière, I., Angot, P. and Belliard, M. (2007) A General Fictitious Domain Method with Immersed Jumps and Multilevel Nested Structured Meshes. Journal of Computational Physics, 225, 1347-1387. https://doi.org/10.1016/j.jcp.2007.01.026
[16]	Sarthou, A., Vincent, S., Angot, P. and Caltagirone, J.-P. (2008) Finite Volumes for Complex Applications V, Chapter the Sub-Mesh-Penalty Method. Wiley, Hoboken, 633-640.
[17]	Sarthou, A., Vincent, S. and Caltagirone, J.-P. (2020) Consistent Velocity-Pressure Coupling for Direct-Forcing and Penalty Methods. Fluids, 5, 92.
[18]	Gibou, F., Fedkiw, R.P., Cheng, L.-T. and Kang, M. (2002) A Second-Order-Accurate Symmetric Discretization of the Poisson Equation on Irregular Domains. Journal of Computational Physics, 176, 205-227. https://doi.org/10.1006/jcph.2001.6977
[19]	Gibou, F. and Fedkiw, R. (2005) A Fourth Order Accurate Discretization for the Laplace and Heat Equations on Arbitrary Domains, with Applications to the Stefan Problem. Journal of Computational Physics, 202, 577-601. https://doi.org/10.1016/j.jcp.2004.07.018
[20]	Leveque, R.J. and Li, Z. (1994) The Immersed Interface Method for Elliptic Equations with Discontinuous Coefficients and Singular Sources. SIAM Journal of Numerical Analysis, 31, 1001-1025. https://doi.org/10.1137/0731054
[21]	Li, Z.L. and Ito, K. (2006) The Immersed Interface Method. Numerical Solutions of PDEs Involving Interfaces and Irregular Domains. SIAM, Philadelphia. https://doi.org/10.1137/1.9780898717464
[22]	Hou, T.Y., Li, Z.L., Osher, S. and Zhao, H. (1997) A Hybrid Method for Moving Interface Problems with Application to the Hele-Shaw Flow. Journal of Computational Physics, 134, 236-252. https://doi.org/10.1006/jcph.1997.5689
[23]	Lai, M.-C. and Tseng, H-.C. (2008) A Simple Implementation of the Immersed Interface Methods for Stokes Flows with Singular Forces. Computers and Fluids, 37, 99-106. https://doi.org/10.1016/j.compfluid.2007.04.003
[24]	Li, Z. (1998) A Fast Iterative Algorithmfor Elliptic Interfaces Problems. SIAM Journal of Numerical Analysis, 35, 230-254. https://doi.org/10.1137/S0036142995291329
[25]	Li, Z.L., Ito, K. and Lai, M.-C. (2007) An Augmented Approach for Stokes Equations with a Discontinuous Viscosity and Singular Forces. Computers & Fluids, 36, 622-635. https://doi.org/10.1016/j.compfluid.2006.03.003
[26]	Ito, K., Lai, M.-C. and Li, Z.L. (2009) A Well-Conditioned Augmented System for Solving Navier-Stokes Equations in Irregular Domains. Journal of Computational Physics, 228, 2616-2628. https://doi.org/10.1016/j.jcp.2008.12.028
[27]	Fedkiw, R.P., Aslam, T., Merriman, B. and Osher, S. (1999) A Non-Oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (The Ghost Fluid Method). Journal of Computational Physics, 152, 457-492. https://doi.org/10.1006/jcph.1999.6236
[28]	Liu, X.-D., Fedkiw, R.P. and Kang, M. (2000) A Boundary Condition Capturing Method for Poisson’s Equation on Irregular Domains. Journal of Computational Physics, 160, 151-178. https://doi.org/10.1006/jcph.2000.6444
[29]	Zhou, Y.C., Zhao, S., Feig, M. and Wei, G.W. (2006) High Order Matched Interface and Boundary Method for Elliptic Equations with Discontinuous Coefficients and Singular Sources. Journal of Computational Physics, 213, 1-30. https://doi.org/10.1016/j.jcp.2005.07.022
[30]	Zhou, Y.C. and Wei, G.W. (2006) On the Fictitious-Domain and Interpolation Formulations of the Matched Interface and Boundary (MIB) Method. Journal of Computational Physics, 219, 228-246. https://doi.org/10.1016/j.jcp.2006.03.027
[31]	Yu, S.N., Zhou, Y.C. and Wei, G.W. (2007) Matched Interface and Boundary (MIB) Method for Elliptic Problems with Sharp-Edged Interfaces. Journal of Computational Physics, 224, 729-756. https://doi.org/10.1016/j.jcp.2006.10.030
[32]	Cisternino, M. and Weynans, L. (2012) A Parallel Second Order Cartesian Method for Elliptic Interface Problems. Communications in Computational Physics, 12, 1562-1587. https://doi.org/10.4208/cicp.160311.090112a
[33]	Sussman, M. and Fatemi, E. (1999) An Efficient, Interface-Preserving Level Set Redistancing Algorithm and Its Application to Interfacial Incompressible Fluid Flow. SIAM Journal of Scientific Computing, 20, 1165-1191. https://doi.org/10.1137/S1064827596298245
[34]	Shin, S. and Juric, D. (2002) Modelling Three-Dimensional Multiphase Flow Using a Level Contour Reconstruction Method for Front Tracking without Connectivity. Journal of Computational Physics, 180, 427-470. https://doi.org/10.1006/jcph.2002.7086
[35]	Sarthou, A., Vincent, S. and Caltagirone, J.P. (2010) A Second-Order Curvilinear to Cartesian Transformation of Immersed Interfaces and Boundaries. Application to Fictitious Domains and Multiphase Flows. Computers & Fluids, 46, 422-428. https://doi.org/10.1016/j.compfluid.2010.11.008
[36]	Sarthou, A. (2009) Fictitious Domain Methods for the Elliptic and Navier-Stokes Equations. Application to Fluid-Structure Coupling. PhD Thesis, Université Bordeaux, Bordeaux.
[37]	Svard, M. and Nordstrom, J. (2006) On the Order of Accuracy for Difference Approximations of Initial-Boundary Value Problems. Journal of Computational Physics, 352, 219-333. https://doi.org/10.1016/j.jcp.2006.02.014
[38]	Schenk, O. and Gartner, K. (2004) Solving Unsymmetric Sparse Systems of Linear Equations with Pardiso. Journal of Future Generation Computing Systems, 20, 475-487. https://doi.org/10.1016/j.future.2003.07.011
[39]	Gustafsson, I. (1978) On First and Second Order Symmetric Factorization Methods for the Solution of Elliptic Difference Equations. Research Report, Department of Computer Sciences, Chalmersuniversity of Technology and the University of Gateborg, Gateborg.
[40]	Saad, Y. and Schultz, M.H. (1986) Gmres: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing, 7, 856-869. https://doi.org/10.1137/0907058
[41]	Angot, P. (2005) A Unified Fictitious Domain Model for General Embedded Boundary Conditions. Comptes Rendus Mathematique, 341, 683-688. https://doi.org/10.1016/j.crma.2005.09.046
[42]	Randrianarivelo, T.N., Pianet, G., Vincent, S. and Caltagirone, J.-P. (2005) Numerical Modelling of Solid Particle Motion Using a New Penalty Method. International Journal for Numerical Methods in Fluids, 47, 1245-1251. https://doi.org/10.1002/fld.914
[43]	Stéphane, V., De Motta, J.C.B., Sarthou, A., Estivalezes, J.-L., Simonin, O. and Climent, E. (2012) A Lagrangian VOF Tensorial Penalty Method for the DNS of Resolved Particle-Laden Flows. Journal of Computational Physics, 256, 582-614.