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.