Monday, November 11, 2013

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.