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.