Dear NGSpy Team,
I was just testing my L[sup]2[/sup] projection, projecting u[sub]ref[/sub]=x[sup]2[/sup] on a finite element space V consisting of linear polynomials, but somehow the error (Pu-u,Pu) is not zero, but of order 10^[sup]-6[/sup].
Here is my code:
uref = CoefficientFunction(x*x)
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
V = L2(mesh, order=1)
u = V.TrialFunction()
v = V.TestFunction()
L = LinearForm(V)
L += uref * v * dx
a = BilinearForm(V, symmetric=True)
a += u*v*dx
Puref = GridFunction(V)
Puref.vec.data = a.mat.Inverse(V.FreeDofs(), inverse='sparsecholesky') * L.vec
Could you point me out to what is wrong?
I think the default integration rule for Integrate is 5.
You can try this to get close to machine precision error:
print(Integrate((uref-Puref)*Puref, mesh, order=2))
I think the problem is not the machine precision error.
If [tex]\Pi_k[/tex] is the L[sup]2[/sup] Projection, then
the following should hold:
[tex](\Pi_k u - u, v_h) = 0 \quad \forall v_h \in P^k(\mathcal T) .[/tex]
This should be true for any Integration formula with algebraic exactness higher than k.
I suspect something in the Integration step.
The linear-form is integrated with twice the order of the finite elements. So an integration order of 2 is used, but order 3 would be needed.
Giving some bonus-order as
L += uref * v * dx(bonus_intorder=1)
gives you the result at machine precision.