Results for 2D Laplace Eqn on the circle, with H=0

h          Nodes   Time    |  L2(Gamma)     H1(Gamma)     Area*St.dev
---------------------------------------------------------------------------------------
0.1000   501    0:0:6s     |  4.19e-02   2.59e-01   9.99e-03
0.0500  1013    0:0:7s    |  5.01e-03   3.06 3.49e-02   2.89 2.10e-03   2.25
0.0250  1997    0:0:8s    |  2.02e-03   1.31 1.33e-02   1.40 3.41e-04   2.62
0.0125  3989   0:0:18s   |  9.38e-04   1.10 5.83e-03   1.19 4.09e-05   3.06
0.0063  7941   0:4:41s   |  2.69e-04   1.80 1.66e-03   1.81 5.01e-06   3.03
0.0031 15961  2:28:17s |  6.94e-05   1.95 4.28e-04   1.95 6.20e-07   3.02
-----------------------------------------------------------------------------------------

As we expected, order of convergence is less than optimal

More results for the Laplace eqn on circle, unfitted with Monte-Carlo method.

Various narrow band widths:

1h strip

h Nodes Time |  L2(Gamma)     H1(Gamma)              Area*St.dev
---------------------------------------------------------------------------------------
0.1000   181    0:0:1s   |  1.33e-03             7.30e-02           5.11e-03
0.0500   357    0:0:5s   |  1.92e-04   2.79    2.03e-02   1.84 1.26e-03   2.02
0.0250   709    0:0:6s   |  2.41e-05   2.99 5.29e-03   1.94 1.68e-04   2.90
0.0125  1381   0:0:14s |  2.97e-06   3.02    1.35e-03   1.97 3.11e-05   2.44
0.0063  2817   0:4:17s |  3.73e-07   2.99    3.38e-04   2.00 3.94e-06   2.98
-----------------------------------------------------------------------------------------
same refinement, several monte-carlo runs:
0.0250   709    0:0:6s |  2.98e-05            5.30e-03            1.61e-04
0.0250   709    0:0:6s |  2.73e-05   0.12 5.30e-03   0.00   1.60e-04   0.01
0.0250   709    0:0:6s |  2.60e-05   0.07   5.29e-03   0.00   1.64e-04   -0.03
0.0250   709    0:0:6s |  2.32e-05   0.16   5.29e-03   -0.00 1.63e-04   0.01
0.0250   709    0:0:6s |  2.31e-05   0.00   5.29e-03   0.00   1.64e-04   -0.01
0.0250   709    0:0:6s |  2.32e-05   -0.00 5.29e-03   -0.00  1.64e-04   0.00
0.0250   709    0:0:6s |  2.35e-05   -0.02 5.29e-03   0.00 1.65e-04   -0.01
0.0250   709    0:0:6s |  2.48e-05   -0.08 5.28e-03   0.00  1.67e-04   -0.02
0.0250   709    0:0:6s |  2.63e-05   -0.08 5.30e-03   -0.00 1.65e-04   0.02
0.0250   709    0:0:6s |  2.43e-05   0.11 5.30e-03   0.00   1.64e-04   0.01
-----------------------------------------------------------------------------------------

Same setup with sub-triangulations:

h Nodes Time |  L2(Gamma)     H1(Gamma)     Area*St.dev
---------------------------------------------------------------------------------------
0.1000   181    0:0:1s   |  1.30e-03           7.36e-02           0.00e+00
0.0500   357    0:0:7s   |  1.75e-04   2.89 2.03e-02   1.86 0.00e+00   NaN
0.0250   709   0:0:16s  |  2.31e-05   2.92 5.29e-03   1.94 0.00e+00   NaN
0.0125  1381   0:0:51s |  2.96e-06   2.96 1.35e-03   1.97 0.00e+00   NaN
0.0063  2817   0:3:16s |  3.73e-07   2.99 3.38e-04   2.00 0.00e+00   NaN
-----------------------------------------------------------------------------------------

3h strip
h        Nodes    Time |  L2(Gamma)       H1(Gamma)        Area*St.dev
---------------------------------------------------------------------------------------
0.1000   353    0:0:6s   |  1.43e-03          8.17e-02            5.38e-03
0.0500   697    0:0:6s   |  1.93e-04  2.89 2.25e-02   1.86 1.39e-03   1.95
0.0250  1361   0:0:7s   |  2.56e-05  2.91 5.84e-03   1.95 1.45e-04   3.27
0.0125  2717   0:0:16s |  3.28e-06  2.96 1.49e-03   1.97 2.20e-05   2.71
0.0063  5457   0:4:34s |  4.16e-07  2.98 3.73e-04   2.00 3.37e-06   2.71
-----------------------------------------------------------------------------------------
same refinement, several monte-carlo runs:
0.0250  1361    0:0:8s |  2.55e-05             5.84e-03            1.48e-04
0.0250  1361    0:0:8s |  2.56e-05   -0.00 5.84e-03   -0.00 1.52e-04   -0.04
0.0250  1361    0:0:7s |  2.55e-05   0.00    5.84e-03   0.00   1.45e-04   0.06
0.0250  1361    0:0:8s |  2.57e-05   -0.01 5.84e-03   -0.00 1.45e-04   -0.00
0.0250  1361    0:0:7s |  2.56e-05   0.01    5.84e-03   0.00   1.50e-04   -0.05
0.0250  1361    0:0:7s |  2.56e-05   -0.00 5.84e-03   -0.00 1.49e-04   0.01
0.0250  1361    0:0:7s |  2.56e-05   0.00    5.84e-03   0.00   1.48e-04   0.00
0.0250  1361    0:0:7s |  2.56e-05   -0.00 5.84e-03   0.00   1.49e-04   -0.00
0.0250  1361    0:0:7s |  2.55e-05   0.00    5.84e-03   -0.00 1.49e-04   -0.00
0.0250  1361    0:0:7s |  2.55e-05   0.00    5.84e-03   0.00   1.50e-04   -0.01
-----------------------------------------------------------------------------------------

5h

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

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

0.1000   501    0:0:6s   |  1.58e-03           8.17e-02           9.69e-03 

0.0500  1013    0:0:6s  |  2.06e-04   2.94 2.25e-02   1.86 2.19e-03   2.14 

0.0250  1997    0:0:8s  |  2.56e-05   3.01 5.84e-03   1.95 3.34e-04   2.71 

0.0125  3989   0:0:17s |  3.29e-06   2.96 1.49e-03   1.97 4.08e-05   3.03 

0.0063  7941   0:4:33s |  4.14e-07   2.99 3.73e-04   2.00 5.01e-06   3.03 

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

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.