diff --git a/src/ReferenceFEs/AWRefFEs.jl b/src/ReferenceFEs/AWRefFEs.jl
index a2c4810e8..72a9a527f 100644
--- a/src/ReferenceFEs/AWRefFEs.jl
+++ b/src/ReferenceFEs/AWRefFEs.jl
@@ -57,11 +57,11 @@ function ArnoldWintherRefFE(::Type{T},p::Polytope,order::Integer) where T
 end
 
 function ReferenceFE(p::Polytope,::ArnoldWinther, order)
-  BDMRefFE(Float64,p,order)
+  ArnoldWintherRefFE(Float64,p,order)
 end
 
 function ReferenceFE(p::Polytope,::ArnoldWinther,::Type{T}, order) where T
-  BDMRefFE(T,p,order)
+  ArnoldWintherRefFE(T,p,order)
 end
 
 function Conformity(reffe::GenericRefFE{ArnoldWinther},sym::Symbol)
diff --git a/src/ReferenceFEs/HHJRefFEs.jl b/src/ReferenceFEs/HHJRefFEs.jl
new file mode 100644
index 000000000..9c0c691e5
--- /dev/null
+++ b/src/ReferenceFEs/HHJRefFEs.jl
@@ -0,0 +1,70 @@
+
+struct HellanHerrmannJhonson <: PushforwardRefFE end
+
+const hhj = HellanHerrmannJhonson()
+
+Pushforward(::Type{<:HellanHerrmannJhonson}) = DoubleContraVariantPiolaMap()
+
+"""
+    struct HellanHerrmannJhonson <: PushforwardRefFE end
+    HellanHerrmannJhonsonRefFE(::Type{T},p::Polytope,order::Integer) where T
+
+Hellan-Herrmann-Jhonson reference finite element.
+
+References:
+
+- `The Hellan-Herrmann-Johnson method with curved elements`, Arnold and Walker (2020)
+
+"""
+function HellanHerrmannJhonsonRefFE(::Type{T},p::Polytope,order::Integer) where T
+  @assert p == TRI "HellanHerrmannJhonson Reference FE only defined for TRIangles"
+  
+  VT = SymTensorValue{2,T}
+  prebasis = MonomialBasis{2}(VT,order,Polynomials._p_filter)
+  fb = MonomialBasis{D-1}(T,order,Polynomials._p_filter)
+  cb = MonomialBasis{2}(VT,order-1,Polynomials._p_filter)
+
+  function cmom(φ,μ,ds) # Cell and Node moment function: σ_K(φ,μ) = ∫(φ:μ)dK
+    Broadcasting(Operation(⊙))(φ,μ)
+  end
+  function fmom(φ,μ,ds) # Face moment function (normal) : σ_F(φ,μ) = ∫((n·φ·n)*μ)dF
+    n = get_facet_normal(ds)
+    φn = Broadcasting(Operation(⋅))(φ,n)
+    nφn = Broadcasting(Operation(⋅))(n,φn)
+    Broadcasting(Operation(*))(nφn,μ)
+  end
+
+  moments = Tuple[
+    (get_dimrange(p,1),fmom,fb), # Face moments
+    (get_dimrange(p,2),cmom,cb)  # Cell moments
+  ]
+
+  return MomentBasedReferenceFE(HellanHerrmannJhonson(),p,prebasis,moments,DivConformity())
+end
+
+function ReferenceFE(p::Polytope,::HellanHerrmannJhonson, order)
+  HellanHerrmannJhonsonRefFE(Float64,p,order)
+end
+
+function ReferenceFE(p::Polytope,::HellanHerrmannJhonson,::Type{T}, order) where T
+  HellanHerrmannJhonsonRefFE(T,p,order)
+end
+
+function Conformity(reffe::GenericRefFE{HellanHerrmannJhonson},sym::Symbol)
+  hdiv = (:Hdiv,:HDiv)
+  if sym == :L2
+    L2Conformity()
+  elseif sym in hdiv
+    DivConformity()
+  else
+    @unreachable """\n
+    It is not possible to use conformity = $sym on a HellanHerrmannJhonson reference FE.
+
+    Possible values of conformity for this reference fe are $((:L2, hdiv...)).
+      """
+  end
+end
+
+function get_face_own_dofs(reffe::GenericRefFE{HellanHerrmannJhonson}, conf::DivConformity)
+  get_face_dofs(reffe)
+end
diff --git a/src/ReferenceFEs/Pullbacks.jl b/src/ReferenceFEs/Pullbacks.jl
index bbdf1ddda..1a825b460 100644
--- a/src/ReferenceFEs/Pullbacks.jl
+++ b/src/ReferenceFEs/Pullbacks.jl
@@ -61,18 +61,6 @@ Represents the inverse of a pushforward map F*, defined as
 where 
   - V̂ is a function space on the reference cell K̂ and 
   - V is a function space on the physical cell K.
-
-# Note: 
-
-Given a pushforward F*, we provide a default implementation for the inverse (F*)^-1 that 
-is not optimal (and I think wrong for nonlinear geometries???). 
-
-For better performance, the user should overload the following methods for each 
-specific pushforward type:  
-
-- `return_cache(k::InversePushforward{PF}, v_phys::Number, args...)`
-- `evaluate!(cache, k::InversePushforward{PF}, v_phys::Number, args...)`
-
 """
 const InversePushforward{PF} = InverseMap{PF} where PF <: Pushforward
 
@@ -82,23 +70,6 @@ function Arrays.lazy_map(
   lazy_map(Broadcasting(Operation(k)), phys_cell_fields, pf_args...)
 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(inverse_map(k),v_ref_basis,args...)
-  return v_ref_basis, pf_cache
-end
-
-function Arrays.evaluate!(
-  cache, k::InversePushforward, v_phys::Number, args...
-)
-  v_ref_basis, pf_cache = cache
-  change = evaluate!(pf_cache,inverse_map(k),v_ref_basis,args...)
-  return v_phys⋅inv(change)
-end
-
 function evaluate!(
   cache, k::InversePushforward, f_phys::AbstractVector{<:Field}, args...
 )
@@ -185,12 +156,6 @@ function evaluate!(
   return v_ref ⋅ (idetJ * Jt)
 end
 
-function return_cache(
-  ::InversePushforward{ContraVariantPiolaMap}, v_phys::Number, Jt::Number
-)
-  nothing
-end
-
 function evaluate!(
   cache, ::InversePushforward{ContraVariantPiolaMap}, v_phys::Number, Jt::Number
 )
@@ -229,17 +194,46 @@ struct CoVariantPiolaMap <: Pushforward end
 function evaluate!(
   cache, ::CoVariantPiolaMap, v_ref::Number, Jt::Number
 )
-  return v_ref ⋅ transpose(inv(Jt))
+  return v_ref ⋅ transpose(pinvJt(Jt))
+end
+
+function evaluate!(
+  cache, ::InversePushforward{CoVariantPiolaMap}, v_phys::Number, Jt::Number
+)
+  return v_phys ⋅ transpose(Jt)
 end
 
-function return_cache(
-  ::InversePushforward{CoVariantPiolaMap}, v_phys::Number, Jt::Number
+# DoubleContraVariantPiolaMap
+
+struct DoubleContraVariantPiolaMap <: Pushforward end
+
+function evaluate!(
+  cache, ::DoubleContraVariantPiolaMap, v_ref::Number, Jt::Number
 )
-  return nothing
+  _Jt = meas(Jt) * Jt
+  return transpose(_Jt) ⋅ v_ref ⋅ _Jt
 end
 
 function evaluate!(
-  cache, ::InversePushforward{CoVariantPiolaMap}, v_phys::Number, Jt::Number
+  cache, ::InversePushforward{DoubleContraVariantPiolaMap}, v_phys::Number, Jt::Number
 )
-  return v_phys ⋅ transpose(Jt)
+  iJt = meas(Jt) * pinvJt(Jt)
+  return transpose(iJt) ⋅ v_phys ⋅ iJt
+end
+
+# DoubleCoVariantPiolaMap
+
+struct DoubleCoVariantPiolaMap <: Pushforward end
+
+function evaluate!(
+  cache, ::DoubleCoVariantPiolaMap, v_ref::Number, Jt::Number
+)
+  iJt = pinvJt(Jt)
+  return iJt ⋅ v_ref ⋅ transpose(iJt)
+end
+
+function evaluate!(
+  cache, ::InversePushforward{DoubleCoVariantPiolaMap}, v_phys::Number, Jt::Number
+)
+  return Jt ⋅ v_phys ⋅ transpose(Jt)
 end