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.