Start page of documentation

Download the code

Documentation: PDF, PS

Correction terms

The points where the 5-point stencil is cut by the interface are called irregular points. There, the approximation (2) has to be corrected by adding so called correction terms. They can be written for each of the appearing derivatives in (1) separately.

Correcting the second order differences

Consider a situation like in figure above and let $({x_\alpha}_1, {y_\alpha}_1)$ be the coordinates of the intersection point $\alpha_1$ and $({x_\alpha}_2, {y_\alpha}_2)$: coordinates of the intersection point $\alpha_2$. Let $h^+:=x_{i+1}-{x_\alpha}_1$, $h^-:=x_i-{x_\alpha}_1$, $k^+:=y_{j+1}-{y_\alpha}_2$, $k^-:=y_j-{y_\alpha}_2$. Then, the corrected differences of second order derivatives are

$$\partial_{xx}u(x_i,y_j) \approx \frac{1}{h_x^2}\left(u_{i+1,j}-2u_{i,j}+u_{i-1,j}\right) -\frac{1... ...partial_x u]_{\alpha_1} + \frac{1}{2}(h^+)^2 [\partial_{xx}u]_{\alpha_1}\right)$$ (3)

$$\partial_{xx}u(x_{i+1},y_j) \approx \frac{1}{h_x^2}\left(u_{i+2,j}-2u_{i+1,j}+u_{i,j}\right) +\frac{1... ...partial_x u]_{\alpha_1} + \frac{1}{2}(h^-)^2 [\partial_{xx}u]_{\alpha_1}\right)$$ (4)

for derivatives in $x$-direction and similarly for the $y$-derivatives

$$\partial_{yy}u(x_i,y_j) \approx \frac{1}{h_y^2}\left(u_{i,j+1}-2u_{i,j}+u_{i,j-1}\right) -\frac{1... ...partial_y u]_{\alpha_2} + \frac{1}{2}(k^+)^2 [\partial_{yy}u]_{\alpha_2}\right)$$ (5)

$$\partial_{yy}u(x_i,y_{j+1}) \approx \frac{1}{h_y^2}\left(u_{i,j+2}-2u_{i,j+1}+u_{i,j}\right) +\frac{1... ...partial_y u]_{\alpha_2} + \frac{1}{2}(k^-)^2 [\partial_{yy}u]_{\alpha_2}\right)$$ (6)

Note, we have introduced new variables into a system. At each intersection point we have 6 new unknowns: $[u]$, $[\partial_x u]$, $[\partial_y u]$, $[\partial_{xx}u]$, $[\partial_{yy}u]$ and $[\partial_{xy}]$, where in the case of Poisson equation some of them are redundant.

In general, jump in the mixed derivative is not used in expressions (3-6) and jumps in $x$-derivatives are not used at $y$-intersections, as well as jumps in $y$-derivatives are not used in $x$-intersections.

All jumps are stored in one supervector $J$ (see ASSEMBLE/run_ejiim.m, variable J) with ordering

$$J = ([u]_{\alpha_1}, [\partial_x u]_{\alpha_1}, [\partial_{y} u... ...al_{xy}u]_{\alpha_1}, [u]_{\alpha_2}, \ldots, [\partial_{xy}u]_{\alpha_l})$$

where $l$ is the number of intersection points $\alpha_1$, $\alpha_2$, ..., $\alpha_l$. Coefficients of $J$ in the correction terms (3-6) are stored in matrix $\boldsymbol{\Psi}$, in code denoted by P (see ASSEMBLE/run_ejiim.m).

With this, approximation of the differential operator $\Delta$ in domain $\Omega^*$ can be written as

$$\mathbf{A} U + \boldsymbol{\Psi} J = F$$ (7)

Details of constructing the P matrix

Matrix $\boldsymbol{\Psi}$, in code stored in variable P is constructed of two parts

$$\boldsymbol{\Psi} = \left( \boldsymbol{\Psi}_D , \boldsymbol{\Psi}_N \right) ,$$

where $\boldsymbol{\Psi}_D$ corresponds to intersection points belonging to $\partial\Omega_D$ and $\boldsymbol{\Psi}_N$ - to intersection points along $\partial\Omega_N$.

Each of $\boldsymbol{\Psi}_D$ and $\boldsymbol{\Psi}_D$ is constructed as sum

$$\boldsymbol{\Psi}_D = \boldsymbol{\Psi}^{xx}_D + \boldsymbol{\Psi}^{yy}_D$$

$$\boldsymbol{\Psi}_N = \boldsymbol{\Psi}^{xx}_N + \boldsymbol{\Psi}^{yy}_N$$

where $\boldsymbol{\Psi}^{xx}$ corrects the approximation of $\partial_{xx}$ derivative according to (3,4) and $\boldsymbol{\Psi}^{yy}$ corrects the approximation of $\partial_{yy}$ derivative according to (5,6). See ASSEMBLE/run_ejiim.m.

Note: Matrix $\Psi$ depends only on the differential operator and geometry. No dependence on boundary conditions!

Construction of the matrix $\boldsymbol{\Psi}$ is done in function DIFFOP/corrections.m. The interface information is stored in the structure variable IX, see Section 3 for details.

Constructing of $\boldsymbol{\Psi}^{xx}$ matrix

Idea: run a loop over all intersection points and check, which grid points are affected.

We know that $\partial_{xx}$ derivative is affected only by $x$-type intersections, so the matrix $\boldsymbol{\Psi}^{xx}$ has zeros in the right block, corresponding to $y$-type intersection points.

In a loop over all $x$-intersections:

1. Determine the coordinates of the anchor point, call it $(x_i, y_j)$.
2. Check if the interface lies left or right from the anchor point.
   • If interface lies LEFT from the anchor, use (4) at point $(x_i,y_j)$ and (3) at $(x_{i-1},y_j)$.
   • If interface lies RIGHT from the anchor, use (3) at the point $(x_i,y_j)$ and (4) at $(x_{i+1},y_j)$

Construction of $\boldsymbol{\Psi}^{yy}$ is done completely analogously, only now zeros are in the left block of $\boldsymbol{\Psi}^{yy}$ and loop has to be done over all intersections.

In the actual version of the code, the EJIIM system is solved by Matlab ` $\backslash$'-operator. For larger problems, especially three dimensional, the Schur-complement approach together with some fast solver for inverting the operator $\mathbf{A}$ is highly suggested [3,1,2].