Split values into geometric mapping and function values

benchmark/helper.jl

@@ -66,17 +66,17 @@ function _generalized_ritz_galerkin_assemble_local_matrix(grid::Ferrite.Abstract
 end
 
 # Minimal Petrov-Galerkin type local assembly loop. We assume that both function spaces share the same integration rule. Test is applied from the left.
-function _generalized_petrov_galerkin_assemble_local_matrix(grid::Ferrite.AbstractGrid, cellvalues_shape::CellValues{<: Ferrite.InterpolationByDim{dim}}, f_shape, cellvalues_test::CellValues{<: Ferrite.InterpolationByDim{dim}}, f_test, op) where {dim}
+function _generalized_petrov_galerkin_assemble_local_matrix(grid::Ferrite.AbstractGrid, cellvalues_shape::CellValues, f_shape, cellvalues_test::CellValues, f_test, op)
     n_basefuncs_shape = getnbasefunctions(cellvalues_shape)
     n_basefuncs_test = getnbasefunctions(cellvalues_test)
     Ke = zeros(n_basefuncs_test, n_basefuncs_shape)
 
     #implicit assumption: Same geometry!
-    X_shape = zeros(Vec{dim,Float64}, Ferrite.getngeobasefunctions(cellvalues_shape))
+    X_shape = zeros(get_coordinate_type(grid), Ferrite.getngeobasefunctions(cellvalues_shape))
     getcoordinates!(X_shape, grid, 1)
     reinit!(cellvalues_shape, X_shape)
 
-    X_test = zeros(Vec{dim,Float64}, Ferrite.getngeobasefunctions(cellvalues_test))
+    X_test = zeros(get_coordinate_type(grid), Ferrite.getngeobasefunctions(cellvalues_test))
     getcoordinates!(X_test, grid, 1)
     reinit!(cellvalues_test, X_test)
 

docs/Project.toml

@@ -16,6 +16,7 @@ Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"

docs/src/devdocs/index.md

@@ -5,5 +5,5 @@ developing the library.
 
 ```@contents
 Depth = 1
-Pages = ["reference_cells.md", "interpolations.md", "elements.md", "dofhandler.md", "performance.md"]
+Pages = ["reference_cells.md", "interpolations.md", "elements.md", "mapping.md", "dofhandler.md", "performance.md"]
 ```

docs/src/devdocs/mapping.md

@@ -0,0 +1,103 @@
+# Mapping of finite elements
+The shape functions and gradients stored in an `FEValues` object, is reinitialized for each cell by calling the `reinit!` function. The main part of this calculation, considers how to map the functions described on the reference cell, to the actual cell.
+
+The geometric mapping of a finite element from the reference coordinates to the real coordinates is shown in the following illustration. 
+
+![mapping_figure](../assets/fe_mapping.svg)
+
+This mapping is given by the geometric shape functions, $\hat{N}_i^g(\boldsymbol{\xi})$, such that 
+```math
+\begin{align*}
+    \boldsymbol{x}(\boldsymbol{\xi}) =& \sum_{i}^N \hat{\boldsymbol{x}}_i \hat{N}_i^g(\boldsymbol{\xi}) \\
+    \boldsymbol{J} :=& \frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}\boldsymbol{\xi}} = \sum_{i}^N \hat{\boldsymbol{x}}_i \otimes \frac{\mathrm{d} \hat{N}_i^g}{\mathrm{d}\boldsymbol{\xi}}\\
+    \boldsymbol{\mathcal{H}} :=&
+    \frac{\mathrm{d} \boldsymbol{J}}{\mathrm{d} \boldsymbol{\xi}} = \sum_{\alpha=1}^N \hat{\boldsymbol{x}}_\alpha \otimes \frac{\mathrm{d}^2 \hat{N}^g_\alpha}{\mathrm{d} \boldsymbol{\xi}^2}
+\end{align*}
+```
+where the defined $\boldsymbol{J}$ is the jacobian of the mapping, and in some cases we will also need the corresponding hessian, $\boldsymbol{\mathcal{H}}$ (3rd order tensor).
+
+We require that the mapping from reference coordinates to real coordinates is [diffeomorphic](https://en.wikipedia.org/wiki/Diffeomorphism), meaning that we can express $\boldsymbol{x} = \boldsymbol{x}(\boldsymbol{\xi}(\boldsymbol{x}))$, such that
+```math
+\begin{align*}
+    \frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}\boldsymbol{x}} = \boldsymbol{I} &= \frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}\boldsymbol{\xi}} \cdot \frac{\mathrm{d}\boldsymbol{\xi}}{\mathrm{d}\boldsymbol{x}} 
+    \quad\Rightarrow\quad 
+    \frac{\mathrm{d}\boldsymbol{\xi}}{\mathrm{d}\boldsymbol{x}} = \left[\frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}\boldsymbol{\xi}}\right]^{-1} = \boldsymbol{J}^{-1}
+\end{align*}
+```
+Depending on the function interpolation, we may want different types of mappings to conserve certain properties of the fields. This results in the different mapping types described below.
+
+## Identity mapping
+`Ferrite.IdentityMapping`
+
+For scalar fields, we always use scalar base functions. For tensorial fields (non-scalar, e.g. vector-fields), the base functions can be constructed from scalar base functions, by using e.g. `VectorizedInterpolation`. From the perspective of the mapping, however, each component is mapped as an individual scalar base function. And for scalar base functions, we only require that the value of the base function is invariant to the element shape (real coordinate), and only depends on the reference coordinate, i.e. 
+```math
+\begin{align*}
+    N(\boldsymbol{x}) &= \hat{N}(\boldsymbol{\xi}(\boldsymbol{x}))\nonumber \\
+    \mathrm{grad}(N(\boldsymbol{x})) &= \frac{\mathrm{d}\hat{N}}{\mathrm{d}\boldsymbol{\xi}} \cdot \boldsymbol{J}^{-1}
+\end{align*}
+```
+
+## Covariant Piola mapping, H(curl)
+`Ferrite.CovariantPiolaMapping`
+
+The covariant Piola mapping of a vectorial base function preserves the tangential components. For the value, the mapping is defined as 
+```math
+\begin{align*}
+    \boldsymbol{N}(\boldsymbol{x}) = \boldsymbol{J}^{-\mathrm{T}} \cdot \hat{\boldsymbol{N}}(\boldsymbol{\xi}(\boldsymbol{x}))
+\end{align*}
+```
+which yields the gradient,
+```math
+\begin{align*}
+    \mathrm{grad}(\boldsymbol{N}(\boldsymbol{x})) &= \boldsymbol{J}^{-T} \cdot \frac{\mathrm{d} \hat{\boldsymbol{N}}}{\mathrm{d} \boldsymbol{\xi}} \cdot \boldsymbol{J}^{-1} - \boldsymbol{J}^{-T} \cdot \left[\hat{\boldsymbol{N}}(\boldsymbol{\xi}(\boldsymbol{x}))\cdot \boldsymbol{J}^{-1} \cdot \boldsymbol{\mathcal{H}}\cdot \boldsymbol{J}^{-1}\right]
+\end{align*}
+```
+
+!!! details "Derivation"
+    Expressing the gradient, $\mathrm{grad}(\boldsymbol{N})$, in index notation,
+    ```math
+    \begin{align*}
+        \frac{\mathrm{d} N_i}{\mathrm{d} x_j} &= \frac{\mathrm{d}}{\mathrm{d} x_j} \left[J^{-\mathrm{T}}_{ik} \hat{N}_k\right] = \frac{\mathrm{d} J^{-\mathrm{T}}_{ik}}{\mathrm{d} x_j} \hat{N}_k + J^{-\mathrm{T}}_{ik}  \frac{\mathrm{d} \hat{N}_k}{\mathrm{d} \xi_l} J_{lj}^{-1}
+    \end{align*}
+    ```
+
+    Except for a few elements, $\boldsymbol{J}$ varies as a function of $\boldsymbol{x}$. The derivative can be calculated as 
+    ```math
+    \begin{align*}
+        \frac{\mathrm{d} J^{-\mathrm{T}}_{ik}}{\mathrm{d} x_j} &= \frac{\mathrm{d} J^{-\mathrm{T}}_{ik}}{\mathrm{d} J_{mn}} \frac{\mathrm{d} J_{mn}}{\mathrm{d} x_j} = - J_{km}^{-1} J_{in}^{-T} \frac{\mathrm{d} J_{mn}}{\mathrm{d} x_j} \nonumber \\
+        \frac{\mathrm{d} J_{mn}}{\mathrm{d} x_j} &= \mathcal{H}_{mno} J_{oj}^{-1}
+    \end{align*}
+    ```
+
+## Contravariant Piola mapping, H(div)
+`Ferrite.ContravariantPiolaMapping`
+
+The covariant Piola mapping of a vectorial base function preserves the normal components. For the value, the mapping is defined as 
+```math
+\begin{align*}
+    \boldsymbol{N}(\boldsymbol{x}) = \frac{\boldsymbol{J}}{\det(\boldsymbol{J})} \cdot \hat{\boldsymbol{N}}(\boldsymbol{\xi}(\boldsymbol{x}))
+\end{align*}
+```
+This gives the gradient
+```math
+\begin{align*}
+    \mathrm{grad}(\boldsymbol{N}(\boldsymbol{x})) = [\boldsymbol{\mathcal{H}}\cdot\boldsymbol{J}^{-1}] : \frac{[\boldsymbol{I} \underline{\otimes} \boldsymbol{I}] \cdot \hat{\boldsymbol{N}}}{\det(\boldsymbol{J})}
+    - \left[\frac{\boldsymbol{J} \cdot \hat{\boldsymbol{N}}}{\det(\boldsymbol{J})}\right] \otimes \left[\boldsymbol{J}^{-T} : \boldsymbol{\mathcal{H}} \cdot \boldsymbol{J}^{-1}\right]
+    + \boldsymbol{J} \cdot \frac{\mathrm{d} \hat{\boldsymbol{N}}}{\mathrm{d} \boldsymbol{\xi}} \cdot \frac{\boldsymbol{J}^{-1}}{\det(\boldsymbol{J})}
+\end{align*}
+```
+
+!!! details "Derivation"
+    Expressing the gradient, $\mathrm{grad}(\boldsymbol{N})$, in index notation,
+    ```math
+    \begin{align*}
+        \frac{\mathrm{d} N_i}{\mathrm{d} x_j} &= \frac{\mathrm{d}}{\mathrm{d} x_j} \left[\frac{J_{ik}}{\det(\boldsymbol{J})} \hat{N}_k\right] =\nonumber\\
+        &= \frac{\mathrm{d} J_{ik}}{\mathrm{d} x_j} \frac{\hat{N}_k}{\det(\boldsymbol{J})} 
+        - \frac{\mathrm{d} \det(\boldsymbol{J})}{\mathrm{d} x_j} \frac{J_{ik} \hat{N}_k}{\det(\boldsymbol{J})^2}
+        + \frac{J_{ik}}{\det(\boldsymbol{J})}  \frac{\mathrm{d} \hat{N}_k}{\mathrm{d} \xi_l} J_{lj}^{-1} \\
+        &= \mathcal{H}_{ikl} J^{-1}_{lj} \frac{\hat{N}_k}{\det(\boldsymbol{J})} 
+        - J^{-T}_{mn} \mathcal{H}_{mnl} J^{-1}_{lj} \frac{J_{ik} \hat{N}_k}{\det(\boldsymbol{J})}
+        + \frac{J_{ik}}{\det(\boldsymbol{J})}  \frac{\mathrm{d} \hat{N}_k}{\mathrm{d} \xi_l} J_{lj}^{-1}
+    \end{align*}
+    ```
+

docs/src/literate-howto/postprocessing.jl

@@ -46,7 +46,7 @@ include("../tutorials/heat_equation.jl");
 # Next we define a function that computes the heat flux for each integration point in the domain.
 # Fourier's law is adopted, where the conductivity tensor is assumed to be isotropic with unit
 # conductivity ``\lambda = 1 ⇒ q = - \nabla u``, where ``u`` is the temperature.
-function compute_heat_fluxes(cellvalues::CellValues{<:ScalarInterpolation}, dh::DofHandler, a::AbstractVector{T}) where T
+function compute_heat_fluxes(cellvalues::CellValues, dh::DofHandler, a::AbstractVector{T}) where T
 
     n = getnbasefunctions(cellvalues)
     cell_dofs = zeros(Int, n)

docs/src/literate-tutorials/incompressible_elasticity.jl

@@ -85,9 +85,9 @@ end
 # use a `PseudoBlockArray` from `BlockArrays.jl`.
 
 function doassemble(
-    cellvalues_u::CellValues{<:VectorInterpolation},
-    cellvalues_p::CellValues{<:ScalarInterpolation},
-    facevalues_u::FaceValues{<:VectorInterpolation},
+    cellvalues_u::CellValues,
+    cellvalues_p::CellValues,
+    facevalues_u::FaceValues,
     K::SparseMatrixCSC, grid::Grid, dh::DofHandler, mp::LinearElasticity
 )
     f = zeros(ndofs(dh))

docs/src/literate-tutorials/linear_shell.jl

@@ -264,6 +264,8 @@ end;
 # ##### Main element routine
 # Below is the main routine that calculates the stiffness matrix of the shell element. 
 # Since it is a so called degenerate shell element, the code is similar to that for an standard continuum element.
+shape_reference_gradient(cv::CellValues, q_point, i) = cv.fun_values.dNdξ[i, q_point]
+
 function integrate_shell!(ke, cv, qr_ooplane, X, data)
     nnodes = getnbasefunctions(cv)
     ndofs = nnodes*5
@@ -281,9 +283,9 @@ function integrate_shell!(ke, cv, qr_ooplane, X, data)
     ef1, ef2, ef3 = fiber_coordsys(p)
 
     for iqp in 1:getnquadpoints(cv)
-        N = cv.N[:,iqp]
-        dNdξ = cv.dNdξ[:,iqp]
-        dNdx = cv.dNdx[:,iqp]
+        N = [shape_value(cv, iqp, i) for i in 1:nnodes]
+        dNdξ = [shape_reference_gradient(cv, iqp, i) for i in 1:nnodes]
+        dNdx = [shape_gradient(cv, iqp, i) for i in 1:nnodes]
 
         for oqp in 1:length(qr_ooplane.weights)
             ζ = qr_ooplane.points[oqp][1]

src/Dofs/DofHandler.jl

@@ -915,7 +915,7 @@ function _evaluate_at_grid_nodes!(data::Union{Vector,Matrix}, sdh::SubDofHandler
         u::Vector{T}, cv::CellValues, drange::UnitRange, ::Type{RT}) where {T, RT}
     ue = zeros(T, length(drange))
     # TODO: Remove this hack when embedding works...
-    if RT <: Vec && cv isa CellValues{<:ScalarInterpolation}
+    if RT <: Vec && get_function_interpolation(cv) isa ScalarInterpolation
         uer = reinterpret(RT, ue)
     else
         uer = ue

src/FEValues/FunctionValues.jl

@@ -0,0 +1,172 @@
+#################################################################
+# Note on dimensions:                                           #
+# sdim = spatial dimension (dimension of the grid nodes)        #
+# rdim = reference dimension (dimension in isoparametric space) #
+# vdim = vector dimension (dimension of the field)              #
+#################################################################
+
+# Helpers to get the correct types for FunctionValues for the given function and, if needed, geometric interpolations.
+struct SInterpolationDims{rdim,sdim} end
+struct VInterpolationDims{rdim,sdim,vdim} end
+function InterpolationDims(::ScalarInterpolation, ip_geo::VectorizedInterpolation{sdim}) where sdim
+    return SInterpolationDims{getdim(ip_geo),sdim}()
+end
+function InterpolationDims(::VectorInterpolation{vdim}, ip_geo::VectorizedInterpolation{sdim}) where {vdim,sdim}
+    return VInterpolationDims{getdim(ip_geo),sdim,vdim}()
+end
+
+typeof_N(::Type{T}, ::SInterpolationDims) where T = T
+typeof_N(::Type{T}, ::VInterpolationDims{<:Any,dim,dim}) where {T,dim} = Vec{dim,T}
+typeof_N(::Type{T}, ::VInterpolationDims{<:Any,<:Any,vdim}) where {T,vdim} = SVector{vdim,T} # Why not ::Vec here?
+
+typeof_dNdx(::Type{T}, ::SInterpolationDims{dim,dim}) where {T,dim} = Vec{dim,T}
+typeof_dNdx(::Type{T}, ::SInterpolationDims{<:Any,sdim}) where {T,sdim} = SVector{sdim,T} # Why not ::Vec here?
+typeof_dNdx(::Type{T}, ::VInterpolationDims{dim,dim,dim}) where {T,dim} = Tensor{2,dim,T}
+typeof_dNdx(::Type{T}, ::VInterpolationDims{<:Any,sdim,vdim}) where {T,sdim,vdim} = SMatrix{vdim,sdim,T} # If vdim=sdim!=rdim Tensor would be possible...
+
+typeof_dNdξ(::Type{T}, ::SInterpolationDims{dim,dim}) where {T,dim} = Vec{dim,T}
+typeof_dNdξ(::Type{T}, ::SInterpolationDims{rdim}) where {T,rdim} = SVector{rdim,T} # Why not ::Vec here?
+typeof_dNdξ(::Type{T}, ::VInterpolationDims{dim,dim,dim}) where {T,dim} = Tensor{2,dim,T}
+typeof_dNdξ(::Type{T}, ::VInterpolationDims{rdim,<:Any,vdim}) where {T,rdim,vdim} = SMatrix{vdim,rdim,T} # If vdim=rdim!=sdim Tensor would be possible...
+
+struct FunctionValues{IP, N_t, dNdx_t, dNdξ_t}
+    ip::IP          # ::Interpolation
+    N_x::N_t        # ::AbstractMatrix{Union{<:Tensor,<:Number}}
+    N_ξ::N_t        # ::AbstractMatrix{Union{<:Tensor,<:Number}}
+    dNdx::dNdx_t    # ::AbstractMatrix{Union{<:Tensor,<:StaticArray}}
+    dNdξ::dNdξ_t    # ::AbstractMatrix{Union{<:Tensor,<:StaticArray}}
+end
+function FunctionValues(::Type{T}, ip::Interpolation, qr::QuadratureRule, ip_geo::VectorizedInterpolation) where T
+    ip_dims = InterpolationDims(ip, ip_geo)
+    N_t = typeof_N(T, ip_dims)
+    dNdx_t = typeof_dNdx(T, ip_dims)
+    dNdξ_t = typeof_dNdξ(T, ip_dims)
+    n_shape = getnbasefunctions(ip)
+    n_qpoints = getnquadpoints(qr)
+
+    N_ξ  = zeros(N_t,    n_shape, n_qpoints)
+    if isa(get_mapping_type(ip), IdentityMapping)
+        N_x = N_ξ
+    else
+        N_x  = zeros(N_t,    n_shape, n_qpoints)
+    end
+    dNdξ = zeros(dNdξ_t, n_shape, n_qpoints)
+    dNdx = fill(zero(dNdx_t) * T(NaN), n_shape, n_qpoints)
+    fv = FunctionValues(ip, N_x, N_ξ, dNdx, dNdξ)
+    precompute_values!(fv, qr) # Precompute N and dNdξ
+    return fv
+end
+function Base.copy(v::FunctionValues)
+    N_ξ_copy = copy(v.N_ξ)
+    N_x_copy = v.N_ξ === v.N_x ? N_ξ_copy : copy(v.N_x) # Preserve aliasing
+    return FunctionValues(copy(v.ip), N_x_copy, N_ξ_copy, copy(v.dNdx), copy(v.dNdξ))
+end
+
+function precompute_values!(fv::FunctionValues, qr::QuadratureRule)
+    n_shape = getnbasefunctions(fv.ip)
+    for (qp, ξ) in pairs(getpoints(qr))
+        for i in 1:n_shape
+            fv.dNdξ[i, qp], fv.N_ξ[i, qp] = shape_gradient_and_value(fv.ip, ξ, i)
+        end
+    end
+end
+
+getnbasefunctions(funvals::FunctionValues) = size(funvals.N_x, 1)
+@propagate_inbounds shape_value(funvals::FunctionValues, q_point::Int, base_func::Int) = funvals.N_x[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_symmetric_gradient(funvals::FunctionValues, q_point::Int, base_func::Int) = symmetric(shape_gradient(funvals, q_point, base_func))
+
+get_function_interpolation(funvals::FunctionValues) = funvals.ip
+
+shape_value_type(funvals::FunctionValues) = eltype(funvals.N_x)
+shape_gradient_type(funvals::FunctionValues) = eltype(funvals.dNdx)
+
+
+# Checks that the user provides the right dimension of coordinates to reinit! methods to ensure good error messages if not
+sdim_from_gradtype(::Type{<:AbstractTensor{<:Any,sdim}}) where sdim = sdim
+sdim_from_gradtype(::Type{<:SVector{sdim}}) where sdim = sdim
+sdim_from_gradtype(::Type{<:SMatrix{<:Any,sdim}}) where sdim = sdim
+
+# For performance, these must be fully inferrable for the compiler.
+# args: valname (:CellValues or :FaceValues), shape_gradient_type, eltype(x)
+function check_reinit_sdim_consistency(valname, gradtype, ::Type{<:Vec{sdim}}) where {sdim}
+    check_reinit_sdim_consistency(valname, Val(sdim_from_gradtype(gradtype)), Val(sdim))
+end
+check_reinit_sdim_consistency(_, ::Val{sdim}, ::Val{sdim}) where sdim = nothing
+function check_reinit_sdim_consistency(valname, ::Val{sdim_val}, ::Val{sdim_x}) where {sdim_val, sdim_x}
+    error("""The $valname and coordinates have different spatial dimensions, $sdim_val and $sdim_x.
+        Could it be that you forgot to vectorize the geometric interpolation when constructing the $valname?""")
+end
+
+# Mapping types 
+struct IdentityMapping end 
+struct CovariantPiolaMapping end 
+struct ContravariantPiolaMapping end 
+# Not yet implemented:
+# struct DoubleCovariantPiolaMapping end 
+# struct DoubleContravariantPiolaMapping end 
+
+get_mapping_type(fv::FunctionValues) = get_mapping_type(fv.ip)
+
+requires_hessian(::IdentityMapping) = false
+requires_hessian(::ContravariantPiolaMapping) = true
+requires_hessian(::CovariantPiolaMapping) = true
+
+# Support for embedded elements
+calculate_Jinv(J::Tensor{2}) = inv(J)
+calculate_Jinv(J::SMatrix) = pinv(J)
+
+# Hotfix to get the dots right for embedded elements until mixed tensors are merged.
+# Scalar/Vector interpolations with sdim == rdim (== vdim)
+@inline dothelper(A, B) = A ⋅ B
+# Vector interpolations with sdim == rdim != vdim
+@inline dothelper(A::SMatrix{vdim, dim}, B::Tensor{2, dim}) where {vdim, dim} = A * SMatrix{dim, dim}(B)
+# Scalar interpolations with sdim > rdim
+@inline dothelper(A::SVector{rdim}, B::SMatrix{rdim, sdim}) where {rdim, sdim} = B' * A
+# Vector interpolations with sdim > rdim
+@inline dothelper(B::SMatrix{vdim, rdim}, A::SMatrix{rdim, sdim}) where {vdim, rdim, sdim} = B * A
+
+@inline function apply_mapping!(funvals::FunctionValues, args...)
+    return apply_mapping!(funvals, get_mapping_type(funvals), args...)
+end
+
+@inline function apply_mapping!(funvals::FunctionValues, ::IdentityMapping, q_point::Int, mapping_values, args...)
+    Jinv = calculate_Jinv(getjacobian(mapping_values))
+    @inbounds for j in 1:getnbasefunctions(funvals)
+        #funvals.dNdx[j, q_point] = funvals.dNdξ[j, q_point] ⋅ Jinv # TODO via Tensors.jl
+        funvals.dNdx[j, q_point] = dothelper(funvals.dNdξ[j, q_point], Jinv)
+    end
+    return nothing
+end
+
+@inline function apply_mapping!(funvals::FunctionValues, ::CovariantPiolaMapping, q_point::Int, mapping_values, cell)
+    H = gethessian(mapping_values)
+    Jinv = inv(getjacobian(mapping_values))
+    @inbounds for j in 1:getnbasefunctions(funvals)
+        d = get_direction(funvals.ip, j, cell)
+        dNdξ = funvals.dNdξ[j, q_point]
+        N_ξ = funvals.N_ξ[j, q_point]
+        funvals.N_x[j, q_point] = d*(N_ξ ⋅ Jinv)
+        funvals.dNdx[j, q_point] = d*(Jinv' ⋅ dNdξ ⋅ Jinv - Jinv' ⋅ (N_ξ ⋅ Jinv ⋅ H ⋅ Jinv))
+    end
+    return nothing
+end
+
+@inline function apply_mapping!(funvals::FunctionValues, ::ContravariantPiolaMapping, q_point::Int, mapping_values, cell)
+    H = gethessian(mapping_values)
+    J = getjacobian(mapping_values)
+    Jinv = inv(J)
+    detJ = det(J)
+    I2 = one(J)
+    H_Jinv = H⋅Jinv
+    A1 = (H_Jinv ⊡ (otimesl(I2,I2))) / detJ
+    A2 = (Jinv' ⊡ H_Jinv) / detJ
+    @inbounds for j in 1:getnbasefunctions(funvals)
+        d = get_direction(funvals.ip, j, cell)
+        dNdξ = funvals.dNdξ[j, q_point]
+        N_ξ = funvals.N_ξ[j, q_point]
+        funvals.N_x[j, q_point] = d*(J ⋅ N_ξ)/detJ
+        funvals.dNdx[j, q_point] = d*(J ⋅ dNdξ ⋅ Jinv/detJ + A1 ⋅ N_ξ - (J ⋅ N_ξ) ⊗ A2)
+    end
+    return nothing
+end