Higher order gradients in CellValues

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

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`
 

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

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)
 
 """