diff --git a/CHANGELOG.md b/CHANGELOG.md index 85fd9cc901..6f3e1d2111 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -321,6 +321,12 @@ more discussion). 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 + constructors, i.e. `CellValues(...; update_hessians=true)`. ([#953][github-938]) + + ### 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 df54430cbc..f4abfdd3cd 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}\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}\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}\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*} + ``` + + The second term can be written 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} + \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..66353e59e5 100644 --- a/src/FEValues/CellValues.jl +++ b/src/FEValues/CellValues.jl @@ -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:** 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) + **Common methods:** * [`reinit!`](@ref) @@ -42,11 +47,18 @@ struct CellValues{FV, GM, QR, detT} <: AbstractCellValues detJdV::detT # AbstractVector{<:Number} or Nothing end function CellValues(::Type{T}, qr::QuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation; - update_gradients::Union{Bool,Nothing} = nothing, update_detJdV::Union{Bool,Nothing} = nothing) where T + update_gradients ::Union{Bool,Nothing} = nothing, #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 + _update_detJdV = update_detJdV === nothing ? true : update_detJdV + _update_hessians && @assert _update_gradients - _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 + 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 @@ -76,9 +88,11 @@ 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) +@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..e4268c1edd 100644 --- a/src/FEValues/FacetValues.jl +++ b/src/FEValues/FacetValues.jl @@ -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:** 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) + **Common methods:** * [`reinit!`](@ref) @@ -41,12 +46,16 @@ struct FacetValues{FV, GM, FQR, detT, nT, V_FV<:AbstractVector{FV}, V_GM<:Abstra 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} + 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 +93,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..a4c2981e17 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) 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, ::Nothing, ::Nothing) where {N_t<:AbstractMatrix} + return new{1, typeof(ip), N_t, typeof(dNdx), typeof(dNdξ), Nothing, Nothing}(ip, Nx, Nξ, dNdx, dNdξ, nothing, nothing) + end + function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, dNdx::AbstractMatrix, dNdξ::AbstractMatrix, d2Ndx2::AbstractMatrix, d2Ndξ2::AbstractMatrix) where {N_t<:AbstractMatrix} + return new{2, typeof(ip), N_t, typeof(dNdx), typeof(dNdξ), typeof(d2Ndx2), typeof(d2Ndξ2)}(ip, Nx, Nξ, dNdx, dNdξ, d2Ndx2, d2Ndξ2) 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ξ) + dNdξ = dNdx = d2Ndξ2 = d2Ndx2 = nothing - if DiffOrder == 0 - dNdξ = dNdx = nothing - elseif DiffOrder == 1 + 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 +@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 + 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..476d3a5f2f 100644 --- a/src/FEValues/GeometryMapping.jl +++ b/src/FEValues/GeometryMapping.jl @@ -14,7 +14,7 @@ struct MappingValues{JT, HT} end @inline getjacobian(mv::MappingValues{<:Union{AbstractTensor, SMatrix}}) = mv.J -# @inline gethessian(mv::MappingValues{<:Any,<:AbstractTensor}) = mv.H # PR798 +@inline gethessian(mv::MappingValues{<:Any,<:AbstractTensor}) = mv.H """ @@ -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 + # ::Nothing (dMdξ or d2Mdξ2 if not required) function GeometryMapping( ip::IP, M::M_t, ::Nothing, ::Nothing ) where {IP <: ScalarInterpolation, M_t<:AbstractMatrix{<:Number}} @@ -46,12 +46,12 @@ struct GeometryMapping{DiffOrder, IP, M_t, dMdξ_t, d2Mdξ2_t} ) where {IP <: ScalarInterpolation, M_t<:AbstractMatrix{<:Number}, dMdξ_t <: AbstractMatrix{<:Vec}} return new{1, IP, M_t, dMdξ_t, Nothing}(ip, M, dMdξ, nothing) end -#= function GeometryMapping( + function GeometryMapping( ip::IP, M::M_t, dMdξ::dMdξ_t, d2Mdξ2::d2Mdξ2_t) where {IP <: ScalarInterpolation, M_t<:AbstractMatrix{<:Number}, dMdξ_t <: AbstractMatrix{<:Vec}, d2Mdξ2_t <: AbstractMatrix{<:Tensor{2}}} return new{2, IP, M_t, dMdξ_t, d2Mdξ2_t}(ip, M, dMdξ, d2Mdξ2) - end =# # PR798 + end end function GeometryMapping{0}(::Type{T}, ip::ScalarInterpolation, qr::QuadratureRule) where T n_shape = getnbasefunctions(ip) @@ -71,7 +71,7 @@ function GeometryMapping{1}(::Type{T}, ip::ScalarInterpolation, qr::QuadratureRu precompute_values!(gm, getpoints(qr)) return gm end -#= function GeometryMapping{2}(::Type{T}, ip::ScalarInterpolation, qr::QuadratureRule) where T +function GeometryMapping{2}(::Type{T}, ip::ScalarInterpolation, qr::QuadratureRule) where T n_shape = getnbasefunctions(ip) n_qpoints = getnquadpoints(qr) @@ -82,7 +82,7 @@ end gm = GeometryMapping(ip, M, dMdξ, d2Mdξ2) precompute_values!(gm, getpoints(qr)) return gm -end =# # PR798 +end function precompute_values!(gm::GeometryMapping{0}, qr_points::Vector{<:Vec}) shape_values!(gm.M, gm.ip, qr_points) @@ -90,9 +90,9 @@ end function precompute_values!(gm::GeometryMapping{1}, qr_points::Vector{<:Vec}) shape_gradients_and_values!(gm.dMdξ, gm.M, gm.ip, qr_points) end -#= function precompute_values!(gm::GeometryMapping{2}, qr_points::Vector{<:Vec}) +function precompute_values!(gm::GeometryMapping{2}, qr_points::Vector{<:Vec}) shape_hessians_gradients_and_values!(gm.d2Mdξ2, gm.dMdξ, gm.M, gm.ip, qr_points) -end =# # PR798 +end function Base.copy(v::GeometryMapping) return GeometryMapping(copy(v.ip), copy(v.M), _copy_or_nothing(v.dMdξ), _copy_or_nothing(v.d2Mdξ2)) @@ -117,9 +117,9 @@ end function otimes_returntype(#=typeof(x)=#::Type{<:Vec{dim,Tx}}, #=typeof(dMdξ)=#::Type{<:Vec{dim,TM}}) where {dim, Tx, TM} return Tensor{2,dim,promote_type(Tx,TM)} end -#= function otimes_returntype(#=typeof(x)=#::Type{<:Vec{dim,Tx}}, #=typeof(d2Mdξ2)=#::Type{<:Tensor{2,dim,TM}}) where {dim, Tx, TM} +function otimes_returntype(#=typeof(x)=#::Type{<:Vec{dim,Tx}}, #=typeof(d2Mdξ2)=#::Type{<:Tensor{2,dim,TM}}) where {dim, Tx, TM} return Tensor{3,dim,promote_type(Tx,TM)} -end =# # PR798 +end @inline function calculate_mapping(::GeometryMapping{0}, q_point, x) return MappingValues(nothing, nothing) @@ -134,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..88a117cc68 100644 --- a/src/FEValues/common_values.jl +++ b/src/FEValues/common_values.jl @@ -211,6 +211,38 @@ 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_symmetric_gradient(fe_v::AbstractValues, q_point::Int, u::AbstractVector, [dof_range]) @@ -308,10 +340,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 7646130fb7..5e794c9c93 100644 --- a/src/Ferrite.jl +++ b/src/Ferrite.jl @@ -17,12 +17,12 @@ using OrderedCollections: using SparseArrays: SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals, sparse using StaticArrays: - StaticArrays, SMatrix, SVector + StaticArrays, MArray, MMatrix, SArray, SMatrix, SVector 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 diff --git a/src/PointEvalHandler.jl b/src/PointEvalHandler.jl index 56f417f292..42c9aac9b2 100644 --- a/src/PointEvalHandler.jl +++ b/src/PointEvalHandler.jl @@ -115,15 +115,15 @@ 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}} +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(V) + local_guess = zero(Vec{dim, T}) max_iters = 10 - tol_norm = 1e-10 converged = false for _ in 1:max_iters - global_guess = zero(V) + 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 diff --git a/src/exports.jl b/src/exports.jl index 681613cff7..a8fe8701f2 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 60f58a4f83..016f2dde08 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) @@ -1531,6 +1529,42 @@ function shape_gradient_and_value(ipv::VectorizedInterpolation{vdim, shape}, ξ: return 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(ξ, ForwardDiff.Tag(f, V)), ForwardDiff.Tag(f, V)) + value_hess = f(ξd) + # Extract the value and gradient + 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 = ForwardDiff.value(vi) + + hess_values_partials = ForwardDiff.partials(hess_values) + for (k, pk) in pairs(hess_values_partials) + grad[i, k] = pk + end + + hess_partials = ForwardDiff.partials(vi) + for (j, partial_j) in pairs(hess_partials) + hess_partials_partials = ForwardDiff.partials(partial_j) + for (k, pk) in pairs(hess_partials_partials) + hess[i, j, k] = pk + end + end + end + return SArray(hess), SMatrix(grad), val +end + reference_coordinates(ip::VectorizedInterpolation) = reference_coordinates(ip.ip) is_discontinuous(::Type{<:VectorizedInterpolation{<:Any, <:Any, <:Any, ip}}) where {ip} = is_discontinuous(ip) diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index 5c1c3cf0a0..05ed1444ff 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,106 @@ (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) + 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) @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) #no value to test against + 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 +123,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 +136,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 @@ -302,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 diff --git a/test/test_facevalues.jl b/test/test_facevalues.jl index 43cff3d366..99568971bc 100644 --- a/test/test_facevalues.jl +++ b/test/test_facevalues.jl @@ -13,70 +13,119 @@ 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) - fv = if VERSION ≥ v"1.9" - @inferred FacetValues(quad_rule, func_interpol, geom_interpol) - 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 = true + update_hessians = (DiffOrder==2 && Ferrite.getorder(func_interpol) > 1) + + 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) @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 - @test function_value(fv, i, u_vector) ≈ (@test_deprecated function_value(fv, i, u)) + @test_deprecated function_value(fv, i, ue_vec) #no value to test against 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 + 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..fc8ec92f5b 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -219,4 +219,32 @@ @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 + + 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 c0fb2e9e83..8bdabf9d20 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -294,3 +294,19 @@ 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 + _, ξ = 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