diff --git a/src/Adaptivity/AdaptivityGlues.jl b/src/Adaptivity/AdaptivityGlues.jl
index 11d830a29..ba417d0ef 100644
--- a/src/Adaptivity/AdaptivityGlues.jl
+++ b/src/Adaptivity/AdaptivityGlues.jl
@@ -235,6 +235,7 @@ function get_d_to_fface_to_cface(glue::AdaptivityGlue{<:RefinementGlue},
   # Local data for each coarse cell, at the RefinementRule level
   rrules = get_old_cell_refinement_rules(glue)
   ccell_to_d_to_faces = lazy_map(rr->map(d->Geometry.get_faces(get_grid_topology(rr.ref_grid),Dc,d),0:Dc),rrules)
+  #ccell_to_d_to_faces = lazy_map(rr->map(d->Geometry.get_faces(get_grid_topology(rr.ref_grid),Dc,d)::Table{Int32, Vector{Int32}, Vector{Int32}},0:Dc),rrules)
   ccell_to_d_to_fface_to_parent_face = lazy_map(get_d_to_face_to_parent_face,rrules)
   fcell_to_child_id = glue.n2o_cell_to_child_id
 
diff --git a/src/Algebra/LinearSolvers.jl b/src/Algebra/LinearSolvers.jl
index 2934f0e24..799d77083 100644
--- a/src/Algebra/LinearSolvers.jl
+++ b/src/Algebra/LinearSolvers.jl
@@ -147,11 +147,11 @@ function test_linear_solver(
   ss = symbolic_setup(ls,A)
   ns = numerical_setup(ss,A)
   numerical_setup!(ns,A)
-  solve!(y,ns,b)
+  @inferred solve!(y,ns,b)
   @test x ≈ y
 
   y .= zero(eltype(y))
-  solve!(y,ls,A,b)
+  @inferred solve!(y,ls,A,b)
   @test x ≈ y
 
   op = AffineOperator(A,b)
diff --git a/src/Algebra/NonlinearSolvers.jl b/src/Algebra/NonlinearSolvers.jl
index 0d01480cf..9363d241f 100644
--- a/src/Algebra/NonlinearSolvers.jl
+++ b/src/Algebra/NonlinearSolvers.jl
@@ -59,15 +59,15 @@ function test_nonlinear_solver(
   pred::Function=isapprox)
 
   x1 = copy(x0)
-  cache = solve!(x1,nls,op)
+  cache = @inferred solve!(x1,nls,op)
   @test pred(x1,x)
 
   x1 = copy(x0)
-  cache = solve!(x1,nls,op,cache)
+  cache = @inferred solve!(x1,nls,op,cache)
   @test pred(x1,x)
 
   x1 = copy(x0)
-  cache = solve!(x1,nls,op,cache)
+  cache = @inferred solve!(x1,nls,op,cache)
   @test pred(x1,x)
 
 end
diff --git a/src/Arrays/Autodiff.jl b/src/Arrays/Autodiff.jl
index b306ab1cd..326d58e19 100644
--- a/src/Arrays/Autodiff.jl
+++ b/src/Arrays/Autodiff.jl
@@ -99,13 +99,15 @@ end
 ConfigMap(f) = ConfigMap(f,nothing)
 
 # TODO Prescribing long chunk size can lead to slow compilation times!
-function return_cache(k::ConfigMap{typeof(ForwardDiff.gradient)},x)
-  cfg = ForwardDiff.GradientConfig(k.tag,x,ForwardDiff.Chunk{length(x)}())
+function return_cache(k::ConfigMap{typeof(ForwardDiff.gradient)}, x)
+  chunk_x = ForwardDiff.Chunk{length(x)}() # TODO fix type instability
+  cfg = ForwardDiff.GradientConfig(k.tag, x, chunk_x)
   return cfg
 end
 
-function return_cache(k::ConfigMap{typeof(ForwardDiff.jacobian)},x)
-  cfg = ForwardDiff.JacobianConfig(k.tag,x,ForwardDiff.Chunk{length(x)}())
+function return_cache(k::ConfigMap{typeof(ForwardDiff.jacobian)}, x)
+  chunk_x = ForwardDiff.Chunk{length(x)}() # TODO fix type instability
+  cfg = ForwardDiff.JacobianConfig(k.tag, x, chunk_x)
   return cfg
 end
 
diff --git a/src/Arrays/KeyToValMaps.jl b/src/Arrays/KeyToValMaps.jl
index b666027ab..212e2f521 100644
--- a/src/Arrays/KeyToValMaps.jl
+++ b/src/Arrays/KeyToValMaps.jl
@@ -13,12 +13,8 @@ function return_cache(m::KeyToValMap,key)
   d = Dict{K,V}()
 end
 
-function evaluate!(cache,m::KeyToValMap,key)
-  if haskey(cache,key)
-    return get(cache,key,nothing)
-  else
-    val = m.key_to_val(key)
-    cache[key] = val
-    return val
-  end
+function evaluate!(cache::Dict{K,V},m::KeyToValMap,key) where {K,V}
+  get!(cache, key) do
+    m.key_to_val(key)
+  end::V
 end
diff --git a/src/CellData/CellStates.jl b/src/CellData/CellStates.jl
index 24c45622f..9b2938d7a 100644
--- a/src/CellData/CellStates.jl
+++ b/src/CellData/CellStates.jl
@@ -1,23 +1,27 @@
 
 """
+    CellState{T,P<:CellPoint,V<:AbstractArray} <: CellField
+
 This can be used as a CellField as long as one evaluates it
 on the stored CellPoint.
 """
-struct CellState{T,P<:CellPoint} <: CellField
+struct CellState{T,P<:CellPoint,V<:AbstractArray} <: CellField
   points::P
-  values::AbstractArray
+  values::V
 
   function CellState{T}(::UndefInitializer,points::CellPoint) where T
     values = _init_values(T,get_data(points))
     P = typeof(points)
-    new{T,P}(points,values)
+    V = typeof(values)
+    new{T,P,V}(points,values)
   end
 
   function CellState(v::Number,points::CellPoint)
     values = _init_values(v,get_data(points))
     T = typeof(v)
     P = typeof(points)
-    new{T,P}(points,values)
+    V = typeof(values)
+    new{T,P,V}(points,values)
   end
 end
 
@@ -50,7 +54,7 @@ function get_data(f::CellState)
 end
 
 get_triangulation(f::CellState) = get_triangulation(f.points)
-DomainStyle(::Type{CellState{T,P}}) where {T,P} = DomainStyle(P)
+DomainStyle(::Type{<:CellState{T,P}}) where {T,P} = DomainStyle(P)
 
 _get_cell_points(a::CellState) = a.points
 
diff --git a/src/FESpaces/CLagrangianFESpaces.jl b/src/FESpaces/CLagrangianFESpaces.jl
index af97a50c6..7c8fbba4f 100644
--- a/src/FESpaces/CLagrangianFESpaces.jl
+++ b/src/FESpaces/CLagrangianFESpaces.jl
@@ -86,7 +86,7 @@ function CLagrangianFESpace(
 end
 
 # Allowed configurations of CLagrangianFESpace constructors
-# We need: 
+# We need:
 #  - H1 conformity
 #  - Both the geometric and FE reffes have to be the same (and linear - this can be relaxed in the future)
 #  - nvertices == nnodes (no periodicity, etc...)
@@ -277,13 +277,13 @@ end
 function _generate_cell_reffe_clagrangian(z,grid)
   ctype_to_reffe = get_reffes(grid)
   cell_to_ctype = get_cell_type(grid)
-  function fun(reffe)
+  function fun(@nospecialize reffe)
     p = get_polytope(reffe)
     orders = get_orders(reffe)
     LagrangianRefFE(typeof(z),p,orders)
   end
-  # using map instead of lazy_map seems to kill the compiler.
-  ctype_to_reffe = collect(lazy_map(fun,ctype_to_reffe))
+
+  ctype_to_reffe = collect(map(fun,ctype_to_reffe))
   cell_to_reffe = expand_cell_data(ctype_to_reffe,cell_to_ctype)
 end
 
diff --git a/src/FESpaces/DiscreteModelWithFEMaps.jl b/src/FESpaces/DiscreteModelWithFEMaps.jl
index e885fdc20..3cde53529 100644
--- a/src/FESpaces/DiscreteModelWithFEMaps.jl
+++ b/src/FESpaces/DiscreteModelWithFEMaps.jl
@@ -33,12 +33,10 @@ function GridWithFEMap(model,orders::AbstractArray; kwargs...)
   ct = get_cell_type(model)
   ps = lazy_map(Reindex(_ps), ct)
 
-  f(a,b) = LagrangianRefFE(T,a,b)
-  reffes = lazy_map(f,ps,os)
+  reffes = map( (a,b) -> LagrangianRefFE(T,a,b), ps,os)
   Vₕ = FESpace(model,reffes;conformity=:H1,kwargs...)
 
-  fs(a,b) = LagrangianRefFE(Ts,a,b)
-  s_reffes = lazy_map(f,ps,os)
+  s_reffes  = map( (a,b) -> LagrangianRefFE(Ts,a,b), ps,os)
   Vₕ_scal = FESpace(model,s_reffes;conformity=:H1)
 
   grid = get_grid(model)
diff --git a/src/FESpaces/Pullbacks.jl b/src/FESpaces/Pullbacks.jl
index 357129804..f0a59363c 100644
--- a/src/FESpaces/Pullbacks.jl
+++ b/src/FESpaces/Pullbacks.jl
@@ -159,7 +159,8 @@ function return_cache(k::NormalSignMap,reffe,facet_own_dofs,cell)
   Dc = num_cell_dims(model)
   topo = get_grid_topology(model)
 
-  cell_facets = get_faces(topo, Dc, Dc-1)
+  # Type anotation needed when get_grid_topology isn't inferrable on model (eg DiscreteModelPortion).
+  cell_facets = get_faces(topo, Dc, Dc-1)::Table{Int32, Vector{Int32}, Vector{Int32}}
   cell_facets_cache = array_cache(cell_facets)
 
   return cell_facets, cell_facets_cache, CachedVector(Float64)
diff --git a/src/Fields/AutoDiff.jl b/src/Fields/AutoDiff.jl
index e467492a4..00f4c2204 100644
--- a/src/Fields/AutoDiff.jl
+++ b/src/Fields/AutoDiff.jl
@@ -1,69 +1,48 @@
 # Number differentiation
 
-function gradient(f::Number)
+function gradient(::V) where V<:Number
   function grad_f(x::Point)
-    zero(return_type(outer,x,f))
+    zero(gradient_type(V,x))
   end
 end
 
-function ∇∇(f::Number)
+function ∇∇(::V) where V<:Number
   function hess_f(x::Point)
-    g = gradient(f)(x)
-    gradient(g)(x)
+    zero(gradient_type(V,x,Val(2)))
   end
 end
 
 # Automatic differentiation of functions
 
-function gradient(f::Function)
-  function _gradient(x)
-    gradient(f,x)
-  end
-end
+# gradient(f) = x -> gradient(f,x)
+gradient(f::Function) = Base.Fix1(gradient,f)
+symmetric_gradient(f::Function) = Base.Fix1(symmetric_gradient,f)
+divergence(f::Function) = Base.Fix1(divergence,f)
+curl(f::Function) = Base.Fix1(curl,f)
+laplacian(f::Function) = Base.Fix1(laplacian,f)
 
-function symmetric_gradient(f::Function)
-  function _symmetric_gradient(x)
-    symmetric_gradient(f,x)
-  end
+function gradient(f::Function, x::Point)
+  fx = return_value(f, x)
+  gradient(f, x, fx)
 end
 
-function divergence(f::Function)
-  function _divergence(x)
-    divergence(f,x)
-  end
+function gradient(f::Function, x::Point{A}, fx::Number) where {A}
+  a = ForwardDiff.gradient(f∘Point, get_array(x))
+  VectorValue{A}(a)
 end
 
-function curl(f::Function)
-  function _curl(x)
-    curl(f,x)
-  end
-end
-
-function laplacian(f::Function)
-  function _laplacian(x)
-    laplacian(f,x)
-  end
-end
-
-# Specialization when x::Point
-
-function gradient(f::Function,x::Point)
-  gradient(f,x,return_value(f,x))
-end
-
-function gradient(f::Function,x::Point,fx::Number)
-  VectorValue(ForwardDiff.gradient(f∘Point,get_array(x)))
-end
-
-function gradient(f::Function,x::Point,fx::VectorValue)
-  TensorValue(transpose(ForwardDiff.jacobian(get_array∘f∘Point,get_array(x))))
+function gradient(f::Function, x::Point{A}, fx::VectorValue{B,T}) where {A,B,T}
+  a = transpose(ForwardDiff.jacobian(get_array∘f∘Point, get_array(x)))
+  # This type assert is necessary because of type inferrability issue in the forwardDiff call, for composed gradients
+  TensorValue{A,B}(a)::TensorValue{A,B,T,A*B}
 end
 
 # Implementation for all second order tensor values
 # Does not exploit possible symmetries
-function gradient(f::Function,x::Point{A},fx::S) where S<:MultiValue{Tuple{B,C}} where {A,B,C}
-  a = transpose(ForwardDiff.jacobian(get_array∘f∘Point,get_array(x)))
-  ThirdOrderTensorValue(reshape(a, (A,B,C)))
+function gradient(f::Function, x::Point{A}, fx::S) where S<:MultiValue{Tuple{B,C}} where {A,B,C}
+  a = transpose(ForwardDiff.jacobian(get_array∘f∘Point, get_array(x)))
+  a = reshape(a, (A,B,C))
+  ThirdOrderTensorValue{A,B,C}(a)
 end
 
 function gradient(f::Function,x::Point,fx::MultiValue)
@@ -108,7 +87,7 @@ function divergence(f::Function,x::Point,fx::TensorValue{3,3})
     a[1,1]+a[2,2]+a[3,3],
     a[4,1]+a[5,2]+a[6,3],
     a[7,1]+a[8,2]+a[9,3],
-   )
+  )
 end
 
 function divergence(f::Function,x::Point{2},fx::AbstractSymTensorValue{2})
@@ -127,7 +106,7 @@ function divergence(f::Function,x::Point{3},fx::AbstractSymTensorValue{3})
     a[1,1]+a[2,2]+a[3,3],
     a[2,1]+a[4,2]+a[5,3],
     a[3,1]+a[5,2]+a[6,3],
-   )
+  )
 end
 
 function divergence(f::Function,x::Point,fx::MultiValue)
@@ -149,7 +128,7 @@ end
 function laplacian(f::Function,x::Point{A},fx::VectorValue{B,T}) where {A,B,T}
   a = MMatrix{A*B,A,T}(undef)
   ForwardDiff.jacobian!(a, y->ForwardDiff.jacobian(get_array∘f∘Point,y), get_array(x))
-  VectorValue{B,T}( ntuple(k -> sum(i-> a[(i-1)*B+k,i], 1:A),B) )
+  VectorValue{B,T}( ntuple(k -> sum(i-> a[(i-1)*B+k,i], 1:A), B) )
 end
 
 function laplacian(f::Function,x::Point{A},fx::S) where S<:MultiValue{Tuple{B,C},T} where {A,B,C,T}
diff --git a/src/Fields/FieldArrays.jl b/src/Fields/FieldArrays.jl
index 30f1489d3..f279d1d58 100644
--- a/src/Fields/FieldArrays.jl
+++ b/src/Fields/FieldArrays.jl
@@ -204,7 +204,11 @@ end
 
 function testvalue(::Type{LinearCombinationField{V,F}}) where {V,F}
   fields = testvalue(F)
-  values = zeros(eltype(V), length(fields), 0)
+  # This method breaks the contract of testvalue for V not subtype of Array
+  #
+  # Ideally, we would extend the API of testvalue with a kwarg that enables
+  # selecting the size of the array when the argument is an array.
+  values = zeros(eltype(V), length(fields), tfill(0,Val(ndims(V)-1))...)
   LinearCombinationField(values,fields,1)
 end
 
@@ -258,12 +262,12 @@ end
 
 function Arrays.testitem(f::LinearCombinationFieldVector{V}) where V
   if !iszero(size(f.values,2))
-    values = f.values
+    values = f.values::V
   elseif V <: Diagonal
     values = Diagonal(zeros(eltype(V), 1))
   else
-    values = zeros(eltype(f.values),size(f.values,1),1)
-  end::V
+    values = zeros(eltype(f.values),size(f.values,1),1)::V
+  end
   LinearCombinationField(values,f.fields,1)
 end
 
@@ -856,3 +860,4 @@ function evaluate!(c,f::ConstantFieldArray{T},x::AbstractArray{<:Point}) where T
   end
   return r
 end
+
diff --git a/src/Fields/Fields.jl b/src/Fields/Fields.jl
index 5f1922440..58d66f5a6 100644
--- a/src/Fields/Fields.jl
+++ b/src/Fields/Fields.jl
@@ -14,7 +14,7 @@ import Gridap.Arrays: testitem
 
 using Gridap.Helpers
 using Gridap.Helpers: @abstractmethod, @notimplemented
-using Gridap.Helpers: @notimplementedif, @unreachable, @check
+using Gridap.Helpers: @notimplementedif, @unreachable, @check, @check_inferred
 using Gridap.Helpers: tfill, first_and_tail
 
 using Gridap.Algebra: mul!
@@ -31,6 +31,7 @@ using Test
 using StaticArrays
 using LinearAlgebra
 using AutoHashEquals: @auto_hash_equals as @ahe
+using Base: Fix1, Fix2
 
 import LinearAlgebra: det, inv, transpose, tr, cross
 import LinearAlgebra: ⋅, dot
diff --git a/src/Fields/FieldsInterfaces.jl b/src/Fields/FieldsInterfaces.jl
index 367bb3202..a991c2f3c 100644
--- a/src/Fields/FieldsInterfaces.jl
+++ b/src/Fields/FieldsInterfaces.jl
@@ -124,13 +124,25 @@ function push_∇∇(∇∇a::Field,ϕ::Field)
 end
 
 """
-    gradient_type(::Type{T}, x::Point) where T
+    gradient_type(::Type{T}, x::Point)
+    gradient_type(::Type{T}, x::Point, ::Val{N})
 
-Tensor type of the gradient of a `T` valued function.
+Tensor type of the gradient (repeated `N` times) of a `T` valued function.
 """
-function gradient_type(::Type{T},x::Point) where T
-  typeof(outer(zero(x),zero(T)))
+gradient_type(::Type{T}, x::Point)           where T = gradient_type(T,x,Val(1))
+gradient_type(::Type{T}, x::Point, ::Val{0}) where T  = typeof(zero(T)*zero(eltype(x)))
+gradient_type(::Type{T}, x::Point, ::Val{1}) where T  = typeof(outer(zero(x),zero(T)))
+function gradient_type(::Type{T}, x::Point, ::Val{2}) where T
+  G = typeof(outer(zero(x),zero(T)))
+  typeof(outer(zero(x),zero(G)))
 end
+function gradient_type(::Type{T}, x::Point, ::Val{3}) where T
+  G = typeof(outer(zero(x),zero(T)))
+  H = typeof(outer(zero(x),zero(G)))
+  typeof(outer(zero(x),zero(H)))
+end
+# I tried a recurive implementation, but that wasn't popely inferrable, so I removed it
+gradient_type(::Type{T}, x::Point, ::Val) where T = @notimplemented
 
 """
     struct FieldGradient{N,F} <: Field
@@ -360,8 +372,11 @@ end
 
 ## Make Function behave like Field
 
-return_cache(f::FieldGradient{N,<:Function},x::Point) where N = gradient(f.object,Val(N))
-evaluate!(c,f::FieldGradient{N,<:Function},x::Point) where N = c(x)
+function return_cache(f::FieldGradient{N,<:Function},x::Point) where N
+  cache = gradient(f.object,Val(N))
+  cache
+end
+evaluate!(cache, f::FieldGradient{N,<:Function}, x::Point) where N = cache(x)
 
 # Operations
 
@@ -393,13 +408,13 @@ function testvalue(::Type{OperationField{O,F}}) where {O,F<:Tuple}
 end
 
 function return_value(c::OperationField,x::Point)
-  fx = map(f -> return_value(f,x),c.fields)
+  fx = map(Fix2(return_value,x), c.fields)
   return return_value(c.op,fx...)
 end
 
 function return_cache(c::OperationField,x::Point)
-  cl = map(fi -> return_cache(fi,x),c.fields)
-  lx = map(fi -> return_value(fi,x),c.fields)
+  cl = map(Fix2(return_cache,x), c.fields)
+  lx = map(Fix2(return_value,x), c.fields)
   ck = return_cache(c.op,lx...)
   return ck, cl
 end
@@ -411,12 +426,12 @@ function evaluate!(cache,c::OperationField,x::Point)
 end
 
 function return_value(c::OperationField,x::AbstractArray{<:Point})
-  fx = map(f -> return_value(f,x),c.fields)
+  fx = map(Fix2(return_value,x), c.fields)
   return return_value(Broadcasting(c.op),fx...)
 end
 
 function return_cache(c::OperationField,x::AbstractArray{<:Point})
-  cf = map(fi -> return_cache(fi,x),c.fields)
+  cf = map(Fix2(return_cache,x), c.fields)
   lx = map((ci,fi) -> evaluate!(ci,fi,x),cf,c.fields)
   ck = return_cache(Broadcasting(c.op),lx...)
   return cf, ck
@@ -793,12 +808,12 @@ argument `grad` can be used. It should contain the result of evaluating `gradien
 Idem for `gradgrad`. The checks are performed with the `@test` macro.
 """
 function test_field(f::Field, x, v, cmp=(==); grad=nothing, gradgrad=nothing)
-  test_map(v,f,x;cmp=cmp)
+  test_map(v,f,x;cmp)
   if grad != nothing
-    test_map(grad,∇(f),x;cmp=cmp)
+    test_map(grad,∇(f),x;cmp)
   end
   if gradgrad != nothing
-    test_map(gradgrad,∇∇(f),x;cmp=cmp)
+    test_map(gradgrad,∇∇(f),x;cmp)
   end
   true
 end
diff --git a/src/Geometry/DiscreteModels.jl b/src/Geometry/DiscreteModels.jl
index 3ae038939..59cc4fbcc 100644
--- a/src/Geometry/DiscreteModels.jl
+++ b/src/Geometry/DiscreteModels.jl
@@ -1,5 +1,5 @@
 """
-    abstract type DiscreteModel{Dc,Dp} <: Grid
+    abstract type DiscreteModel{Dc,Dp} <: Grid{Dc,Dp}
 
 Abstract type holding information about a physical grid,
 the underlying grid topology, and a labeling of
@@ -73,21 +73,6 @@ get_reffes(g::DiscreteModel) = get_reffes(get_grid(g))
 
 # Default API
 
-"""
-    num_dims(model::DiscreteModel)
-"""
-num_dims(model::DiscreteModel) = num_dims(get_grid_topology(model))
-
-"""
-    num_cell_dims(model::DiscreteModel)
-"""
-num_cell_dims(model::DiscreteModel) = num_cell_dims(get_grid_topology(model))
-
-"""
-    num_point_dims(model::DiscreteModel)
-"""
-num_point_dims(model::DiscreteModel) = num_point_dims(get_grid_topology(model))
-
 """
     num_faces(g::DiscreteModel,d::Integer)
     num_faces(g::DiscreteModel)
diff --git a/src/Geometry/FaceLabelings.jl b/src/Geometry/FaceLabelings.jl
index f54ff6f0a..80ee155d5 100644
--- a/src/Geometry/FaceLabelings.jl
+++ b/src/Geometry/FaceLabelings.jl
@@ -545,7 +545,7 @@ function face_labeling_from_cell_tags(
   # Face entities:
   n_entities = n_tags
   to_key(e) = sort!(unique!(filter!(!isequal(UNSET), collect(Int32,e))))
-  entities = Dict{UInt64,Int32}([hash(to_key([i])) => i for i in 1:n_tags])
+  entities = Dict{UInt,Int32}([hash(to_key([i])) => i for i in 1:n_tags])
   for d in D-1:-1:0
     if split_dimensions
       empty!(entities) # Reset entities
diff --git a/src/Helpers/Helpers.jl b/src/Helpers/Helpers.jl
index bc5dacdec..661703295 100644
--- a/src/Helpers/Helpers.jl
+++ b/src/Helpers/Helpers.jl
@@ -5,6 +5,7 @@ $(public_names_in_md(@__MODULE__))
 """
 module Helpers
 using DocStringExtensions
+using Test
 
 @static if VERSION >= v"1.6"
   using Preferences
@@ -18,6 +19,7 @@ export @notimplemented
 export @notimplementedif
 export @unreachable
 export @check
+export @check_inferred
 export tfill
 export get_val_parameter
 export first_and_tail
diff --git a/src/Helpers/Macros.jl b/src/Helpers/Macros.jl
index 8ba21d6dc..78c166f7f 100644
--- a/src/Helpers/Macros.jl
+++ b/src/Helpers/Macros.jl
@@ -66,3 +66,67 @@ macro check(test,msg="A check failed")
     end
   end
 end
+
+"""
+    @check_inferred expression
+
+Macro used unsure that `expression` is inferrable (passes [`Test.@inferred`](@ref)).
+The check gets deactivated when running Gridap in performance mode.
+"""
+macro check_inferred(ex)
+  ex.head != :call && error("@check_inferred requires a call expression, got $(string(ex))")
+  ex_str = string(ex)
+  args_symbs = map(ex.args[2:end]) do arg_ex
+    string(arg_ex)
+  end
+  args_ex = map(ex.args[2:end]) do arg_ex
+    esc(arg_ex)
+  end
+
+  @static if execution_mode == "debug"
+    quote
+      try
+        $(esc( :(Helpers.Test.@inferred $(ex)) ))
+      catch e
+        if e isa InferrabilityException
+          # This ensures that we track the innermost inferrability failure of the tree
+          rethrow(e)
+
+        elseif e isa ErrorException && contains(e.msg , "does not match inferred return type")
+          args_symbs = $args_symbs
+          args_types = map(string∘typeof, tuple($(args_ex...)))
+          args_vals = tuple($(args_ex...))
+          ids = 1:length(args_symbs)
+          args_str = map(ids, args_symbs, args_types, args_vals) do i, symb, T, val
+            "    $symb::$T = $val\n"
+          end
+          ret_type = typeof($(esc(ex)))
+
+          msg = """
+          InferrabilityException: The following call is not type-inferrable:
+              $($ex_str)
+
+          Returned value:
+          $ret_type
+
+          Call arguments:
+          $(join(args_str, "\n"))
+          """
+          rethrow(InferrabilityException(msg))
+
+        else
+          rethrow()
+        end
+      end # catch
+    end # quote
+  else
+    esc(ex)
+  end
+end
+
+struct InferrabilityException <: Exception
+  msg
+end
+
+Base.showerror(io::IO, e::InferrabilityException) = print(io, e.msg)
+
diff --git a/src/ODEs/ODEs.jl b/src/ODEs/ODEs.jl
index 2be3e347d..760a3542d 100644
--- a/src/ODEs/ODEs.jl
+++ b/src/ODEs/ODEs.jl
@@ -13,6 +13,7 @@ using SparseArrays
 using BlockArrays
 using NLsolve
 using ForwardDiff
+using Base: Fix1
 
 using Gridap.Helpers
 using Gridap.Algebra
diff --git a/src/ODEs/TimeDerivatives.jl b/src/ODEs/TimeDerivatives.jl
index 78f7608d7..503958091 100644
--- a/src/ODEs/TimeDerivatives.jl
+++ b/src/ODEs/TimeDerivatives.jl
@@ -105,20 +105,24 @@ end
 # Specialisation for `Function` #
 #################################
 function time_derivative(f::Function)
+  # dfdt = t ->  ( x -> _time_derivative(f, t, x) )
+  #dfdt(t) = Fix1(Fix1(_time_derivative, f), t)
   function dfdt(t)
-    ft = f(t)
-    function dfdt_t(x)
-      T = return_type(ft, x)
-      _time_derivative(T, f, t, x)
-    end
+    # dfdt_t(x) = _time_derivative(f, t, x)
+    dfdt_t(x) = Fix1(Fix1(_time_derivative, f), t)(x)
     dfdt_t
   end
   dfdt
 end
 
+function _time_derivative(f, t, x)
+  T = return_type(f(t), x)
+  _time_derivative(T, f, t, x)
+end
+
 function _time_derivative(T::Type{<:Real}, f, t, x)
   partial(t) = f(t)(x)
-  ForwardDiff.derivative(partial, t)
+  T(ForwardDiff.derivative(partial, t))::T
 end
 
 function _time_derivative(T::Type{<:MultiValue}, f, t, x)
@@ -126,6 +130,9 @@ function _time_derivative(T::Type{<:MultiValue}, f, t, x)
   T(ForwardDiff.derivative(partial, t))
 end
 
+#Temporary for Julia 1.10 and 1.11 (Base.Fix1 did support 3 arg functions)
+(dfdt::Base.Fix1{typeof(_time_derivative), <:Function})(t, x) = _time_derivative(dfdt.x, t, x)
+
 ##########################################
 # Specialisation for `TimeSpaceFunction` #
 ##########################################
diff --git "a/src/Polynomials/BarycentricP\316\233Bases.jl" "b/src/Polynomials/BarycentricP\316\233Bases.jl"
index 9151cfe5e..c53cf17f1 100644
--- "a/src/Polynomials/BarycentricP\316\233Bases.jl"
+++ "b/src/Polynomials/BarycentricP\316\233Bases.jl"
@@ -437,11 +437,8 @@ function _return_cache(b::_BaryPΛBasis,x,::Type{G},::Val{N_deriv}) where {G,N_d
   cB = CachedVector(zeros(T,ndof_bernstein))
   # Cache for derivatives of Bα (∇Bα or HBα)
   if N_deriv > 0
-    DB = T
     xi = testitem(x)
-    for _ in 1:N_deriv
-      DB = gradient_type(DB,xi)
-    end
+    DB = gradient_type(T,xi,Val(N_deriv))
     # Cache for N_deriv's derivatives all of scalar nD-Bernstein polynomials
     t = (( nothing for _ in 2:N_deriv)..., CachedArray(zeros(DB,(1,ndof_bernstein))))
     s = MArray{Tuple{size(DB)...},T}(undef)
diff --git a/src/Polynomials/ModalC0Bases.jl b/src/Polynomials/ModalC0Bases.jl
index d8aa5a7f8..5cfbb3ce6 100644
--- a/src/Polynomials/ModalC0Bases.jl
+++ b/src/Polynomials/ModalC0Bases.jl
@@ -167,7 +167,7 @@ function _sort_by_nfaces!(terms::Vector{CartesianIndex{D}},orders) where D
   rang_nfaces = map(f,bin_ids_nfaces)
   rang_own_dofs = map(g,bin_ids_nfaces)
 
-  P = Int64[]
+  P = Int[]
   for i = 1:length(bin_ids_nfaces)
     cis_nfaces = CartesianIndices(rang_nfaces[i])
     cis_own_dofs = CartesianIndices(rang_own_dofs[i])
diff --git a/src/Polynomials/PolynomialInterfaces.jl b/src/Polynomials/PolynomialInterfaces.jl
index c58433a33..97ab66611 100644
--- a/src/Polynomials/PolynomialInterfaces.jl
+++ b/src/Polynomials/PolynomialInterfaces.jl
@@ -165,10 +165,8 @@ function return_cache(
   @assert D == length(eltype(x)) "Incorrect number of point components"
   f = fg.fa
   xi = testitem(x)
-  G = _return_val_eltype(f,x)
-  for _ in 1:N
-    G = gradient_type(G,xi)
-  end
+  VV = _return_val_eltype(f,x)
+  G = gradient_type(VV,xi,Val(N))
   _return_cache(f,x,G,Val(N))
 end
 
diff --git a/src/ReferenceFEs/Dofs.jl b/src/ReferenceFEs/Dofs.jl
index b5fdab866..03492c258 100644
--- a/src/ReferenceFEs/Dofs.jl
+++ b/src/ReferenceFEs/Dofs.jl
@@ -115,14 +115,10 @@ function Base.getindex(b::ConcatenatedDofVector, i::Integer)
 end
 
 function Arrays.return_cache(b::ConcatenatedDofVector, field)
-  args_caches = ()
-  Ts = ()
-  for dofs_i in b.args
-    ci = return_cache(dofs_i, field)
+  args_caches = map(dofs_i -> return_cache(dofs_i, field), b.args)
+  Ts = map(b.args, args_caches) do dofs_i, ci
     vals = evaluate!(ci, dofs_i, field)
-    Ti = eltype(vals)
-    Ts = (Ts..., Ti)
-    args_caches = (args_caches..., ci)
+    eltype(vals)
   end
 
   r_size = field isa AbstractVector ? (length(b),length(field)) : (length(b),)
diff --git a/src/ReferenceFEs/GeneralPolytopes.jl b/src/ReferenceFEs/GeneralPolytopes.jl
index 5c810d87e..140260382 100644
--- a/src/ReferenceFEs/GeneralPolytopes.jl
+++ b/src/ReferenceFEs/GeneralPolytopes.jl
@@ -503,7 +503,7 @@ function get_dimranges(p::GeneralPolytope)
   p.dimranges
 end
 
-function get_dimrange(p::GeneralPolytope,d::Integer)
+function get_dimrange(p::GeneralPolytope,d::Integer)::UnitRange{Int}
   setup_dimranges!(p)
   p.dimranges[d+1]
 end
diff --git a/src/ReferenceFEs/PolytopalQuadratures.jl b/src/ReferenceFEs/PolytopalQuadratures.jl
index 96421d59b..aff99b16f 100644
--- a/src/ReferenceFEs/PolytopalQuadratures.jl
+++ b/src/ReferenceFEs/PolytopalQuadratures.jl
@@ -49,10 +49,10 @@ end
 
 struct PQuadCoordsMap <: Map end
 
-function Arrays.return_cache(::PQuadCoordsMap, q::PolytopalQuadrature)
+function Arrays.return_cache(::PQuadCoordsMap, q::PolytopalQuadrature{D}) where D
   conn, coords = q.conn, q.coords
-  cmap = affine_map(Tuple(coords[first(conn)]))
-  
+  cmap = affine_map(NTuple{D+1}(coords[first(conn)]))
+
   x = get_coordinates(q.quad)
   cmap_cache = return_cache(cmap,x)
   y = evaluate!(cmap_cache,cmap,x)
@@ -61,7 +61,7 @@ function Arrays.return_cache(::PQuadCoordsMap, q::PolytopalQuadrature)
   return y_cache, cmap_cache
 end
 
-function Arrays.evaluate!(cache,::PQuadCoordsMap, q::PolytopalQuadrature)
+function Arrays.evaluate!(cache,::PQuadCoordsMap, q::PolytopalQuadrature{D}) where D
   y_cache, cmap_cache = cache
   conn, coords = q.conn, q.coords
 
@@ -71,7 +71,7 @@ function Arrays.evaluate!(cache,::PQuadCoordsMap, q::PolytopalQuadrature)
 
   nx = length(x)
   for (k,verts) in enumerate(conn)
-    cmap = affine_map(Tuple(coords[verts]))
+    cmap = affine_map(NTuple{D+1}(coords[verts]))
     y[(k-1)*nx+1:k*nx] .= evaluate!(cmap_cache, cmap, x)
   end
   return y
@@ -85,7 +85,7 @@ function Arrays.return_cache(::PQuadWeightsMap, q::PolytopalQuadrature)
   return CachedArray(similar(w,(length(w)*length(conn))))
 end
 
-function Arrays.evaluate!(cache,::PQuadWeightsMap, q::PolytopalQuadrature)
+function Arrays.evaluate!(cache,::PQuadWeightsMap, q::PolytopalQuadrature{D}) where D
   conn, coords = q.conn, q.coords
 
   w = get_weights(q.quad)
@@ -94,7 +94,7 @@ function Arrays.evaluate!(cache,::PQuadWeightsMap, q::PolytopalQuadrature)
 
   nw = length(w)
   for (k,verts) in enumerate(conn)
-    dV = meas(affine_map(Tuple(coords[verts])).gradient)
+    dV = meas(affine_map(NTuple{D+1}(coords[verts])).gradient)
     z[(k-1)*nw+1:k*nw] .= w .* dV
   end
   return z
diff --git a/src/ReferenceFEs/Polytopes.jl b/src/ReferenceFEs/Polytopes.jl
index 9c85e52ed..da9a3486d 100644
--- a/src/ReferenceFEs/Polytopes.jl
+++ b/src/ReferenceFEs/Polytopes.jl
@@ -118,12 +118,12 @@ function get_dimranges(p::Polytope)
 end
 
 """
-    get_dimrange(p::Polytope, d::Integer)
+    get_dimrange(p::Polytope, d::Integer)::UnitRange{Int64}
 
 Indices range of the `d`-dimensional faces of `p`. Equivalent to
 [`get_dimranges(p)[d+1]`](@ref get_dimrange).
 """
-function get_dimrange(p::Polytope, d::Integer)
+function get_dimrange(p::Polytope, d::Integer)::UnitRange{Int}
   get_dimranges(p)[d+1]
 end
 
diff --git a/test/AdaptivityTests/ComplexChangeDomainTests.jl b/test/AdaptivityTests/ComplexChangeDomainTests.jl
index 275932ec5..50cff4d1f 100644
--- a/test/AdaptivityTests/ComplexChangeDomainTests.jl
+++ b/test/AdaptivityTests/ComplexChangeDomainTests.jl
@@ -193,9 +193,9 @@ function get_model3(model2)
         Gridap.Geometry.Oriented(),
         has_affine_map=true))
 
-    n2o_faces_map = Vector{Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int64}}}(undef, 3)
-    n2o_faces_map[3] = Gridap.Arrays.Table(Vector{Int64}[[1], [2, 3, 4, 5]])
-    n2o_cell_to_child_id = Gridap.Arrays.Table(Vector{Int64}[[1], [1, 2, 3, 4]])
+    n2o_faces_map = Vector{Gridap.Arrays.Table{Int, Vector{Int}, Vector{Int}}}(undef, 3)
+    n2o_faces_map[3] = Gridap.Arrays.Table(Vector{Int}[[1], [2, 3, 4, 5]])
+    n2o_cell_to_child_id = Gridap.Arrays.Table(Vector{Int}[[1], [1, 2, 3, 4]])
 
 
     ref_rule_1, ref_rule_2 = get_ref_rules()
@@ -241,4 +241,4 @@ uh2 = interpolate(uh1, Uh2)
 dΩ = Measure(t3, 2)
 @test sum(∫((uh2-u)⋅(uh2-u))dΩ) < 1e-10
 
-end
\ No newline at end of file
+end
diff --git a/test/ArraysTests/KeyToValMapsTests.jl b/test/ArraysTests/KeyToValMapsTests.jl
index 2e7c6847d..bda052c85 100644
--- a/test/ArraysTests/KeyToValMapsTests.jl
+++ b/test/ArraysTests/KeyToValMapsTests.jl
@@ -20,7 +20,7 @@ _r = key_to_val.(keys)
 @test all(_r .== collect(r))
 
 
-key_to_val(i) = LagrangianRefFE(Float64,HEX,i)
+key_to_val(i) = LagrangianRefFE(Float64,HEX,i)::GenericLagrangianRefFE{GradConformity, 3}
 
 m = KeyToValMap(key_to_val)
 keys = [1,2,3,4,5,4,3,2,1]
diff --git a/test/CellDataTests/DiracDeltasTests.jl b/test/CellDataTests/DiracDeltasTests.jl
index c71d0b07d..816c8b375 100644
--- a/test/CellDataTests/DiracDeltasTests.jl
+++ b/test/CellDataTests/DiracDeltasTests.jl
@@ -48,14 +48,12 @@ p = Point(0.2,0.3)
 δ = DiracDelta(model,p)
 @test sum(δ(v)) ≈ v(p)
 
-@show δ = DiracDelta(Ω,p)
 @test sum(δ(v)) ≈ v(p)
 
 pvec = [p,π*p]
 δ = DiracDelta(model,pvec)
 @test sum(δ(v)) ≈ sum(v(pvec))
 
-@show δ = DiracDelta(Ω,pvec)
 @test sum(δ(v)) ≈ sum(v(pvec))
 
 p = Point(0.2,0.3)
diff --git a/test/FESpacesTests/CellFieldsTests.jl b/test/FESpacesTests/CellFieldsTests.jl
index 0c966ea7d..bbbb7df66 100644
--- a/test/FESpacesTests/CellFieldsTests.jl
+++ b/test/FESpacesTests/CellFieldsTests.jl
@@ -34,7 +34,9 @@ Random.seed!(0)
 
     coeff0 = rand(Float64)
     coeffs = rand(SVector{D,Float64})
-    f(x) = coeffs ⋅ SVector(Tuple(x)) + coeff0
+    f(x) = let coeffs=coeffs, coeff0=coeff0
+      coeffs ⋅ SVector(Tuple(x)) + coeff0
+    end
     # TODO: use this mechanism instead to project
     # Francesc Verdugo @fverdugo 13:11
     # a(u,v) = ∫( u*v )dΩ
@@ -124,7 +126,7 @@ source_model = CartesianDiscreteModel((0,1,0,1),(2,2))
   end
 
   # VectorValued Lagrangian
-  fᵥ(x) = VectorValue([x[1], x[1]+x[2]])
+  fᵥ(x) = VectorValue(x[1], x[1]+x[2])
   reffe = ReferenceFE(lagrangian, VectorValue{2,et}, 1)
   V₁ = FESpace(source_model, reffe, conformity=:H1)
   fh = interpolate_everywhere(fᵥ, V₁)
@@ -160,7 +162,7 @@ end
 
 @testset "Test interpolation RT" begin
   # RT Space -> RT Space
-  f(x) = VectorValue([x[1], x[2]])
+  f(x) = VectorValue(x[1], x[2])
   reffe = RaviartThomasRefFE(et, p, 0)
   V₁ = FESpace(source_model, reffe, conformity=:HDiv)
   fh = interpolate_everywhere(f, V₁);
@@ -186,7 +188,7 @@ source_model = CartesianDiscreteModel((0,1,0,1),(2,2)) |> simplexify
 @testset "Test interpolation BDM" begin
   # BDM Space -> BDM Space
 
-  f(x) = VectorValue([x[1], x[2]])
+  f(x) = VectorValue(x[1], x[2])
   reffe = BDMRefFE(et, p, 1)
   V₁ = FESpace(source_model, reffe, conformity=:HDiv)
   fh = interpolate_everywhere(f, V₁);
diff --git a/test/FESpacesTests/ConformingFESpacesTests.jl b/test/FESpacesTests/ConformingFESpacesTests.jl
index 7a03a2040..837a7537f 100644
--- a/test/FESpacesTests/ConformingFESpacesTests.jl
+++ b/test/FESpacesTests/ConformingFESpacesTests.jl
@@ -202,20 +202,33 @@ end
 # - (0/D)-form mapped: ~h⁰ (D-form same as 1-form because we don't use the broken Piola map in Gridap)
 # - 1-form mapped: ~h¹
 # - (D-1)-form mapped: ~hᴰ⁻¹
-vec_value(D) = (D > 1 ? x->VectorValue{D}(tfill(1.0,Val(D))) : x->1.0)
+function one_function(D)
+  if     D == 1
+    return x -> 1.
+  elseif D == 2
+    return x->VectorValue(1.,1.)
+  elseif D == 3
+    return x->VectorValue(1.,1.,1.)
+  elseif D == 4
+    return x->VectorValue(1.,1.,1.,1.)
+  else
+    @unreachable
+  end
+end
 
 reffe = ReferenceFE(nedelec, Float64, 3)
-_test_FESpace_dof_scaling(reffe, vec_value)
+_test_FESpace_dof_scaling(reffe, one_function)
+
 
 reffe = ReferenceFE(nedelec2, Float64, 3)
-_test_FESpace_dof_scaling(reffe, vec_value; n_cube=false)
+_test_FESpace_dof_scaling(reffe, one_function; n_cube=false)
 
 reffe = ReferenceFE(raviart_thomas, Float64, 3)
-_test_FESpace_dof_scaling(reffe, vec_value)
-_test_FESpace_dof_scaling(reffe, vec_value, dims=4:4, n_cube=false)
+_test_FESpace_dof_scaling(reffe, one_function)
+_test_FESpace_dof_scaling(reffe, one_function, dims=4:4, n_cube=false)
 
 reffe = ReferenceFE(bdm, Float64, 3)
-_test_FESpace_dof_scaling(reffe, vec_value; n_cube=false)
+_test_FESpace_dof_scaling(reffe, one_function; n_cube=false)
 
 # All trivial elements
 #
diff --git a/test/FESpacesTests/ConstantFESpaceTests.jl b/test/FESpacesTests/ConstantFESpaceTests.jl
index ab838f4b5..7f57326c6 100644
--- a/test/FESpacesTests/ConstantFESpaceTests.jl
+++ b/test/FESpacesTests/ConstantFESpaceTests.jl
@@ -11,7 +11,8 @@ model = CartesianDiscreteModel(domain,partition)
 Gridap.FESpaces.test_fe_space(Λ)
 M = TrialFESpace(Λ)
 
-order = 2
+let order = 2
+
 u((x,y)) = (x+y)^order
 f(x) = -Δ(u,x)
 reffe = ReferenceFE(lagrangian,Float64,order)
@@ -39,6 +40,8 @@ op2 = AffineFEOperator(a2,l2,M2,Λ2)
 μ2h = solve(op2)
 @assert sum(∫(μ2h⋅μ2h)dΩ) < 1.0e-12
 
+end # let order
+
 trian = Triangulation(model,[1,2,3,4])
 Λ3 = ConstantFESpace(trian,field_type=VectorValue{2,Float64})
 Gridap.FESpaces.test_fe_space(Λ3)
diff --git a/test/FESpacesTests/PolytopalFESpacesTests.jl b/test/FESpacesTests/PolytopalFESpacesTests.jl
index 512d18a15..e33a87a13 100644
--- a/test/FESpacesTests/PolytopalFESpacesTests.jl
+++ b/test/FESpacesTests/PolytopalFESpacesTests.jl
@@ -16,10 +16,10 @@ function test_l2_proj(model,V,order,u_exact)
 
   mass_lhs(u,v) = ∫(u⋅v)dΩ
   mass_rhs(v) = ∫(v⋅u_exact)dΩ
-  
+
   op = AffineFEOperator(mass_lhs,mass_rhs,U,V)
   uh = solve(op)
-  
+
   eh = uh - u_exact
   @test sum(∫(eh⋅eh)dΩ) < 1e-10
 
@@ -41,31 +41,31 @@ function test_dg_lap(model,V,order,u_exact)
   Ω = Triangulation(model)
   Γ = Boundary(model)
   Λ = Skeleton(model)
-  
+
   dΩ = Measure(Ω,2*order)
   dΓ = Measure(Γ,2*order)
   dΛ = Measure(Λ,2*order)
-  
+
   β = 100.0
   f(x) = -Δ(u_exact)(x)
   nΛ = get_normal_vector(Λ)
   βΛ = CellField(β ./ get_array(∫(1)dΛ),Λ)
   nΓ = get_normal_vector(Γ)
   βΓ = CellField(β ./ get_array(∫(1)dΓ),Γ)
-  
-  lap_lhs(u,v) = ∫(∇(u)⋅∇(v))dΩ + 
-                 ∫(βΛ*jump(u*nΛ)⋅jump(v*nΛ) - mean(∇(u))⋅jump(v*nΛ) - mean(∇(v))⋅jump(u*nΛ))dΛ + 
+
+  lap_lhs(u,v) = ∫(∇(u)⋅∇(v))dΩ +
+                 ∫(βΛ*jump(u*nΛ)⋅jump(v*nΛ) - mean(∇(u))⋅jump(v*nΛ) - mean(∇(v))⋅jump(u*nΛ))dΛ +
                  ∫(βΓ*u*v - (∇(u)⋅nΓ)*v - (∇(v)⋅nΓ)*u)dΓ
   lap_rhs(v) = ∫(v*f)dΩ + ∫(βΓ*u_exact*v - (∇(v)⋅nΓ)*u_exact)dΓ
-  
+
   op = AffineFEOperator(lap_lhs,lap_rhs,U,V)
   uh = solve(op)
-  
+
   eh = uh - u_exact
   @test sum(∫(eh⋅eh)dΩ) < 1e-10
 end
 
-order = 2
+let order = 2
 
 # 2D bulk
 u_exact_2d(x) = x[1]^order + x[2]^order
@@ -94,7 +94,7 @@ VΓ = FESpaces.PolytopalFESpace(Γ,VectorValue{2,Float64},order,space=:P,dirichl
 @test any(get_cell_dof_ids(VΓ).data .< 0)
 test_l2_proj(pmodel,VΓ,order,u_exact_2d_vec)
 
-# TODO: Does this make sense? 
+# TODO: Does this make sense?
 # VΓ = FESpaces.PolytopalFESpace(Γ,VectorValue{2,Float64},order,space=:P,dirichlet_tags=["boundary"],dirichlet_masks=[true,false])
 # @test any(get_cell_dof_ids(VΓ).data .< 0)
 # test_l2_proj(pmodel,VΓ,order,u_exact_2d_vec)
@@ -132,4 +132,6 @@ ptopo = Geometry.PatchTopology(get_grid_topology(model), Table([[1,2],[3,4]]))
 
 pspace = FESpaces.PatchFESpace(model,ptopo,Float64,order,space=:P)
 
-end
\ No newline at end of file
+end # let order
+
+end
diff --git a/test/FESpacesTests/ZeroMeanFESpacesTests.jl b/test/FESpacesTests/ZeroMeanFESpacesTests.jl
index 46b298c0e..fc9910ec0 100644
--- a/test/FESpacesTests/ZeroMeanFESpacesTests.jl
+++ b/test/FESpacesTests/ZeroMeanFESpacesTests.jl
@@ -38,13 +38,13 @@ uh = interpolate(f, U)
 @test abs(sum(∫(uh)*dΩ)) < 1.0e-10
 
 # Interpolate a function with zero mean
-ĝ(x) = x[1]^2 + x[2]
-g_mean = sum(∫(ĝ)*dΩ)/sum(∫(1)*dΩ)
-g(x) = ĝ(x) - g_mean
-vh = interpolate(g, U)
-eh = vh - g
-@test abs(sum(∫(vh)*dΩ)) < 1.0e-10
-@test sum(∫(eh*eh)*dΩ) < 1.0e-10
+let ĝ(x) = x[1]^2 + x[2], g_mean = sum(∫(ĝ)*dΩ)/sum(∫(1)*dΩ)
+  g(x) = ĝ(x) - g_mean
+  vh = interpolate(g, U)
+  eh = vh - g
+  @test abs(sum(∫(vh)*dΩ)) < 1.0e-10
+  @test sum(∫(eh*eh)*dΩ) < 1.0e-10
+end
 
 # Check Stokes/Navier-Stokes problem properties
 
diff --git a/test/FieldsTests/DiffOperatorsTests.jl b/test/FieldsTests/DiffOperatorsTests.jl
index 29c1ec011..066a62939 100644
--- a/test/FieldsTests/DiffOperatorsTests.jl
+++ b/test/FieldsTests/DiffOperatorsTests.jl
@@ -9,6 +9,7 @@ using Gridap.Fields: skew_symmetric_gradient, ShiftedNabla
 
 using LinearAlgebra
 using FillArrays
+using StaticArrays: @MMatrix
 
 np = 4
 p = Point(1,2)
@@ -66,7 +67,10 @@ f = Fill(f,l)
 @test Broadcasting(curl)(f) == Broadcasting(Operation(grad2curl))(Broadcasting(∇)(f))
 @test Broadcasting(ε)(f) == Broadcasting(Operation(symmetric_part))(Broadcasting(∇)(f))
 
-@test evaluate(Broadcasting(∇+p)(f),x) == evaluate( Broadcasting(Operation((g,f)->g+p⊗f))(Broadcasting(∇)(f),f)  ,x)
+let p=p
+  g_plus_pf = (g,f)->g+p⊗f
+  @test evaluate(Broadcasting(∇+p)(f),x) == evaluate( Broadcasting(Operation(g_plus_pf))(Broadcasting(∇)(f),f)  ,x)
+end
 
 # Test automatic differentiation
 
@@ -106,7 +110,7 @@ u_qten(x)  = SymTracelessTensorValue( x[1]^3 + 2x[2]^3, 5*x[1]^3 - 7x[2]^3)
 u_ten3(x)  = ThirdOrderTensorValue{2,1,2}( x[1]^2 + x[2], 4*x[1] - x[2]^2, -x[1]^2 - x[2], -4*x[1] + x[2]^2 )
 Δu_ten3(x) = ThirdOrderTensorValue{2,1,2}( 2, -2, -2, 2 )
 
-xs = [ Point(1.,1.), Point(2.,0.), Point(0.,3.), Point(-1.,3.)]
+xs = [ Point(1.,2.), Point(1.,1.), Point(2.,0.), Point(0.,3.), Point(-1.,3.)]
 for x in xs
   @test  ∇(u_scal)(x) == ∇u_scal(x)
   @test  Δ(u_scal)(x) == Δu_scal(x)
@@ -169,27 +173,28 @@ for x in xs
   @test (∇⋅εu)(x) == VectorValue(4.,0.)
 end
 
-x0 = VectorValue(-0.3, 0.5)
-A = TensorValue(0.0, -1.0, 1.0, 0.0)
-u(x) = (x - x0) ⋅ A
-∇u(x) = A
-Δu(x) = VectorValue(0.0, 0.0)
+let x0=VectorValue(-0.3, 0.5), A=TensorValue(0.0, -1.0, 1.0, 0.0)
+  u(x) = (x - x0) ⋅ A
+  ∇u(x) = A
+  Δu(x) = VectorValue(0.0, 0.0)
 
-for x in xs
-  @test ∇(u)(x) == ∇u(x)
-  @test Δ(u)(x) == Δu(x)
-  @test Δ(u)(x) == (∇⋅∇u)(x)
+  for x in xs
+    @test ∇(u)(x) == ∇u(x)
+    @test Δ(u)(x) == Δu(x)
+    @test Δ(u)(x) == (∇⋅∇u)(x)
+  end
 end
 
-x0 = VectorValue(1.0, 0.5)
-u(x) = exp(-(x - x0) ⋅ (x - x0))
-∇u(x) = -2(x - x0) * u(x)
-Δu(x) = 4u(x) * ((x - x0) ⋅ (x - x0) - 1)
+let x0=VectorValue(1.0, 0.5)
+  u(x) = exp(-(x - x0) ⋅ (x - x0))
+  ∇u(x) = -2(x - x0) * u(x)
+  Δu(x) = 4u(x) * ((x - x0) ⋅ (x - x0) - 1)
 
-for x in xs
-  @test ∇(u)(x) == ∇u(x)
-  @test Δ(u)(x) == Δu(x)
-  @test Δ(u)(x) == (∇⋅∇u)(x)
+  for x in xs
+    @test ∇(u)(x) == ∇u(x)
+    @test Δ(u)(x) == Δu(x)
+    @test Δ(u)(x) == (∇⋅∇u)(x)
+  end
 end
 
 u(x) = VectorValue( x[1]^2 + 2*x[3]^2, -x[1]^2, -x[2]^2 + x[3]^2 )
@@ -210,26 +215,28 @@ for x in xs
   @test (∇⋅εu)(x) == VectorValue(4.,-1.,1.)
 end
 
-v = VectorValue(0.5, -1.0, 2.0)
-x0 = VectorValue(0.3, 1.0, -0.2)
-u(x) = (x - x0) × v + v ⋅ (x ⊗ x)
-∇u(x) = (
-  TensorValue([0.0 -v[3] v[2]; v[3] 0.0 -v[1]; -v[2] v[1] 0.0])
-  + v ⊗ x + v ⋅ x * TensorValue(Matrix(I,3,3))
-)
-Δu(x) = 2v
-εu(x) = symmetric_part(∇u(x))
-νu(x) = TensorValues.skew_symmetric_part(∇u(x))
-
-xs = [ Point(1.,1.,2.0), Point(2.,0.,1.), Point(0.,3.,0.), Point(-1.,3.,2.)]
-for x in xs
-  @test ∇(u)(x) == ∇u(x)
-  @test (∇⋅u)(x) == tr(∇u(x))
-  @test (∇×u)(x) == grad2curl(∇u(x))
-  @test Δ(u)(x) == Δu(x)
-  @test ε(u)(x) == εu(x)
-  @test Fields.skew_symmetric_gradient(u)(x) == νu(x)
-  @test Δ(u)(x) == (∇⋅∇u)(x)
+let v=VectorValue(0.5, -1.0, 2.0), x0=VectorValue(0.3, 1.0, -0.2)
+  u(x) = (x - x0) × v + v ⋅ (x ⊗ x)
+  function ∇u(x)
+      a = TensorValue(@MMatrix [0.0 -v[3] v[2]; v[3] 0.0 -v[1]; -v[2] v[1] 0.0] )
+      b =  v ⊗ x + v ⋅ x * one(TensorValue{3,3,Float64})
+      a+b
+  end
+  Δu(x) = 2v
+  εu(x) = symmetric_part(∇u(x))
+  νu(x) = TensorValues.skew_symmetric_part(∇u(x))
+
+  xs = [ Point(1.,1.,2.0), Point(2.,0.,1.), Point(0.,3.,0.), Point(-1.,3.,2.)]
+  x = first(xs)
+  for x in xs
+    @test ∇(u)(x) == ∇u(x)
+    @test (∇⋅u)(x) == tr(∇u(x))
+    @test (∇×u)(x) == grad2curl(∇u(x))
+    @test Δ(u)(x) == Δu(x)
+    @test ε(u)(x) == εu(x)
+    @test Fields.skew_symmetric_gradient(u)(x) == νu(x)
+    @test Δ(u)(x) == (∇⋅∇u)(x)
+  end
 end
 
 # FieldArray diff ops
diff --git a/test/FieldsTests/FieldInterfacesTests.jl b/test/FieldsTests/FieldInterfacesTests.jl
index 214b0fbf9..d8da7bde8 100644
--- a/test/FieldsTests/FieldInterfacesTests.jl
+++ b/test/FieldsTests/FieldInterfacesTests.jl
@@ -404,4 +404,19 @@ test_field(∇f, p, ∇cp, grad=∇∇cp)
 test_field(∇f, x, ∇(f).(x), grad=∇∇(f).(x))
 test_field(∇f, z, ∇(f).(z), grad=∇∇(f).(z))
 
+p = Point(1., 1.)
+V = VectorValue{3,Float32}
+@test gradient_type(V,p) == TensorValue{2,3,Float64,6}
+@test gradient_type(V,p,Val(0)) == VectorValue{3,Float64}
+@test gradient_type(V,p,Val(1)) == TensorValue{2,3,Float64,6}
+@test gradient_type(V,p,Val(2)) == ThirdOrderTensorValue{2,2,3,Float64,12}
+@test gradient_type(V,p,Val(3)) == HighOrderTensorValue{Tuple{2,2,2,3},Float64,4,24}
+@test_throws "not yet implemented" gradient_type(V,p,Val(4))
+
+G = gradient_type(V,p)
+H = gradient_type(G,p)
+@test H == ThirdOrderTensorValue{2,2,3,Float64,12}
+J = gradient_type(H,p)
+@test J  == HighOrderTensorValue{Tuple{2,2,2,3},Float64,4,24}
+
 end # module
diff --git a/test/HelpersTests/MacrosTests.jl b/test/HelpersTests/MacrosTests.jl
index 99b136460..dc4012922 100644
--- a/test/HelpersTests/MacrosTests.jl
+++ b/test/HelpersTests/MacrosTests.jl
@@ -1,3 +1,19 @@
 module MacrosTests
 
+using Test
+using Gridap.Helpers
+
+y = 0
+f(x) = y+x
+@test_throws Helpers.InferrabilityException @check_inferred f(0)
+
+let y=y
+g(x) = y+x
+@test_nowarn @check_inferred g(0)
+end
+
+const z = 0
+h(x) = z+x
+@test_nowarn @check_inferred h(0)
+
 end # module
diff --git a/test/ODEsTests/TransientFEOperatorsSolutionsTests.jl b/test/ODEsTests/TransientFEOperatorsSolutionsTests.jl
index e176928fb..7c7f81cf7 100644
--- a/test/ODEsTests/TransientFEOperatorsSolutionsTests.jl
+++ b/test/ODEsTests/TransientFEOperatorsSolutionsTests.jl
@@ -104,7 +104,7 @@ end
 # First-order #
 ###############
 order = 1
-ft(t) = x -> ∂t(u)(t, x) - Δ(u)(t, x)
+ft(t) = let u=u; x -> ∂t(u)(t, x) - Δ(u)(t, x) end
 f = TimeSpaceFunction(ft)
 
 mass(t, ∂ₜu, v) = ∫(∂ₜu ⋅ v) * dΩ
@@ -200,7 +200,7 @@ end
 # Second-order #
 ################
 order = 2
-ft(t) = x -> ∂tt(u)(t, x) + ∂t(u)(t, x) - Δ(u)(t, x)
+ft(t) = let u=u; x -> ∂tt(u)(t, x) + ∂t(u)(t, x) - Δ(u)(t, x) end
 f = TimeSpaceFunction(ft)
 
 mass(t, ∂ₜₜu, v) = ∫(∂ₜₜu ⋅ v) * dΩ
diff --git a/test/ODEsTests/TransientFEProblemsTests/FreeSurfacePotentialFlowTests.jl b/test/ODEsTests/TransientFEProblemsTests/FreeSurfacePotentialFlowTests.jl
index 9da67f79d..3224ef617 100644
--- a/test/ODEsTests/TransientFEProblemsTests/FreeSurfacePotentialFlowTests.jl
+++ b/test/ODEsTests/TransientFEProblemsTests/FreeSurfacePotentialFlowTests.jl
@@ -6,21 +6,21 @@ using Gridap
 using Gridap.Geometry
 
 # Parameters
-L = 2 * π
-H = 1.0
-n = 8
-order = 2
-g = 9.81
-ξ = 0.1
-λ = L / 2
-k = 2 * π / L
-h = L / n
-ω = √(g * k * tanh(k * H))
-t0 = 0.0
-dt = h / (2 * λ * ω)
-tF = 10 * dt # 2 * π
-α = 2 / dt
-tol = 1.0e-2
+const L = 2 * π
+const H = 1.0
+const n = 8
+const order = 2
+const g = 9.81
+const ξ = 0.1
+const λ = L / 2
+const k = 2 * π / L
+const h = L / n
+const ω = √(g * k * tanh(k * H))
+const t0 = 0.0
+const dt = h / (2 * λ * ω)
+const tF = 10 * dt # 2 * π
+const α = 2 / dt
+const tol = 1.0e-2
 
 # Exact solution
 ϕₑ(t) = x -> ω / k * ξ * (cosh(k * (x[2]))) / sinh(k * H) * sin(k * x[1] - ω * t)
diff --git a/test/ODEsTests/TransientFEProblemsTests/HeatEquationDGTests.jl b/test/ODEsTests/TransientFEProblemsTests/HeatEquationDGTests.jl
index 1c0aa3167..5d977d34d 100644
--- a/test/ODEsTests/TransientFEProblemsTests/HeatEquationDGTests.jl
+++ b/test/ODEsTests/TransientFEProblemsTests/HeatEquationDGTests.jl
@@ -38,7 +38,7 @@ dΛ = Measure(Λ, degree)
 nΛ = get_normal_vector(Λ)
 
 # FE operator
-ft(t) = x -> ∂t(u)(t, x) - Δ(u)(t, x)
+ft(t) = let u=u; x -> ∂t(u)(t, x) - Δ(u)(t, x) end
 f = TimeSpaceFunction(ft)
 h = 1 / 5
 γ = order * (order + 1)
diff --git a/test/ODEsTests/TransientFEProblemsTests/HeatEquationMultiFieldTests.jl b/test/ODEsTests/TransientFEProblemsTests/HeatEquationMultiFieldTests.jl
index ed7c0616b..90322a63d 100644
--- a/test/ODEsTests/TransientFEProblemsTests/HeatEquationMultiFieldTests.jl
+++ b/test/ODEsTests/TransientFEProblemsTests/HeatEquationMultiFieldTests.jl
@@ -33,7 +33,7 @@ degree = 2 * order
 dΩ = Measure(Ω, degree)
 
 # FE operator
-ft(t) = x -> ∂t(u)(t, x) - Δ(u)(t, x)
+ft(t) = let u=u; x -> ∂t(u)(t, x) - Δ(u)(t, x) end
 f = TimeSpaceFunction(ft)
 _mass(t, ∂ₜu, v) = ∫(∂ₜu ⋅ v) * dΩ
 _mass(t, u, ∂ₜu, v) = _mass(t, ∂ₜu, v)
diff --git a/test/ODEsTests/TransientFEProblemsTests/HeatEquationNeumannTests.jl b/test/ODEsTests/TransientFEProblemsTests/HeatEquationNeumannTests.jl
index 167b77d03..3159c1785 100644
--- a/test/ODEsTests/TransientFEProblemsTests/HeatEquationNeumannTests.jl
+++ b/test/ODEsTests/TransientFEProblemsTests/HeatEquationNeumannTests.jl
@@ -36,7 +36,7 @@ dΓ = Measure(Γ, degree)
 nΓ = get_normal_vector(Γ)
 
 # FE operator
-ft(t) = x -> ∂t(u)(t, x) - Δ(u)(t, x)
+ft(t) = let u=u; x -> ∂t(u)(t, x) - Δ(u)(t, x) end
 f = TimeSpaceFunction(ft)
 mass(t, ∂ₜu, v) = ∫(∂ₜu ⋅ v) * dΩ
 mass(t, u, ∂ₜu, v) = mass(t, ∂ₜu, v)
diff --git a/test/ODEsTests/TransientFEProblemsTests/HeatEquationScalarTests.jl b/test/ODEsTests/TransientFEProblemsTests/HeatEquationScalarTests.jl
index d7d0ddd57..88465923e 100644
--- a/test/ODEsTests/TransientFEProblemsTests/HeatEquationScalarTests.jl
+++ b/test/ODEsTests/TransientFEProblemsTests/HeatEquationScalarTests.jl
@@ -30,7 +30,7 @@ degree = 2 * order
 dΩ = Measure(Ω, degree)
 
 # FE operator
-ft(t) = x -> ∂t(u)(t, x) - Δ(u)(t, x)
+ft(t) = let u=u; x -> ∂t(u)(t, x) - Δ(u)(t, x) end
 f = TimeSpaceFunction(ft)
 mass(t, ∂ₜu, v) = ∫(∂ₜu ⋅ v) * dΩ
 mass(t, u, ∂ₜu, v) = mass(t, ∂ₜu, v)
diff --git a/test/ODEsTests/TransientFEProblemsTests/HeatEquationVectorTests.jl b/test/ODEsTests/TransientFEProblemsTests/HeatEquationVectorTests.jl
index 86c75ab06..46c6f3426 100644
--- a/test/ODEsTests/TransientFEProblemsTests/HeatEquationVectorTests.jl
+++ b/test/ODEsTests/TransientFEProblemsTests/HeatEquationVectorTests.jl
@@ -30,7 +30,7 @@ degree = 2 * order
 dΩ = Measure(Ω, degree)
 
 # FE operator
-ft(t) = x -> ∂t(u)(t, x) - Δ(u)(t, x)
+ft(t) = let u=u; x -> ∂t(u)(t, x) - Δ(u)(t, x) end
 f = TimeSpaceFunction(ft)
 mass(t, ∂ₜu, v) = ∫(∂ₜu ⋅ v) * dΩ
 mass(t, u, ∂ₜu, v) = mass(t, ∂ₜu, v)
diff --git a/test/ODEsTests/TransientFEProblemsTests/SecondOrderEquationTests.jl b/test/ODEsTests/TransientFEProblemsTests/SecondOrderEquationTests.jl
index 6cc04f36f..80a26382d 100644
--- a/test/ODEsTests/TransientFEProblemsTests/SecondOrderEquationTests.jl
+++ b/test/ODEsTests/TransientFEProblemsTests/SecondOrderEquationTests.jl
@@ -30,7 +30,7 @@ dΩ = Measure(Ω, degree)
 
 # FE operator
 order = 2
-ft(t) = x -> ∂tt(u)(t, x) + ∂t(u)(t, x) - Δ(u)(t, x)
+ft(t) = let u=u; x -> ∂tt(u)(t, x) + ∂t(u)(t, x) - Δ(u)(t, x) end
 f = TimeSpaceFunction(ft)
 mass(t, ∂ₜₜu, v) = ∫(∂ₜₜu ⋅ v) * dΩ
 mass(t, u, ∂ₜₜu, v) = mass(t, ∂ₜₜu, v)
diff --git a/test/ODEsTests/TransientFEProblemsTests/StokesEquationTests.jl b/test/ODEsTests/TransientFEProblemsTests/StokesEquationTests.jl
index c78227991..be86b256b 100644
--- a/test/ODEsTests/TransientFEProblemsTests/StokesEquationTests.jl
+++ b/test/ODEsTests/TransientFEProblemsTests/StokesEquationTests.jl
@@ -40,7 +40,7 @@ degree = 2 * order
 dΩ = Measure(Ω, degree)
 
 # FE operator
-ft(t) = x -> ∂t(u)(t, x) - Δ(u)(t, x) + ∇(p)(t, x)
+ft(t) = let u=u, p=p; x -> ∂t(u)(t, x) - Δ(u)(t, x) + ∇(p)(t, x) end
 f = TimeSpaceFunction(ft)
 g = ∇ ⋅ u
 mass(t, ∂ₜu, v) = ∫(∂ₜu ⋅ v) * dΩ
diff --git a/test/ReferenceFEsTests/BDMRefFEsTests.jl b/test/ReferenceFEsTests/BDMRefFEsTests.jl
index f799800d1..e4d1a6d42 100644
--- a/test/ReferenceFEsTests/BDMRefFEsTests.jl
+++ b/test/ReferenceFEsTests/BDMRefFEsTests.jl
@@ -34,8 +34,7 @@ face_dofs = Vector{Int}[[],[],[],[1,2],[3,4],[5,6],[1,2,3,4,5,6]]
 prebasis = get_prebasis(reffe)
 dof_basis = get_dof_basis(reffe)
 
-v = VectorValue(3.0,0.0)
-field = GenericField(x->v*x[1])
+field = GenericField(x->VectorValue(3.0,0.0)*x[1])
 
 cache = return_cache(dof_basis,field)
 r = evaluate!(cache, dof_basis, field)
@@ -58,8 +57,7 @@ reffe = BDMRefFE(et,p,order)
 prebasis = get_prebasis(reffe)
 dof_basis = get_dof_basis(reffe)
 
-v = VectorValue(3.0,0.0)
-field = GenericField(x->v*x[1])
+field = GenericField(x->VectorValue(3.0,0.0)*x[1])
 
 cache = return_cache(dof_basis,field)
 r = evaluate!(cache, dof_basis, field)
@@ -97,8 +95,7 @@ test_reference_fe(reffe)
 prebasis = get_prebasis(reffe)
 dof_basis = get_dof_basis(reffe)
 
-v = VectorValue(0.0,3.0,0.0)
-field = GenericField(x->v)
+field = GenericField(x->VectorValue(0.0,3.0,0.0))
 
 cache = return_cache(dof_basis,field)
 r = evaluate!(cache, dof_basis, field)
@@ -120,8 +117,7 @@ test_reference_fe(reffe)
 prebasis = get_prebasis(reffe)
 dof_basis = get_dof_basis(reffe)
 
-v = VectorValue(0.0,3.0,0.0)
-field = GenericField(x->v)
+field = GenericField(x->VectorValue(0.0,3.0,0.0))
 
 cache = return_cache(dof_basis,field)
 r = evaluate!(cache, dof_basis, field)
diff --git a/test/ReferenceFEsTests/CrouzeixRaviartFEsTests.jl b/test/ReferenceFEsTests/CrouzeixRaviartFEsTests.jl
index a8f3e51ca..ce038770d 100644
--- a/test/ReferenceFEsTests/CrouzeixRaviartFEsTests.jl
+++ b/test/ReferenceFEsTests/CrouzeixRaviartFEsTests.jl
@@ -144,9 +144,7 @@ end
 
 
 partition = (0,1,0,1)
-
-# Stokes
-u_exact(x) = VectorValue( [sin(pi*x[1])^2*sin(pi*x[2])*cos(pi*x[2]), -sin(pi*x[2])^2*sin(pi*x[1])*cos(pi*x[1])] )
+u_exact(x) = VectorValue( sin(pi*x[1])^2*sin(pi*x[2])*cos(pi*x[2]), -sin(pi*x[2])^2*sin(pi*x[1])*cos(pi*x[1]) )
 p_exact(x) = sin(2*pi*x[1])*sin(2*pi*x[2])
 
 # Scalar Poisson
diff --git a/test/ReferenceFEsTests/GeometricDecompositionTests.jl b/test/ReferenceFEsTests/GeometricDecompositionTests.jl
index 99c4a1369..3a8cab23e 100644
--- a/test/ReferenceFEsTests/GeometricDecompositionTests.jl
+++ b/test/ReferenceFEsTests/GeometricDecompositionTests.jl
@@ -46,6 +46,8 @@ If put in src, this should be a supertype of `Dof`, the latter being a *linear*
 """
 abstract type Form <: Map end
 
+struct FaceIntegralForm <: Form end
+
 """
     FaceIntegralFormVector(
       p::Polytope{D},
@@ -66,22 +68,23 @@ In the final basis, the forms are ordered by moment descriptor, then by face.
 
 We are assuming that all the faces in a moment are of the same type.
 """
-struct FaceIntegralFormVector{T,FI,DS} <: AbstractVector{Form}
-  form_faces::Vector{Vector{NTuple{2,Int}}} # f -> faces over which f is integrated computed
+struct FaceIntegralFormVector{T,N,FI,DS} <: AbstractVector{FaceIntegralForm}
+  form_faces::NTuple{N,Vector{NTuple{2,Int}}} # f -> faces over which f is integrated computed
   form_integrands::FI # i -> σᵢ
   form_measures::DS # tuple of `FaceMeasure`s
   face_own_forms::Vector{Vector{Int}} # "inverse" of form_faces
 
   function FaceIntegralFormVector(faces,integrands,measures,f_own_forms)
+    N = length(faces)
+    @check N == length(integrands) == length(measures)
     T = Float64 # TODO
     FI = typeof(integrands)
     DS = typeof(measures)
-    new{T,FI,DS}(faces,integrands,measures,f_own_forms)
+    new{T,N,FI,DS}(faces,integrands,measures,f_own_forms)
   end
 
   function FaceIntegralFormVector(p::Polytope,order::Int,forms)
     n_faces = num_faces(p)
-    n_form = length(forms)
     face_dims = get_facedims(p)
     face_offsets = get_offsets(p)
     reffaces, face_types = ReferenceFEs._compute_reffaces_and_face_types(p)
@@ -89,28 +92,32 @@ struct FaceIntegralFormVector{T,FI,DS} <: AbstractVector{Form}
     # Create facemeasures and integrand for each integral form
     # Count number of forms per face
     # compute local to global form index
+    face_form_ids = ()
     integrands = ()
     measures = ()
     face_n_forms = zeros(Int,n_faces)
-    face_form_ids = Vector{Vector{NTuple{2,Int}}}(undef,n_form)
-    form_I = 0
-    for (k,(faces,σ)) in enumerate(forms)
+    form_I = [0]
+
+    face_form_ids = map(forms) do (faces, σ)
       face_ids = Vector{NTuple{2,Int}}(undef, length(faces))
       for (i,face) in enumerate(faces)
         d = face_dims[face]
         lface = face - face_offsets[d+1]
-        face_ids[i] = (form_I + i, lface)
+        face_ids[i] = (form_I[1] + i, lface)
       end
-      face_form_ids[k] = face_ids
-      form_I += length(faces)
-
+      form_I[1] += length(faces)
       face_n_forms[faces] .+= 1
-      integrands = (integrands..., σ)
+      face_ids
+    end
+    integrands = map(forms) do (faces, σ)
+      σ
+    end
 
+    measures = map(forms) do (faces, σ)
       ftype = face_types[first(faces)]
       @check all(isequal(ftype), face_types[faces])
       fp = reffaces[ftype]
-      measures = (measures..., FaceMeasure(p,fp,order) )
+      FaceMeasure(p,fp,order)
     end
 
     # Compute face forms and nodes indices
@@ -122,20 +129,19 @@ struct FaceIntegralFormVector{T,FI,DS} <: AbstractVector{Form}
       n_forms += n_forms_i
     end
 
-    FaceIntegralFormVector(face_form_ids,integrands,measures,face_own_forms)
+    FaceIntegralFormVector(Tuple(face_form_ids),integrands,measures,face_own_forms)
   end
 end
 
 Base.size(a::FaceIntegralFormVector) = (num_forms(a),)
 Base.axes(a::FaceIntegralFormVector) = (Base.OneTo(num_forms(a)),)
-Base.getindex(::FaceIntegralFormVector,::Integer) = Form()
+Base.getindex(::FaceIntegralFormVector,::Integer) = FaceIntegralForm()
 Base.IndexStyle(::FaceIntegralFormVector) = IndexLinear()
 num_forms(b::FaceIntegralFormVector) = mapreduce(length, +, b.face_own_forms)
 
-function return_cache(b::FaceIntegralFormVector{T}, f) where {T}
-  cms = ()
-  for (_,σ,ds) in zip(b.form_faces, b.form_integrands, b.form_measures)
-    cms = ( cms..., return_cache(σ,f,_μ,ds))
+function return_cache(b::FaceIntegralFormVector{T,N,FI,DS}, f) where {T,N,FI,DS}
+  cms = map(b.form_faces, b.form_integrands, b.form_measures) do faces, σ, ds
+    return_cache(σ,f,_μ,ds)
   end
 
   r = Array{T}(undef, (num_forms(b), size(f)...))
@@ -160,7 +166,6 @@ function evaluate!(cache, b::FaceIntegralFormVector, f::AbstractVector{<:Field})
   forms_vals
 end
 
-
 ###########################################################################################
 # Traces L2-norm descriptors and geometric decompositions trace-forms for each conformity #
 ###########################################################################################
@@ -210,8 +215,10 @@ end
 # norm form on each boundary face, zero inside
 function get_trace_forms(p::Polytope{D}, field, ::GradConformity) where D
   face_ranges = get_dimranges(p)
-  trace_forms = Tuple[ (face_ranges[d], HGrad_face_tr_norm) for d in 1:D ]
-  push!(trace_forms, (last(face_ranges), zero_moment))
+  trace_forms = ntuple(Val(D)) do d
+    (face_ranges[d], HGrad_face_tr_norm)
+  end
+  trace_forms = (trace_forms..., (last(face_ranges), zero_moment))
   qorder = 2*get_order(field)+1
   FaceIntegralFormVector(p,qorder,trace_forms)
 end
@@ -223,7 +230,7 @@ function get_trace_forms(p::Polytope{D}, field, ::CurlConformity) where D
   trace_forms = Tuple[ (first(face_ranges), zero_moment) ]
   D≥2 && push!(trace_forms, (face_ranges[2], HCurl_edge_tr_norm))
   D≥3 && push!(trace_forms, (face_ranges[D], HCurl_facet_tr_norm))
-  push!(  trace_forms, (last( face_ranges), zero_moment))
+  trace_forms = (trace_forms..., (last(face_ranges), zero_moment))
   qorder = 2*get_order(field)+1
   FaceIntegralFormVector(p,qorder,trace_forms)
 end
@@ -232,15 +239,15 @@ end
 function get_trace_forms(p::Polytope{D}, field, ::DivConformity) where D
   face_ranges = get_dimranges(p)
   trace_forms = Tuple[ (face_ranges[d], zero_moment) for d in 1:D-1 ]
-  push!(trace_forms, (face_ranges[D], HDiv_facet_tr_norm))
-  push!(trace_forms, (last( face_ranges), zero_moment))
+  trace_forms = (trace_forms..., (face_ranges[D], HDiv_facet_tr_norm))
+  trace_forms = (trace_forms..., (last( face_ranges), zero_moment))
   qorder = 2*get_order(field)+1
   FaceIntegralFormVector(p,qorder,trace_forms)
 end
 
 function get_normal_flux_forms(p::Polytope{D}, field) where D
   facet_range = get_dimrange(p,D-1)
-  fun_flux_form = Tuple[ (facet_range, HDiv_facet_flux)  ]
+  fun_flux_form = Tuple[ (facet_range, HDiv_facet_flux) ]
   qorder = get_order(field)+2
   FaceIntegralFormVector(p,qorder,fun_flux_form)
 end
diff --git a/test/ReferenceFEsTests/LagrangianDofBasesTests.jl b/test/ReferenceFEsTests/LagrangianDofBasesTests.jl
index 55c071467..b9189bf90 100644
--- a/test/ReferenceFEsTests/LagrangianDofBasesTests.jl
+++ b/test/ReferenceFEsTests/LagrangianDofBasesTests.jl
@@ -13,14 +13,14 @@ x4 = Point(1,1)
 x = [x1,x2,x3,x4]
 np = length(x)
 
-db = LagrangianDofBasis(Float64,x)
+let v = 3.0
 
+db = LagrangianDofBasis(Float64,x)
 @test db.nodes === x
 @test db.node_and_comp_to_dof == collect(1:np)
 @test db.dof_to_node == collect(1:np)
 @test db.dof_to_comp == fill(1,np)
 
-v = 3.0
 f = GenericField(x->v*x[1])
 fx = evaluate(f,x)
 test_dof_array(db,f,fx)
@@ -29,6 +29,12 @@ ndof = 8
 b = fill(f,ndof)
 bx = evaluate(b,x)
 test_dof_array(db,b,bx)
+end # let v
+
+let v = VectorValue(1,2)
+
+f = GenericField(x->v*x[1])
+dbf = [0, 1, 0, 1, 0, 2, 0, 2]
 
 T = VectorValue{2,Float64}
 db = LagrangianDofBasis(T,x)
@@ -36,11 +42,6 @@ db = LagrangianDofBasis(T,x)
 @test db.node_and_comp_to_dof == VectorValue{2,Int}[(1, 5), (2, 6), (3, 7), (4, 8)]
 @test db.dof_to_node == [1, 2, 3, 4, 1, 2, 3, 4]
 @test db.dof_to_comp == [1, 1, 1, 1, 2, 2, 2, 2]
-
-v = VectorValue(1,2)
-f = GenericField(x->v*x[1])
-dbf = [0, 1, 0, 1, 0, 2, 0, 2]
-
 test_dof_array(db,f,dbf)
 
 ndof = 8
@@ -58,11 +59,14 @@ db = LagrangianDofBasis(T,x)
 @test db.node_and_comp_to_dof == TensorValue{2,2,Int}[(1,5,9,13), (2,6,10,14), (3,7,11,15), (4,8,12,16)]
 @test db.dof_to_node == [1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4]
 @test db.dof_to_comp == [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4]
+end # let v
 
-v = TensorValue(1,2,3,4)
+let v = TensorValue(1,2,3,4)
 f = GenericField(x->v*x[1])
 dbf = [0, 1, 0, 1, 0, 2, 0, 2, 0, 3, 0, 3, 0, 4, 0, 4]
 
+T = TensorValue{2,2,Float64}
+db = LagrangianDofBasis(T,x)
 test_dof_array(db,f,dbf)
 
 ndof = 16
@@ -75,7 +79,9 @@ dbb = [
   0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4;]
 
 test_dof_array(db,b,dbb)
+end # let v
 
+let v = SymTensorValue(1,2,3)
 
 T = SymTensorValue{2,Float64}
 db = LagrangianDofBasis(T,x)
@@ -84,7 +90,6 @@ db = LagrangianDofBasis(T,x)
 @test db.dof_to_node == [1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4]
 @test db.dof_to_comp == [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3]
 
-v = SymTensorValue(1,2,3)
 f = GenericField(x->v*x[1])
 dbf = [0, 1, 0, 1, 0, 2, 0, 2, 0, 3, 0, 3]
 
@@ -99,7 +104,11 @@ dbb = [
   0 0 0 0 0 0 0 0 0 0 0 0; 3 3 3 3 3 3 3 3 3 3 3 3; 0 0 0 0 0 0 0 0 0 0 0 0; 3 3 3 3 3 3 3 3 3 3 3 3;]
 
 test_dof_array(db,b,dbb)
+end # let v
 
+let v = SymTracelessTensorValue(1,2)
+f = GenericField(x->v*x[1])
+dbf = [0, 1, 0, 1, 0, 2, 0, 2]
 
 T = SymTracelessTensorValue{2,Float64}
 db = LagrangianDofBasis(T,x)
@@ -108,10 +117,6 @@ db = LagrangianDofBasis(T,x)
 @test db.dof_to_node == [1, 2, 3, 4, 1, 2, 3, 4]
 @test db.dof_to_comp == [1, 1, 1, 1, 2, 2, 2, 2]
 
-v = SymTracelessTensorValue(1,2)
-f = GenericField(x->v*x[1])
-dbf = [0, 1, 0, 1, 0, 2, 0, 2]
-
 test_dof_array(db,f,dbf)
 
 ndof = 8
@@ -123,4 +128,6 @@ dbb = [
 
 test_dof_array(db,b,dbb)
 
+end # let v
+
 end # module
diff --git a/test/ReferenceFEsTests/NedelecRefFEsTests.jl b/test/ReferenceFEsTests/NedelecRefFEsTests.jl
index e1d50b8d9..82927431c 100644
--- a/test/ReferenceFEsTests/NedelecRefFEsTests.jl
+++ b/test/ReferenceFEsTests/NedelecRefFEsTests.jl
@@ -47,8 +47,7 @@ test_reference_fe(reffe)
 prebasis = get_prebasis(reffe)
 dof_basis = get_dof_basis(reffe)
 
-v = VectorValue(3.0,0.0)
-field = GenericField(x->v*x[1])
+field = GenericField(x->VectorValue(3.0,0.0)*x[1])
 
 cache = return_cache(dof_basis,field)
 r = evaluate!(cache, dof_basis, field)
@@ -76,8 +75,7 @@ test_reference_fe(reffe)
 prebasis = get_prebasis(reffe)
 dof_basis = get_dof_basis(reffe)
 
-v = VectorValue(3.0,1.0,-2.)
-field = GenericField(x->v*x[1]*x[2]*x[3])
+field = GenericField(x->VectorValue(3.0,1.0,-2.)*x[1]*x[2]*x[3])
 
 cache = return_cache(dof_basis,field)
 r = evaluate!(cache, dof_basis, field)
@@ -241,18 +239,18 @@ dof_basis = get_dof_basis(reffe)
 
 face_odofs_L2 = get_face_own_dofs(reffe,L2Conformity())
 
-@test face_odofs_L2 == [Int64[], Int64[], Int64[], Int64[],
-                     Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[],
+@test face_odofs_L2 == [Int[], Int[], Int[], Int[],
+                     Int[], Int[], Int[], Int[], Int[], Int[], Int[], Int[], Int[], Int[],
                     collect(1:20)]
 
 face_odofs    = get_face_own_dofs(reffe)
 face_cdofs    = get_face_dofs(reffe)
 
-@test face_odofs == [Int64[], Int64[], Int64[], Int64[],
+@test face_odofs == [Int[], Int[], Int[], Int[],
                     [1, 2], [3, 4], [5, 6], [7, 8], [9, 10], [11, 12], [13, 14], [15, 16], [17, 18], [19, 20],
-                    Int64[]]
+                    Int[]]
 
-@test face_cdofs == [Int64[], Int64[], Int64[], Int64[],
+@test face_cdofs == [Int[], Int[], Int[], Int[],
                      [1, 2], [3, 4], [5, 6], [7, 8], [9, 10], [11, 12],
                      [1, 2, 3, 4, 5, 6, 13, 14], [1, 2, 7, 8, 9, 10, 15, 16], [3, 4, 7, 8, 11, 12, 17, 18], [5, 6, 9, 10, 11, 12, 19, 20],
                      [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]]
@@ -271,14 +269,14 @@ dof_basis = get_dof_basis(reffe)
 
 face_odofs_L2 = get_face_own_dofs(reffe,L2Conformity())
 
-@test face_odofs_L2 == [Int64[], Int64[], Int64[], Int64[],
-                     Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[],
+@test face_odofs_L2 == [Int[], Int[], Int[], Int[],
+                     Int[], Int[], Int[], Int[], Int[], Int[], Int[], Int[], Int[], Int[],
                     collect(1:84)]
 
 face_odofs    = get_face_own_dofs(reffe)
 face_cdofs    = get_face_dofs(reffe)
 
-@test face_odofs == [Int64[], Int64[], Int64[], Int64[],
+@test face_odofs == [Int[], Int[], Int[], Int[],
   [1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16], [17, 18, 19, 20], [21, 22, 23, 24],
   [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36],  # facet 123
   [37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48],  # facet 124
diff --git a/test/ReferenceFEsTests/RaviartThomasRefFEsTests.jl b/test/ReferenceFEsTests/RaviartThomasRefFEsTests.jl
index b2f4e1e81..898082ef3 100644
--- a/test/ReferenceFEsTests/RaviartThomasRefFEsTests.jl
+++ b/test/ReferenceFEsTests/RaviartThomasRefFEsTests.jl
@@ -50,8 +50,7 @@ test_reference_fe(reffe)
 prebasis = get_prebasis(reffe)
 dof_basis = get_dof_basis(reffe)
 
-v = VectorValue(3.0,0.0)
-field = GenericField(x->v*x[1])
+field = GenericField(x->VectorValue(3.0,0.0)*x[1])
 
 cache = return_cache(dof_basis,field)
 r = evaluate!(cache, dof_basis, field)
@@ -85,8 +84,7 @@ reffe = RaviartThomasRefFE(et,p,order)
 prebasis = get_prebasis(reffe)
 dof_basis = get_dof_basis(reffe)
 
-v = VectorValue(3.0,0.0)
-field = GenericField(x->v*x[1])
+field = GenericField(x->VectorValue(3.0,0.0)*x[1])
 
 predofs = get_dof_basis(reffe).predofs
 nodes, nf_nodes, nf_moments =  get_nodes(predofs), get_face_nodes_dofs(predofs), get_face_moments(predofs)
@@ -128,8 +126,7 @@ test_reference_fe(reffe)
 prebasis = get_prebasis(reffe)
 dof_basis = get_dof_basis(reffe)
 
-v = VectorValue(0.0,3.0,0.0)
-field = GenericField(x->v)
+field = GenericField(x->VectorValue(0.0,3.0,0.0))
 
 cache = return_cache(dof_basis,field)
 r = evaluate!(cache, dof_basis, field)
diff --git a/test/TensorValuesTests/OperationsTests.jl b/test/TensorValuesTests/OperationsTests.jl
index 4fb325431..a94b7feec 100644
--- a/test/TensorValuesTests/OperationsTests.jl
+++ b/test/TensorValuesTests/OperationsTests.jl
@@ -1195,7 +1195,7 @@ c = HighOrderTensorValue{Tuple{3, 3, 3, 3}}([16058 57615 111407; 57615 167523 22
 a = SymFourthOrderTensorValue{4}(primes[1:100]...)
 b = SymFourthOrderTensorValue{4}(primes[101:200]...)
 #c = MultiValue(SArray(a)) ⋅² MultiValue(SArray(b))
-c = HighOrderTensorValue{Tuple{4, 4, 4, 4}, Int64}([188413 712136 1304296 2026734; 712136 2748622 3488732 4312968; 1304296 3488732 5129196 5935236; 2026734 4312968 5935236 6749138;;; 189575 717188 1313908 2041950; 717188 2769358 3515192 4345748; 1313908 3515192 5168440 5980840; 2041950 4345748 5980840 6800930;;; 190827 722090 1322958 2056060; 722090 2788564 3539610 4375934; 1322958 3539610 5204362 6022438; 2056060 4375934 6022438 6848208;;; 191791 725884 1330012 2067110; 725884 2803590 3558720 4399596; 1330012 3558720 5232580 6055128; 2067110 4399596 6055128 6885390;;;; 189575 717188 1313908 2041950; 717188 2769358 3515192 4345748; 1313908 3515192 5168440 5980840; 2041950 4345748 5980840 6800930;;; 193313 732352 1342192 2086054; 732352 2829214 3591684 4440152; 1342192 3591684 5281152 6111468; 2086054 4440152 6111468 6949414;;; 194655 738228 1353364 2103654; 738228 2853238 3622320 4478092; 1353364 3622320 5326516 6164152; 2103654 4478092 6164152 7009246;;; 196095 743898 1363934 2120148; 743898 2875580 3650802 4513262; 1363934 3650802 5368478 6212798; 2120148 4513262 6212798 7064468;;;; 190827 722090 1322958 2056060; 722090 2788564 3539610 4375934; 1322958 3539610 5204362 6022438; 2056060 4375934 6022438 6848208;;; 194655 738228 1353364 2103654; 738228 2853238 3622320 4478092; 1353364 3622320 5326516 6164152; 2103654 4478092 6164152 7009246;;; 197749 750626 1376534 2139800; 750626 2902328 3684882 4555438; 1376534 3684882 5418786 6271086; 2139800 4555438 6271086 7130732;;; 199085 756340 1387284 2156610; 756340 2925218 3714192 4591660; 1387284 3714192 5462008 6321248; 2156610 4591660 6321248 7187638;;;; 191791 725884 1330012 2067110; 725884 2803590 3558720 4399596; 1330012 3558720 5232580 6055128; 2067110 4399596 6055128 6885390;;; 196095 743898 1363934 2120148; 743898 2875580 3650802 4513262; 1363934 3650802 5368478 6212798; 2120148 4513262 6212798 7064468;;; 199085 756340 1387284 2156610; 756340 2925218 3714192 4591660; 1387284 3714192 5462008 6321248; 2156610 4591660 6321248 7187638;;; 200937 763194 1399646 2175840; 763194 2951440 3747226 4632710; 1399646 3747226 5510550 6377398; 2175840 4632710 6377398 7251424])
+c = HighOrderTensorValue{Tuple{4, 4, 4, 4}, Int}([188413 712136 1304296 2026734; 712136 2748622 3488732 4312968; 1304296 3488732 5129196 5935236; 2026734 4312968 5935236 6749138;;; 189575 717188 1313908 2041950; 717188 2769358 3515192 4345748; 1313908 3515192 5168440 5980840; 2041950 4345748 5980840 6800930;;; 190827 722090 1322958 2056060; 722090 2788564 3539610 4375934; 1322958 3539610 5204362 6022438; 2056060 4375934 6022438 6848208;;; 191791 725884 1330012 2067110; 725884 2803590 3558720 4399596; 1330012 3558720 5232580 6055128; 2067110 4399596 6055128 6885390;;;; 189575 717188 1313908 2041950; 717188 2769358 3515192 4345748; 1313908 3515192 5168440 5980840; 2041950 4345748 5980840 6800930;;; 193313 732352 1342192 2086054; 732352 2829214 3591684 4440152; 1342192 3591684 5281152 6111468; 2086054 4440152 6111468 6949414;;; 194655 738228 1353364 2103654; 738228 2853238 3622320 4478092; 1353364 3622320 5326516 6164152; 2103654 4478092 6164152 7009246;;; 196095 743898 1363934 2120148; 743898 2875580 3650802 4513262; 1363934 3650802 5368478 6212798; 2120148 4513262 6212798 7064468;;;; 190827 722090 1322958 2056060; 722090 2788564 3539610 4375934; 1322958 3539610 5204362 6022438; 2056060 4375934 6022438 6848208;;; 194655 738228 1353364 2103654; 738228 2853238 3622320 4478092; 1353364 3622320 5326516 6164152; 2103654 4478092 6164152 7009246;;; 197749 750626 1376534 2139800; 750626 2902328 3684882 4555438; 1376534 3684882 5418786 6271086; 2139800 4555438 6271086 7130732;;; 199085 756340 1387284 2156610; 756340 2925218 3714192 4591660; 1387284 3714192 5462008 6321248; 2156610 4591660 6321248 7187638;;;; 191791 725884 1330012 2067110; 725884 2803590 3558720 4399596; 1330012 3558720 5232580 6055128; 2067110 4399596 6055128 6885390;;; 196095 743898 1363934 2120148; 743898 2875580 3650802 4513262; 1363934 3650802 5368478 6212798; 2120148 4513262 6212798 7064468;;; 199085 756340 1387284 2156610; 756340 2925218 3714192 4591660; 1387284 3714192 5462008 6321248; 2156610 4591660 6321248 7187638;;; 200937 763194 1399646 2175840; 763194 2951440 3747226 4632710; 1399646 3747226 5510550 6377398; 2175840 4632710 6377398 7251424])
 @test SArray(c) == SArray(a ⋅² b)
 
 a = ThirdOrderTensorValue{3,3,3}(primes[1:27]...)