P2 moment fitting

The problems appear to be with boundary triangles that are very small compared to mesh size.

Old method with $h^{-1}$ subtriangles fitting the boundary

h        Nodes   Time      |   L2(Gamma)            H1(Gamma)    
---------------------------------------------------------------------------------------
0.1000      10    0:0:1s    |  1.44e-03                8.17e-02
0.0500      20    0:0:4s    |  1.93e-04   2.91      2.25e-02   1.86
0.0250      40   0:0:16s   |  2.56e-05   2.91      5.84e-03   1.95
0.0125      80    0:1:2s    |  3.28e-06   2.96      1.49e-03   1.97
0.0063     160   0:4:25s  |  4.14e-07   2.99      3.73e-04   2.00
0.0031     320  0:22:13s |  5.21e-08   2.99      9.38e-05   1.99
----------------------------------------------------------------------------------

Moments fitting P2
\begin{array}[h] {lrr|rr|rr} h & d.o.f & time & L_2(\Gamma) & & H_1(\Gamma) & \\ \hline 0.1000 &    10 &    0:0:1s &  2.58e-03  & & 8.16e-02 & \\0.05000 &    20 &    0:0:2s &  7.33e-04 & 1.82 & 2.29e-02 & 1.84 \\0.02500 &    40 &    0:0:4s &  1.47e-04 & 2.32 & 9.77e-03 & 1.23 \\0.01250 &    80 &    0:0:8s &  5.25e-05 & 1.49 & 7.82e-03 & 0.32 \\0.00625 &   160 &   0:0:15s &  8.28e-06 & 2.66 & 1.66e-03 & 2.24 \\0.00313 &   320 &   0:0:31s &  2.47e-06 & 1.74 & 1.05e-03 & 0.67 \\0.00156 &   640 &    0:1:2s &  1.82e-06 & 0.44 & 1.75e-03 & -0.74 \\0.00078 &  1280 &    0:2:6s &  2.34e-07 & 2.96 & 4.21e-04 & 2.06 \\0.00039 &  2560 &   0:4:14s &  2.81e-07 & -0.27 & 9.64e-04 & -1.20 \\0.00020 &  5120 &   0:8:42s &  2.47e-06 & -3.14 & 2.11e-02 & -4.45 \\0.00010 & 10240 &  0:19:22s &  1.82e-06 & 0.45 & 4.93e-02 & -1.22 \\\end{array}

P1 FEM with moments fitting integration on boundary

To figure out the problem with P2 MF scheme, I went back to P1 elements and use the moments fitting scheme of h^2 precision. The idea works, but there appears to be some numerical problem when I get to higher precision. The L2 error stagnates at 1E-8.

Update: this is not an error. We have seen this problem before - the strip width of 2h is just too small. 5h width has better results.

Strip width = 2h (surface $\pm h$)
Old method
\begin{array}[h] {lrr|rr|rr} h & d.o.f & time & L_2(\Gamma) & & H_1(\Gamma) & \\ \hline 0.1000 &    10 &    0:0:0s &  5.50e-02  & & 1.23e+00 & \\0.05000 &    20 &    0:0:1s &  1.55e-02 & 1.83 & 6.19e-01 & 0.99 \\0.02500 &    40 &    0:0:1s &  4.28e-03 & 1.86 & 3.19e-01 & 0.96 \\0.01250 &    80 &    0:0:2s &  1.08e-03 & 1.99 & 1.60e-01 & 0.99 \\0.00625 &   160 &    0:0:3s &  2.64e-04 & 2.02 & 7.98e-02 & 1.00 \\0.00313 &   320 &    0:0:6s &  6.68e-05 & 1.98 & 3.99e-02 & 1.00 \\0.00156 &   640 &   0:0:13s &  1.59e-05 & 2.07 & 1.99e-02 & 1.00 \\0.00078 &  1280 &   0:0:26s &  4.22e-06 & 1.92 & 9.90e-03 & 1.01 \\0.00039 &  2560 &   0:0:50s &  1.36e-06 & 1.63 & 4.91e-03 & 1.01 \\0.00020 &  5120 &   0:1:40s &  3.10e-07 & 2.14 & 2.41e-03 & 1.03 \\0.00010 & 10240 &   0:3:22s &  8.48e-08 & 1.87 & 1.16e-03 & 1.05 \\0.00005 & 20480 &   0:6:51s &  3.69e-08 & 1.20 & 5.35e-04 & 1.12 \\0.00002 & 40960 &   0:14:6s &  8.42e-09 & 2.13 & 2.23e-04 & 1.26 \\\end{array}

Moments fitting Strip width = 2h (surface $\pm h$)
\begin{array}[h] {lrr|rr|rr} h & d.o.f & time & L_2(\Gamma) & & H_1(\Gamma) & \\ \hline 0.1000 &    10 &   0:0:1s &  5.74e-02  & & 1.26e+00 & \\0.05000 &    20 &    0:0:1s &  1.40e-02 & 2.04 & 6.29e-01 & 1.01 \\0.02500 &    40 &    0:0:2s &  3.68e-03 & 1.92 & 3.20e-01 & 0.97 \\0.01250 &    80 &    0:0:3s &  8.99e-04 & 2.03 & 1.60e-01 & 1.00 \\0.00625 &   160 &    0:0:6s &  2.39e-04 & 1.91 & 7.98e-02 & 1.00 \\0.00313 &   320 &   0:0:11s &  5.82e-05 & 2.04 & 3.99e-02 & 1.00 \\0.00156 &   640 &   0:0:23s &  1.65e-05 & 1.82 & 1.99e-02 & 1.00 \\0.00078 &  1280 &   0:0:45s &  4.83e-06 & 1.77 & 9.90e-03 & 1.01 \\0.00039 &  2560 &   0:1:30s &  1.24e-06 & 1.96 & 4.91e-03 & 1.01 \\0.00020 &  5120 &    0:3:5s &  4.10e-07 & 1.60 & 2.41e-03 & 1.03 \\0.00010 & 10240 &   0:6:40s &  1.28e-07 & 1.68 & 1.16e-03 & 1.05 \\0.00005 & 20480 &   0:15:6s &  4.65e-08 & 1.46 & 5.35e-04 & 1.12 \\0.00002 & 40960 &  0:36:21s &  1.49e-08 & 1.64 & 2.23e-04 & 1.26 \\0.00001 & 81920 &  1:37:48s &  2.78e-08 & -0.90 & 6.74e-05 & 1.72 \\0.00001 & 163840 &   4:54:3s &  3.20e-08 & -0.20 & 1.25e-05 & 2.43 \\\end{array}

Old method Strip width = 5h (surface $\pm 2.5 h$)
\begin{array}[h] {lrr|rr|rr} h & d.o.f & time & L_2(\Gamma) & & H_1(\Gamma) & \\ \hline 0.1000 &    10 &    0:0:3s &  6.19e-02  & & 1.27e+00 & \\0.05000 &    20 &    0:0:4s &  1.51e-02 & 2.03 & 6.29e-01 & 1.01 \\0.02500 &    40 &    0:0:6s &  3.78e-03 & 2.00 & 3.20e-01 & 0.97 \\0.01250 &    80 &    0:0:5s &  9.32e-04 & 2.02 & 1.60e-01 & 1.00 \\0.00625 &   160 &    0:0:7s &  2.34e-04 & 1.99 & 7.99e-02 & 1.00 \\0.00313 &   320 &   0:0:12s &  5.88e-05 & 2.00 & 3.99e-02 & 1.00 \\0.00156 &   640 &   0:0:21s &  1.46e-05 & 2.01 & 2.00e-02 & 1.00 \\0.00078 &  1280 &   0:0:36s &  3.65e-06 & 2.00 & 9.98e-03 & 1.00 \\0.00039 &  2560 &    0:1:6s &  8.82e-07 & 2.05 & 4.99e-03 & 1.00 \\0.00020 &  5120 &    0:2:8s &  2.34e-07 & 1.91 & 2.49e-03 & 1.00 \\0.00010 & 10240 &   0:4:15s &  5.41e-08 & 2.11 & 1.24e-03 & 1.01 \\0.00005 & 20480 &   0:8:20s &  1.34e-08 & 2.02 & 6.16e-04 & 1.01 \\\end{array}

Moments fitting Strip width = 5h (surface $\pm 2.5 h$)
\begin{array}[h] {lrr|rr|rr} h & d.o.f & time & L_2(\Gamma) & & H_1(\Gamma) & \\ \hline 0.1000 &    10 &    0:0:4s &  6.21e-02  & & 1.27e+00 & \\0.05000 &    20 &    0:0:5s &  1.51e-02 & 2.04 & 6.30e-01 & 1.01 \\0.02500 &    40 &    0:0:8s &  3.73e-03 & 2.02 & 3.20e-01 & 0.98 \\0.01250 &    80 &    0:0:8s &  9.34e-04 & 2.00 & 1.60e-01 & 1.00 \\0.00625 &   160 &   0:0:13s &  2.28e-04 & 2.03 & 7.99e-02 & 1.00 \\0.00313 &   320 &   0:0:22s &  5.74e-05 & 1.99 & 3.99e-02 & 1.00 \\0.00156 &   640 &   0:0:40s &  1.44e-05 & 2.00 & 2.00e-02 & 1.00 \\0.00078 &  1280 &   0:1:15s &  3.61e-06 & 1.99 & 9.98e-03 & 1.00 \\0.00039 &  2560 &   0:2:28s &  9.15e-07 & 1.98 & 4.99e-03 & 1.00 \\0.00020 &  5120 &   0:4:51s &  2.26e-07 & 2.02 & 2.49e-03 & 1.00 \\0.00010 & 10240 &  0:10:18s &  5.62e-08 & 2.01 & 1.24e-03 & 1.01 \\0.00005 & 20480 &  0:22:33s &  1.43e-08 & 1.97 & 6.16e-04 & 1.01 \\\end{array}

Moments fitting, H=0
\begin{array}[h] {lrr|rr|rr} h & d.o.f & time & L_2(\Gamma) & & H_1(\Gamma) & \\ \hline 0.1000 &    10 &    0:0:1s &  5.27e-02  & & 1.23e+00 & \\0.0500 &    20 &    0:0:1s &  1.46e-02 & 1.85 & 6.19e-01 & 0.99 \\0.0250 &    40 &    0:0:2s &  4.11e-03 & 1.83 & 3.19e-01 & 0.96 \\0.0125 &    80 &    0:0:3s &  1.02e-03 & 2.01 & 1.60e-01 & 1.00 \\0.0063 &   160 &    0:0:5s &  2.69e-04 & 1.92 & 7.97e-02 & 1.00 \\0.0031 &   320 &   0:0:11s &  6.59e-05 & 2.03 & 3.99e-02 & 1.00 \\0.0016 &   640 &   0:0:21s &  1.84e-05 & 1.84 & 1.99e-02 & 1.00 \\0.0008 &  1280 &   0:0:43s &  5.21e-06 & 1.82 & 9.90e-03 & 1.01 \\0.0004 &  2560 &   0:1:25s &  1.32e-06 & 1.98 & 4.91e-03 & 1.01 \\0.0002 &  5120 &   0:2:58s &  4.30e-07 & 1.61 & 2.41e-03 & 1.03 \\\end{array}

Moment Fitting quadrature on implicit surface (Muller)

I found an error in my Surface integration weights calculations. In some triangles the direction of the normal was found wrong.

Instead of using distance function at location of the normal, I should make sure dot product of normal vector at each quadrature point with normal of interface edge is positive.

Update: fixed above error.

Description. 

I'm integrating f(x) over the zero level set of a narrow-band region using the quadrature method by Muller. f(x) is either a constant function or $cos(9 \theta)$ in 1st quadrant.

The scheme produced should be precise for polynomials up to 2nd order in each triangle, so the total error should be of order O(h^2). 

Nodes are taken from a 13-point Gaussian quadrature over triangles. 1-D integration over the edges is done using Gaussian rule precise for polynomials of order 9.

Below table is for mesh refinements.

h Time       |  f=1 error                             f=$cos(9 \theta)$ error 

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

0.100000    0:0:0s  |  2.95e-06                   9.92e-05    

0.050000    0:0:1s  |  4.22e-09         9.45 3.87e-06 4.68 

0.025000    0:0:1s  |  1.01e-08        -1.25 5.29e-07 2.87 

0.012500    0:0:2s  |  6.17e-10         4.03 7.47e-09 6.15 

0.006250    0:0:4s  |  1.17e-10         2.40 2.38e-09 1.65 

0.003125    0:0:7s  |  4.90e-12         4.57 3.44e-10 2.79 

0.001563   0:0:15s |  3.48e-13         3.82 1.28e-11 4.74 

0.000781   0:0:31s |  7.99e-15         5.44 2.11e-13 5.93 

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

I'm somewhat puzzled how to interpret these results. ...

More on previous post about H=0 run:

The table in previous post actually comes from H=0, 5h strip width experiment.
:
5h strip, H=0
h Nodes Time |  L2(Gamma)     H1(Gamma)     Area*St.dev
---------------------------------------------------------------------------------------
0.1000   501    0:0:6s   |  4.20e-02   2.60e-01   9.82e-03
0.0500  1013    0:0:6s  |  5.01e-03   3.07 3.49e-02   2.89 2.32e-03   2.08
0.0250  1997    0:0:8s  |  2.02e-03   1.31 1.33e-02   1.40 3.40e-04   2.77
0.0125  3989   0:0:16s |  9.38e-04   1.10 5.83e-03   1.19 4.10e-05   3.05
0.0063  7941   0:4:21s |  2.69e-04   1.80 1.66e-03   1.81 5.00e-06   3.04
-----------------------------------------------------------------------------------------

5h strip, exact H
h Nodes Time |  L2(Gamma)     H1(Gamma)     Area*St.dev
---------------------------------------------------------------------------------------
0.1000   501    0:0:6s   |  2.61e-03   8.17e-02   1.00e-02
0.0500  1013    0:0:6s  |  2.08e-04   3.65 2.25e-02   1.86 2.11e-03   2.25
0.0250  1997    0:0:8s  |  2.56e-05   3.02 5.84e-03   1.95 3.44e-04   2.62
0.0125  3989   0:0:16s |  3.28e-06   2.96 1.49e-03   1.97 4.10e-05   3.07
0.0063  7941   0:4:21s |  4.14e-07   2.99 3.73e-04   2.00 5.00e-06   3.04
-----------------------------------------------------------------------------------------

Surprisingly, the narrower strip width gives the expected convergence of P1 method:
strip width 1h, H=0
h Nodes Time |  L2(Gamma)     H1(Gamma)     Area*St.dev
---------------------------------------------------------------------------------------
0.1000   181    0:0:1s |  1.06e-02   7.28e-02   7.08e-03
0.0500   357    0:0:5s |  7.10e-04   3.89 2.10e-02   1.79 1.21e-03   2.54
0.0250   709    0:0:6s |  1.79e-04   1.98 5.58e-03   1.92 1.66e-04   2.87
0.0125  1381   0:0:15s |  4.49e-05   2.00 1.43e-03   1.97 3.10e-05   2.42
0.0063  2817   0:4:13s |  1.13e-05   1.99 3.58e-04   1.99 3.95e-06   2.97
-----------------------------------------------------------------------------------------

I am not sure why this would happen.

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.

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.

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: