Laplace equation in 2D with Dirichlet BC

This example is exactly the same as the 2D Laplace-Beltrami equation on curve y=sqrt(x) 2 posts ago.

The only difference is instead of Neumann boundary condition at the two edges of the curve, I used Dirichlet BC. This is slightly cheating since I did not actually fit the mesh at the right side boundary, instead assigning correct values to mesh nodes for which curve arclength > 1. In figure below these nodes are shown in green.

Convergence:

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.

Deckelnik-Eliott unfitted scheme tests

 I tested the unfitted FE scheme from Deckelnik paper on a unit circle, solving 

$-\Delta u + u = f$
on unit circle,

& Neumann BC

for a known
$u = cos(\pi r^2)$

Here are my error convergence results:

It appears there is no problems with my implementation.

---

Unfitted finite elements approach from Deckelnick / Eliott / Dziuk paper


imajna.oxfordjournals.org/cgi/reprint/drn049v1.pdf

This paper describes almost the same approach that we are using - obtaining solution on surface by extending to unfitted mesh containing h-narrow band around the surface.

Only difference seems that they are using orthogonal projection from earlier Dziuk's paper to extend the domain, and they do mention Greer's idea as well.

In essence, the unfitted mesh is restricted to a fitted mesh. Elements that fall across the boundary level sets are chopped to only the part inside the band. This creates problems when the restricted triangle comes out very small or narrow.














In the paper, diagonal preconditioning is used. Here are my convergence results with no preconditioning:

       

L2 convergence of restrictions to bands around surface

1. We expected the L2 error on Surface to level off with increasing width of band after a certain value. The following plot confirms this:

(The series represent meshes of edge sizes 0.002, 0.005, 0.01 restricted to bands of increasing width around the surface)

2. Next I refine the mesh and keep same strip width from above runs.

 

Conclusion:

While it is true that sufficient strip width (0.5) gives O(1/nodes) convergence [same as O(1/h^2)], it seems discouraging as we have to take quite a number of triangles that are not crossed by the surface. Figure below shows that the jagged boundary causes the method to behave worse than the case with simple square extension.  am currently looking for ways to improve this, as well as searching how did other people handle such problems.

Related papers I am looking at:

http://www.sciencedirect.com/science/article/pii/0898122185900586#

The problem is to solve a 2D elliptic equation with moving interface using a fixed unfitted mesh. Vertical derivatives of solution across the interface are prescribed.

 

Their approach is also to split the interface gradient into 2 components, one of which is given by oblique condition and the other component is added to the equation matrix. The ellipticity of the problem is proved.

Feb 1 2013 (Friday) meeting

1. Look at paper Chernyshenko - Non-degenerate Eulerian FEM for PDEs on surfaces

• implement simple benchmark test cases for known problems
• extend to parabolic equations (?)
• Look at texture on surface as solution (steady equilibrium state) of parabolic PDE
• Mean curvature flow (Find more references about it!)
• Look at paper Greer-Improvement of Eulerian Method for PDE on general geometries

2. Look at vizualisation soft