Skip to content

Implement Hessians and Identity Mapping for embedded elements #1257

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

Draft
henrij22 wants to merge 11 commits into
base: master
Choose a base branch
from
Draft

Implement Hessians and Identity Mapping for embedded elements #1257

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
0eeecd9
wip first draft
henrij22 Dec 9, 2025
06ac744
vector-valued + more tests
henrij22 Dec 9, 2025
53e0939
more tests + vector-valued
henrij22 Dec 9, 2025
7edb56b
fixes for vector valued
henrij22 Dec 9, 2025
8737dd3
update
henrij22 Dec 10, 2025
831fd7f
Merge remote-tracking branch 'origin/HEAD' into dev/embeddedhessians
henrij22 Dec 10, 2025
158cc2b
fix computation of H in calculate_mapping(gip::ScalarInterpolation)
henrij22 Dec 15, 2025
8011c2f
more performance
henrij22 Dec 15, 2025
6f7fbbd
fix
henrij22 Dec 15, 2025
7260a95
fix2
henrij22 Dec 15, 2025
d3214ec
update
henrij22 Dec 27, 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
1 change: 1 addition & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,3 +24,4 @@ docs/src/changelog.md
LocalPreferences.toml
docs/LocalPreferences.toml
benchmark/tune.json
.vscode/
50 changes: 43 additions & 7 deletions src/FEValues/FunctionValues.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -174,6 +174,39 @@ required_geo_diff_order(::CovariantPiolaMapping, fun_diff_order::Int) = 1 + fun_
# Vector interpolations with sdim > rdim
@inline dothelper(B::SMatrix{vdim, rdim}, A::SMatrix{rdim, sdim}) where {vdim, rdim, sdim} = B * A

# aᵢBᵢⱼₖ = Cⱼₖ
@inline function dothelper(A::SVector{sdim}, B::SArray{Tuple{sdim, rdim, rdim}}) where {rdim, sdim}
return SMatrix{rdim, rdim}((dot(A, B[:, j, k]) for j in 1:rdim, k in 1:rdim))
end

@inline function dothelper(A::SMatrix{vdim, sdim}, B::SArray{Tuple{sdim, rdim, rdim}}) where {vdim, rdim, sdim}
return SArray{Tuple{vdim, rdim, rdim}}(
(dothelper(A[i, :], B)[j, k] for i in 1:vdim, j in 1:rdim, k in 1:rdim)
)
Copy link
Copy Markdown
Collaborator

@lijas lijas Dec 14, 2025

Choose a reason for hiding this comment

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

I think this type of construction of SArrays is fine, but I dont know how efficiently the compiler can optimize it, maybe someone else can review this part :)

Copy link
Copy Markdown
Contributor Author

@henrij22 henrij22 Dec 15, 2025

Choose a reason for hiding this comment

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

You were right, there are a lot of calls in there.
This implementation takes 24 ns on my system, I found a better version which only takes about 2 ns

end

@inline dothelper(A::SMatrix{sdim, rdim}, B::SMatrix{rdim, rdim}) where {rdim, sdim} = A * B
@inline dothelper(A::SMatrix{sdim, rdim}, B::SMatrix{rdim, sdim}) where {rdim, sdim} = A * B

@inline otimesu_helper(A, B) = otimesu(A, B)

# Cᵢⱼₖₗ = AᵢₖBⱼₗ
@inline function otimesu_helper(A::SMatrix{rdim, sdim}, B::SMatrix{rdim, sdim}) where {rdim, sdim}
return SArray{Tuple{rdim, rdim, sdim, sdim}}(
(A[i, k] * B[j, l] for i in 1:rdim, j in 1:rdim, k in 1:sdim, l in 1:sdim)
)
end

@inline dcontract_helper(A, B) = A ⊡ B
@inline dcontract_helper(A::SMatrix{rdim, rdim}, B::SMatrix{rdim, rdim}) where {rdim} = Tensor{2, rdim}(A) ⊡ Tensor{2, rdim}(B)

# Cᵢₗₘ = AᵢⱼₖBⱼₖₗₘ
@inline function dcontract_helper(A::SArray{Tuple{vdim, rdim, rdim}}, B::SArray{Tuple{rdim, rdim, sdim, sdim}}) where {vdim, sdim, rdim}
return SArray{Tuple{vdim, sdim, sdim}}(
(dcontract_helper(A[i, :, :], B[:, :, l, m]) for i in 1:vdim, l in 1:sdim, m in 1:sdim)
)
end

# =============
# Apply mapping
# =============
Expand All @@ -198,18 +231,21 @@ end
@inline function apply_mapping!(funvals::FunctionValues{2}, ::IdentityMapping, q_point::Int, mapping_values, args...)
Jinv = calculate_Jinv(getjacobian(mapping_values))

sdim, rdim = size(Jinv)
(rdim != sdim) && error("apply_mapping! for second order gradients and embedded elements not implemented")

H = gethessian(mapping_values)
is_vector_valued = first(funvals.Nx) isa Vec
Jinv_otimesu_Jinv = is_vector_valued ? otimesu(Jinv, Jinv) : nothing
is_vector_valued = first(funvals.Nx) isa Union{<:Vec, <:SVector}
@inbounds for j in 1:getnbasefunctions(funvals)
dNdx = dothelper(funvals.dNdξ[j, q_point], Jinv)
if is_vector_valued
d2Ndx2 = (funvals.d2Ndξ2[j, q_point] - dNdx ⋅ H) ⊡ Jinv_otimesu_Jinv
t = (funvals.d2Ndξ2[j, q_point] - dothelper(dNdx, H))
Jinv_otimesu_Jinv = otimesu_helper(Jinv, Jinv)
d2Ndx2 = dcontract_helper(t, Jinv_otimesu_Jinv)

# d2Ndx2 = (funvals.d2Ndξ2[j, q_point] - dNdx ⋅ H) ⊡ Jinv_otimesu_Jinv
else
d2Ndx2 = Jinv' ⋅ (funvals.d2Ndξ2[j, q_point] - dNdx ⋅ H) ⋅ Jinv
t = dothelper(Jinv', (funvals.d2Ndξ2[j, q_point] - dothelper(dNdx, H)))
d2Ndx2 = dothelper(t, Jinv)

# d2Ndx2 = Jinv' ⋅ (funvals.d2Ndξ2[j, q_point] - dothelper(dNdx, H)) ⋅ Jinv
end

funvals.dNdx[j, q_point] = dNdx
Expand Down
29 changes: 22 additions & 7 deletions src/FEValues/GeometryMapping.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ end

@inline getjacobian(mv::MappingValues{<:Union{AbstractTensor, SMatrix}}) = mv.J
@inline gethessian(mv::MappingValues{<:Any, <:AbstractTensor}) = mv.H
@inline gethessian(mv::MappingValues{<:SMatrix, <:SArray}) = mv.H


"""
Expand Down Expand Up @@ -112,6 +113,12 @@ geometric_interpolation(geo_mapping::GeometryMapping) = geo_mapping.ip
@inline function otimes_helper(x::Vec{sdim}, dMdξ::Vec{rdim}) where {sdim, rdim}
return SMatrix{sdim, rdim}((x[i] * dMdξ[j] for i in 1:sdim, j in 1:rdim)...)
end

@inline otimes_helper(x::Vec{dim}, d2Mdξ2::Tensor{2, dim}) where {dim} = x ⊗ d2Mdξ2
@inline function otimes_helper(x::Vec{sdim}, d2Mdξ2::Tensor{2, rdim}) where {sdim, rdim}
return SArray{Tuple{sdim, rdim, rdim}}((x[i] * d2Mdξ2[j, k] for i in 1:sdim, j in 1:rdim, k in 1:rdim)...)
end

# End of embedded hot-fixes

# For creating initial value
Expand All @@ -124,6 +131,10 @@ end
function otimes_returntype(#=typeof(x)=# ::Type{<:Vec{dim, Tx}}, #=typeof(d2Mdξ2)=# ::Type{<:Tensor{2, dim, TM}}) where {dim, Tx, TM}
return Tensor{3, dim, promote_type(Tx, TM)}
end
function otimes_returntype(::Type{<:Vec{sdim, Tx}}, #=typeof(d2Mdξ2)=# ::Type{<:Tensor{2, rdim, TM}}) where {sdim, rdim, Tx, TM}
return typeof(StaticArrays.@SArray zeros(promote_type(Tx, TM), sdim, rdim, rdim))
end


@inline function calculate_mapping(::GeometryMapping{0}, q_point::Int, x::AbstractVector{<:Vec})
return MappingValues(nothing, nothing)
Expand All @@ -140,12 +151,13 @@ end

@inline function calculate_mapping(geo_mapping::GeometryMapping{2}, q_point::Int, x::AbstractVector{<:Vec})
J = zero(otimes_returntype(eltype(x), eltype(geo_mapping.dMdξ)))
sdim, rdim = size(J)
(rdim != sdim) && error("hessian for embedded elements not implemented (rdim=$rdim, sdim=$sdim)")
H = zero(otimes_returntype(eltype(x), eltype(geo_mapping.d2Mdξ2)))
@inbounds for j in 1:getngeobasefunctions(geo_mapping)
J += x[j] ⊗ geo_mapping.dMdξ[j, q_point]
H += x[j] ⊗ geo_mapping.d2Mdξ2[j, q_point]
J += otimes_helper(x[j], geo_mapping.dMdξ[j, q_point])
H += otimes_helper(x[j], geo_mapping.d2Mdξ2[j, q_point])

# J += x[j] ⊗ geo_mapping.dMdξ[j, q_point]
# H += x[j] ⊗ geo_mapping.d2Mdξ2[j, q_point]
end
return MappingValues(J, H)
end
Expand All @@ -170,13 +182,16 @@ end
@inline function calculate_mapping(gip::ScalarInterpolation, ξ::Vec{rdim, T}, x::AbstractVector{<:Vec{sdim}}, ::Val{2}) where {T, rdim, sdim}
n_basefuncs = getnbasefunctions(gip)
@boundscheck checkbounds(x, Base.OneTo(n_basefuncs))
(rdim != sdim) && error("hessian for embedded elements not implemented (rdim=$rdim, sdim=$sdim)")
# (rdim != sdim) && error("hessian for embedded elements not implemented (rdim=$rdim, sdim=$sdim)")
Comment thread
henrij22 marked this conversation as resolved.
Outdated
J = zero(otimes_returntype(Vec{sdim, T}, Vec{rdim, T}))
H = zero(otimes_returntype(eltype(x), typeof(J)))
@inbounds for j in 1:n_basefuncs
d2Mdξ2, dMdξ, _ = reference_shape_hessian_gradient_and_value(gip, ξ, j)
J += x[j] ⊗ dMdξ
H += x[j] ⊗ d2Mdξ2

J += otimes_helper(x[j], dMdξ)

# J += x[j] ⊗ dMdξ
# H += x[j] ⊗ d2Mdξ2
Copy link
Copy Markdown
Collaborator

@lijas lijas Dec 14, 2025

Choose a reason for hiding this comment

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

where did the computation of the hessian go?

Copy link
Copy Markdown
Contributor Author

@henrij22 henrij22 Dec 15, 2025

Choose a reason for hiding this comment

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

Aha, I only added it in above. This code then might not be tested atm. Also, what is the difference between calculate_mapping(gip::ScalarInterpolation) and calculate_mapping(geo_mapping::GeometryMapping). They seem quite similar.

end
return MappingValues(J, H)
end
Expand Down
80 changes: 77 additions & 3 deletions test/test_cellvalues.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -441,9 +441,83 @@
cv_vector = CellValues(qr, ip^2, ip^3; update_hessians = true)
cv_scalar = CellValues(qr, ip, ip^3; update_hessians = true)

coords = [Vec{3}((x[1], x[2], 0.0)) for x in Ferrite.reference_coordinates(ip)]
@test_throws ErrorException reinit!(cv_vector, coords) # Not implemented for embedded elements
@test_throws ErrorException reinit!(cv_scalar, coords)
coords = [Vec{3}((2x[1], 0.5x[2], rand())) for x in Ferrite.reference_coordinates(ip)]

# Test Scalar H
reinit!(cv_scalar, coords)

qp = 1
H = Ferrite.calculate_mapping(cv_scalar.geo_mapping, qp, coords).H

d2Ndξ2 = cv_scalar.fun_values.d2Ndξ2
∂A₁∂₁ = Vec{3}(i -> getindex.(d2Ndξ2[:, qp], 1, 1) ⋅ getindex.(coords, i))
∂A₂∂₂ = Vec{3}(i -> getindex.(d2Ndξ2[:, qp], 2, 2) ⋅ getindex.(coords, i))
∂A₁∂₂ = Vec{3}(i -> getindex.(d2Ndξ2[:, qp], 1, 2) ⋅ getindex.(coords, i)) # = ∂A₂∂₁
∂A₂∂₁ = Vec{3}(i -> getindex.(d2Ndξ2[:, qp], 2, 1) ⋅ getindex.(coords, i))

@test H[:, 1, 1] ≈ ∂A₁∂₁
@test H[:, 2, 2] ≈ ∂A₂∂₂
@test H[:, 1, 2] ≈ ∂A₁∂₂
@test H[:, 2, 1] ≈ ∂A₂∂₁


# Test Vector H
reinit!(cv_vector, coords)

H = Ferrite.calculate_mapping(cv_vector.geo_mapping, qp, coords).H

d2Ndξ2 = cv_vector.fun_values.d2Ndξ2
∂A₁∂₁ = Vec{3}(i -> getindex.(d2Ndξ2[1:2:18, qp], 1, 1, 1) ⋅ getindex.(coords, i))
∂A₂∂₂ = Vec{3}(i -> getindex.(d2Ndξ2[1:2:18, qp], 1, 2, 2) ⋅ getindex.(coords, i))
∂A₁∂₂ = Vec{3}(i -> getindex.(d2Ndξ2[1:2:18, qp], 1, 1, 2) ⋅ getindex.(coords, i)) # = ∂A₂∂₁
∂A₂∂₁ = Vec{3}(i -> getindex.(d2Ndξ2[1:2:18, qp], 1, 2, 1) ⋅ getindex.(coords, i))


@test H[:, 1, 1] ≈ ∂A₁∂₁
@test H[:, 2, 2] ≈ ∂A₂∂₂
@test H[:, 1, 2] ≈ ∂A₁∂₂
@test H[:, 2, 1] ≈ ∂A₂∂₁

# Test Mapping Scalar
coords_scaled = [Vec{3}((2x[1], 0.5x[2], 0.0)) for x in Ferrite.reference_coordinates(ip)]
reinit!(cv_scalar, coords_scaled)

scale_x = 2.0
scale_y = 0.5

coords_ref = [Vec{3}((x[1], x[2], 0.0)) for x in Ferrite.reference_coordinates(ip)]
cv_ref = CellValues(qr, ip, ip^3; update_hessians = true)
reinit!(cv_ref, coords_ref)

# Test scaling embedded vs embedded
@test shape_hessian(cv_scalar, qp, 1)[1, 1] * scale_x^2 ≈ shape_hessian(cv_ref, qp, 1)[1, 1]
@test shape_hessian(cv_scalar, qp, 1)[2, 2] * scale_y^2 ≈ shape_hessian(cv_ref, qp, 1)[2, 2]
@test shape_hessian(cv_scalar, qp, 1)[3, 3] ≈ shape_hessian(cv_ref, qp, 1)[3, 3]
@test shape_hessian(cv_scalar, qp, 1)[1, 2] * scale_x * scale_y ≈ shape_hessian(cv_ref, qp, 1)[1, 2]
@test shape_hessian(cv_scalar, qp, 1)[2, 1] * scale_x * scale_y ≈ shape_hessian(cv_ref, qp, 1)[2, 1]

# Test Mapping Vector
cv_ref2 = CellValues(qr, ip^2, ip^3; update_hessians = true)
reinit!(cv_vector, coords_scaled)
reinit!(cv_ref2, coords_ref)

# Test scaling embedded vs embedded
@test shape_hessian(cv_vector, qp, 1)[1, 1, 1] * scale_x^2 ≈ shape_hessian(cv_ref2, qp, 1)[1, 1, 1]
@test shape_hessian(cv_vector, qp, 1)[1, 2, 2] * scale_y^2 ≈ shape_hessian(cv_ref2, qp, 1)[1, 2, 2]
@test shape_hessian(cv_vector, qp, 1)[1, 3, 3] ≈ shape_hessian(cv_ref2, qp, 1)[1, 3, 3]
@test shape_hessian(cv_vector, qp, 1)[1, 1, 2] * scale_x * scale_y ≈ shape_hessian(cv_ref2, qp, 1)[1, 1, 2]
@test shape_hessian(cv_vector, qp, 1)[1, 2, 1] * scale_x * scale_y ≈ shape_hessian(cv_ref2, qp, 1)[1, 2, 1]


# Test embedded vs non-embedded
cv_ref2n = CellValues(qr, ip^2, ip^2; update_hessians = true)
coords_refn = [Vec{2}((x[1], x[2])) for x in Ferrite.reference_coordinates(ip)]
reinit!(cv_ref2n, coords_refn)

@test shape_hessian(cv_ref2, qp, 1)[1, 1, 1] ≈ shape_hessian(cv_ref2n, qp, 1)[1, 1, 1]
@test shape_hessian(cv_ref2, qp, 1)[1, 2, 1] ≈ shape_hessian(cv_ref2n, qp, 1)[1, 2, 1]
@test shape_hessian(cv_ref2, qp, 1)[2, 1, 1] ≈ shape_hessian(cv_ref2n, qp, 1)[2, 1, 1]
@test shape_hessian(cv_ref2, qp, 1)[2, 2, 1] ≈ shape_hessian(cv_ref2n, qp, 1)[2, 2, 1]
Copy link
Copy Markdown
Collaborator

@lijas lijas Dec 14, 2025 •
edited
Loading

Choose a reason for hiding this comment

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

I think we can add tests for the function_hessian also, as we do here:

                coords_nl = [x + rand(x) * 0.01 for x in coords] # add some displacement to nodes
                reinit!(cv, coords_nl)

                _ue_nl = [u_funk(coords_nl[i], V, G, H) for i in 1:n_basefunc_base]
                ue_nl = reinterpret(Float64, _ue_nl)

                for i in 1:getnquadpoints(cv)
                    xqp = spatial_coordinate(cv, i, coords_nl)
                    Hqp, Gqp, Vqp = Tensors.hessian(x -> function_value_from_physical_coord(func_interpol, coords_nl, x, ue_nl), xqp, :all)
                    @test function_value(cv, i, ue_nl) ≈ Vqp
                    @test function_gradient(cv, i, ue_nl) ≈ Gqp
                    if update_hessians
                        @test Ferrite.function_hessian(cv, i, ue_nl) ≈ Hqp
                    end
                end

Can you add a similar test in this @testset for hessians? (I am not sure if function_value_from_physical_coord works with embedded elements)

Copy link
Copy Markdown
Collaborator

@lijas lijas Dec 15, 2025

Choose a reason for hiding this comment

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

This might be more difficult than i first thought. For the embedded case, the function_value_from_physical_coord would need some signed distance search to find the parametric coord on the surface

end
end

Expand Down
Loading