diff --git a/src/FESpaces/CurlConformingFESpaces.jl b/src/FESpaces/CurlConformingFESpaces.jl
index 0bf9e45cc..55b3d484a 100644
--- a/src/FESpaces/CurlConformingFESpaces.jl
+++ b/src/FESpaces/CurlConformingFESpaces.jl
@@ -45,7 +45,7 @@ function get_cell_shapefuns(
 
   cell_reffe_shapefuns = lazy_map(get_shapefuns,cell_reffe)
   cell_map = get_cell_map(Triangulation(model))
+  cell_Jt  = lazy_map(Broadcasting(∇),cell_map)
   k = ReferenceFEs.CoVariantPiolaMap()
-  lazy_map(k,cell_reffe_shapefuns,cell_map)
+  lazy_map(k,cell_reffe_shapefuns,cell_Jt)
 end
-
diff --git a/src/FESpaces/DivConformingFESpaces.jl b/src/FESpaces/DivConformingFESpaces.jl
index db08c8ba2..e72724247 100644
--- a/src/FESpaces/DivConformingFESpaces.jl
+++ b/src/FESpaces/DivConformingFESpaces.jl
@@ -47,9 +47,10 @@ function get_cell_shapefuns(
   sign_flip = get_sign_flip(model, cell_reffe)
 )
   cell_map = get_cell_map(get_grid(model))
+  cell_Jt  = lazy_map(Broadcasting(∇),cell_map)
   cell_shapefuns = lazy_map(get_shapefuns,cell_reffe)
   k = ContraVariantPiolaMap()
-  lazy_map(k,cell_shapefuns,cell_map,lazy_map(Broadcasting(constant_field),sign_flip))
+  lazy_map(k,cell_shapefuns,cell_Jt,lazy_map(Broadcasting(constant_field),sign_flip))
 end
 
 struct SignFlipMap{T} <: Map
diff --git a/src/ReferenceFEs/MomentBasedReferenceFEs.jl b/src/ReferenceFEs/MomentBasedReferenceFEs.jl
index 27e551928..72abb3af9 100644
--- a/src/ReferenceFEs/MomentBasedReferenceFEs.jl
+++ b/src/ReferenceFEs/MomentBasedReferenceFEs.jl
@@ -36,7 +36,7 @@ struct MomentBasedDofBasis{P,V} <: AbstractVector{Moment}
 end
 
 Base.size(a::MomentBasedDofBasis) = (num_dofs(a),)
-Base.axes(a::MomentBasedDofBasis) = (Base.OneTo(num_dofs),)
+Base.axes(a::MomentBasedDofBasis) = (Base.OneTo(num_dofs(a)),)
 Base.getindex(a::MomentBasedDofBasis,i::Integer) = Moment()
 Base.IndexStyle(::MomentBasedDofBasis) = IndexLinear()
 
diff --git a/src/ReferenceFEs/NedelecRefFEs.jl b/src/ReferenceFEs/NedelecRefFEs.jl
index 3faa07361..985984b92 100644
--- a/src/ReferenceFEs/NedelecRefFEs.jl
+++ b/src/ReferenceFEs/NedelecRefFEs.jl
@@ -115,43 +115,3 @@ function get_face_dofs(reffe::GenericRefFE{Nedelec,Dc}) where Dc
   end
   face_dofs
 end
-
-struct CoVariantPiolaMap <: Map end
-
-function evaluate!(
-  cache,
-  ::Broadcasting{typeof(∇)},
-  a::Fields.BroadcastOpFieldArray{CoVariantPiolaMap})
-  v, Jt = a.args
-  # Assuming J comes from an affine map
-  ∇v = Broadcasting(∇)(v)
-  k = CoVariantPiolaMap()
-  Broadcasting(Operation(k))(∇v,Jt)
-end
-
-function lazy_map(
-  ::Broadcasting{typeof(gradient)},
-  a::LazyArray{<:Fill{Broadcasting{Operation{CoVariantPiolaMap}}}})
-  v, Jt = a.args
-  ∇v = lazy_map(Broadcasting(∇),v)
-  k = CoVariantPiolaMap()
-  lazy_map(Broadcasting(Operation(k)),∇v,Jt)
-end
-
-function evaluate!(cache,::CoVariantPiolaMap,v::Number,Jt::Number)
-  v⋅transpose(inv(Jt)) # we multiply by the right side to compute the gradient correctly
-end
-
-function evaluate!(cache,k::CoVariantPiolaMap,v::AbstractVector{<:Field},phi::Field)
-  Jt = ∇(phi)
-  Broadcasting(Operation(k))(v,Jt)
-end
-
-function lazy_map(
-  k::CoVariantPiolaMap,
-  cell_ref_shapefuns::AbstractArray{<:AbstractArray{<:Field}},
-  cell_map::AbstractArray{<:Field})
-
-  cell_Jt = lazy_map(∇,cell_map)
-  lazy_map(Broadcasting(Operation(k)),cell_ref_shapefuns,cell_Jt)
-end
diff --git a/src/ReferenceFEs/Pullbacks.jl b/src/ReferenceFEs/Pullbacks.jl
index 9c6477238..d339d446b 100644
--- a/src/ReferenceFEs/Pullbacks.jl
+++ b/src/ReferenceFEs/Pullbacks.jl
@@ -11,7 +11,7 @@ where
 abstract type Pushforward <: Map end
 
 function Arrays.lazy_map(
-  k::Pushforward, ref_cell_fields, pf_args...
+  k::Pushforward, ref_cell_fields::AbstractArray, pf_args::AbstractArray...
 )
   lazy_map(Broadcasting(Operation(k)), ref_cell_fields, pf_args...)
 end
@@ -22,6 +22,12 @@ function Arrays.evaluate!(
   @abstractmethod
 end
 
+function evaluate!(
+  cache, k::Pushforward, f_ref::AbstractVector{<:Field}, args...
+)
+  Broadcasting(Operation(k))(f_ref,args...)
+end
+
 function Arrays.lazy_map(
   ::Broadcasting{typeof(gradient)}, a::LazyArray{<:Fill{Broadcasting{Operation{<:Pushforward}}}}
 )
@@ -30,6 +36,17 @@ function Arrays.lazy_map(
   return lazy_map(a.maps.value,cell_ref_gradient,args...)
 end
 
+function Arrays.evaluate!(
+  cache,
+  ::Broadcasting{typeof(∇)},
+  a::Fields.BroadcastOpFieldArray{<:Pushforward}
+)
+  v, Jt, sign_flip = a.args
+  ∇v = Broadcasting(∇)(v)
+  k = ContraVariantPiolaMap()
+  Broadcasting(Operation(k))(∇v,Jt,sign_flip)
+end
+
 # InversePushforward
 
 """
@@ -53,7 +70,7 @@ Arrays.inverse_map(pf::Pushforward) = InversePushforward(pf)
 Arrays.inverse_map(ipf::InversePushforward) = ipf.pushforward
 
 function Arrays.lazy_map(
-  k::InversePushforward, phys_cell_fields, pf_args...
+  k::InversePushforward, phys_cell_fields::AbstractArray, pf_args::AbstractArray...
 )
   lazy_map(Broadcasting(Operation(k)), phys_cell_fields, pf_args...)
 end
@@ -61,6 +78,7 @@ end
 function Arrays.return_cache(
   k::InversePushforward, v_phys::Number, args...
 )
+  mock_basis(::VectorValue{D,T}) where {D,T} = one(TensorValue{D,D,T})
   v_ref_basis = mock_basis(v_phys)
   pf_cache = return_cache(k.pushforward,v_ref_basis,args...)
   return v_ref_basis, pf_cache
@@ -71,10 +89,16 @@ function Arrays.evaluate!(
 )
   v_ref_basis, pf_cache = cache
   change = evaluate!(pf_cache,k.pushforward,v_ref_basis,args...)
-  return inv(change)⋅v_phys
+  return v_phys⋅inv(change)
+end
+
+function evaluate!(
+  cache, k::InversePushforward, f_phys::AbstractVector{<:Field}, args...
+)
+  Broadcasting(Operation(k))(f_phys,args...)
 end
 
-# Pushforward
+# Pullback
 
 """
     struct Pullback{PF <: Pushforward} <: Map end
@@ -96,7 +120,7 @@ struct Pullback{PF} <: Map
 end
 
 function Arrays.lazy_map(
-  ::typeof{evaluate},k::LazyArray{<:Fill{<:Pullback}},ref_cell_fields
+  ::typeof(evaluate),k::LazyArray{<:Fill{<:Pullback}},ref_cell_fields::AbstractArray
 )
   pb = k.maps.value
   phys_cell_dofs, pf_args = k.args
@@ -104,6 +128,12 @@ function Arrays.lazy_map(
   return lazy_map(evaluate,phys_cell_dofs,phys_cell_fields)
 end
 
+function evaluate!(
+  cache, k::PullBack, σ_phys::MomentBasedDofBasis, args...
+)
+  Broadcasting(Operation(k))(f_phys,args...)
+end
+
 # InversePullback
 
 """
@@ -129,7 +159,7 @@ Arrays.inverse_map(pb::Pullback) = InversePullback(pb.pushforward)
 Arrays.inverse_map(ipb::InversePullback) = Pullback(ipb.pushforward)
 
 function Arrays.lazy_map(
-  ::typeof{evaluate},k::LazyArray{<:Fill{<:InversePullback}},phys_cell_fields
+  ::typeof(evaluate),k::LazyArray{<:Fill{<:InversePullback}},phys_cell_fields::AbstractArray
 )
   pb = inverse_map(k.maps.value)
   ref_cell_dofs, pf_args = k.args
@@ -142,11 +172,20 @@ end
 struct ContraVariantPiolaMap <: Pushforward end
 
 function evaluate!(
-  cache,::ContraVariantPiolaMap,
-  v::Number,
-  Jt::Number,
-  sign_flip::Bool
+  cache, ::ContraVariantPiolaMap, v_ref::Number, Jt::Number, sign_flip::Bool
+)
+  sign  = (-1)^sign_flip
+  idetJ = 1. / meas(Jt)
+  return (sign*v_ref)⋅(idetJ*Jt)
+end
+
+# CoVariantPiolaMap
+
+struct CoVariantPiolaMap <: Pushforward end
+
+function evaluate!(
+  cache, ::CoVariantPiolaMap, v_ref::Number, Jt::Number
 )
-  idetJ = 1/meas(Jt)
-  ((-1)^sign_flip*v)⋅(idetJ*Jt)
+  # we right-multiply to compute the gradient correctly
+  return v_ref⋅transpose(inv(Jt))
 end
diff --git a/src/ReferenceFEs/RaviartThomasRefFEs.jl b/src/ReferenceFEs/RaviartThomasRefFEs.jl
index 29844461b..b507293b8 100644
--- a/src/ReferenceFEs/RaviartThomasRefFEs.jl
+++ b/src/ReferenceFEs/RaviartThomasRefFEs.jl
@@ -88,59 +88,3 @@ function compute_jacobi_basis(::Type{T},p::ExtrusionPolytope{D},orders) where {D
   terms = _monomial_terms(extrusion,orders)
   JacobiPolynomialBasis{D}(T,orders,terms)
 end
-
-# ContraVariantPiolaMap
-
-struct ContraVariantPiolaMap <: Map end
-
-function evaluate!(
-  cache,
-  ::Broadcasting{typeof(∇)},
-  a::Fields.BroadcastOpFieldArray{ContraVariantPiolaMap}
-)
-  v, Jt, sign_flip = a.args
-  ∇v = Broadcasting(∇)(v)
-  k = ContraVariantPiolaMap()
-  Broadcasting(Operation(k))(∇v,Jt,sign_flip)
-end
-
-function lazy_map(
-  ::Broadcasting{typeof(gradient)},
-  a::LazyArray{<:Fill{Broadcasting{Operation{ContraVariantPiolaMap}}}}
-)
-  v, Jt, sign_flip = a.args
-  ∇v = lazy_map(Broadcasting(∇),v)
-  k = ContraVariantPiolaMap()
-  lazy_map(Broadcasting(Operation(k)),∇v,Jt,sign_flip)
-end
-
-function lazy_map(
-  k::ContraVariantPiolaMap,
-  cell_ref_shapefuns::AbstractArray{<:AbstractArray{<:Field}},
-  cell_map::AbstractArray{<:Field},
-  sign_flip::AbstractArray{<:AbstractArray{<:Field}}
-)
-  cell_Jt = lazy_map(∇,cell_map)
-  lazy_map(Broadcasting(Operation(k)),cell_ref_shapefuns,cell_Jt,sign_flip)
-end
-
-function evaluate!(
-  cache,::ContraVariantPiolaMap,
-  v::Number,
-  Jt::Number,
-  sign_flip::Bool
-)
-  idetJ = 1/meas(Jt)
-  ((-1)^sign_flip*v)⋅(idetJ*Jt)
-end
-
-function evaluate!(
-  cache,
-  k::ContraVariantPiolaMap,
-  v::AbstractVector{<:Field},
-  phi::Field,
-  sign_flip::AbstractVector{<:Field}
-)
-  Jt = ∇(phi)
-  Broadcasting(Operation(k))(v,Jt,sign_flip)
-end
diff --git a/src/ReferenceFEs/ReferenceFEs.jl b/src/ReferenceFEs/ReferenceFEs.jl
index abfc4766d..513ce5ada 100644
--- a/src/ReferenceFEs/ReferenceFEs.jl
+++ b/src/ReferenceFEs/ReferenceFEs.jl
@@ -248,6 +248,8 @@ include("StrangQuadratures.jl")
 
 include("XiaoGimbutasQuadratures.jl")
 
+include("Pullbacks.jl")
+
 include("MomentBasedReferenceFEs.jl")
 
 include("RaviartThomasRefFEs.jl")
diff --git a/test/pullbacks.jl b/test/pullbacks.jl
new file mode 100644
index 000000000..a59f851b3
--- /dev/null
+++ b/test/pullbacks.jl
@@ -0,0 +1,43 @@
+
+using FillArrays
+using Gridap
+using Gridap.ReferenceFEs, Gridap.FESpaces
+
+using Gridap.ReferenceFEs: Pullback
+
+model = CartesianDiscreteModel((0,2,0,4),(2,2))
+Ω = Triangulation(model)
+dΩ = Measure(Ω,2)
+pts = CellData.get_data(get_cell_points(dΩ))
+
+reffe = RaviartThomasRefFE(Float64,QUAD,1)
+
+V = FESpace(model,reffe)
+
+u(x) = VectorValue(x[1], -x[2])
+uh = interpolate(u,V)
+
+φ_phys = get_fe_basis(V).cell_basis
+φ_ref  = φ_phys.args[1]
+Jt, sign = φ_phys.args[2:end]
+
+σ_phys = get_fe_dof_basis(V).cell_dof
+σ_ref  = Fill(get_dof_basis(reffe),num_cells(model))
+
+App = lazy_map(evaluate,σ_phys,φ_phys)[1]
+Arr = lazy_map(evaluate,σ_ref,φ_ref)[1]
+
+pf  = ContraVariantPiolaMap()
+
+ipf = inverse_map(pf)
+f_ref = lazy_map(ipf,φ_phys,Jt,sign)
+Brr = lazy_map(evaluate,σ_ref,f_ref)[1]
+Brr == Arr
+φ_ref_x = lazy_map(evaluate,φ_ref,pts)[1]
+f_ref_x = lazy_map(evaluate,f_ref,pts)[1]
+f_ref_x == φ_ref_x
+
+pb  = Pullback(pf)
+θ_ref = lazy_map(pb,σ_phys,Jt,sign)
+
+ipb = inverse_map(pb)