Skip to content

Commit

Permalink
add higher order gradients in fe values
Browse files Browse the repository at this point in the history
  • Loading branch information
@lijas
lijas committed May 22, 2024
1 parent c015df0 commit a429acd
Show file tree
Hide file tree
Showing 12 changed files with 444 additions and 121 deletions.
43 changes: 43 additions & 0 deletions docs/src/topics/FEValues.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,49 @@ For scalar fields, we always use scalar base functions. For tensorial fields (no
\end{align*}
```

Second order gradients of the shape functions are computed as

```math
\begin{align*}
\mathrm{grad}(\mathrm{grad}(N(\boldsymbol{x}))) = \frac{\mathrm{d}^2 N}{\mathrm{d}\boldsymbol{x}^2} = \boldsymbol{J}^{-T} \cdot \frac{\mathrm{d}^2\hat{N}}{\mathrm{d}\boldsymbol{\xi}^2} \cdot \boldsymbol{J}^{-1} - \boldsymbol{J}^{-T} \cdot\mathrm{grad}(N) \cdot \boldsymbol{\mathcal{H}} \cdot \boldsymbol{J}^{-1}
\end{align*}
```
!!! details "Derivation"
The gradient of the shape functions is obtained using the chain rule:

\begin{align*}
\frac{\mathrm{d} N}{\mathrm{d}x_i} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\frac{\mathrm{d} \xi_r}{\mathrm{d} x_i} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r} J^{-1}_{ri}
\end{align*}

For the second order gradients, we first use the product rule on the equation above:

\begin{align*}
\frac{\mathrm{d}^2 N}{\mathrm{d}x_i \mathrm{d}x_j} = \frac{\mathrm{d}}{\mathrm{d}x_j}(\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}) J^{-1}_{ri} + \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r} \frac{\mathrm{d}}{\mathrm{d}x_j}(J^{-1}_{ri})
\end{align*}

The first term can be computed as:

\begin{align*}
\frac{\mathrm{d}}{\mathrm{d}x_j}(\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}) J^{-1}_{ri} = J^{-1}_{si}\frac{\mathrm{d}}{\mathrm{d}\xi_s}(\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}) J^{-1}_{ri} = J^{-1}_{si}\frac{\mathrm{d}^2 \hat N}{\mathrm{d} \xi_s\mathrm{d} \xi_r} J^{-1}_{ri}
\end{align*}

The second term can be written as:

\begin{align*}
\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\frac{\mathrm{d}}{\mathrm{d}x_j}(J^{-1}_{ri}) = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}[\frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}\xi_s}]\frac{\mathrm{d} \xi_s}{\mathrm{d}x_j} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}[ - J^{-1}_{ri}\mathcal{H}_{ris} J^{-1}_{li}]J^{-1}_{sj} = - \frac{\mathrm{d} \hat N}{\mathrm{d} x_i}\mathcal{H}_{ris} J^{-1}_{li}J^{-1}_{sj}
\end{align*}

where we have used that the inverse of the jacobian can be computed as:

\begin{align*}
0 = \frac{\mathrm{d}}{\mathrm{d}\xi_s} (J_{rl} J^{-1}_{li} ) = \frac{\mathrm{d}J_{rl}}{\mathrm{d}\xi_s} J^{-1}_{li} + J_{ri} \frac{\mathrm{d}J^{-1}_{li}}{\mathrm{d}\xi_s} = 0 \quad \Rightarrow \\
\end{align*}

\begin{align*}
\frac{\mathrm{d}J^{-1}_{li}}{\mathrm{d}\xi_s} = - J^{-1}_{ri}\frac{\mathrm{d}J_{rl}}{\mathrm{d}\xi_s} J^{-1}_{li} = - J^{-1}_{ri}\mathcal{H}_{ris} J^{-1}_{li}\\
\end{align*}


### Covariant Piola mapping, H(curl)
`Ferrite.CovariantPiolaMapping`

Expand Down
18 changes: 12 additions & 6 deletions src/FEValues/CellValues.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,15 +41,20 @@ 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

_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
function CellValues(::Type{T}, qr::QuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation;
update_gradients::Bool = true, update_hessians = false, update_detJdV::Bool = true) 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

FunDiffOrder = convert(Int, _update_gradients)
FunDiffOrder += convert(Int, _update_hessians)
GeoDiffOrder = max(required_geo_diff_order(mapping_type(ip_fun), FunDiffOrder), _update_detJdV)
geo_mapping = GeometryMapping{GeoDiffOrder}(T, ip_geo.ip, qr)
fun_values = FunctionValues{FunDiffOrder}(T, ip_fun, qr, ip_geo)
detJdV = _update_detJdV ? fill(T(NaN), length(getweights(qr))) : nothing
detJdV = update_detJdV ? fill(T(NaN), length(getweights(qr))) : nothing
return CellValues(fun_values, geo_mapping, qr, detJdV)
end

Expand Down Expand Up @@ -79,6 +84,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
Expand Down
19 changes: 13 additions & 6 deletions src/FEValues/FacetValues.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,10 +40,16 @@ 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 = convert(Int, _update_gradients)
FunDiffOrder += convert(Int, _update_hessians)
GeoDiffOrder = max(required_geo_diff_order(mapping_type(ip_fun), FunDiffOrder), 1)
# Not sure why the type-asserts are required here but not for CellValues,
# but they solve the type-instability for FacetValues
Expand All @@ -55,9 +61,9 @@ function FacetValues(::Type{T}, fqr::FacetQuadratureRule, ip_fun::Interpolation,
return FacetValues(fun_values, geo_mapping, fqr, detJdV, normals, ScalarWrapper(1))
end

FacetValues(qr::FacetQuadratureRule, ip::Interpolation, args...; kwargs...) = FacetValues(Float64, qr, ip, args...; kwargs...)
function FacetValues(::Type{T}, qr::FacetQuadratureRule, ip::Interpolation, ip_geo::ScalarInterpolation; kwargs...) where T
return FacetValues(T, qr, ip, VectorizedInterpolation(ip_geo); kwargs...)
FaceValues(qr::FacetQuadratureRule, ip::Interpolation, args...; kwargs...) = FaceValues(Float64, qr, ip, args...; kwargs...)
function FaceValues(::Type{T}, qr::FacetQuadratureRule, ip::Interpolation, ip_geo::ScalarInterpolation; kwargs...) where T
return FaceValues(T, qr, ip, VectorizedInterpolation(ip_geo); kwargs...)
end

function Base.copy(fv::FacetValues)
Expand All @@ -84,6 +90,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)

"""
Expand Down
95 changes: 74 additions & 21 deletions src/FEValues/FunctionValues.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -33,36 +43,51 @@ for both the reference cell (precalculated) and the real cell (updated in `reini
"""
FunctionValues

struct FunctionValues{DiffOrder, IP, N_t, dNdx_t, dNdξ_t}
struct FunctionValues{DiffOrder, IP, N_t, dNdx_t, dNdξ_t, d2Ndx2_t, d2Ndξ2_t}
ip::IP # ::Interpolation
Nx::N_t # ::AbstractMatrix{Union{<:Tensor,<:Number}}
Nξ::N_t # ::AbstractMatrix{Union{<:Tensor,<:Number}}
dNdx::dNdx_t # ::AbstractMatrix{Union{<:Tensor,<:StaticArray}} or Nothing
dNdξ::dNdξ_t # ::AbstractMatrix{Union{<:Tensor,<:StaticArray}} or Nothing
function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, ::Nothing, ::Nothing) where {N_t<:AbstractMatrix}
return new{0, typeof(ip), N_t, Nothing, Nothing}(ip, Nx, Nξ, nothing, nothing)
d2Ndx2::d2Ndx2_t # ::AbstractMatrix{<:Tensor{2}} Hessians of geometric shape functions in ref-domain
d2Ndξ2::d2Ndξ2_t # ::AbstractMatrix{<:Tensor{2}} Hessians of geometric shape functions in ref-domain
function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, ::Nothing, ::Nothing, ::Nothing, ::Nothing) where {N_t<:AbstractMatrix}
return new{0, typeof(ip), N_t, Nothing, Nothing, Nothing, Nothing}(ip, Nx, Nξ, nothing, nothing, nothing, nothing)
end
function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, dNdx::AbstractMatrix, dNdξ::AbstractMatrix) 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}
sdim = n_components(ip_geo)
rdim = getdim(ip)
((rdim != sdim) && DiffOrder>1) && error("CellValues/FunctionValues for embedded elements and higher order gradient not implemented")

n_shape = getnbasefunctions(ip)
n_qpoints = getnquadpoints(qr)

Nξ = zeros(typeof_N(T, ip, ip_geo), n_shape, n_qpoints)
Nx = isa(mapping_type(ip), IdentityMapping) ? Nξ : similar(Nξ)

if DiffOrder == 0
dNdξ = dNdx = nothing
elseif DiffOrder == 1
dNdξ = dNdx = d2Ndξ2 = d2Ndx2 = nothing

if DiffOrder >= 1
dNdξ = zeros(typeof_dNdξ(T, ip, ip_geo), n_shape, n_qpoints)
dNdx = fill(zero(typeof_dNdx(T, ip, ip_geo)) * T(NaN), n_shape, n_qpoints)
else
end

if DiffOrder >= 2
d2Ndξ2 = zeros(typeof_d2Ndξ2(T, ip, ip_geo), n_shape, n_qpoints)
d2Ndx2 = fill(zero(typeof_d2Ndx2(T, ip, ip_geo)) * T(NaN), n_shape, n_qpoints)
end

if DiffOrder > 2
throw(ArgumentError("Currently only values and gradients can be updated in FunctionValues"))
end

fv = FunctionValues(ip, Nx, Nξ, dNdx, dNdξ)
fv = FunctionValues(ip, Nx, Nξ, dNdx, dNdξ, d2Ndx2, d2Ndξ2)
precompute_values!(fv, getpoints(qr)) # Separate function for qr point update in PointValues
return fv
end
Expand All @@ -73,25 +98,34 @@ 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
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
Expand Down Expand Up @@ -167,6 +201,25 @@ 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))
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
Loading

0 comments on commit a429acd

Please sign in to comment.