# A correction function-based kernel-free boundary integral method for
elliptic PDEs with implicitly defined interfaces

Han Zhou  Wenjun Ying<EMAIL_ADDRESS>School of Mathematical Sciences, MOE-
LSC and Institute of Natural Sciences, Shanghai Jiao Tong University, Minhang,
Shanghai 200240, P. R. China

###### Abstract

This work addresses a novel version of the kernel-free boundary integral
(KFBI) method for solving elliptic PDEs with implicitly defined irregular
boundaries and interfaces. We focus on boundary value problems and interface
problems, which are reformulated into boundary integral equations and solved
with the matrix-free GMRES method. In the KFBI method, evaluating boundary and
volume integrals only requires solving equivalent but much simpler interface
problems in a bounding box, for which fast solvers such as FFTs and geometric
multigrid methods are applicable. For the simple interface problem, a
correction function is introduced for both the evaluation of right-hand side
correction terms and the interpolation of a non-smooth potential function. A
mesh-free collocation method is proposed to compute the correction function
near the interface. The new method avoids complicated derivation for
derivative jumps of the solution and is easy to implement, especially for the
fourth-order method in three space dimensions. Various numerical examples are
presented, including challenging cases such as high-contrast coefficients,
arbitrarily close interfaces and heterogeneous interface problems. The
reported numerical results verify that the proposed method is both accurate
and efficient.

###### keywords:

Elliptic PDEs; Interface problems; Jump conditions; Cartesian grid-based
method; Compact finite difference method

##  1 Introduction

Boundary value problems and interface problems of elliptic partial
differential equations (PDEs) draw much attention due to their wide scientific
and industrial applications, such as viscous incompressible flow [1, 2, 3, 4],
heat transfer [5, 6], biomolecular electrostatics [7, 8], electromagnetics [9,
10], and many others. In practical situations, domain boundaries and material
interfaces are complex and even move with time, making it challenging to
design accurate and efficient numerical methods for these problems.

Body-fitted discretization approaches, such as finite element methods [11, 12,
13, 14], approximate the computational domain with an unstructured mesh, which
conforms with the geometry of boundaries and interfaces to achieve high-order
accuracy. However, it is always difficult and time-consuming to generate high-
quality body-fitted meshes for complex geometries, especially when the
boundary or interface has a substantial motion with time evolution. In
addition, the linear systems generated from the discretization of the PDE on
body-fitted meshes are less structured than that from a Cartesian grid, and
fast solvers such as FFTs and geometric multigrid methods can not be applied.

Immersed methods have been prevalent in recent decades, in which the complex
boundary or interface is immersed into a fixed grid. The pioneering work of
immersed methods is the immersed boundary method (IBM) [15, 16, 17] that was
initially proposed by C. S. Peskin for simulations of cardiac mechanics and
blood flows. In IBM, Peskin uses Lagrangian marker points on the boundary and
regularized Dirac delta functions to approximate the singular force and spread
it into the Eulerian grid. The IBM is quite robust but is restricted to first-
order accuracy due to the non-smoothness of the solution in the vicinity of
the boundary. Motivated by IBM, a number of immersed-type approaches have also
been developed to improve the performance of conventional IBM. Among them are
the immersed interface method (IIM) [18, 19, 20, 21, 22, 23], the ghost-fluid
method (GFM) [24, 25, 26, 27, 28], the matched interface and boundary (MIB)
method [29, 30, 31, 32], the correction function method (CFM) [33, 34, 35],
the Immersed Boundary Smooth Extension (IBSE) method [36, 37]. The methods
mentioned above are mainly based on finite difference discretizations. Since
the finite element method may provide more rigorous convergence analysis,
similar ideas have also been used to develop finite element-based immersed
methods, such as the extended finite element method (XFEM) [38] and the
immersed finite element method (IFEM) [39, 40, 41, 42].

The kernel-free boundary integral (KFBI) is a potential theory-based Cartesian
grid method, which was initially proposed by W. Ying and C. S. Henriquez [43]
as an extension of Mayo’s method [44, 45, 46]. The KFBI method is also an
immersed approach. Unlike traditional boundary integral methods (BIMs)/
boundary element methods (BEMs) [47, 48, 49, 50, 51, 52, 53, 54, 55, 56],
layer and volume potentials are computed by solving equivalent but much
simpler interface problems on a Cartesian grid, and the linear system can be
efficiently solved with FFTs or geometric multigrid methods. Therefore, the
KFBI method has several attractive advantages: (a) no analytical expression of
Green’s function is needed for solving the boundary integral equation; (b)
singular and nearly singular integrals are avoided ;(c) it can be applied to
variable coefficient problems. In the KFBI method, solving the constant
coefficient interface problem is a fundamental building block. In previous
works [57, 58, 59, 43], the simple interface problem is discretized with
standard finite difference methods with a modified right-hand side. The
correction terms for the right-hand side are linear combinations of derivative
jumps $[u],[u_{x}],[u_{y}].[u_{z}],[u_{xx}],\cdots$, which are computed by
repeatedly taking tangential derivatives of the jump values and applying the
local coordinate transformation. The coordinate-transformation method for
derivative jumps is accurate yet complicated when too many derivative terms
are needed, such as for high-order schemes and in three space dimensions, for
example, in [57].

In this work, we present a novel KFBI method that is both simple and accurate
for two and three space dimensional BVPs and interface problems. Motivated by
the correction function method (CFM) [33, 34, 35], we introduce a correction
function in the vicinity of the interface to derive correction terms of the
right-hand side for the constant coefficient interface problem. In order to
solve the local Cauchy problem for the correction function, we proposed a
mesh-free collocation method based on an overlapping surface decomposition for
the interface. Unlike the original CFM [33, 34, 35], no surface quadrature is
required since the collocation method works with the strong form of the Cauchy
problem. The overlapping surface decomposition representation of the interface
also provides a good choice of collocation points such that the resulting
collocation problem is accurate and stable. Another property of the
collocation method is that the discrete system of the collocation problem is a
square one and can be solved accurately, which is different from the original
CFM in that the linear system is overdetermined and needs to be solved in the
least-square sense. The new approach for the constant coefficient interface
problem is built into the KFBI framework to accommodate elliptic BVPs and more
general interface problems. The resulting method is named the correction
function-based KFBI method.

The paper is organized as follows. The governing equations and their boundary
integral equations are described in section 2 and 3. In section 4, the main
idea of the KFBI method is described. The details of the numerical method for
the constant coefficient interface problem are described in section 5. The
algorithm is summarized in section 6. In section 7, numerical results
demonstrating the method with examples are presented. Finally, we discuss the
improvement and advantages of the proposed method in section 8.

##  2 Governing equations

###  2.1 Boundary value problem

Let $\Omega\subset\mathbb{R}^{d},d=2,3$ be a complex domain with smooth
boundary $\Gamma=\partial\Omega$, as illustrated in Figure 1. The BVP of an
elliptic PDE is given by

$\nabla\cdot(\sigma\nabla u)-\kappa u=f,\quad\text{in }\Omega,$ (1)

subject to either the Dirichlet boundary condition or the Neumann boundary
condition

$u=g_{D},\quad\text{or}\quad\sigma\partial_{\boldsymbol{n}}u=g_{N},\quad\text{on
}\Gamma,$ (2)

where $\sigma>0$ is diffusitivity and $\kappa\geq 0$ is reaction coefficient.
In this paper, we assume that $\sigma$ and $\kappa$ are constants.

###  2.2 Interface problem

Let $\Gamma\subset\mathbb{R}^{d},d=2,3$ be a sharp interface that separates a
larger domain $\mathcal{B}\subset\mathbb{R}^{d}$ into two subdomains
$\Omega_{1}$ and $\Omega_{2}$, as illustrated in Figure 1. The interface
problem of an elliptic PDE is given by

$\nabla\cdot(\sigma_{i}\nabla u)-\kappa_{i}u=f_{i},\quad\text{in
}\Omega_{i},\quad i=1,2,$ (3)

subject to two interface jump conditions

$[u]=g_{1},\quad[\sigma\partial_{\boldsymbol{n}}u]=g_{2},\quad\text{on
}\Gamma,$ (4)

and a homogeneous Dirichlet boundary condition on the outer boundary

$u=0,\quad\text{on }\partial\mathcal{B},$ (5)

where $\sigma_{i}>0,i=1,2$ are diffusitivities and $\kappa_{i}\geq 0,i=1,2$
are reaction coefficients. Similarly, we only consider the case that
$\sigma_{i}$ and $\kappa_{i}$ are constants.  Note that the Dirichlet boundary
condition (5) is chosen only for simplicity, since the treatment for boundary
conditions on $\partial\mathcal{B}$ only depends on the finite difference
scheme and is much simpler. Different boundary conditions, such as Neumann and
periodic ones, can also be used.

Here, the boundary/interface $\Gamma$ is assumed to be implicitly defined as
the zero level set of a function. In the case that $\Gamma$ is defined by a
parametric surface or spline, it can also be transformed into an implicit
form.

Figure 1: A schematic of the (a) boundary value problem and (b) interface
problem. Irregular domains and interfaces are embedded into a larger bounding
box, in which a uniform Cartesian grid is used for computation.

##  3 Boundary integral equations

Both the boundary value problem (1) and (2) and the interface problem (3) and
(4) are solved by reformulating them as boundary integral equations.

###  3.1 Boundary value problem

Let $G(\boldsymbol{q},\boldsymbol{p})$ be Green’s function such that for each
fixed $\boldsymbol{p}\in\mathcal{B}$,

$\displaystyle\nabla_{\boldsymbol{q}}\cdot(\sigma(\boldsymbol{q})\nabla_{\boldsymbol{q}}G(\boldsymbol{q},\boldsymbol{p}))-\kappa(\boldsymbol{q})G(\boldsymbol{q};\boldsymbol{p})$
$\displaystyle=\delta(\boldsymbol{q}-\boldsymbol{p}),\quad$
$\displaystyle\text{in }\mathcal{B},$ (6) $\displaystyle
G(\boldsymbol{q},\boldsymbol{p})$ $\displaystyle=0,\quad$
$\displaystyle\text{on }\partial\mathcal{B}.$

Let $\varphi,\psi$ be two density functions. Define the single layer, double
layer, adjoint double layer and hyper-singular integrals, respectively, by

$\displaystyle\mathcal{S}\psi(\boldsymbol{p})$
$\displaystyle=\int_{\Gamma}G(\boldsymbol{q};\boldsymbol{p})\psi(\boldsymbol{q})\,d\boldsymbol{s}_{\boldsymbol{q}},\quad$
$\displaystyle\boldsymbol{p}\in\Gamma,$ (7)
$\displaystyle\mathcal{K}\varphi(\boldsymbol{p})$
$\displaystyle=\int_{\Gamma}\sigma(\boldsymbol{q})\dfrac{\partial
G(\boldsymbol{q};\boldsymbol{p})}{\partial\boldsymbol{n}_{\boldsymbol{q}}}\varphi(\boldsymbol{q})\,d\boldsymbol{s}_{\boldsymbol{q}},\quad$
$\displaystyle\boldsymbol{p}\in\Gamma,$ (8)
$\displaystyle\mathcal{K}^{\prime}\psi(\boldsymbol{p})$
$\displaystyle=\int_{\Gamma}\sigma(\boldsymbol{p})\dfrac{\partial
G(\boldsymbol{q};\boldsymbol{p})}{\partial\boldsymbol{n}_{\boldsymbol{p}}}\psi(\boldsymbol{q})\,d\boldsymbol{s}_{\boldsymbol{q}},\quad$
$\displaystyle\boldsymbol{p}\in\Gamma,$ (9)
$\displaystyle\mathcal{D}\varphi(\boldsymbol{p})$
$\displaystyle=\int_{\Gamma}\sigma(\boldsymbol{p})\sigma(\boldsymbol{q})\dfrac{\partial^{2}G(\boldsymbol{q};\boldsymbol{p})}{\partial\boldsymbol{n}_{\boldsymbol{q}}\partial\boldsymbol{n}_{\boldsymbol{p}}}\varphi(\boldsymbol{q})\,d\boldsymbol{s}_{\boldsymbol{q}},\quad$
$\displaystyle\boldsymbol{p}\in\Gamma.$ (10)

Define the volume integrals by

$\displaystyle\mathcal{G}f(\boldsymbol{p})$
$\displaystyle=\int_{\Omega}G(\boldsymbol{q};\boldsymbol{p})f(\boldsymbol{q})\,d\boldsymbol{q},\quad\boldsymbol{p}\in\Gamma,$
(11) $\displaystyle\partial_{\mathbf{n}}\mathcal{G}f(\boldsymbol{p})$
$\displaystyle=\dfrac{\partial}{\partial\mathbf{n}_{\mathbf{p}}}\int_{\Omega}G(\boldsymbol{q};\boldsymbol{p})f(\boldsymbol{q})\,d\boldsymbol{q},\quad\boldsymbol{p}\in\Gamma.$
(12)

The Dirichlet and Neumann BVP (1) and (2) can be, respectively, reformulated
as boundary integral equations

$\displaystyle(\dfrac{1}{2}+\mathcal{K})\varphi$
$\displaystyle=g_{D}-\mathcal{G}f,$ (13)
$\displaystyle(\dfrac{1}{2}-\mathcal{K}^{\prime})\psi$
$\displaystyle=g_{N}-\sigma\partial_{\boldsymbol{n}}\mathcal{G}f,$ (14)

which are both Fredholm integral equations of the second kind and well-
conditioned.

###  3.2 Interface problem

Let $G_{i}(\boldsymbol{q},\boldsymbol{p}),i=1,2$ be Green’s functions such
that for each fixed $\boldsymbol{p}\in\mathcal{B}$,

$\displaystyle\nabla_{\boldsymbol{q}}\cdot(\sigma_{i}(\boldsymbol{q})\nabla_{\boldsymbol{q}}G_{i}(\boldsymbol{q},\boldsymbol{p}))-\kappa_{i}(\boldsymbol{q})G_{i}(\boldsymbol{q};\boldsymbol{p})$
$\displaystyle=\delta(\boldsymbol{q}-\boldsymbol{p}),\quad$
$\displaystyle\text{in }\mathcal{B},$ (15) $\displaystyle
G_{i}(\boldsymbol{q},\boldsymbol{p})$ $\displaystyle=0,\quad$
$\displaystyle\text{on }\partial\mathcal{B}.$

Similarly, we can define the single layer, double layer, adjoint double layer,
hyper-singular and volume integral operators
$\mathcal{S}_{i},\mathcal{K}_{i},\mathcal{K}^{\prime}_{i},\mathcal{D}_{i},\mathcal{G}_{i},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\partial_{\mathbf{n}}\mathcal{G}_{i}},i=1,2$.
By introducing two unknown densities functions $\varphi=u_{1}$ and
$\psi=\sigma_{2}\partial_{\boldsymbol{n}}u_{2}$, the interface problem (3) can
be reformulated as a system of boundary integral equations

$\displaystyle\varphi-(\mathcal{K}_{1}-\mathcal{K}_{2})\varphi+(\mathcal{S}_{1}-\mathcal{S}_{2})\psi$
$\displaystyle=\dfrac{1}{2}g_{1}+\mathcal{G}_{1}f_{1}+\mathcal{G}_{2}f_{2}+\mathcal{K}_{2}g_{1}-\mathcal{S}_{1}g_{2},$
(16)
$\displaystyle\psi-(\mathcal{D}_{1}-\mathcal{D}_{2})\varphi+(\mathcal{K}^{\prime}_{1}-\mathcal{K}^{\prime}_{2})\psi$
$\displaystyle=-\dfrac{1}{2}g_{2}+\sigma_{1}\partial_{\boldsymbol{n}}\mathcal{G}_{1}f_{1}+\sigma_{2}\partial_{\boldsymbol{n}}\mathcal{G}_{2}f_{2}+\mathcal{D}_{2}g_{1}-\mathcal{K}^{\prime}_{1}g_{2}.$

In the case of $\kappa_{i}=0,i=1,2$ or
$\sigma_{1}/\sigma_{2}=\kappa_{1}/\kappa_{2}$, by dividing the two equations
in (3) by $\sigma_{i}$, respectively, it yields

$\Delta u-\tilde{\kappa}u=\left\\{\begin{aligned}
f_{1}/\sigma_{1},\quad\text{in }\Omega_{1},\\\ f_{2}/\sigma_{2},\quad\text{in
}\Omega_{2},\end{aligned}\right.$ (17)

where $\tilde{\kappa}=0$ or
$\tilde{\kappa}=\kappa_{1}/\sigma_{1}=\kappa_{2}/\sigma_{2}$. Let
$\psi=[\partial_{\boldsymbol{n}}u]$, we may also get a simpler boundary
integral equation

$\dfrac{1}{2}\psi+\mu\mathcal{K}^{\prime}\psi=\dfrac{g_{2}}{\sigma_{1}+\sigma_{2}}+\mu(\mathcal{D}g_{1}+\partial_{\boldsymbol{n}}\mathcal{G}f),$
(18)

where $\mu=(\sigma_{2}-\sigma_{1})/(\sigma_{2}+\sigma_{1})\in(-1,1)$ and
$\mathcal{K},\mathcal{D}$ and $\mathcal{G}$ are the integral operators
associated with the Green’s function of the operator $\Delta-\tilde{\kappa}$.
One may refer to the reference [58] for detailed derivations of the boundary
integral equations.

##  4 Kernel-free boundary integral method

In the kernel-free boundary integral method, values of the boundary integrals
and volume integral at boundary or interface $\Gamma$ are not evaluated with
quadrature methods. Instead, they are evaluated by solving equivalent but much
simpler interface problems for boundary and volume potentials.

For Green’s function $G(\boldsymbol{q},\boldsymbol{p})$, which is associated
with the elliptic operator $\sigma\Delta-\kappa$, define the single layer
potential $-S\psi$, the double layer potential $D\varphi$ and the Newtonian
potential $Nf$ by

$\displaystyle-S\psi(\boldsymbol{p})$
$\displaystyle=-\int_{\Gamma}G(\boldsymbol{q},\boldsymbol{p})\psi(\boldsymbol{q})\,d\boldsymbol{s}_{\boldsymbol{q}},\quad$
$\displaystyle\boldsymbol{p}\in\mathcal{B},$ (19) $\displaystyle
D\varphi(\boldsymbol{p})$
$\displaystyle=\int_{\Gamma}\sigma(\boldsymbol{q})\dfrac{\partial
G(\boldsymbol{q},\boldsymbol{p})}{\partial\boldsymbol{n}_{\boldsymbol{q}}}\varphi(\boldsymbol{q})\,d\boldsymbol{s}_{\boldsymbol{q}},\quad$
$\displaystyle\boldsymbol{p}\in\mathcal{B},$ $\displaystyle
Nf(\boldsymbol{p})$
$\displaystyle=\int_{\Omega}G(\boldsymbol{q};\boldsymbol{p})f(\boldsymbol{q})\,d\boldsymbol{q},\quad$
$\displaystyle\boldsymbol{p}\in\mathcal{B}.$

Then the boundary integrals $\mathcal{S}\psi$, $\mathcal{K}\varphi$,
$\mathcal{K}^{\prime}\psi$, $\mathcal{D}\varphi$ and the volume integrals
$\mathcal{G}f,\partial_{\mathbf{n}}\mathcal{G}f$ coincide with boundary values
or normal derivatives of the potentials $S\psi$, $D\varphi$ and $Nf$. The
above three potential functions are not smooth at $\Gamma$ and they satisfy
equivalent interface problems, by classical potential theory (see [43, 57]).
The equivalent interface problems for the single layer potential $-S\psi$, the
double layer potential $D\varphi$ and the Newton potential $Nf$ can be unified
as

$\left\\{\begin{aligned} &\nabla\cdot(\sigma\nabla u)-\kappa
u=F,\quad&\text{in }\mathcal{B}\setminus\Gamma,\\\ &[u]=\Phi,\quad&\text{on
}\Gamma,\\\ &[\partial_{\boldsymbol{n}}u]=\Psi,\quad&\text{on }\Gamma,\\\
&u=0,\quad&\text{on }\partial\mathcal{B}.\end{aligned}\right.$ (20)

The functions $\Phi$, $\Psi$ and $F$ are specified for each potential by

  * 1.

$-S\psi$: $\Phi=F=0$, $\Psi=\psi$.

  * 2.

$D\psi$: $\Phi=\varphi$, $\Psi=F=0$.

  * 3.

$Nf$: $\Phi=\Psi=0$. $F$ is an arbitrary extension of $f$ to the whole box
$\mathcal{B}$. For simplicity, we set the extended value as zero.

Once the interface problem (20) is solved for the potentials, the boundary
integrals $\mathcal{S}\psi$, $\mathcal{K}\varphi$, $\mathcal{K}^{\prime}\psi$,
$\mathcal{D}\varphi$ and the volume integrals
$\mathcal{G}f,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\partial_{\mathbf{n}}\mathcal{G}f}$
can be obtained from the grid data of these potentials with an interpolation
method.

##  5 Equivalent simple interface problem

Solving the constant coefficient interface problem (20) is an essential part
of the KFBI method. For simplicity, we drop the constant $\sigma$ and proceed
with the following problem

$\left\\{\begin{aligned} &\Delta u-\kappa u=f,\quad&\text{in
}\mathcal{B}\setminus\Gamma,\\\ &[u]=a,\quad&\text{on }\Gamma,\\\
&[\partial_{\boldsymbol{n}}u]=b,\quad&\text{on }\Gamma,\\\ &u=0,\quad&\text{on
}\partial\mathcal{B}.\end{aligned}\right.$ (21)

where $a$, $b$ and $f$ are given data and $\kappa$ is a constant. The right-
hand side $f$ is possibly discontinuous across the interface $\Gamma$. The
constant coefficient interface problem is a much simpler case of the more
general interface problem (3).

###  5.1 Interface representation

In this work, the interface $\Gamma$ is implicitly defined by the level set
function $H$ in the following way

$\Gamma=\\{\mathbf{x}\in\mathbb{R}^{d}|H(\mathbf{x})=0\\},\quad\text{for
}d=2,3.$ (22)

We assume the level set function $H$ is at least $C^{4}$ (for the fourth-order
method) and $|\nabla H|>c_{0}$ for some $c_{0}>0$ near the interface $\Gamma$.
The level set function allows us to easily determine the intersection points
of the surface with grid lines. For instance, if we have an intersection point
on the line segment between two grid nodes $\mathbf{x}_{1}$ and
$\mathbf{x}_{2}$, we expect the values $H(\mathbf{x}_{1})$ and
$H(\mathbf{x}_{2})$ to have opposite signs. By solving the scalar algebraic
equation for $t$ as follows:

$H(t\mathbf{x}_{1}+(1-t)\mathbf{x}_{2})=0,\quad t\in[0,1],$ (23)

using methods such as the Newton method or the bisection method, one can
obtain the coordinates of the intersection point. To compute the unit outward
normal at a surface point, we utilize the gradient of the level set function.
The unit outward normal vector $\mathbf{n}(\mathbf{x})$ at a point
$\mathbf{x}$ on the surface is given by:

$\mathbf{n}(\mathbf{x})=\dfrac{\nabla H(\mathbf{x})}{|\nabla H(\mathbf{x})|}.$
(24)

The method described in [59] is employed in this work for the representation
of the interface using only a subset of intersection points. For each
$r=1,\ldots,d$, $\mathbf{e}_{r}$ represents the $r$-th Cartesian basis vector
in $\mathbb{R}^{d}$, and $\alpha\in(\cos^{-1}(1/\sqrt{d}),\pi/2)$ is a fixed
angle. We define the subset:

$\Gamma_{r}=\\{\mathbf{x}\in\Gamma:|\mathbf{n}(\mathbf{x})\cdot\mathbf{e}_{r}|>\cos\alpha\\},$
(25)

which forms an overlapping surface decomposition of $\Gamma$. The discrete
representation of the interface $\Gamma_{r}$ only considers the intersection
points between $\Gamma_{r}$ and the grid lines aligned with the
$\mathbf{e}_{r}$ direction (refer to Figure 2). We denote the set of these
intersection points as $\Gamma_{r}^{h}$. The union of all sets
$\Gamma_{r}^{h}$ for $r=1,\ldots,d$ is denoted as $\Gamma^{h}$, which
represents the discrete set of points used to approximate $\Gamma$ and
allocate surface degrees of freedom. For more detailed information on the
surface discretization algorithm, please refer to [59].

Figure 2: Illustrations of surface points: (a) points in $\Gamma_{1}^{h}$ (in
2D); (b) points in $\Gamma_{2}^{h}$ (in 2D); (c) points in $\Gamma_{1}^{h}$
(in 3D).

With the help of the overlapping surface decomposition-based discretization,
the interface $\Gamma$ can be locally parameterized by a reference coordinate
plane. Candidate reference planes are

$\displaystyle\Pi_{i}:\\{(x_{1},x_{2},\cdots,x_{d})\in\mathbb{R}^{d}|x_{i}=0\\},\quad
i=1,\cdots,d,$ (26)

for $d=2\text{ or }3$. Suppose at a point $\boldsymbol{x}\in\Gamma$, if the
$i$-th component of the local normal $\boldsymbol{n}(\boldsymbol{x})$ has the
largest absolute value, then we choose $\Pi_{i}$ as the reference plane of
$\Gamma$ near $\boldsymbol{x}$. In such a way, the interpolation stencils on
$\Gamma$ can be easily found with the help of the Cartesian grid on the
reference plane. Numerical integration and interpolation on $\Gamma$ can be
done in a similar way as those on a planer domain as well. We remark that, in
principle, the Cartesian grid used for the representation of $\Gamma$ is not
necessarily chosen as the same one for solving PDEs. In this work, we use the
same Cartesian grid only for simplicity.

###  5.2 Corrected finite difference scheme

For simplicity, the bounding box is assumed to be a unit cube, $i.e.$
$\mathcal{B}=(0,1)^{3}$. Given a positive integer $N$, the domain
$\mathcal{B}$ is uniformly partitioned into a Cartesian grid with mesh
parameter $h=1/N$. Let $P_{i,j,k}$ denote the grid node
$(x_{i},y_{j},z_{k}),i,j,k=0,1,\cdots,$, where $x_{i}=ih$, $y_{j}=jh$ and
$z_{k}=kh$ are node coordinates. For an irregular domain
$\Omega\subset\mathcal{B}$, the interior and exterior grid nodes are defined
as $\Omega_{h}$ and $\Omega_{h}^{C}$, respectively,

$\displaystyle\Omega_{h}$
$\displaystyle=\\{P_{i,j,k}|(x_{i},y_{j},z_{k})\in\Omega,\quad
i,j,k=1,\cdots-1\\},$ (27) $\displaystyle\Omega_{h}^{C}$
$\displaystyle=\\{P_{i,j,k}|(x_{i},y_{j},z_{k})\in\mathcal{B}\setminus\Omega,\quad
i,j,k=1,\cdots-1\\}.$ (28)

In the absence of interfaces, it is known that the following two compact
finite difference schemes (29) and (30) are fourth-order accurate for 2D and
3D cases, respectively.

$\displaystyle-(\dfrac{10}{3h^{2}}+\dfrac{2}{3}\kappa)u_{i,j}+(\dfrac{2}{3h^{2}}-\dfrac{1}{12}\kappa)\sum_{|r|+|s|=1}u_{i+r,j+s}+\dfrac{1}{6h^{2}}\sum_{|r|+|s|=2}u_{i+r,j+s}=f_{i,j}+\dfrac{h^{2}}{12}\Delta
f_{i,j}.$ (29)
$\displaystyle-(\dfrac{25}{6h^{2}}+\dfrac{1}{2}\kappa)u_{i,j,k}+(\dfrac{5}{12h^{2}}-\dfrac{1}{12}\kappa)\sum_{|r|+|s|+|t|=1}u_{i+r,j+s,k+t}+\dfrac{1}{8h^{2}}\sum_{|r|+|s|+|t|=2}u_{i+r,j+s,k+t}$
$\displaystyle+\dfrac{1}{48h^{2}}\sum_{|r|+|s|+|t|=3}u_{i+r,j+s,k+t}=f_{i,j,k}+\dfrac{h^{2}}{12}\Delta
f_{i,j,k}.$ (30)

The two schemes are adopted to derive the corresponding corrected finite
difference schemes for the interface problem (21). We write the finite
difference schemes as a general form

$\sum_{P_{i+r,j+s,k+t}\in\mathcal{S}_{i,j,k}}c_{r,s,t}u_{i+r,j+s,k+t}=F_{i,j,k},$
(31)

where $c_{r,s,t}$ is the coefficient of $u_{i+r,j+s,k+t}$, $F_{i,j,k}$ is the
right-hand side of the finite difference equation and $\mathcal{S}_{i,j,k}$ is
the node set that contains all grid nodes with $c_{r,s,t}\neq 0$ at
$P_{i,j,k}$. Then we define regular nodes $\mathcal{R}_{h}$ and irregular
nodes $\mathcal{I}_{h}$ as follows,

$\displaystyle\mathcal{R}_{h}=\\{P_{i,j,k}|S_{i,j,k}\cap\Omega_{h}=\emptyset\text{
or }S_{i,j,k}\cap\Omega_{h}^{C}=\emptyset\\},$ (32)
$\displaystyle\mathcal{I}_{h}=\\{P_{i,j,k}|S_{i,j,k}\cap\Omega_{h}\neq\emptyset\text{
and }S_{i,j,k}\cap\Omega_{h}^{C}\neq\emptyset\\},$ (33)

At irregular nodes, since the finite difference approximation is taken across
the discontinuity at the interface, large local truncation errors may occur
and result in inaccurate or even divergent results. Precisely, let
$\mathcal{A}_{h}$ denote the difference operator in the finite difference
scheme (31). Suppose the local truncation error is on the order of
$\mathcal{O}(h^{p})$ at a regular node. The local truncation error at an
irregular node is given by

$\displaystyle
E_{h}(x_{i},y_{j},z_{k})=\mathcal{A}_{h}u(x_{i},y_{j},z_{k})-F_{i,j,k}$ (34)
$\displaystyle=$ $\displaystyle\left\\{\begin{aligned}
&\sum_{P_{i+r,j+s,k+t}\in\Omega_{h}\cap\mathcal{S}_{i,j,k}}c_{r,s,t}u^{+}(x_{i+r},y_{j+s},z_{k+t})\\\
&+\sum_{P_{i+r,j+s,k+t}\in\Omega_{h}^{C}\cap\mathcal{S}_{i,j,k}}c_{r,s,t}u^{-}(x_{i+r},y_{j+s},z_{k+t})-F_{i,j,k},\quad
P_{i,j,k}\in\Omega_{h},\\\
&\sum_{P_{i+r,j+s,k+t}\in\Omega_{h}^{C}\cap\mathcal{S}_{i,j,k}}c_{r,s,t}u^{-}(x_{i+r},y_{j+s},z_{k+t})\\\
&+\sum_{P_{i+r,j+s,k+t}\in\Omega_{h}\cap\mathcal{S}_{i,j,k}}c_{r,s,t}u^{+}(x_{i+r},y_{j+s},z_{k+t})-F_{i,j,k},\quad
P_{i,j,k}\in\Omega^{C}_{h},\\\ \end{aligned}\right.$ $\displaystyle=$
$\displaystyle\left\\{\begin{aligned}
&\sum_{P_{i+r,j+s,k+t}\in\mathcal{S}_{i,j,k}}c_{r,s,t}u^{+}(x_{i+r},y_{j+s},z_{k+t})\\\
&+\sum_{P_{i+r,j+s,k+t}\in\Omega_{h}^{C}\cap\mathcal{S}_{i,j,k}}c_{r,s,t}(u^{-}-u^{+})(x_{i+r},y_{j+s},z_{k+t})-F_{i,j,k},\quad
P_{i,j,k}\in\Omega_{h},\\\
&\sum_{P_{i+r,j+s,k+t}\in\mathcal{S}_{i,j,k}}c_{r,s,t}u^{-}(x_{i+r},y_{j+s},z_{k+t})\\\
&+\sum_{P_{i+r,j+s,k+t}\in\Omega_{h}\cap\mathcal{S}_{i,j,k}}c_{r,s,t}(u^{+}-u^{-})(x_{i+r},y_{j+s},z_{k+t})-F_{i,j,k},\quad
P_{i,j,k}\in\Omega^{C}_{h},\\\ \end{aligned}\right.$ $\displaystyle=$
$\displaystyle\left\\{\begin{aligned}
\dfrac{1}{h^{2}}\sum_{P_{i+r,j+s,k+t}\in\Omega_{h}^{C}\cap\mathcal{S}_{i,j,k}}c_{r,s,t}(u^{-}-u^{+})(x_{i+r},y_{j+s},z_{k+t})+\mathcal{O}(h^{p}),\quad
P_{i,j,k}\in\Omega_{h},\\\
\dfrac{1}{h^{2}}\sum_{P_{i+r,j+s,k+t}\in\Omega_{h}\cap\mathcal{S}_{i,j,k}}c_{r,s,t}(u^{+}-u^{-})(x_{i+r},y_{j+s},z_{k+t})+\mathcal{O}(h^{p}),\quad
P_{i,j,k}\in\Omega^{C}_{h},\\\ \end{aligned}\right.$

where $u^{+}$ and $u^{-}$ are two smooth functions that coincide with $u$ in
the domain $\Omega$ and $\Omega^{C}$, respectively. It can be found that the
leading term in the local truncation error at an irregular node is on the
order of $\mathcal{O}(h^{-2})$, which is not acceptable for the sake of
accuracy. The problem can be fixed by including the leading terms of the local
truncation error, as correction terms, into the final finite difference
equations. Define the correction function
$C(\boldsymbol{x})=u^{+}(\boldsymbol{x})-u^{-}(\boldsymbol{x})$. Then, the
corrected finite difference scheme can be written as

$\sum_{P_{i+r,j+s,k+t}\in\mathcal{S}_{i,j,k}}c_{r,s,t}u_{i+r,j+s,k+t}=F_{i,j,k}+C_{i,j,k}.$
(35)

where the correction term $C_{i,j,k}$ is given by

$C_{i,j,k}=\left\\{\begin{aligned} &0,\quad&P_{i,j,k}\in\mathcal{R}_{h},\\\
&-\dfrac{1}{h^{2}}\sum_{P_{i+r,j+s,k+t}\in\Omega_{h}^{C}\cap\mathcal{S}_{i,j,k}}c_{r,s,t}C(x_{i+r},y_{j+s},z_{k+t}),\quad&P_{i,j,k}\in\Omega_{h}\cap\mathcal{I}_{h},\\\
&\dfrac{1}{h^{2}}\sum_{P_{i+r,j+s,k+t}\in\Omega_{h}\cap\mathcal{S}_{i,j,k}}c_{r,s,t}C(x_{i+r},y_{j+s},z_{k+t}),\quad&P_{i,j,k}\in\Omega^{C}_{h}\cap\mathcal{I}_{h}.\\\
\end{aligned}\right.$ (36)

######  Remark 5.1.

If exact values of the correction function $C(\boldsymbol{x})$ are given, then
the local truncation error of the corrected finite difference scheme (35) is
on the order $\mathcal{O}(h^{p})$ at each node. However, it happens only when
the interface coincides with grid nodes and $C(\boldsymbol{x})$ equals the
Dirichlet jump condition $[u]$. In practice, approximate values of the
correction function $C(\boldsymbol{x})$ are used. For the fourth-order method
in this work, the correction function only needs to be approximated with an
error on the order of $\mathcal{O}(h^{5})$ such that the local truncation
error becomes $\mathcal{O}(h^{3})$ at irregular nodes and $\mathcal{O}(h^{4})$
elsewhere.

######  Remark 5.2.

If two interfaces are arbitrarily close, the line segment between two grid
nodes may intersect interfaces more than once (see Figure 3). Let $u^{(i)}$ be
restrictions of the piece-wise smooth solution $u$ in $\Omega^{(i)}$ for
$i=0,1,2$. Denote by $C^{(1)}=u^{(1)}-u^{(0)}$ and $C^{(2)}=u^{(2)}-u^{(0)}$
two correction functions that are solved near $\Gamma^{(1)}$ and
$\Gamma^{(2)}$. For the correction term $C_{i,j}$, the value
$C(x_{i+r},y_{j+s})$ in (36) is computed by

$\displaystyle C(x_{i+r},y_{j+s})$
$\displaystyle=u^{(1)}(x_{i+r},y_{j+s})-u^{(2)}(x_{i+r},y_{j+s})$ (37)
$\displaystyle=u^{(1)}(x_{i+r},y_{j+s})-u^{(0)}(x_{i+r},y_{j+s})+u^{(0)}(x_{i+r},y_{j+s})-u^{(2)}(x_{i+r},y_{j+s})$
$\displaystyle=C^{(1)}(x_{i+r},y_{j+s})-C^{(2)}(x_{i+r},y_{j+s}).$

This is simply adding and subtracting a middle term and is similar to the
technique used in [60].

\begin{overpic}[scale={0.7}]{fig/double_cross.png}
\put(130.0,450.0){$P_{i,j}$} \put(800.0,450.0){$P_{i+r,j+s}$}
\put(200.0,700.0){$\Gamma^{(1)}$} \put(700.0,700.0){$\Gamma^{(2)}$}
\put(200.0,350.0){$\Omega^{(1)}$} \put(700.0,350.0){$\Omega^{(2)}$}
\put(450.0,250.0){$\Omega^{(0)}$} \end{overpic} Figure 3: An illustration of a
line segment intersecting two interfaces.

###  5.3 Local Cauchy problem

Suppose $\Gamma$ is sufficiently smooth and the right-hand side $f$ is also
piece-wise smooth. Denote by $\Omega_{\Gamma}$ a narrow band around $\Gamma$
that covers all irregular nodes. Let $f^{+}$ and $f^{-}$ be two smooth
extension functions of $f$ in $\Omega_{\Gamma}$ from two different sides
$\Omega_{i}$ and $\Omega_{e}$, respectively. Then the function
$\tilde{f}=f^{+}-f^{-}$ is also smooth in $\Omega_{\Gamma}$. The smoothness of
$\tilde{f}$ is relevant to the accuracy of $C(x)$, see 5.3.  Notice that the
correction function $C(\boldsymbol{x})$ satisfies the Cauchy problem

$\displaystyle\Delta C(\boldsymbol{x})-\kappa C(\boldsymbol{x})$
$\displaystyle=\tilde{f}(\boldsymbol{x}),\quad$
$\displaystyle\boldsymbol{x}\in\Omega_{\Gamma},$ (38) $\displaystyle
C(\boldsymbol{x})$ $\displaystyle=a(\boldsymbol{x}),\quad$
$\displaystyle\boldsymbol{x}\in\Gamma,$
$\displaystyle\partial_{\boldsymbol{n}}C(\boldsymbol{x})$
$\displaystyle=b(\boldsymbol{x}),\quad$
$\displaystyle\boldsymbol{x}\in\Gamma.$

Even though the Cauchy problem is known as ill-posed in the sense of Hadmond
that small perturbation in the boundary data will grow exponentially away from
the boundary and it is difficult to obtain a global numerical solution. Since
the correction function is only required at irregular nodes that are close to
the boundary $\Gamma$, we are only interested in the local solution of the
Cauchy problem. In that case, numerical errors can be bounded from above. The
localness of the Cauchy problem also suggests that numerical schemes with a
small stencil, such as compact finite difference schemes, are preferred to
work with for the correction function method.

To locally solve the Cauchy problem (38), we approximate the local solution in
the narrow band $\Omega_{\Gamma}$ with a partition of unity approach. Using
the quasi-uniform point set
$\left\\{\boldsymbol{p}_{i}\right\\}_{i=1}^{N_{p}}\subset\Gamma$ that are
primary points on the boundary $\Gamma$. Let $\Omega_{\Gamma,i}$ be a
neighborhood of the point $\boldsymbol{p}_{i}$. Define $\Omega_{\Gamma}$ as
the union of the neighborhoods

$\Omega_{\Gamma}=\bigcup_{i=1}^{N_{p}}\Omega_{\Gamma,i}.$ (39)

Then $\\{\Omega_{\Gamma,i}\\}_{i=1}^{N_{p}}$ forms an overlapping
decomposition of $\Omega_{\Gamma}$. Note that each $\Omega_{\Gamma,i}$ should
be chosen such that $\Omega_{\Gamma}$ covers all irregular grid nodes.  Unlike
the original CFM [33], where $\Omega_{\Gamma}$ is defined as some particular
grid patches relying on the cut pattern of $\Gamma$ with grid cells, the
current definition of $\Omega_{\Gamma}$ is flexible since it only depends on
the location of surface points. This decomposition gives us a simple way to
represent $C(\boldsymbol{x})$ in $\Omega_{\Gamma}$.

For the partitions $\Omega_{\Gamma,i},i=1,2,\cdots,N_{p}$, define the
compactly supported weight functions $\omega_{i}(\boldsymbol{x})$ such that
$\text{supp}(\omega_{i})=\Omega_{\Gamma,i}$ and

$\sum_{i=1}^{N_{p}}\omega_{i}(\boldsymbol{y})\equiv
1,\quad\boldsymbol{y}\in\Omega_{\Gamma}=\bigcup_{i=1}^{N_{p}}\Omega_{\Gamma,i}.$
(40)

In practice, the weight function $\omega_{i}$ can be constructed in many ways,
such as Shepard’s method [61]. In this work, we use a simple non-smooth weight
function,

$\omega_{i}(\boldsymbol{x})=\left\\{\begin{aligned} 1,&\quad\text{if
$\boldsymbol{p}_{i}$ is the closest point to $\boldsymbol{x}$ for
$i=1,2,\cdots,N_{p}$,}\\\ 0,&\quad\text{otherwise}.\end{aligned}\right.$ (41)

We remark that the smoothness of the weight function has a neglectable effect
on the algorithm. The above simple weight function works very well for all
numerical experiments. Suppose $C_{h,i}(\boldsymbol{x})$ is an approximation
to $C(\boldsymbol{x})$ for $\boldsymbol{x}\in\Omega_{\Gamma,i}$. With the
partition of unity, the complete approximate solution $C_{h}(\boldsymbol{x})$
for $\boldsymbol{x}\in\Omega_{\Gamma}$ is constructed by the linear
combination of local solutions $C_{h,i}$,

$C_{h}(\boldsymbol{x})=\sum_{i=1}^{N_{p}}\omega_{i}(\boldsymbol{x})C_{h,i}(\boldsymbol{x}),\quad\boldsymbol{x}\in\Omega_{\Gamma}.$
(42)

To this end, we restrict the Cauchy problem (38) to the partition
$\Omega_{\Gamma,i}$ and consider numerically solving a sequence of subproblems
for $i=1,2,\cdots N_{p}$,

$\displaystyle\Delta C_{i}(\boldsymbol{x})-\kappa C_{i}(\boldsymbol{x})$
$\displaystyle=\tilde{f}(\boldsymbol{x}),\quad$
$\displaystyle\boldsymbol{x}\in\Omega_{\Gamma,i},$ (43) $\displaystyle
C_{i}(\boldsymbol{x})$ $\displaystyle=a(\boldsymbol{x}),\quad$
$\displaystyle\boldsymbol{x}\in\Gamma\cap\Omega_{\Gamma,i},$
$\displaystyle\partial_{\boldsymbol{n}}C_{i}(\boldsymbol{x})$
$\displaystyle=b(\boldsymbol{x}),\quad$
$\displaystyle\boldsymbol{x}\in\Gamma\cap\Omega_{\Gamma,i}.$

to obtain numerical solutions $C_{h,i}(\boldsymbol{x})$. The restricted
problems (43) are both temporally and spatially local, which explains the
terminology “local Cauchy problem”.

The method is more understandable if one regards the normal direction of
$\Gamma$ as a time variable and the problems (43) as initial-boundary value
problems (IBVPs). Solving the restricted problems for the full Cauchy problem
(38) refers to the explicit method for time-dependent PDEs. In the correction
function method, one does not need to be concerned with the stability of the
explicit method since the solution is only solved by one step away from the
boundary $\Gamma$.

######  Remark 5.3.

To obtain an accurate correction function $C(\boldsymbol{x})$, the right-hand
side $\tilde{f}(\boldsymbol{x})$ should be sufficiently smooth. For a fourth-
order method, $C(\boldsymbol{x})$ is required to be at least $C^{4}$, and,
consequently, $\tilde{f}(\boldsymbol{x})$ is required to be at least $C^{2}$.
Numerically, we can use the same partition of unity approach to represent
$\tilde{f}$ in $\Omega_{\Gamma}$. In each $\Omega_{\Gamma,i}$, $\tilde{f}$ is
replaced by a simple quadratic function using the jump information of $f$ (for
example, in 2D, we use $[f]$, $[f_{x}]$, $[f_{y}]$, $[f_{xx}]$, $[f_{yy}]$,
and $[f_{xy}]$). There are also several different ways to obtain smooth
$f^{+}$ and $f^{-}$, such as the PDE-based method [62] and the partition of
unity extension (PUX) method [63].

####  5.3.1 A mesh-free collocation method

Let $\\{\phi_{l,m,n}(\boldsymbol{x})\\}_{l+m+n\leq p}$ denote the basis of
Taylor polynomial of degree no more than $p$, where the subscripts $l$, $m$
and $n$ are non-negative integers. The elements of the basis are given by, for
example,

$\displaystyle\phi_{0,0,0}(x,y,z)=1,$ (44)
$\displaystyle\phi_{1,0,0}(x,y,z)=x,\quad\phi_{0,1,0}(x,y,z)=y,\quad\phi_{0,0,1}(x,y,z)=z,$
$\displaystyle\phi_{2,0,0}(x,y,z)=x^{2},\quad\phi_{0,2,0}(x,y,z)=y^{2},\quad\phi_{0,0,2}(x,y,z)=z^{2},$
$\displaystyle\phi_{1,1,0}(x,y,z)=xy,\quad\phi_{1,0,1}(x,y,z)=xz,\quad\phi_{0,1,1}(x,y,z)=yz,$
$\displaystyle\cdots.$

The approximate solution $C_{h,i}(\boldsymbol{x})$ is expressed as the linear
combination of the basis

$C_{h,i}(\boldsymbol{x})=\sum_{l+m+n\leq
p}d_{l,m,n}\phi_{l,m,n}(\xi,\eta,\zeta),\quad\boldsymbol{x}\in\Omega_{\Gamma,i},$
(45)

where $\xi$, $\eta$ and $\zeta$ are scaled local coordinate of
$\boldsymbol{x}$. Suppose $\boldsymbol{p}_{i}=(x^{(i)},y^{(i)},z^{(i)})$ is
the center point of the local domain $\Omega_{\Gamma,i}$. The scaled local
coordinate $\tilde{\boldsymbol{x}}=(\xi,\eta,\zeta)$ of
$\boldsymbol{x}=(x,y,z)$ is defined as

$\xi=(x-x^{(i)})/h,\quad\eta=(y-y^{(i)})/h,\quad\zeta=(z-z^{(i)})/h,$ (46)

where $h$ is the mesh parameter. To determine the coefficients $d_{l,m,n}$, we
replace $C_{i}$ with $C_{h,i}$ in the problem (43) and let the equations be
exactly satisfied at multiple points. The resulting method is essentially
mesh-free and falls into the category of collocation methods. Then the chosen
points are called ”collocation points”. Since the problem (43) involves both
bulk PDE and boundary conditions, it involves collocation points both in
$\Omega_{\Gamma,i}$ and $\Gamma\cap\Omega_{\Gamma,i}$. Collocation points can
be classified into three types based on the equations at which they are
satisfied. Let $\boldsymbol{x}_{j}^{pde},j=1,2,\cdots,m_{1}$ be the points in
$\Omega_{\Gamma,i}$ where the PDE is satisfied. Let
$\boldsymbol{x}_{j}^{D},j=1,2,\cdots,m_{2}$ and
$\boldsymbol{x}_{j}^{N},j=1,2,\cdots,m_{3}$ be the points on
$\Gamma\cap\Omega_{\Gamma,i}$ where the Dirichlet and Neumann conditions are
satisfied, respectively. The problem (43) is approximated by the finite-
dimensional problem

$\displaystyle\sum_{l+m+n\leq
p}(\Delta-\kappa)\phi_{l,m,n}(\tilde{\boldsymbol{x}}^{pde}_{j})d_{l,m,n}$
$\displaystyle=\tilde{f}(\boldsymbol{x}_{j}^{pde}),\quad$
$\displaystyle\text{for }j=1,2,\cdots,m_{1}$ (47)
$\displaystyle\sum_{l+m+n\leq
p}\phi_{l,m,n}(\tilde{\boldsymbol{x}}^{D}_{j})d_{l,m,n}$
$\displaystyle=a(\boldsymbol{x}_{j}^{D}),\quad$ $\displaystyle\text{for
}j=1,2,\cdots,m_{2},$ $\displaystyle\sum_{l+m+n\leq
p}\boldsymbol{n}(\boldsymbol{x}^{N}_{j})\cdot\nabla\phi_{l,m,n}(\tilde{\boldsymbol{x}}^{N}_{j})d_{l,m,n}$
$\displaystyle=b(\boldsymbol{x}_{j}^{N}),\quad$ $\displaystyle\text{for
}j=1,2,\cdots,m_{3}.$

The approximate problem (47) forms a linear system

$\mathbf{M}\mathbf{U}=\mathbf{Q},$ (48)

where the unknown vector $\mathbf{U}$ consists of the coefficients
$d_{l,m,n}$.

######  Remark 5.4.

The collocation method is closely related to the local coordinate-
transformation approach used in previous works [43, 59, 58, 57]. The
coordinate-transformation approach can also be viewed as a method for solving
the local Cauchy problem (43) since the correction function
$C(\boldsymbol{x})$ can also be approximated with the jump of derivatives
$[u],[u_{x}],[u_{y}],[u_{z}],[u_{xx}]\cdots$ in terms of Taylor polynomial.
However, the derivation of derivative jumps in the coordinate-transformation
approach involves repeatedly taking tangential derivatives and applying the
chains rule, which requires tedious calculation, especially for high-order and
3D cases. The collocation introduced here is much simpler since applying the
chains rule is not required.

###  5.4 Selection of collocation points

Selecting collocation points is an essential part of the mesh-free collocation
method to ensure accuracy and stability of the algorithm. Different selection
procedures for collocation points result in different systems (48) and
different results. For example, one can choose many collocation points such
that their number is much more than the number of unknowns. In that case, the
linear system (48) becomes overdetermined and can be solved in the least-
square sense, which is similar to the method in [33].  Here, an interpolation-
type method is employed so that each equation in the system (48) is accurately
satisfied. An advantage of using an interpolation-type method is that when a
boundary point coincides with a grid node, the correction function is accurate
at the point since the Dirichlet jump condition $[u]$ is enforced accurately.

Before describing the selection procedure of collocation points, we emphasize
a few key rules:

  1. (a)

Collocation points should be chosen in $\Omega_{\Gamma,i}$ for the PDE and on
$\Gamma\cap\Omega_{\Gamma,i}$ for boundary conditions.

  2. (b)

Collocation points of the same type should be well-separated such that the
resulting linear system is non-singular.

  3. (c)

For each equation in (43), the number of collocation points should be chosen
to meet the formal accuracy requirement.

Rule (a) is a basic requirement for consistency of the collocation method.
Rule (b) is addressed to avoid nearly singular or rank-deficient matrix
$\mathbf{M}$ and to ensure stability of the method. For collocation points of
the same type to be well-separated, the distance between two different points
should have a positive lower bound. Moreover, the number of projections of
these points onto each spatial direction should be sufficiently large such
that the interpolation bases associated with the points can span the
polynomial space. Rule (c) ensures accuracy of the collocation method. Note
that the three equations in (43) have different times of derivatives of
$C_{i}$, and thus a polynomial approximation of $C_{i}$ results in different
orders of accuracy for each equation. Since the equations in (43) are enforced
accurately at collocation points, the collocation problem is also referred to
an interpolation problem. With the error estimation of polynomial
interpolation, one can find that the approximation errors at a point
$\boldsymbol{x}$, away from collocation points, satisfy

$\displaystyle(\Delta-\kappa)C_{h,i}(\boldsymbol{x})-\tilde{f}(\boldsymbol{x})$
$\displaystyle=\mathcal{O}(h^{p-1}),\quad$
$\displaystyle\boldsymbol{x}\in\Omega_{\Gamma,i},$ (49) $\displaystyle
C_{h,i}(\boldsymbol{x})-a(\boldsymbol{x})$
$\displaystyle=\mathcal{O}(h^{p+1}),\quad$
$\displaystyle\boldsymbol{x}\in\Gamma\cap\Omega_{\Gamma,i},$
$\displaystyle\partial_{\boldsymbol{n}}C_{h,i}(\boldsymbol{x})-b(\boldsymbol{x})$
$\displaystyle=\mathcal{O}(h^{p}),\quad$
$\displaystyle\boldsymbol{x}\in\Gamma\cap\Omega_{\Gamma,i}.$

To take into account the consistency and stability requirement and to balance
the approximation errors, we choose collocation points as interpolation points
such that the corresponding Lagrange interpolant on these points has the same
order of accuracy as shown in (49). Precisely, collocation points are chosen
as interpolation points of a polynomial of degree (i) $(p-2)$ for the PDE ;
(ii) $p$ for the Dirichlet boundary condition; and (iii) $(p-1)$ for the
Neumann boundary condition. It should be mentioned that the Lagrangian
interpolant associated with the PDE is in $d$ space dimensions and those for
the boundary conditions are in the $(d-1)$ space dimensions. Therefore, to
choose collocation points for boundary conditions, we first project the
boundary $\Gamma$ into its reference plane locally such that we can find the
local stencil by working with the Cartesian grid on the planer domain. A good
choice of the distribution of collocation points is illustrated in Figure 4.
Similar point selection strategies for multivariate interpolation are used in
[57, 58, 59].

Figure 4: Schematics of collocation points in 3D for (a)Dirichlet boundary
condition; (b)Neumann boundary condition; (c)PDE. For center points of the
local Cauchy problem that are located in the shaded region, collocation points
are marked as black circles. Figures (a) and (b) show the projections of
collocation points on the reference plane.

If collocation points are chosen as above, the number of collocation points
equals the number of degree of freedoms. For example, in three space
dimensions, the collocation equation numbers for the PDE, and Dirichlet and
Neumann boundary condition are $\sum_{i=0}^{p-2}(i+1)(i+2)/2$, $(p+1)(p+2)/2$
and $p(p+1)/2$, respectively. Obviously, it yields

$N^{eqn}=\sum_{i=0}^{p-2}\dfrac{(i+1)(i+2)}{2}+\dfrac{(p+1)(p+2)}{2}+\dfrac{p(p+1)}{2}=\sum_{i=0}^{p}\dfrac{(i+1)(i+2)}{2}=N^{dof}.$
(50)

Then the system (48) is a square one. One can easily verify that similar
results hold for the two space dimensional case as well. The invertibility of
the matrix $\mathbf{M}$ is difficult to prove since it depends on the geometry
of $\Gamma$. Nevertheless, if the collocation points are chosen as
aforementioned, the linear system is always uniquely solvable with a standard
decomposition method, such as the QR decomposition method.

######  Remark 5.5.

We suggest using the scaled local coordinate $(\xi,\eta,\zeta)$ instead of the
original coordinate $(x,y,z)$ for solving the problem (47). It is equivalent
to rescaling the local Cauchy problem such that its characteristic length
changes from $\mathcal{O}(h)$ to $\mathcal{O}(1)$. Thus, the condition number
of the problem $\eqref{eqn:cau-ls}$ is essentially independent of the grid
size $h$ and the accuracy and the stability of the algorithm can be improved
to reduce the effect of round-off error. In the numerical experiments, the
condition number $\text{cond}(\mathbf{M})$ is always on the order of
$10^{2}\sim 10^{3}$ regardless of how small the grid size is.

###  5.5 Extracting boundary data

After solving the linear system of the corrected finite difference scheme, one
can obtain the numerical solution at Cartesian grid nodes. However, in the
KFBI method, one needs to frequently use boundary/interface data, such as
boundary value or normal derivative of the solution, at boundary nodes rather
than Cartesian grid nodes. In order to extract boundary data of the numerical
solution, Lagrange interpolation is used to compute off-grid data. One should
also take into account the jump values of the potential function such that the
Lagrange interpolation has high-order accuracy. The correction function
$C(\boldsymbol{x})$ introduced before now offers a suitable way to take into
account the non-smoothness of the solution. With the correction function, it
is simple to reconstruct smooth data for interpolation using the piece-wise
smooth grid value.

For example, given a boundary point $\boldsymbol{p}\in\Gamma$, we try to
obtain the one-sided limit boundary data of the numerical solution $v_{h}$ in
$\Omega^{+}$. Let $\boldsymbol{q}_{i},i=1,2,\cdots$ be a grid node in the
interpolation stencil near $\boldsymbol{p}$. Suppose the numerical solution
$v_{h}$ is piece-wise smooth and coincides with the smooth functions
$v_{h}^{+}$ and $v_{h}^{-}$ in $\Omega^{+}$ and $\Omega^{-}$, respectively. We
add the correction function $C(\boldsymbol{q}_{i})$ to the grid value
$v_{h}(\boldsymbol{q}_{i})$ if $\boldsymbol{q}_{i}\in\Omega^{-}$ so that the
interpolation data is smooth. By making Taylor expansion at $\boldsymbol{p}$,
it yields

$\displaystyle v_{h}(\boldsymbol{q}_{i})$ $\displaystyle=\sum_{l+m+n\leq
p}\dfrac{p!}{l!m!n!}\xi^{l}\eta^{m}\zeta^{n}\dfrac{\partial^{l+m+n}}{\partial
x^{l}\partial y^{m}\partial
z^{n}}v_{h}^{+}(\boldsymbol{p})+\mathcal{O}(|\boldsymbol{q}_{i}-\boldsymbol{p}|^{p+1}),\quad\text{if
}\boldsymbol{q}_{i}\in\Omega^{+},$ (51) $\displaystyle
v_{h}(\boldsymbol{q}_{i})+C(\boldsymbol{q}_{i})$
$\displaystyle=\sum_{l+m+n\leq
p}\dfrac{p!}{l!m!n!}\xi^{l}\eta^{m}\zeta^{n}\dfrac{\partial^{l+m+n}}{\partial
x^{l}\partial y^{m}\partial
z^{n}}v_{h}^{+}(\boldsymbol{p})+\mathcal{O}(|\boldsymbol{q}_{i}-\boldsymbol{p}|^{p+1}),\quad\text{if
}\boldsymbol{q}_{i}\in\Omega^{-}.$ (52)

where $(\xi,\eta,\zeta)^{T}=\boldsymbol{q}_{i}-\boldsymbol{p}$. Now, by
solving the interpolation problem, the function value and derivatives of
$v_{h}^{+}(\boldsymbol{p})$ are obtained.

##  6 Algorithm Summary

In this section, we summarize the proposed method. We take the boundary
integral equation (16) as an example. The algorithms for the boundary integral
equations (13),(14),(17) are similar. The algorithm is summarized in Algorithm
1.

Algorithm 1 Correction function-based KFBI method

  1. 1.

Compute the right-hand side of (16), in which the integral operators
$\mathcal{S}_{i},\mathcal{K}_{i},\mathcal{K}^{\prime}_{i},\mathcal{D}_{i},\mathcal{G}_{i},\partial_{\mathbf{n}}\mathcal{G}f,i=1,2$
are computed using the same approach in Step 3;

  2. 2.

Given an initial guess for $\varphi$ and $\psi$;

  3. 3.

Compute the integral operators
$\mathcal{S}_{i},\mathcal{K}_{i},\mathcal{K}^{\prime}_{i},\mathcal{D}_{i},i=1,2$
by solving the equivalent interface problem (20);

    1. (a)

Computes the correction function in $\Omega_{\Gamma}$ by solving the local
Cauchy problem (38);

    2. (b)

Compute the correction terms in the right-hand side of (35);

    3. (c)

Solve the linear system of the finite difference scheme (35) with FFT;

    4. (d)

Compute the integral operator values by interpolation from the grid solution;

  4. 4.

Generate the next $\varphi$ and $\psi$ using the GMRES method and repeat Step
3 until the residual is less than a given tolerance.

In each iteration, individually computing the integral operators would require
a total of eight calls of the FFT solver. We stress that the number can be
reduced to two since the terms $\mathcal{K}_{i}\varphi-\mathcal{S}_{i}\psi$
and $\mathcal{D}_{i}\varphi-\mathcal{K}^{\prime}_{i}\psi$ for $i=1$ or $2$ can
be computed by calling the FFT solver only once. By the principle of linear
superposition, one only needs to solve the interface problem for the potential
$D\varphi-S\psi$ and interpolate the function value and normal derivative on
$\Gamma$ to obtain the terms. In this way, it only requires two calls of the
FFT solver in each GMRES iteration.

##  7 Numerical results

In this section, numerical results for boundary and interface problems in both
two and three space dimensions are presented. In the following examples,
irregular domains and interfaces are given in their level-set forms, which
will be specified for each case. Irregular domains and interfaces are embedded
into a bounding box $\mathcal{B}$, which is chosen as a square in 2D and a
cube in 3D. The box $\mathcal{B}$ is uniformly partitioned into $N$ intervals
in each direction for simplicity. The total number of primary boundary points
representing the interface $\Gamma$ is denoted by $N_{b}$.

The following numerical experiments are performed on a personal computer with
3.80 GHz Intel Core i7 processors. The codes for conducting the numerical
experiments are written in C++ computer language. The tolerance in the GMRES
method is fixed as $10^{-10}$. GMRES iteration numbers and CPU times (in
seconds) are reported to quantify the computational complexity. Numerical
errors on the grid node set $\Omega^{h}$ in the $L_{2}$ and maximum norms are
defined as

$\|e\|_{2}=\sqrt{\dfrac{\sum_{\boldsymbol{x}\in\Omega^{h}}|v(\boldsymbol{x})-u(\boldsymbol{x})|^{2}}{M}},\quad\|e\|_{\infty}=\max_{\boldsymbol{x}\in\Omega^{h}}|v(\boldsymbol{x})-u(\boldsymbol{x})|,$
(53)

where $M$ is the number of grid nodes in $\Omega^{h}$, $v$ and $u$ are the
numerical and exact solution, respectively.

###  7.1 Two space dimensional examples

####  7.1.1 Boundary value problem

In the first example, we solve the 2D Dirichlet BVP of the Poisson equation on
a rotated ellipse-shaped domain $\Omega$

$\Omega=\left\\{(x,y)\in\mathbb{R}^{2}:\dfrac{(x\cos\theta+y\sin\theta)^{2}}{a^{2}}+\dfrac{(y\cos\theta-x\sin\theta)^{2}}{b^{2}}<1\right\\},$
(54)

with $a=1,b=0.5,\theta=-\pi/6$. The ellipse is embedded into the bounding box
$\mathcal{B}=[-1.2,1.2]^{2}$. Boundary condition and right-hand side are taken
such that the exact solution satisfies

$u(x,y)=\exp(x)\sin(\cos(\pi/3)x+\sin(\pi/3)y).$ (55)

Numerical results are summarized in Table 1. Nearly fifth-order accuracy in
both the $L_{2}$ and maximum norms can be observed. The increase in
convergence order may be caused by the error of quartic polynomial
interpolation, which is fifth-order accurate and dominates the numerical error
in the vicinity of the boundary. As the grid refines, the GMRES iteration
number is essentially independent of grid size, which is a main advantage of
the present method. Taking into account FFT solvers and boundary operations in
each iteration, the overall computational complexity of the method is given by
$\mathcal{O}(N^{2}\log N+N_{b})$ in two space dimensions.  On coarse grids,
the CPU time scaling is close to $\mathcal{O}(N_{b})$, implying that boundary
operations dominate the computational cost. On finer grids, the CPU time is
roughly linearly proportional to $\mathcal{O}(N^{2}\log N)$, which implies
that the computational cost is dominated by the FFT solver. Isocontours of the
numerical solution are also presented in Figure 6.

Table 1: Numerical results for the Dirichlet BVP of the Poisson equation on an ellipse-shaped domain. grid size | 64$\times$64 | 128$\times$128 | 256$\times$256 | 512$\times$512 | 1024$\times$1024
---|---|---|---|---|---
$N_{b}$ | 116 | 230 | 460 | 918 | 1838
itr no. | 10 | 10 | 9 | 9 | 9
$\|e\|_{2}$ | 7.40E-06 | 1.12E-07 | 3.03E-09 | 6.86E-11 | 2.31E-12
$\|e\|_{\infty}$ | 1.31E-04 | 3.69E-06 | 1.03E-07 | 3.56E-09 | 1.24E-10
CPU time | 3.91E-03 | 6.35E-03 | 1.86E-02 | 5.81E-02 | 2.33E-01

####  7.1.2 Interface problem with multiple interfaces

In the second example, we solve the 2D Poisson interface problem with multiple
disjoint interfaces, which are eight circles and a five-fold star, on the
domain $\mathcal{B}=[-1.7,1.7]^{2}$. The circles are given by

$\Gamma_{m}^{cir}=\left\\{(x,y)\in\mathbb{R}^{2}:(x-\cos(m\pi/4))^{2}+(y-\sin(m\pi/4))^{2}=r^{2}\right\\},\quad
m=1,2,\cdots,8,$ (56)

with $r=0.383$. The five-fold star is given by

$\Gamma^{star}=\left\\{(x,y)\in\mathbb{R}^{2}:\dfrac{x^{2}}{a^{2}}+\dfrac{y^{2}}{b^{2}}=(1.0+\varepsilon\sin(m\arctan(\dfrac{y}{x})))^{2}\right\\},$
(57)

with $a=b=0.514,\varepsilon=0.2,m=5$. Two adjacent interfaces may become very
close to each other, and, as a result, there may be more than one intersection
point between two adjacent grid nodes. Boundary condition, interface condition
and right-hand side are chosen such that the exact solution is given by

$u(x,y)=\left\\{\begin{aligned} &\exp(0.6x+0.8y),\quad&\text{in
}\Omega_{i},\\\ &\sin(\pi(x+1)/2)\sin(\pi(y+1)/2),\quad&\text{in
}\Omega_{e},\end{aligned}\right.$ (58)

where $\Omega_{i}$ denotes the union of the interiors of the circles and the
star and $\Omega_{e}$ denotes the exterior domain. The diffusion coefficients
are chosen as $\sigma_{i}=1$ in $\Omega_{i}$ and $\sigma_{e}=3$ in
$\Omega_{e}$. For this and the following examples, the subscripts $i$ and $e$
represent variables in the interior and exterior regions, respectively.

Numerical results are summarized in Table 2. The solution in both interior and
exterior domains has fourth-order accuracy. The GMRES iteration number is
essentially independent of grid size, even if there are arbitrarily close
interfaces. It can be observed that the iteration number is slightly larger on
the coarsest grid $N=64$. A coarse Cartesian grid may not be able to
accurately capture the geometry of complex interfaces. It affects the well-
conditioned property of the discrete boundary integral equation and causes the
increase in iteration number. Isocontours of the numerical solution are shown
in Figure 6.

Table 2: Numerical results for the Poisson interface problem with multiple touching interfaces. grid size | 64$\times$64 | 128$\times$128 | 256$\times$256 | 512$\times$512 | 1024$\times$1024
---|---|---|---|---|---
$N_{b}$ | 392 | 784 | 1568 | 3136 | 6282
itr no. | 29 | 20 | 20 | 19 | 19
$\|e\|_{\infty,\Omega_{i}}$ | 5.21E-05 | 4.65E-07 | 3.14E-09 | 1.90E-10 | 3.51E-11
$\|e\|_{\infty,\Omega_{e}}$ | 3.94E-05 | 5.73E-07 | 3.19E-09 | 1.89E-10 | 3.51E-11
CPU time | 7.81E-02 | 8.59E-02 | 1.48E-01 | 2.58E-01 | 5.86E-01

Figure 5: Numerical solution to the Dirichlet BVP of the Poisson equation on
an ellipse-shaped domain.

Figure 6: Numerical solution to the Poisson interface problem with multiple
touching interfaces.

###  7.2 Three space dimensional examples

To demonstrate the applicability of the present method, we consider solving
three space dimensional problems.

####  7.2.1 Poisson BVP

This example is the Neumann BVP of the Poisson equation on a torus in 3D. The
torus is given by

$\Omega=\left\\{(x,y,z)\in\mathbb{R}^{3}:(1-\sqrt{x^{2}+y^{2}})^{2}+z^{2}<0.4^{2}\right\\}.$
(59)

The bounding box is taken as $\mathcal{B}=[-1.5,1.5]^{3}$. The boundary
condition and right-hand side are taken such that the exact solution satisfies

$u(x,y,z)=\exp(z)(\cos(2x)+\cos(3y)).$ (60)

Note that the solution to the Poisson Neumann BVP is only determined up to an
additive constant. We first subtract a constant from the right-hand side of
the linear system such that it has zero mean. At the same time, the matrix-
vector products in the GMRES iterations are subtracted by a constant such that
their means are zero. To compute numerical errors, we need to add a constant
to the numerical solution such that it matches the exact solution at a point.

Numerical results and the numerical solution are presented in Table 3 and
Figure 8, respectively. Fourth-order accuracy in both the $L_{2}$ and maximum
norms is reached for the Neumann BVP. In this example, the GMRES iteration
number decreases slightly as the grid refines. Since the discrete linear
system mimics the original well-conditioned BIE and the approximation with a
fine grid is more accurate. We believe that the better approximation property
with a fine grid gives a linear system with a better condition number and is
the main reason that is responsible for the faster convergence of the GMRES
method.

Theoretically, the computational complexity in three space dimensions is
$\mathcal{O}(N^{3}\log N+N_{b})$. The cost of boundary operations is more
important than that in two space dimensions since the polynomial approximation
for the correction function needs more terms in this case. As a result, the
overall computational cost is closer to $\mathcal{O}(N^{2})$ since we have
$N_{b}=\mathcal{O}(N^{2})$.

Table 3: Numerical results for the Neumann BVP of the Poisson equation on a torus. grid size | 64$\times$64$\times$64 | 128$\times$128$\times$128 | 256$\times$256$\times$256 | 512$\times$512$\times$512
---|---|---|---|---
$N_{b}$ | 6168 | 24656 | 98668 | 394548
itr no. | 23 | 21 | 18 | 17
$\|e\|_{2}$ | 2.39E-04 | 2.76E-05 | 2.38E-06 | 1.68E-07
$\|e\|_{\infty}$ | 1.18E-03 | 7.97E-05 | 5.98E-06 | 4.01E-07
CPU time | 1.62E+00 | 6.35E+00 | 2.50E+01 | 1.35E+02

####  7.2.2 Modified Helmholtz BVP

As in the preceding example, we solve the Dirichlet BVP of the modified
Helmholtz equation with $\kappa=100$ on the domain $\Omega$, which is given by

$\Omega=\left\\{(x,y,z)\in\mathbb{R}^{3}:(1+4x^{2})(1+4y^{2})(1+4z^{2})+64xyz+4x^{2}+4y^{2}+4z^{2}<3\right\\},$
(61)

This domain has relatively large curvature and is difficult to be captured by
a coarse grid. The bounding box is taken as $[-0.7,0.7]^{3}$. Boundary
condition and right-hand side are chosen such that the exact solution
satisfies

$u(x,y,z)=\exp(z)(\cos(5x)+\cos(2y)).$ (62)

Numerical results are summarized in Table 4. The numerical solution is
presented in Figure 8. One can observe that the numerical error is large on
the grid $N=64$ and decreases rapidly when the grid is refined to $N=128$. It
can be explained by the fact that the coarse grid $N=64$ may not be able to
fully capture the fast changes of the boundary and cause large errors for
near-interface corrections and surface interpolations. As the grid refines,
the decrease in numerical errors matches the fourth-order accuracy, as
anticipated. The coarse grid with $N=64$ also requires more GMRES iterations
to converge. In each iteration, the CPU time scaling is close to
$\mathcal{O}(N_{b})$ due to the dominance of boundary operations.

Table 4: Numerical results for the Dirichlet BVP of the modified Helmholtz equation.  grid size | 64$\times$64$\times$64 | 128$\times$128$\times$128 | 256$\times$256$\times$256 | 512$\times$512$\times$512
---|---|---|---|---
$N_{b}$ | 7082 | 28482 | 114018 | 455450
itr no. | 14 | 9 | 9 | 9
$\|e\|_{2}$ | 1.14E-05 | 7.34E-08 | 2.82E-09 | 1.20E-10
$\|e\|_{\infty}$ | 1.80E-03 | 6.34E-06 | 5.06E-07 | 2.17E-08
CPU time | 2.98E+00 | 8.52E+00 | 3.82E+01 | 1.75E+02

Figure 7: Numerical solution to the Neumann BVP of the Poisson equation on a
torus.

Figure 8: Numerical solution to the Dirichlet BVP of the modified Helmholtz
equation.

####  7.2.3 Interface problem with high-contrast coefficients

In this example, we solve the Poisson interface equation with a four-atom
molecular-shaped interface in the domain $\mathcal{B}=[-1.2,1.2]^{3}$. The
interface $\Gamma$ is given by

$\Gamma=\left\\{\boldsymbol{x}=(x,y,z)\in\mathbb{R}^{3}:\sum_{k=1}^{4}\exp(-\dfrac{|\boldsymbol{x}-\boldsymbol{x}_{k}|^{2}}{r^{2}})=0.6\right\\},$
(63)

with $\boldsymbol{x}_{1}=(\sqrt{3}/3,0,-\sqrt{6}/12)$,
$\boldsymbol{x}_{2}=(-\sqrt{3}/6,0.5,-\sqrt{6}/12)$,
$\boldsymbol{x}_{3}=(-\sqrt{3}/6,-0.5,-\sqrt{6}/12)$ and
$\boldsymbol{x}_{4}=(0,0,\sqrt{6}/4)$.

$u(x,y,z)=\left\\{\begin{aligned}
&\sin^{2}(2x)\cos^{2}(2y)\cos(z),\quad&\text{in the interior }\Omega_{i},\\\
&\cos(x)\cos(y)\cos(z),\quad&\text{in the exterior
}\Omega_{e}.\end{aligned}\right.$ (64)

The coefficient ratio $\sigma_{e}/\sigma_{i}$ varies from $10$ to $10^{4}$,
and its effect on the performance of the present method is studied in this
example. This effect was also studied by [29, 14, 35]. The numerical solution
is shown in Figure 10. According to the numerical results presented in Table
5, high-contrast coefficients only have a small effect on the numerical
accuracy, even for the extreme case $\sigma_{e}/\sigma_{i}=10^{4}$. The GMRES
iteration number is slightly affected by the coefficient ratio on coarse
grids. As the grid refines, the GMRES iteration number is rather stable and is
independent of the coefficient ratio. This is also due to the fact that a fine
grid has a better approximation property, as aforementioned.

Table 5: Numerical results for the Poisson interface problems with varying coefficient ratios.  $\sigma_{e}:\sigma_{i}$ | N | itr no. | $\|e\|_{2,\Omega_{i}}$ | $\|e\|_{\infty,\Omega_{i}}$ | $\|e\|_{2,\Omega_{e}}$ | $\|e\|_{\infty,\Omega_{e}}$
---|---|---|---|---|---|---
$10:1$ | 128 | 11 | 5.15E-08 | 2.18E-07 | 1.43E-08 | 2.33E-07
256 | 10 | 3.25E-09 | 1.22E-08 | 9.03E-10 | 1.33E-08
512 | 10 | 2.04E-10 | 7.59E-10 | 5.76E-11 | 8.08E-10
$10^{2}:1$ | 128 | 13 | 5.78E-08 | 2.74E-07 | 1.89E-08 | 2.92E-07
256 | 11 | 3.63E-09 | 1.61E-08 | 1.18E-09 | 1.71E-08
512 | 10 | 2.29E-10 | 1.02E-09 | 7.59E-11 | 1.06E-09
$10^{4}:1$ | 128 | 14 | 5.86E-08 | 2.81E-07 | 1.94E-08 | 2.99E-07
256 | 11 | 3.68E-09 | 1.66E-08 | 1.21E-09 | 1.75E-08
512 | 10 | 2.33E-10 | 1.05E-09 | 7.83E-11 | 1.09E-09

####  7.2.4 Interface problem with arbitrarily close interfaces

In this case, we solve the Poisson interface problem with the presence of
arbitrarily close interfaces in three space dimensions. Interfaces are taken
as a torus and an ellipsoid. The torus-shaped interface $\Gamma^{tor}$ is
given by the boundary of the domain $\Omega$ defined in (59). The ellipsoid-
shaped interface is given by

$\Gamma^{ell}=\left\\{(x,y,z)\in\mathbb{R}^{3}:\dfrac{x^{2}}{a^{2}}+\dfrac{y^{2}}{b^{2}}+\dfrac{z^{2}}{c^{2}}=1\right\\},$
(65)

with $a=b=0.6,c=1$. The two interfaces are very close to each other near the
curve

$S=\left\\{(x,y,z)\in\mathbb{R}^{3}:\dfrac{x^{2}}{a^{2}}+\dfrac{y^{2}}{b^{2}}=1,\quad
z=0\right\\}.$ (66)

In this configuration, since the curve $S$ is a dimension one object, the
number of multi-intersection grid line segments, a grid line segment that
intersects with interfaces multiple times, is on the order of
$\mathcal{O}(N)$. The problem is challenging for classical body-fitted
approaches, since it is nearly impossible to be resolved by a body-fitted mesh
since the two interface is too close. The bounding box $\mathcal{B}$ is taken
as $[-1.5,1.5]^{3}$. Boundary condition, interface condition and right-hand
side are chosen such that the exact solution reads

$u(x,y,z)=\left\\{\begin{aligned}
&\sin^{2}(2x)\cos^{2}(2y)\cos(z),\quad&\text{in torus }\Omega_{i,1},\\\
&\exp(z)(\cos(2x)+\cos(3y)),\quad&\text{in ellipsoid }\Omega_{i,2},\\\
&\cos(x)\cos(y)\cos(z),\quad&\text{in exterior region
}\Omega_{e}.\end{aligned}\right.$ (67)

The coefficients are chosen as $\sigma_{i}=1$ in
$\Omega_{i,1}\cup\Omega_{i,2}$ and $\sigma_{e}=3$ in $\Omega_{e}$. Numerical
results are summarized in Table 6. The numerical solution is visualized and
shown in Figure 10. It is observed that fourth-order accuracy is achieved in
all regions, except for an accuracy loss on the coarsest grid $N=64$ due to
similar reasons that were mentioned before.

Table 6: Numerical results for the Poisson interface problem with a touching curve. grid size | 64$\times$64$\times$64 | 128$\times$128$\times$128 | 256$\times$256$\times$256 | 512$\times$512$\times$512
---|---|---|---|---
$N_{b}$ | 8738 | 34866 | 139454 | 557794
itr no. | 22 | 21 | 19 | 17
$\|e\|_{\infty,\Omega_{i,1}}$ | 2.33E-03 | 6.75E-07 | 4.38E-08 | 2.13E-09
$\|e\|_{\infty,\Omega_{i,2}}$ | 5.56E-05 | 7.72E-07 | 4.39E-08 | 2.11E-09
$\|e\|_{\infty,\Omega_{e}}$ | 3.24E-03 | 5.97E-07 | 3.65E-08 | 2.06E-09
CPU time | 5.55E+00 | 2.08E+01 | 7.07E+01 | 2.86E+02

Figure 9: Numerical solution to the Poisson interface problem with varying
coefficient ratios.

Figure 10: Numerical solution to the Poisson interface problem with a touching
curve.

####  7.2.5 Heterogeneous interface problem

In the final example, we consider the heterogeneous interface problem in three
space dimensions. Interfaces are taken as three spheres with radius $r=0.7$
whose centers are chosen as $\boldsymbol{x}_{1}=(0.5,0.5,0.5)$,
$\boldsymbol{x}_{2}=(-0.5,-0.5,0.5)$ and $\boldsymbol{x}_{3}=(0.5,-0.5,-0.5)$,
respectively. The coefficients on each side of the interfaces are given as

$\sigma_{i}=1,\quad\kappa_{i}=0,\quad\sigma_{e}=4,\quad\kappa_{e}=10,$ (68)

such that the unknown function $u$ satisfies the Poisson equation in the
interior region and the modified Helmholtz equation in the exterior region. It
is called a heterogeneous interface problem since the elliptic differential
operator on two sides of the interfaces are of different types. The
heterogeneous interface problem is a linearized version of the Poisson-
Boltzmann equation, which appears in the Poisson-Boltzmann theory in
biophysics for modeling solvated biomolecular systems. Table 7 and Figure 11
show the numerical results and the visualization of the numerical solution,
respectively. Once again we observe the fourth-order convergence in both
regions. The number of GMRES iterations is essentially independent of the grid
size.

Table 7: Numerical results for the heterogeneous interface problem. grid size | 64$\times$64$\times$64 | 128$\times$128$\times$128 | 256$\times$256$\times$256 | 512$\times$512$\times$512
---|---|---|---|---
$N_{b}$ | 6987 | 27939 | 111774 | 447264
itr no. | 21 | 20 | 20 | 20
$\|e\|_{\infty,\Omega_{i}}$ | 6.71E-06 | 4.79E-07 | 3.32E-08 | 2.26E-09
$\|e\|_{\infty,\Omega_{e}}$ | 3.89E-06 | 2.65E-07 | 1.73E-08 | 1.15E-09
CPU time | 1.34E+00 | 6.59E+00 | 4.08E+01 | 3.20E+02
Figure 11: Numerical solution to the heterogeneous interface problem.

##  8 Discussion

This work proposes a new version of the kernel-free boundary integral method
for solving elliptic partial differential equations in two and three space
dimensions with high accuracy. The KFBI method solves boundary and interface
problems with their boundary integral formulations. It computes boundary and
volume integrals by solving equivalent interface problems with fast PDE
solvers and then obtains boundary values by interpolation.

The equivalent interface problems are simpler than the original problem and
are essential for the KFBI method. To accommodate the jump conditions in the
interface, a correction function is introduced in the vicinity of the
interface to derive corrected finite difference schemes and the boundary
interpolation scheme. Unlike the original KFBI method, which applies a local
coordinate transformation to calculate correction terms, the new approach
obtains correction terms by solving a local Cauchy problem for the correction
function. The local Cauchy problem is solved with a mesh-free collocation
method, for which we also proposed a strategy to choose collocation points
such that the resulting linear system is accurate and stable. The resulting
method avoids repeatedly taking tangent derivatives on the jump condition and
significantly simplifies the derivation procedure.

The presented method is efficient and accurate, which is demonstrated through
several challenging numerical experiments. The efficiency of the method relies
on the well-conditionness of the boundary integral equations and the
applicability of fast PDE solvers (FFT, geometric multigrid methods) on a
Cartesian grid. Even though the presented numerical results are based on a
fourth-order implementation of the method, the method can be extended to
arbitrary accuracy in principle [33].

Finally, we emphasize that the present method is designed for implicitly
defined interfaces with level-set formulations. Although this work uses an
analytic expression of the level-set function, extending the method to cases
when the level-set function is only given at Cartesian grid nodes is
straightforward. It may have advantages for solving moving interface problems
and free boundary problems when combined with the level-set method [64, 65].

## Acknowledgments

This work is financially supported by the Shanghai Science and Technology
Innovation Action Plan in Basic Research Area (Project No. 22JC1401700). It is
also partially supported by the National Key R&D Program of China (Project No.
2020YFA0712000), the Strategic Priority Research Program of Chinese Academy of
Sciences (Grant No. XDA25010405) and the National Natural Science Foundation
of China (Grant No. DMS-11771290).

## References

  * Quartapelle [1993] L. Quartapelle, Numerical solution of the incompressible Navier-Stokes equations, volume 113, Springer Science & Business Media, 1993.
  * Chorin [1968] A. J. Chorin,  Numerical solution of the navier-stokes equations,  Mathematics of computation 22 (1968) 745–762.
  * Calhoun [2002] D. Calhoun,  A cartesian grid method for solving the two-dimensional streamfunction-vorticity equations in irregular regions,  Journal of Computational Physics 176 (2002) 231–275.
  * Russell and Jane Wang [2003] D. Russell, Z. Jane Wang,  A cartesian grid method for modeling multiple moving objects in 2d incompressible viscous flow,  Journal of Computational Physics 191 (2003) 177–205.
  * Bergman et al. [2011] T. L. Bergman, A. S. Lavine, F. P. Incropera, D. P. DeWitt, Introduction to heat transfer, John Wiley & Sons, 2011.
  * Patankar [2018] S. V. Patankar, Numerical heat transfer and fluid flow, CRC press, 2018.
  * Davis and McCammon [1990] M. E. Davis, J. A. McCammon,  Electrostatics in biomolecular structure and dynamics,  Chemical Reviews 90 (1990) 509–521.
  * Honig and Nicholls [1995] B. Honig, A. Nicholls,  Classical electrostatics in biology and chemistry,  Science 268 (1995) 1144–1149.
  * Rothwell and Cloud [2018] E. J. Rothwell, M. J. Cloud, Electromagnetics, CRC press, 2018.
  * Bondeson et al. [2012] A. Bondeson, T. Rylander, P. Ingelström, Computational electromagnetics, Springer, 2012.
  * Babuška [1970] I. Babuška,  The finite element method for elliptic equations with discontinuous coefficients,  Computing 5 (1970) 207–213.
  * Guyomarc’h et al. [2009] G. Guyomarc’h, C.-O. Lee, K. Jeon,  A discontinuous galerkin method for elliptic interface problems with application to electroporation,  Communications in Numerical Methods in Engineering 25 (2009) 991–1008.
  * Bramble and King [1996] J. H. Bramble, J. T. King,  A finite element method for interface problems in domains with smooth boundaries and interfaces,  Advances in Computational Mathematics 6 (1996) 109–138.
  * Wang [2004] W. C. Wang,  A jump condition capturing finite difference scheme for elliptic interface problems,  SIAM Journal on Scientific Computing 25 (2004) 1479–1496.
  * Peskin [1972] C. S. Peskin,  Flow patterns around heart valves: A numerical method,  Journal of Computational Physics 10 (1972) 252–271.
  * Peskin [2002] C. S. Peskin,  The immersed boundary method,  Acta Numerica 11 (2002) 479–517.
  * Peskin [1977] C. S. Peskin,  Numerical analysis of blood flow in the heart,  Journal of Computational Physics 25 (1977) 220–252.
  * Tan et al. [2008] Z. Tan, D. V. Le, Z. Li, K. M. Lim, B. C. Khoo,  An immersed interface method for solving incompressible viscous flows with piecewise constant viscosity across a moving elastic membrane,  Journal of Computational Physics 227 (2008) 9955–9983.
  * Deng et al. [2003] S. Deng, K. Ito, Z. Li,  Three-dimensional elliptic solvers for interface problems and applications,  Journal of Computational Physics 184 (2003) 215–243.
  * Zhilin [1998] L. I. Zhilin,  A fast iterative algorithm for elliptic interface problems,  SIAM Journal on Numerical Analysis 35 (1998) 230–254.
  * Leveque and Li [1994] R. J. Leveque, Z. Li,  Immersed interface method for elliptic equations with discontinuous coefficients and singular sources,  SIAM Journal on Numerical Analysis 31 (1994) 1019–1044.
  * Leveque and Li [1997] R. J. Leveque, Z. Li,  Immersed interface methods for stokes flow with elastic boundaries or surface tension,  SIAM Journal of Scientific Computing 18 (1997) 709–735.
  * Li and Lai [2001] Z. Li, M. C. Lai,  The immersed interface method for the navier-stokes equations with singular forces,  Journal of Computational Physics 171 (2001) 822–842.
  * Fedkiw et al. [1999a] R. P. Fedkiw, T. Aslam, S. Xu,  The ghost fluid method for deflagration and detonation discontinuities,  Journal of Computational Physics 154 (1999a) 393–427.
  * Fedkiw et al. [1999b] R. P. Fedkiw, T. Aslam, B. Merriman, S. Osher,  A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method),  Journal of Computational Physics 152 (1999b) 457–492.
  * Nguyen et al. [2001] D. Q. Nguyen, R. P. Fedkiw, M. Kang,  A boundary condition capturing method for incompressible flame discontinuities,  Journal of Computational Physics 172 (2001) 71–98.
  * Liu et al. [2000] X. D. Liu, R. P. Fedkiw, M. Kang,  A boundary condition capturing method for poisson’s equation on irregular domains,  Journal of Computational Physics 160 (2000) 151–178.
  * Luo et al. [2006] X. Y. Luo, M. J. Ni, A. Ying, M. Abdou,  A boundary condition capturing method for multiphase flow with phase change,  Collection of Technical Papers - 44th AIAA Aerospace Sciences Meeting 20 (2006) 15443–15451.
  * Zhou et al. [2006] Y. C. Zhou, S. Zhao, M. Feig, G. W. Wei,  High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources,  Journal of Computational Physics 213 (2006) 1–30.
  * Yu et al. [2007] S. Yu, Y. Zhou, G. W. Wei,  Matched interface and boundary (mib) method for elliptic problems with sharp-edged interfaces,  Journal of Computational Physics 224 (2007) 729–756.
  * Wang et al. [2015] B. Wang, K. Xia, G. W. Wei,  Matched interface and boundary method for elasticity interface problems,  Journal of Computational and Applied Mathematics 285 (2015) 203–225.
  * Feng et al. [2019] H. Feng, G. Long, S. Zhao,  An augmented matched interface and boundary (mib)method for solving elliptic interface problem,  Journal of Computational and Applied Mathematics 361 (2019) 426–443.
  * Marques et al. [2011] A. N. Marques, J. C. Nave, R. R. Rosales,  A correction function method for poisson problems with interface jump conditions,  Journal of Computational Physics 230 (2011) 7567–7597.
  * Marques et al. [2019] A. N. Marques, J. C. Nave, R. R. Rosales,  Imposing jump conditions on nonconforming interfaces for the correction function method: A least squares approach,  Journal of Computational Physics 397 (2019) 108869.
  * Marques et al. [2017] A. N. Marques, J. C. Nave, R. R. Rosales,  High order solution of poisson problems with piecewise constant coefficients and interface jumps,  Journal of Computational Physics 335 (2017) 497–515.
  * Stein et al. [2017] D. B. Stein, R. D. Guy, B. Thomases,  Immersed boundary smooth extension (ibse): A high-order method for solving incompressible flows in arbitrary smooth domains,  Journal of Computational Physics 335 (2017) 155–178.
  * Stein et al. [2016] D. B. Stein, R. D. Guy, B. Thomases,  Immersed boundary smooth extension: A high-order method for solving pde on arbitrary smooth domains using fourier spectral methods,  Journal of Computational Physics 304 (2016) 252–274.
  * Moës et al. [1999] N. Moës, J. Dolbow, T. Belytschko,  A finite element method for crack growth without remeshing,  International Journal for Numerical Methods in Engineering 46 (1999) 131–150.
  * Li et al. [2003] Z. Li, T. Lin, X. Wu,  New cartesian grid methods for interface problems using the finite element formulation,  Numerische Mathematik 96 (2003) 61–98.
  * Gong et al. [2008] Y. Gong, B. Li, Z. Li,  Immersed-interface finite-element methods for elliptic interface problems with nonhomogeneous jump conditions,  SIAM Journal on Numerical Analysis 46 (2008) 472–495.
  * Hou and Liu [2005] S. Hou, X. D. Liu,  A numerical method for solving variable coefficient elliptic equation with interfaces,  Journal of Computational Physics 202 (2005) 411–445.
  * Guo and Lin [2020] R. Guo, T. Lin,  An immersed finite element method for elliptic interface problems in three dimensions,  Journal of Computational Physics 414 (2020) 109478.
  * Ying and Henriquez [2007] W. Ying, C. S. Henriquez,  A kernel-free boundary integral method for elliptic boundary value problems,  Journal of Computational Physics 227 (2007) 1046–1074.
  * Mayo [1984] A. Mayo,  The fast solution of poisson’s and the biharmonic equations on irregular regions,  SIAM Journal on Numerical Analysis 21 (1984) 285–299.
  * Mayo [1985] A. Mayo,  Fast high order accurate solution of laplace’s equation on irregular regions,  SIAM Journal on Scientific and Statistical Computing 6 (1985) 144–157.
  * McKenney et al. [1995] A. McKenney, L. Greengard, A. Mayo,  A fast poisson solver for complex geometries,  Journal of Computational Physics 118 (1995) 348–355.
  * Beale [2004] J. T. Beale,  A grid-based boundary integral method for elliptic problems in three dimensions,  SIAM Journal on Numerical Analysis 42 (2004) 599–620.
  * Beale and Lai [2001] J. T. Beale, M.-C. Lai,  A method for computing nearly singular integrals,  SIAM Journal on Numerical Analysis 38 (2001) 1902–1925.
  * Greengard and Kropinski [1998] L. Greengard, M. C. Kropinski,  An integral equation approach to the incompressible navier-stokes equations in two dimensions,  SIAM Journal of Scientific Computing 20 (1998) 318–336.
  * Carrier et al. [1988a] J. Carrier, L. GreenGardi, V. Rokhlin,  a fast adaptive multipole algorithm,  SIAM J. SCl. STAT. COMPUT. 9 (1988a) 669–686.
  * Carrier et al. [1988b] J. Carrier, L. Greengard, V. Rokhlin,  A fast adaptive multipole algorithm for particle simulations,  SIAM Journal on Scientific and Statistical Computing 9 (1988b) 669–686.
  * Greengard and Rokhlin [1987] L. Greengard, V. Rokhlin,  A fast algorithm for particle simulations,  Journal of Computational Physics 73 (1987) 325–348.
  * Liu et al. [1981] P. L.-F. Liu, A. H.-D. Cheng, J. A. Liggett, J. H. Lee,  Boundary integral equation solutions to moving interface between two fluids in porous media,  Water Resources Research 17 (1981) 1445–1452.
  * Ying et al. [2006] L. Ying, G. Biros, D. Zorin,  A high-order 3d boundary integral equation solver for elliptic pdes in smooth domains,  Journal of Computational Physics 219 (2006) 247–275.
  * Klaseboer and Khoo [2004] E. Klaseboer, B. C. Khoo,  Boundary integral equations as applied to an oscillating bubble near a fluid-fluid interface,  Computational Mechanics 33 (2004) 129–138.
  * Tlupova and Cortez [2009] S. Tlupova, R. Cortez,  Boundary integral solutions of coupled stokes and darcy flows,  Journal of Computational Physics 228 (2009) 158–179.
  * Xie and Ying [2020] Y. Xie, W. Ying,  A fourth-order kernel-free boundary integral method for implicitly defined surfaces in three space dimensions,  Journal of Computational Physics 415 (2020) 109526.
  * Ying and Wang [2014] W. Ying, W. C. Wang,  A kernel-free boundary integral method for variable coefficients elliptic pdes,  Communications in Computational Physics 15 (2014) 1108–1140.
  * Ying and Wang [2013] W. Ying, W. C. Wang,  A kernel-free boundary integral method for implicitly defined surfaces,  Journal of Computational Physics 252 (2013) 606–624.
  * Ying [2018] W. Ying,  A cartesian grid-based boundary integral method for an elliptic interface problem on closely packed cells,  Communications in Computational Physics 24 (2018) 1196–1220.
  * Shepard [1968] D. Shepard,  Two- dimensional interpolation function for irregularly- spaced data,  Proc 23rd Nat Conf (1968) 517–524.
  * Aslam [2004] T. D. Aslam,  A partial differential equation approach to multidimensional extrapolation,  Journal of Computational Physics 193 (2004) 349–355.
  * Fryklund et al. [2018] F. Fryklund, E. Lehto, A.-K. K. Tornberg,  Partition of unity extension of functions on complex domains,  Journal of Computational Physics 375 (2018) 57–79.
  * Osher and Sethian [1988] S. Osher, J. A. Sethian,  Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations,  Journal of Computational Physics 79 (1988) 12–49.
  * Sethian [1996] J. A. Sethian,  A fast marching level set method for monotonically advancing fronts,  Proceedings of the National Academy of Sciences of the United States of America 93 (1996) 1591–1595.