ElementLinePp(2) plotting incorrect

[Vinc0110]

Comparing a one-dimensional simulation with higher order elements, I noticed a strange behavior of ElementLinePp() when used with p=2. Other values for p seem to be fine. Of course, in this special case it would anyways be better to use ElementLineP2(), as stated in the warning message.

Example code:

import skfem
from skfem.models.poisson import laplace
from skfem.visuals.matplotlib import plot, show

mesh = skfem.MeshLine(p=[0, 0.2, 0.5, 1.0])

#element = skfem.ElementLineP2()
element = skfem.ElementLinePp(2)
basis = skfem.Basis(mesh, element)

A = skfem.asm(laplace, basis)
b = basis.zeros()
phi = basis.zeros()
phi[basis.get_dofs(lambda x: x == 0)] = -0.5
phi[basis.get_dofs(lambda x: x == 1.0)] = 2.0

phi = skfem.solve(*skfem.condense(A, b, x=phi, D=basis.get_dofs()))

ax = plot(basis, phi)
ax.set_title(f'{type(element).__name__}')
show()

[kinnala]

It remainds me of #969. Perhaps related or not.

[kinnala]

Yes, it's again something to do with 1D plotting. Here is a workaround:

ax = plot(basis, phi, nrefs=0)

[Vinc0110]

With nrefs=0 the plot does look correct, thanks for that hint!

Regardless of the plotting, the calculated values are still not identical:

With ElementLineP2():

phi = array([-5.00000000e-01,  8.32667268e-17,  7.50000000e-01,  2.00000000e+00, 
-2.50000000e-01,  3.75000000e-01,  1.37500000e+00])

With ElementLinePp(2):

phi = array([-5.00000000e-01,  5.55111512e-17,  7.50000000e-01,  2.00000000e+00,
0.00000000e+00,  0.00000000e+00,  0.00000000e+00])

By the way, my code above was slightly wrong: I messed up the boundary identification.
Replace

phi[basis.get_dofs(lambda x: x == 0)] = -0.5
phi[basis.get_dofs(lambda x: x == 1.0)] = 2.0

with

phi[basis.get_dofs(lambda x: x[0] == 0)] = -0.5
phi[basis.get_dofs(lambda x: x[0] == 1.0)] = 2.0

[kinnala]

I believe that the difference is because ElementLineP2 uses nodal basis and ElementLinePp uses hierarchical basis.

[Vinc0110]

I'm not entirely sure what that means. Is ElementLinePp(2) simply different from ElementLineP2() and cannot be used as a universal (but slower) substitute?
In a larger simulation, where I'm trying to calculate eigenmodes in 1d, I also get different results when using ElementLinePp(2) instead of ElementLineP2().

[kinnala]

In most applications you can use either one of those but they don't give exactly the same result because the basis functions are different. There are infinite different bases you can use to represent quadratic polynomials. There are many different choices that are used in FEM depending on the use case.

[Vinc0110]

Understood. Thanks for the help!

[kinnala]

We can leave this open because of the bug in skfem.visuals.matplotlib.plot for MeshLine1.