Higher order gradients in CellValues

docs/src/topics/FEValues.md

@@ -42,6 +42,49 @@ For scalar fields, we always use scalar base functions. For tensorial fields (no
 \end{align*}
 ```
 
+Second order gradients of the shape functions are computed as
+
+```math
+\begin{align*} 
+    \mathrm{grad}(\mathrm{grad}(N(\boldsymbol{x}))) = \frac{\mathrm{d}^2 N}{\mathrm{d}\boldsymbol{x}^2} = \boldsymbol{J}^{-T} \cdot \frac{\mathrm{d}^2\hat{N}}{\mathrm{d}\boldsymbol{\xi}^2} \cdot \boldsymbol{J}^{-1} -  \boldsymbol{J}^{-T} \cdot\mathrm{grad}(N) \cdot \boldsymbol{\mathcal{H}}  \cdot \boldsymbol{J}^{-1}
+\end{align*}
+```
+!!! details "Derivation"
+    The gradient of the shape functions is obtained using the chain rule:
+
+    \begin{align*} 
+        \frac{\mathrm{d} N}{\mathrm{d}x_i} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\frac{\mathrm{d} \xi_r}{\mathrm{d} x_i} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r} J^{-1}_{ri}
+    \end{align*}
+
+    For the second order gradients, we first use the product rule on the equation above:
+
+    \begin{align*} 
+        \frac{\mathrm{d}^2 N}{\mathrm{d}x_i \mathrm{d}x_j} = \frac{\mathrm{d}}{\mathrm{d}x_j}(\frac{\mathrm{d} \hat N}{\mathrm{d}   \xi_r}) J^{-1}_{ri} + \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r} \frac{\mathrm{d}}{\mathrm{d}x_j}(J^{-1}_{ri}) 
+    \end{align*}
+
+    The first term can be computed as:
+
+    \begin{align*} 
+        \frac{\mathrm{d}}{\mathrm{d}x_j}(\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}) J^{-1}_{ri} = J^{-1}_{si}\frac{\mathrm{d}}{\mathrm{d}\xi_s}(\frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}) J^{-1}_{ri} = J^{-1}_{si}\frac{\mathrm{d}^2 \hat N}{\mathrm{d} \xi_s\mathrm{d} \xi_r} J^{-1}_{ri}
+    \end{align*}
+
+    The second term can be written as:
+
+    \begin{align*} 
+        \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}\frac{\mathrm{d}}{\mathrm{d}x_j}(J^{-1}_{ri}) = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}[\frac{\mathrm{d}J^{-1}_{ri}}{\mathrm{d}\xi_s}]\frac{\mathrm{d} \xi_s}{\mathrm{d}x_j} = \frac{\mathrm{d} \hat N}{\mathrm{d} \xi_r}[ - J^{-1}_{ri}\mathcal{H}_{ris} J^{-1}_{li}]J^{-1}_{sj} = - \frac{\mathrm{d} \hat N}{\mathrm{d} x_i}\mathcal{H}_{ris} J^{-1}_{li}J^{-1}_{sj} 
+    \end{align*}
+
+    where we have used that the inverse of the jacobian can be computed as:
+
+    \begin{align*} 
+    0 = \frac{\mathrm{d}}{\mathrm{d}\xi_s} (J_{rl} J^{-1}_{li} ) = \frac{\mathrm{d}J_{rl}}{\mathrm{d}\xi_s} J^{-1}_{li}  + J_{ri} \frac{\mathrm{d}J^{-1}_{li}}{\mathrm{d}\xi_s} = 0 \quad \Rightarrow \\
+    \end{align*}
+
+    \begin{align*} 
+    \frac{\mathrm{d}J^{-1}_{li}}{\mathrm{d}\xi_s} = - J^{-1}_{ri}\frac{\mathrm{d}J_{rl}}{\mathrm{d}\xi_s} J^{-1}_{li} = - J^{-1}_{ri}\mathcal{H}_{ris} J^{-1}_{li}\\
+    \end{align*}
+
+
 ### Covariant Piola mapping, H(curl)
 `Ferrite.CovariantPiolaMapping`
 

src/FEValues/CellValues.jl

@@ -42,10 +42,17 @@ struct CellValues{FV, GM, QR, detT} <: AbstractCellValues
     detJdV::detT   # AbstractVector{<:Number} or Nothing
 end
 function CellValues(::Type{T}, qr::QuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation; 
-        update_gradients::Union{Bool,Nothing} = nothing, update_detJdV::Union{Bool,Nothing} = nothing) where T 
+        update_gradients::Union{Bool,Nothing} = nothing, 
+        update_hessians::Union{Bool,Nothing} = nothing, 
+        update_detJdV::Union{Bool,Nothing} = nothing) where T 
 
-    _update_detJdV = update_detJdV === nothing ? true : update_detJdV
-    FunDiffOrder = update_gradients === nothing ? 1 : convert(Int, update_gradients) # Logic must change when supporting update_hessian kwargs
+    _update_gradients = update_gradients === nothing ? true : update_gradients
+    _update_hessians  = update_hessians  === nothing ? false : update_hessians
+    _update_detJdV    = update_detJdV    === nothing ? true : update_detJdV
+    _update_hessians && @assert _update_gradients
+
+    FunDiffOrder  = convert(Int, _update_gradients) 
+    FunDiffOrder += convert(Int, _update_hessians) 
     GeoDiffOrder = max(required_geo_diff_order(mapping_type(ip_fun), FunDiffOrder), _update_detJdV)
     geo_mapping = GeometryMapping{GeoDiffOrder}(T, ip_geo.ip, qr)
     fun_values = FunctionValues{FunDiffOrder}(T, ip_fun, qr, ip_geo)
@@ -79,6 +86,7 @@ shape_gradient_type(cv::CellValues) = shape_gradient_type(cv.fun_values)
 
 @propagate_inbounds shape_value(cv::CellValues, q_point::Int, i::Int) = shape_value(cv.fun_values, q_point, i)
 @propagate_inbounds shape_gradient(cv::CellValues, q_point::Int, i::Int) = shape_gradient(cv.fun_values, q_point, i)
+@propagate_inbounds shape_hessian(cv::CellValues, q_point::Int, i::Int) = shape_hessian(cv.fun_values, q_point, i)
 @propagate_inbounds shape_symmetric_gradient(cv::CellValues, q_point::Int, i::Int) = shape_symmetric_gradient(cv.fun_values, q_point, i)
 
 # Access quadrature rule values 

src/FEValues/FaceValues.jl

@@ -41,9 +41,15 @@ struct FaceValues{FV, GM, FQR, detT, nT, V_FV<:AbstractVector{FV}, V_GM<:Abstrac
 end
 
 function FaceValues(::Type{T}, fqr::FaceQuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation{sdim} = default_geometric_interpolation(ip_fun); 
-        update_gradients::Union{Bool,Nothing} = nothing) where {T,sdim} 
-
-    FunDiffOrder = update_gradients === nothing ? 1 : convert(Int, update_gradients) # Logic must change when supporting update_hessian kwargs
+    update_gradients::Union{Bool,Nothing} = nothing, 
+    update_hessians::Union{Bool,Nothing} = nothing) where {T,sdim}
+
+    _update_gradients = update_gradients === nothing ? true : update_gradients
+    _update_hessians  = update_hessians  === nothing ? false : update_hessians
+    _update_hessians && @assert _update_gradients
+
+    FunDiffOrder  = convert(Int, _update_gradients) 
+    FunDiffOrder += convert(Int, _update_hessians) 
     GeoDiffOrder = max(required_geo_diff_order(mapping_type(ip_fun), FunDiffOrder), 1)
     # Not sure why the type-asserts are required here but not for CellValues, 
     # but they solve the type-instability for FaceValues
@@ -84,6 +90,7 @@ get_fun_values(fv::FaceValues) = @inbounds fv.fun_values[getcurrentface(fv)]
 
 @propagate_inbounds shape_value(fv::FaceValues, q_point::Int, i::Int) = shape_value(get_fun_values(fv), q_point, i)
 @propagate_inbounds shape_gradient(fv::FaceValues, q_point::Int, i::Int) = shape_gradient(get_fun_values(fv), q_point, i)
+@propagate_inbounds shape_hessian(fv::FaceValues, q_point::Int, i::Int) = shape_hessian(get_fun_values(fv), q_point, i)
 @propagate_inbounds shape_symmetric_gradient(fv::FaceValues, q_point::Int, i::Int) = shape_symmetric_gradient(get_fun_values(fv), q_point, i)
 
 """

src/FEValues/FunctionValues.jl

@@ -6,24 +6,34 @@
 #################################################################
 
 # Scalar, sdim == rdim                                                 sdim                     rdim
-typeof_N(   ::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = T
-typeof_dNdx(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Vec{dim, T}
-typeof_dNdξ(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Vec{dim, T}
+typeof_N(     ::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = T
+typeof_dNdx(  ::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Vec{dim, T}
+typeof_dNdξ(  ::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Vec{dim, T}
+typeof_d2Ndx2(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Tensor{2, dim, T}
+typeof_d2Ndξ2(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Tensor{2, dim, T}
 
 # Vector, vdim == sdim == rdim              vdim                            sdim                     rdim
-typeof_N(   ::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Vec{dim, T}
-typeof_dNdx(::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Tensor{2, dim, T}
-typeof_dNdξ(::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Tensor{2, dim, T}
+typeof_N(     ::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Vec{dim, T}
+typeof_dNdx(  ::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Tensor{2, dim, T}
+typeof_dNdξ(  ::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Tensor{2, dim, T}
+typeof_d2Ndx2(::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Tensor{3, dim, T}
+typeof_d2Ndξ2(::Type{T}, ::VectorInterpolation{dim}, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Tensor{3, dim, T}
 
 # Scalar, sdim != rdim (TODO: Use Vec if (s|r)dim <= 3?)
 typeof_N(   ::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, sdim, rdim} = T
 typeof_dNdx(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, sdim, rdim} = SVector{sdim, T}
 typeof_dNdξ(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, sdim, rdim} = SVector{rdim, T}
+typeof_d2Ndx2(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, sdim, rdim} = SMatrix{sdim, sdim, T, sdim*sdim}
+typeof_d2Ndξ2(::Type{T}, ::ScalarInterpolation, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, sdim, rdim} = SMatrix{rdim, rdim, T, rdim*rdim}
+
 
 # Vector, vdim != sdim != rdim (TODO: Use Vec/Tensor if (s|r)dim <= 3?)
 typeof_N(   ::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, vdim, sdim, rdim} = SVector{vdim, T}
-typeof_dNdx(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, vdim, sdim, rdim} = SMatrix{vdim, sdim, T}
-typeof_dNdξ(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, vdim, sdim, rdim} = SMatrix{vdim, rdim, T}
+typeof_dNdx(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, vdim, sdim, rdim} = SMatrix{vdim, sdim, T, vdim*sdim}
+typeof_dNdξ(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, vdim, sdim, rdim} = SMatrix{vdim, rdim, T, vdim*rdim}
+typeof_d2Ndx2(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, vdim, sdim, rdim} = SArray{Tuple{vdim, sdim, sdim}, T, 3, vdim*sdim*sdim}
+typeof_d2Ndξ2(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, vdim, sdim, rdim} = SArray{Tuple{vdim, rdim, rdim}, T, 3, vdim*rdim*rdim}
+
 
 """
     FunctionValues{DiffOrder}(::Type{T}, ip_fun, qr::QuadratureRule, ip_geo::VectorizedInterpolation)
@@ -33,36 +43,51 @@ for both the reference cell (precalculated) and the real cell (updated in `reini
 """
 FunctionValues
 
-struct FunctionValues{DiffOrder, IP, N_t, dNdx_t, dNdξ_t}
+struct FunctionValues{DiffOrder, IP, N_t, dNdx_t, dNdξ_t, d2Ndx2_t, d2Ndξ2_t}
     ip::IP          # ::Interpolation
     Nx::N_t         # ::AbstractMatrix{Union{<:Tensor,<:Number}}
     Nξ::N_t         # ::AbstractMatrix{Union{<:Tensor,<:Number}}
     dNdx::dNdx_t    # ::AbstractMatrix{Union{<:Tensor,<:StaticArray}} or Nothing
     dNdξ::dNdξ_t    # ::AbstractMatrix{Union{<:Tensor,<:StaticArray}} or Nothing
-    function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, ::Nothing, ::Nothing) where {N_t<:AbstractMatrix}
-        return new{0, typeof(ip), N_t, Nothing, Nothing}(ip, Nx, Nξ, nothing, nothing)
+    d2Ndx2::d2Ndx2_t   # ::AbstractMatrix{<:Tensor{2}}  Hessians of geometric shape functions in ref-domain
+    d2Ndξ2::d2Ndξ2_t   # ::AbstractMatrix{<:Tensor{2}}  Hessians of geometric shape functions in ref-domain
+    function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, ::Nothing, ::Nothing, ::Nothing, ::Nothing) where {N_t<:AbstractMatrix}
+        return new{0, typeof(ip), N_t, Nothing, Nothing, Nothing, Nothing}(ip, Nx, Nξ, nothing, nothing, nothing, nothing)
+    end
+    function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, dNdx::AbstractMatrix, dNdξ::AbstractMatrix, ::Nothing, ::Nothing) where {N_t<:AbstractMatrix}
+        return new{1, typeof(ip), N_t, typeof(dNdx), typeof(dNdξ), Nothing, Nothing}(ip, Nx, Nξ, dNdx, dNdξ, nothing, nothing)
     end
-    function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, dNdx::AbstractMatrix, dNdξ::AbstractMatrix) where {N_t<:AbstractMatrix}
-        return new{1, typeof(ip), N_t, typeof(dNdx), typeof(dNdξ)}(ip, Nx, Nξ, dNdx, dNdξ)
+    function FunctionValues(ip::Interpolation, Nx::N_t, Nξ::N_t, dNdx::AbstractMatrix, dNdξ::AbstractMatrix, d2Ndx2::AbstractMatrix, d2Ndξ2::AbstractMatrix) where {N_t<:AbstractMatrix}
+        return new{2, typeof(ip), N_t, typeof(dNdx), typeof(dNdξ), typeof(d2Ndx2), typeof(d2Ndξ2)}(ip, Nx, Nξ, dNdx, dNdξ, d2Ndx2, d2Ndξ2)
     end
 end
 function FunctionValues{DiffOrder}(::Type{T}, ip::Interpolation, qr::QuadratureRule, ip_geo::VectorizedInterpolation) where {DiffOrder, T}
+    sdim = n_components(ip_geo)
+    rdim = getdim(ip)
+    ((rdim != sdim) && DiffOrder>1) && error("CellValues/FunctionValues for embedded elements and higher order gradient not implemented")
+
     n_shape = getnbasefunctions(ip)
     n_qpoints = getnquadpoints(qr)
 
     Nξ = zeros(typeof_N(T, ip, ip_geo), n_shape, n_qpoints)
     Nx = isa(mapping_type(ip), IdentityMapping) ? Nξ : similar(Nξ)
-
-    if DiffOrder == 0
-        dNdξ = dNdx = nothing
-    elseif DiffOrder == 1
+    dNdξ = dNdx = d2Ndξ2 = d2Ndx2 = nothing
+
+    if DiffOrder >= 1
         dNdξ = zeros(typeof_dNdξ(T, ip, ip_geo),               n_shape, n_qpoints)
         dNdx = fill(zero(typeof_dNdx(T, ip, ip_geo)) * T(NaN), n_shape, n_qpoints)
-    else
+    end
+
+    if DiffOrder >= 2
+        d2Ndξ2 = zeros(typeof_d2Ndξ2(T, ip, ip_geo),               n_shape, n_qpoints)
+        d2Ndx2 = fill(zero(typeof_d2Ndx2(T, ip, ip_geo)) * T(NaN), n_shape, n_qpoints)
+    end
+
+    if DiffOrder > 2
         throw(ArgumentError("Currently only values and gradients can be updated in FunctionValues"))
     end
 
-    fv = FunctionValues(ip, Nx, Nξ, dNdx, dNdξ)
+    fv = FunctionValues(ip, Nx, Nξ, dNdx, dNdξ, d2Ndx2, d2Ndξ2)
     precompute_values!(fv, getpoints(qr)) # Separate function for qr point update in PointValues
     return fv
 end
@@ -73,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
@@ -167,6 +201,25 @@ end
     return nothing
 end
 
+@inline function apply_mapping!(funvals::FunctionValues{2}, ::IdentityMapping, q_point::Int, mapping_values, args...)
+    Jinv = calculate_Jinv(getjacobian(mapping_values))
+    H    = gethessian(mapping_values)
+    is_vector_valued = first(funvals.Nx) isa Vec
+    Jinv_otimesu_Jinv = is_vector_valued ? otimesu(Jinv, Jinv) : nothing
+    @inbounds for j in 1:getnbasefunctions(funvals)
+        dNdx = dothelper(funvals.dNdξ[j, q_point], Jinv)
+        if is_vector_valued #TODO - combine with helper function ? 
+            d2Ndx2 = (funvals.d2Ndξ2[j, q_point] - dNdx⋅H) ⊡ Jinv_otimesu_Jinv
+        else
+            d2Ndx2 = Jinv'⋅(funvals.d2Ndξ2[j, q_point] - dNdx⋅H)⋅Jinv
+        end
+
+        funvals.dNdx[j, q_point]   = dNdx
+        funvals.d2Ndx2[j, q_point] = d2Ndx2
+    end
+    return nothing
+end
+
 # TODO in PR798, apply_mapping! for 
 # * CovariantPiolaMapping
 # * ContravariantPiolaMapping