diff --git a/src/ReferenceFEs/CRRefFEs.jl b/src/ReferenceFEs/CRRefFEs.jl
new file mode 100644
index 000000000..261fcd64b
--- /dev/null
+++ b/src/ReferenceFEs/CRRefFEs.jl
@@ -0,0 +1,59 @@
+struct CR <: ReferenceFEName end   
+const cr = CR()
+
+"""
+CRRefFE(::Type{et},p::Polytope,order::Integer) where et
+
+The `order` argument has the following meaning: the divergence of the  functions in this basis
+is in the P space of degree `order-1`.
+
+"""
+function CRRefFE(::Type{T},p::Polytope,order::Integer) where T 
+  D = num_dims(p)
+
+  if is_simplex(p) && order == 1
+    prebasis = MonomialBasis{D}(T,order,Polynomials._p_filter)
+    fb = MonomialBasis{D-1}(T,0,Polynomials._p_filter)
+  else
+    @notimplemented "CR Reference FE only available for simplices and lowest order"
+  end
+
+  function fmom(φ,μ,ds) # Face moment function : σ_F(φ,μ) = 1/|F| ( ∫((φ)*μ)dF )
+    D = num_dims(ds.cpoly)
+    facet_measure = get_facet_measure(ds.cpoly, D-1)  
+    facet_measure_1 = Gridap.Fields.ConstantField(1 / facet_measure[ds.face])
+    φμ = Broadcasting(Operation(⋅))(φ,μ)
+    Broadcasting(Operation(*))(φμ,facet_measure_1)
+  end
+
+  moments = Tuple[
+    (get_dimrange(p,D-1),fmom,fb), # Face moments
+  ]
+
+  return Gridap.ReferenceFEs.MomentBasedReferenceFE(CR(),p,prebasis,moments,L2Conformity())
+end
+
+
+function ReferenceFE(p::Polytope,::CR, order)
+  CRRefFE(Float64,p,order)
+end
+
+function ReferenceFE(p::Polytope,::CR,::Type{T}, order) where T
+  CRRefFE(T,p,order)
+end
+
+function Conformity(reffe::GenericRefFE{CR},sym::Symbol)
+  if sym == :L2
+    L2Conformity()
+  else
+    @unreachable """\n
+    It is not possible to use conformity = $sym on a CR reference FE.
+
+    Possible values of conformity for this reference fe are $((:L2, hdiv...)).
+      """
+  end
+end
+
+function get_face_own_dofs(reffe::GenericRefFE{CR}, conf::L2Conformity)
+  get_face_dofs(reffe)
+end
diff --git a/src/ReferenceFEs/MomentBasedReferenceFEs.jl b/src/ReferenceFEs/MomentBasedReferenceFEs.jl
index c032e4c89..7974016d2 100644
--- a/src/ReferenceFEs/MomentBasedReferenceFEs.jl
+++ b/src/ReferenceFEs/MomentBasedReferenceFEs.jl
@@ -166,6 +166,49 @@ function get_extension(m::FaceMeasure{Df,Dc}) where {Df,Dc}
   return ConstantField(TensorValue(hcat([vs[2]-vs[1]...],[vs[3]-vs[1]...])))
 end
 
+# TO DO: Bug in 3D, when n==2; n==4 and D==3
+function get_facet_measure(p::Polytope{D}, face::Int) where D
+  measures = Float64[]
+  facet_entities = get_face_coordinates(p)
+  for entity in facet_entities   
+    n = length(entity)
+    if n == 1 
+       push!(measures, 0.0)  # A point has zero measure
+    elseif n == 2
+        # Length of an edge
+        p1, p2 = entity
+        push!(measures, norm(p2-p1))
+    elseif n == 3 && D == 2
+      # Perimeter of the closed polygon
+      n = length(entity)
+      perimeter = 0.0
+      for i in 1:n
+          p1, p2 = entity[i], entity[mod1(i+1, n)]  # cyclic indices
+          perimeter += norm([p2[1] - p1[1], p2[2] - p1[2]])
+      end
+      push!(measures, perimeter)
+    elseif n == 3 && D == 3
+      # Area of a simplex 
+      p1, p2, p3 = entity
+      v1 = [p2[i] - p1[i] for i in 1:D]
+      v2 = [p3[i] - p1[i] for i in 1:D]
+      area = 0.5 * norm(cross(v1, v2))
+      push!(measures, area)
+    elseif n == 4 && D == 3
+      # Volume of a tetrahedron ( To do: Should be perimeter of the tetrahedron.)
+      p1, p2, p3, p4 = entity
+      v1 = [p2[i] - p1[i] for i in 1:D]
+      v2 = [p3[i] - p1[i] for i in 1:D]
+      v3 = [p4[i] - p1[i] for i in 1:D]
+      volume = abs(dot(v1, cross(v2, v3))) / 6
+      push!(measures, volume)
+    end
+
+  end
+  dim = get_dimranges(p)[face+1]
+  return measures[dim]
+end
+
 function Arrays.return_cache(
   σ::Function,                # σ(φ,μ,ds) -> Field/Array{Field}
   φ::AbstractArray{<:Field},  # φ: prebasis (defined on the cell)
diff --git a/src/ReferenceFEs/ReferenceFEs.jl b/src/ReferenceFEs/ReferenceFEs.jl
index 44cbc6bee..4799cb396 100644
--- a/src/ReferenceFEs/ReferenceFEs.jl
+++ b/src/ReferenceFEs/ReferenceFEs.jl
@@ -56,6 +56,7 @@ export get_dimranges
 export get_dimrange
 export get_vertex_coordinates
 export get_facet_normal
+export get_facet_measure
 export get_facet_orientations
 export get_edge_tangent
 export get_vertex_permutations
@@ -182,6 +183,7 @@ export BDMRefFE
 export NedelecRefFE
 export BezierRefFE
 export ModalC0RefFE
+export CRRefFE
 
 export Lagrangian
 export DivConforming
@@ -197,6 +199,7 @@ export bdm
 export nedelec
 export bezier
 export modalC0
+export cr
 
 export Quadrature
 export QuadratureName
@@ -250,6 +253,8 @@ include("BDMRefFEs.jl")
 
 include("NedelecRefFEs.jl")
 
+include("CRRefFEs.jl")
+
 include("MockDofs.jl")
 
 include("BezierRefFEs.jl")
diff --git a/test/ReferenceFEsTests/BDMRefFEsTests.jl b/test/ReferenceFEsTests/BDMRefFEsTests.jl
index f695d2c19..7b124ae18 100644
--- a/test/ReferenceFEsTests/BDMRefFEsTests.jl
+++ b/test/ReferenceFEsTests/BDMRefFEsTests.jl
@@ -1,4 +1,4 @@
-module BDMRefFEsTest
+# module BDMRefFEsTest
 
 using Test
 using Gridap.Polynomials
@@ -34,7 +34,7 @@ field = GenericField(x->v*x[1])
 
 cache = return_cache(dof_basis,field)
 r = evaluate!(cache, dof_basis, field)
-test_dof_array(dof_basis,field,r)
+@enter test_dof_array(dof_basis,field,r)
 
 cache = return_cache(dof_basis,prebasis)
 r = evaluate!(cache, dof_basis, prebasis)
@@ -98,4 +98,4 @@ reffe = ReferenceFE(TET,bdm,Float64,1)
 
 @test BDM() == bdm
 
-end # module
+# end # module