L2Projector not workingwith `Pyramid`s

[AbdAlazezAhmed]

using Ferrite
grid = generate_grid(Pyramid, (1,1,1));
ip = Lagrange{RefPyramid, 1}()
proj = L2Projector(ip, grid)

ERROR: PosDefException: matrix is not positive definite; Factorization failed.
Stacktrace:
 [1] cholesky!(F::SparseArrays.CHOLMOD.Factor{Float64, Int64}, A::SparseArrays.CHOLMOD.Sparse{Float64, Int64}; shift::Float64, check::Bool)
   @ SparseArrays.CHOLMOD ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/SparseArrays/src/solvers/cholmod.jl:1462
 [2] cholesky!
   @ ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/SparseArrays/src/solvers/cholmod.jl:1455 [inlined]
 [3] #cholesky#14
   @ ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/SparseArrays/src/solvers/cholmod.jl:1497 [inlined]
 [4] cholesky
   @ ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/SparseArrays/src/solvers/cholmod.jl:1490 [inlined]
 [5] cholesky
   @ ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/SparseArrays/src/solvers/cholmod.jl:1609 [inlined]
 [6] close!(proj::L2Projector)
   @ Ferrite ~/.julia/dev/Ferrite/src/L2_projection.jl:135
 [7] #L2Projector#507
   @ ~/.julia/dev/Ferrite/src/L2_projection.jl:90 [inlined]
 [8] L2Projector(ip::Lagrange{RefPyramid, 1}, grid::Grid{3, Pyramid, Float64})
   @ Ferrite ~/.julia/dev/Ferrite/src/L2_projection.jl:80
 [9] top-level scope

[lijas]

Does it work by just increasing the number of quadrature points?

[AbdAlazezAhmed]

Interesting, yes starting from order 3, haven't checked the tests tho

proj = L2Projector(ip, grid; qr_lhs = QuadratureRule{RefPyramid}(3))

L2Projector
  projection on:           6/6 cells in grid
  function interpolation:  Lagrange{RefPyramid, 1}()
  geometric interpolation: Lagrange{RefPyramid, 1}()

[lijas]

For second order Pyramids, "5th order" integration rule seems to work.

function _mass_qr(::Lagrange{shape, 1}) where {shape <: RefPyramid}
    return QuadratureRule{shape}(3)
end
function _mass_qr(::Lagrange{shape, 2}) where {shape <: RefPyramid}
    return QuadratureRule{shape}(5)
end

Can you add this to your PR #1215 that will test for this?

[termi-official]

We probably need different quadrature rules here. The Ansatz functions on pyramids are simply not polynomials (and the name Lagrange for the interpolation here is very misleading, as there are no Lagrange polynomials involved).
See

Ferrite.jl/src/interpolations.jl

Lines 1239 to 1242 in 0e96cf1

	i == 1 && return zzero ? zero(T) : (-x * y + (z - 1) * (-x - y - z + 1)) / (z - 1)
	i == 2 && return zzero ? zero(T) : x * (y + z - 1) / (z - 1)
	i == 3 && return zzero ? zero(T) : y * (x + z - 1) / (z - 1)
	i == 4 && return zzero ? zero(T) : -x * y / (z - 1)

We should probably deprecate the usage with a warning and rename it. And in addition we can provide a different approach (e.g. https://onlinelibrary.wiley.com/doi/abs/10.1002/cnm.1070) and provide corresponding quadrature rules for these.

[lijas]

What exactly makes an element "Lagrage"? The shape functions are from here which clasifies it as Lagrange: https://defelement.org/elements/examples/pyramid-lagrange-equispaced-1.html

[termi-official]

Good question. After some search there seem to be two different definitons around.

1. A finite element using Lagrange interpolation polynomials as the basis functions
2. A nodal finite element with only evaluations of the basis functions (i.e. no derivative evaluations of the basis functions as e.g. in Hermite elements)

After reading up more I kinda lean towards the latter type of definition, although I would restrict it a bit more to separate this class from e.g. Crouzeix-Raviart type elements. However, I still think it is misleading, as most users would probably also assume that the functions are just polynomials. This also raises the question for me why we should prefer this rational ansatz space instead of the (piecewise) polynomial ansatz space (as we can quite easily construct composed quadrature rules).

[KnutAM]

although I would restrict it a bit more to separate this class from e.g. Crouzeix-Raviart type elements

We also have Serendipity...