Lagrange elements of higher order

[edljk]

I wonder if there are some plans to generalize interpolation.jl (and whatever needed) to higher order Lagrange elements on simplices (2D/3D)?
Following Getem++ Appendix, a naive implementation in 2D for instance would be the following functions.
Trying to test these "new" elements (like P3), I obtained an assertion error ERROR: AssertionError: ip_info.nedgedofs <= 1, etc..
Do you have any idea of the required steps to make these new elements available?
Thanks

function Ferrite.getnbasefunctions(::Lagrange{2,RefTetrahedron,order}) where order 
    return (order + 1) * (order + 2) ÷ 2
end
Ferrite.nvertexdofs(::Lagrange{2,RefTetrahedron,order}) where order = 1
Ferrite.nfacedofs(::Lagrange{2,RefTetrahedron,order}) where order = order - 1
function Ferrite.vertices(::Lagrange{2,RefTetrahedron,order}) where order 
    return (1, 1 + order, (order + 1) * (order + 2) ÷ 2)
end
function Ferrite.faces(::Lagrange{2,RefTetrahedron,order}) where order     
    return (((1:(order + 1))...,), (cumsum((order + 1):-1:1)...,),
            ([1, (1 .+ cumsum((order + 1):-1:2))...]...,))
end
function Ferrite.reference_coordinates(::Lagrange{2,RefTetrahedron,order}) where order
    coordpts = Vector{Ferrite.Vec{2, Float64}}()
    for k = 0:order
        for l = 0:(order - k)
            push!(coordpts, Ferrite.Vec{2, Float64}((l / order, k / order)))
        end
    end
    return coordpts
end
function Ferrite.value(ip::Lagrange{2,RefTetrahedron,order}, i::Int, 
                           ξ::Ferrite.Vec{2}) where order
    ξ1, ξ2, ξ3= ξ[1], ξ[2], 1. - sum(ξ)
    i1, i2, i3 = _numlin_basis2DT(i, order) 
    val = one(typeof(ξ1))
    i1 ≥ 1 && (val *= prod((order * ξ1 - j) / (j + 1) for j = 0:(i1 - 1))) 
    i2 ≥ 1 && (val *= prod((order * ξ2 - j) / (j + 1) for j = 0:(i2 - 1)))
    i3 ≥ 1 && (val *= prod((order * ξ3 - j) / (j + 1) for j = 0:(i3 - 1)))
    return val
end
function _numlin_basis2DT(i, order) 
    c, j1, j2, j3 = 0, 0, 0, 0
    for k = 0:order
        if i <= c + (order + 1 - k)
            j2 = i - c - 1
            break
        else
            j3 += 1
            c += order + 1 - k
        end
    end
    j1 = order - j2 -j3
    return j1, j2, j3
end

[termi-official]

Hi @edljk ,
adding higher order elements is on my todo, but I cannot find time to implement them. Your implementation for higher order Lagrangian basis functions looks fine. However, the problem is that the current dof assignment algorithm does not play well with higher order interpolations. What you need is correct orientation information of your dof assignment along edges and faces.
To elaborate on th eproblem, think about having two triangles forming a quad, then the shared edge has a different orientation for each triangle. In this case, there the current dof assignment algorithm would match the wrong dofs, because the orientation information is neglected. As far as I am aware, this should be the "only" thing that has to be fixed. Or have I missed something @koehlerson ?

[edljk]

Hi @termi-official,
thank you very much for the feedback.
Does the following faces definition fix the orientation or something has still to be fixed somewhere else?

function Ferrite.faces(::Lagrange{2,RefTetrahedron,order}) where order     
    return (((1:(order + 1))...,), (cumsum((order + 1):-1:1)...,),
    ((order + 1) * (order + 2) ÷ 2 .- cumsum([0, collect(2:(order + 1))...])...,))
end

I am not sure to understand what is expected for the dof assignment along faces in 2D? Was it a 3D observation?

[termi-official]

Sorry, I think I was not specific enough on this matter. I do not think that we can fix this simply by reorienting the faces, because this will break other stuff related to face orientations. If I remember correctly, then the problem has to be fixed in the __close! function.

[fredrikekre]

However, the problem is that the current dof assignment algorithm does not play well with higher order interpolations. What you need is correct orientation information of your dof assignment along edges and faces.

Yes, and this is exactly why the assertion exist:

I obtained an assertion error ERROR: AssertionError: ip_info.nedgedofs <= 1, etc..

because as long as there is just one dof for an edge the orientation doesn't matter.

[edljk]

Thank you fredrikekre.
OK, I understood that a key point is the orientation that should be considered more attentively in some funtion of Ferrite.jl.
Do you have an idea of which part of the code would require a more specific attention to make higher order elements work? What should be the first steps for that improvement?

[termi-official]

Let me comment on this since I put some time into this. If I see it correctly, then, as long as you are only interested in H1 spaces, fixing the assignment in __close! should be everything. What can be done here is a check if an edge or face already has dofs assigned - in this case the dof assignment has to be permutated correctly. The simplest approach to find the correct permutations is to define tables, if you just need some special cases. IIRC then more elaborate algorithms should be defined in the FEM book by Pavel Solin (Higher-Order Finite Element Methods). Hope that I could help here.

[edljk]

Thank you termi-official, the situation is clear now.
I commented a few assertions and gave it a try. I wanted to test previous Pk implementation to solve a simple H1 eigenvalue problem (Neumann eigenvalues of the square..)
I paste the code below.
At least the number of degrees of freedom seems to be ok but there is still something wrong. Most probably in my implementation or perhaps with the compatibility of the numbering of dofs and vertices.. Do not really know. I will continue to investigate.

julia> test_H1Pk(degfem = 1)
real.(λ) = [-5.3759690384243275e-15, 2.4724574365104006, 2.4724575455509052, 4.965112609288587, 9.950481029670694, 9.950806702685107, 12.466769445826522, 12.540827949304449, 20.220092594807863, 22.61794567596597]
error = 0.20382244794275017
nbdof = 441

julia> test_H1Pk(degfem = 2)
real.(λ) = [6.052240143512161e-14, 2.46740315865208, 2.467403161776321, 4.934830898153275, 9.869735676495825, 9.869735679677525, 12.33727069978334, 12.337492721885768, 19.741019778974213, 22.208093643681096]
error = 0.000487220524069798
nbdof = 1681

julia> test_H1Pk(degfem = 3)
real.(λ) = [75.29284195065094, 81.6602541282197, 90.55757001630815, 93.29508764352964, 106.32266184039761, 118.07857760635417, 122.39257092374741, 138.9103367965305, 139.3015386333036, 148.14236328416197]
error = 126.57333129516879
nbdof = 2921

using Ferrite, SparseArrays, Arpack

function test_H1Pk(; degfem::Int64 = 2, degcurve::Int64 = 1, degquad::Int64 = 6,
                     dim::Int64 = 2)
    typecellq = degcurve == 2 ? QuadraticTriangle : Ferrite.Triangle
    grid = generate_grid(typecellq, (20, 20))
    typecellq = dim == 2 ? QuadraticTriangle : QuadraticTetrahedron
    ip = Ferrite.Lagrange{dim, RefTetrahedron, degfem}()
    qr = Ferrite.QuadratureRule{dim, RefTetrahedron}(:legendre, degquad)
    ip_geom = Ferrite.Lagrange{dim, RefTetrahedron, degcurve}() 
    cellvalues = Ferrite.CellScalarValues(qr, ip, ip_geom)
    dh = Ferrite.DofHandler(grid)
    Ferrite.push!(dh, :u, 1, ip)
    Ferrite.close!(dh)

    K = create_sparsity_pattern(dh)
    doassemble_K_ferrite!(K, cellvalues, dh)
    M = create_sparsity_pattern(dh)
    doassemble_M_ferrite!(M, cellvalues, dh)

    # compute eigenvalues
    λ, U = eigs(K, M, nev = 10, which = :LM, tol = 0., maxiter = 30_000, 
                sigma = rand() / 1e4) 
    # display
    error = maximum(abs.([0, 1/4, 1/4, 1/2, 1, 1, 5/4, 5/4] * pi ^ 2 - λ[1:8]))
    @show real.(λ)
    @show error
    nbdof = size(K, 1)
    @show nbdof
    return 
end

function doassemble_K_ferrite!(K::SparseMatrixCSC, 
                       cellvalues::CellScalarValues{dim}, 
                       dh::DofHandler) where {dim}
    n_basefuncs = getnbasefunctions(cellvalues)
    Ke = zeros(n_basefuncs, n_basefuncs)
    assembler = start_assemble(K)
    @inbounds for cell in CellIterator(dh)
        fill!(Ke, 0)
        reinit!(cellvalues, cell)
        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            for i in 1:n_basefuncs
                v  = shape_value(cellvalues, q_point, i)
                ∇v = shape_gradient(cellvalues, q_point, i)
                for j in 1:n_basefuncs
                    ∇u = shape_gradient(cellvalues, q_point, j)
                    Ke[i, j] += (∇v ⋅ ∇u) * dΩ
                end
            end
        end
        assemble!(assembler, celldofs(cell), Ke)
    end
    return K
end
function doassemble_M_ferrite!(M::SparseMatrixCSC, 
                       cellvalues::CellScalarValues{dim}, dh::DofHandler;
                       ρ::Union{Function, Nothing} = nothing) where {dim}
    n_basefuncs = getnbasefunctions(cellvalues)
    Me = zeros(n_basefuncs, n_basefuncs)
    assembler = start_assemble(M)
    valρ = 1.
    @inbounds for cell in CellIterator(dh)
        fill!(Me, 0)
        reinit!(cellvalues, cell)
        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            if ρ != nothing 
                pq = spatial_coordinate(cellvalues, q_point, 
                                        getcoordinates(cell))
                valρ = ρ(pq)
            end
            for i in 1:n_basefuncs
                v  = shape_value(cellvalues, q_point, i) 
                for j in 1:n_basefuncs
                    u = shape_value(cellvalues, q_point, j)                   
                    Me[i, j] += (v * u) * valρ * dΩ
                end
            end
        end
        assemble!(assembler, celldofs(cell), Me)
    end
    return M
end

#--------------------------------------------------------------------------------------------------------
#----------------------------------------- Pk elements --------------------------------------------------
#--------------------------------------------------------------------------------------------------------
function Ferrite.getnbasefunctions(::Lagrange{2,RefTetrahedron,order}) where order 
    return (order + 1) * (order + 2) ÷ 2
end

Ferrite.nvertexdofs(::Lagrange{2,RefTetrahedron,order}) where order = 1
#Ferrite.nedgedofs(::Lagrange{2,RefTetrahedron,order}) where order = order + 1
Ferrite.nfacedofs(::Lagrange{2,RefTetrahedron,order}) where order = order - 1

function Ferrite.vertices(::Lagrange{2,RefTetrahedron,order}) where order 
    return (1 + order, (order + 1) * (order + 2) ÷ 2, 1)
end

function Ferrite.faces(::Lagrange{2,RefTetrahedron,order}) where order     
    return (((1:(order + 1))...,), (cumsum((order + 1):-1:1)...,),
            ((order + 1) * (order + 2) ÷ 2 .- cumsum([0, collect(2:(order + 1))...])...,))
    #=
    return (reverse(((1:(order + 1))...,)), 
            reverse(((order + 1) * (order + 2) ÷ 2 .- cumsum([0, collect(2:(order + 1))...])...,)),
            reverse((cumsum((order + 1):-1:1)...,)))
    =#
end

function Ferrite.reference_coordinates(::Lagrange{2,RefTetrahedron,order}) where order
    coordpts = Vector{Ferrite.Vec{2, Float64}}()
    for k = 0:order
        for l = 0:(order - k)
            push!(coordpts, Ferrite.Vec{2, Float64}((l / order, k / order)))
        end
    end
    return coordpts
end

function Ferrite.value(ip::Lagrange{2,RefTetrahedron,order}, i::Int, 
                           ξ::Ferrite.Vec{2}) where order
    ξ1, ξ2, ξ3= 1. - sum(ξ), ξ[1], ξ[2] 
    i1, i2, i3 = _numlin_basis2DT(i, order) 
    val = one(typeof(ξ1))
    i1 ≥ 1 && (val *= prod((order * ξ1 - j) / (j + 1) for j = 0:(i1 - 1))) 
    i2 ≥ 1 && (val *= prod((order * ξ2 - j) / (j + 1) for j = 0:(i2 - 1)))
    i3 ≥ 1 && (val *= prod((order * ξ3 - j) / (j + 1) for j = 0:(i3 - 1)))
    return val
end
function _numlin_basis2DT(i, order) 
    c, j1, j2, j3 = 0, 0, 0, 0
    for k = 0:order
        if i <= c + (order + 1 - k)
            j2 = i - c - 1
            break
        else
            j3 += 1
            c += order + 1 - k
        end
    end
    j1 = order - j2 -j3
    return j1, j2, j3
end

[edljk]

I identified one mistake in previous code (the barycentric coordinates ξ1, ξ2, ξ3 were not correctly defined). It is updated.
I tried to reproduce the P2 (correct) results of the original Ferrite implementation with this "new" implementation.
I observed that the result depends on the ordering of the reference_coordinates (which of course should be compatible with the value function).
I didn't understand what are the constraints that this ordering should satisfy to be compatible with the assembling procedures in an H1 discretization.. Does someone has an idea? Thanks.

[fredrikekre]

1. Ordering of reference_coordinates must match the ordering of the shape functions (i.e. in the value function)
2. Ordering of shape functions must be
   1. vertices
   2. edges
   3. faces
   4. internal

In particular this means that the first N dofs are the same for linear and quadratic, and first M dofs the same for quad and third order etc

[fredrikekre]

You can test 1. with:

julia> using Ferrite, Test

julia> interpolation = Lagrange{3,RefTetrahedron,2}()
Lagrange{3, RefTetrahedron, 2}()

julia> for i in 1:getnbasefunctions(interpolation)
           for (j, x) in enumerate(Ferrite.reference_coordinates(interpolation))
               val = Ferrite.value(interpolation, i, x)
               if i == j
                   @test val == 1
               else
                   @test val == 0
               end
           end
       end

[edljk]

Thank you @fredrikekre for the comments and the test function. It seems that my previous piece of code passes the tests.
I am still not sure to understand the ordering:

• I see that the dofs of order 1 are also dofs of order 2 elements. But this is not true between order 2 and 3, is it?..
• To be perhaps more explicit. Does the following ordering look correct to you for order 3?

julia> ip = Lagrange{2, RefTetrahedron, 3}()
julia> Ferrite.getnbasefunctions(ip)
10
julia> Ferrite.nvertexdofs(ip)
1
julia> Ferrite.nedgedofs(ip) # correct?
2
julia> Ferrite.nfacedofs(ip) # correct?
2
julia> Ferrite.vertices(ip)
(1, 2, 3)
julia> Ferrite.faces(ip)
((1, 2, 4, 5), (2, 3, 6, 7), (3, 1, 8, 9))

julia> Ferrite.reference_coordinates(ip)

10×1 Matrix{Vec{2, Float64}}:
 [1.0, 0.0]
 [0.0, 1.0]
 [0.0, 0.0]
 [0.6666666666666666, 0.3333333333333333]
 [0.3333333333333333, 0.6666666666666666]
 [0.0, 0.6666666666666666]
 [0.0, 0.3333333333333333]
 [0.3333333333333333, 0.0]
 [0.6666666666666666, 0.0]
 [0.3333333333333333, 0.3333333333333333]

[fredrikekre]

But this is not true between order 2 and 3, is it?

No, probably this can not hold for order 3

To be perhaps more explicit. Does the following ordering look correct to you for order 3?

Yea, I believe that all looks fine.

Can you share your implementation, maybe open it as a WIP pull request? That might make it easier to communicate and play around with it.

[fredrikekre]

WIP PR opened here for reference: #482

[edljk]

Curved mesh of order > 2 not implemented yet. Not sure what should be done. Will try to have a look. Thanks.

[termi-official]

That should not be an issue. The geometrical approximation is decoupled from the algebraic approximation. With push! you can simply provide your favorite approximation (see https://github.com/Ferrite-FEM/Ferrite.jl/blob/master/src/Dofs/DofHandler.jl#L115) while you provide the geometric interpolator separately to the Values object (see https://github.com/Ferrite-FEM/Ferrite.jl/blob/master/src/FEValues/cell_values.jl#L51-L54).

[edljk]

Perhaps I am missing something. Here is a simple example with a circular boundary of order 3 which produces an error.

p = [-0.222521     -0.974928
     4.15312e-17   4.93183e-17
     -0.900969     -0.433884
     -0.148347     -0.649952
     -0.0741736    -0.324976
     -0.300323     -0.144628
     -0.600646     -0.289256
     -0.739526     -0.673128
     -0.491692     -0.870769
     -0.62349      -0.781831
     0.62349      -0.781831
     0.0651862    -0.997873
     0.37423      -0.927336
     0.41566      -0.521221
     0.20783      -0.26061
     0.222521     -0.974928
     1.0           2.03666e-16
     0.820812     -0.571199
     0.958349     -0.2856
     0.666667     -5.52021e-17
     0.333333      3.57829e-17
     0.900969     -0.433884
     -0.222521      0.974928
     -0.900969      0.433884
     -0.491692      0.870769
     -0.739526      0.673128
     -0.600646      0.289256
     -0.300323      0.144628
     -0.0741736     0.324976
     -0.148347      0.649952
     -0.62349       0.781831
     0.62349       0.781831
     0.20783       0.26061
     0.41566       0.521221
     0.37423       0.927336
     0.0651862     0.997873
     0.222521      0.974928
     -0.98736       0.158496
     -0.98736      -0.158496
     -1.0           1.25146e-16
     0.958349      0.2856
     0.820812      0.571199
     0.900969      0.433884]
np = size(p, 1)
nodes = [Ferrite.Node((p[k, :]...,)) for k = 1:np]
cells = [Cell{2, 10, 3}((1, 2, 3, 4, 5, 6, 7, 8, 9, 10)),
         Cell{2, 10, 3}((2, 1, 11, 5, 4, 12, 13, 14, 15, 16)),
         Cell{2, 10, 3}((11, 17, 2, 18, 19, 20, 21, 15, 14, 22)),
         Cell{2, 10, 3}((23, 24, 2, 25, 26, 27, 28, 29, 30, 31)),
         Cell{2, 10, 3}((23, 2, 32, 30, 29, 33, 34, 35, 36, 37)),
         Cell{2, 10, 3}((24, 3, 2, 38, 39, 7, 6, 28, 27, 40)),
         Cell{2, 10, 3}((2, 17, 32, 21, 20, 41, 42, 34, 33, 43))]
grid = Ferrite.Grid(cells, nodes)
@show grid
ip = Lagrange{2, RefTetrahedron, 1}()
qr = QuadratureRule{2, RefTetrahedron}(:legendre, 6)
ip_geom = Lagrange{2, RefTetrahedron, 3}() 
cellvalues = CellScalarValues(qr, ip, ip_geom)
dh = DofHandler(grid)
push!(dh, :u, 1, ip)
close!(dh)

ERROR: LoadError: MethodError: no method matching default_interpolation(::Type{Cell{2, 10, 3}})
Closest candidates are:
  default_interpolation(::Union{Type{Ferrite.Line}, Type{Line2D}, Type{Line3D}}) at .julia/dev/Ferrite/src/Grid/grid.jl:738
  default_interpolation(::Union{Type{Quadrilateral}, Type{Quadrilateral3D}}) at .julia/dev/Ferrite/src/Grid/grid.jl:742
  default_interpolation(::Type{QuadraticLine}) at .julia/dev/Ferrite/src/Grid/grid.jl:739
  ...
Stacktrace:
 [1] push!(dh::DofHandler{2, Float64, Grid{2, Cell{2, 10, 3}, Float64}}, name::Symbol, dim::Int64, ip::Lagrange{2, RefTetrahedron, 1})
   @ Ferrite  .julia/dev/Ferrite/src/Dofs/DofHandler.jl:121
 [2] top-level scope
   @ Julia/experiments/ferrite/test_Ferritecurve.jl:62
 [3] include(fname::String)
   @ Base.MainInclude ./client.jl:451
 [4] top-level scope
   @ REPL[27]:1