diff --git a/Project.toml b/Project.toml
index c18ee1634a..5da10c89b7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ version = "0.3.14"
 
 [deps]
 EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
diff --git a/src/FEValues/common_values.jl b/src/FEValues/common_values.jl
index c3a42cdb19..34950300b8 100644
--- a/src/FEValues/common_values.jl
+++ b/src/FEValues/common_values.jl
@@ -296,13 +296,13 @@ Compute the spatial coordinate in a given quadrature point. `cell_coordinates` c
 The coordinate is computed, using the geometric interpolation, as
 ``\\mathbf{x} = \\sum\\limits_{i = 1}^n M_i (\\mathbf{\\xi}) \\mathbf{\\hat{x}}_i``
 """
-spatial_coordinate(ip::VectorizedInterpolation, ξ::Vec{<:Any,T}, cell_coordinates::AbstractVector{<:Vec{sdim, T}}) where {T, sdim} = spatial_coordinate(ip, ξ, cell_coordinates)
+spatial_coordinate(ip::VectorizedInterpolation, ξ::Vec, cell_coordinates::AbstractVector{<:Vec}) = spatial_coordinate(ip, ξ, cell_coordinates)
 
-function spatial_coordinate(interpolation::ScalarInterpolation, ξ::Vec{<:Any,T}, cell_coordinates::AbstractVector{<:Vec{sdim, T}}) where {T, sdim}
+function spatial_coordinate(interpolation::ScalarInterpolation, ξ::Vec, cell_coordinates::AbstractVector{<:Vec{sdim, Tx}}) where {sdim, Tx}
     n_basefuncs = getnbasefunctions(interpolation)
     @boundscheck checkbounds(cell_coordinates, Base.OneTo(n_basefuncs))
 
-    x = zero(Vec{sdim, T})
+    x = zero(Vec{sdim, promote_type(eltype(ξ), Tx)})
     @inbounds for j in 1:n_basefuncs
         M = shape_value(interpolation, ξ, j)
         x += M * cell_coordinates[j]
diff --git a/src/Ferrite.jl b/src/Ferrite.jl
index 013421f21f..e99388a324 100644
--- a/src/Ferrite.jl
+++ b/src/Ferrite.jl
@@ -9,6 +9,7 @@ using StaticArrays
 using Base: @propagate_inbounds
 using NearestNeighbors
 using EnumX
+import ForwardDiff
 
 include("exports.jl")
 
diff --git a/src/PointEvalHandler.jl b/src/PointEvalHandler.jl
index a59c9939bd..d29779884e 100644
--- a/src/PointEvalHandler.jl
+++ b/src/PointEvalHandler.jl
@@ -109,6 +109,24 @@ cellcenter(::Type{<:RefSimplex{dim}}, _::Type{T}) where {dim, T} = Vec{dim, T}((
 _solve_helper(A::Tensor{2,dim}, b::Vec{dim}) where {dim} = inv(A) ⋅ b
 _solve_helper(A::SMatrix{idim, odim}, b::Vec{idim,T}) where {odim, idim, T} = Vec{odim,T}(pinv(A) * b)
 
+
+# Temporary solution until we have MixedTensors
+function _spatial_coordinate(interpolation::ScalarInterpolation, ξ::SVector, cell_coordinates::AbstractVector{<:Vec{sdim, Tx}}) where {sdim, Tx}
+    n_basefuncs = getnbasefunctions(interpolation)
+    @boundscheck checkbounds(cell_coordinates, Base.OneTo(n_basefuncs))
+
+    x = zero(SVector{sdim, promote_type(eltype(ξ), Tx)})
+    @inbounds for j in 1:n_basefuncs
+        M = shape_value(interpolation, Vec(ξ.data), j)
+        x += M * SVector(cell_coordinates[j].data)
+    end
+    return x
+end
+tensor_or_svec_gradient_value(f::F, ξ::Vec{dim}, ::Type{<:Vec{dim}}) where {F, dim} = gradient(f, ξ, :all)
+function tensor_or_svec_gradient_value(f::F, ξ::Vec{rdim}, ::Type{<:Vec{sdim}}) where {F, rdim, sdim}
+    return ForwardDiff.jacobian(f, SVector(ξ.data)), f(ξ)
+end
+
 # find_local_coordinate(interpolation::IP, cell_coordinates::Vector{V}, global_coordinate::V, warn::Bool, linesearch_max_substeps::Int , max_iters::Int, tol_norm::T) where {dim, T, V <: Vec{dim, T}, ref_shape, IP <: Interpolation{ref_shape}}
 #
 # 
@@ -122,11 +140,12 @@ function find_local_coordinate(interpolation::IP, cell_coordinates::Vector{V}, g
     for iter in 1:max_iters
         # Check if still inside element
         check_isoparametric_boundaries(ref_shape, local_guess, sqrt(tol_norm)) || break
-        # Setup J(ξ) and x(ξ)
-        mapping, global_guess = calculate_mapping_and_spatial_coordinate(interpolation, local_guess, cell_coordinates)
-        J = getjacobian(mapping)
+        # Setup J(ξ) and r(ξ) (special handling until MixedTensors)
+        rf(ξ::Vec) = spatial_coordinate(interpolation, ξ, cell_coordinates) - global_coordinate
+        rf(ξ::SVector) = _spatial_coordinate(interpolation, ξ, cell_coordinates) - SVector(global_coordinate.data)
+        J, residual = tensor_or_svec_gradient_value(rf, local_guess, V)
+        
         # Check if converged
-        residual = global_guess - global_coordinate
         best_residual_norm = norm(residual) # for line search below
         if best_residual_norm ≤ tol_norm
             converged = true