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.