Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

L2 projection

More
2 years 7 months ago - 2 years 7 months ago #3897 by Macs
L2 projection was created by Macs
Dear NGSpy Team,

I was just testing my L2 projection, projecting uref=x2 on a finite element space V consisting of linear polynomials, but somehow the error (Pu-u,Pu) is not zero, but of order 10^-6.

Here is my code:
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!
Last edit: 2 years 7 months ago by Macs.
More
2 years 7 months ago #3898 by Guosheng Fu
Replied by Guosheng Fu on topic L2 projection
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))
More
2 years 7 months ago #3899 by sneidehx
Replied by sneidehx on topic L2 projection
Hey together,

I think the problem is not the machine precision error.
If [tex]\Pi_k[/tex] is the L2 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
More
2 years 7 months ago #3900 by joachim
Replied by joachim on topic L2 projection
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
Code:
L += uref * v * dx(bonus_intorder=1)

gives you the result at machine precision.

Joachim
The following user(s) said Thank You: Macs, sneidehx
More
2 years 7 months ago #3901 by Macs
Replied by Macs on topic L2 projection
Thank you!
Time to create page: 0.120 seconds