Implement hermitian polynomials

src/FEValues/CellValues.jl

@@ -108,7 +108,7 @@ function reinit!(cv::CellValues, cell::Union{AbstractCell, Nothing}, x::Abstract
     n_geom_basefuncs = getngeobasefunctions(geo_mapping)
 
     check_reinit_sdim_consistency(:CellValues, shape_gradient_type(cv), eltype(x))
-    if cell === nothing && !isa(mapping_type(fun_values), IdentityMapping)
+    if cell === nothing && !isa(mapping_type(fun_values), IdentityMapping) && !isa(mapping_type(fun_values), HermiteMapping)
         throw(ArgumentError("The cell::AbstractCell input is required to reinit! non-identity function mappings"))
     end
     if !checkbounds(Bool, x, 1:n_geom_basefuncs) || length(x) != n_geom_basefuncs

src/FEValues/FunctionValues.jl

@@ -142,6 +142,7 @@
 
 # Mapping types
 struct IdentityMapping end
+struct HermiteMapping end
 # Not yet implemented:
 # struct CovariantPiolaMapping end # PR798
 # struct ContravariantPiolaMapping end # PR798
@@ -157,7 +158,8 @@
 the function values and gradients from the reference cell
 to the physical cell geometry.
 """
-required_geo_diff_order(::IdentityMapping,           fun_diff_order::Int) = fun_diff_order
+required_geo_diff_order(::IdentityMapping, fun_diff_order::Int)   = fun_diff_order
+required_geo_diff_order(::HermiteMapping,  fun_diff_order::Int64) = fun_diff_order
 #required_geo_diff_order(::ContravariantPiolaMapping, fun_diff_order::Int) = 1 + fun_diff_order # PR798
 #required_geo_diff_order(::CovariantPiolaMapping,     fun_diff_order::Int) = 1 + fun_diff_order # PR798
 
@@ -220,6 +222,43 @@
     return nothing
 end
 
+@inline function apply_mapping!(funvals::FunctionValues{2}, ::HermiteMapping, q_point::Int, mapping_values, args...)
+    J = getjacobian(mapping_values)
+    Jinv = calculate_Jinv(J)
+    sdim, rdim = size(Jinv)
+    #TODO Generalize to 2d/3d
+    @assert sdim == 1 && rdim == 1
+
+    (rdim != sdim) && error("apply_mapping! for second order gradients and embedded elements not implemented")
+
+    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:2:getnbasefunctions(funvals)
+        dNdx = dothelper(funvals.dNdξ[j, q_point], Jinv)
+        if is_vector_valued
+            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.Nx[j, q_point] = funvals.Nξ[j, q_point]
+        funvals.dNdx[j, q_point]   = dNdx
+        funvals.d2Ndx2[j, q_point] = d2Ndx2
+    end
+
+    @inbounds for j in 2:2:getnbasefunctions(funvals)
+        N = funvals.Nξ[j, q_point]*det(J)
+        d2Ndx2 = funvals.d2Ndξ2[j, q_point] ⋅ Jinv
+
+        funvals.Nx[j, q_point]   = N
+        funvals.dNdx[j, q_point] = funvals.dNdξ[j, q_point]
+        funvals.d2Ndx2[j, q_point] = d2Ndx2
+    end
+
+    return nothing
+end
+
 # TODO in PR798, apply_mapping! for
 # * CovariantPiolaMapping
 # * ContravariantPiolaMapping

src/interpolations.jl

@@ -1594,6 +1594,31 @@
     throw(ArgumentError("no shape function $i for interpolation $ip"))
 end
 
+#################################
+# Hermitian dim 1 order 3 #
+#################################
+struct Hermitian{shape, order} <: ScalarInterpolation{shape, order} end 
+
+adjust_dofs_during_distribution(::Hermitian{RefLine,3}) = false
+vertexdof_indices(::Hermitian{RefLine,3}) = ((1,2),(3,4))
+getnbasefunctions(::Hermitian{RefLine,3}) = 4
+edgedof_indices(::Hermitian{RefLine,3}) = ((1,2,3,4),)
+
+function reference_coordinates(::Hermitian{RefLine,1})
+return [Vec{1, Float64}((-1.0,)),
+        Vec{1, Float64}(( 1.0,))] #Not sure
+end
+
+function Ferrite.reference_shape_value(ip::Hermitian{RefLine,3}, ξ::Vec{1}, i::Int)
+    ξ_x = (ξ[1]+1)/2
+    i == 1 && return 2*ξ_x^3 - 3*ξ_x^2 + 1
+    i == 2 && return 2ξ_x*(ξ_x^2 - 2*ξ_x + 1) #note scale with 2
+    i == 3 && return ξ_x^2*(3-2*ξ_x)
+    i == 4 && return 2ξ_x^2*(ξ_x-1) #note scale with 2
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+mapping_type(::Hermitian) = HermiteMapping()
+
 ##################################################
 # VectorizedInterpolation{<:ScalarInterpolation} #
 ##################################################

test/test_interpolations.jl

@@ -265,5 +265,54 @@ end
     @test_throws ArgumentError Ferrite.facetdof_interior_indices(ip)
 end
 
+@testset "Test Hermitian" begin
+
+    grid = generate_grid(Line, (1,), Vec(0.0), Vec(5.0))
+
+    ip = Ferrite.Hermitian{RefLine,3}()
+
+    dh = DofHandler(grid)
+    add!(dh, :u, ip)
+    close!(dh)
+
+    qr = QuadratureRule{RefLine}(5)
+    cv = CellValues(qr, ip; update_hessians = true)
+
+    function ebbeam(cv)
+        Ke = zeros(4,4)
+        fe = zeros(4)
+        for qp in 1:getnquadpoints(cv)
+            dV = getdetJdV(cv,qp)
+            for i in 1:getnbasefunctions(cv)
+                dNi = shape_hessian(cv,qp,i)
+                fe[i] += shape_value(cv, qp, i) * dV
+                for j in 1:getnbasefunctions(cv)
+                    dNj = shape_hessian(cv,qp,j)
+                    Ke[i,j] += dNi⊡dNj*dV
+                end
+            end
+        end
+        return Ke, fe
+    end
+    function ebbeam_anlytical(L)
+        ke = [12 6L -12 6L;
+         6L 4L^2 -6L 2L^2;
+         -12 -6L 12 -6L;
+         6L 2L^2 -6L 4L^2] * (1/L^3)
+        fe = [0.5L, 1/12*L^2, 0.5L, -1/12*L^2]
+        return ke, fe
+    end
+
+    coords = getcoordinates(grid, 1)
+    reinit!(cv, coords)
+    L = norm(coords[2] - coords[1])
+
+    ke, fe = ebbeam(cv)
+    kea, fea = ebbeam_anlytical(L)
+
+    @test ke ≈ kea
+    @test fe ≈ fea
+
+end
 
 end # testset

tmp.jl

@@ -0,0 +1,18 @@
+using Ferrite
+
+function get_2d_grid()
+    # GIVEN: two cells, a quad and a triangle sharing one face
+    cells = [
+        Quadrilateral((1, 2, 3, 4)),
+        Triangle((3, 2, 5))
+        ]
+    coords = zeros(Vec{2,Float64}, 5)
+    nodes = [Node(coord) for coord in zeros(Vec{2,Float64}, 5)]
+    return Grid(cells,nodes)
+end
+
+grid = get_2d_grid()
+dh = DofHandler(grid);
+sdh1 = SubDofHandler(dh, Set([1,2]))
+add!(sdh1, :u, Lagrange{Ferrite.RefAnyshape{2},1}())
+close!(dh)