diff --git a/src/FEValues/cell_values.jl b/src/FEValues/cell_values.jl
index 2c9c55d8bb..41288a6d1b 100644
--- a/src/FEValues/cell_values.jl
+++ b/src/FEValues/cell_values.jl
@@ -1,7 +1,7 @@
 # Defines CellScalarValues and CellVectorValues and common methods
 """
-    CellScalarValues([::Type{T}], quad_rule::QuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation])
-    CellVectorValues([::Type{T}], quad_rule::QuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation])
+    CellScalarValues([::Type{T}], quad_rule::QuadratureRule, func_interp::Interpolation, [geo_interp::Interpolation])
+    CellVectorValues([::Type{T}], quad_rule::QuadratureRule, func_interp::Interpolation, [geo_interp::Interpolation])
 
 A `CellValues` object facilitates the process of evaluating values of shape functions, gradients of shape functions,
 values of nodal functions, gradients and divergences of nodal functions etc. in the finite element cell. There are
@@ -12,8 +12,8 @@ utilizes scalar shape functions and `CellVectorValues` utilizes vectorial shape
 **Arguments:**
 * `T`: an optional argument (default to `Float64`) to determine the type the internal data is stored as.
 * `quad_rule`: an instance of a [`QuadratureRule`](@ref)
-* `func_interpol`: an instance of an [`Interpolation`](@ref) used to interpolate the approximated function
-* `geom_interpol`: an optional instance of a [`Interpolation`](@ref) which is used to interpolate the geometry
+* `func_interp`: an instance of an [`Interpolation`](@ref) used to interpolate the approximated function
+* `geo_interp`: an optional instance of a [`Interpolation`](@ref) which is used to interpolate the geometry
 
 **Common methods:**
 * [`reinit!`](@ref)
@@ -34,7 +34,7 @@ utilizes scalar shape functions and `CellVectorValues` utilizes vectorial shape
 CellValues
 
 # CellScalarValues
-struct CellScalarValues{dim,T<:Real,refshape<:AbstractRefShape} <: CellValues{dim,T,refshape}
+struct CellScalarValues{dim,T<:Real,refshape<:AbstractRefShape,FI,GI} <: CellValues{dim,T,refshape,FI,GI}
     N::Matrix{T}
     dNdx::Matrix{Vec{dim,T}}
     dNdξ::Matrix{Vec{dim,T}}
@@ -42,47 +42,49 @@ struct CellScalarValues{dim,T<:Real,refshape<:AbstractRefShape} <: CellValues{di
     M::Matrix{T}
     dMdξ::Matrix{Vec{dim,T}}
     qr_weights::Vector{T}
+    func_interp::FI
+    geo_interp::GI
 end
 
-function CellScalarValues(quad_rule::QuadratureRule, func_interpol::Interpolation,
-        geom_interpol::Interpolation=func_interpol)
-    CellScalarValues(Float64, quad_rule, func_interpol, geom_interpol)
+function CellScalarValues(quad_rule::QuadratureRule, func_interp::Interpolation,
+        geo_interp::Interpolation=func_interp)
+    CellScalarValues(Float64, quad_rule, func_interp, geo_interp)
 end
 
-function CellScalarValues(::Type{T}, quad_rule::QuadratureRule{dim,shape}, func_interpol::Interpolation,
-        geom_interpol::Interpolation=func_interpol) where {dim,T,shape<:AbstractRefShape}
+function CellScalarValues(::Type{T}, quad_rule::QuadratureRule{dim,shape,}, func_interp::Interpolation,
+        geo_interp::Interpolation=func_interp) where {dim,T,shape<:AbstractRefShape}
 
-    @assert getdim(func_interpol) == getdim(geom_interpol)
-    @assert getrefshape(func_interpol) == getrefshape(geom_interpol) == shape
+    @assert getdim(func_interp) == getdim(geo_interp)
+    @assert getrefshape(func_interp) == getrefshape(geo_interp) == shape
     n_qpoints = length(getweights(quad_rule))
 
     # Function interpolation
-    n_func_basefuncs = getnbasefunctions(func_interpol)
+    n_func_basefuncs = getnbasefunctions(func_interp)
     N    = fill(zero(T)          * T(NaN), n_func_basefuncs, n_qpoints)
     dNdx = fill(zero(Vec{dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
     dNdξ = fill(zero(Vec{dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
 
     # Geometry interpolation
-    n_geom_basefuncs = getnbasefunctions(geom_interpol)
+    n_geom_basefuncs = getnbasefunctions(geo_interp)
     M    = fill(zero(T)          * T(NaN), n_geom_basefuncs, n_qpoints)
     dMdξ = fill(zero(Vec{dim,T}) * T(NaN), n_geom_basefuncs, n_qpoints)
 
     for (qp, ξ) in enumerate(quad_rule.points)
         for i in 1:n_func_basefuncs
-            dNdξ[i, qp], N[i, qp] = gradient(ξ -> value(func_interpol, i, ξ), ξ, :all)
+            dNdξ[i, qp], N[i, qp] = gradient(ξ -> value(func_interp, i, ξ), ξ, :all)
         end
         for i in 1:n_geom_basefuncs
-            dMdξ[i, qp], M[i, qp] = gradient(ξ -> value(geom_interpol, i, ξ), ξ, :all)
+            dMdξ[i, qp], M[i, qp] = gradient(ξ -> value(geo_interp, i, ξ), ξ, :all)
         end
     end
 
     detJdV = fill(T(NaN), n_qpoints)
 
-    CellScalarValues{dim,T,shape}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule.weights)
+    CellScalarValues{dim,T,shape,typeof(func_interp),typeof(geo_interp)}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule.weights, func_interp, geo_interp)
 end
 
 # CellVectorValues
-struct CellVectorValues{dim,T<:Real,refshape<:AbstractRefShape,M} <: CellValues{dim,T,refshape}
+struct CellVectorValues{dim,T<:Real,refshape<:AbstractRefShape,FI,GI,M} <: CellValues{dim,T,refshape,FI,GI}
     N::Matrix{Vec{dim,T}}
     dNdx::Matrix{Tensor{2,dim,T,M}}
     dNdξ::Matrix{Tensor{2,dim,T,M}}
@@ -90,34 +92,36 @@ struct CellVectorValues{dim,T<:Real,refshape<:AbstractRefShape,M} <: CellValues{
     M::Matrix{T}
     dMdξ::Matrix{Vec{dim,T}}
     qr_weights::Vector{T}
+    func_interp::FI
+    geo_interp::GI
 end
 
-function CellVectorValues(quad_rule::QuadratureRule, func_interpol::Interpolation, geom_interpol::Interpolation=func_interpol)
-    CellVectorValues(Float64, quad_rule, func_interpol, geom_interpol)
+function CellVectorValues(quad_rule::QuadratureRule, func_interp::Interpolation, geo_interp::Interpolation=func_interp)
+    CellVectorValues(Float64, quad_rule, func_interp, geo_interp)
 end
 
-function CellVectorValues(::Type{T}, quad_rule::QuadratureRule{dim,shape}, func_interpol::Interpolation,
-        geom_interpol::Interpolation=func_interpol) where {dim,T,shape<:AbstractRefShape}
+function CellVectorValues(::Type{T}, quad_rule::QuadratureRule{dim,shape}, func_interp::Interpolation,
+        geo_interp::Interpolation=func_interp) where {dim,T,shape<:AbstractRefShape}
 
-    @assert getdim(func_interpol) == getdim(geom_interpol)
-    @assert getrefshape(func_interpol) == getrefshape(geom_interpol) == shape
+    @assert getdim(func_interp) == getdim(geo_interp)
+    @assert getrefshape(func_interp) == getrefshape(geo_interp) == shape
     n_qpoints = length(getweights(quad_rule))
 
     # Function interpolation
-    n_func_basefuncs = getnbasefunctions(func_interpol) * dim
+    n_func_basefuncs = getnbasefunctions(func_interp) * dim
     N    = fill(zero(Vec{dim,T})      * T(NaN), n_func_basefuncs, n_qpoints)
     dNdx = fill(zero(Tensor{2,dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
     dNdξ = fill(zero(Tensor{2,dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
 
     # Geometry interpolation
-    n_geom_basefuncs = getnbasefunctions(geom_interpol)
+    n_geom_basefuncs = getnbasefunctions(geo_interp)
     M    = fill(zero(T)          * T(NaN), n_geom_basefuncs, n_qpoints)
     dMdξ = fill(zero(Vec{dim,T}) * T(NaN), n_geom_basefuncs, n_qpoints)
 
     for (qp, ξ) in enumerate(quad_rule.points)
         basefunc_count = 1
-        for basefunc in 1:getnbasefunctions(func_interpol)
-            dNdξ_temp, N_temp = gradient(ξ -> value(func_interpol, basefunc, ξ), ξ, :all)
+        for basefunc in 1:getnbasefunctions(func_interp)
+            dNdξ_temp, N_temp = gradient(ξ -> value(func_interp, basefunc, ξ), ξ, :all)
             for comp in 1:dim
                 N_comp = zeros(T, dim)
                 N_comp[comp] = N_temp
@@ -130,14 +134,14 @@ function CellVectorValues(::Type{T}, quad_rule::QuadratureRule{dim,shape}, func_
             end
         end
         for basefunc in 1:n_geom_basefuncs
-            dMdξ[basefunc, qp], M[basefunc, qp] = gradient(ξ -> value(geom_interpol, basefunc, ξ), ξ, :all)
+            dMdξ[basefunc, qp], M[basefunc, qp] = gradient(ξ -> value(geo_interp, basefunc, ξ), ξ, :all)
         end
     end
 
     detJdV = fill(T(NaN), n_qpoints)
     MM = Tensors.n_components(Tensors.get_base(eltype(dNdx)))
 
-    CellVectorValues{dim,T,shape,MM}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule.weights)
+    CellVectorValues{dim,T,shape,typeof(func_interp),typeof(geo_interp),MM}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule.weights, func_interp, geo_interp)
 end
 
 function reinit!(cv::CellValues{dim}, x::AbstractVector{Vec{dim,T}}) where {dim,T}
@@ -154,7 +158,7 @@ function reinit!(cv::CellValues{dim}, x::AbstractVector{Vec{dim,T}}) where {dim,
             fecv_J += x[j] ⊗ cv.dMdξ[j, i]
         end
         detJ = det(fecv_J)
-        detJ > 0.0 || throw(ArgumentError("det(J) is not positive: det(J) = $(detJ)"))
+        detJ > 0.0 || throw_detJ_not_pos(detJ)
         cv.detJdV[i] = detJ * w
         Jinv = inv(fecv_J)
         for j in 1:n_func_basefuncs
diff --git a/src/FEValues/common_values.jl b/src/FEValues/common_values.jl
index d69328485f..cdac8f45d7 100644
--- a/src/FEValues/common_values.jl
+++ b/src/FEValues/common_values.jl
@@ -2,14 +2,16 @@
 
 using Base: @propagate_inbounds
 
-const ScalarValues{dim,T,shape} = Union{CellScalarValues{dim,T,shape},FaceScalarValues{dim,T,shape}}
-const VectorValues{dim,T,shape} = Union{CellVectorValues{dim,T,shape},FaceVectorValues{dim,T,shape}}
+@noinline throw_detJ_not_pos(detJ) = throw(ArgumentError("det(J) is not positive: det(J) = $(detJ)"))
 
-getnbasefunctions(cv::Values) = size(cv.N, 1)
-getngeobasefunctions(cv::Values) = size(cv.M, 1)
+const ScalarValues{dim,T,shape,func_interp,geo_interp} = Union{CellScalarValues{dim,T,shape},FaceScalarValues{dim,T,shape,func_interp,geo_interp}}
+const VectorValues{dim,T,shape,func_interp,geo_interp} = Union{CellVectorValues{dim,T,shape},FaceVectorValues{dim,T,shape,func_interp,geo_interp}}
 
-getn_scalarbasefunctions(cv::ScalarValues) = size(cv.N, 1)
-getn_scalarbasefunctions(cv::VectorValues{dim}) where {dim} = size(cv.N, 1) ÷ dim
+getnbasefunctions(fe::Values{dim}) where {dim} = getnbasefunctions(fe.func_interp) * (fe isa VectorValues ? dim : 1)
+getngeobasefunctions(fe::Values) = getnbasefunctions(fe.geo_interp)
+
+getn_scalarbasefunctions(cv::ScalarValues) = getnbasefunctions(cv)
+getn_scalarbasefunctions(cv::VectorValues{dim}) where {dim} = getnbasefunctions(cv) ÷ dim
 
 """
     reinit!(cv::CellValues, x::Vector)
@@ -104,8 +106,7 @@ where ``u_i`` are the value of ``u`` in the nodes. For a vector valued function
 nodal values of ``\\mathbf{u}``.
 """
 function function_value(fe_v::Values{dim}, q_point::Int, u::AbstractVector{T}, dof_range = eachindex(u)) where {dim,T}
-    n_base_funcs = getn_scalarbasefunctions(fe_v)
-    isa(fe_v, VectorValues) && (n_base_funcs *= dim)
+    n_base_funcs = getnbasefunctions(fe_v)
     @assert length(dof_range) == n_base_funcs
     @boundscheck checkbounds(u, dof_range)
     val = zero(_valuetype(fe_v, u))
@@ -150,8 +151,7 @@ For a vector valued function the gradient is computed as
 where ``\\mathbf{u}_i`` are the nodal values of ``\\mathbf{u}``.
 """
 function function_gradient(fe_v::Values{dim}, q_point::Int, u::AbstractVector{T}, dof_range = eachindex(u)) where {dim,T}
-    n_base_funcs = getn_scalarbasefunctions(fe_v)
-    isa(fe_v, VectorValues) && (n_base_funcs *= dim)
+    n_base_funcs = getnbasefunctions(fe_v)
     @assert length(dof_range) == n_base_funcs
     @boundscheck checkbounds(u, dof_range)
     grad = zero(_gradienttype(fe_v, u))
diff --git a/src/FEValues/face_values.jl b/src/FEValues/face_values.jl
index 428494c601..5c60726baa 100644
--- a/src/FEValues/face_values.jl
+++ b/src/FEValues/face_values.jl
@@ -1,7 +1,7 @@
 # Defines FaceScalarValues and FaceVectorValues and common methods
 """
-    FaceScalarValues([::Type{T}], quad_rule::QuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation])
-    FaceVectorValues([::Type{T}], quad_rule::QuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation])
+    FaceScalarValues([::Type{T}], quad_rule::QuadratureRule, func_interp::Interpolation, [geo_interp::Interpolation])
+    FaceVectorValues([::Type{T}], quad_rule::QuadratureRule, func_interp::Interpolation, [geo_interp::Interpolation])
 
 A `FaceValues` object facilitates the process of evaluating values of shape functions, gradients of shape functions,
 values of nodal functions, gradients and divergences of nodal functions etc. on the faces of finite elements. There are
@@ -17,8 +17,8 @@ For a scalar field, the `FaceScalarValues` type should be used. For vector field
 
 * `T`: an optional argument to determine the type the internal data is stored as.
 * `quad_rule`: an instance of a [`QuadratureRule`](@ref)
-* `func_interpol`: an instance of an [`Interpolation`](@ref) used to interpolate the approximated function
-* `geom_interpol`: an optional instance of an [`Interpolation`](@ref) which is used to interpolate the geometry
+* `func_interp`: an instance of an [`Interpolation`](@ref) used to interpolate the approximated function
+* `geo_interp`: an optional instance of an [`Interpolation`](@ref) which is used to interpolate the geometry
 
 **Common methods:**
 
@@ -40,7 +40,7 @@ For a scalar field, the `FaceScalarValues` type should be used. For vector field
 FaceValues
 
 # FaceScalarValues
-struct FaceScalarValues{dim,T<:Real,refshape<:AbstractRefShape} <: FaceValues{dim,T,refshape}
+struct FaceScalarValues{dim,T<:Real,refshape<:AbstractRefShape,FI,GI} <: FaceValues{dim,T,refshape,FI,GI}
     N::Array{T,3}
     dNdx::Array{Vec{dim,T},3}
     dNdξ::Array{Vec{dim,T},3}
@@ -50,54 +50,56 @@ struct FaceScalarValues{dim,T<:Real,refshape<:AbstractRefShape} <: FaceValues{di
     dMdξ::Array{Vec{dim,T},3}
     qr_weights::Vector{T}
     current_face::ScalarWrapper{Int}
+    func_interp::FI
+    geo_interp::GI
 end
 
-function FaceScalarValues(quad_rule::QuadratureRule, func_interpol::Interpolation,
-                          geom_interpol::Interpolation=func_interpol)
-    FaceScalarValues(Float64, quad_rule, func_interpol, geom_interpol)
+function FaceScalarValues(quad_rule::QuadratureRule, func_interp::Interpolation,
+                          geo_interp::Interpolation=func_interp)
+    FaceScalarValues(Float64, quad_rule, func_interp, geo_interp)
 end
 
-function FaceScalarValues(::Type{T}, quad_rule::QuadratureRule{dim_qr,shape}, func_interpol::Interpolation,
-        geom_interpol::Interpolation=func_interpol) where {dim_qr,T,shape<:AbstractRefShape}
+function FaceScalarValues(::Type{T}, quad_rule::QuadratureRule{dim_qr,shape}, func_interp::Interpolation,
+        geo_interp::Interpolation=func_interp) where {dim_qr,T,shape<:AbstractRefShape}
 
-    @assert getdim(func_interpol) == getdim(geom_interpol)
-    @assert getrefshape(func_interpol) == getrefshape(geom_interpol) == shape
+    @assert getdim(func_interp) == getdim(geo_interp)
+    @assert getrefshape(func_interp) == getrefshape(geo_interp) == shape
     n_qpoints = length(getweights(quad_rule))
     dim = dim_qr + 1
 
-    face_quad_rule = create_face_quad_rule(quad_rule, func_interpol)
+    face_quad_rule = create_face_quad_rule(quad_rule, func_interp)
     n_faces = length(face_quad_rule)
 
     # Normals
     normals = zeros(Vec{dim,T}, n_qpoints)
 
     # Function interpolation
-    n_func_basefuncs = getnbasefunctions(func_interpol)
+    n_func_basefuncs = getnbasefunctions(func_interp)
     N =    fill(zero(T)          * T(NaN), n_func_basefuncs, n_qpoints, n_faces)
     dNdx = fill(zero(Vec{dim,T}) * T(NaN), n_func_basefuncs, n_qpoints, n_faces)
     dNdξ = fill(zero(Vec{dim,T}) * T(NaN), n_func_basefuncs, n_qpoints, n_faces)
 
     # Geometry interpolation
-    n_geom_basefuncs = getnbasefunctions(geom_interpol)
+    n_geom_basefuncs = getnbasefunctions(geo_interp)
     M =    fill(zero(T)          * T(NaN), n_geom_basefuncs, n_qpoints, n_faces)
     dMdξ = fill(zero(Vec{dim,T}) * T(NaN), n_geom_basefuncs, n_qpoints, n_faces)
 
     for face in 1:n_faces, (qp, ξ) in enumerate(face_quad_rule[face].points)
         for i in 1:n_func_basefuncs
-            dNdξ[i, qp, face], N[i, qp, face] = gradient(ξ -> value(func_interpol, i, ξ), ξ, :all)
+            dNdξ[i, qp, face], N[i, qp, face] = gradient(ξ -> value(func_interp, i, ξ), ξ, :all)
         end
         for i in 1:n_geom_basefuncs
-            dMdξ[i, qp, face], M[i, qp, face] = gradient(ξ -> value(geom_interpol, i, ξ), ξ, :all)
+            dMdξ[i, qp, face], M[i, qp, face] = gradient(ξ -> value(geo_interp, i, ξ), ξ, :all)
         end
     end
 
     detJdV = fill(T(NaN), n_qpoints, n_faces)
 
-    FaceScalarValues{dim,T,shape}(N, dNdx, dNdξ, detJdV, normals, M, dMdξ, quad_rule.weights, ScalarWrapper(0))
+    FaceScalarValues{dim,T,shape,typeof(func_interp),typeof(geo_interp)}(N, dNdx, dNdξ, detJdV, normals, M, dMdξ, quad_rule.weights, ScalarWrapper(0), func_interp, geo_interp)
 end
 
 # FaceVectorValues
-struct FaceVectorValues{dim,T<:Real,refshape<:AbstractRefShape,M} <: FaceValues{dim,T,refshape}
+struct FaceVectorValues{dim,T<:Real,refshape<:AbstractRefShape,FI,GI,M} <: FaceValues{dim,T,refshape,FI,GI}
     N::Array{Vec{dim,T},3}
     dNdx::Array{Tensor{2,dim,T,M},3}
     dNdξ::Array{Tensor{2,dim,T,M},3}
@@ -107,41 +109,43 @@ struct FaceVectorValues{dim,T<:Real,refshape<:AbstractRefShape,M} <: FaceValues{
     dMdξ::Array{Vec{dim,T},3}
     qr_weights::Vector{T}
     current_face::ScalarWrapper{Int}
+    func_interp::FI
+    geo_interp::GI
 end
 
-function FaceVectorValues(quad_rule::QuadratureRule, func_interpol::Interpolation, geom_interpol::Interpolation=func_interpol)
-    FaceVectorValues(Float64, quad_rule, func_interpol, geom_interpol)
+function FaceVectorValues(quad_rule::QuadratureRule, func_interp::Interpolation, geo_interp::Interpolation=func_interp)
+    FaceVectorValues(Float64, quad_rule, func_interp, geo_interp)
 end
 
-function FaceVectorValues(::Type{T}, quad_rule::QuadratureRule{dim_qr,shape}, func_interpol::Interpolation,
-        geom_interpol::Interpolation=func_interpol) where {dim_qr,T,shape<:AbstractRefShape}
+function FaceVectorValues(::Type{T}, quad_rule::QuadratureRule{dim_qr,shape}, func_interp::Interpolation,
+        geo_interp::Interpolation=func_interp) where {dim_qr,T,shape<:AbstractRefShape}
 
-    @assert getdim(func_interpol) == getdim(geom_interpol)
-    @assert getrefshape(func_interpol) == getrefshape(geom_interpol) == shape
+    @assert getdim(func_interp) == getdim(geo_interp)
+    @assert getrefshape(func_interp) == getrefshape(geo_interp) == shape
     n_qpoints = length(getweights(quad_rule))
     dim = dim_qr + 1
 
-    face_quad_rule = create_face_quad_rule(quad_rule, func_interpol)
+    face_quad_rule = create_face_quad_rule(quad_rule, func_interp)
     n_faces = length(face_quad_rule)
 
     # Normals
     normals = zeros(Vec{dim,T}, n_qpoints)
 
     # Function interpolation
-    n_func_basefuncs = getnbasefunctions(func_interpol) * dim
+    n_func_basefuncs = getnbasefunctions(func_interp) * dim
     N    = fill(zero(Vec{dim,T})      * T(NaN), n_func_basefuncs, n_qpoints, n_faces)
     dNdx = fill(zero(Tensor{2,dim,T}) * T(NaN), n_func_basefuncs, n_qpoints, n_faces)
     dNdξ = fill(zero(Tensor{2,dim,T}) * T(NaN), n_func_basefuncs, n_qpoints, n_faces)
 
     # Geometry interpolation
-    n_geom_basefuncs = getnbasefunctions(geom_interpol)
+    n_geom_basefuncs = getnbasefunctions(geo_interp)
     M    = fill(zero(T)          * T(NaN), n_geom_basefuncs, n_qpoints, n_faces)
     dMdξ = fill(zero(Vec{dim,T}) * T(NaN), n_geom_basefuncs, n_qpoints, n_faces)
 
     for face in 1:n_faces, (qp, ξ) in enumerate(face_quad_rule[face].points)
         basefunc_count = 1
-        for basefunc in 1:getnbasefunctions(func_interpol)
-            dNdξ_temp, N_temp = gradient(ξ -> value(func_interpol, basefunc, ξ), ξ, :all)
+        for basefunc in 1:getnbasefunctions(func_interp)
+            dNdξ_temp, N_temp = gradient(ξ -> value(func_interp, basefunc, ξ), ξ, :all)
             for comp in 1:dim
                 N_comp = zeros(T, dim)
                 N_comp[comp] = N_temp
@@ -154,21 +158,20 @@ function FaceVectorValues(::Type{T}, quad_rule::QuadratureRule{dim_qr,shape}, fu
             end
         end
         for basefunc in 1:n_geom_basefuncs
-            dMdξ[basefunc, qp, face], M[basefunc, qp, face] = gradient(ξ -> value(geom_interpol, basefunc, ξ), ξ, :all)
+            dMdξ[basefunc, qp, face], M[basefunc, qp, face] = gradient(ξ -> value(geo_interp, basefunc, ξ), ξ, :all)
         end
     end
 
     detJdV = fill(T(NaN), n_qpoints, n_faces)
     MM = Tensors.n_components(Tensors.get_base(eltype(dNdx)))
 
-    FaceVectorValues{dim,T,shape,MM}(N, dNdx, dNdξ, detJdV, normals, M, dMdξ, quad_rule.weights, ScalarWrapper(0))
+    FaceVectorValues{dim,T,shape,typeof(func_interp),typeof(geo_interp),MM}(N, dNdx, dNdξ, detJdV, normals, M, dMdξ, quad_rule.weights, ScalarWrapper(0), func_interp, geo_interp)
 end
 
 function reinit!(fv::FaceValues{dim}, x::AbstractVector{Vec{dim,T}}, face::Int) where {dim,T}
     n_geom_basefuncs = getngeobasefunctions(fv)
-    n_func_basefuncs = getn_scalarbasefunctions(fv)
+    n_func_basefuncs = getnbasefunctions(fv)
     @assert length(x) == n_geom_basefuncs
-    isa(fv, FaceVectorValues) && (n_func_basefuncs *= dim)
 
     fv.current_face[] = face
     cb = getcurrentface(fv)
@@ -183,7 +186,7 @@ function reinit!(fv::FaceValues{dim}, x::AbstractVector{Vec{dim,T}}, face::Int)
         fv.normals[i] = weight_norm / norm(weight_norm)
         detJ = norm(weight_norm)
 
-        detJ > 0.0 || throw(ArgumentError("det(J) is not positive: det(J) = $(detJ)"))
+        detJ > 0.0 || throw_detJ_not_pos(detJ)
         fv.detJdV[i, cb] = detJ * w
         Jinv = inv(fefv_J)
         for j in 1:n_func_basefuncs
@@ -215,16 +218,16 @@ struct BCValues{T}
     current_face::ScalarWrapper{Int}
 end
 
-BCValues(func_interpol::Interpolation, geom_interpol::Interpolation) =
-    BCValues(Float64, func_interpol, geom_interpol)
+BCValues(func_interp::Interpolation, geo_interp::Interpolation) =
+    BCValues(Float64, func_interp, geo_interp)
 
-function BCValues(::Type{T}, func_interpol::Interpolation{dim,refshape}, geom_interpol::Interpolation{dim,refshape}) where {T,dim,refshape}
+function BCValues(::Type{T}, func_interp::Interpolation{dim,refshape}, geo_interp::Interpolation{dim,refshape}) where {T,dim,refshape}
     # set up quadrature rules for each face with dof-positions
-    # (determined by func_interpol) as the quadrature points
-    interpolation_coords = reference_coordinates(func_interpol)
+    # (determined by func_interp) as the quadrature points
+    interpolation_coords = reference_coordinates(func_interp)
 
     qrs = QuadratureRule{dim,refshape,T}[]
-    for face in faces(func_interpol)
+    for face in faces(func_interp)
         dofcoords = Vec{dim,T}[]
         for facedof in face
             push!(dofcoords, interpolation_coords[facedof])
@@ -235,12 +238,12 @@ function BCValues(::Type{T}, func_interpol::Interpolation{dim,refshape}, geom_in
 
     n_qpoints = length(getweights(qrs[1])) # assume same in all
     n_faces = length(qrs)
-    n_geom_basefuncs = getnbasefunctions(geom_interpol)
+    n_geom_basefuncs = getnbasefunctions(geo_interp)
     M =    fill(zero(T)           * T(NaN), n_geom_basefuncs, n_qpoints, n_faces)
 
     for face in 1:n_faces, (qp, ξ) in enumerate(qrs[face].points)
         for i in 1:n_geom_basefuncs
-            M[i, qp, face] = value(geom_interpol, i, ξ)
+            M[i, qp, face] = value(geo_interp, i, ξ)
         end
     end
 
diff --git a/src/Ferrite.jl b/src/Ferrite.jl
index e4e67d2715..51e9e81ae3 100644
--- a/src/Ferrite.jl
+++ b/src/Ferrite.jl
@@ -23,9 +23,9 @@ struct RefCube <: AbstractRefShape end
 """
 Abstract type which has `CellValues` and `FaceValues` as subtypes
 """
-abstract type Values{dim,T,refshape} end
-abstract type CellValues{dim,T,refshape} <: Values{dim,T,refshape} end
-abstract type FaceValues{dim,T,refshape} <: Values{dim,T,refshape} end
+abstract type Values{dim,T,refshape,FI,GI} end
+abstract type CellValues{dim,T,refshape,FI,GI} <: Values{dim,T,refshape,FI,GI} end
+abstract type FaceValues{dim,T,refshape,FI,GI} <: Values{dim,T,refshape,FI,GI} end
 
 include("utils.jl")
 
diff --git a/src/interpolations.jl b/src/interpolations.jl
index 9fdab87371..bcdc450c09 100644
--- a/src/interpolations.jl
+++ b/src/interpolations.jl
@@ -34,17 +34,17 @@ abstract type Interpolation{dim,shape,order} end
 """
 Return the dimension of an `Interpolation`
 """
-@inline getdim(ip::Interpolation{dim}) where {dim} = dim
+@inline getdim(::Interpolation{dim}) where {dim} = dim
 
 """
 Return the reference shape of an `Interpolation`
 """
-@inline getrefshape(ip::Interpolation{dim,shape}) where {dim,shape} = shape
+@inline getrefshape(::Interpolation{dim,shape}) where {dim,shape} = shape
 
 """
 Return the polynomial order of the `Interpolation`
 """
-@inline getorder(ip::Interpolation{dim,shape,order}) where {dim,shape,order} = order
+@inline getorder(::Interpolation{dim,shape,order}) where {dim,shape,order} = order
 
 """
 Compute the value of the shape functions at a point ξ for a given interpolation
@@ -60,6 +60,8 @@ function derivative(ip::Interpolation{dim}, ξ::Vec{dim,T}) where {dim,T}
     [gradient(ξ -> value(ip, i, ξ), ξ) for i in 1:getnbasefunctions(ip)]
 end
 
+Base.copy(ip::Interpolation) = ip
+
 #####################
 # Utility functions #
 #####################
diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl
index 3c2ea0c44d..cee4102cf0 100644
--- a/test/test_cellvalues.jl
+++ b/test/test_cellvalues.jl
@@ -1,5 +1,5 @@
 @testset "CellValues" begin
-for (func_interpol, quad_rule) in  (
+for (func_interp, quad_rule) in  (
                                     (Lagrange{1, RefCube, 1}(), QuadratureRule{1, RefCube}(2)),
                                     (Lagrange{1, RefCube, 2}(), QuadratureRule{1, RefCube}(2)),
                                     (Lagrange{2, RefCube, 1}(), QuadratureRule{2, RefCube}(2)),
@@ -13,14 +13,14 @@ for (func_interpol, quad_rule) in  (
                                    )
 
     for fe_valtype in (CellScalarValues, CellVectorValues)
-        cv = fe_valtype(quad_rule, func_interpol)
-        ndim = Ferrite.getdim(func_interpol)
-        n_basefuncs = getnbasefunctions(func_interpol)
+        cv = fe_valtype(quad_rule, func_interp)
+        ndim = Ferrite.getdim(func_interp)
+        n_basefuncs = getnbasefunctions(func_interp)
 
         fe_valtype == CellScalarValues && @test getnbasefunctions(cv) == n_basefuncs
-        fe_valtype == CellVectorValues && @test getnbasefunctions(cv) == n_basefuncs * Ferrite.getdim(func_interpol)
+        fe_valtype == CellVectorValues && @test getnbasefunctions(cv) == n_basefuncs * Ferrite.getdim(func_interp)
 
-        x, n = valid_coordinates_and_normals(func_interpol)
+        x, n = valid_coordinates_and_normals(func_interp)
         reinit!(cv, x)
 
         # We test this by applying a given deformation gradient on all the nodes.
@@ -58,16 +58,16 @@ for (func_interpol, quad_rule) in  (
         for i in 1:getnquadpoints(cv)
             vol += getdetJdV(cv,i)
         end
-        @test vol ≈ calculate_volume(func_interpol, x)
+        @test vol ≈ calculate_volume(func_interp, x)
 
         # Test quadrature rule after reinit! with ref. coords
-        x = Ferrite.reference_coordinates(func_interpol)
+        x = Ferrite.reference_coordinates(func_interp)
         reinit!(cv, x)
         vol = 0.0
         for i in 1:getnquadpoints(cv)
             vol += getdetJdV(cv,i)
         end
-        @test vol ≈ reference_volume(func_interpol)
+        @test vol ≈ reference_volume(func_interp)
 
         # Test spatial coordinate (after reinit with ref.coords we should get back the quad_points)
         for (i, qp_x) in enumerate(getpoints(quad_rule))
@@ -78,8 +78,12 @@ for (func_interpol, quad_rule) in  (
         cvc = copy(cv)
         for fname in fieldnames(typeof(cv))
             @test typeof(cv) == typeof(cvc)
-            @test pointer(getfield(cv, fname)) != pointer(getfield(cvc, fname))
-            @test getfield(cv, fname) == getfield(cvc, fname)
+            v = getfield(cv, fname)
+            vc = getfield(cvc, fname)
+            if hasmethod(pointer, Tuple{typeof(v)})
+                @test pointer(v) != pointer(vc)
+            end
+            @test v == vc
         end
     end
 end
diff --git a/test/test_facevalues.jl b/test/test_facevalues.jl
index f775f31fbc..20fdb337e7 100644
--- a/test/test_facevalues.jl
+++ b/test/test_facevalues.jl
@@ -1,5 +1,5 @@
 @testset "FaceValues" begin
-for (func_interpol, quad_rule) in  (
+for (func_interp, quad_rule) in  (
                                     (Lagrange{1, RefCube, 1}(), QuadratureRule{0, RefCube}(2)),
                                     (Lagrange{1, RefCube, 2}(), QuadratureRule{0, RefCube}(2)),
                                     (Lagrange{2, RefCube, 1}(), QuadratureRule{1, RefCube}(2)),
@@ -13,15 +13,15 @@ for (func_interpol, quad_rule) in  (
                                    )
 
     for fe_valtype in (FaceScalarValues, FaceVectorValues)
-        fv = fe_valtype(quad_rule, func_interpol)
-        ndim = Ferrite.getdim(func_interpol)
-        n_basefuncs = getnbasefunctions(func_interpol)
+        fv = fe_valtype(quad_rule, func_interp)
+        ndim = Ferrite.getdim(func_interp)
+        n_basefuncs = getnbasefunctions(func_interp)
 
         fe_valtype == FaceScalarValues && @test getnbasefunctions(fv) == n_basefuncs
-        fe_valtype == FaceVectorValues && @test getnbasefunctions(fv) == n_basefuncs * Ferrite.getdim(func_interpol)
+        fe_valtype == FaceVectorValues && @test getnbasefunctions(fv) == n_basefuncs * Ferrite.getdim(func_interp)
 
-        xs, n = valid_coordinates_and_normals(func_interpol)
-        for face in 1:getnfaces(func_interpol)
+        xs, n = valid_coordinates_and_normals(func_interp)
+        for face in 1:getnfaces(func_interp)
             reinit!(fv, xs, face)
             @test Ferrite.getcurrentface(fv) == face
 
@@ -61,17 +61,17 @@ for (func_interpol, quad_rule) in  (
             for i in 1:getnquadpoints(fv)
                 vol += getdetJdV(fv,i)
             end
-            x_face = xs[[Ferrite.faces(func_interpol)[face]...]]
-            @test vol ≈ calculate_volume(Ferrite.getlowerdim(func_interpol), x_face)
+            x_face = xs[[Ferrite.faces(func_interp)[face]...]]
+            @test vol ≈ calculate_volume(Ferrite.getlowerdim(func_interp), x_face)
 
             # Test quadrature rule after reinit! with ref. coords
-            x = Ferrite.reference_coordinates(func_interpol)
+            x = Ferrite.reference_coordinates(func_interp)
             reinit!(fv, x, face)
             vol = 0.0
             for i in 1:getnquadpoints(fv)
                 vol += getdetJdV(fv, i)
             end
-            @test vol ≈ reference_volume(func_interpol, face)
+            @test vol ≈ reference_volume(func_interp, face)
 
             # Test spatial coordinate (after reinit with ref.coords we should get back the quad_points)
             # TODO: Renable somehow after quad rule is no longer stored in FaceValues
@@ -84,11 +84,13 @@ for (func_interpol, quad_rule) in  (
         # test copy
         fvc = copy(fv)
         for fname in fieldnames(typeof(fv))
-            @test typeof(fv) == typeof(fvc)
-            if !isa(getfield(fv, fname), Ferrite.ScalarWrapper)
-                @test pointer(getfield(fv, fname)) != pointer(getfield(fvc, fname))
-                @test getfield(fv, fname) == getfield(fvc, fname)
+            isa(getfield(fv, fname), Ferrite.ScalarWrapper) && continue
+            v = getfield(fv, fname)
+            vc = getfield(fvc, fname)
+            if hasmethod(pointer, Tuple{typeof(v)})
+                @test pointer(v) != pointer(vc)
             end
+            @test v == vc
         end
     end
 end