From b2c4758736e77c28b2441944831f4761b69e0c6f Mon Sep 17 00:00:00 2001 From: lijas Date: Fri, 24 May 2024 13:51:36 +0200 Subject: [PATCH 01/16] Add higher order gradients in Fe-values --- docs/src/topics/FEValues.md | 53 +++++++++++++++ src/FEValues/CellValues.jl | 16 +++-- src/FEValues/FacetValues.jl | 15 +++-- src/FEValues/FunctionValues.jl | 97 ++++++++++++++++++++------ src/FEValues/GeometryMapping.jl | 30 +++++---- src/FEValues/common_values.jl | 37 +++++++++- src/Ferrite.jl | 2 +- src/exports.jl | 1 + src/interpolations.jl | 42 ++++++++++-- test/runtests.jl | 3 +- test/test_cellvalues.jl | 116 ++++++++++++++++++++++---------- test/test_facevalues.jl | 113 ++++++++++++++++++++++--------- test/test_interpolations.jl | 15 +++++ test/test_utils.jl | 38 +++++++++++ 14 files changed, 457 insertions(+), 121 deletions(-) diff --git a/docs/src/topics/FEValues.md b/docs/src/topics/FEValues.md index df54430cbc..a6a48c0147 100644 --- a/docs/src/topics/FEValues.md +++ b/docs/src/topics/FEValues.md @@ -42,6 +42,59 @@ 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: + ```math + \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: + + ```math + \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} + ``` + + Using the fact that $\frac{\mathrm{d}}{\mathrm{d}x_j} = J^{-1}_{sj}\frac{\mathrm{d}}{\mathrm{d}\xi_s}$, the first term in equation the equation above can be expressed as: + + ```math + \begin{align*} + \frac{\mathrm{d}}{\mathrm{d}x_j}(\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}) J^{-1}_{ri} = J^{-1}_{sj}\frac{\mathrm{d}}{\mathrm{d}\xi_s}(\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}) J^{-1}_{ri} = J^{-1}_{sj}\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: + + ```math + \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}]J^{-1}_{sj} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}[ - J^{-1}_{rk}\mathcal{H}_{kps} J^{-1}_{pi}]J^{-1}_{sj} = - \frac{\mathrm{d} \hat N}{\mathrm{d} x_k}\mathcal{H}_{kps} J^{-1}_{pi}J^{-1}_{sj} + \end{align*} + ``` + + where we have used that the inverse of the jacobian can be computed as: + + ```math + \begin{align*} + 0 = \frac{\mathrm{d}}{\mathrm{d}\xi_s} (J_{kr} J^{-1}_{ri} ) = \frac{\mathrm{d}J_{kp}}{\mathrm{d}\xi_s} J^{-1}_{pi} + J_{kr} \frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}\xi_s} = 0 \quad \Rightarrow \\ + \end{align*} + ``` + + ```math + \begin{align*} + \frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}\xi_s} = - J^{-1}_{rk}\frac{\mathrm{d}J_{kp}}{\mathrm{d}\xi_s} J^{-1}_{pi} = - J^{-1}_{rk}\mathcal{H}_{kps} J^{-1}_{pi}\\ + \end{align*} + ``` + ### Covariant Piola mapping, H(curl) `Ferrite.CovariantPiolaMapping` diff --git a/src/FEValues/CellValues.jl b/src/FEValues/CellValues.jl index 209a36d7d2..b00b9bcd4c 100644 --- a/src/FEValues/CellValues.jl +++ b/src/FEValues/CellValues.jl @@ -41,12 +41,19 @@ struct CellValues{FV, GM, QR, detT} <: AbstractCellValues qr::QR # QuadratureRule 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 +function CellValues(::Type{T}, qr::QuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation; + 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 = _update_hessians ? 2 : (_update_gradients ? 1 : 0) 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) detJdV = _update_detJdV ? fill(T(NaN), length(getweights(qr))) : nothing @@ -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/FacetValues.jl b/src/FEValues/FacetValues.jl index 52d315cba8..ed1184f0c8 100644 --- a/src/FEValues/FacetValues.jl +++ b/src/FEValues/FacetValues.jl @@ -40,13 +40,17 @@ struct FacetValues{FV, GM, FQR, detT, nT, V_FV<:AbstractVector{FV}, V_GM<:Abstra current_facet::ScalarWrapper{Int} end -function FacetValues(::Type{T}, fqr::FacetQuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation{sdim} = default_geometric_interpolation(ip_fun); - update_gradients::Union{Bool,Nothing} = nothing) where {T,sdim} +function FacetValues(::Type{T}, fqr::FacetQuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation{sdim} = default_geometric_interpolation(ip_fun); + update_gradients::Union{Bool,Nothing} = nothing, + update_hessians ::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 = update_gradients === nothing ? true : update_gradients + _update_hessians = update_hessians === nothing ? false : update_hessians + _update_hessians && @assert _update_gradients + + FunDiffOrder = _update_hessians ? 2 : (_update_gradients ? 1 : 0) 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 FacetValues + geo_mapping = [GeometryMapping{GeoDiffOrder}(T, ip_geo.ip, qr) for qr in fqr.face_rules]::Vector{<:GeometryMapping{GeoDiffOrder}} fun_values = [FunctionValues{FunDiffOrder}(T, ip_fun, qr, ip_geo) for qr in fqr.face_rules]::Vector{<:FunctionValues{FunDiffOrder}} max_nquadpoints = maximum(qr->length(getweights(qr)), fqr.face_rules) @@ -84,6 +88,7 @@ get_fun_values(fv::FacetValues) = @inbounds fv.fun_values[getcurrentfacet(fv)] @propagate_inbounds shape_value(fv::FacetValues, q_point::Int, i::Int) = shape_value(get_fun_values(fv), q_point, i) @propagate_inbounds shape_gradient(fv::FacetValues, q_point::Int, i::Int) = shape_gradient(get_fun_values(fv), q_point, i) +@propagate_inbounds shape_hessian(fv::FacetValues, q_point::Int, i::Int) = shape_hessian(get_fun_values(fv), q_point, i) @propagate_inbounds shape_symmetric_gradient(fv::FacetValues, 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 0e664dacea..64c9d868ef 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,17 +43,22 @@ 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} @@ -52,17 +67,23 @@ function FunctionValues{DiffOrder}(::Type{T}, ip::Interpolation, qr::QuadratureR 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 - throw(ArgumentError("Currently only values and gradients can be updated in FunctionValues")) end - fv = FunctionValues(ip, Nx, Nξ, dNdx, dNdξ) + 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, gradients, and hessians can be updated in FunctionValues")) + end + + 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 +94,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 +119,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 +197,29 @@ end return nothing end -# TODO in PR798, apply_mapping! for +@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 + @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 5821d611b4..09af19a962 100644 --- a/src/FEValues/GeometryMapping.jl +++ b/src/FEValues/GeometryMapping.jl @@ -13,8 +13,8 @@ struct MappingValues{JT, HT} H::HT # dJ/dξ # Hessian end -@inline getjacobian(mv::MappingValues{<:Union{AbstractTensor, SMatrix}}) = mv.J -# @inline gethessian(mv::MappingValues{<:Any,<:AbstractTensor}) = mv.H # PR798 +@inline getjacobian(mv::MappingValues{<:Union{AbstractTensor, SMatrix}}) = mv.J +@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( - ip::IP, M::M_t, dMdξ::dMdξ_t, d2Mdξ2::d2Mdξ2_t) where - {IP <: ScalarInterpolation, M_t<:AbstractMatrix{<:Number}, + 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,15 +134,17 @@ 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ξ))) + 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] 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 c7a694a748..841d239910 100644 --- a/src/FEValues/common_values.jl +++ b/src/FEValues/common_values.jl @@ -211,6 +211,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]) @@ -308,10 +343,8 @@ 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/Ferrite.jl b/src/Ferrite.jl index 2b0a286522..32504f07cb 100644 --- a/src/Ferrite.jl +++ b/src/Ferrite.jl @@ -17,7 +17,7 @@ using OrderedCollections: using SparseArrays: SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals, sparse, spzeros using StaticArrays: - StaticArrays, MMatrix, SMatrix, SVector + StaticArrays, MArray, MMatrix, SArray, SMatrix, SVector using WriteVTK: WriteVTK, VTKCellTypes using Tensors: diff --git a/src/exports.jl b/src/exports.jl index 0e834f88a6..95029d1010 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 4cad56d16b..4545ac601e 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -182,7 +182,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) @@ -197,7 +196,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) @@ -232,7 +231,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) @@ -242,7 +240,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) @@ -1540,6 +1538,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/runtests.jl b/test/runtests.jl index ba74624411..665503d941 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -34,9 +34,10 @@ end include("test_utils.jl") # Unit tests -include("test_interpolations.jl") +#include("test_interpolations.jl") include("test_cellvalues.jl") include("test_facevalues.jl") +asdf include("test_interfacevalues.jl") include("test_quadrules.jl") include("test_assemble.jl") diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index 5c1c3cf0a0..0d1ec1bb0a 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) rdim = Ferrite.getrefdim(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 = zeros(Vec{rdim, Float64}, n_basefunc_base) - u_scal = zeros(n_basefunc_base) - H = rand(Tensor{2, rdim}) - V = rand(Tensor{1, rdim}) + V, G, H = if func_interpol isa Ferrite.ScalarInterpolation + (rand(), rand(Tensor{1, rdim}), Tensor{2, rdim}((i,j)-> i==j ? rand() : 0.0)) + else + (rand(Tensor{1, rdim}), rand(Tensor{2, rdim}), Tensor{3, rdim}((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) + rdim == 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{rdim,Float64}) for i in 1:n_basefunc_base] + G_vector = rand(Tensor{2, rdim}) 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 - rdim == 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 rdim == 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 43cff3d366..d6b6f525c0 100644 --- a/test/test_facevalues.jl +++ b/test/test_facevalues.jl @@ -13,70 +13,117 @@ for (scalar_interpol, quad_rule) in ( (Lagrange{RefPyramid, 2}(), FacetQuadratureRule{RefPyramid}(2)), (Lagrange{RefPrism, 2}(), FacetQuadratureRule{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 FacetValues(quad_rule, func_interpol, geom_interpol) + FacetValues(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 - FacetValues(quad_rule, func_interpol, geom_interpol) + FacetValues(quad_rule, func_interpol, geom_interpol; update_gradients, update_hessians) end rdim = Ferrite.getrefdim(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.nfacets(func_interpol) - reinit!(fv, xs, face) + reinit!(fv, coords, face) @test Ferrite.getcurrentfacet(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 = zeros(Vec{rdim, Float64}, n_basefunc_base) - u_scal = zeros(n_basefunc_base) - H = rand(Tensor{2, rdim}) - V = rand(Tensor{1, rdim}) + V, G, H = if func_interpol isa Ferrite.ScalarInterpolation + (rand(), rand(Tensor{1, rdim}), Tensor{2, rdim}((i,j)-> i==j ? rand() : 0.0)) + else + (rand(Tensor{1, rdim}), rand(Tensor{2, rdim}), Tensor{3, rdim}((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) + rdim == 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{rdim,Float64}) for i in 1:n_basefunc_base] + G_vector = rand(Tensor{2, rdim}) 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 - rdim == 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 rdim == 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.facetdof_indices(ip_base)[face]...]] + x_face = coords[[Ferrite.facetdof_indices(ip_base)[face]...]] @test vol ≈ calculate_facet_area(ip_base, x_face, face) end diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index 63b0b5c61c..d3bc685c93 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -219,4 +219,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 c0fb2e9e83..1ed3db5462 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 From a5851676f079396b403e83a3971a6904cd812379 Mon Sep 17 00:00:00 2001 From: lijas Date: Fri, 24 May 2024 14:24:49 +0200 Subject: [PATCH 02/16] some clean up and update changelog --- CHANGELOG.md | 5 +++++ docs/src/topics/FEValues.md | 2 +- src/FEValues/GeometryMapping.jl | 2 +- test/runtests.jl | 1 - 4 files changed, 7 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 456316478c..5df24694ee 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -317,6 +317,11 @@ more discussion). a given grid (based on its node coordinates), and returns the minimum and maximum vertices of the bounding box. ([#880][github-880]) +- CellValues and FacetValues can now store and map second order gradients (Hessians). The number + of gradients computed in CellValues/FacetValues is specified using the keyword arguments + `update_gradients::Bool` (default true) and `update_hessians::Bool` (default false) in the + constructors, i.e. `CellValues(...; update_hessians=true)`. + ### Changed - `create_sparsity_pattern` now supports cross-element dof coupling by passing kwarg diff --git a/docs/src/topics/FEValues.md b/docs/src/topics/FEValues.md index a6a48c0147..9c6aaee84c 100644 --- a/docs/src/topics/FEValues.md +++ b/docs/src/topics/FEValues.md @@ -65,7 +65,7 @@ Second order gradients of the shape functions are computed as \end{align} ``` - Using the fact that $\frac{\mathrm{d}}{\mathrm{d}x_j} = J^{-1}_{sj}\frac{\mathrm{d}}{\mathrm{d}\xi_s}$, the first term in equation the equation above can be expressed as: + Using the fact that $\frac{\mathrm{d}}{\mathrm{d}x_j} = J^{-1}_{sj}\frac{\mathrm{d}}{\mathrm{d}\xi_s}$, the first term in the equation above can be expressed as: ```math \begin{align*} diff --git a/src/FEValues/GeometryMapping.jl b/src/FEValues/GeometryMapping.jl index 09af19a962..7a711838b1 100644 --- a/src/FEValues/GeometryMapping.jl +++ b/src/FEValues/GeometryMapping.jl @@ -35,7 +35,7 @@ struct GeometryMapping{DiffOrder, IP, M_t, dMdξ_t, d2Mdξ2_t} M::M_t # ::AbstractMatrix{<:Number} Values of geometric shape functions dMdξ::dMdξ_t # ::AbstractMatrix{<:Vec} Gradients of geometric shape functions in ref-domain d2Mdξ2::d2Mdξ2_t # ::AbstractMatrix{<:Tensor{2}} Hessians of geometric shape functions in ref-domain - # ::Nothing When not required + function GeometryMapping( ip::IP, M::M_t, ::Nothing, ::Nothing ) where {IP <: ScalarInterpolation, M_t<:AbstractMatrix{<:Number}} diff --git a/test/runtests.jl b/test/runtests.jl index 665503d941..b5898bbc2f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -37,7 +37,6 @@ include("test_utils.jl") #include("test_interpolations.jl") include("test_cellvalues.jl") include("test_facevalues.jl") -asdf include("test_interfacevalues.jl") include("test_quadrules.jl") include("test_assemble.jl") From 6bcf7d4565c5e01600b053a09e7bd896f6fd71ea Mon Sep 17 00:00:00 2001 From: lijas Date: Fri, 24 May 2024 15:40:47 +0200 Subject: [PATCH 03/16] derp --- src/FEValues/FunctionValues.jl | 2 +- test/runtests.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/FEValues/FunctionValues.jl b/src/FEValues/FunctionValues.jl index 64c9d868ef..8629b71838 100644 --- a/src/FEValues/FunctionValues.jl +++ b/src/FEValues/FunctionValues.jl @@ -201,7 +201,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 diff --git a/test/runtests.jl b/test/runtests.jl index b5898bbc2f..ba74624411 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -34,7 +34,7 @@ end include("test_utils.jl") # Unit tests -#include("test_interpolations.jl") +include("test_interpolations.jl") include("test_cellvalues.jl") include("test_facevalues.jl") include("test_interfacevalues.jl") From e77886f62b896967877f4d3d45a26a5fe03c4013 Mon Sep 17 00:00:00 2001 From: lijas Date: Mon, 27 May 2024 14:34:12 +0200 Subject: [PATCH 04/16] Update docs/src/topics/FEValues.md Co-authored-by: Knut Andreas Meyer --- docs/src/topics/FEValues.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/topics/FEValues.md b/docs/src/topics/FEValues.md index 9c6aaee84c..5703cd7681 100644 --- a/docs/src/topics/FEValues.md +++ b/docs/src/topics/FEValues.md @@ -77,7 +77,7 @@ Second order gradients of the shape functions are computed as ```math \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}]J^{-1}_{sj} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}[ - J^{-1}_{rk}\mathcal{H}_{kps} J^{-1}_{pi}]J^{-1}_{sj} = - \frac{\mathrm{d} \hat N}{\mathrm{d} x_k}\mathcal{H}_{kps} J^{-1}_{pi}J^{-1}_{sj} + \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}x_j} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\left[\frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}\xi_s}\right]J^{-1}_{sj} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\left[- J^{-1}_{rk}\mathcal{H}_{kps} J^{-1}_{pi}\right] J^{-1}_{sj} = - \frac{\mathrm{d} \hat N}{\mathrm{d} x_k}\mathcal{H}_{kps} J^{-1}_{pi}J^{-1}_{sj} \end{align*} ``` From 0c4c63dc395da8a33b506379d3d7b2b94dbb4b36 Mon Sep 17 00:00:00 2001 From: lijas Date: Mon, 27 May 2024 14:39:37 +0200 Subject: [PATCH 05/16] Apply suggestions from code review Co-authored-by: Knut Andreas Meyer --- docs/src/topics/FEValues.md | 6 +++--- src/FEValues/GeometryMapping.jl | 2 +- src/FEValues/common_values.jl | 3 --- test/test_cellvalues.jl | 4 ++++ test/test_interpolations.jl | 13 +++++++++++++ test/test_utils.jl | 2 +- 6 files changed, 22 insertions(+), 8 deletions(-) diff --git a/docs/src/topics/FEValues.md b/docs/src/topics/FEValues.md index 5703cd7681..8d71609927 100644 --- a/docs/src/topics/FEValues.md +++ b/docs/src/topics/FEValues.md @@ -61,15 +61,15 @@ Second order gradients of the shape functions are computed as ```math \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}) + \frac{\mathrm{d}^2 N}{\mathrm{d}x_i \mathrm{d}x_j} = \frac{\mathrm{d}}{\mathrm{d}x_j}\left[\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\right] J^{-1}_{ri} + \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r} \frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}x_j} \end{align} ``` - Using the fact that $\frac{\mathrm{d}}{\mathrm{d}x_j} = J^{-1}_{sj}\frac{\mathrm{d}}{\mathrm{d}\xi_s}$, the first term in the equation above can be expressed as: + Using the fact that $\frac{\mathrm{d}\hat{f}(\boldsymbol{\xi})}{\mathrm{d}x_j} = \frac{\mathrm{d}\hat{f}(\boldsymbol{\xi})}{\mathrm{d}\xi_s} J^{-1}_{sj}$, the first term in the equation above can be expressed as: ```math \begin{align*} - \frac{\mathrm{d}}{\mathrm{d}x_j}(\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}) J^{-1}_{ri} = J^{-1}_{sj}\frac{\mathrm{d}}{\mathrm{d}\xi_s}(\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}) J^{-1}_{ri} = J^{-1}_{sj}\frac{\mathrm{d}^2 \hat N}{\mathrm{d} \xi_s\mathrm{d} \xi_r} J^{-1}_{ri} + \frac{\mathrm{d}}{\mathrm{d}x_j}\left[\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\right] J^{-1}_{ri} = J^{-1}_{sj}\frac{\mathrm{d}}{\mathrm{d}\xi_s}\left[\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\right] J^{-1}_{ri} = J^{-1}_{sj}\left[\frac{\mathrm{d}^2 \hat N}{\mathrm{d} \xi_s\mathrm{d} \xi_r}\right] J^{-1}_{ri} \end{align*} ``` diff --git a/src/FEValues/GeometryMapping.jl b/src/FEValues/GeometryMapping.jl index 7a711838b1..74b86cdd94 100644 --- a/src/FEValues/GeometryMapping.jl +++ b/src/FEValues/GeometryMapping.jl @@ -35,7 +35,7 @@ struct GeometryMapping{DiffOrder, IP, M_t, dMdξ_t, d2Mdξ2_t} M::M_t # ::AbstractMatrix{<:Number} Values of geometric shape functions dMdξ::dMdξ_t # ::AbstractMatrix{<:Vec} Gradients of geometric shape functions in ref-domain d2Mdξ2::d2Mdξ2_t # ::AbstractMatrix{<:Tensor{2}} Hessians of geometric shape functions in ref-domain - + # ::Nothing (dMdξ or d2Mdξ2 if not required) function GeometryMapping( ip::IP, M::M_t, ::Nothing, ::Nothing ) where {IP <: ScalarInterpolation, M_t<:AbstractMatrix{<:Number}} diff --git a/src/FEValues/common_values.jl b/src/FEValues/common_values.jl index 841d239910..59894f7198 100644 --- a/src/FEValues/common_values.jl +++ b/src/FEValues/common_values.jl @@ -242,9 +242,6 @@ 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]) diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index 0d1ec1bb0a..90e7422752 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -23,6 +23,10 @@ update_gradients = true update_hessians = (DiffOrder==2 && Ferrite.getorder(func_interpol) > 1) cv = CellValues(quad_rule, func_interpol, geom_interpol; update_gradients, update_hessians) + if update_gradients && !update_hessians # Check correct and type-stable default constructor + cv_default = @inferred CellValues(quad_rule, func_interpol, geom_interpol) + @test typeof(cv) === typeof(cv_default) + end rdim = Ferrite.getrefdim(func_interpol) n_basefuncs = getnbasefunctions(func_interpol) diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index d3bc685c93..1c1d1bec09 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -232,6 +232,19 @@ end @test G ≈ G_sa @test H ≈ H_sa end + + ips = Lagrange{RefQuadrilateral,2}() + vdim = 3 + ipv = ips^vdim + ξ = rand(Vec{2, Float64}) + for ipv_ind in 1:getnbasefunctions(ipv) + ips_ind, v_ind = fldmod1(ipv_ind, vdim) + H, G, V = Ferrite.shape_hessian_gradient_and_value(ipv, ξ, ipv_ind) + h, g, v = Ferrite.shape_hessian_gradient_and_value(ips, ξ, ips_ind) + @test h ≈ H[v_ind, :, :] + @test g ≈ G[v_ind, :] + @test v ≈ V[v_ind] + end end end # testset diff --git a/test/test_utils.jl b/test/test_utils.jl index 1ed3db5462..399204d6b1 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -303,7 +303,7 @@ function function_value_from_physical_coord(interpolation::Interpolation, cell_c 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? + u = zero(shape_value(interpolation, ξ, 1)) for j in 1:n_basefuncs N = shape_value(interpolation, ξ, j) u += N * ue[j] From 7224a6bfa392a02fad6e56fada11dd2e7d9074a4 Mon Sep 17 00:00:00 2001 From: lijas Date: Mon, 27 May 2024 15:38:32 +0200 Subject: [PATCH 06/16] update some doc-strings --- src/FEValues/CellValues.jl | 9 +++++++-- src/FEValues/FacetValues.jl | 7 ++++++- test/test_facevalues.jl | 10 ++++++---- 3 files changed, 19 insertions(+), 7 deletions(-) diff --git a/src/FEValues/CellValues.jl b/src/FEValues/CellValues.jl index b00b9bcd4c..484fa9dc3f 100644 --- a/src/FEValues/CellValues.jl +++ b/src/FEValues/CellValues.jl @@ -1,5 +1,5 @@ """ - CellValues([::Type{T},] quad_rule::QuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation]) + CellValues([::Type{T},] quad_rule::QuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation]; update_detJdV=true, update_gradients=true, update_hessians=false) A `CellValues` object facilitates the process of evaluating values of shape functions, gradients of shape functions, values of nodal functions, gradients and divergences of nodal functions etc. in the finite element cell. @@ -12,6 +12,11 @@ values of nodal functions, gradients and divergences of nodal functions etc. in By default linear Lagrange interpolation is used. For embedded elements the geometric interpolations should be vectorized to the spatial dimension. +**Keyword arguments:** +* `update_gradients`: Specifies if the gradients of the shape functions should be updated (default true) +* `update_hessians`: Specifies if the hessians of the shape functions should be updated (default false) +* `update_detJdV`: Specifies if the volume associated with each quadrature point should be updated (default true) + **Common methods:** * [`reinit!`](@ref) @@ -42,7 +47,7 @@ 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_gradients ::Union{Bool,Nothing} = nothing, #Use Union{Bool,Nothing} to get type-stable code update_hessians ::Union{Bool,Nothing} = nothing, update_detJdV ::Union{Bool,Nothing} = nothing) where T diff --git a/src/FEValues/FacetValues.jl b/src/FEValues/FacetValues.jl index ed1184f0c8..212f0692c4 100644 --- a/src/FEValues/FacetValues.jl +++ b/src/FEValues/FacetValues.jl @@ -1,5 +1,5 @@ """ - FacetValues([::Type{T}], quad_rule::FacetQuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation]) + FacetValues([::Type{T}], quad_rule::FacetQuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation]; ; update_gradients=true, update_hessians=false, ) A `FacetValues` object facilitates the process of evaluating values of shape functions, gradients of shape functions, values of nodal functions, gradients and divergences of nodal functions etc. on the facets of finite elements. @@ -12,6 +12,11 @@ values of nodal functions, gradients and divergences of nodal functions etc. on * `geom_interpol`: an optional instance of an [`Interpolation`](@ref) which is used to interpolate the geometry. By default linear Lagrange interpolation is used. +**Keyword arguments:** + +* `update_gradients`: Specifies if the gradients of the shape functions should be updated (default true) +* `update_hessians`: Specifies if the hessians of the shape functions should be updated (default false) + **Common methods:** * [`reinit!`](@ref) diff --git a/test/test_facevalues.jl b/test/test_facevalues.jl index d6b6f525c0..32fbbbdf9d 100644 --- a/test/test_facevalues.jl +++ b/test/test_facevalues.jl @@ -19,11 +19,13 @@ for (scalar_interpol, quad_rule) in ( n_basefunc_base = getnbasefunctions(scalar_interpol) update_gradients = true update_hessians = (DiffOrder==2 && Ferrite.getorder(func_interpol) > 1) - fv = if VERSION ≥ v"1.9" - FacetValues(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 - FacetValues(quad_rule, func_interpol, geom_interpol; update_gradients, update_hessians) + + fv = FacetValues(quad_rule, func_interpol, geom_interpol; update_gradients, update_hessians) + if update_gradients && !update_hessians && VERSION ≥ v"1.9" # Check correct and type-stable default constructor + fv_default = FacetValues(quad_rule, func_interpol, geom_interpol) + @test typeof(fv) === typeof(fv_default) end + rdim = Ferrite.getrefdim(func_interpol) n_basefuncs = getnbasefunctions(func_interpol) From 8c14a53b389908f648df7dd8234365260be83860 Mon Sep 17 00:00:00 2001 From: lijas Date: Mon, 27 May 2024 16:04:34 +0200 Subject: [PATCH 07/16] fix linting errors --- CHANGELOG.md | 2 +- docs/src/topics/FEValues.md | 14 +++++++------- src/FEValues/CellValues.jl | 6 +++--- src/FEValues/FacetValues.jl | 4 ++-- src/FEValues/FunctionValues.jl | 4 ++-- src/FEValues/GeometryMapping.jl | 8 ++++---- src/FEValues/common_values.jl | 2 +- test/test_cellvalues.jl | 4 ++-- test/test_facevalues.jl | 2 +- test/test_utils.jl | 2 +- 10 files changed, 24 insertions(+), 24 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5df24694ee..59d350ae5e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -320,7 +320,7 @@ more discussion). - CellValues and FacetValues can now store and map second order gradients (Hessians). The number of gradients computed in CellValues/FacetValues is specified using the keyword arguments `update_gradients::Bool` (default true) and `update_hessians::Bool` (default false) in the - constructors, i.e. `CellValues(...; update_hessians=true)`. + constructors, i.e. `CellValues(...; update_hessians=true)`. ### Changed diff --git a/docs/src/topics/FEValues.md b/docs/src/topics/FEValues.md index 8d71609927..47764c3d67 100644 --- a/docs/src/topics/FEValues.md +++ b/docs/src/topics/FEValues.md @@ -45,14 +45,14 @@ For scalar fields, we always use scalar base functions. For tensorial fields (no Second order gradients of the shape functions are computed as ```math -\begin{align*} +\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: ```math - \begin{align*} + \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*} ``` @@ -60,7 +60,7 @@ Second order gradients of the shape functions are computed as For the second order gradients, we first use the product rule on the equation above: ```math - \begin{align} + \begin{align} \frac{\mathrm{d}^2 N}{\mathrm{d}x_i \mathrm{d}x_j} = \frac{\mathrm{d}}{\mathrm{d}x_j}\left[\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\right] J^{-1}_{ri} + \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r} \frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}x_j} \end{align} ``` @@ -68,7 +68,7 @@ Second order gradients of the shape functions are computed as Using the fact that $\frac{\mathrm{d}\hat{f}(\boldsymbol{\xi})}{\mathrm{d}x_j} = \frac{\mathrm{d}\hat{f}(\boldsymbol{\xi})}{\mathrm{d}\xi_s} J^{-1}_{sj}$, the first term in the equation above can be expressed as: ```math - \begin{align*} + \begin{align*} \frac{\mathrm{d}}{\mathrm{d}x_j}\left[\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\right] J^{-1}_{ri} = J^{-1}_{sj}\frac{\mathrm{d}}{\mathrm{d}\xi_s}\left[\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\right] J^{-1}_{ri} = J^{-1}_{sj}\left[\frac{\mathrm{d}^2 \hat N}{\mathrm{d} \xi_s\mathrm{d} \xi_r}\right] J^{-1}_{ri} \end{align*} ``` @@ -76,7 +76,7 @@ Second order gradients of the shape functions are computed as The second term can be written as: ```math - \begin{align*} + \begin{align*} \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}x_j} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\left[\frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}\xi_s}\right]J^{-1}_{sj} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\left[- J^{-1}_{rk}\mathcal{H}_{kps} J^{-1}_{pi}\right] J^{-1}_{sj} = - \frac{\mathrm{d} \hat N}{\mathrm{d} x_k}\mathcal{H}_{kps} J^{-1}_{pi}J^{-1}_{sj} \end{align*} ``` @@ -84,13 +84,13 @@ Second order gradients of the shape functions are computed as where we have used that the inverse of the jacobian can be computed as: ```math - \begin{align*} + \begin{align*} 0 = \frac{\mathrm{d}}{\mathrm{d}\xi_s} (J_{kr} J^{-1}_{ri} ) = \frac{\mathrm{d}J_{kp}}{\mathrm{d}\xi_s} J^{-1}_{pi} + J_{kr} \frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}\xi_s} = 0 \quad \Rightarrow \\ \end{align*} ``` ```math - \begin{align*} + \begin{align*} \frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}\xi_s} = - J^{-1}_{rk}\frac{\mathrm{d}J_{kp}}{\mathrm{d}\xi_s} J^{-1}_{pi} = - J^{-1}_{rk}\mathcal{H}_{kps} J^{-1}_{pi}\\ \end{align*} ``` diff --git a/src/FEValues/CellValues.jl b/src/FEValues/CellValues.jl index 484fa9dc3f..e29e5e42b1 100644 --- a/src/FEValues/CellValues.jl +++ b/src/FEValues/CellValues.jl @@ -47,9 +47,9 @@ 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, #Use Union{Bool,Nothing} to get type-stable code - update_hessians ::Union{Bool,Nothing} = nothing, - update_detJdV ::Union{Bool,Nothing} = nothing) where T + update_gradients ::Union{Bool,Nothing} = nothing, #Use Union{Bool,Nothing} to get type-stable code + update_hessians ::Union{Bool,Nothing} = nothing, + update_detJdV ::Union{Bool,Nothing} = nothing) where T _update_gradients = update_gradients === nothing ? true : update_gradients _update_hessians = update_hessians === nothing ? false : update_hessians diff --git a/src/FEValues/FacetValues.jl b/src/FEValues/FacetValues.jl index 212f0692c4..d712c0de5b 100644 --- a/src/FEValues/FacetValues.jl +++ b/src/FEValues/FacetValues.jl @@ -45,8 +45,8 @@ struct FacetValues{FV, GM, FQR, detT, nT, V_FV<:AbstractVector{FV}, V_GM<:Abstra current_facet::ScalarWrapper{Int} end -function FacetValues(::Type{T}, fqr::FacetQuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation{sdim} = default_geometric_interpolation(ip_fun); - update_gradients::Union{Bool,Nothing} = nothing, +function FacetValues(::Type{T}, fqr::FacetQuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation{sdim} = default_geometric_interpolation(ip_fun); + update_gradients::Union{Bool,Nothing} = nothing, update_hessians ::Union{Bool,Nothing} = nothing) where {T,sdim} _update_gradients = update_gradients === nothing ? true : update_gradients diff --git a/src/FEValues/FunctionValues.jl b/src/FEValues/FunctionValues.jl index 8629b71838..86f9b8b0dd 100644 --- a/src/FEValues/FunctionValues.jl +++ b/src/FEValues/FunctionValues.jl @@ -208,7 +208,7 @@ end 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 ? + 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 @@ -220,6 +220,6 @@ end return nothing end -# TODO in PR798, apply_mapping! for +# TODO in PR798, apply_mapping! for # * CovariantPiolaMapping # * ContravariantPiolaMapping diff --git a/src/FEValues/GeometryMapping.jl b/src/FEValues/GeometryMapping.jl index 74b86cdd94..476d3a5f2f 100644 --- a/src/FEValues/GeometryMapping.jl +++ b/src/FEValues/GeometryMapping.jl @@ -13,8 +13,8 @@ struct MappingValues{JT, HT} H::HT # dJ/dξ # Hessian end -@inline getjacobian(mv::MappingValues{<:Union{AbstractTensor, SMatrix}}) = mv.J -@inline gethessian(mv::MappingValues{<:Any,<:AbstractTensor}) = mv.H +@inline getjacobian(mv::MappingValues{<:Union{AbstractTensor, SMatrix}}) = mv.J +@inline gethessian(mv::MappingValues{<:Any,<:AbstractTensor}) = mv.H """ @@ -47,8 +47,8 @@ struct GeometryMapping{DiffOrder, IP, M_t, dMdξ_t, d2Mdξ2_t} return new{1, IP, M_t, dMdξ_t, Nothing}(ip, M, dMdξ, nothing) end function GeometryMapping( - ip::IP, M::M_t, dMdξ::dMdξ_t, d2Mdξ2::d2Mdξ2_t) where - {IP <: ScalarInterpolation, M_t<:AbstractMatrix{<:Number}, + 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 diff --git a/src/FEValues/common_values.jl b/src/FEValues/common_values.jl index 59894f7198..88a117cc68 100644 --- a/src/FEValues/common_values.jl +++ b/src/FEValues/common_values.jl @@ -215,7 +215,7 @@ 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. + 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) diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index 90e7422752..2b29b44f73 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -45,7 +45,7 @@ (rand(Tensor{1, rdim}), rand(Tensor{2, rdim}), Tensor{3, rdim}((i,j,k)-> i==j==k ? rand() : 0.0)) end - u_funk(x,V,G,H) = begin + u_funk(x,V,G,H) = begin if update_hessians 0.5*x⋅H⋅x + G⋅x + V else @@ -63,7 +63,7 @@ @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) + #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 diff --git a/test/test_facevalues.jl b/test/test_facevalues.jl index 32fbbbdf9d..6321d5ff0a 100644 --- a/test/test_facevalues.jl +++ b/test/test_facevalues.jl @@ -63,7 +63,7 @@ for (scalar_interpol, quad_rule) in ( @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) + #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 diff --git a/test/test_utils.jl b/test/test_utils.jl index 399204d6b1..9bff00058e 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -296,7 +296,7 @@ module DummyRefShapes end ############################################################ -# Inverse parametric mapping ξ = ϕ(x) for testing hessians # +# 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) From 59deb76b72e12b225539dfbe8101bac08b913492 Mon Sep 17 00:00:00 2001 From: lijas Date: Mon, 27 May 2024 16:17:50 +0200 Subject: [PATCH 08/16] fix all linting errors (using pre-commit) --- docs/src/topics/FEValues.md | 4 ++-- src/FEValues/CellValues.jl | 2 +- src/FEValues/FunctionValues.jl | 2 +- src/interpolations.jl | 4 ++-- test/test_cellvalues.jl | 2 +- test/test_facevalues.jl | 4 ++-- test/test_interpolations.jl | 2 +- test/test_utils.jl | 2 +- 8 files changed, 11 insertions(+), 11 deletions(-) diff --git a/docs/src/topics/FEValues.md b/docs/src/topics/FEValues.md index 47764c3d67..f4abfdd3cd 100644 --- a/docs/src/topics/FEValues.md +++ b/docs/src/topics/FEValues.md @@ -77,12 +77,12 @@ Second order gradients of the shape functions are computed as ```math \begin{align*} - \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}x_j} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\left[\frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}\xi_s}\right]J^{-1}_{sj} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\left[- J^{-1}_{rk}\mathcal{H}_{kps} J^{-1}_{pi}\right] J^{-1}_{sj} = - \frac{\mathrm{d} \hat N}{\mathrm{d} x_k}\mathcal{H}_{kps} J^{-1}_{pi}J^{-1}_{sj} + \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}x_j} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\left[\frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}\xi_s}\right]J^{-1}_{sj} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\left[- J^{-1}_{rk}\mathcal{H}_{kps} J^{-1}_{pi}\right] J^{-1}_{sj} = - \frac{\mathrm{d} \hat N}{\mathrm{d} x_k}\mathcal{H}_{kps} J^{-1}_{pi}J^{-1}_{sj} \end{align*} ``` where we have used that the inverse of the jacobian can be computed as: - + ```math \begin{align*} 0 = \frac{\mathrm{d}}{\mathrm{d}\xi_s} (J_{kr} J^{-1}_{ri} ) = \frac{\mathrm{d}J_{kp}}{\mathrm{d}\xi_s} J^{-1}_{pi} + J_{kr} \frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}\xi_s} = 0 \quad \Rightarrow \\ diff --git a/src/FEValues/CellValues.jl b/src/FEValues/CellValues.jl index e29e5e42b1..e839489a83 100644 --- a/src/FEValues/CellValues.jl +++ b/src/FEValues/CellValues.jl @@ -46,7 +46,7 @@ struct CellValues{FV, GM, QR, detT} <: AbstractCellValues qr::QR # QuadratureRule detJdV::detT # AbstractVector{<:Number} or Nothing end -function CellValues(::Type{T}, qr::QuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation; +function CellValues(::Type{T}, qr::QuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation; update_gradients ::Union{Bool,Nothing} = nothing, #Use Union{Bool,Nothing} to get type-stable code update_hessians ::Union{Bool,Nothing} = nothing, update_detJdV ::Union{Bool,Nothing} = nothing) where T diff --git a/src/FEValues/FunctionValues.jl b/src/FEValues/FunctionValues.jl index 86f9b8b0dd..9745bc5d52 100644 --- a/src/FEValues/FunctionValues.jl +++ b/src/FEValues/FunctionValues.jl @@ -68,7 +68,7 @@ function FunctionValues{DiffOrder}(::Type{T}, ip::Interpolation, qr::QuadratureR Nξ = zeros(typeof_N(T, ip, ip_geo), n_shape, n_qpoints) Nx = isa(mapping_type(ip), IdentityMapping) ? Nξ : similar(Nξ) 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) diff --git a/src/interpolations.jl b/src/interpolations.jl index 4545ac601e..5496b7a457 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -1557,12 +1557,12 @@ function _shape_hessian_gradient_and_value_static_array(ipv::VectorizedInterpola 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) diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index 2b29b44f73..5289517b69 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -105,7 +105,7 @@ _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) diff --git a/test/test_facevalues.jl b/test/test_facevalues.jl index 6321d5ff0a..cfb9d31cf1 100644 --- a/test/test_facevalues.jl +++ b/test/test_facevalues.jl @@ -45,7 +45,7 @@ for (scalar_interpol, quad_rule) in ( (rand(Tensor{1, rdim}), rand(Tensor{2, rdim}), Tensor{3, rdim}((i,j,k)-> i==j==k ? rand() : 0.0)) end - u_funk(x,V,G,H) = begin + u_funk(x,V,G,H) = begin if update_hessians 0.5*x⋅H⋅x + G⋅x + V else @@ -105,7 +105,7 @@ for (scalar_interpol, quad_rule) in ( _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) diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index 1c1d1bec09..fc8ec92f5b 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -232,7 +232,7 @@ end @test G ≈ G_sa @test H ≈ H_sa end - + ips = Lagrange{RefQuadrilateral,2}() vdim = 3 ipv = ips^vdim diff --git a/test/test_utils.jl b/test/test_utils.jl index 9bff00058e..1b71a03665 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -331,4 +331,4 @@ function MAPPING(interpolation, cell_coordinates, global_coordinate::Vec{dim,T}) ξ -= inv(J) ⋅ residual end return ξ -end \ No newline at end of file +end From ed578c51f53a3a1e9926ee1e981b067176b059fc Mon Sep 17 00:00:00 2001 From: lijas Date: Tue, 28 May 2024 16:26:51 +0200 Subject: [PATCH 09/16] add shape_hessian_type for cellvalues --- src/FEValues/CellValues.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/FEValues/CellValues.jl b/src/FEValues/CellValues.jl index e839489a83..d1f2187739 100644 --- a/src/FEValues/CellValues.jl +++ b/src/FEValues/CellValues.jl @@ -88,6 +88,7 @@ function_interpolation(cv::CellValues) = function_interpolation(cv.fun_values) function_difforder(cv::CellValues) = function_difforder(cv.fun_values) shape_value_type(cv::CellValues) = shape_value_type(cv.fun_values) shape_gradient_type(cv::CellValues) = shape_gradient_type(cv.fun_values) +shape_hessian_type(cv::CellValues) = shape_hessian_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) From 9f205f1bdaa9bccae4d27ea6146159a6f5751dec Mon Sep 17 00:00:00 2001 From: lijas Date: Tue, 28 May 2024 17:08:20 +0200 Subject: [PATCH 10/16] test constructor for embedded cellvalues with update_hessians=true --- test/test_cellvalues.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index 5289517b69..8eba22906b 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -352,6 +352,18 @@ end @test zeros(vdim) == function_gradient(csv3, 1, ue)[:, 3] end end + + @testset "CellValues with hessians" begin + ip = Lagrange{RefQuadrilateral,2}() + qr = QuadratureRule{RefQuadrilateral}(2) + + 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) + end end @testset "CellValues constructor entry points" begin From 49288485cd1f9eb16ecef3869823c096942512ab Mon Sep 17 00:00:00 2001 From: lijas Date: Tue, 28 May 2024 17:17:11 +0200 Subject: [PATCH 11/16] use Ferrite.find_local_coordiante in test --- src/PointEvalHandler.jl | 11 +++++------ test/test_utils.jl | 26 ++------------------------ 2 files changed, 7 insertions(+), 30 deletions(-) diff --git a/src/PointEvalHandler.jl b/src/PointEvalHandler.jl index 304a425843..7252eccd71 100644 --- a/src/PointEvalHandler.jl +++ b/src/PointEvalHandler.jl @@ -114,17 +114,16 @@ function _check_isoparametric_boundaries(::Type{RefSimplex{dim}}, x_local::Vec{d end # See https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975/2 -# TODO: should we make iteration params optional keyword arguments? -function find_local_coordinate(interpolation, cell_coordinates::Vector{V}, global_coordinate::V) where {dim, T, V <: Vec{dim, T}} +# TODO: should we make iteration params optional keyword arguments? +function find_local_coordinate(interpolation, cell_coordinates::Vector{V}, global_coordinate::V2; tol_norm = 1e-10) where {dim, T, T2, V <: Vec{dim, T}, V2 <: Vec{dim, T2}} n_basefuncs = getnbasefunctions(interpolation) @assert length(cell_coordinates) == n_basefuncs - local_guess = zero(V) + local_guess = zero(V2) max_iters = 10 - tol_norm = 1e-10 converged = false for _ in 1:max_iters - global_guess = zero(V) - J = zero(Tensor{2, dim, T}) + global_guess = zero(V2) + J = zero(Tensor{2, dim, T2}) # TODO batched eval after 764 is merged. for j in 1:n_basefuncs dNdξ, N = shape_gradient_and_value(interpolation, local_guess, j) diff --git a/test/test_utils.jl b/test/test_utils.jl index 1b71a03665..4e7d814318 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -302,33 +302,11 @@ function function_value_from_physical_coord(interpolation::Interpolation, cell_c n_basefuncs = getnbasefunctions(interpolation) scalar_ip = interpolation isa Ferrite.ScalarInterpolation ? interpolation : interpolation.ip @assert length(ue) == n_basefuncs - ξ = MAPPING(scalar_ip, cell_coordinates, X) + _, ξ = Ferrite.find_local_coordinate(scalar_ip, cell_coordinates, X; tol_norm=1e-16) u = zero(shape_value(interpolation, ξ, 1)) 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 +end \ No newline at end of file From 5e1ec33f7aad71353f81da08f37b99010c023ae1 Mon Sep 17 00:00:00 2001 From: lijas Date: Tue, 28 May 2024 17:17:41 +0200 Subject: [PATCH 12/16] precommit --- src/PointEvalHandler.jl | 2 +- test/test_cellvalues.jl | 2 +- test/test_utils.jl | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/PointEvalHandler.jl b/src/PointEvalHandler.jl index 7252eccd71..fd534d3cea 100644 --- a/src/PointEvalHandler.jl +++ b/src/PointEvalHandler.jl @@ -114,7 +114,7 @@ function _check_isoparametric_boundaries(::Type{RefSimplex{dim}}, x_local::Vec{d end # See https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975/2 -# TODO: should we make iteration params optional keyword arguments? +# TODO: should we make iteration params optional keyword arguments? function find_local_coordinate(interpolation, cell_coordinates::Vector{V}, global_coordinate::V2; tol_norm = 1e-10) where {dim, T, T2, V <: Vec{dim, T}, V2 <: Vec{dim, T2}} n_basefuncs = getnbasefunctions(interpolation) @assert length(cell_coordinates) == n_basefuncs diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index 8eba22906b..1d7cc8d79d 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -356,7 +356,7 @@ end @testset "CellValues with hessians" begin ip = Lagrange{RefQuadrilateral,2}() qr = QuadratureRule{RefQuadrilateral}(2) - + cv_vector = CellValues(qr, ip^2, ip^3; update_hessians = true) cv_scalar = CellValues(qr, ip, ip^3; update_hessians = true) diff --git a/test/test_utils.jl b/test/test_utils.jl index 4e7d814318..8bdabf9d20 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -309,4 +309,4 @@ function function_value_from_physical_coord(interpolation::Interpolation, cell_c u += N * ue[j] end return u -end \ No newline at end of file +end From 59815143a08e3897d8846c30b6e8658d50d13200 Mon Sep 17 00:00:00 2001 From: lijas Date: Wed, 29 May 2024 10:35:11 +0200 Subject: [PATCH 13/16] Apply suggestions from code review Co-authored-by: Knut Andreas Meyer --- src/FEValues/CellValues.jl | 2 +- src/FEValues/FacetValues.jl | 4 ++-- src/FEValues/FunctionValues.jl | 4 ++-- src/PointEvalHandler.jl | 9 +++++---- 4 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/FEValues/CellValues.jl b/src/FEValues/CellValues.jl index d1f2187739..00d8eb298b 100644 --- a/src/FEValues/CellValues.jl +++ b/src/FEValues/CellValues.jl @@ -12,7 +12,7 @@ values of nodal functions, gradients and divergences of nodal functions etc. in By default linear Lagrange interpolation is used. For embedded elements the geometric interpolations should be vectorized to the spatial dimension. -**Keyword arguments:** +**Keyword arguments:** The following keyword arguments are experimental and may change in future minor releases * `update_gradients`: Specifies if the gradients of the shape functions should be updated (default true) * `update_hessians`: Specifies if the hessians of the shape functions should be updated (default false) * `update_detJdV`: Specifies if the volume associated with each quadrature point should be updated (default true) diff --git a/src/FEValues/FacetValues.jl b/src/FEValues/FacetValues.jl index d712c0de5b..e4268c1edd 100644 --- a/src/FEValues/FacetValues.jl +++ b/src/FEValues/FacetValues.jl @@ -1,5 +1,5 @@ """ - FacetValues([::Type{T}], quad_rule::FacetQuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation]; ; update_gradients=true, update_hessians=false, ) + FacetValues([::Type{T}], quad_rule::FacetQuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation]) A `FacetValues` object facilitates the process of evaluating values of shape functions, gradients of shape functions, values of nodal functions, gradients and divergences of nodal functions etc. on the facets of finite elements. @@ -12,7 +12,7 @@ values of nodal functions, gradients and divergences of nodal functions etc. on * `geom_interpol`: an optional instance of an [`Interpolation`](@ref) which is used to interpolate the geometry. By default linear Lagrange interpolation is used. -**Keyword arguments:** +**Keyword arguments:** The following keyword arguments are experimental and may change in future minor releases * `update_gradients`: Specifies if the gradients of the shape functions should be updated (default true) * `update_hessians`: Specifies if the hessians of the shape functions should be updated (default false) diff --git a/src/FEValues/FunctionValues.jl b/src/FEValues/FunctionValues.jl index 9745bc5d52..a4c2981e17 100644 --- a/src/FEValues/FunctionValues.jl +++ b/src/FEValues/FunctionValues.jl @@ -203,12 +203,12 @@ end sdim, rdim = size(Jinv) (rdim != sdim) && error("apply_mapping! for second order gradients and embedded elements not implemented") - H = gethessian(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 ? + 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 diff --git a/src/PointEvalHandler.jl b/src/PointEvalHandler.jl index 3fe09e4481..42c9aac9b2 100644 --- a/src/PointEvalHandler.jl +++ b/src/PointEvalHandler.jl @@ -115,15 +115,16 @@ end # See https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975/2 # TODO: should we make iteration params optional keyword arguments? -function find_local_coordinate(interpolation, cell_coordinates::Vector{V}, global_coordinate::V2; tol_norm = 1e-10) where {dim, T, T2, V <: Vec{dim, T}, V2 <: Vec{dim, T2}} +function find_local_coordinate(interpolation, cell_coordinates::Vector{<:Vec{dim}}, global_coordinate::Vec{dim}; tol_norm = 1e-10) where dim + T = promote_type(eltype(cell_coordinates[1]), eltype(global_coordinate)) n_basefuncs = getnbasefunctions(interpolation) @assert length(cell_coordinates) == n_basefuncs - local_guess = zero(V2) + local_guess = zero(Vec{dim, T}) max_iters = 10 converged = false for _ in 1:max_iters - global_guess = zero(V2) - J = zero(Tensor{2, dim, T2}) + global_guess = zero(Vec{dim, T}) + J = zero(Tensor{2, dim, T}) # TODO batched eval after 764 is merged. for j in 1:n_basefuncs dNdξ, N = shape_gradient_and_value(interpolation, local_guess, j) From 89d059029ee0f9a383407dee2de22ecfecc74027 Mon Sep 17 00:00:00 2001 From: lijas Date: Wed, 29 May 2024 10:42:53 +0200 Subject: [PATCH 14/16] add excplicit imports and run pre-commit --- CHANGELOG.md | 2 +- src/FEValues/CellValues.jl | 2 +- src/Ferrite.jl | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d79453f347..6f3e1d2111 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -320,7 +320,7 @@ more discussion). - A new function, `geometric_interpolation`, is exported, which gives the geometric interpolation for each cell type. This is equivalent to the deprecated `Ferrite.default_interpolation` function. ([#953][github-953]) - + - CellValues and FacetValues can now store and map second order gradients (Hessians). The number of gradients computed in CellValues/FacetValues is specified using the keyword arguments `update_gradients::Bool` (default true) and `update_hessians::Bool` (default false) in the diff --git a/src/FEValues/CellValues.jl b/src/FEValues/CellValues.jl index 00d8eb298b..66353e59e5 100644 --- a/src/FEValues/CellValues.jl +++ b/src/FEValues/CellValues.jl @@ -1,5 +1,5 @@ """ - CellValues([::Type{T},] quad_rule::QuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation]; update_detJdV=true, update_gradients=true, update_hessians=false) + CellValues([::Type{T},] quad_rule::QuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation]) A `CellValues` object facilitates the process of evaluating values of shape functions, gradients of shape functions, values of nodal functions, gradients and divergences of nodal functions etc. in the finite element cell. diff --git a/src/Ferrite.jl b/src/Ferrite.jl index a766c11fc9..5e794c9c93 100644 --- a/src/Ferrite.jl +++ b/src/Ferrite.jl @@ -22,7 +22,7 @@ using WriteVTK: WriteVTK, VTKCellTypes using Tensors: Tensors, AbstractTensor, SecondOrderTensor, SymmetricTensor, Tensor, Vec, gradient, - rotation_tensor, symmetric, tovoigt! + rotation_tensor, symmetric, tovoigt!, hessian, otimesu using ForwardDiff: ForwardDiff From df0d06d7f60aae190742c01be1b0c87b9650548c Mon Sep 17 00:00:00 2001 From: lijas Date: Wed, 29 May 2024 11:08:14 +0200 Subject: [PATCH 15/16] never ending ci fixes --- src/interpolations.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/interpolations.jl b/src/interpolations.jl index e7c65d52ec..016f2dde08 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -1540,23 +1540,23 @@ 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)) + ξd = Tensors._load(Tensors._load(ξ, ForwardDiff.Tag(f, V)), ForwardDiff.Tag(f, V)) value_hess = f(ξd) # Extract the value and gradient - val = Vec{vdim, T}(i -> Tensors.value(Tensors.value(value_hess[i]))) + val = Vec{vdim, T}(i -> ForwardDiff.value(ForwardDiff.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 = ForwardDiff.value(vi) - hess_values_partials = Tensors.partials(hess_values) + hess_values_partials = ForwardDiff.partials(hess_values) for (k, pk) in pairs(hess_values_partials) grad[i, k] = pk end - hess_partials = Tensors.partials(vi) + hess_partials = ForwardDiff.partials(vi) for (j, partial_j) in pairs(hess_partials) - hess_partials_partials = Tensors.partials(partial_j) + hess_partials_partials = ForwardDiff.partials(partial_j) for (k, pk) in pairs(hess_partials_partials) hess[i, j, k] = pk end From 7de217655abde63ad0ad46e6234a54b2af476324 Mon Sep 17 00:00:00 2001 From: lijas Date: Wed, 29 May 2024 15:05:30 +0200 Subject: [PATCH 16/16] re-add test for deprecated function --- test/test_cellvalues.jl | 2 +- test/test_facevalues.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index 1d7cc8d79d..05ed1444ff 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -93,7 +93,7 @@ if rdim == 3 @test (@test_deprecated function_curl(cv, i, ue_vec)) ≈ Ferrite.curl_from_gradient(G_vector) end - @test_deprecated function_value(cv, i, ue_vec) + @test_deprecated function_value(cv, i, ue_vec) #no value to test against end end diff --git a/test/test_facevalues.jl b/test/test_facevalues.jl index cfb9d31cf1..99568971bc 100644 --- a/test/test_facevalues.jl +++ b/test/test_facevalues.jl @@ -93,7 +93,7 @@ for (scalar_interpol, quad_rule) in ( if rdim == 3 @test (@test_deprecated function_curl(fv, i, ue_vec)) ≈ Ferrite.curl_from_gradient(G_vector) end - function_value(fv, i, ue_vec) + @test_deprecated function_value(fv, i, ue_vec) #no value to test against end end