- Thank you received: 0
L2 projection
3 years 3 months ago - 3 years 3 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:
Could you point me out to what is wrong?
Thank you!
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!
Attachments:
Last edit: 3 years 3 months ago by Macs.
- Guosheng Fu
-
- Offline
- Elite Member
-
Less
More
- Thank you received: 6
3 years 3 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))
You can try this to get close to machine precision error:
print(Integrate((uref-Puref)*Puref, mesh, order=2))
3 years 3 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
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
3 years 3 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
gives you the result at machine precision.
Joachim
Giving some bonus-order as
Code:
L += uref * v * dx(bonus_intorder=1)
gives you the result at machine precision.
Joachim
Time to create page: 0.104 seconds