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)
a.Assemble()
L.Assemble()
Puref.vec.data = a.mat.Inverse(V.FreeDofs(), inverse='sparsecholesky') * L.vec
print(Integrate((uref-Puref)*Puref, mesh))

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.