diff --git a/docs/src/devdocs/FEValues.md b/docs/src/devdocs/FEValues.md
index e329e25c86..627c11b2c6 100644
--- a/docs/src/devdocs/FEValues.md
+++ b/docs/src/devdocs/FEValues.md
@@ -22,7 +22,7 @@ Ferrite.BCValues
 ## Internal utilities
 ```@docs
 Ferrite.embedding_det
-Ferrite.shape_value_type
+Ferrite.shape_value_type(::Ferrite.AbstractValues)
 Ferrite.shape_gradient_type
 Ferrite.ValuesUpdateFlags
 ```
diff --git a/docs/src/devdocs/interpolations.md b/docs/src/devdocs/interpolations.md
index 22c45ac3ef..5c3f441bc4 100644
--- a/docs/src/devdocs/interpolations.md
+++ b/docs/src/devdocs/interpolations.md
@@ -19,6 +19,7 @@ Ferrite.reference_shape_values!
 Ferrite.reference_shape_gradients!
 Ferrite.reference_shape_gradients_and_values!
 Ferrite.reference_shape_hessians_gradients_and_values!
+Ferrite.shape_value_type(ip::Interpolation, ::Type{T}) where T<:Number
 ```
 
 ### Required methods to implement for all subtypes of `Interpolation` to define a new finite element
diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl
index a89ffef1cb..84f59e62b7 100644
--- a/src/Dofs/DofHandler.jl
+++ b/src/Dofs/DofHandler.jl
@@ -930,13 +930,13 @@ function evaluate_at_grid_nodes(dh::DofHandler, u::AbstractVector, fieldname::Sy
 end
 
 # Internal method that have the vtk option to allocate the output differently
-function _evaluate_at_grid_nodes(dh::DofHandler, u::AbstractVector{T}, fieldname::Symbol, ::Val{vtk} = Val(false)) where {T, vtk}
+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 = ip isa ScalarInterpolation ? T : Vec{n_components(ip), T}
+    RT = shape_value_type(ip, T)
     if vtk
         # VTK output of solution field (or L2 projected scalar data)
         n_c = n_components(ip)
@@ -957,15 +957,10 @@ function _evaluate_at_grid_nodes(dh::DofHandler, u::AbstractVector{T}, fieldname
         ip_geo = geometric_interpolation(CT)
         local_node_coords = reference_coordinates(ip_geo)
         qr = QuadratureRule{getrefshape(ip)}(zeros(length(local_node_coords)), local_node_coords)
-        if ip isa VectorizedInterpolation
-            # TODO: Remove this hack when embedding works...
-            cv = CellValues(qr, ip.ip, ip_geo)
-        else
-            cv = CellValues(qr, ip, ip_geo)
-        end
+        cv = CellValues(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, RT)
+        _evaluate_at_grid_nodes!(data, sdh, u, cv, drange)
     end
     return data
 end
@@ -973,15 +968,9 @@ end
 # Loop over the cells and use shape functions to compute the value
 function _evaluate_at_grid_nodes!(
         data::Union{Vector, Matrix}, sdh::SubDofHandler,
-        u::AbstractVector{T}, cv::CellValues, drange::UnitRange, ::Type{RT}
-    ) where {T, RT}
+        u::AbstractVector{T}, cv::CellValues, drange::UnitRange
+    ) where {T}
     ue = zeros(T, length(drange))
-    # TODO: Remove this hack when embedding works...
-    if RT <: Vec && function_interpolation(cv) isa ScalarInterpolation
-        uer = reinterpret(RT, ue)
-    else
-        uer = ue
-    end
     for cell in CellIterator(sdh)
         # Note: We are only using the shape functions: no reinit!(cv, cell) necessary
         @assert getnquadpoints(cv) == length(cell.nodes)
@@ -989,7 +978,7 @@ function _evaluate_at_grid_nodes!(
             ue[i] = u[cell.dofs[I]]
         end
         for (qp, nodeid) in pairs(cell.nodes)
-            val = function_value(cv, qp, uer)
+            val = function_value(cv, qp, ue)
             if data isa Matrix # VTK
                 data[1:length(val), nodeid] .= val
                 data[(length(val) + 1):end, nodeid] .= 0 # purge the NaN
diff --git a/src/interpolations.jl b/src/interpolations.jl
index afc2119340..8bf9edcbda 100644
--- a/src/interpolations.jl
+++ b/src/interpolations.jl
@@ -57,6 +57,17 @@ n_components(::VectorInterpolation{vdim}) where {vdim} = vdim
 # Number of components that are allowed to prescribe in e.g. Dirichlet BC
 n_dbc_components(ip::Interpolation) = n_components(ip)
 
+"""
+    shape_value_type(ip::Interpolation, ::Type{T}) where T<:Number
+
+Return the type of `shape_value(ip::Interpolation, ξ::Vec, ib::Int)`.
+"""
+shape_value_type(::Interpolation, ::Type{T}) where {T <: Number}
+
+shape_value_type(::ScalarInterpolation, ::Type{T}) where {T <: Number} = T
+shape_value_type(::VectorInterpolation{vdim}, ::Type{T}) where {vdim, T <: Number} = Vec{vdim, T}
+#shape_value_type(::MatrixInterpolation, T::Type) = Tensor  #958
+
 # TODO: Add a fallback that errors if there are multiple dofs per edge/face instead to force
 #       interpolations to opt-out instead of silently do nothing.
 """