Skip to content

Hessians 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

Open
henrij22 wants to merge 10 commits into
base: master
Choose a base branch
from
+177 −51
Open

Hessians for embedded elements #1257

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 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
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/
118 changes: 77 additions & 41 deletions src/FEValues/FunctionValues.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,33 +6,33 @@
#################################################################

# Scalar, sdim == rdim sdim rdim
typeof_N(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim, <:AbstractRefShape{dim}}) where {T, dim} = T
typeof_dNdx(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim, <:AbstractRefShape{dim}}) where {T, dim} = Vec{dim, T}
typeof_dNdξ(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim, <:AbstractRefShape{dim}}) where {T, dim} = Vec{dim, T}
typeof_d2Ndx2(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim, <:AbstractRefShape{dim}}) where {T, dim} = Tensor{2, dim, T}
typeof_d2Ndξ2(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim, <:AbstractRefShape{dim}}) where {T, dim} = Tensor{2, dim, T}
typeof_N(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim,<:AbstractRefShape{dim}}) where {T,dim} = T
typeof_dNdx(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim,<:AbstractRefShape{dim}}) where {T,dim} = Vec{dim,T}
typeof_dNdξ(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim,<:AbstractRefShape{dim}}) where {T,dim} = Vec{dim,T}
typeof_d2Ndx2(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim,<:AbstractRefShape{dim}}) where {T,dim} = Tensor{2,dim,T}
typeof_d2Ndξ2(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim,<:AbstractRefShape{dim}}) where {T,dim} = Tensor{2,dim,T}

# Vector, vdim == sdim == rdim vdim sdim rdim
typeof_N(::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim, <:AbstractRefShape{dim}}) where {T, dim} = Vec{dim, T}
typeof_dNdx(::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim, <:AbstractRefShape{dim}}) where {T, dim} = Tensor{2, dim, T}
typeof_dNdξ(::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim, <:AbstractRefShape{dim}}) where {T, dim} = Tensor{2, dim, T}
typeof_d2Ndx2(::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim, <:AbstractRefShape{dim}}) where {T, dim} = Tensor{3, dim, T}
typeof_d2Ndξ2(::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim, <:AbstractRefShape{dim}}) where {T, dim} = Tensor{3, dim, T}
typeof_N(::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim,<:AbstractRefShape{dim}}) where {T,dim} = Vec{dim,T}
typeof_dNdx(::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim,<:AbstractRefShape{dim}}) where {T,dim} = Tensor{2,dim,T}
typeof_dNdξ(::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim,<:AbstractRefShape{dim}}) where {T,dim} = Tensor{2,dim,T}
typeof_d2Ndx2(::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim,<:AbstractRefShape{dim}}) where {T,dim} = Tensor{3,dim,T}
typeof_d2Ndξ2(::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim,<:AbstractRefShape{dim}}) where {T,dim} = Tensor{3,dim,T}

# Scalar, sdim != rdim (TODO: Use Vec if (s|r)dim <= 3?)
typeof_N(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{sdim, <:AbstractRefShape{rdim}}) where {T, sdim, rdim} = T
typeof_dNdx(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{sdim, <:AbstractRefShape{rdim}}) where {T, sdim, rdim} = SVector{sdim, T}
typeof_dNdξ(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{sdim, <:AbstractRefShape{rdim}}) where {T, sdim, rdim} = SVector{rdim, T}
typeof_d2Ndx2(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{sdim, <:AbstractRefShape{rdim}}) where {T, sdim, rdim} = SMatrix{sdim, sdim, T, sdim * sdim}
typeof_d2Ndξ2(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{sdim, <:AbstractRefShape{rdim}}) where {T, sdim, rdim} = SMatrix{rdim, rdim, T, rdim * rdim}
typeof_N(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{sdim,<:AbstractRefShape{rdim}}) where {T,sdim,rdim} = T
typeof_dNdx(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{sdim,<:AbstractRefShape{rdim}}) where {T,sdim,rdim} = SVector{sdim,T}
typeof_dNdξ(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{sdim,<:AbstractRefShape{rdim}}) where {T,sdim,rdim} = SVector{rdim,T}
typeof_d2Ndx2(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{sdim,<:AbstractRefShape{rdim}}) where {T,sdim,rdim} = SMatrix{sdim,sdim,T,sdim * sdim}
typeof_d2Ndξ2(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{sdim,<:AbstractRefShape{rdim}}) where {T,sdim,rdim} = SMatrix{rdim,rdim,T,rdim * rdim}


# Vector, vdim != sdim != rdim (TODO: Use Vec/Tensor if (s|r)dim <= 3?)
typeof_N(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim, <:AbstractRefShape{rdim}}) where {T, vdim, sdim, rdim} = SVector{vdim, T}
typeof_dNdx(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim, <:AbstractRefShape{rdim}}) where {T, vdim, sdim, rdim} = SMatrix{vdim, sdim, T, vdim * sdim}
typeof_dNdξ(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim, <:AbstractRefShape{rdim}}) where {T, vdim, sdim, rdim} = SMatrix{vdim, rdim, T, vdim * rdim}
typeof_d2Ndx2(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim, <:AbstractRefShape{rdim}}) where {T, vdim, sdim, rdim} = SArray{Tuple{vdim, sdim, sdim}, T, 3, vdim * sdim * sdim}
typeof_d2Ndξ2(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim, <:AbstractRefShape{rdim}}) where {T, vdim, sdim, rdim} = SArray{Tuple{vdim, rdim, rdim}, T, 3, vdim * rdim * rdim}
typeof_N(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim,<:AbstractRefShape{rdim}}) where {T,vdim,sdim,rdim} = SVector{vdim,T}
typeof_dNdx(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim,<:AbstractRefShape{rdim}}) where {T,vdim,sdim,rdim} = SMatrix{vdim,sdim,T,vdim * sdim}
typeof_dNdξ(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim,<:AbstractRefShape{rdim}}) where {T,vdim,sdim,rdim} = SMatrix{vdim,rdim,T,vdim * rdim}
typeof_d2Ndx2(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim,<:AbstractRefShape{rdim}}) where {T,vdim,sdim,rdim} = SArray{Tuple{vdim,sdim,sdim},T,3,vdim * sdim * sdim}
typeof_d2Ndξ2(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim,<:AbstractRefShape{rdim}}) where {T,vdim,sdim,rdim} = SArray{Tuple{vdim,rdim,rdim},T,3,vdim * rdim * rdim}


"""
Expand All @@ -43,25 +43,25 @@ for both the reference cell (precalculated) and the real cell (updated in `reini
"""
FunctionValues

struct FunctionValues{DiffOrder, IP, N_t, dNdx_t, dNdξ_t, d2Ndx2_t, d2Ndξ2_t}
struct FunctionValues{DiffOrder,IP,N_t,dNdx_t,dNdξ_t,d2Ndx2_t,d2Ndξ2_t}
ip::IP # ::Interpolation
Nx::N_t # ::AbstractMatrix{Union{<:Tensor,<:Number}}
Nξ::N_t # ::AbstractMatrix{Union{<:Tensor,<:Number}}
dNdx::dNdx_t # ::AbstractMatrix{Union{<:Tensor,<:StaticArray}} or Nothing
dNdξ::dNdξ_t # ::AbstractMatrix{Union{<:Tensor,<:StaticArray}} or Nothing
d2Ndx2::d2Ndx2_t # ::AbstractMatrix{<:Tensor{2}} Hessians of geometric shape functions in ref-domain
d2Ndξ2::d2Ndξ2_t # ::AbstractMatrix{<:Tensor{2}} Hessians of geometric shape functions in ref-domain
function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, ::Nothing, ::Nothing, ::Nothing, ::Nothing) where {N_t <: AbstractMatrix}
return new{0, typeof(ip), N_t, Nothing, Nothing, Nothing, Nothing}(ip, Nx, Nξ, nothing, nothing, nothing, nothing)
function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, ::Nothing, ::Nothing, ::Nothing, ::Nothing) where {N_t<:AbstractMatrix}
return new{0,typeof(ip),N_t,Nothing,Nothing,Nothing,Nothing}(ip, Nx, Nξ, nothing, nothing, nothing, nothing)
end
function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, dNdx::AbstractMatrix, dNdξ::AbstractMatrix, ::Nothing, ::Nothing) where {N_t <: AbstractMatrix}
return new{1, typeof(ip), N_t, typeof(dNdx), typeof(dNdξ), Nothing, Nothing}(ip, Nx, Nξ, dNdx, dNdξ, nothing, nothing)
function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, dNdx::AbstractMatrix, dNdξ::AbstractMatrix, ::Nothing, ::Nothing) where {N_t<:AbstractMatrix}
return new{1,typeof(ip),N_t,typeof(dNdx),typeof(dNdξ),Nothing,Nothing}(ip, Nx, Nξ, dNdx, dNdξ, nothing, nothing)
end
function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, dNdx::AbstractMatrix, dNdξ::AbstractMatrix, d2Ndx2::AbstractMatrix, d2Ndξ2::AbstractMatrix) where {N_t <: AbstractMatrix}
return new{2, typeof(ip), N_t, typeof(dNdx), typeof(dNdξ), typeof(d2Ndx2), typeof(d2Ndξ2)}(ip, Nx, Nξ, dNdx, dNdξ, d2Ndx2, d2Ndξ2)
function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, dNdx::AbstractMatrix, dNdξ::AbstractMatrix, d2Ndx2::AbstractMatrix, d2Ndξ2::AbstractMatrix) where {N_t<:AbstractMatrix}
return new{2,typeof(ip),N_t,typeof(dNdx),typeof(dNdξ),typeof(d2Ndx2),typeof(d2Ndξ2)}(ip, Nx, Nξ, dNdx, dNdξ, d2Ndx2, d2Ndξ2)
end
end
function FunctionValues{DiffOrder}(::Type{T}, ip::Interpolation, qr::QuadratureRule, ip_geo::VectorizedInterpolation) where {DiffOrder, T}
function FunctionValues{DiffOrder}(::Type{T}, ip::Interpolation, qr::QuadratureRule, ip_geo::VectorizedInterpolation) where {DiffOrder,T}
assert_same_refshapes(qr, ip, ip_geo)
n_shape = getnbasefunctions(ip)
n_qpoints = getnquadpoints(qr)
Expand Down Expand Up @@ -126,9 +126,9 @@ shape_hessian_type(::FunctionValues{1}) = nothing


# Checks that the user provides the right dimension of coordinates to reinit! methods to ensure good error messages if not
sdim_from_gradtype(::Type{<:AbstractTensor{<:Any, sdim}}) where {sdim} = sdim
sdim_from_gradtype(::Type{<:AbstractTensor{<:Any,sdim}}) where {sdim} = sdim
sdim_from_gradtype(::Type{<:SVector{sdim}}) where {sdim} = sdim
sdim_from_gradtype(::Type{<:SMatrix{<:Any, sdim}}) where {sdim} = sdim
sdim_from_gradtype(::Type{<:SMatrix{<:Any,sdim}}) where {sdim} = sdim

# For performance, these must be fully inferable for the compiler.
# args: valname (:CellValues or :FacetValues), shape_gradient_type, eltype(x)
Expand All @@ -138,7 +138,7 @@ function check_reinit_sdim_consistency(valname, gradtype::Type, ::Type{<:Vec{sdi
end
check_reinit_sdim_consistency(_, ::Nothing, ::Type{<:Vec}) = nothing # gradient not stored, cannot check
check_reinit_sdim_consistency(_, ::Val{sdim}, ::Val{sdim}) where {sdim} = nothing
function check_reinit_sdim_consistency(valname, ::Val{sdim_val}, ::Val{sdim_x}) where {sdim_val, sdim_x}
function check_reinit_sdim_consistency(valname, ::Val{sdim_val}, ::Val{sdim_x}) where {sdim_val,sdim_x}
throw(ArgumentError("The $valname (sdim=$sdim_val) and coordinates (sdim=$sdim_x) have different spatial dimensions."))
end

Expand Down Expand Up @@ -168,11 +168,44 @@ required_geo_diff_order(::CovariantPiolaMapping, fun_diff_order::Int) = 1 + fun_
# Scalar/Vector interpolations with sdim == rdim (== vdim)
@inline dothelper(A, B) = A ⋅ B
# Vector interpolations with sdim == rdim != vdim
@inline dothelper(A::SMatrix{vdim, dim}, B::Tensor{2, dim}) where {vdim, dim} = A * SMatrix{dim, dim}(B)
@inline dothelper(A::SMatrix{vdim,dim}, B::Tensor{2,dim}) where {vdim,dim} = A * SMatrix{dim,dim}(B)
# Scalar interpolations with sdim > rdim
@inline dothelper(A::SVector{rdim}, B::SMatrix{rdim, sdim}) where {rdim, sdim} = B' * A
@inline dothelper(A::SVector{rdim}, B::SMatrix{rdim,sdim}) where {rdim,sdim} = B' * A
# Vector interpolations with sdim > rdim
@inline dothelper(B::SMatrix{vdim, rdim}, A::SMatrix{rdim, sdim}) where {vdim, rdim, sdim} = B * A
@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}}(
(sum(a -> A[i, a] * B[a, j, k], 1:sdim) for i in 1:vdim, j in 1:rdim, k in 1:rdim)
)
Copy link
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
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}}(
(sum(A[i, j, k] * B[j, k, l, m] for j in 1:rdim, k in 1:rdim) 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)")
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ξ)
H += otimes_helper(x[j], d2Mdξ2)

# J += x[j] ⊗ dMdξ
# H += x[j] ⊗ d2Mdξ2
Copy link
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
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
Loading
Loading