Fix wrong use of shape_value_type in evaluate_at_grid_nodes

Minor additions to interpolation testing
Ferrite-FEM · Feb 24, 2025 · 048cc38 · 048cc38
1 parent 6eead25
commit 048cc38
Show file tree

Hide file tree

Showing 2 changed files with 14 additions and 9 deletions.
diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl
@@ -933,14 +933,17 @@ function evaluate_at_grid_nodes(dh::DofHandler, u::AbstractVector, fieldname::Sy
     return _evaluate_at_grid_nodes(dh, u, fieldname)
 end
 
+_function_value_type(::ScalarInterpolation, ::AbstractVector{T}) where {T} = T
+_function_value_type(::VectorInterpolation{vdim}, ::AbstractVector{T}) where {vdim, T <: Number} = Vec{vdim, T}
+
 # Internal method that have the vtk option to allocate the output differently
 function _evaluate_at_grid_nodes(dh::DofHandler{sdim}, u::AbstractVector{T}, fieldname::Symbol, ::Val{vtk} = Val(false)) where {T, vtk, sdim}
     # Make sure the field exists
     fieldname ∈ getfieldnames(dh) || error("Field $fieldname not found.")
     # Figure out the return type (scalar or vector)
     field_idx = find_field(dh, fieldname)
     ip = getfieldinterpolation(dh, field_idx)
-    RT = shape_value_type(ip, T)
+    RT = _function_value_type(ip, typeof(u))
     if vtk
         # VTK output of solution field (or L2 projected scalar data)
         n_c = n_components(ip)
@@ -960,8 +963,8 @@ function _evaluate_at_grid_nodes(dh::DofHandler{sdim}, u::AbstractVector{T}, fie
         ip = getfieldinterpolation(sdh, field_idx)
         ip_geo = geometric_interpolation(CT)
         local_node_coords = reference_coordinates(ip_geo)
-        qr = QuadratureRule{getrefshape(ip)}(zeros(length(local_node_coords)), local_node_coords)
-        cv = CellValues(qr, ip, ip_geo^sdim; update_gradients = false, update_hessians = false, update_detJdV = false)
+        qr = QuadratureRule{getrefshape(ip)}(T, zeros(length(local_node_coords)), local_node_coords)
+        cv = CellValues(T, qr, ip, ip_geo^sdim; update_gradients = false, update_hessians = false, update_detJdV = false)
         drange = dof_range(sdh, field_idx)
         # Function barrier
         _evaluate_at_grid_nodes!(data, sdh, u, cv, drange)

diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl
@@ -52,14 +52,16 @@ function test_interpolation_properties(ip::Interpolation{RefShape, FunOrder}) wh
         # Test A-D
         _test_interpolation_properties(dof_data, refshape_data)
 
-        # Test E: All base functions implemented.
+        # Test E: All base functions implemented and infers correct types
         # Argument errors for 0th and n+1 indices.
-        ξ = zero(Vec{Ferrite.getrefdim(ip)})
-        @test_throws ArgumentError Ferrite.reference_shape_value(ip, ξ, 0)
-        for i in 1:getnbasefunctions(ip)
-            @test Ferrite.reference_shape_value(ip, ξ, i) isa Ferrite.shape_value_type(ip, Float64)
+        for T in (Float64, Float32)
+            ξ = zero(Vec{Ferrite.getrefdim(ip), T})
+            @test_throws ArgumentError Ferrite.reference_shape_value(ip, ξ, 0)
+            for i in 1:getnbasefunctions(ip)
+                @test (@inferred Ferrite.reference_shape_value(ip, ξ, i)) isa Ferrite.shape_value_type(ip, T)
+            end
+            @test_throws ArgumentError Ferrite.reference_shape_value(ip, ξ, getnbasefunctions(ip) + 1)
         end
-        @test_throws ArgumentError Ferrite.reference_shape_value(ip, ξ, getnbasefunctions(ip) + 1)
     end
 end