L2 projection

Main Forum Old-Forum - do not post here
1

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[1]-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))

Could you point me out to what is wrong?

Thank you!

Attachment: testL2_2021-08-12.py

  1. sup ↩︎

2

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))

3

Hey together,

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.

Best

4

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.

Joachim

5

Thank you!