Skip to content

Add Lagrange{RefTetrahedron,3} #1255

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
termi-official wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Add Lagrange{RefTetrahedron,3} #1255

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
77a5ee8
Add Lagrange{RefTetrahedron,3}
termi-official Nov 26, 2025
23b31ed
Runic
termi-official Nov 26, 2025
ca070ed
Fix oopsie.
termi-official Nov 26, 2025
b01d1b6
Fix CI warning
termi-official Nov 26, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
141 changes: 141 additions & 0 deletions src/interpolations.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ The following interpolations are implemented:
* `Lagrange{RefHexahedron, 2}`
* `Lagrange{RefTetrahedron, 1}`
* `Lagrange{RefTetrahedron, 2}`
* `Lagrange{RefTetrahedron, 3}`
* `Lagrange{RefPrism, 1}`
* `Lagrange{RefPrism, 2}`
* `Lagrange{RefPyramid, 1}`
Expand Down Expand Up @@ -932,6 +933,80 @@ function reference_shape_value(ip::Lagrange{RefTetrahedron, 2}, ξ::Vec{3}, i::I
throw(ArgumentError("no shape function $i for interpolation $ip"))
end

###################################
# Lagrange RefTetrahedron order 3 #
###################################
getnbasefunctions(::Lagrange{RefTetrahedron, 3}) = 20
Copy link
Collaborator

@AbdAlazezAhmed AbdAlazezAhmed Dec 1, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe

Suggested change
getnbasefunctions(::Lagrange{RefTetrahedron, 3}) = 20
getnbasefunctions(::Lagrange{RefTetrahedron, order}) where order = (order+1)*(order+2)*(order+3)÷6

would be beneficial in the future


facedof_indices(::Lagrange{RefTetrahedron, 3}) = ((1, 3, 2, 10, 9, 8, 7, 6, 5, 17), (1, 2, 4, 5, 6, 13, 14, 11, 12, 18), (2, 3, 4, 7, 8, 15, 16, 13, 14, 19), (1, 4, 3, 11, 12, 15, 16, 9, 10, 20))
facedof_interior_indices(::Lagrange{RefTetrahedron, 3}) = ((17,), (18,), (19,), (20,))
edgedof_indices(::Lagrange{RefTetrahedron, 3}) = ((1, 2, 5, 6), (2, 3, 7, 8), (3, 1, 9, 10), (1, 4, 11, 12), (2, 4, 13, 14), (3, 4, 15, 16))
edgedof_interior_indices(::Lagrange{RefTetrahedron, 3}) = ((5, 6), (7, 8), (9, 10), (11, 12), (13, 14), (15, 16))

function reference_coordinates(::Lagrange{RefTetrahedron, 3})
return [
# Vertices
Vec{3, Float64}((0.0, 0.0, 0.0)),
Vec{3, Float64}((1.0, 0.0, 0.0)),
Vec{3, Float64}((0.0, 1.0, 0.0)),
Vec{3, Float64}((0.0, 0.0, 1.0)),
# Edges
#1
Vec{3, Float64}((1 / 3, 0.0, 0.0)),
Vec{3, Float64}((2 / 3, 0.0, 0.0)),
#2
Vec{3, Float64}((2 / 3, 1 / 3, 0.0)),
Vec{3, Float64}((1 / 3, 2 / 3, 0.0)),
#3
Vec{3, Float64}((0.0, 2 / 3, 0.0)),
Vec{3, Float64}((0.0, 1 / 3, 0.0)),
#4
Vec{3, Float64}((0.0, 0.0, 1 / 3)),
Vec{3, Float64}((0.0, 0.0, 2 / 3)),
#5
Vec{3, Float64}((2 / 3, 0.0, 1 / 3)),
Vec{3, Float64}((1 / 3, 0.0, 2 / 3)),
#6
Vec{3, Float64}((0.0, 2 / 3, 1 / 3)),
Vec{3, Float64}((0.0, 1 / 3, 2 / 3)),
# Faces
Vec{3, Float64}((1 / 3, 1 / 3, 0.0)),
Vec{3, Float64}((1 / 3, 0.0, 1 / 3)),
Vec{3, Float64}((1 / 3, 1 / 3, 1 / 3)),
Vec{3, Float64}((0.0, 1 / 3, 1 / 3)),
]
end

# http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch09.d/AFEM.Ch09.pdf
# http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch10.d/AFEM.Ch10.pdf
function reference_shape_value(ip::Lagrange{RefTetrahedron, 3}, ξ::Vec{3, T}, i::Int) where {T}
ξ_x = ξ[1]
ξ_y = ξ[2]
ξ_z = ξ[3]
# TODO implement as tensor product form after https://github.com/Ferrite-FEM/Ferrite.jl/pull/1188 is merged
i == 1 && return T(1.0) + T(-5.5) * ξ_x + T(9.0) * ξ_x^2 + T(-4.5) * ξ_x^3 + T(-5.5) * ξ_y + T(18.0) * ξ_x * ξ_y + T(-13.5) * ξ_x^2 * ξ_y + T(9.0) * ξ_y^2 + T(-13.5) * ξ_x * ξ_y^2 + T(-4.5) * ξ_y^3 + T(-5.5) * ξ_z + T(18.0) * ξ_x * ξ_z + T(-13.5) * ξ_x^2 * ξ_z + T(18.0) * ξ_y * ξ_z + T(-27.0) * ξ_x * ξ_y * ξ_z + T(-13.5) * ξ_y^2 * ξ_z + T(9.0) * ξ_z^2 + T(-13.5) * ξ_x * ξ_z^2 + T(-13.5) * ξ_y * ξ_z^2 + T(-4.5) * ξ_z^3
i == 2 && return ξ_x + T(-4.5) * ξ_x^2 + T(4.5) * ξ_x^3
i == 3 && return ξ_y + T(-4.5) * ξ_y^2 + T(4.5) * ξ_y^3
i == 4 && return ξ_z + T(-4.5) * ξ_z^2 + T(4.5) * ξ_z^3
i == 5 && return T(9.0) * ξ_x + T(-22.5) * ξ_x^2 + T(13.5) * ξ_x^3 + T(-22.5) * ξ_x * ξ_y + T(27.0) * ξ_x^2 * ξ_y + T(13.5) * ξ_x * ξ_y^2 + T(-22.5) * ξ_x * ξ_z + T(27.0) * ξ_x^2 * ξ_z + T(27.0) * ξ_x * ξ_y * ξ_z + T(13.5) * ξ_x * ξ_z^2
i == 6 && return T(-4.5) * ξ_x + T(18.0) * ξ_x^2 + T(-13.5) * ξ_x^3 + T(4.5) * ξ_x * ξ_y + T(-13.5) * ξ_x^2 * ξ_y + T(4.5) * ξ_x * ξ_z + T(-13.5) * ξ_x^2 * ξ_z
i == 7 && return T(-4.5) * ξ_x * ξ_y + T(13.5) * ξ_x^2 * ξ_y
i == 8 && return T(-4.5) * ξ_x * ξ_y + T(13.5) * ξ_x * ξ_y^2
i == 9 && return T(-4.5) * ξ_y + T(4.5) * ξ_x * ξ_y + T(18) * ξ_y^2 + T(-13.5) * ξ_x * ξ_y^2 + T(-13.5) * ξ_y^3 + T(4.5) * ξ_y * ξ_z + T(-13.5) * ξ_y^2 * ξ_z
i == 10 && return T(9.0) * ξ_y + T(-22.5) * ξ_x * ξ_y + T(13.5) * ξ_x^2 * ξ_y + T(-22.5) * ξ_y^2 + T(27.0) * ξ_x * ξ_y^2 + T(13.5) * ξ_y^3 + T(-22.5) * ξ_y * ξ_z + T(27.0) * ξ_x * ξ_y * ξ_z + T(27.0) * ξ_y^2 * ξ_z + T(13.5) * ξ_y * ξ_z^2
i == 11 && return T(9.0) * ξ_z + T(-22.5) * ξ_x * ξ_z + T(13.5) * ξ_x^2 * ξ_z + T(-22.5) * ξ_y * ξ_z + T(27.0) * ξ_x * ξ_y * ξ_z + T(13.5) * ξ_y^2 * ξ_z + T(-22.5) * ξ_z^2 + T(27.0) * ξ_x * ξ_z^2 + T(27.0) * ξ_y * ξ_z^2 + T(13.5) * ξ_z^3
i == 12 && return T(-4.5) * ξ_z + T(4.5) * ξ_x * ξ_z + T(4.5) * ξ_y * ξ_z + T(18.0) * ξ_z^2 + T(-13.5) * ξ_x * ξ_z^2 + T(-13.5) * ξ_y * ξ_z^2 + T(-13.5) * ξ_z^3
i == 13 && return T(-4.5) * ξ_x * ξ_z + T(13.5) * ξ_x^2 * ξ_z
i == 14 && return T(-4.5) * ξ_x * ξ_z + T(13.5) * ξ_x * ξ_z^2
i == 15 && return T(-4.5) * ξ_y * ξ_z + T(13.5) * ξ_y^2 * ξ_z
i == 16 && return T(-4.5) * ξ_y * ξ_z + T(13.5) * ξ_y * ξ_z^2
i == 17 && return T(27.0) * ξ_x * ξ_y + T(-27.0) * ξ_x^2 * ξ_y + T(-27.0) * ξ_x * ξ_y^2 + T(-27.0) * ξ_x * ξ_y * ξ_z
i == 18 && return T(27.0) * ξ_x * ξ_z + T(-27.0) * ξ_x^2 * ξ_z + T(-27.0) * ξ_x * ξ_y * ξ_z + T(-27.0) * ξ_x * ξ_z^2
i == 19 && return T(27.0) * ξ_x * ξ_y * ξ_z
i == 20 && return T(27.0) * ξ_y * ξ_z + T(-27.0) * ξ_x * ξ_y * ξ_z + T(-27.0) * ξ_y^2 * ξ_z + T(-27.0) * ξ_y * ξ_z^2
throw(ArgumentError("no shape function $i for interpolation $ip"))
end

##################################
# Lagrange RefHexahedron order 1 #
##################################
Expand Down Expand Up @@ -2132,3 +2207,69 @@ adjust_dofs_during_distribution(::Nedelec{RefHexahedron, 1}) = false
function get_direction(::Nedelec{RefHexahedron, 1}, shape_nr, cell)
return get_edge_direction(cell, shape_nr) # shape_nr = edge_nr
end

# TODO where to put this? I am a lazy man and I won't compute this by hand.
Copy link
Member

@KnutAM KnutAM Dec 1, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems like this should be a macro or generated function?

Copy link
Member Author

@termi-official termi-official Dec 1, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not 100% sure yet. This is for now just a template to get things started on higher order element generators, as we can simply generate these in a systematic way.

Copy link
Collaborator

@AbdAlazezAhmed AbdAlazezAhmed Dec 1, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe

@generated function do_magic(::Lagrange{RefTetrahedron, order}, ξ::Vec{3, Float64}, i::Int) where {order}

    N = 0
    for p₃ in 0:order
        for p₂ in 0:(order - p₃)
            for p₁ in 0:(order - p₃ - p₂)
                N += 1
            end
        end
    end

    @assert N == 20
    xs = reference_coordinates(Lagrange{RefTetrahedron, order})
    V = zeros(N, N)
    j = 0
    for p₃ in 0:order
        for p₂ in 0:(order - p₃)
            for p₁ in 0:(order - p₃ - p₂)
                j += 1
                for i in 1:N
                    x = xs[i]
                    V[i, j] = x[1]^p₁ * x[2]^p₂ * x[3]^p₃
                end
            end
        end
    end
    C = inv(V)' # These are our coefficients
    ex = :(    ξ_x = ξ[1];
    ξ_y = ξ[2];
    ξ_z = ξ[3];)
    for i in 1:N
        j = 0
        x = xs[i]
        ex_i = :(0.0)
        for p₃ in 0:order
            for p₂ in 0:(order - p₃)
                for p₁ in 0:(order - p₃ - p₂)
                    j += 1
                    ex_i = :($ex_i + $(C[i, j]) * ξ_x^$p₁ * ξ_y^$p₂ * ξ_z^$p₃)
                end
            end
        end
        ex = :($ex; i == $i && return $ex_i)
    end

    return ex
end
julia> @benchmark Ferrite.do_magic($ip, $(Vec(-0.0,-0.0,-0.0)),$(1))
BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
 Range (min … max):  4.882 ns … 20.719 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     5.008 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   5.030 ns ±  0.389 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▅     █▂                                                   
  ▅█▆▃▂▂▄██▄▂▂▂▂▁▂▂▂▁▂▄▄▃▂▂▂▂▂▂▂▂▁▁▂▂▂▂▁▁▁▂▂▂▂▁▁▂▂▁▁▁▁▁▂▂▁▁▂ ▃
  4.88 ns        Histogram: frequency by time        5.86 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark Ferrite.reference_shape_value($ip, $(Vec(-0.0,-0.0,-0.0)),$(1))
BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
 Range (min … max):  4.581 ns … 9.742 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     4.619 ns             ┊ GC (median):    0.00%
 Time  (mean ± σ):   4.673 ns ± 0.209 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▅█                                                        
  ▂██▂▂▂▂▄▄▂▂▂▂▂▂▂▂▂▂▂▄▂▂▂▂▂▂▂▁▂▂▂▁▂▂▂▂▁▁▁▁▁▁▁▂▁▂▂▂▂▂▁▁▁▁▁▂ ▂
  4.58 ns        Histogram: frequency by time       5.53 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

compiler vroom

function tabulate_lagrange(interpolation::Lagrange{RefTetrahedron, order}; tol = 1.0e-8, with_T = false) where {order}
N = 0
for p₃ in 0:order
for p₂ in 0:(order - p₃)
for p₁ in 0:(order - p₃ - p₂)
N += 1
end
end
end

@assert N == getnbasefunctions(interpolation)
xs = reference_coordinates(interpolation)
V = zeros(N, N)
j = 0
for p₃ in 0:order
for p₂ in 0:(order - p₃)
for p₁ in 0:(order - p₃ - p₂)
j += 1
for i in 1:N
x = xs[i]
V[i, j] = x[1]^p₁ * x[2]^p₂ * x[3]^p₃
end
end
end
end
C = inv(V)' # These are our coefficients

for i in 1:N
j = 0
x = xs[i]
isfirst = true
for p₃ in 0:order
for p₂ in 0:(order - p₃)
for p₁ in 0:(order - p₃ - p₂)
j += 1
(abs(C[i, j]) < tol) && continue
if isfirst
print("i == $i && return")
else
print(" +")
end
isfirst = false
if with_T
print(" T($(C[i, j]))")
else
print(" $(C[i, j])")
end
if p₁ > 0
p₁ > 1 ? print(" * ξ_x^$p₁") : print(" * ξ_x")
end
if p₂ > 0
p₂ > 1 ? print(" * ξ_y^$p₂") : print(" * ξ_y")
end
if p₃ > 0
p₃ > 1 ? print(" * ξ_z^$p₃") : print(" * ξ_z")
end
end
end
end
println()
end

return C
end
3 changes: 2 additions & 1 deletion test/integration/test_simple_scalar_convergence.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ module ConvergenceTestHelper
get_geometry(::Ferrite.Interpolation{RefPyramid}) = Pyramid

get_quadrature_order(::Lagrange{shape, order}) where {shape, order} = max(2 * order - 1, 2)
get_quadrature_order(::Lagrange{RefTriangle, 5}) where {shape, order} = 8
get_quadrature_order(::Lagrange{RefTriangle, 5}) = 8
get_quadrature_order(::Lagrange{RefPrism, order}) where {order} = 2 * order # Don't know why
get_quadrature_order(::Serendipity{shape, order}) where {shape, order} = max(2 * order - 1, 2)
get_quadrature_order(::CrouzeixRaviart{shape, order}) where {shape, order} = max(2 * order - 1, 2)
Expand Down Expand Up @@ -195,6 +195,7 @@ end
Lagrange{RefTriangle, 2}(),
Lagrange{RefHexahedron, 2}(),
Lagrange{RefTetrahedron, 2}(),
Lagrange{RefTetrahedron, 3}(),
Lagrange{RefPrism, 2}(),
CrouzeixRaviart{RefTriangle, 1}(),
CrouzeixRaviart{RefTetrahedron, 1}(),
Expand Down
6 changes: 4 additions & 2 deletions test/test_interpolations.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,6 +157,7 @@ end
Serendipity{RefHexahedron, 2}(),
Lagrange{RefTetrahedron, 1}(),
Lagrange{RefTetrahedron, 2}(),
Lagrange{RefTetrahedron, 3}(),
Lagrange{RefPrism, 1}(),
Lagrange{RefPrism, 2}(),
Lagrange{RefPyramid, 1}(),
Expand Down Expand Up @@ -219,14 +220,15 @@ end
end
end

# Check for dirac delta property of interpolation
# Check for dirac delta property of interpolatio
@testset "dirac delta property of dof $dof" for dof in 1:n_basefuncs
for k in 1:n_basefuncs
N_dof = reference_shape_value(interpolation, coords[dof], k)
if k == dof
@test N_dof ≈ 1.0
else
factor = interpolation isa Lagrange{RefQuadrilateral, 3} ? 200 : 4
# High order elements are right now not implemented in factorized form, so there is some small numerical noise in the evaluation.
factor = Ferrite.getorder(interpolation) > 2 ? 200 : 4
@test N_dof ≈ 0.0 atol = factor * eps(typeof(N_dof))
end
end
Expand Down
Loading