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 e5ea8c3fe4..eea714ccbd 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,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) + ) +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 @@ -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 diff --git a/src/FEValues/GeometryMapping.jl b/src/FEValues/GeometryMapping.jl index 7443821c20..fe47168e36 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) @@ -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 @@ -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 end return MappingValues(J, H) end diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index daffd081b4..ca1e943d1d 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -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] end end