Applying previous method to the Laplace equation on the unit circle with Monte-Carlo integration, I got following results:

h        Nodes      Time       |  L2(Gamma)            H1(Gamma)              Area*St.dev
---------------------------------------------------------------------------------------
0.1000     353     0:0:1s     |  1.71e-03                8.15e-02                   3.01e-03
0.0500     697     0:0:1s     |  3.14e-04   2.45      2.25e-02   1.86         5.58e-04        2.43
0.0250    1361    0:0:4s     |  3.39e-05   3.21      5.84e-03   1.95         9.86e-05        2.50
0.0125    2717    0:0:19s   |  3.29e-06   3.37      1.49e-03   1.97         1.66e-05        2.57
0.0063    5457    0:5:16s   |  4.14e-07   2.99      3.73e-04   2.00         3.07e-06        2.44
0.0031   10837   2:46:24s |  5.21e-08   2.99      9.38e-05   1.99         4.37e-07        2.81
-----------------------------------------------------------------------------------------

Another way to integrate irregular areas in unfitted triangles is to place ${h^-1}$ points on the domain boundary, triangulate and integrate.

h        Nodes       Time      |  L2(Gamma)            H1(Gamma)               Area*St.dev

---------------------------------------------------------------------------------------

0.1000     353     0:0:1s     |  1.43e-03                8.14e-02                   0.00e+00 

0.0500     697     0:0:4s     |  1.93e-04   2.89      2.25e-02   1.86         0.00e+00        NaN 

0.0250    1361    0:0:17s   |  2.55e-05   2.91      5.83e-03   1.94         0.00e+00        NaN 

0.0125    2717    0:1:5s     |  3.28e-06   2.96      1.49e-03   1.97         0.00e+00        NaN 

0.0063    5457    0:4:53s   |  4.14e-07   2.99      3.73e-04   2.00         0.00e+00        NaN 

0.0031   10837   0:26:53s |  5.21e-08   2.99      9.38e-05   1.99         0.00e+00        NaN 

-----------------------------------------------------------------------------------------

I am surprised to get very similar results for the Monte-Carlo and triangulation method. This seems a bit suspicious.

Instead of using uniform random points over entire triangle, integrate linear approximation by splitting into triangles and using Monte-Carlo integration in the thin rectangular strip shown in following figure:
(Click the figure to view)

To generate uniform random points in the thin strip:

1. For triangle K, take $x_1$, $x_2$ to be intersection of boundary $\delta \Omega$ and edges of K ($\phi(x_i) = 0$).

2. Approximate m = max($\phi(x)$) on the segment between $x_1$ and $x_2$. This gives the thin rectangle containing leftover region to be integrated.

3. Generate uniform random points in (0,1)x(0,1), scale rotate and translate them to the thin rectangle.

Results for 2D Laplace equation inside the unit circle:

Results:
h            Nodes     Time     |  L2                          H1                      MaxVariance
---------------------------------------------------------------------------------------
0.2500     261     0:0:0s     |  1.61e-02                3.00e-01             6.30e-03
0.1250     937     0:0:1s     |  1.56e-03   3.36      4.03e-02   2.90   8.59e-04    2.87
0.0625    3499    0:0:6s     |  8.20e-05   4.25      5.78e-03   2.80   1.14e-04    2.92
0.0312   13477   0:2:7s     |  6.49e-06   3.66      1.01e-03   2.51   2.32e-05    2.29
0.0156   52707   1:63:34s |  6.62e-07   3.29      1.53e-04   2.73   2.43e-06    3.26
-----------------------------------------------------------------------------------------

Convergence appears to be as expected or better. I used 6th order quadrature scheme with 9 points to integrate most of the stiffness matrix terms exactly, which probably is causing better than expected orders of convergence. Should change that to a 2nd order scheme.

P2 Finite elements with Monte-Carlo integration

Before trying to solve equation on surface, I try to make sure the method works with a simple 2D example,
\begin{eqnarray*}
& -\Delta u + u = f  \\
& u = \cos(\pi r^2)  \\
\end{eqnarray*}
inside the unit circle in 2D.

using P2 elements and Monte-Carlo integration. For P2 method, I need each integral in assembly to be approximated to order $h^2$. Therefore, I use $M := C_1h^{-4},\ C_1 = 1$ uniform random points in the triangle, generated once and reused in every triangle.

For each triangle, the Monte-Carlo variance is made sure to be less than $C_2h^4$. If any of the 36+6 integrals do not satisfy the variance criteria, M new points are generated and added to the integration. This causes most of the integrals per triangle to be very over-integrated.

Unfitted boundary triangles: 46 - 102 - 210 - 430 - 866
Monte-Carlo points per triangle: 256 - 4096 - 65536 - 1048576 - 16777216

h             Nodes Time        |  L2                          H1                            MaxVariance
---------------------------------------------------------------------------------------
0.2500     261    0:0:1s      |  1.56e-02                2.75e-01                   1.24e-02
0.1250     937    0:0:4s      |  1.32e-03   3.57      3.94e-02   2.80         3.12e-03   1.99
0.0625    3499   0:0:26s    |  1.87e-04   2.82      7.00e-03   2.50         7.78e-04   2.01
0.0312   13477  0:11:51s  |  8.71e-06   4.42      9.85e-04   2.83         1.93e-04   2.01
0.0156   52707 2:177:28s |  6.22e-07   3.81      1.68e-04   2.55         3.74e-05   2.37
-----------------------------------------------------------------------------------------

L_2 convergence order lower than optimal order of 3 are caused by the following issue:
$|E[f(x_i)] - I(f)| < 5\sqrt{Var[f(x_i)]/M}$ holds with probablility $1-10^{-5}$. So $C_2$ must be 1. Also, when the number of triangles increases, this probability may grow.

The better than optimal order of convergence in H1 norm is probably caused by over-integration.

Repeat of 2D Laplace experiment on circle with Lagrange P2 elements.

1. Fitted mesh

Done.
Since the mesh approximates the region with line segments, the L2 error can not go better than O(h^2).

2. Unfitted mesh, linear interface approximation on boundary triangles

Done.
Again, L2 error can not decrease faster than O(h^2) because of linear boundary approximation.

3. Unfitted mesh with Monte-Carlo integration on boundaries

Done (?)


4. Unfitted mesh with linear interface approximation and Monte-Carlo integration for leftover irregular regions.



5. Unfitted mesh, Moments fitting Gaussian integration on boundary triangles

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.