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.