diff --git a/docs/src/topics/FEValues.md b/docs/src/topics/FEValues.md index 1a4f836be8..3d0ebc662a 100644 --- a/docs/src/topics/FEValues.md +++ b/docs/src/topics/FEValues.md @@ -42,6 +42,49 @@ For scalar fields, we always use scalar base functions. For tensorial fields (no \end{align*} ``` +Second order gradients of the shape functions are computed as + +```math +\begin{align*} + \mathrm{grad}(\mathrm{grad}(N(\boldsymbol{x}))) = \frac{\mathrm{d}^2 N}{\mathrm{d}\boldsymbol{x}^2} = \boldsymbol{J}^{-T} \cdot \frac{\mathrm{d}^2\hat{N}}{\mathrm{d}\boldsymbol{\xi}^2} \cdot \boldsymbol{J}^{-1} - \boldsymbol{J}^{-T} \cdot\mathrm{grad}(N) \cdot \boldsymbol{\mathcal{H}} \cdot \boldsymbol{J}^{-1} +\end{align*} +``` +!!! details "Derivation" + The gradient of the shape functions is obtained using the chain rule: + + \begin{align*} + \frac{\mathrm{d} N}{\mathrm{d}x_i} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\frac{\mathrm{d} \xi_r}{\mathrm{d} x_i} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r} J^{-1}_{ri} + \end{align*} + + For the second order gradients, we first use the product rule on the equation above: + + \begin{align*} + \frac{\mathrm{d}^2 N}{\mathrm{d}x_i \mathrm{d}x_j} = \frac{\mathrm{d}}{\mathrm{d}x_j}(\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}) J^{-1}_{ri} + \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r} \frac{\mathrm{d}}{\mathrm{d}x_j}(J^{-1}_{ri}) + \end{align*} + + The first term can be computed as: + + \begin{align*} + \frac{\mathrm{d}}{\mathrm{d}x_j}(\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}) J^{-1}_{ri} = J^{-1}_{si}\frac{\mathrm{d}}{\mathrm{d}\xi_s}(\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}) J^{-1}_{ri} = J^{-1}_{si}\frac{\mathrm{d}^2 \hat N}{\mathrm{d} \xi_s\mathrm{d} \xi_r} J^{-1}_{ri} + \end{align*} + + The second term can be written as: + + \begin{align*} + \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\frac{\mathrm{d}}{\mathrm{d}x_j}(J^{-1}_{ri}) = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}[\frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}\xi_s}]\frac{\mathrm{d} \xi_s}{\mathrm{d}x_j} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}[ - J^{-1}_{ri}\mathcal{H}_{ris} J^{-1}_{li}]J^{-1}_{sj} = - \frac{\mathrm{d} \hat N}{\mathrm{d} x_i}\mathcal{H}_{ris} J^{-1}_{li}J^{-1}_{sj} + \end{align*} + + where we have used that the inverse of the jacobian can be computed as: + + \begin{align*} + 0 = \frac{\mathrm{d}}{\mathrm{d}\xi_s} (J_{rl} J^{-1}_{li} ) = \frac{\mathrm{d}J_{rl}}{\mathrm{d}\xi_s} J^{-1}_{li} + J_{ri} \frac{\mathrm{d}J^{-1}_{li}}{\mathrm{d}\xi_s} = 0 \quad \Rightarrow \\ + \end{align*} + + \begin{align*} + \frac{\mathrm{d}J^{-1}_{li}}{\mathrm{d}\xi_s} = - J^{-1}_{ri}\frac{\mathrm{d}J_{rl}}{\mathrm{d}\xi_s} J^{-1}_{li} = - J^{-1}_{ri}\mathcal{H}_{ris} J^{-1}_{li}\\ + \end{align*} + + ### Covariant Piola mapping, H(curl) `Ferrite.CovariantPiolaMapping` diff --git a/src/FEValues/CellValues.jl b/src/FEValues/CellValues.jl index 5383381ed7..0b907bcc2b 100644 --- a/src/FEValues/CellValues.jl +++ b/src/FEValues/CellValues.jl @@ -42,10 +42,17 @@ struct CellValues{FV, GM, QR, detT} <: AbstractCellValues detJdV::detT # AbstractVector{<:Number} or Nothing end function CellValues(::Type{T}, qr::QuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation; - update_gradients::Union{Bool,Nothing} = nothing, update_detJdV::Union{Bool,Nothing} = nothing) where T + update_gradients::Union{Bool,Nothing} = nothing, + update_hessians::Union{Bool,Nothing} = nothing, + update_detJdV::Union{Bool,Nothing} = nothing) where T - _update_detJdV = update_detJdV === nothing ? true : update_detJdV - FunDiffOrder = update_gradients === nothing ? 1 : convert(Int, update_gradients) # Logic must change when supporting update_hessian kwargs + _update_gradients = update_gradients === nothing ? true : update_gradients + _update_hessians = update_hessians === nothing ? false : update_hessians + _update_detJdV = update_detJdV === nothing ? true : update_detJdV + _update_hessians && @assert _update_gradients + + FunDiffOrder = convert(Int, _update_gradients) + FunDiffOrder += convert(Int, _update_hessians) GeoDiffOrder = max(required_geo_diff_order(mapping_type(ip_fun), FunDiffOrder), _update_detJdV) geo_mapping = GeometryMapping{GeoDiffOrder}(T, ip_geo.ip, qr) fun_values = FunctionValues{FunDiffOrder}(T, ip_fun, qr, ip_geo) @@ -79,6 +86,7 @@ shape_gradient_type(cv::CellValues) = shape_gradient_type(cv.fun_values) @propagate_inbounds shape_value(cv::CellValues, q_point::Int, i::Int) = shape_value(cv.fun_values, q_point, i) @propagate_inbounds shape_gradient(cv::CellValues, q_point::Int, i::Int) = shape_gradient(cv.fun_values, q_point, i) +@propagate_inbounds shape_hessian(cv::CellValues, q_point::Int, i::Int) = shape_hessian(cv.fun_values, q_point, i) @propagate_inbounds shape_symmetric_gradient(cv::CellValues, q_point::Int, i::Int) = shape_symmetric_gradient(cv.fun_values, q_point, i) # Access quadrature rule values diff --git a/src/FEValues/FaceValues.jl b/src/FEValues/FaceValues.jl index d9353fc7aa..d38f155196 100644 --- a/src/FEValues/FaceValues.jl +++ b/src/FEValues/FaceValues.jl @@ -41,9 +41,15 @@ struct FaceValues{FV, GM, FQR, detT, nT, V_FV<:AbstractVector{FV}, V_GM<:Abstrac end function FaceValues(::Type{T}, fqr::FaceQuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation{sdim} = default_geometric_interpolation(ip_fun); - update_gradients::Union{Bool,Nothing} = nothing) where {T,sdim} - - FunDiffOrder = update_gradients === nothing ? 1 : convert(Int, update_gradients) # Logic must change when supporting update_hessian kwargs + update_gradients::Union{Bool,Nothing} = nothing, + update_hessians::Union{Bool,Nothing} = nothing) where {T,sdim} + + _update_gradients = update_gradients === nothing ? true : update_gradients + _update_hessians = update_hessians === nothing ? false : update_hessians + _update_hessians && @assert _update_gradients + + FunDiffOrder = convert(Int, _update_gradients) + FunDiffOrder += convert(Int, _update_hessians) GeoDiffOrder = max(required_geo_diff_order(mapping_type(ip_fun), FunDiffOrder), 1) # Not sure why the type-asserts are required here but not for CellValues, # but they solve the type-instability for FaceValues @@ -84,6 +90,7 @@ get_fun_values(fv::FaceValues) = @inbounds fv.fun_values[getcurrentface(fv)] @propagate_inbounds shape_value(fv::FaceValues, q_point::Int, i::Int) = shape_value(get_fun_values(fv), q_point, i) @propagate_inbounds shape_gradient(fv::FaceValues, q_point::Int, i::Int) = shape_gradient(get_fun_values(fv), q_point, i) +@propagate_inbounds shape_hessian(fv::FaceValues, q_point::Int, i::Int) = shape_hessian(get_fun_values(fv), q_point, i) @propagate_inbounds shape_symmetric_gradient(fv::FaceValues, q_point::Int, i::Int) = shape_symmetric_gradient(get_fun_values(fv), q_point, i) """ diff --git a/src/FEValues/FunctionValues.jl b/src/FEValues/FunctionValues.jl index 361f9a6619..974a331408 100644 --- a/src/FEValues/FunctionValues.jl +++ b/src/FEValues/FunctionValues.jl @@ -6,24 +6,34 @@ ################################################################# # 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_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_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} + # 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} -typeof_dNdξ(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, vdim, sdim, rdim} = SMatrix{vdim, rdim, 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} + """ FunctionValues{DiffOrder}(::Type{T}, ip_fun, qr::QuadratureRule, ip_geo::VectorizedInterpolation) @@ -33,36 +43,51 @@ for both the reference cell (precalculated) and the real cell (updated in `reini """ FunctionValues -struct FunctionValues{DiffOrder, IP, N_t, dNdx_t, dNdξ_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 - function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, ::Nothing, ::Nothing) where {N_t<:AbstractMatrix} - return new{0, typeof(ip), N_t, Nothing, Nothing}(ip, Nx, Nξ, nothing, 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) + 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) end - function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, dNdx::AbstractMatrix, dNdξ::AbstractMatrix) where {N_t<:AbstractMatrix} - return new{1, typeof(ip), N_t, typeof(dNdx), typeof(dNdξ)}(ip, Nx, Nξ, dNdx, dNdξ) + 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} + sdim = n_components(ip_geo) + rdim = getdim(ip) + ((rdim != sdim) && DiffOrder>1) && error("CellValues/FunctionValues for embedded elements and higher order gradient not implemented") + n_shape = getnbasefunctions(ip) n_qpoints = getnquadpoints(qr) Nξ = zeros(typeof_N(T, ip, ip_geo), n_shape, n_qpoints) Nx = isa(mapping_type(ip), IdentityMapping) ? Nξ : similar(Nξ) - - if DiffOrder == 0 - dNdξ = dNdx = nothing - elseif DiffOrder == 1 + dNdξ = dNdx = d2Ndξ2 = d2Ndx2 = nothing + + if DiffOrder >= 1 dNdξ = zeros(typeof_dNdξ(T, ip, ip_geo), n_shape, n_qpoints) dNdx = fill(zero(typeof_dNdx(T, ip, ip_geo)) * T(NaN), n_shape, n_qpoints) - else + end + + if DiffOrder >= 2 + d2Ndξ2 = zeros(typeof_d2Ndξ2(T, ip, ip_geo), n_shape, n_qpoints) + d2Ndx2 = fill(zero(typeof_d2Ndx2(T, ip, ip_geo)) * T(NaN), n_shape, n_qpoints) + end + + if DiffOrder > 2 throw(ArgumentError("Currently only values and gradients can be updated in FunctionValues")) end - fv = FunctionValues(ip, Nx, Nξ, dNdx, dNdξ) + fv = FunctionValues(ip, Nx, Nξ, dNdx, dNdξ, d2Ndx2, d2Ndξ2) precompute_values!(fv, getpoints(qr)) # Separate function for qr point update in PointValues return fv end @@ -73,18 +98,24 @@ end function precompute_values!(fv::FunctionValues{1}, qr_points::Vector{<:Vec}) shape_gradients_and_values!(fv.dNdξ, fv.Nξ, fv.ip, qr_points) end +function precompute_values!(fv::FunctionValues{2}, qr_points::Vector{<:Vec}) + shape_hessians_gradients_and_values!(fv.d2Ndξ2, fv.dNdξ, fv.Nξ, fv.ip, qr_points) +end function Base.copy(v::FunctionValues) Nξ_copy = copy(v.Nξ) Nx_copy = v.Nξ === v.Nx ? Nξ_copy : copy(v.Nx) # Preserve aliasing dNdx_copy = _copy_or_nothing(v.dNdx) dNdξ_copy = _copy_or_nothing(v.dNdξ) - return FunctionValues(copy(v.ip), Nx_copy, Nξ_copy, dNdx_copy, dNdξ_copy) + d2Ndx2_copy = _copy_or_nothing(v.d2Ndx2) + d2Ndξ2_copy = _copy_or_nothing(v.d2Ndξ2) + return FunctionValues(copy(v.ip), Nx_copy, Nξ_copy, dNdx_copy, dNdξ_copy, d2Ndx2_copy, d2Ndξ2_copy) end getnbasefunctions(funvals::FunctionValues) = size(funvals.Nx, 1) @propagate_inbounds shape_value(funvals::FunctionValues, q_point::Int, base_func::Int) = funvals.Nx[base_func, q_point] @propagate_inbounds shape_gradient(funvals::FunctionValues, q_point::Int, base_func::Int) = funvals.dNdx[base_func, q_point] +@propagate_inbounds shape_hessian(funvals::FunctionValues{2}, q_point::Int, base_func::Int) = funvals.d2Ndx2[base_func, q_point] @propagate_inbounds shape_symmetric_gradient(funvals::FunctionValues, q_point::Int, base_func::Int) = symmetric(shape_gradient(funvals, q_point, base_func)) function_interpolation(funvals::FunctionValues) = funvals.ip @@ -92,6 +123,9 @@ function_difforder(::FunctionValues{DiffOrder}) where DiffOrder = DiffOrder shape_value_type(funvals::FunctionValues) = eltype(funvals.Nx) shape_gradient_type(funvals::FunctionValues) = eltype(funvals.dNdx) shape_gradient_type(::FunctionValues{0}) = nothing +shape_hessian_type(funvals::FunctionValues) = eltype(funvals.d2Ndx2) +shape_hessian_type(::FunctionValues{0}) = nothing +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 @@ -167,6 +201,25 @@ end return nothing end +@inline function apply_mapping!(funvals::FunctionValues{2}, ::IdentityMapping, q_point::Int, mapping_values, args...) + Jinv = calculate_Jinv(getjacobian(mapping_values)) + H = gethessian(mapping_values) + is_vector_valued = first(funvals.Nx) isa Vec + Jinv_otimesu_Jinv = is_vector_valued ? otimesu(Jinv, Jinv) : nothing + @inbounds for j in 1:getnbasefunctions(funvals) + dNdx = dothelper(funvals.dNdξ[j, q_point], Jinv) + if is_vector_valued #TODO - combine with helper function ? + d2Ndx2 = (funvals.d2Ndξ2[j, q_point] - dNdx⋅H) ⊡ Jinv_otimesu_Jinv + else + d2Ndx2 = Jinv'⋅(funvals.d2Ndξ2[j, q_point] - dNdx⋅H)⋅Jinv + end + + funvals.dNdx[j, q_point] = dNdx + funvals.d2Ndx2[j, q_point] = d2Ndx2 + end + return nothing +end + # TODO in PR798, apply_mapping! for # * CovariantPiolaMapping # * ContravariantPiolaMapping diff --git a/src/FEValues/GeometryMapping.jl b/src/FEValues/GeometryMapping.jl index b6346b7550..fdf5798ef3 100644 --- a/src/FEValues/GeometryMapping.jl +++ b/src/FEValues/GeometryMapping.jl @@ -14,7 +14,7 @@ struct MappingValues{JT, HT} end @inline getjacobian(mv::MappingValues{<:Union{AbstractTensor, SMatrix}}) = mv.J -# @inline gethessian(mv::MappingValues{<:Any,<:AbstractTensor}) = mv.H # PR798 +@inline gethessian(mv::MappingValues{<:Any,<:AbstractTensor}) = mv.H """ @@ -46,12 +46,12 @@ struct GeometryMapping{DiffOrder, IP, M_t, dMdξ_t, d2Mdξ2_t} ) where {IP <: ScalarInterpolation, M_t<:AbstractMatrix{<:Number}, dMdξ_t <: AbstractMatrix{<:Vec}} return new{1, IP, M_t, dMdξ_t, Nothing}(ip, M, dMdξ, nothing) end -#= function GeometryMapping( + function GeometryMapping( ip::IP, M::M_t, dMdξ::dMdξ_t, d2Mdξ2::d2Mdξ2_t) where {IP <: ScalarInterpolation, M_t<:AbstractMatrix{<:Number}, dMdξ_t <: AbstractMatrix{<:Vec}, d2Mdξ2_t <: AbstractMatrix{<:Tensor{2}}} return new{2, IP, M_t, dMdξ_t, d2Mdξ2_t}(ip, M, dMdξ, d2Mdξ2) - end =# # PR798 + end end function GeometryMapping{0}(::Type{T}, ip::ScalarInterpolation, qr::QuadratureRule) where T n_shape = getnbasefunctions(ip) @@ -71,7 +71,7 @@ function GeometryMapping{1}(::Type{T}, ip::ScalarInterpolation, qr::QuadratureRu precompute_values!(gm, getpoints(qr)) return gm end -#= function GeometryMapping{2}(::Type{T}, ip::ScalarInterpolation, qr::QuadratureRule) where T +function GeometryMapping{2}(::Type{T}, ip::ScalarInterpolation, qr::QuadratureRule) where T n_shape = getnbasefunctions(ip) n_qpoints = getnquadpoints(qr) @@ -82,7 +82,7 @@ end gm = GeometryMapping(ip, M, dMdξ, d2Mdξ2) precompute_values!(gm, getpoints(qr)) return gm -end =# # PR798 +end function precompute_values!(gm::GeometryMapping{0}, qr_points::Vector{<:Vec}) shape_values!(gm.M, gm.ip, qr_points) @@ -90,9 +90,9 @@ end function precompute_values!(gm::GeometryMapping{1}, qr_points::Vector{<:Vec}) shape_gradients_and_values!(gm.dMdξ, gm.M, gm.ip, qr_points) end -#= function precompute_values!(gm::GeometryMapping{2}, qr_points::Vector{<:Vec}) +function precompute_values!(gm::GeometryMapping{2}, qr_points::Vector{<:Vec}) shape_hessians_gradients_and_values!(gm.d2Mdξ2, gm.dMdξ, gm.M, gm.ip, qr_points) -end =# # PR798 +end function Base.copy(v::GeometryMapping) return GeometryMapping(copy(v.ip), copy(v.M), _copy_or_nothing(v.dMdξ), _copy_or_nothing(v.d2Mdξ2)) @@ -117,9 +117,9 @@ end function otimes_returntype(#=typeof(x)=#::Type{<:Vec{dim,Tx}}, #=typeof(dMdξ)=#::Type{<:Vec{dim,TM}}) where {dim, Tx, TM} return Tensor{2,dim,promote_type(Tx,TM)} end -#= function otimes_returntype(#=typeof(x)=#::Type{<:Vec{dim,Tx}}, #=typeof(d2Mdξ2)=#::Type{<:Tensor{2,dim,TM}}) where {dim, Tx, TM} +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 =# # PR798 +end @inline function calculate_mapping(::GeometryMapping{0}, q_point, x) return MappingValues(nothing, nothing) @@ -134,7 +134,7 @@ end return MappingValues(fecv_J, nothing) end -#= @inline function calculate_mapping(geo_mapping::GeometryMapping{2}, q_point, x) +@inline function calculate_mapping(geo_mapping::GeometryMapping{2}, q_point, x) J = zero(otimes_returntype(eltype(x), eltype(geo_mapping.dMdξ))) H = zero(otimes_returntype(eltype(x), eltype(geo_mapping.d2Mdξ2))) @inbounds for j in 1:getngeobasefunctions(geo_mapping) @@ -142,7 +142,7 @@ end H += x[j] ⊗ geo_mapping.d2Mdξ2[j, q_point] end return MappingValues(J, H) -end =# # PR798 +end calculate_detJ(J::Tensor{2}) = det(J) calculate_detJ(J::SMatrix) = embedding_det(J) diff --git a/src/FEValues/common_values.jl b/src/FEValues/common_values.jl index 181e54b03c..088bbae7ff 100644 --- a/src/FEValues/common_values.jl +++ b/src/FEValues/common_values.jl @@ -205,6 +205,41 @@ function function_gradient_init(cv::AbstractValues, ::AbstractVector{T}) where { return zero(T) ⊗ zero(shape_gradient_type(cv)) end +""" + function_hessian(fe_v::AbstractValues{dim}, q_point::Int, u::AbstractVector{<:AbstractFloat}, [dof_range]) + + Compute the hessian of the function in a quadrature point. `u` is a vector with values + for the degrees of freedom. +""" +function function_hessian(fe_v::AbstractValues, q_point::Int, u::AbstractVector, dof_range = eachindex(u)) + n_base_funcs = getnbasefunctions(fe_v) + length(dof_range) == n_base_funcs || throw_incompatible_dof_length(length(dof_range), n_base_funcs) + @boundscheck checkbounds(u, dof_range) + @boundscheck checkquadpoint(fe_v, q_point) + hess = function_hessian_init(fe_v, u) + @inbounds for (i, j) in pairs(dof_range) + hess += shape_hessian(fe_v, q_point, i) * u[j] + end + return hess +end + +""" + shape_hessian_type(fe_v::AbstractValues) + +Return the type of `shape_hessian(fe_v, q_point, base_function)` +""" +function shape_hessian_type(fe_v::AbstractValues) + # Default fallback + return typeof(shape_hessian(fe_v, 1, 1)) +end + +function function_hessian_init(cv::AbstractValues, ::AbstractVector{T}) where {T} + return zero(shape_hessian_type(cv)) * zero(T) +end +function function_hessian_init(cv::AbstractValues, ::AbstractVector{T}) where {T <: AbstractVector} + return zero(T) ⊗ zero(shape_hessian_type(cv)) +end + """ function_symmetric_gradient(fe_v::AbstractValues, q_point::Int, u::AbstractVector, [dof_range]) @@ -302,11 +337,9 @@ function shape_gradients_and_values!(gradients::AbstractMatrix, values::Abstract end end -#= PR798 function shape_hessians_gradients_and_values!(hessians::AbstractMatrix, gradients::AbstractMatrix, values::AbstractMatrix, ip, qr_points::Vector{<:Vec}) for (qp, ξ) in pairs(qr_points) shape_hessians_gradients_and_values!(@view(hessians[:, qp]), @view(gradients[:, qp]), @view(values[:, qp]), ip, ξ) end end -=# diff --git a/src/exports.jl b/src/exports.jl index dab7e368fb..9ee5f580f7 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -33,6 +33,7 @@ export reinit!, shape_value, shape_gradient, + shape_hessian, shape_symmetric_gradient, shape_divergence, shape_curl, diff --git a/src/interpolations.jl b/src/interpolations.jl index 537759b15c..17e1194c6d 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -229,7 +229,6 @@ function shape_gradients_and_values!(gradients::GAT, values::SAT, ip::IP, ξ::Ve end end -#= PR798 """ shape_hessians_gradients_and_values!(hessians::AbstractVector, gradients::AbstractVector, values::AbstractVector, ip::Interpolation, ξ::Vec) @@ -244,7 +243,7 @@ and store them in `hessians`, `gradients`, and `values`. hessians[i], gradients[i], values[i] = shape_hessian_gradient_and_value(ip, ξ, i) end end -=# + """ shape_value(ip::Interpolation, ξ::Vec, i::Int) @@ -279,7 +278,6 @@ function shape_gradient_and_value(ip::Interpolation, ξ::Vec, i::Int) return gradient(x -> shape_value(ip, x, i), ξ, :all) end -#= PR798 """ shape_hessian_gradient_and_value(ip::Interpolation, ξ::Vec, i::Int) @@ -289,7 +287,7 @@ Optimized version combining the evaluation [`Ferrite.shape_value(::Interpolation function shape_hessian_gradient_and_value(ip::Interpolation, ξ::Vec, i::Int) return hessian(x -> shape_value(ip, x, i), ξ, :all) end -=# + """ reference_coordinates(ip::Interpolation) @@ -1560,6 +1558,42 @@ function shape_gradient_and_value(ipv::VectorizedInterpolation{vdim, shape}, ξ: return SMatrix(grad), val end +# vdim == refdim +function shape_hessian_gradient_and_value(ipv::VectorizedInterpolation{dim, shape}, ξ::Vec{dim}, I::Int) where {dim, shape <: AbstractRefShape{dim}} + return invoke(shape_hessian_gradient_and_value, Tuple{Interpolation, Vec, Int}, ipv, ξ, I) +end +# vdim != refdim +function shape_hessian_gradient_and_value(ipv::VectorizedInterpolation{vdim, shape}, ξ::V, I::Int) where {vdim, refdim, shape <: AbstractRefShape{refdim}, T, V <: Vec{refdim, T}} + _shape_hessian_gradient_and_value_static_array(ipv, ξ, I) +end +function _shape_hessian_gradient_and_value_static_array(ipv::VectorizedInterpolation{vdim, shape}, ξ::V, I::Int) where {vdim, refdim, shape <: AbstractRefShape{refdim}, T, V <: Vec{refdim, T}} + # Load with dual numbers and compute the value + f = x -> shape_value(ipv, x, I) + ξd = Tensors._load(Tensors._load(ξ, Tensors.Tag(f, V)), Tensors.Tag(f, V)) + value_hess = f(ξd) + # Extract the value and gradient + val = Vec{vdim, T}(i -> Tensors.value(Tensors.value(value_hess[i]))) + grad = zero(MMatrix{vdim, refdim, T}) + hess = zero(MArray{Tuple{vdim, refdim, refdim}, T}) + for (i, vi) in pairs(value_hess) + hess_values = Tensors.value(vi) + + hess_values_partials = Tensors.partials(hess_values) + for (k, pk) in pairs(hess_values_partials) + grad[i, k] = pk + end + + hess_partials = Tensors.partials(vi) + for (j, partial_j) in pairs(hess_partials) + hess_partials_partials = Tensors.partials(partial_j) + for (k, pk) in pairs(hess_partials_partials) + hess[i, j, k] = pk + end + end + end + return SArray(hess), SMatrix(grad), val +end + reference_coordinates(ip::VectorizedInterpolation) = reference_coordinates(ip.ip) is_discontinuous(::Type{<:VectorizedInterpolation{<:Any, <:Any, <:Any, ip}}) where {ip} = is_discontinuous(ip) diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index 731060f682..f544b54fcc 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -1,5 +1,5 @@ @testset "CellValues" begin -@testset "ip=$scalar_interpol quad_rule=$(typeof(quad_rule))" for (scalar_interpol, quad_rule) in ( +@testset "ip=$scalar_interpol" for (scalar_interpol, quad_rule) in ( (Lagrange{RefLine, 1}(), QuadratureRule{RefLine}(2)), (Lagrange{RefLine, 2}(), QuadratureRule{RefLine}(2)), (Lagrange{RefQuadrilateral, 1}(), QuadratureRule{RefQuadrilateral}(2)), @@ -16,56 +16,102 @@ (Lagrange{RefPrism, 2}(), QuadratureRule{RefPrism}(2)), (Lagrange{RefPyramid, 2}(), QuadratureRule{RefPyramid}(2)), ) - - for func_interpol in (scalar_interpol, VectorizedInterpolation(scalar_interpol)) + for func_interpol in (scalar_interpol, VectorizedInterpolation(scalar_interpol)), DiffOrder in 1:2 + (DiffOrder==2 && Ferrite.getorder(func_interpol)==1) && continue #No need to test linear interpolations again geom_interpol = scalar_interpol # Tests below assume this n_basefunc_base = getnbasefunctions(scalar_interpol) - cv = @inferred CellValues(quad_rule, func_interpol, geom_interpol) + update_gradients = true + update_hessians = (DiffOrder==2 && Ferrite.getorder(func_interpol) > 1) + cv = CellValues(quad_rule, func_interpol, geom_interpol; update_gradients, update_hessians) ndim = Ferrite.getdim(func_interpol) n_basefuncs = getnbasefunctions(func_interpol) @test getnbasefunctions(cv) == n_basefuncs - x, n = valid_coordinates_and_normals(func_interpol) - reinit!(cv, x) - @test_call reinit!(cv, x) + coords, n = valid_coordinates_and_normals(func_interpol) + reinit!(cv, coords) + @test_call reinit!(cv, coords) # We test this by applying a given deformation gradient on all the nodes. # Since this is a linear deformation we should get back the exact values # from the interpolation. - u = Vec{ndim, Float64}[zero(Tensor{1,ndim}) for i in 1:n_basefunc_base] - u_scal = zeros(n_basefunc_base) - H = rand(Tensor{2, ndim}) - V = rand(Tensor{1, ndim}) + V, G, H = if func_interpol isa Ferrite.ScalarInterpolation + (rand(), rand(Tensor{1, ndim}), Tensor{2, ndim}((i,j)-> i==j ? rand() : 0.0)) + else + (rand(Tensor{1, ndim}), rand(Tensor{2, ndim}), Tensor{3, ndim}((i,j,k)-> i==j==k ? rand() : 0.0)) + end + + u_funk(x,V,G,H) = begin + if update_hessians + 0.5*x⋅H⋅x + G⋅x + V + else + G⋅x + V + end + end + + _ue = [u_funk(coords[i],V,G,H) for i in 1:n_basefunc_base] + ue = reinterpret(Float64, _ue) + + for i in 1:getnquadpoints(cv) + xqp = spatial_coordinate(cv, i, coords) + Hqp, Gqp, Vqp = Tensors.hessian(x -> u_funk(x,V,G,H), xqp, :all) + + @test function_value(cv, i, ue) ≈ Vqp + @test function_gradient(cv, i, ue) ≈ Gqp + if update_hessians + #Note, the jacobian of the element is constant, which makes the hessian (of the mapping) + #zero. So this is not the optimal test + @test Ferrite.function_hessian(cv, i, ue) ≈ Hqp + end + if func_interpol isa Ferrite.VectorInterpolation + @test function_symmetric_gradient(cv, i, ue) ≈ 0.5(Gqp + Gqp') + @test function_divergence(cv, i, ue) ≈ tr(Gqp) + ndim == 3 && @test function_curl(cv, i, ue) ≈ Ferrite.curl_from_gradient(Gqp) + else + @test function_divergence(cv, i, ue) ≈ sum(Gqp) + end + end + + #Test CellValues when input is a ::Vector{<:Vec} (most of which is deprecated) + ue_vec = [zero(Vec{ndim,Float64}) for i in 1:n_basefunc_base] + G_vector = rand(Tensor{2, ndim}) for i in 1:n_basefunc_base - u[i] = H ⋅ x[i] - u_scal[i] = V ⋅ x[i] + ue_vec[i] = G_vector ⋅ coords[i] end - u_vector = reinterpret(Float64, u) for i in 1:getnquadpoints(cv) if func_interpol isa Ferrite.ScalarInterpolation - @test function_gradient(cv, i, u) ≈ H - @test function_symmetric_gradient(cv, i, u) ≈ 0.5(H + H') - @test function_divergence(cv, i, u_scal) ≈ sum(V) - @test function_divergence(cv, i, u) ≈ tr(H) - @test function_gradient(cv, i, u_scal) ≈ V - ndim == 3 && @test function_curl(cv, i, u) ≈ Ferrite.curl_from_gradient(H) - function_value(cv, i, u) - function_value(cv, i, u_scal) + @test function_gradient(cv, i, ue_vec) ≈ G_vector else# func_interpol isa Ferrite.VectorInterpolation - @test function_gradient(cv, i, u_vector) ≈ H - @test (@test_deprecated function_gradient(cv, i, u)) ≈ H - @test function_symmetric_gradient(cv, i, u_vector) ≈ 0.5(H + H') - @test (@test_deprecated function_symmetric_gradient(cv, i, u)) ≈ 0.5(H + H') - @test function_divergence(cv, i, u_vector) ≈ tr(H) - @test (@test_deprecated function_divergence(cv, i, u)) ≈ tr(H) + @test (@test_deprecated function_gradient(cv, i, ue_vec)) ≈ G_vector + @test (@test_deprecated function_symmetric_gradient(cv, i, ue_vec)) ≈ 0.5(G_vector + G_vector') + @test (@test_deprecated function_divergence(cv, i, ue_vec)) ≈ tr(G_vector) if ndim == 3 - @test function_curl(cv, i, u_vector) ≈ Ferrite.curl_from_gradient(H) - @test (@test_deprecated function_curl(cv, i, u)) ≈ Ferrite.curl_from_gradient(H) + @test (@test_deprecated function_curl(cv, i, ue_vec)) ≈ Ferrite.curl_from_gradient(G_vector) + end + @test_deprecated function_value(cv, i, ue_vec) + end + end + + #Check if the non-linear mapping is correct + #Only do this for one interpolation becuase it relise on AD on "iterative function" + if scalar_interpol === Lagrange{RefQuadrilateral, 2}() + 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 - @test function_value(cv, i, u_vector) ≈ (@test_deprecated function_value(cv, i, u)) end + reinit!(cv, coords) # reinit back to old coords end # Test of volume @@ -73,11 +119,11 @@ for i in 1:getnquadpoints(cv) vol += getdetJdV(cv,i) end - @test vol ≈ calculate_volume(func_interpol, x) + @test vol ≈ calculate_volume(func_interpol, coords) # Test quadrature rule after reinit! with ref. coords - x = Ferrite.reference_coordinates(func_interpol) - reinit!(cv, x) + coords = Ferrite.reference_coordinates(func_interpol) + reinit!(cv, coords) vol = 0.0 for i in 1:getnquadpoints(cv) vol += getdetJdV(cv,i) @@ -86,7 +132,7 @@ # Test spatial coordinate (after reinit with ref.coords we should get back the quad_points) for (i, qp_x) in pairs(Ferrite.getpoints(quad_rule)) - @test spatial_coordinate(cv, i, x) ≈ qp_x + @test spatial_coordinate(cv, i, coords) ≈ qp_x end @testset "copy(::CellValues)" begin diff --git a/test/test_facevalues.jl b/test/test_facevalues.jl index 056cc542e2..f70d624ebe 100644 --- a/test/test_facevalues.jl +++ b/test/test_facevalues.jl @@ -13,70 +13,117 @@ for (scalar_interpol, quad_rule) in ( (Lagrange{RefPyramid, 2}(), FaceQuadratureRule{RefPyramid}(2)), (Lagrange{RefPrism, 2}(), FaceQuadratureRule{RefPrism}(2)), ) - - for func_interpol in (scalar_interpol, VectorizedInterpolation(scalar_interpol)) + for func_interpol in (scalar_interpol, VectorizedInterpolation(scalar_interpol)), DiffOrder in 1:2 + (DiffOrder==2 && Ferrite.getorder(func_interpol)==1) && continue #No need to test linear interpolations again geom_interpol = scalar_interpol # Tests below assume this n_basefunc_base = getnbasefunctions(scalar_interpol) + update_gradients = true + update_hessians = (DiffOrder==2 && Ferrite.getorder(func_interpol) > 1) fv = if VERSION ≥ v"1.9" - @inferred FaceValues(quad_rule, func_interpol, geom_interpol) + FaceValues(quad_rule, func_interpol, geom_interpol; update_gradients, update_hessians) else # Type unstable on 1.6, but works at least for 1.9 and later. PR882 - FaceValues(quad_rule, func_interpol, geom_interpol) + FaceValues(quad_rule, func_interpol, geom_interpol; update_gradients, update_hessians) end ndim = Ferrite.getdim(func_interpol) n_basefuncs = getnbasefunctions(func_interpol) @test getnbasefunctions(fv) == n_basefuncs - xs, n = valid_coordinates_and_normals(func_interpol) + coords, n = valid_coordinates_and_normals(func_interpol) for face in 1:Ferrite.nfaces(func_interpol) - reinit!(fv, xs, face) + reinit!(fv, coords, face) @test Ferrite.getcurrentface(fv) == face # We test this by applying a given deformation gradient on all the nodes. # Since this is a linear deformation we should get back the exact values # from the interpolation. - u = Vec{ndim, Float64}[zero(Tensor{1,ndim}) for i in 1:n_basefunc_base] - u_scal = zeros(n_basefunc_base) - H = rand(Tensor{2, ndim}) - V = rand(Tensor{1, ndim}) + V, G, H = if func_interpol isa Ferrite.ScalarInterpolation + (rand(), rand(Tensor{1, ndim}), Tensor{2, ndim}((i,j)-> i==j ? rand() : 0.0)) + else + (rand(Tensor{1, ndim}), rand(Tensor{2, ndim}), Tensor{3, ndim}((i,j,k)-> i==j==k ? rand() : 0.0)) + end + + u_funk(x,V,G,H) = begin + if update_hessians + 0.5*x⋅H⋅x + G⋅x + V + else + G⋅x + V + end + end + + _ue = [u_funk(coords[i],V,G,H) for i in 1:n_basefunc_base] + ue = reinterpret(Float64, _ue) + + for i in 1:getnquadpoints(fv) + xqp = spatial_coordinate(fv, i, coords) + Hqp, Gqp, Vqp = Tensors.hessian(x -> u_funk(x,V,G,H), xqp, :all) + + @test function_value(fv, i, ue) ≈ Vqp + @test function_gradient(fv, i, ue) ≈ Gqp + if update_hessians + #Note, the jacobian of the element is constant, which makes the hessian (of the mapping) + #zero. So this is not the optimal test + @test Ferrite.function_hessian(fv, i, ue) ≈ Hqp + end + if func_interpol isa Ferrite.VectorInterpolation + @test function_symmetric_gradient(fv, i, ue) ≈ 0.5(Gqp + Gqp') + @test function_divergence(fv, i, ue) ≈ tr(Gqp) + ndim == 3 && @test function_curl(fv, i, ue) ≈ Ferrite.curl_from_gradient(Gqp) + else + @test function_divergence(fv, i, ue) ≈ sum(Gqp) + end + end + + #Test CellValues when input is a ::Vector{<:Vec} (most of which is deprecated) + ue_vec = [zero(Vec{ndim,Float64}) for i in 1:n_basefunc_base] + G_vector = rand(Tensor{2, ndim}) for i in 1:n_basefunc_base - u[i] = H ⋅ xs[i] - u_scal[i] = V ⋅ xs[i] + ue_vec[i] = G_vector ⋅ coords[i] end - u_vector = reinterpret(Float64, u) + for i in 1:getnquadpoints(fv) - @test getnormal(fv, i) ≈ n[face] if func_interpol isa Ferrite.ScalarInterpolation - @test function_gradient(fv, i, u) ≈ H - @test function_symmetric_gradient(fv, i, u) ≈ 0.5(H + H') - @test function_divergence(fv, i, u_scal) ≈ sum(V) - @test function_divergence(fv, i, u) ≈ tr(H) - @test function_gradient(fv, i, u_scal) ≈ V - ndim == 3 && @test function_curl(fv, i, u) ≈ Ferrite.curl_from_gradient(H) - function_value(fv, i, u) - function_value(fv, i, u_scal) - else # func_interpol isa Ferrite.VectorInterpolation - @test function_gradient(fv, i, u_vector) ≈ H - @test (@test_deprecated function_gradient(fv, i, u)) ≈ H - @test function_symmetric_gradient(fv, i, u_vector) ≈ 0.5(H + H') - @test (@test_deprecated function_symmetric_gradient(fv, i, u)) ≈ 0.5(H + H') - @test function_divergence(fv, i, u_vector) ≈ tr(H) - @test (@test_deprecated function_divergence(fv, i, u)) ≈ tr(H) + @test function_gradient(fv, i, ue_vec) ≈ G_vector + else# func_interpol isa Ferrite.VectorInterpolation + @test (@test_deprecated function_gradient(fv, i, ue_vec)) ≈ G_vector + @test (@test_deprecated function_symmetric_gradient(fv, i, ue_vec)) ≈ 0.5(G_vector + G_vector') + @test (@test_deprecated function_divergence(fv, i, ue_vec)) ≈ tr(G_vector) if ndim == 3 - @test function_curl(fv, i, u_vector) ≈ Ferrite.curl_from_gradient(H) - @test (@test_deprecated function_curl(fv, i, u)) ≈ Ferrite.curl_from_gradient(H) + @test (@test_deprecated function_curl(fv, i, ue_vec)) ≈ Ferrite.curl_from_gradient(G_vector) + end + function_value(fv, i, ue_vec) + end + end + + #Check if the non-linear mapping is correct + #Only do this for one interpolation becuase it relise on AD on "iterative function" + if scalar_interpol === Lagrange{RefQuadrilateral, 2}() + coords_nl = [x+rand(x)*0.01 for x in coords] #add some displacement to nodes + reinit!(fv, coords_nl, face) + + _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(fv) + xqp = spatial_coordinate(fv, 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(fv, i, ue_nl) ≈ Vqp + @test function_gradient(fv, i, ue_nl) ≈ Gqp + if update_hessians + @test Ferrite.function_hessian(fv, i, ue_nl) ≈ Hqp end - @test function_value(fv, i, u_vector) ≈ (@test_deprecated function_value(fv, i, u)) end + reinit!(fv, coords, face) # reinit back to old coords end + # Test of volume vol = 0.0 for i in 1:getnquadpoints(fv) vol += getdetJdV(fv,i) end let ip_base = func_interpol isa VectorizedInterpolation ? func_interpol.ip : func_interpol - x_face = xs[[Ferrite.facedof_indices(ip_base)[face]...]] + x_face = coords[[Ferrite.facedof_indices(ip_base)[face]...]] @test vol ≈ calculate_face_area(ip_base, x_face, face) end diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index 25c6905154..41fbc539dc 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -225,4 +225,19 @@ @test Ferrite.is_discontinuous(d_ip_t) == true end +@testset "Correctness of AD of embedded interpolations" begin + ip = Lagrange{RefHexahedron,2}()^3 + ξ = rand(Vec{3,Float64}) + for I in 1:getnbasefunctions(ip) + #Call StaticArray-version + H_sa, G_sa, V_sa = Ferrite._shape_hessian_gradient_and_value_static_array(ip, ξ, I) + #Call tensor AD version + H, G, V = Ferrite.shape_hessian_gradient_and_value(ip, ξ, I) + + @test V ≈ V_sa + @test G ≈ G_sa + @test H ≈ H_sa + end +end + end # testset diff --git a/test/test_utils.jl b/test/test_utils.jl index 140954df8b..7305eaa411 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -294,3 +294,41 @@ module DummyRefShapes ) end end + +############################################################ +# Inverse parametric mapping ξ = ϕ(x) for testing hessians # +############################################################ +function function_value_from_physical_coord(interpolation::Interpolation, cell_coordinates, X::Vec{dim,T}, ue) where {dim,T} + n_basefuncs = getnbasefunctions(interpolation) + scalar_ip = interpolation isa Ferrite.ScalarInterpolation ? interpolation : interpolation.ip + @assert length(ue) == n_basefuncs + ξ = MAPPING(scalar_ip, cell_coordinates, X) + u = zero(typeof(shape_value(interpolation, ξ, 1))) #Is there a utility function for this init? + for j in 1:n_basefuncs + N = shape_value(interpolation, ξ, j) + u += N * ue[j] + end + return u +end + +function MAPPING(interpolation, cell_coordinates, global_coordinate::Vec{dim,T}) where {dim,T} + ξ = zero(global_coordinate) + n_basefuncs = getnbasefunctions(interpolation) + max_iters = 10 + tol_norm = 1e-16 + for _ in 1:max_iters + global_guess = zero(global_coordinate) + J = zero(Tensor{2,dim,T}) + for j in 1:n_basefuncs + dNdξ, N = Ferrite.shape_gradient_and_value(interpolation, ξ, j) + global_guess += N * cell_coordinates[j] + J += cell_coordinates[j] ⊗ dNdξ + end + residual = global_guess - global_coordinate + if norm(residual) <= tol_norm + break + end + ξ -= inv(J) ⋅ residual + end + return ξ +end \ No newline at end of file