Immersed hybrid difference methods for elliptic boundary value problems by artificial interface conditions

1.   Introduction

In this paper we propose the immersed hybrid difference (IHD) methods for elliptic boundary value problems. Let us consider an elliptic boundary value problem on a multiply connected domain such that

Referring to Fig. 1 let ${ {\Omega} } = { {\Omega}^-}\cup { {\Omega}^+}$, ${\Gamma} = {\partial} {\Omega}^+ \cap {\partial} {\Omega}^-$ and ${\Gamma}_0 = {\partial} {\Omega}$. The vector ${\boldsymbol{\nu}}^\pm$ denotes the unit outward normal vector on ${\Gamma}$ from ${\Omega}^\pm$, respectively. We assume the source $f$ is smooth enough in ${\Omega}^+$ and the inner boundary ${\Gamma}$ is smooth for the sake of concise regularity analyses. The case $(\lambda, \mu) = (1, 0)$ corresponds to the Dirichlet condition and $(\lambda, \mu) = (0, 1)$ corresponds to the mixed boundary value problem with the Neumann interior boundary condition. It will be called the Neumann condition from here on since we are interested in treating the interior boundary condition. When ${\lambda} \mu \neq 0$ it is called the Robin condition. The boundary condition on the outer boundary is assumed to be the homogeneous Dirichlet condition for simplicity. For exterior domain problems the outer boundary condition can be replaced by some absorbing boundary condition, which has many physical applications such as the heat conduction and wave propagation problems.

The finite element and finite difference methods work very well for the boundary value problems if a mesh is constructed properly. Many researchers have studied and developed numerical methods for various problem by using finite difference (FD) and finite element (FE) approaches, to name a few, the (compact) finite difference [20,30], conforming and non-conforming FE (see [4] and references therein), mixed FE [1,6], discontinuous Galerkin (DG) [2,7,10] and hybridized DG (HDG) [9]. For these methods, boundary fitting mesh generation is the most effort taking part in implementation when the boundary of a domain is geometrically complicated. For curved interfaces high order methods fail to achieve the optimal order of convergence since the boundary conditions are imposed on a non-physical boundary with the usual elements. Thus, we need additional treatments (e.g. isoparametric curvilinear elements [3,5,19]) to recover the optimal order of convergence. On the other hand, various immersed methods for elliptic boundary value problems and relate problems have been introduced and developed to handle the computational domain involving complex geometries or curved interfaces/boundaries. For example, the boundary integral method (BIM) was introduced by Mayo [24,25]. Recent work by Marques et al. [23] and the kernel-free boundary integral method [32,34,33] are inspired by the BIM. Some of other well established methods are the immersed boundary method [27], the immersed interface method [21,22], the explicit-jump immersed interface method [31], the ghost fluid method [11,12], the immersed boundary smooth extension method [29], the matched interface and boundary method, and other methods [8,28,13]. These methods commonly reduce the efforts of mesh generation for the interface problem, because robust and efficient structured grids are used instead of a fitted mesh generation. The differences are how they enforce the boundary conditions and how they produce the solution in the domain.

The immersed hybrid difference method was developed by the first author [14] for the elliptic interface problems, and this idea is extended to solve boundary value problems. In the immersed boundary approach, instead of solving the equation on a boundary fitting mesh, we consider the extended problem of it on the whole domain by treating the inner boundary as an immersed interface. To do that, we introduced artificial interface conditions maintaining ellipticity of the weak formulation and imposing the given boundary condition simultaneously. Then, the immersed interface type numerical method solves the problem on the uniform rectangular mesh. Similar to the existing immersed methods, our approach reduces mesh generation efforts compared to the standard finite element or finite difference approaches. Also, there is no numerical integration to solve the problem and the number of global degrees of freedom can be reduced via embedded static condensation since our approach is based on the hybrid difference method. On the other hand, there are drawbacks in our approach compared to the existing immersed methods. The condition number of the local system can be large when the penalty parameter is very small, and the extended problem requires extra computational cost. If the size of ${\Omega}^-$ is small compared to that of ${\Omega}^+$ the additional computational cost can be ignored, and this kind of situation arises commonly in physical applications such as wave scattering problems.

The hybrid difference (HD) method is a finite difference version of the hybridized discontinuous Galerkin method and it was introduced by the authors and their colleagues for the elliptic, Stokes and Navier-Stokes equations [15,16,17,18]. Recently, the immersed boundary approach was applied to the HD method to develop the IHD method by the first author [14] for elliptic interface problems. The key feature of the IHD method is the introduction of the virtual to real transformation (VR-T) to find the finite difference stencil on interface cells, which refer to the cells that contain a portion of the given interface. The VR-T defines a relation between two locally extended sub-solutions on interface cells. In this paper we treat boundary ${\Gamma}$ as an immersed interface, therefore, we need a mesh generation of the whole domain ${\Omega}^+ \cup {\Omega}^-$, and the boundary conditions are imposed weakly by artificial interface conditions which contain a penalty parameter. The technique reminds us the Nitsche trick somehow.

The immersed hybrid difference method is consisting of two kinds of finite differences. Let us consider a decomposition of the domain into rectangular cells, ${\mathcal{T}}_h$ so that ${\Omega} = \cup_{R \in {\mathcal{T}}_h} R$. The cell difference approximates the local PDE on $R \in {\mathcal{T}}_h$,

The intercell difference patches the local solutions together on intercell edges by approximating the flux continuity,

The continuity of $u$ is guaranteed by assuming $u$ is single valued on cell edges. Here, ${\boldsymbol{\nu}}$ and ${\boldsymbol{\nu}}'$ are the outward unit normal vectors from $R$ and $R'$, respectively. In the IHD the boundary conditions on ${\Gamma}_0$ appears explicitly, however, the boundary condition on ${\Gamma}$ is included in the VR-T and it does appear explicitly.

The paper is organized as follows. In §2 the artificial interface conditions are derived, which will replace the boundary conditions on ${\Gamma}$ in (1.1), and some related analyses are provided. In §3 one dimensional immersed hybrid difference is introduced and the artificial interface conditions in §2 is used in deriving the virtual to real transformation. The two dimensional extension of the immersed hybrid difference are discussed in §4. In §5, numerical examples are provided to justify the theory and to demonstrate efficiency of the proposed method. Lastly, some concluding remarks are given at the end of the paper.

2.   The immersed boundary value problem by artificial interface conditions

In this section we extend the given boundary value problem (1.1) in a multiply connected domain to the immersed boundary value problem in a simply connected domain by introducing artificial interface conditions. For simplicity of discussion on solution regularity we assume that the interior boundary is smooth enough.

2.1.   Dirichlet condition

Let us consider the Dirichlet condition.

By solving another Dirichlet condition in ${\Omega}^-$

one obtains a solution

By the theory of elliptic operators in [26] we have the following lemmas.

Lemma 2.1. Suppose ${\Omega}^- \in C^{r+1, 1}$ for $r \ge 0$. Then,

with $-r-1 \le s \le r+1$.

Lemma 2.2. Suppose ${\Omega}^+$ is Lipschitz. Then,

and

We convert the above decoupled problems (2.1) and (2.2) into an approximating coupled problem with suitable interface conditions, and it will be called the immersed boundary value problem from here on.

Let

be the solution of the coupled problem;

with the exterior boundary condition $u^+ = 0$ on ${\Gamma}_0$, and with the interface conditions

Here, ${\boldsymbol{\nu}}^+$ and ${\boldsymbol{\nu}}^-$ are the unit outward normal vectors on ${\Gamma}$ from ${\Omega}^+$ and ${\Omega}^-$, respectively, and $0 < \epsilon < 1$.

Let us introduce function spaces for the extended solution $u = u^+ \cup u^-$.

The extended solution $u \in H^1_0 ({\Omega}^- \cup {\Omega}^+)$ of (2.1) and (2.2) satisfies the variational form, for $v \in H^1_0({\Omega}^- \cup {\Omega}^+)$,

Moreover, the solution $u_{{\epsilon}} \in H^1_0 ({\Omega}^- \cup {\Omega}^+)$ of the immersed boundary value problem (2.3) and (2.4) satisfies the variational form, for $v \in H^1_0({\Omega}^- \cup {\Omega}^+)$,

Subtracting (2.5) from (2.6) we have

Taking $v = u-u_{{\epsilon}}$ one obtains that

2.2.   Neumann condition

Let us consider the the problem with the Neumann boundary condition only on the interior boundary.

The extension of $u^+$ to the whole domain is made by considering an additional Dirichlet condition.

Then, the total solution $u \in H^1_0({\Omega}^- \cup {\Omega}^+)$ satisfies the following variational problem, for $v \in H^1_0({\Omega}^- \cup {\Omega}^+)$,

The immersed boundary approximation of the problems (2.8) and (2.9) is obtained by introducing appropriate interface conditions as follows.

with the interface conditions

Then, $u_{{\epsilon}} \in H^1_0({\Omega}^- \cup {\Omega}^+)$ satisfies the variational problem, for $v \in H^1_0({\Omega}^- \cup {\Omega}^+)$,

Subtracting (2.13) from (2.10) one obtains that

Take $v = u - u_{{\epsilon}}$, and using the duality of the Sobolev space, Lemma 2.1 and the trace theorem we have

Notice that $u^-$ is independent of ${{\epsilon}}$. The Poincaré inequality yields that

2.3.   Robin condition

Let us consider the the problem with the Robin boundary condition on the interior boundary (we consider only the case $\mu {\lambda} \ge 0$ and $\mu \neq 0$ in (1.1b)). Without loss of generality assume $\mu = 1$ and ${\lambda} \ge 0$.

As in the Neumann case the extension of $u^+$ to the whole domain is made by considering an additional Dirichlet condition.

Then, $u \in H^1_0({\Omega}^- \cup {\Omega}^+)$ satisfies the following variational problem, for $v \in H^1_0({\Omega}^- \cup {\Omega}^+)$,

The immersed boundary approximation of the problems (2.16) and (2.17) is obtained by introducing appropriate interface conditions as follows.

with the interface conditions

Then, $u_{{\epsilon}} \in H^1_0({\Omega}^- \cup {\Omega}^+)$ satisfies the variational problem, for $v \in H^1_0({\Omega}^- \cup {\Omega}^+)$,

Subtracting (2.19) from (2.18) one obtains that

By a similar argument as in the Neumann condition

3.   One dimensional immersed hybrid difference method

In this section we introduce the immersed hybrid difference method for the immersed boundary value problems.

3.1.   The 2nd order method for the Dirichlet condition

We begin with the simplest one dimensional method. Consider the domain $[0, 1]$ with two interface points, ${\Gamma}_1 \in I_{k-1}$ and ${\Gamma}_2 \in I_l$, where it is divided into $N$-cells so that $[0, 1] = \cup_{j = 0}^{N-1} I_j$ with $I_j = [\eta_{2j}, \eta_{2j+2}]$. The non-interface cell contains one cell-point and two intercell points, and the interface cells $I_{k-1}$ and $I_l$ contain an additional point that is the interface point as in Fig. 2.

We consider the one dimensional elliptic Dirichlet condition on two disjoint intervals ${\Omega}^+ = (0, {\Gamma}_1) \cup ({\Gamma}_2, 1)$: find $u^+$ such that

In addition to the above problem we consider an auxiliary problem such that

From here on we call $u = u^+ \cup u^-$ the exact solution.

The reason to consider the auxiliary problem is that one wishes to solve the problem (3.1) on a non-matching grid with an extra cost of solving the problem in (3.2). The numerical method on non-matching grids will be called the immersed (interface) hybrid difference method. By following the argument in the previous section we are going to solve the following coupled problem instead of (3.1) and (3.2). Let $u_{{\epsilon}} = u^+_{{\epsilon}} \cup u^- _{{\epsilon}}$ be the solution of the coupled problem;

with the exterior boundary condition $u^+ (0) = u^+(1) = 0$, and the interface conditions

Then the equations (3.3) with the interface conditions (3.4) can be solved by an immersed (interface) numerical methods. The advantage of this immersed approach will become more clearer for higher dimensional problems since the problem can be solved on a uniform mesh, which is independent of the geometric shape of an interface. In our approach the interior boundaries ${\Gamma}_1$ and ${\Gamma}_2$ are treated as immersed interface.

For the sake of simplicity we use $u$ for $u_{{\epsilon}}$ from here on as long as there is no danger of being confused. The total exact solution $u = u^+ \cup u^-$ satisfies

We are going to use the finite difference method. Then, the approximation

with $h_{k-1} = \eta_{2k}- \eta_{2k-2}$ is not appropriate for approximation of $u_{xx}$ on the interface cell since $u$ is not smooth enough on the interface cell $I_{k-1} = [\eta_{2k-2}, \eta_{2k}]$. It is inevitable to consider extensions of $u^{\pm}$ on interface cells. The interface cell $I_{k-1}$ is meant by the cell with which $I_{k-1} \cap {\Omega}^+ \neq \emptyset$ and $I_{k-1} \cap {\Omega}^- \neq \emptyset$.

Now, we consider an extension $\widetilde{ u}$ of $u$ on interface cells.

with reference to Fig. 2. Here, $\widetilde{ u}^+ = u^+$ on ${\Omega}^+$ and $\widetilde{ u}^- = u^-$ on ${\Omega}^-$, and $\widetilde{ u}$ is multiply valued on the interface cells $I_{k-1} = [\eta_{2k-2}, \eta_{2k}]$ and $I_l = [\eta_{2l}, \eta_{2l+2}]$. Then, the finite differences based on $\widetilde{ u}^{\pm}$ yield good approximations for derivatives as follows; on $I_{k-1}$

on the interface cell $I_{k-1}$. Here, the unit normal vector ${\boldsymbol{\nu}}^+$ on the cell interface represents the normal vector from the neighboring right cell and ${\boldsymbol{\nu}}^-$ represents that from the left cell. Notice that we use the notations, ${\partial}_{{\boldsymbol{\nu}}^-}$ and ${\partial}_{{\boldsymbol{\nu}}^+}$ for the normal derivatives on the genuine interface.

Now, we consider how to obtain $\widetilde{ u}$ from $u$, that is, how to obtain the extended values $\{ \widetilde{ u}^-_{2k-2}, \widetilde{ u}^+_{2k-1}, \widetilde{ u}^+_{2k} \}$ from the existing values $\{ \widetilde{ u}^+_{2k-2}, \widetilde{ u}^-_{2k-1}, \widetilde{ u}^-_{2k} \}$ (see Fig. 3). From here on, $\{ \widetilde{ u}^-_{2k-2}, \widetilde{ u}^+_{2k-1}, \widetilde{ u}^+_{2k} \}$ are called the virtual values, and $\{ \widetilde{ u}^+_{2k-2}, \widetilde{ u}^-_{2k-1}, \widetilde{ u}^-_{2k} \}$ are called the real values. It is time to introduce a relation between the real and virtual values, which is called the virtual to real transformation (the VR-T). To obtain the VR-T we use the artificial interface conditions in (3.4) and the governing differential equation

To impose the conditions in (3.6) the Lagrange interpolations, $\widetilde{ u}^- = \widetilde{ u}^-_{2k-2} \phi_1 + \widetilde{ u}^-_{2k-1} \phi_2 + \widetilde{ u}^-_{2k} \phi_3$ and $\widetilde{ u}^+ = \widetilde{ u}^+_{2k-2} \phi_1 + \widetilde{ u}^+_{2k-1} \phi_2 + \widetilde{ u}^+_{2k} \phi_3$ are used, where $\phi_1, \phi_2, \phi_3 \in P_2 (I_k) = span \{ 1, x, x^2 \}$ are the Lagrange basis. In the VR-T the first two conditions corresponds to the $1$-d version of the interface conditions in (2.4), and the third condition is called the consistency condition, which is derived from the governing differential equation. In view of Fig. 3 the VR-T yields the following linear relation between the virtual values and the real values

Then, the virtual values in the finite differences (3.5) can be rewritten in terms of all real values via the VR-T. For example, let us look at $-D_{xx}^h \widetilde{ u}^-_{2k-1}$, ${\partial}_{{\boldsymbol{\nu}}^+}^h \widetilde{ u}_{2k-2}^+$ and ${\partial}_{{\boldsymbol{\nu}}^-}^h \widetilde{ u}_{2k}^-$ in (3.5a) and (3.5b). The finite differences in real values are obtained as follows.

where $m_{ij}$ and $l_{ij}$ are the $(i, j)$ components of matrices $\mathcal{M}$ and $\mathcal{L}$, respectively.

Theorem 3.1. Let ${\Gamma}_1 \in [\eta_{2k-2}, \eta_{2k}]$ and $\delta = 2 \frac{{\Gamma}_1 - \eta_{2k-1}}{h_{k-1}}$. The matrix $\mathcal M$ for the Dirichlet condition is nonsingular if $|\delta^3- \delta | \neq \frac{4 {{\epsilon}}}{h_{k-1}}$, where ${{\epsilon}}$ is the penalty parameter, $0 < {{\epsilon}} < 1$.

Proof. The matrix $\mathcal M$ is derived by solving the homogeneous jump and consistency conditions, that are,

Consider the function $d(x) = u^+(x) - u^- (x)$. Using $d \in P_2(I_{k-1})$, the conditions in (3.7) yield that $d_{xx} = 0$, $d_x = K$ and $d = K(x-{\Gamma}_1)$ for a constant $K = -\frac{u^+({\Gamma}_1)}{{{\epsilon}}}$ on $I_{k-1}$.

To show that $\mathcal M$ is an injection let us assume that the virtual values are zero so that $u^-_{2k-2} = u^+_{2k-1} = u^+_{2k} = 0$ (that is, $-1 < \delta \le 0$ with ${\Gamma}_1 = \eta_{2k-1} + \delta \frac{h_{k-1}}{2}$ as in Fig. 3). Then, $d = K(x-{\Gamma}_1)$ yields that

The Taylor expansion yields that

where $\delta = \frac{{\Gamma}_1-\eta_{2k-1}}{h_{k-1}/2}$. The condition, $K = -\frac{u^+({\Gamma}_1)}{{{\epsilon}}}$ yields that $K( h_{k-1}(\delta-\delta^3) + 4 {{\epsilon}}) = 0$. If $\delta-\delta^3 \neq - \frac{ 4 {{\epsilon}}}{h_{k-1}}$, $K = 0$. Then, the real values satisfy $u^+_{2k-2} = u^-_{2k-1} = u^-_{2k} = 0$, and $\mathcal M$ is injective.

Now, let us consider the other case so that the virtual values are $\{ u^-_{2k-2}, u^-_{2k-1}, u^+_{2k} \}$ (that is, $0 \le \delta < 1$). Let $u^-_{2k-2} = u^-_{2k-1} = u^+_{2k} = 0$. Then, $d = K(x-{\Gamma}_1)$ yields that

The Taylor expansion yields that

where $\delta = \frac{{\Gamma}_1-\eta_{2k-1}}{h_{k-1}/2}$. The condition, $K = -\frac{u^-({\Gamma}_1)}{{{\epsilon}}}$ yields that $K( - h_{k-1}(\delta-\delta^3) + 4 {{\epsilon}}) = 0$. Then, $\mathcal M$ is injective under the assumption.

Theorem 3.2. Let ${\Gamma}_1 \in [\eta_{2k-2}, \eta_{2k}]$ and $\delta = 2 \frac{{\Gamma}_1 - \eta_{2k-1}}{h_{k-1}}$. The matrix $\mathcal M$ for the Neumann condition is nonsingular with the penalty parameter $0 < {{\epsilon}} < 1$.

Proof. The matrix $\mathcal M$ is obtained by solving

Let us consider $d(x) = u^+(x) - {{\epsilon}} u^- (x)$. Using $d \in P_2(I_{k-1})$, the conditions in (3.4) yield that $d_{xx} = 0$, $d_x = 0$, hence $u^+ = {{\epsilon}} u^- + K$ for some constant $K$. Now, using $(u^+ - u^- )({\Gamma}_1) = 0$ we have $u^+({\Gamma}_1) = u^-({\Gamma}_1) = \frac{K}{1-{{\epsilon}}}$.

To show that $\mathcal M$ is an injection let us assume that the virtual values are zero so that $u^-_{2k-2} = u^+_{2k-1} = u^+_{2k} = 0$ (that is, $-1 < \delta \le 0$ with ${\Gamma}_1 = \eta_{2k-1} + \delta \frac{h_{k-1}}{2}$ as in Fig. 3). Then,

The Taylor expansion yields that

Using $u^+({\Gamma}_1) = \frac{K}{1-{{\epsilon}}}$, we have $K(\delta^2-\delta - \frac{2}{1-{{\epsilon}}} ) = 0$. Since $\delta^2-\delta - \frac{2}{1-{{\epsilon}}} \neq 0$, $K = 0$, which implies $u^+_{2k-2} = u^-_{2k-1} = u^-_{2k} = 0$. The injectiveness of $\mathcal M$ is proved. The proof for the remaining case, $0 \le \delta < 1$, can be obtained by a similar way.

Remark 1. In our approach we use the penalty parameter ${{\epsilon}} = h^\alpha$ with $\alpha > 1$. Therefore, the case $|\delta^3- \delta| = 4 h^{\alpha-1}$ (for Theorem 3.1) can be avoided if the genuine interface ${\Gamma}_1$ is not located somewhere very close to the cell center or the cell interfaces.

Now, let us introduce the hybrid difference method. Let $U$ be the approximate solution of $\widetilde{ u}_{{\epsilon}}$ and $U^+$ and $U^-$ are those of $\widetilde{ u}^+_{{\epsilon}}$ and $\widetilde{ u}^-_{{\epsilon}}$ by the immersed hybrid difference method. Then the hybrid difference method can written as follows. The PDE (3.1) at the cell interior point and the flux continuity at the intercell points derives the following hybrid difference equations:

and on the interface cell $I_{k-1}$ and $I_l$ (see that $\eta_{2k-1} \in {\Omega}^-$ and $\eta_{2l+1} \in {\Omega}^+$ in Fig. 2),

Note that the intercell equations related to (3.11) are contained in (3.9b) and (3.10b). The VR-T does not appear explicitly in the hybrid difference formulation, and it is contained implicitly when evaluating the finite differences in the interface cells.

3.2.   A higher order method

Now we consider a higher order method based on the cell configuration in Fig. 4. It should be noted that the interior cell points $\{ \eta_{4k+1}, \eta_{4k+2}\}$ of the cell $[\eta_{4k}, \eta_{4k+3}]$ are the Gaussian points.

In this case we have the four real and four virtual values in the interface cell. The real values are $\{ U^-_{4k}, U^-_{4k+1}, U^+_{4k+2}, U^+_{4k+3} \}$ and the virtual values are $\{ U^+_{4k}, U^+_{4k+1}, U^-_{4k+2}, U^-_{4k+3} \}$. The VR-T is derived by solving the following four independent equations by considering $U^-, U^+ \in P_3(I_k) = span\{1, x, x^2, x^3 \}$.

For the Neumann condition we use the following VR-T,

For the Neumann VR-T we can use $-U_{xx}^+ (\eta) = -U_{xx}^-(\eta) +f$ instead of $-U_{xx}^+ (\eta) = -{{\epsilon}} U_{xx}^-(\eta) +f$ since $-D_{xx} u^-_{{{\epsilon}}} = 0$, and the former makes the convergence behavior more regular according to our numerical test in Section 5. For the Robin condition only the equation (3.12b) will be changed to

4.   Two dimensional IHD method

In this section the one dimensional immersed method is revised to accommodate the two dimensional case. Since the immersed hybrid difference method for the non-interface cell is trivial, we mainly discuss on the IHD on interface cells and related two dimensional VR transformation.

4.1.   Hybrid difference method

The hybrid difference method for the equation (1.1) is composed of two kinds of finite differences, the cell finite difference and the intercell finite difference as follows. With reference to the cell configuration in Fig. 5 the cell finite difference is

on the cell $R_1$, and the interface finite differences are

In the above finite differences $U^+(\eta_1)$ and $U^+(\eta_3)$ are virtual values, and they should be replaced by the real values $\{ U^-(\eta_1), U^-(\eta_3), U^+(\eta_4), U^+(\eta_5), U^+(\eta_8) \}$ via the VR transformation.

4.2.   Two dimensional VR transformations

In Fig. 6, two grid types are displayed. The left one is the $Q_2^*$-grid cell, the others are the $Q_3^*$-grid cells. Note that the interior grid points are the Gauss points on the corresponding dashed line. If an interface passed through a cell, this cell is called the interface cell. The interface cells are divided into two types. It is called the coherent interface cell if the interface and the lines of integration (the red lines) intersect, and the incoherent interface cell has no intersection. In Fig. 6, the left and center ones are coherent and the right one is incoherent. In the IHD method, an incoherent cell is regarded as a regular cell which an interface does not pass through at all. Now we introduce two dimensional VR transformations for interface coherent cells. Firstly, we consider the constitutive equations on the $Q_2^*$-coherent grid cell. Since the $Q_2^*$-grid cell contains five real and virtual values, respectively, we need the constitutive equations composed of five equations as follows.

The VR-T for the Dirichlet condition (2.3):

The VR-T for the Neumann condition (2.11):

These are simple and natural generalizations of the one dimensional VR transformations.

Since there are twelve real and virtual variables, respectively, (see, the center in Fig. 6) in the coherent $Q_3^*$-cell grid, we need a set of $12$-constitutive equations as follows.

Here, $p = \tau_1, \tau_2, \tau_3, \tau_4$ and $q = \eta_4, \eta_5, \eta_8, \eta_9$. The conditions (4.5) and (4.6) are used to derive the VR transformations for numerical computation in Section 5.

In order to define the VR transformation, the number of equations required is $((k+1)^2-4)$ which is equal to the number of grid points in a $Q_k^*$-grid cell. It is natural to consider for the $4(k-1)$ distinct interface points and the $(k-1)^2$ interior Gaussian points of a cell.

4.3.   The binomial space $Q_k^*(R)$

In order to impose the constitutive equations, we need binomial spaces for $U^-$ and $U^+$ satisfying the following properties:

(P1): The degree of freedom of $Q^*_k (R) = (k+1)^2-4$.

(P2): $P_k(R) \subset Q_k^*(R) \subset Q_k(R)$

(P3): $Q_{k-2} (R) \subset \Delta Q_k^*(R)$.

Here, $Q_k(R)$ is the space of binomials of degree $\le k$ for each variable, and $P_k(R)$ is that of the total degree $\le k$. In order to obtain the optimal convergence rate, the condition (P2) is required. The condition (P3) guarantees independent consistency equations. Since $U$ and ${\partial}_{{\boldsymbol{\nu}}} U$ become functions which have the number of degrees of freedom $2k$ and $(2k-1)$, respectively, in the corresponding variable on a slanted interface, the $(2k-2)$-collocations on $U$ and ${\partial}_{{\boldsymbol{\nu}}} U$ on the interface are independent. If the interface is parallel to the coordinate axes then both $U$ and ${\partial}_{{\boldsymbol{\nu}}} U$ degenerate to polynomials of degree $k$. Since $2k-2 \le k+1$ for $k\le 3$, they are independent only when $k \le 3$.

The following two binomial spaces are likely to be most commonly used,

● $Q_2^*(R) = span \{ 1, x, y, x^2, y^2 \}$,

● $Q_3^* (R) = span \{ 1, x, y, x^2, xy, y^2, x^3, x^2y, xy^2, y^3, x^3 y, xy^3 \}$.

It is easy to see that $Q_3^*$ satisfy the properties (P1)–(P3) and $Q_2^*$ satisfies $P_1(R) \subset Q_2^*(R)$ instead of (P2).

Remark 2. Let the extended solutions satisfy $\widetilde{ u}^\pm \in C^{k+1}(R)$ for an interface cell $R$ with a smooth enough interface. Then,

for $k \ge 3$, and

5.   Numerical experiments

In this section, we present numerical experiments in one and two dimensions to illustrate efficiency of the IHD method. We also performed numerical experiments by changing $\alpha$ which determines the penalty parameter $\epsilon$ (see, Remark 1).

5.1.   One-dimensional numerical examples

We perform numerical tests for a one-dimensional problem.

Ex. 1. Let us consider the differential equation, $-u_{xx} = f$ on ${\Omega}^+ = [0, \frac{4}{7}] \cup [\frac{5}{7}, 1]$, and the exact solution satisfies

Then, the extended solution satisfies $u^-(x) = (e^{\frac{4}{7}} -1) (5-7x)+\sin(\frac{5}{7} \pi)(7x-4)$ on $[\frac{4}{7}, \frac{5}{7}]$.

The immersed boundary value problem solves the extended problem,

with $u(0) = u(1) = 0$. The interface conditions are imposed as

for the Dirichlet condition, and the interface conditions for the Neumann condition are

The penalty parameter is taken to be ${{\epsilon}} = h^\alpha$ with $\alpha = 1, ..., 4$.

The computation is made on the cell configuration in Fig. 4. According to the standard theory of approximation the best possible order of convergence is of $O(h^4 + h^\alpha)$ when $u \in C^4({\Omega}^+)$ with ${{\epsilon}} = h^\alpha$. Here, $h$ represents the scale of a mesh. Fig. 7 shows that all the numerical results for the Dirichlet condition concord well with the expected convergence. For the Neumann problem convergence behaves a bit erratically, and convergence stays around $O(h^3)$ while the expected convergence is of $O(h^4)$ when $\alpha = 4$ (see Fig. 8). It is interesting to see the scaled condition number ($(\text{condition number})\times h^3$) stays constant with respect to the penalty parameter ${{\epsilon}} = h^\alpha$. It is clear that the condition number for the VR-T heavily dependent on the penalty parameter ${{\epsilon}}$, however, this ill-conditioning in the VR-T does not affect on the condition number of the global stiffness system.

5.2.   Two-dimensional numerical examples

We provide two-dimensional numerical examples for the Poisson equation on three multiply connected domains with a curved inner boundary. Numerical experiments were performed only on the $Q_3^*$-grids with the uniform mesh grids. The errors are measured in the discrete $L_2$-norm,

for a mesh ${\mathcal{T}}_h$ of ${\Omega} = {\Omega}^+ \cup {\Omega}^-$. Fig. 9 represents three different immersed boundaries and corresponding mesh grids that are used in our computation. Since the inner boundary is treated as an interface in our approach we call it frequently as interface. A mesh grid is called interface coherent if all interface cells are coherent and it is called interface incoherent if there is at least one incoherent interface cell. The incoherent interface cells are treated as a regular cell in our computation.

In the immersed methods one wishes to use a fixed mesh (usually, the uniform mesh), independently of the shape of interface or immersed boundary. Then, it is inevitable to have incoherent interface cells for curved interfaces. Let us introduce a measure that indicates the portion of coherent interface cells among all interface cells, which is called the coherence ratio:

In the Fig. 9 the immersed boundaries satisfy

respectively. The the five-leafed flower is expressed in the polar coordinates.

Fig. 10 represents the coherence ratio for each immersed boundary for the mesh configurations in Fig. 9. The rotated ellipse has the least coherence ratio among them. The five-leafed flower has the less oscillatory coherence ratio.

Ex. 2. We consider the Poisson equation on the multiply connected domains as shown in Fig. 9.

where the exact solution is given as $u(x, y) = e^x \cos(y)+x^2$ on the square $[-1, 1]\times[-1, 1]$. The Dirichlet and Neumann condition are considered. As shown in Section 2 we solve the extended problem

with the interface conditions

with $g_D = u|_{\Gamma}$ and $g_N = {{\partial}_{{\boldsymbol{\nu}}^+}} u|_{\Gamma}$ by varying ${{\epsilon}} = h^\alpha$, $\alpha = 1, ..., 4$. The computation is performed on the (quasi-) uniform $Q_3^*$-mesh grid.

Fig. 11 represents numerical results for the Dirichlet problem on the domain $[-1, 1]\times [-1, 1] \setminus D$, where $D$ is a disk with radius $0.3$, and Fig. 12 represents that for the Neumann condition on the same domain. Both show the optimally expected convergence, $O(h^\alpha + h^4)$. Unlikely to the one dimensional case the Neumann condition also achieves the optimal convergence.

Fig. 13 represents the convergence history for the Dirichlet problem on the domain with the rotated-ellipse inner boundary. Even though the coherence rate is relatively lower than those of others, convergence did not seem to be influenced much by the lower coherence ratio. Fig. 14 represents the convergence history for the five-leafed flower. The boundary is non-convex and a poorer or oscillatory convergence is observed with the uniform mesh. Therefore, we adopt a quasi-uniform mesh, which is obtained by modifying the uniform mesh slightly on the places that the interface becomes horizontally or vertically flat. By doing that we can reduce the number of interface cells with reentrant interface boundary, otherwise cells with reentrant interface boundary are recognized as interface incoherent cells in our programing. The cell with reentrant interface is meant by that interface comes in and goes out of a cell through the same edge of a cell. For convex interfaces as in the first two figures in Fig. 9 the four interface collocation points as in the middle of Fig. 6 are used to imposed artificial interface conditions. For the five-leaped flower case we picked two sets of collocation points, each of which are consisting of four collocation points. One set is used to impose $[[ U]] = 0$, and the other set is used for the remaining interface condition. This seems to enhance the stability of the VR transformation. The convergences are optimal, however it is much lower than the theoretically expected one when $\alpha = 4$. The bad conditioning by a very small penalty parameter in the VR transformation and non-convexity of the domain seems to have some negative influence in convergence.

The condition number of the VR transformation increases in the order of $O(h^\alpha)$ as $h$ decreases, however the (scaled) condition number for the global stiffness system is independent of the penalty parameter $h^\alpha$ as mentioned in the one-dimensional case.

According to our numerical experiments non-convex interface boundaries induce some difficulties in mesh generation if one wish to reduce the number of cells with reentrant interface, which is an important part for a more regular convergence of numerical solutions. It is intuitively trivial to see that the higher degree methods yields higher coherence ratios, hence it is desired to use at least the $Q_3^*$-IHD method. Moreover, complexity of the boundary shapes seems to have more influence on convergence than the coherence ratio.

6.   Concluding remarks

In this paper the immersed hybrid difference method for elliptic boundary value problems on multiply connected domains is presented. The artificial interface conditions combined with the virtual to real transformation introduced in [14] enable us to solve the boundary value problems on general rectangular meshes which are constructed independently of the geometric shape of the boundaries. This approach dramatically reduces mesh generation efforts when a problem domain has a boundary of complicated shapes. Although we need additional degrees of freedom due to the approximation in $\Omega^-$, the additional computational cost can be negligible when the size of $\Omega^-$ is relatively small, compared with that of $\Omega^+$. The numerical results indicate that the IHD method can be promising in handling problems on domains with a complicatedly shaped boundary. The optimal rate of convergence can be achieved when the parameter $\alpha$ is chosen properly.

Acknowledgments

We would like to express sincere thanks to anonymous referees, whose invaluable comments improved the quality of the paper a lot.