Add direct computation of mapping values (Ferrite-FEM#1018)

src/FEValues/GeometryMapping.jl

@@ -126,15 +126,15 @@ end
 end
 
 @inline function calculate_mapping(geo_mapping::GeometryMapping{1}, q_point::Int, x::AbstractVector{<:Vec})
-    fecv_J = zero(otimes_returntype(eltype(x), eltype(geo_mapping.dMdξ)))
+    J = zero(otimes_returntype(eltype(x), eltype(geo_mapping.dMdξ)))
     @inbounds for j in 1:getngeobasefunctions(geo_mapping)
-        #fecv_J += x[j] ⊗ geo_mapping.dMdξ[j, q_point]
-        fecv_J += otimes_helper(x[j], geo_mapping.dMdξ[j, q_point])
+        # J += x[j] ⊗ geo_mapping.dMdξ[j, q_point]
+        J += otimes_helper(x[j], geo_mapping.dMdξ[j, q_point])
     end
-    return MappingValues(fecv_J, nothing)
+    return MappingValues(J, nothing)
 end
 
-@inline function calculate_mapping(geo_mapping::GeometryMapping{2}, q_point, x)
+@inline function calculate_mapping(geo_mapping::GeometryMapping{2}, q_point::Int, x::AbstractVector{<:Vec})
     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)")
@@ -146,6 +146,37 @@ end
     return MappingValues(J, H)
 end
 
+@inline function calculate_mapping(gip::ScalarInterpolation, ξ::Vec, x::AbstractVector{<:Vec}, ::Val{0})
+    return MappingValues(nothing, nothing)
+end
+
+@inline function calculate_mapping(gip::ScalarInterpolation, ξ::Vec{rdim,T}, x::AbstractVector{<:Vec{sdim}}, ::Val{1}) where {T,rdim, sdim}
+    n_basefuncs = getnbasefunctions(gip)
+    @boundscheck checkbounds(x, Base.OneTo(n_basefuncs))
+
+    J = zero(otimes_returntype(Vec{sdim,T}, Vec{rdim,T}))
+    @inbounds for j in 1:n_basefuncs
+        dMdξ = reference_shape_gradient(gip, ξ, j)
+        # J += x[j] ⊗ dMdξ # https://github.com/Ferrite-FEM/Tensors.jl/pull/188
+        J += otimes_helper(x[j], dMdξ)
+    end
+    return MappingValues(J, nothing)
+end
+
+@inline function calculate_mapping(gip::ScalarInterpolation, ξ::Vec{rdim,T}, x::AbstractVector{<:Vec{sdim}}, ::Val{2}) where {T,rdim, sdim}
+    n_basefuncs = getnbasefunctions(gip)
+    @boundscheck checkbounds(x, Base.OneTo(n_basefuncs))
+    (rdim != sdim) && error("hessian for embedded elements not implemented (rdim=$rdim, sdim=$sdim)")
+    J = zero(otimes_returntype(Vec{sdim,T}, Vec{rdim,T}))
+    H = zero(otimes_returntype(eltype(x), typeof(J)))
+    @inbounds for j in 1:n_basefuncs
+        d2Mdξ2, dMdξ, _ = reference_shape_hessian_gradient_and_value(gip, ξ, j)
+        J += x[j] ⊗ dMdξ
+        H += x[j] ⊗ d2Mdξ2
+    end
+    return MappingValues(J, H)
+end
+
 calculate_detJ(J::Tensor{2}) = det(J)
 calculate_detJ(J::SMatrix) = embedding_det(J)
 

test/test_cellvalues.jl

@@ -167,6 +167,27 @@
     end
 end
 
+@testset "GeometryMapping" begin
+    grid = generate_grid(Quadrilateral, (1,1))
+    cc = first(CellIterator(grid))
+
+    qr = QuadratureRule{RefQuadrilateral}(1)
+    ξ = first(Ferrite.getpoints(qr))
+    ip = Lagrange{RefQuadrilateral,1}()
+
+    cv0 = CellValues(Float64, qr, ip, ip^2; update_detJdV=false, update_gradients=false, update_hessians=false)
+    reinit!(cv0, cc)
+    @test Ferrite.calculate_mapping(cv0.geo_mapping, 1, cc.coords) == Ferrite.calculate_mapping(ip, ξ, cc.coords, Val(0))
+
+    cv1 = CellValues(Float64, qr, ip, ip^2; update_detJdV=false, update_gradients=true,  update_hessians=false)
+    reinit!(cv1, cc)
+    @test Ferrite.calculate_mapping(cv1.geo_mapping, 1, cc.coords) == Ferrite.calculate_mapping(ip, ξ, cc.coords, Val(1))
+
+    cv2 = CellValues(Float64, qr, ip, ip^2; update_detJdV=false, update_gradients=false, update_hessians=true)
+    reinit!(cv2, cc)
+    @test Ferrite.calculate_mapping(cv2.geo_mapping, 1, cc.coords) == Ferrite.calculate_mapping(ip, ξ, cc.coords, Val(2))
+end
+
 @testset "#265: error message for incompatible geometric interpolation" begin
     dim = 1
     deg = 1