Laplace eq on a sphere in 3D

I use example from Chernyshenko-Olsh. paper for the 3D case:

$u=12(3x_1^2 x_2- x_3^3)/||x||^3)$

$f = -(72(x_2^3-x_2^2 x_3+x_2 x_3^2+x_3^3-x_1^2(5x_2+x_3)))/||x||^5) + u$

Details
5-point scheme of second order was used to integrate inside the tetrahedrals.

Due to memory/processing time I was not able to use $h< 0.05$

Result

Convergence plots

band of width 3h

2D Laplace equation on $y=sqrt(x)$, Neumann BC

Example 2:

Let $\Gamma$ denote curve $y=\sqrt{x}$ in 2D, $s$ be the arclength of  $\Gamma$ starting at the origin, and $\phi(x,y)$ denote signed distance function from $\Gamma$.

I look at problem of solving Laplace-Beltrami equation on $\Gamma$ for s = 0 to 1, for the known solution function $u=cos(4 \pi s)$. At the two endpoints, Neumann boundary condition is used.

As in previous work, signed distance function is used to define a narrow band extension of width $kh$ around $\Gamma$ and an unfitted finite element scheme[Deckelnik/Eliott] with element size h is applied to solve extended problem [Chern./Olsh.].

Arclength, distance, closest point transform, hessian
$\phi(x,y)$ is computed up to error of order h^2 using Matlab distance transform as follows: binary image of size sqrt(2N) by 2N is generated, for N > h^-2. Pixels at (x+1, Round(sqrt(x))+1) are set as foreground and the distance image is obtained using matlab bwdist function, which uses algorithm by Maurer. However this method becomes difficult to use for h < 1e-2 since Matlab runs out of memory.

For the surface y=sqrt(x) distance and closest point can also be found as roots of polynomial $2y^3 + (-2x_0+1)y -y_0=0$ where $x_0$ and $y_0$ are given point and $y$ gives the closest point on the curve.

From the values of distance function, Hessian is approximated at each point of the grid using centered differences. Cubic interpolation is used to get values at mesh points.

Arclength is approximated to h^2 using trapezoid method.

Unfitted mesh problem with Neumann boundary
The Deckelnik/Eliott unfitted mesh method needed some extra effort to be applied to this problem because of neumann boundary at the 2 sides of curve $\Gamma$.

Convergence

Clearly I am doing something wrong since the error does not follow expected convergence rates.

Narrow band with unfitted element scheme.

The width of band depends on mesh size h. I got following convergence plots:

To try strip width less than 3h, I would need to add some extra code to handle triangles that get cut off twice.

Example 5.1 from Eliott's paper:

$\Delta u = f$ on [-2,2]x[-2,2] \ Unit circle
$u = g_D$ on the sides of the square
${\partial u \over \partial n} = g_N$ on unit circle boundary
with
$f = -r^2 \\
u = (r^2+r^{-2})cos(2 \theta) + {r^4 \over 12}(sin^4(\theta)+cos^4(\theta)) \\
g_N = -(x^4+y^4)/3$


I get O(h^2) convergence in L2 norm:

 

Next: Apply the now working unfitted fem to narrow band

Found an error: My code occasionally subdivided quadrilaterals incorrectly.

Blue: interior
Green: boundary triangles
Red: part of boundary triangles that is actually used

After fixing above problem, I get following convergence plot.

I have also made the grid symmetric wrt axes and removed the code to zero-mean the solution.

Error plots
h = 0.05

 

0.0125

Deckelnik/Eliott method implementation (2)

I was able to significantly improve and speed up my code, which was needed to run more tests with finer meshes:

- Added ILU preconditioner and changed iterative solver from GMRES to CGS (conjugate gradient squared)

- Switched mesh generation to 'dumb' gridding with squares or hexagons. Since we don't actually need a mesh for any specific domain, I generate a very simple grid that contains required region, and restrict to that grid. Another benefit is that all triangles come out exactly same, so I simplified a number of matrix assembly steps by computing stiffness and mass sub-matrices once.

- Profiled my code and removed several bottlenecks. Most important of those was fixing the ways I handled the sparse matrix

This allowed me to run a couple more refinements for last week's results. Unfortunately it looks like I was wrong, I do not have O(h^2) convergence:

Update: I found an error in my reference triangle calculations (item 2 in list above).
Update 2: Above error with indexing only affected square meshes. Equilateral triangles in hex meshes are not affected by this error.

I have also checked again if my overall FEM implementation was ok, using a simple example u=cos(pi x)cos(pi y) on a [-1 1]x[-1 1] square. It was correct:

 

 So, I am still unclear why I did not get the Deckelnik unfitted method to work.