From 0eeecd9b47db5bac310c0224eb803663c9adcb2f Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Tue, 9 Dec 2025 07:55:49 +0000 Subject: [PATCH 1/9] wip first draft --- .vscode/settings.json | 6 +++++ src/FEValues/FunctionValues.jl | 19 +++++++++++++-- src/FEValues/GeometryMapping.jl | 29 ++++++++++++++++++----- test/test_cellvalues.jl | 41 ++++++++++++++++++++++++++++++--- 4 files changed, 84 insertions(+), 11 deletions(-) create mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000000..43356975b0 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,6 @@ +{ + "[julia]": { + "editor.defaultFormatter": "jkillian.custom-local-formatters" + }, + "julia.environmentPath": "/home/USADR/ac136110/.julia/dev/Ferrite", +} \ No newline at end of file diff --git a/src/FEValues/FunctionValues.jl b/src/FEValues/FunctionValues.jl index e5ea8c3fe4..eca7a28b94 100644 --- a/src/FEValues/FunctionValues.jl +++ b/src/FEValues/FunctionValues.jl @@ -174,6 +174,18 @@ 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 +# Result: aᵢBᵢⱼₖ = Cⱼₖ (?) +function dothelper(A::SVector{sdim}, B::SArray{Tuple{sdim, rdim, rdim}}) where {rdim, sdim} + C = zeros(rdim, rdim) + for i in 1:sdim, j in 1:rdim, k in 1:rdim + C[j, k] += A[i] * B[i, j, k] + end + return SMatrix{rdim, rdim}(C) +end + +dothelper(A::SMatrix{sdim, rdim}, B::SMatrix{rdim, rdim}) where {rdim, sdim} = A * B +dothelper(A::SMatrix{sdim, rdim}, B::SMatrix{rdim, sdim}) where {rdim, sdim} = A * B + # ============= # Apply mapping # ============= @@ -199,7 +211,7 @@ end Jinv = calculate_Jinv(getjacobian(mapping_values)) sdim, rdim = size(Jinv) - (rdim != sdim) && error("apply_mapping! for second order gradients and embedded elements not implemented") + # (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 @@ -209,7 +221,10 @@ end if is_vector_valued 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 diff --git a/src/FEValues/GeometryMapping.jl b/src/FEValues/GeometryMapping.jl index 7443821c20..bb28dd96ce 100644 --- a/src/FEValues/GeometryMapping.jl +++ b/src/FEValues/GeometryMapping.jl @@ -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 """ @@ -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 @@ -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) @@ -141,11 +152,14 @@ 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)") + # (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 @@ -170,13 +184,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)") 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 end return MappingValues(J, H) end diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index daffd081b4..5c81ae824f 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -441,9 +441,44 @@ 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)] + # TODO + # @test_throws ErrorException reinit!(cv_vector, coords) # Not implemented for embedded elements + + # 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 ∂A₁∂₁ ≈ H[:, 1, 1] + @test H[:, 2, 2] ≈ ∂A₂∂₂ + @test H[:, 1, 2] ≈ ∂A₁∂₂ + @test H[:, 2, 1] ≈ ∂A₂∂₁ + + # Test Mapping + 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 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] + end end From 06ac744d985c7f703250e5419dbb19dfab41a332 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Tue, 9 Dec 2025 09:52:01 +0000 Subject: [PATCH 2/9] vector-valued + more tests --- .vscode/settings.json | 6 ------ 1 file changed, 6 deletions(-) delete mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json deleted file mode 100644 index 43356975b0..0000000000 --- a/.vscode/settings.json +++ /dev/null @@ -1,6 +0,0 @@ -{ - "[julia]": { - "editor.defaultFormatter": "jkillian.custom-local-formatters" - }, - "julia.environmentPath": "/home/USADR/ac136110/.julia/dev/Ferrite", -} \ No newline at end of file From 53e0939d8f6612a8577b6b60deb8e80f95a0b942 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Tue, 9 Dec 2025 09:52:23 +0000 Subject: [PATCH 3/9] more tests + vector-valued --- src/FEValues/FunctionValues.jl | 36 ++++++++++++++++++++------ test/test_cellvalues.jl | 47 +++++++++++++++++++++++++++++++--- 2 files changed, 71 insertions(+), 12 deletions(-) diff --git a/src/FEValues/FunctionValues.jl b/src/FEValues/FunctionValues.jl index eca7a28b94..04486571ad 100644 --- a/src/FEValues/FunctionValues.jl +++ b/src/FEValues/FunctionValues.jl @@ -176,16 +176,34 @@ required_geo_diff_order(::CovariantPiolaMapping, fun_diff_order::Int) = 1 + fun_ # Result: aᵢBᵢⱼₖ = Cⱼₖ (?) function dothelper(A::SVector{sdim}, B::SArray{Tuple{sdim, rdim, rdim}}) where {rdim, sdim} - C = zeros(rdim, rdim) - for i in 1:sdim, j in 1:rdim, k in 1:rdim - C[j, k] += A[i] * B[i, j, k] - end - return SMatrix{rdim, rdim}(C) + return SMatrix{rdim, rdim}((dot(A, B[:, j, k]) for j in 1:rdim, k in 1:rdim)) +end + +# For Vector Valued +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) + ) end dothelper(A::SMatrix{sdim, rdim}, B::SMatrix{rdim, rdim}) where {rdim, sdim} = A * B dothelper(A::SMatrix{sdim, rdim}, B::SMatrix{rdim, sdim}) where {rdim, sdim} = A * B +@inline dothelper2(A, B) = A ⋅ B + +# dothelper2 helps in disambiguating calls to dothelper +function dothelper2(A::SMatrix{sdim, rdim}, B::SArray{Tuple{vdim, rdim, rdim}}) where {vdim, rdim, sdim} + return SArray{Tuple{vdim, sdim, rdim}}( + ((A * B[i, :, :])[j, k] for i in 1:vdim, j in 1:sdim, k in 1:rdim) + ) +end + +function dothelper2(A::SArray{Tuple{vdim, sdim, rdim}}, B::SMatrix{rdim, sdim}) where {vdim, rdim, sdim} + return SArray{Tuple{vdim, sdim, sdim}}( + ((A[i, :, :] * B)[j, k] for i in 1:vdim, j in 1:sdim, k in 1:sdim) + ) +end + # ============= # Apply mapping # ============= @@ -214,12 +232,14 @@ end # (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 = dothelper2(Jinv', (funvals.d2Ndξ2[j, q_point] - dothelper(dNdx, H))) + d2Ndx2 = dothelper2(t, Jinv) + + # d2Ndx2 = (funvals.d2Ndξ2[j, q_point] - dNdx ⋅ H) ⊡ Jinv_otimesu_Jinv else t = dothelper(Jinv', (funvals.d2Ndξ2[j, q_point] - dothelper(dNdx, H))) d2Ndx2 = dothelper(t, Jinv) diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index 5c81ae824f..ca1e943d1d 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -442,8 +442,6 @@ cv_scalar = CellValues(qr, ip, ip^3; update_hessians = true) coords = [Vec{3}((2x[1], 0.5x[2], rand())) for x in Ferrite.reference_coordinates(ip)] - # TODO - # @test_throws ErrorException reinit!(cv_vector, coords) # Not implemented for embedded elements # Test Scalar H reinit!(cv_scalar, coords) @@ -457,12 +455,30 @@ ∂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 ∂A₁∂₁ ≈ H[:, 1, 1] + @test H[:, 1, 1] ≈ ∂A₁∂₁ @test H[:, 2, 2] ≈ ∂A₂∂₂ @test H[:, 1, 2] ≈ ∂A₁∂₂ @test H[:, 2, 1] ≈ ∂A₂∂₁ - # Test Mapping + + # 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) @@ -473,12 +489,35 @@ 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] end end From 7edb56b5b274a60613ce321b8ecdb565ab20d6ce Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Tue, 9 Dec 2025 13:51:30 +0000 Subject: [PATCH 4/9] fixes for vector valued --- src/FEValues/FunctionValues.jl | 36 +++++++++++++++++---------------- src/FEValues/GeometryMapping.jl | 2 -- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/FEValues/FunctionValues.jl b/src/FEValues/FunctionValues.jl index 04486571ad..62e3772d73 100644 --- a/src/FEValues/FunctionValues.jl +++ b/src/FEValues/FunctionValues.jl @@ -174,33 +174,37 @@ 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 -# Result: aᵢBᵢⱼₖ = Cⱼₖ (?) -function dothelper(A::SVector{sdim}, B::SArray{Tuple{sdim, rdim, rdim}}) where {rdim, sdim} +# 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 # For Vector Valued -function dothelper(A::SMatrix{vdim, sdim}, B::SArray{Tuple{sdim, rdim, rdim}}) where {vdim, rdim, sdim} +@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) ) end -dothelper(A::SMatrix{sdim, rdim}, B::SMatrix{rdim, rdim}) where {rdim, sdim} = A * B -dothelper(A::SMatrix{sdim, rdim}, B::SMatrix{rdim, sdim}) where {rdim, sdim} = A * B +@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 dothelper2(A, B) = A ⋅ B +@inline otimesu_helper(A, B) = otimesu(A, B) -# dothelper2 helps in disambiguating calls to dothelper -function dothelper2(A::SMatrix{sdim, rdim}, B::SArray{Tuple{vdim, rdim, rdim}}) where {vdim, rdim, sdim} - return SArray{Tuple{vdim, sdim, rdim}}( - ((A * B[i, :, :])[j, k] for i in 1:vdim, j in 1:sdim, k in 1:rdim) +# 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 -function dothelper2(A::SArray{Tuple{vdim, sdim, rdim}}, B::SMatrix{rdim, sdim}) where {vdim, rdim, sdim} +@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}}( - ((A[i, :, :] * B)[j, k] for i in 1:vdim, j in 1:sdim, k in 1:sdim) + (dcontract_helper(A[i, :, :], B[:, :, l, m]) for i in 1:vdim, l in 1:sdim, m in 1:sdim) ) end @@ -228,16 +232,14 @@ 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 Union{<:Vec, <:SVector} @inbounds for j in 1:getnbasefunctions(funvals) dNdx = dothelper(funvals.dNdξ[j, q_point], Jinv) if is_vector_valued - t = dothelper2(Jinv', (funvals.d2Ndξ2[j, q_point] - dothelper(dNdx, H))) - d2Ndx2 = dothelper2(t, 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 diff --git a/src/FEValues/GeometryMapping.jl b/src/FEValues/GeometryMapping.jl index bb28dd96ce..0707530f86 100644 --- a/src/FEValues/GeometryMapping.jl +++ b/src/FEValues/GeometryMapping.jl @@ -151,8 +151,6 @@ 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 += otimes_helper(x[j], geo_mapping.dMdξ[j, q_point]) From 8737dd3f72e3160e3e5ac906ed5ae12cc63e80d6 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Wed, 10 Dec 2025 13:35:54 +0000 Subject: [PATCH 5/9] update --- .gitignore | 1 + src/FEValues/FunctionValues.jl | 1 - 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 07fd400412..5773c480ef 100644 --- a/.gitignore +++ b/.gitignore @@ -24,3 +24,4 @@ docs/src/changelog.md LocalPreferences.toml docs/LocalPreferences.toml benchmark/tune.json +.vscode/ diff --git a/src/FEValues/FunctionValues.jl b/src/FEValues/FunctionValues.jl index 62e3772d73..44c8f0258a 100644 --- a/src/FEValues/FunctionValues.jl +++ b/src/FEValues/FunctionValues.jl @@ -179,7 +179,6 @@ required_geo_diff_order(::CovariantPiolaMapping, fun_diff_order::Int) = 1 + fun_ return SMatrix{rdim, rdim}((dot(A, B[:, j, k]) for j in 1:rdim, k in 1:rdim)) end -# For Vector Valued @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) From 158cc2bf2dbe28854d5316d9670e098bb04f142e Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Mon, 15 Dec 2025 08:19:01 +0000 Subject: [PATCH 6/9] fix computation of H in calculate_mapping(gip::ScalarInterpolation) --- src/FEValues/GeometryMapping.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/FEValues/GeometryMapping.jl b/src/FEValues/GeometryMapping.jl index 0707530f86..fe47168e36 100644 --- a/src/FEValues/GeometryMapping.jl +++ b/src/FEValues/GeometryMapping.jl @@ -182,13 +182,13 @@ 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 += otimes_helper(x[j], dMdξ) + H += otimes_helper(x[j], d2Mdξ2) # J += x[j] ⊗ dMdξ # H += x[j] ⊗ d2Mdξ2 From 8011c2f5c0dbf3f62556dbbf679b065a937741c7 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Mon, 15 Dec 2025 09:36:46 +0000 Subject: [PATCH 7/9] more performance --- src/FEValues/FunctionValues.jl | 96 +++++++++++++++++----------------- 1 file changed, 48 insertions(+), 48 deletions(-) diff --git a/src/FEValues/FunctionValues.jl b/src/FEValues/FunctionValues.jl index 44c8f0258a..83a5e5a08a 100644 --- a/src/FEValues/FunctionValues.jl +++ b/src/FEValues/FunctionValues.jl @@ -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} """ @@ -43,7 +43,7 @@ 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}} @@ -51,17 +51,17 @@ struct FunctionValues{DiffOrder, IP, N_t, dNdx_t, dNdξ_t, d2Ndx2_t, d2Ndξ2_t} 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) @@ -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) @@ -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 @@ -168,42 +168,42 @@ 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)) +@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) +@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) ) 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 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}}( +@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) +@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) +@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 a in j:rdim, b in k:rdim) for i in 1:vdim, l in 1:sdim, m in 1:sdim) ) end @@ -232,7 +232,7 @@ end Jinv = calculate_Jinv(getjacobian(mapping_values)) H = gethessian(mapping_values) - is_vector_valued = first(funvals.Nx) isa Union{<:Vec, <:SVector} + 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 From 6f7fbbdfa1d8c53e874580b8f183d071889100d7 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Mon, 15 Dec 2025 09:43:21 +0000 Subject: [PATCH 8/9] fix --- src/FEValues/FunctionValues.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/FEValues/FunctionValues.jl b/src/FEValues/FunctionValues.jl index 83a5e5a08a..83f91f2657 100644 --- a/src/FEValues/FunctionValues.jl +++ b/src/FEValues/FunctionValues.jl @@ -203,7 +203,7 @@ end # 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 a in j:rdim, b in k:rdim) for i in 1:vdim, l in 1:sdim, m in 1:sdim) + (sum(A[i, j, k] * B[j, k, l, m] for a in j:rdim, k in k:rdim) for i in 1:vdim, l in 1:sdim, m in 1:sdim) ) end From 7260a95f84158b430537947225ec274bc94457b4 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Mon, 15 Dec 2025 09:48:21 +0000 Subject: [PATCH 9/9] fix2 --- src/FEValues/FunctionValues.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/FEValues/FunctionValues.jl b/src/FEValues/FunctionValues.jl index 83f91f2657..eea714ccbd 100644 --- a/src/FEValues/FunctionValues.jl +++ b/src/FEValues/FunctionValues.jl @@ -203,7 +203,7 @@ end # 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 a in j:rdim, k in k:rdim) for i in 1:vdim, l in 1:sdim, m in 1: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