diff --git a/docs/src/ReferenceFEs.md b/docs/src/ReferenceFEs.md
index 86f9102b2..33a85105f 100644
--- a/docs/src/ReferenceFEs.md
+++ b/docs/src/ReferenceFEs.md
@@ -162,10 +162,10 @@ get_face_own_nodes_permutations(reffe::NodalReferenceFE,d::Integer)
 get_face_nodes(reffe::NodalReferenceFE,d::Integer)
 ```
 
-### GenericNodalRefFE
+### GenericNodalCartesianRefFE
 
 ```@docs
-GenericNodalRefFE
+GenericNodalCartesianRefFE
 ```
 
 ### Lagrangian reference elements
diff --git a/docs/src/TensorValues.md b/docs/src/TensorValues.md
index e47f09c4e..80e7678df 100644
--- a/docs/src/TensorValues.md
+++ b/docs/src/TensorValues.md
@@ -34,6 +34,6 @@ n_components
 inner
 outer
 meas
-symmetic_part
+symmetric_part
 ```
 
diff --git a/src/Arrays/Apply.jl b/src/Arrays/Apply.jl
index 2de0c12a2..1d7cb1eb4 100644
--- a/src/Arrays/Apply.jl
+++ b/src/Arrays/Apply.jl
@@ -118,7 +118,7 @@ end
 """
     apply_all(f::Tuple,a::AbstractArray...) -> Tuple
 
-Numerically equivalent to 
+Numerically equivalent to
 
     tuple( ( apply(fi, a...) for fi in f)... )
 
@@ -370,4 +370,3 @@ Base.IndexStyle(::Type{<:ArrayWithCounter{T,N,A}}) where {T,A,N} = IndexStyle(A)
 function reset_counter!(a::ArrayWithCounter)
   a.counter[:] .= 0
 end
-
diff --git a/src/Arrays/Tables.jl b/src/Arrays/Tables.jl
index 367fb434a..d9964317e 100644
--- a/src/Arrays/Tables.jl
+++ b/src/Arrays/Tables.jl
@@ -6,7 +6,7 @@
       ptrs::Vector{P}
     end
 
-Type representing a list of lists (i.e., a table) in 
+Type representing a list of lists (i.e., a table) in
 compressed format.
 """
 struct Table{T,P} <: AbstractVector{Vector{T}}
@@ -337,7 +337,7 @@ end
 """
     collect1d(a)
 
-Equivalent to 
+Equivalent to
 
     [a[i] for in 1:length(a)]
 """
@@ -442,5 +442,3 @@ function from_dict(::Type{Table{T,P}}, dict::Dict{Symbol,Any}) where {T,P}
   ptrs::Vector{P} = dict[:ptrs]
   Table(data,ptrs)
 end
-
-
diff --git a/src/FESpaces/CLagrangianFESpaces.jl b/src/FESpaces/CLagrangianFESpaces.jl
index ddc8187a5..810c0e4b5 100644
--- a/src/FESpaces/CLagrangianFESpaces.jl
+++ b/src/FESpaces/CLagrangianFESpaces.jl
@@ -9,7 +9,7 @@
     end
 """
 struct CLagrangianFESpace{S} <: SingleFieldFESpace
-  space::UnsconstrainedFESpace
+  space::UnconstrainedFESpace
   grid::Grid
   dof_to_node::Vector{Int}
   dof_to_comp::Vector{Int8}
@@ -19,7 +19,7 @@ struct CLagrangianFESpace{S} <: SingleFieldFESpace
       CLagrangianFESpace(::Type{T},grid::Grid) where T
   """
   function CLagrangianFESpace(::Type{T},grid::Grid) where T
-    space, dof_to_node, dof_to_comp, node_and_comp_to_dof = _generate_clargangian_fespace(T,grid)
+    space, dof_to_node, dof_to_comp, node_and_comp_to_dof = _generate_clagrangian_fespace(T,grid)
     S = eltype(node_and_comp_to_dof)
     new{S}(space,grid,dof_to_node,dof_to_comp,node_and_comp_to_dof)
   end
@@ -81,7 +81,7 @@ end
 
 # Helpers
 
-function _generate_clargangian_fespace(T,grid)
+function _generate_clagrangian_fespace(T,grid)
 
     grid_reffes = get_reffes(grid)
     function fun(reffe)
@@ -114,7 +114,7 @@ function _generate_clargangian_fespace(T,grid)
     dirichlet_cells = Int[]
     ntags = 0
 
-    space = UnsconstrainedFESpace(
+    space = UnconstrainedFESpace(
       nfree,
       ndirichlet,
       cell_dofs,
@@ -235,4 +235,3 @@ function _fill_cell_dofs_clagrangian_fespace!(
   end
 
 end
-
diff --git a/src/FESpaces/CellBases.jl b/src/FESpaces/CellBases.jl
index c25da70e4..7bdf2937f 100644
--- a/src/FESpaces/CellBases.jl
+++ b/src/FESpaces/CellBases.jl
@@ -41,10 +41,20 @@ end
 
 # Define how the metadata is preserved
 
+function change_ref_style(cf::CellBasis)
+  ref_sty = RefStyle(cf)
+  bool = !get_val_parameter(ref_sty)
+  new_sty = Val{bool}()
+  trial_style = TrialStyle(cf)
+  ar = get_array(cf)
+  cm = get_cell_map(cf)
+  GenericCellBasis(trial_style,ar,cm,new_sty)
+end
+
 function similar_object(cf::CellBasis,array::AbstractArray)
   cm = get_cell_map(cf)
   trial_style = TrialStyle(cf)
-  GenericCellBasis(trial_style,array,cm)
+  GenericCellBasis(trial_style,array,cm,RefStyle(cf))
 end
 
 function similar_object(a::CellBasis,b::CellField,v::AbstractArray)
@@ -73,7 +83,8 @@ end
 
 function _similar_cell_basis_test_trial(v,a,b)
   cm = get_cell_map(a)
-  GenericCellMatrixField(v,cm)
+  @assert is_in_ref_space(a) == is_in_ref_space(b)
+  GenericCellMatrixField(v,cm,RefStyle(a))
 end
 
 # Operations
@@ -114,13 +125,13 @@ end
 
 function operate(op,cf1::CellBasis,object)
   cm = get_cell_map(cf1)
-  cf2 = convert_to_cell_field(object,cm)
+  cf2 = convert_to_cell_field(object,cm,RefStyle(cf1))
   operate(op,cf1,cf2)
 end
 
 function operate(op,object,cf2::CellBasis)
   cm = get_cell_map(cf2)
-  cf1 = convert_to_cell_field(object,cm)
+  cf1 = convert_to_cell_field(object,cm,RefStyle(cf2))
   operate(op,cf1,cf2)
 end
 
@@ -191,17 +202,16 @@ end
 
 """
 """
-struct GenericCellBasis{T} <: CellBasis
+struct GenericCellBasis{T,R} <: CellBasis
   trial_style::Val{T}
   array
   cell_map
+  ref_trait::Val{R}
 end
 
-"""
-"""
-function GenericCellBasis(array::AbstractArray,cell_map::AbstractArray)
-  trial_style = Val{false}()
-  GenericCellBasis(trial_style,array,cell_map)
+function GenericCellBasis(trial_style::Val{T},array::AbstractArray,cell_map::AbstractArray) where T
+  ref_trait = Val{true}()
+  GenericCellBasis(trial_style,array,cell_map,ref_trait)
 end
 
 get_array(a::GenericCellBasis) = a.array
@@ -212,7 +222,9 @@ function TrialStyle(::Type{<:GenericCellBasis{T}}) where T
   Val{T}()
 end
 
-# CellMatrixField
+function RefStyle(::Type{<:GenericCellBasis{T,R}}) where {T,R}
+  Val{R}()
+end
 
 """
 """
@@ -226,9 +238,18 @@ end
 
 # Define how the metadata is preserved
 
+function change_ref_style(cf::CellMatrixField)
+  ref_sty = RefStyle(cf)
+  bool = !get_val_parameter(ref_sty)
+  new_sty = Val{bool}()
+  ar = get_array(cf)
+  cm = get_cell_map(cf)
+  GenericCellMatrixField(ar,cm,new_sty)
+end
+
 function similar_object(cf::CellMatrixField,a::AbstractArray)
   cm = get_cell_map(cf)
-  GenericCellMatrixField(a,cm)
+  GenericCellMatrixField(a,cm,RefStyle(cf))
 end
 
 function similar_object(a::CellMatrixField,b::CellField,v::AbstractArray)
@@ -274,25 +295,29 @@ end
 
 function operate(op,a::CellMatrixField,b)
   cm = get_cell_map(a)
-  _b = convert_to_cell_field(b,cm)
+  _b = convert_to_cell_field(b,cm,RefStyle(a))
   operate(op,a,_b)
 end
 
 function operate(op,a,b::CellMatrixField)
   cm = get_cell_map(b)
-  _a = convert_to_cell_field(a,cm)
+  _a = convert_to_cell_field(a,cm,RefStyle(b))
   operate(op,_a,b)
 end
 
+
 # Concrete CellMatrixField
 
 """
 """
-struct GenericCellMatrixField{A,B} <: CellMatrixField
+struct GenericCellMatrixField{A,B,R} <: CellMatrixField
   array::A
   cell_map::B
+  ref_style::Val{R}
 end
 
+RefStyle(::Type{GenericCellMatrixField{A,B,R}}) where {A,B,R} = Val{R}()
+
 get_array(a::GenericCellMatrixField) = a.array
 
 get_cell_map(a::GenericCellMatrixField) = a.cell_map
@@ -300,7 +325,8 @@ get_cell_map(a::GenericCellMatrixField) = a.cell_map
 # Restrictions
 
 function restrict(cf::CellBasis,trian::Triangulation)
-  a = get_array(cf)
+  _cf = to_ref_space(cf)
+  a = get_array(_cf)
   r = restrict(a,trian)
   trial_style = TrialStyle(cf)
   _restrict_cell_basis(trial_style,r,trian)
@@ -398,13 +424,13 @@ end
 
 function operate(op,cf1::SkeletonCellBasis,object)
   cm = get_cell_map(cf1)
-  cf2 = convert_to_cell_field(object,cm)
+  cf2 = convert_to_cell_field(object,cm,RefStyle(cf1.left))
   operate(op,cf1,cf2)
 end
 
 function operate(op,object,cf2::SkeletonCellBasis)
   cm = get_cell_map(cf2)
-  cf1 = convert_to_cell_field(object,cm)
+  cf1 = convert_to_cell_field(object,cm,RefStyle(cf2.left))
   operate(op,cf1,cf2)
 end
 
@@ -450,13 +476,13 @@ end
 
 function operate(op,cf1::ReducedSkeletonCellBasis,object)
   cm = get_cell_map(cf1)
-  cf2 = convert_to_cell_field(object,cm)
+  cf2 = convert_to_cell_field(object,cm,RefStyle(cf1.left))
   operate(op,cf1,cf2)
 end
 
 function operate(op,object,cf2::ReducedSkeletonCellBasis)
   cm = get_cell_map(cf2)
-  cf1 = convert_to_cell_field(object,cm)
+  cf1 = convert_to_cell_field(object,cm,RefStyle(cf2.left))
   operate(op,cf1,cf2)
 end
 
@@ -533,13 +559,13 @@ end
 
 function operate(op,object,cf2::SkeletonCellMatrixField)
   cm = get_cell_map(cf2)
-  cf1 = convert_to_cell_field(object,cm)
+  cf1 = convert_to_cell_field(object,cm,RefStyle(cf2.ll))
   operate(op,cf1,cf2)
 end
 
 function operate(op,cf1::SkeletonCellMatrixField,object)
   cm = get_cell_map(cf1)
-  cf2 = convert_to_cell_field(object,cm)
+  cf2 = convert_to_cell_field(object,cm,RefStyle(cf1.ll))
   operate(op,cf1,cf2)
 end
 
@@ -645,4 +671,3 @@ end
   add_to_array!(vecd,vec)
   (mat, vecd)
 end
-
diff --git a/src/FESpaces/CellDofBases.jl b/src/FESpaces/CellDofBases.jl
new file mode 100644
index 000000000..9a4118f98
--- /dev/null
+++ b/src/FESpaces/CellDofBases.jl
@@ -0,0 +1,55 @@
+"""
+   abstract type CellDofBasis end
+
+Abstract type that represents a cell array of `Dof`. The main motivation for
+its definition is to provide a trait that informs whether the `Dof` entries are
+defined for functions in the reference or physical space
+"""
+abstract type CellDofBasis end
+
+RefStyle(::Type{<:CellDofBasis}) = @abstractmethod
+get_array(::CellDofBasis) = @abstractmethod
+
+function test_cell_dof_basis(cf::CellDofBasis,f::CellFieldLike)
+  ar = get_array(cf)
+  @test isa(ar,AbstractArray)
+  a = evaluate(cf,f)
+  _ = collect(a)
+  rs = RefStyle(cf)
+  @test isa(get_val_parameter(rs),Bool)
+end
+
+"""
+    evaluate(dof_array::CellDofBasis,field_array::AbstractArray)
+
+Evaluates the `CellDofBasis` for the `Field` objects
+at the array `field` element by element.
+
+The result is numerically equivalent to
+
+    map(evaluate_dof, dof_array.array, field_array)
+
+but it is described with a more memory-friendly lazy type.
+"""
+function evaluate(cell_dofs::CellDofBasis,cell_field::CellFieldLike)
+ _evaluate_cell_dofs(cell_dofs,cell_field,RefStyle(cell_dofs))
+end
+
+function  _evaluate_cell_dofs(cell_dofs,cell_field,ref_trait::Val{true})
+  evaluate_dof_array(get_array(cell_dofs),get_array(to_ref_space(cell_field)))
+end
+
+function  _evaluate_cell_dofs(cell_dofs,cell_field,ref_trait::Val{false})
+  evaluate_dof_array(get_array(cell_dofs),get_array(to_physical_space(cell_field)))
+end
+
+"""
+"""
+struct GenericCellDofBasis{R} <: CellDofBasis
+   ref_trait::Val{R}
+   array::AbstractArray{<:Dof}
+end
+
+RefStyle(::Type{<:GenericCellDofBasis{R}}) where R = Val{R}()
+
+get_array(a::GenericCellDofBasis) = a.array
diff --git a/src/FESpaces/ConformingFESpaces.jl b/src/FESpaces/ConformingFESpaces.jl
index 0f0e28fd4..3b2434cdd 100644
--- a/src/FESpaces/ConformingFESpaces.jl
+++ b/src/FESpaces/ConformingFESpaces.jl
@@ -5,12 +5,13 @@ function GradConformingFESpace(
   reffes::Vector{<:ReferenceFE},
   model::DiscreteModel,
   dirichlet_tags,
-  dirichlet_components=nothing)
+  dirichlet_components=nothing,
+  is_ref=true)
 
-  face_labeing = get_face_labeling(model)
+  face_labeling = get_face_labeling(model)
 
   GradConformingFESpace(
-    reffes,model,face_labeing,dirichlet_tags,dirichlet_components)
+    reffes,model,face_labeling,dirichlet_tags,dirichlet_components,is_ref)
 
 end
 
@@ -19,11 +20,12 @@ function GradConformingFESpace(
   model::DiscreteModel,
   order::Integer,
   dirichlet_tags,
-  dirichlet_components=nothing) where T
+  dirichlet_components=nothing,
+  is_ref=true) where T
 
-  face_labeing = get_face_labeling(model)
+  face_labeling = get_face_labeling(model)
 
-  GradConformingFESpace(T,model,order,face_labeing,dirichlet_tags,dirichlet_components)
+  GradConformingFESpace(T,model,order,face_labeling,dirichlet_tags,dirichlet_components,is_ref)
 
 end
 
@@ -31,9 +33,10 @@ function GradConformingFESpace(
   ::Type{T},
   model::DiscreteModel,
   order::Integer,
-  face_labeing::FaceLabeling,
+  face_labeling::FaceLabeling,
   dirichlet_tags,
-  dirichlet_components=nothing) where T
+  dirichlet_components=nothing,
+  is_ref=true) where T
 
   grid_topology = get_grid_topology(model)
   polytopes = get_polytopes(grid_topology)
@@ -41,7 +44,7 @@ function GradConformingFESpace(
   reffes = [ LagrangianRefFE(T,p,order) for p in polytopes ]
 
   GradConformingFESpace(
-    reffes,model,face_labeing,dirichlet_tags,dirichlet_components)
+    reffes,model,face_labeling,dirichlet_tags,dirichlet_components,is_ref)
 
 end
 
@@ -50,16 +53,17 @@ end
 function GradConformingFESpace(
   reffes::Vector{<:LagrangianRefFE},
   model::DiscreteModel,
-  face_labeing::FaceLabeling,
+  face_labeling::FaceLabeling,
   dirichlet_tags,
-  dirichlet_components=nothing)
+  dirichlet_components=nothing,
+  is_ref=true)
 
   grid_topology = get_grid_topology(model)
 
   _dirichlet_components = _convert_dirichlet_components(dirichlet_tags,dirichlet_components)
 
   cell_dofs, nfree, ndirichlet, dirichlet_dof_tag, dirichlet_cells = compute_conforming_cell_dofs(
-    reffes,grid_topology,face_labeing,dirichlet_tags,_dirichlet_components)
+    reffes,grid_topology,face_labeling,dirichlet_tags,_dirichlet_components)
 
   ntags = length(dirichlet_tags)
 
@@ -67,9 +71,13 @@ function GradConformingFESpace(
   cell_to_ctype = get_cell_type(grid_topology)
   cell_map = get_cell_map(grid)
 
-  cell_shapefuns, cell_dof_basis = compute_cell_space(reffes, cell_to_ctype, cell_map)
+  if is_ref
+    cell_shapefuns, cell_dof_basis = compute_cell_space(reffes, cell_to_ctype, cell_map)
+  else
+    cell_shapefuns, cell_dof_basis = compute_cell_space_physical(reffes, cell_to_ctype, cell_map)
+  end
 
-  UnsconstrainedFESpace(
+  UnconstrainedFESpace(
     nfree,
     ndirichlet,
     cell_dofs,
@@ -88,25 +96,94 @@ function compute_cell_space(reffes, cell_to_ctype, cell_map)
 
   dof_basis = map(get_dof_basis,reffes)
   cell_dof_basis = CompressedArray(dof_basis,cell_to_ctype)
+  cell_dof_basis = GenericCellDofBasis(Val{true}(),cell_dof_basis)
 
   shapefuns =  map(get_shapefuns,reffes)
   refshapefuns = CompressedArray(shapefuns,cell_to_ctype)
   cell_shapefuns = attachmap(refshapefuns,cell_map)
+  cell_shapefuns = GenericCellBasis(Val{false}(),cell_shapefuns,cell_map,Val{true}())
+
+  (cell_shapefuns, cell_dof_basis)
+end
+
+"""
+It creates the cell-wise DOF basis and shape functions, when the DOFs
+are evaluated at the physical space. The DOFs (moments) for the prebasis
+are assumed to be computable at a reference FE space.
+"""
+function compute_cell_space_physical(reffes, cell_to_ctype, cell_map)
+# E.g., if one has to implement $\int_F q ϕ_h(x)$ for $q \in P(F)$, we can
+# assume that it can be written as $\sum_{x_{gp} \in Q}_{\hat{F}}
+# \hat{q}(\tilde{x}_{gp}) ϕ(x_{gp})$ for $\hat{q} \in P(\hat{F})$.
+
+  dof_bases = map(get_dof_basis,reffes)
+
+  cell_dof_basis = _cell_dof_basis_physical_space(dof_bases,cell_to_ctype,cell_map)
+  ref_style = Val{false}()
+  cell_dof_basis = GenericCellDofBasis(ref_style,cell_dof_basis)
+
+  prebasis =  map(get_prebasis,reffes)
+  cell_prebasis = CompressedArray(prebasis,cell_to_ctype)
+  cell_prebasis = GenericCellBasis(Val{false}(),cell_prebasis,cell_map,Val{false}())
+  cell_shapefuns = _cell_shape_functions_physical_space(cell_prebasis,cell_dof_basis,cell_map)
 
   (cell_shapefuns, cell_dof_basis)
 end
 
+function _cell_dof_basis_physical_space(dof_bases,cell_to_ctype,cell_map)
+  @notimplemented
+end
+
+function _cell_dof_basis_physical_space(
+  dof_bases::Vector{MomentBasedDofBasis{T,V}},
+  cell_to_ctype,cell_map) where {T,V}
+
+  ctype_to_refnodes = map(get_nodes,dof_bases)
+  cell_to_refnodes = CompressedArray(ctype_to_refnodes,cell_to_ctype)
+  cell_physnodes = evaluate(cell_map,cell_to_refnodes)
+
+  # Not efficient, create a Kernel
+  ct_face_moments = map(get_face_moments,dof_bases)
+  c_face_moments = CompressedArray(ct_face_moments,cell_to_ctype)
+  ct_face_nodes_dofs = map(get_face_nodes_dofs,dof_bases)
+  c_face_nodes_dofs = CompressedArray(ct_face_nodes_dofs,cell_to_ctype)
+  cell_dof_basis = apply( (n,m,nd) -> MomentBasedDofBasis(n,m,nd),
+                          cell_physnodes, c_face_moments, c_face_nodes_dofs)
+end
+
+function _cell_dof_basis_physical_space(
+  dof_bases::Vector{LagrangianDofBasis{T,V}},
+  cell_to_ctype,cell_map) where {T,V}
+
+  # Create new dof_basis with nodes in the physical space
+  # ctype_to_refnodes= map(get_node_coordinates,reffes)
+  ctype_to_refnodes= map(get_nodes,dof_bases)
+  cell_to_refnodes = CompressedArray(ctype_to_refnodes,cell_to_ctype)
+  cell_physnodes = evaluate(cell_map,cell_to_refnodes)
+  cell_dof_basis = apply( nodes -> LagrangianDofBasis(V,nodes), cell_physnodes )
+end
+
+function _cell_shape_functions_physical_space(cell_prebasis,cell_dof_basis,cell_map)
+  cell_matrix = evaluate(cell_dof_basis,cell_prebasis)
+  cell_matrix_inv = apply(inv,cell_matrix)
+  cell_shapefuns_phys = apply(change_basis,get_array(cell_prebasis),cell_matrix_inv)
+  ref_style = Val{false}()
+  trial_style = Val{false}()
+  GenericCellBasis(trial_style,cell_shapefuns_phys,cell_map,ref_style)
+   # @santiagobadia : better implementation in the future...
+end
+
 """
   compute_conforming_cell_dofs(
     reffes,
     grid_topology,
-    face_labeing,
+    face_labeling,
     dirichlet_tags)
 
   compute_conforming_cell_dofs(
     reffes,
     grid_topology,
-    face_labeing,
+    face_labeling,
     dirichlet_tags,
     dirichlet_components)
 
@@ -121,7 +198,7 @@ If `dirichlet_components`  is given, then `get_dof_to_comp` has to be defined
 for the reference elements in `reffes`.
 """
 function compute_conforming_cell_dofs(
-  reffes,grid_topology,face_labeing,dirichlet_tags,dirichlet_components=nothing)
+  reffes,grid_topology,face_labeling,dirichlet_tags,dirichlet_components=nothing)
 
   D = num_cell_dims(grid_topology)
   n_faces = num_faces(grid_topology)
@@ -140,7 +217,7 @@ function compute_conforming_cell_dofs(
      d_to_offset,
      d_to_ctype_to_ldface_to_own_ldofs)
 
-  d_to_dface_to_tag = [ get_face_tag_index(face_labeing,dirichlet_tags,d)  for d in 0:D]
+  d_to_dface_to_tag = [ get_face_tag_index(face_labeling,dirichlet_tags,d)  for d in 0:D]
   cell_to_faces = Table(get_cell_faces(grid_topology))
 
   nfree, ndiri, diri_dof_tag = _split_face_own_dofs_into_free_and_dirichlet!(
@@ -526,4 +603,3 @@ end
 ###  getindex!(cache,a,cell)
 ###end
 ###
-
diff --git a/src/FESpaces/CurlConformingFESpaces.jl b/src/FESpaces/CurlConformingFESpaces.jl
index d1461a1c2..167ac0b07 100644
--- a/src/FESpaces/CurlConformingFESpaces.jl
+++ b/src/FESpaces/CurlConformingFESpaces.jl
@@ -4,19 +4,20 @@
     CurlConformingFESpace(
       reffes::Vector{<:ReferenceFE},
       model::DiscreteModel,
-      face_labeing::FaceLabeling,
+      face_labeling::FaceLabeling,
       dirichlet_tags)
 """
 function CurlConformingFESpace(
   reffes::Vector{<:ReferenceFE},
   model::DiscreteModel,
-  face_labeing::FaceLabeling,
-  dirichlet_tags)
+  face_labeling::FaceLabeling,
+  dirichlet_tags,
+  is_ref)
 
   grid_topology = get_grid_topology(model)
 
   cell_dofs, nfree, ndirichlet, dirichlet_dof_tag, dirichlet_cells = compute_conforming_cell_dofs(
-    reffes,grid_topology,face_labeing,dirichlet_tags)
+    reffes,grid_topology,face_labeling,dirichlet_tags)
 
   ntags = length(dirichlet_tags)
 
@@ -24,9 +25,14 @@ function CurlConformingFESpace(
   cell_to_ctype = get_cell_type(grid_topology)
   cell_map = get_cell_map(grid)
 
-  cell_shapefuns, cell_dof_basis = _compute_hcurl_cell_space(reffes, cell_to_ctype, cell_map)
+  # cell_shapefuns, cell_dof_basis = _compute_hcurl_cell_space(reffes, cell_to_ctype, cell_map)
+  if is_ref
+    cell_shapefuns, cell_dof_basis = compute_cell_space(reffes,cell_to_ctype,cell_map)
+  else
+    cell_shapefuns, cell_dof_basis = compute_cell_space_physical(reffes,cell_to_ctype,cell_map)
+  end
 
-  UnsconstrainedFESpace(
+  UnconstrainedFESpace(
     nfree,
     ndirichlet,
     cell_dofs,
@@ -41,5 +47,5 @@ end
 
 function _compute_hcurl_cell_space(reffes, cell_to_ctype, cell_map)
   #TODO: fine for structured hex meshes, but not otherwise
-  compute_cell_space(reffes,cell_to_ctype,cell_map)
+  # compute_cell_space(reffes,cell_to_ctype,cell_map)
 end
diff --git a/src/FESpaces/DiscontinuousFESpaces.jl b/src/FESpaces/DiscontinuousFESpaces.jl
index 15677f5af..d0c2565f1 100644
--- a/src/FESpaces/DiscontinuousFESpaces.jl
+++ b/src/FESpaces/DiscontinuousFESpaces.jl
@@ -2,7 +2,7 @@
 """
     DiscontinuousFESpace(reffes::Vector{<:ReferenceFE}, trian::Triangulation)
 """
-function DiscontinuousFESpace(reffes::Vector{<:ReferenceFE}, trian::Triangulation)
+function DiscontinuousFESpace(reffes::Vector{<:ReferenceFE}, trian::Triangulation, is_ref=true)
 
   cell_to_ctype = get_cell_type(trian)
   cell_map = get_cell_map(trian)
@@ -14,9 +14,13 @@ function DiscontinuousFESpace(reffes::Vector{<:ReferenceFE}, trian::Triangulatio
   dirichlet_cells = Int[]
   ntags = 0
 
-  cell_shapefuns, cell_dof_basis = compute_cell_space(reffes, cell_to_ctype, cell_map)
+  if is_ref
+    cell_shapefuns, cell_dof_basis = compute_cell_space(reffes, cell_to_ctype, cell_map)
+  else
+    cell_shapefuns, cell_dof_basis = compute_cell_space_physical(reffes, cell_to_ctype, cell_map)
+  end
 
-  UnsconstrainedFESpace(
+  UnconstrainedFESpace(
     nfree,
     ndirichlet,
     cell_dofs,
diff --git a/src/FESpaces/DivConformingFESpaces.jl b/src/FESpaces/DivConformingFESpaces.jl
index e93351c2a..8517cdcdd 100644
--- a/src/FESpaces/DivConformingFESpaces.jl
+++ b/src/FESpaces/DivConformingFESpaces.jl
@@ -3,19 +3,20 @@
     DivConformingFESpace(
       reffes::Vector{<:ReferenceFE},
       model::DiscreteModel,
-      face_labeing::FaceLabeling,
+      face_labeling::FaceLabeling,
       dirichlet_tags)
 """
 function DivConformingFESpace(
   reffes::Vector{<:ReferenceFE},
   model::DiscreteModel,
-  face_labeing::FaceLabeling,
-  dirichlet_tags)
+  face_labeling::FaceLabeling,
+  dirichlet_tags,
+  is_ref)
 
   grid_topology = get_grid_topology(model)
 
   cell_dofs, nfree, ndirichlet, dirichlet_dof_tag, dirichlet_cells = compute_conforming_cell_dofs(
-    reffes,grid_topology,face_labeing,dirichlet_tags)
+    reffes,grid_topology,face_labeling,dirichlet_tags)
 
   ntags = length(dirichlet_tags)
 
@@ -23,9 +24,15 @@ function DivConformingFESpace(
   cell_to_ctype = get_cell_type(grid_topology)
   cell_map = get_cell_map(grid)
 
-  cell_shapefuns, cell_dof_basis = _compute_hdiv_cell_space(reffes, cell_to_ctype, cell_map)
+  # cell_shapefuns, cell_dof_basis = _compute_hdiv_cell_space(reffes, cell_to_ctype, cell_map)
+  if is_ref
+    cell_shapefuns, cell_dof_basis = compute_cell_space(reffes,cell_to_ctype,cell_map)
+  else
+    cell_shapefuns, cell_dof_basis = compute_cell_space_physical(reffes,cell_to_ctype,cell_map)
+  end
 
-  UnsconstrainedFESpace(
+
+  UnconstrainedFESpace(
     nfree,
     ndirichlet,
     cell_dofs,
@@ -40,5 +47,5 @@ end
 
 function _compute_hdiv_cell_space(reffes, cell_to_ctype, cell_map)
   #TODO: fine for structured hex meshes, but not otherwise
-  compute_cell_space(reffes,cell_to_ctype,cell_map)
+  # compute_cell_space(reffes,cell_to_ctype,cell_map)
 end
diff --git a/src/FESpaces/ExtendedFESpaces.jl b/src/FESpaces/ExtendedFESpaces.jl
index c7b5c49e3..29b6583a4 100644
--- a/src/FESpaces/ExtendedFESpaces.jl
+++ b/src/FESpaces/ExtendedFESpaces.jl
@@ -180,23 +180,27 @@ function get_cell_basis(f::ExtendedFESpace)
 
   cm = get_cell_map(f.trian.oldtrian)
   trial_style = TrialStyle(cell_basis)
-  GenericCellBasis(trial_style,array,cm)
+  GenericCellBasis(trial_style,array,cm,RefStyle(cell_basis))
 end
 
 function get_cell_dof_basis(f::ExtendedFESpace)
 
   cell_to_val = get_cell_dof_basis(f.space)
+  ref_trait = RefStyle(cell_to_val)
+  cell_to_val = get_array(cell_to_val)
 
   D = num_dims(f.trian)
   T = Float64 # TODO
   void_to_val = Fill(LagrangianDofBasis(T,Point{D,T}[]),length(f.trian.void_to_oldcell))
 
-  ExtendedVector(
-    void_to_val,
-    cell_to_val,
-    f.trian.oldcell_to_cell,
-    f.trian.void_to_oldcell,
-    f.trian.cell_to_oldcell)
+  eb = ExtendedVector(
+         void_to_val,
+         cell_to_val,
+         f.trian.oldcell_to_cell,
+         f.trian.void_to_oldcell,
+         f.trian.cell_to_oldcell)
+
+  cell_dof_basis = GenericCellDofBasis(ref_trait,eb)
 
 end
 
@@ -249,4 +253,3 @@ function TrialFESpace(f::ExtendedFESpace)
   U = TrialFESpace(f.space)
   ExtendedFESpace(U,f.trian)
 end
-
diff --git a/src/FESpaces/FESpaceFactories.jl b/src/FESpaces/FESpaceFactories.jl
index eee302959..1d4a91cc8 100644
--- a/src/FESpaces/FESpaceFactories.jl
+++ b/src/FESpaces/FESpaceFactories.jl
@@ -4,6 +4,7 @@
 function FESpace(;kwargs...)
 
   constraint = _get_kwarg(:constraint,kwargs,nothing)
+  # constraint = nothing
   reffe = _get_kwarg(:reffe,kwargs)
   @notimplementedif !isa(reffe,Symbol) "For the moment, reffe can only be a symbol"
 
@@ -90,6 +91,11 @@ function _setup_hdiv_space(kwargs)
   conformity = _get_kwarg(:conformity,kwargs,true)
   diritags = _get_kwarg(:dirichlet_tags,kwargs,Int[])
   order = _get_kwarg(:order,kwargs,nothing)
+  dofspace = _get_kwarg(:dof_space,kwargs,:reference)
+  ( dofspace == :reference ? true : false )
+
+  is_ref = (dofspace==:reference)
+
   Tf = _get_kwarg(:valuetype,kwargs,VectorValue{1,Float64})
   T = eltype(Tf)
 
@@ -101,7 +107,7 @@ function _setup_hdiv_space(kwargs)
   reffes = [RaviartThomasRefFE(T,p,order) for p in polytopes]
 
   if conformity in [true, :default, :HDiv, :Hdiv]
-      V =  DivConformingFESpace(reffes,model,labels,diritags)
+      V =  DivConformingFESpace(reffes,model,labels,diritags,is_ref)
   else
     s = "Conformity $conformity not implemented for $reffe reference FE on polytopes $(polytopes...)"
     @unreachable s
@@ -118,9 +124,12 @@ function _setup_hcurl_space(kwargs)
   conformity = _get_kwarg(:conformity,kwargs,true)
   diritags = _get_kwarg(:dirichlet_tags,kwargs,Int[])
   order = _get_kwarg(:order,kwargs,nothing)
+  dofspace = _get_kwarg(:dof_space,kwargs,:reference)
   Tf = _get_kwarg(:valuetype,kwargs,VectorValue{1,Float64})
   T = eltype(Tf)
 
+  is_ref = dofspace==:reference
+
   if order == nothing
     @unreachable "order is a mandatory keyword argument in FESpace constructor for Nedelec reference FEs"
   end
@@ -129,7 +138,7 @@ function _setup_hcurl_space(kwargs)
   reffes = [NedelecRefFE(T,p,order) for p in polytopes]
 
   if conformity in [true, :default, :HCurl, :Hcurl]
-      V =  CurlConformingFESpace(reffes,model,labels,diritags)
+      V =  CurlConformingFESpace(reffes,model,labels,diritags,is_ref)
   else
     s = "Conformity $conformity not implemented for $reffe reference FE on polytopes $(polytopes...)"
     @unreachable s
@@ -146,9 +155,14 @@ function _setup_lagrange_spaces(kwargs)
   T = _get_kwarg(:valuetype,kwargs,nothing)
   diritags = _get_kwarg(:dirichlet_tags,kwargs,Int[])
   dirimasks = _get_kwarg(:dirichlet_masks,kwargs,nothing)
+  # dirimasks = nothing
+  dofspace = _get_kwarg(:dof_space,kwargs,:reference)
   labels = _get_kwarg(:labels,kwargs,nothing)
+  # labels = nothing
   model = _get_kwarg(:model,kwargs,nothing)
 
+  is_ref = (dofspace==:reference)
+
   if T == nothing
     @unreachable "valuetype is a mandatory keyword argument in FESpace constructor for Lagrangian reference FEs"
   end
@@ -188,7 +202,7 @@ function _setup_lagrange_spaces(kwargs)
       end
     end
 
-    return  DiscontinuousFESpace(_reffes,trian)
+    return  DiscontinuousFESpace(_reffes,trian,is_ref)
 
   elseif conformity in [true, :default, :H1, :C0]
 
@@ -207,9 +221,9 @@ function _setup_lagrange_spaces(kwargs)
       _reffes = [LagrangianRefFE(T,p,order) for p in polytopes]
     end
     if labels == nothing
-      return GradConformingFESpace(_reffes,model,diritags,dirimasks)
+      return GradConformingFESpace(_reffes,model,diritags,dirimasks,is_ref)
     else
-      return GradConformingFESpace(_reffes,model,labels,diritags,dirimasks)
+      return GradConformingFESpace(_reffes,model,labels,diritags,dirimasks,is_ref)
     end
 
   else
diff --git a/src/FESpaces/FESpaces.jl b/src/FESpaces/FESpaces.jl
index cae949496..1579bc414 100644
--- a/src/FESpaces/FESpaces.jl
+++ b/src/FESpaces/FESpaces.jl
@@ -22,9 +22,12 @@ using Gridap.Algebra
 using Gridap.Polynomials
 using Gridap.TensorValues
 
+using Gridap.ReferenceFEs: evaluate_dof_array
+
 using Gridap.Geometry: CellFieldLike
 using Gridap.Geometry: UnimplementedField
 using Gridap.Geometry: test_cell_field_like
+
 using Gridap.Arrays: _split
 using Gridap.Arrays: Reindexed
 using Gridap.Arrays: IdentityVector
@@ -41,6 +44,8 @@ import Gridap.Geometry: get_cell_map
 import Gridap.Geometry: get_cell_shapefuns
 import Gridap.Geometry: get_reffes
 import Gridap.Geometry: get_cell_type
+import Gridap.Geometry: RefStyle
+import Gridap.Geometry: change_ref_style
 import Gridap.Helpers: operate
 import Gridap.Geometry: similar_object
 import Gridap.Geometry: jump
@@ -131,7 +136,7 @@ export get_cell_dof_basis
 
 export SingleFieldFEFunction
 
-export UnsconstrainedFESpace
+export UnconstrainedFESpace
 export GradConformingFESpace
 export DiscontinuousFESpace
 
@@ -146,6 +151,10 @@ export is_trial
 export is_test
 export attach_dirichlet_bcs
 
+export CellDofBasis
+export GenericCellDofBasis
+export test_cell_dof_basis
+
 export FECellBasisStyle
 export is_a_fe_cell_basis
 
@@ -205,6 +214,8 @@ export update_state_variables!
 
 include("CellBases.jl")
 
+include("CellDofBases.jl")
+
 include("Law.jl")
 
 include("FEFunctions.jl")
diff --git a/src/FESpaces/FESpacesInterfaces.jl b/src/FESpaces/FESpacesInterfaces.jl
index 7f44bc8f6..68fad9a7e 100644
--- a/src/FESpaces/FESpacesInterfaces.jl
+++ b/src/FESpaces/FESpacesInterfaces.jl
@@ -268,4 +268,3 @@ end
   a = apply_kernel!(cmat,k.kmat,mat,cellid)
   (a,vec)
 end
-
diff --git a/src/FESpaces/SingleFieldFEFunctions.jl b/src/FESpaces/SingleFieldFEFunctions.jl
index 51b413a77..669e5ee1f 100644
--- a/src/FESpaces/SingleFieldFEFunctions.jl
+++ b/src/FESpaces/SingleFieldFEFunctions.jl
@@ -1,12 +1,13 @@
 
 """
 """
-struct SingleFieldFEFunction <: CellField
+struct SingleFieldFEFunction{R} <: CellField
   array
   cell_vals
   free_values
   dirichlet_values
   fe_space
+  ref_style::Val{R}
   @doc """
   """
   function SingleFieldFEFunction(
@@ -15,7 +16,10 @@ struct SingleFieldFEFunction <: CellField
     free_values::AbstractVector,
     dirichlet_values::AbstractVector,
     fe_space::SingleFieldFESpace)
-    new(array,cell_vals,free_values,dirichlet_values,fe_space)
+
+    ref_style = RefStyle(get_cell_dof_basis(fe_space))
+    R = get_val_parameter(ref_style)
+    new{R}(array,cell_vals,free_values,dirichlet_values,fe_space,ref_style)
   end
 end
 
@@ -33,3 +37,4 @@ get_cell_map(f::SingleFieldFEFunction) = get_cell_map(f.fe_space)
 
 get_cell_values(f::SingleFieldFEFunction) = f.cell_vals
 
+RefStyle(::Type{SingleFieldFEFunction{R}}) where R = Val{R}()
diff --git a/src/FESpaces/SingleFieldFESpaces.jl b/src/FESpaces/SingleFieldFESpaces.jl
index 935a0cd3b..d9380e8eb 100644
--- a/src/FESpaces/SingleFieldFESpaces.jl
+++ b/src/FESpaces/SingleFieldFESpaces.jl
@@ -135,63 +135,38 @@ function gather_free_values(f::SingleFieldFESpace,cell_vals)
   free_values
 end
 
-"""
-cell_field defined in the reference space with derivatives in the physical one
-"""
-function compute_free_and_dirichlet_values(f::SingleFieldFESpace, cell_field::CellField)
-  cell_vals = _compute_cell_vals(f, cell_field)
-  gather_free_and_dirichlet_values(f,cell_vals)
-end
-
-"""
-"""
-function compute_dirichlet_values(f::SingleFieldFESpace,cell_field::CellField)
-  cell_vals = _compute_cell_vals(f, cell_field)
-  gather_dirichlet_values(f,cell_vals)
-end
-
-"""
-"""
-function compute_free_values(f::SingleFieldFESpace,cell_field::CellField)
-  cell_vals = _compute_cell_vals(f, cell_field)
-  gather_free_values(f,cell_vals)
-end
-
-function _compute_cell_vals(f,cell_field)
-  cell_dof_basis = get_cell_dof_basis(f)
-  cell_dofs = get_cell_dofs(f)
-  cell_vals = evaluate_dof_array(cell_dof_basis,get_array(cell_field))
-  cell_vals
-end
-
 """
 The resulting FE function is in the space (in particular it fulfills Dirichlet BCs
 even in the case that the given cell field does not fulfill them)
 """
 function interpolate(fs::SingleFieldFESpace,object)
-  cell_map = get_cell_map(fs)
-  cell_field = convert_to_cell_field(object,cell_map)
-  free_values = compute_free_values(fs,cell_field)
+  cell_vals = _cell_vals(fs,object)
+  free_values = gather_free_values(fs,cell_vals)
   FEFunction(fs,free_values)
 end
 
+function _cell_vals(fs::SingleFieldFESpace,object)
+  cdb = get_cell_dof_basis(fs)
+  cm = get_cell_map(fs)
+  cf = convert_to_cell_field(object,cm,RefStyle(cdb))
+  cell_vals = evaluate(cdb,cf)
+end
+
 """
 like interpolate, but also compute new degrees of freedom for the dirichlet component.
 The resulting FEFunction does not necessary belongs to the underlying space
 """
 function interpolate_everywhere(fs::SingleFieldFESpace,object)
-  cell_map = get_cell_map(fs)
-  cell_field = convert_to_cell_field(object,cell_map)
-  free_values, dirichlet_values = compute_free_and_dirichlet_values(fs,cell_field)
-  FEFunction(fs,free_values, dirichlet_values)
+  cell_vals = _cell_vals(fs,object)
+  free_values, dirichlet_values = gather_free_and_dirichlet_values(fs,cell_vals)
+  FEFunction(fs,free_values,dirichlet_values)
 end
 
 """
 """
 function interpolate_dirichlet(fs::SingleFieldFESpace,object)
-  cell_map = get_cell_map(fs)
-  cell_field = convert_to_cell_field(object,cell_map)
-  dirichlet_values = compute_dirichlet_values(fs,cell_field)
+  cell_vals = _cell_vals(fs,object)
+  dirichlet_values = gather_dirichlet_values(fs,cell_vals)
   free_values = zero_free_values(fs)
   FEFunction(fs,free_values, dirichlet_values)
 end
@@ -200,12 +175,13 @@ end
 """
 function compute_dirichlet_values_for_tags(f::SingleFieldFESpace,tag_to_object)
   dirichlet_dof_to_tag = get_dirichlet_dof_tag(f)
+  cdb = get_cell_dof_basis(f)
   cell_map = get_cell_map(f)
   dirichlet_values = zero_dirichlet_values(f)
   _tag_to_object = _convert_to_collectable(tag_to_object,num_dirichlet_tags(f))
   for (tag, object) in enumerate(_tag_to_object)
-    cell_field = convert_to_cell_field(object,cell_map)
-    dv = compute_dirichlet_values(f,cell_field)
+    cell_vals = _cell_vals(f,object)
+    dv = gather_dirichlet_values(f,cell_vals)
     _fill_dirichlet_values_for_tag!(dirichlet_values,dv,tag,dirichlet_dof_to_tag)
   end
   dirichlet_values
@@ -228,7 +204,7 @@ function _convert_to_collectable(object::Function,ntags)
   _convert_to_collectable(fill(object,ntags),ntags)
 end
 
+
 function _convert_to_collectable(object::Number,ntags)
   _convert_to_collectable(fill(object,ntags),ntags)
 end
-
diff --git a/src/FESpaces/TrialFESpaces.jl b/src/FESpaces/TrialFESpaces.jl
index 6d93703a9..59bdb1ffa 100644
--- a/src/FESpaces/TrialFESpaces.jl
+++ b/src/FESpaces/TrialFESpaces.jl
@@ -31,7 +31,7 @@ function  _prepare_trial_cell_basis(space)
   a = get_array(cb)
   cm = get_cell_map(cb)
   trial_style = Val{true}()
-  cell_basis = GenericCellBasis(trial_style,a,cm)
+  cell_basis = GenericCellBasis(trial_style,a,cm,RefStyle(cb))
 end
 
 # Genuine functions
diff --git a/src/FESpaces/UnconstrainedFESpaces.jl b/src/FESpaces/UnconstrainedFESpaces.jl
index 5a40e124e..8cdf8b046 100644
--- a/src/FESpaces/UnconstrainedFESpaces.jl
+++ b/src/FESpaces/UnconstrainedFESpaces.jl
@@ -3,7 +3,7 @@
 Generic implementation of an unconstrained single-field FE space
 Private fields and type parameters
 """
-struct UnsconstrainedFESpace{A,B,C} <: SingleFieldFESpace
+struct UnconstrainedFESpace{A,B,C} <: SingleFieldFESpace
   nfree::Int
   ndirichlet::Int
   cell_dofs::A
@@ -15,19 +15,17 @@ struct UnsconstrainedFESpace{A,B,C} <: SingleFieldFESpace
 
   @doc """
   """
-  function UnsconstrainedFESpace(
+  function UnconstrainedFESpace(
     nfree::Int,
     ndirichlet::Int,
     cell_dofs::AbstractArray,
-    cell_shapefuns::AbstractArray,
-    cell_dof_basis::AbstractArray,
+    cell_basis::CellBasis,
+    cell_dof_basis::CellDofBasis,
     cell_map::AbstractArray,
     dirichlet_dof_tag::Vector{Int8},
     dirichlet_cells::Vector{Int},
     ntags) where T
 
-    cell_basis = GenericCellBasis(cell_shapefuns,cell_map)
-
     A = typeof(cell_dofs)
     B = typeof(cell_basis)
     C = typeof(cell_dof_basis)
@@ -46,53 +44,53 @@ end
 
 # FESpace interface
 
-constraint_style(::Type{<:UnsconstrainedFESpace}) = Val{false}()
+constraint_style(::Type{<:UnconstrainedFESpace}) = Val{false}()
 
-function num_free_dofs(f::UnsconstrainedFESpace)
+function num_free_dofs(f::UnconstrainedFESpace)
   f.nfree
 end
 
-function get_cell_basis(f::UnsconstrainedFESpace)
+function get_cell_basis(f::UnconstrainedFESpace)
   f.cell_basis
 end
 
-function zero_free_values(::Type{T},f::UnsconstrainedFESpace) where T
+function zero_free_values(::Type{T},f::UnconstrainedFESpace) where T
   zeros(T,num_free_dofs(f))
 end
 
 # SingleFieldFESpace interface
 
-function get_cell_dofs(f::UnsconstrainedFESpace)
+function get_cell_dofs(f::UnconstrainedFESpace)
   f.cell_dofs
 end
 
-function get_cell_dof_basis(f::UnsconstrainedFESpace)
+function get_cell_dof_basis(f::UnconstrainedFESpace)
   f.cell_dof_basis
 end
 
-function num_dirichlet_dofs(f::UnsconstrainedFESpace)
+function num_dirichlet_dofs(f::UnconstrainedFESpace)
   f.ndirichlet
 end
 
-function num_dirichlet_tags(f::UnsconstrainedFESpace)
+function num_dirichlet_tags(f::UnconstrainedFESpace)
   f.ntags
 end
 
-function zero_dirichlet_values(f::UnsconstrainedFESpace)
+function zero_dirichlet_values(f::UnconstrainedFESpace)
   T = Float64 # TODO
   zeros(T,num_dirichlet_dofs(f))
 end
 
-function get_dirichlet_dof_tag(f::UnsconstrainedFESpace)
+function get_dirichlet_dof_tag(f::UnconstrainedFESpace)
   f.dirichlet_dof_tag
 end
 
-function scatter_free_and_dirichlet_values(f::UnsconstrainedFESpace,free_values,dirichlet_values)
+function scatter_free_and_dirichlet_values(f::UnconstrainedFESpace,free_values,dirichlet_values)
   cell_dofs = get_cell_dofs(f)
   LocalToGlobalPosNegArray(cell_dofs,free_values,dirichlet_values)
 end
 
-function gather_free_and_dirichlet_values(f::UnsconstrainedFESpace,cell_vals)
+function gather_free_and_dirichlet_values(f::UnconstrainedFESpace,cell_vals)
 
   cell_dofs = get_cell_dofs(f)
   cache_vals = array_cache(cell_vals)
@@ -113,7 +111,7 @@ function gather_free_and_dirichlet_values(f::UnsconstrainedFESpace,cell_vals)
   (free_vals, dirichlet_vals)
 end
 
-function gather_dirichlet_values(f::UnsconstrainedFESpace,cell_vals)
+function gather_dirichlet_values(f::UnconstrainedFESpace,cell_vals)
 
   cell_dofs = get_cell_dofs(f)
   cache_vals = array_cache(cell_vals)
diff --git a/src/Fields/DiffOperators.jl b/src/Fields/DiffOperators.jl
index 8f499d945..d39582f8f 100644
--- a/src/Fields/DiffOperators.jl
+++ b/src/Fields/DiffOperators.jl
@@ -9,7 +9,7 @@ function symmetric_gradient end
 """
     symmetric_gradient(f)
 """
-symmetric_gradient(f) = symmetic_part(gradient(f))
+symmetric_gradient(f) = symmetric_part(gradient(f))
 
 """
     const ε = symmetric_gradient
diff --git a/src/Fields/FieldOperations.jl b/src/Fields/FieldOperations.jl
index 67420ec92..5e35ce939 100644
--- a/src/Fields/FieldOperations.jl
+++ b/src/Fields/FieldOperations.jl
@@ -40,7 +40,7 @@ end
 
 # Unary operations on fields and arrays of fields
 
-for op in (:+,:-,:tr, :transpose, :adjoint, :symmetic_part)
+for op in (:+,:-,:tr, :transpose, :adjoint, :symmetric_part)
   @eval begin
 
     function ($op)(f::Field)
diff --git a/src/Fields/Fields.jl b/src/Fields/Fields.jl
index f65e4ab19..ce2116f95 100644
--- a/src/Fields/Fields.jl
+++ b/src/Fields/Fields.jl
@@ -69,12 +69,14 @@ export AffineMap
 export field_operation
 export field_array_operation
 
+export function_field
+
 import Gridap.Arrays: kernel_cache
 import Gridap.Arrays: apply_kernel!
 import Gridap.Arrays: kernel_return_type
 import Gridap.TensorValues: outer
 import Gridap.TensorValues: inner
-import Gridap.TensorValues: symmetic_part
+import Gridap.TensorValues: symmetric_part
 import Base: +, - , *
 import LinearAlgebra: cross
 import LinearAlgebra: tr
@@ -85,6 +87,8 @@ include("FieldInterface.jl")
 
 include("MockFields.jl")
 
+include("FunctionFields.jl")
+
 include("ConstantFields.jl")
 
 include("Homothecies.jl")
diff --git a/src/Fields/FunctionFields.jl b/src/Fields/FunctionFields.jl
new file mode 100644
index 000000000..bb324c94e
--- /dev/null
+++ b/src/Fields/FunctionFields.jl
@@ -0,0 +1,25 @@
+struct FunctionField{F} <: Field
+  f::F
+end
+
+function_field(f::Function) =  FunctionField(f)
+
+function field_cache(f::FunctionField,x)
+  nx = length(x)
+  Te = eltype(x)
+  c = zeros(return_type(f.f,Te),nx)
+  CachedArray(c)
+end
+
+function evaluate_field!(c,f::FunctionField,x)
+  nx = length(x)
+  setsize!(c,(nx,))
+  for i in eachindex(x)
+    c[i] = f.f(x[i])
+  end
+  c
+end
+
+function field_gradient(f::FunctionField)
+  FunctionField(gradient(f.f))
+end
diff --git a/src/Geometry/CellFields.jl b/src/Geometry/CellFields.jl
index 687153a55..dee177ac3 100644
--- a/src/Geometry/CellFields.jl
+++ b/src/Geometry/CellFields.jl
@@ -18,6 +18,53 @@ function get_cell_map(cf::CellFieldLike)
   @abstractmethod
 end
 
+"""
+This trait returns `Val{true}()` when the `CellFieldLike` is defined in a
+reference finite element space, and `Val{false}()` when it is defined in the
+physical space
+"""
+RefStyle(::Type{<:CellFieldLike}) = @notimplemented
+
+# We use duck typing here for all types marked with the RefStyle
+# @santiagobadia : That is dangerous, it creates overflows if called with
+RefStyle(::T) where T = RefStyle(T)
+is_in_ref_space(::Type{T}) where T = get_val_parameter(RefStyle(T))
+is_in_ref_space(::T) where T = is_in_ref_space(T)
+is_in_physical_space(::Type{T}) where T = !is_in_ref_space(T)
+is_in_physical_space(a::T) where T = !is_in_ref_space(T)
+
+to_ref_space(a::CellFieldLike) = _to_ref_space(a,RefStyle(a))
+_to_ref_space(a,::Val{true}) = a
+function _to_ref_space(a,::Val{false})
+  cell_map = get_cell_map(a)
+  array = compose(  get_array(a), cell_map  )
+  no = similar_object(a,array)
+  change_ref_style(no)
+end
+
+to_physical_space(a::CellFieldLike) = _to_physical_space(a,RefStyle(a))
+_to_physical_space(a,::Val{true}) = @notimplemented # and probably not doable in some cases
+_to_physical_space(a,::Val{false}) = a
+
+# Assumption : x ALWAIS defined in the reference space
+# In the future we can also add the RefStyle to x
+"""
+    evaluate(cf::CellFieldLike,x)
+"""
+function evaluate(cf::CellFieldLike,x::AbstractArray)
+  _evaluate(cf,x,RefStyle(cf))
+end
+
+function _evaluate(cf::CellFieldLike,x::AbstractArray,::Val{true})
+  evaluate_field_array(get_array(cf),x)
+end
+
+function _evaluate(cf::CellFieldLike,x::AbstractArray,::Val{false})
+  cm = get_cell_map(cf)
+  _x = evaluate(cm,x)
+  evaluate_field_array(get_array(cf),_x)
+end
+
 """
     similar_object(cf::CellFieldLike,array::AbstractArray)
 """
@@ -32,6 +79,15 @@ function similar_object(cf1::CellFieldLike,cf2::CellFieldLike,array::AbstractArr
   @abstractmethod
 end
 
+"""
+   change_ref_style(cf::CellFieldLike)
+
+Return an object with the same array and metadata as in `cf`, except for `RefStyle` which is changed.
+"""
+function change_ref_style(cf::CellFieldLike)
+  @abstractmethod
+end
+
 """
     gradient(cf::CellFieldLike)
 """
@@ -63,14 +119,12 @@ function test_cell_field_like(cf::CellFieldLike,x::AbstractArray,b::AbstractArra
     g = evaluate(gradient(cf),x)
     test_array(g,grad,pred)
   end
-end
-
-"""
-    evaluate(cf::CellFieldLike,x)
-"""
-function evaluate(cf::CellFieldLike,x)
-  a = get_array(cf)
-  evaluate_field_array(a,x)
+  rs = RefStyle(cf)
+  @test isa(get_val_parameter(rs),Bool)
+  _cf = change_ref_style(cf)
+  @test get_array(_cf) === get_array(cf)
+  @test is_in_ref_space(cf) == !is_in_ref_space(_cf)
+  @test is_in_physical_space(cf) == !is_in_physical_space(_cf)
 end
 
 """
@@ -97,12 +151,22 @@ end
 
 function similar_object(cf::CellField,array::AbstractArray)
   cm = get_cell_map(cf)
-  GenericCellField(array,cm)
+  GenericCellField(array,cm,RefStyle(cf))
 end
 
 function similar_object(cf1::CellField,cf2::CellField,array::AbstractArray)
   cm = get_cell_map(cf1)
-  GenericCellField(array,cm)
+  @assert is_in_ref_space(cf1) == is_in_ref_space(cf2)
+  GenericCellField(array,cm,RefStyle(cf1))
+end
+
+function change_ref_style(cf::CellField)
+  ref_sty = RefStyle(cf)
+  bool = !get_val_parameter(ref_sty)
+  new_sty = Val{bool}()
+  ar = get_array(cf)
+  cm = get_cell_map(cf)
+  GenericCellField(ar,cm,new_sty)
 end
 
 # Diff operations
@@ -139,51 +203,74 @@ end
 
 function operate(op,cf1::CellField,object)
   cm = get_cell_map(cf1)
-  cf2 = convert_to_cell_field(object,cm)
+  cf2 = convert_to_cell_field(object,cm,RefStyle(cf1))
   operate(op,cf1,cf2)
 end
 
 function operate(op,object,cf2::CellField)
   cm = get_cell_map(cf2)
-  cf1 = convert_to_cell_field(object,cm)
+  cf1 = convert_to_cell_field(object,cm,RefStyle(cf2))
   operate(op,cf1,cf2)
 end
 
 # Conversions
 
+function convert_to_cell_field(object,cell_map)
+  ref_style = Val{true}()
+  convert_to_cell_field(object,cell_map,ref_style)
+end
+
 """
-    convert_to_cell_field(object::CellField,cell_map)
+    convert_to_cell_field(object,cell_map,ref_style)
 """
-function convert_to_cell_field(object::CellField,cell_map)
+function convert_to_cell_field(object::CellField,cell_map,ref_style::Val)
+  @assert RefStyle(object) == ref_style
   object
 end
 
-function convert_to_cell_field(object::AbstractArray,cell_map)
+function convert_to_cell_field(object::CellField,cell_map)
+  ref_style = RefStyle(object)
+  convert_to_cell_field(object,cell_map,ref_style)
+end
+
+function convert_to_cell_field(object::AbstractArray,cell_map,ref_style::Val)
   @assert length(object) == length(cell_map)
-  GenericCellField(object,cell_map)
+  GenericCellField(object,cell_map,ref_style)
 end
 
-function convert_to_cell_field(object::Function,cell_map)
+function convert_to_cell_field(object::Function,cell_map,ref_style::Val{true})
   b = compose(object,cell_map)
-  GenericCellField(b,cell_map)
+  GenericCellField(b,cell_map,Val{true}())
 end
 
-function convert_to_cell_field(object::Number,cell_map)
+function convert_to_cell_field(fun::Function,cell_map,ref_style::Val{false})
+  field = function_field(fun)
+  cell_field = Fill(field,length(cell_map))
+  GenericCellField(cell_field,cell_map,Val{false}())
+end
+
+function convert_to_cell_field(object::Number,cell_map,ref_style::Val)
   a = Fill(object,length(cell_map))
-  GenericCellField(a,cell_map)
+  GenericCellField(a,cell_map,ref_style)
 end
 
 # Concrete implementation
 
 """
-    struct GenericCellField <: CellField
+struct GenericCellField{R} <: CellField
       array::AbstractArray
       cell_map::AbstractArray
+      ref_trait::Val{R}
     end
 """
-struct GenericCellField <: CellField
+struct GenericCellField{R} <: CellField
   array::AbstractArray
   cell_map::AbstractArray
+  ref_trait::Val{R}
+end
+
+function GenericCellField(array::AbstractArray,cell_map::AbstractArray)
+  GenericCellField(array,cell_map,Val{true}())
 end
 
 function get_array(cf::GenericCellField)
@@ -194,6 +281,9 @@ function get_cell_map(cf::GenericCellField)
   cf.cell_map
 end
 
+function RefStyle(::Type{<:GenericCellField{R}}) where {R}
+  Val{R}()
+end
 
 # Skeleton related
 
@@ -281,13 +371,12 @@ end
 
 function operate(op,cf1::SkeletonCellField,object)
   cm = get_cell_map(cf1)
-  cf2 = convert_to_cell_field(object,cm)
+  cf2 = convert_to_cell_field(object,cm,RefStyle(cf1.left))
   operate(op,cf1,cf2)
 end
 
 function operate(op,object,cf2::SkeletonCellField)
   cm = get_cell_map(cf2)
-  cf1 = convert_to_cell_field(object,cm)
+  cf1 = convert_to_cell_field(object,cm,RefStyle(cf2.left))
   operate(op,cf1,cf2)
 end
-
diff --git a/src/Geometry/CellQuadratures.jl b/src/Geometry/CellQuadratures.jl
index 4435ec260..907d6be20 100644
--- a/src/Geometry/CellQuadratures.jl
+++ b/src/Geometry/CellQuadratures.jl
@@ -101,7 +101,8 @@ function integrate(cell_field,trian::Triangulation,quad::CellQuadrature)
   q = get_coordinates(quad)
   w = get_weights(quad)
   j = gradient(cell_map)
-  f = convert_to_cell_field(cell_field,cell_map)
+  _f = CellField(cell_field,trian)
+  f = to_ref_space(_f)
   @assert length(f) == length(cell_map) "Are you using the right triangulation to integrate?"
   @assert length(f) == length(w) "Are you using the right quadrature to integrate?"
   integrate(get_array(f),q,w,j)
diff --git a/src/Geometry/Geometry.jl b/src/Geometry/Geometry.jl
index 4c1ea2355..859d152ce 100644
--- a/src/Geometry/Geometry.jl
+++ b/src/Geometry/Geometry.jl
@@ -88,9 +88,16 @@ import Gridap.Arrays: apply_kernel!
 export CellField
 export GenericCellField
 export SkeletonCellField
+export RefStyle
 export similar_object
+export change_ref_style
 export test_cell_field
 export convert_to_cell_field
+export to_ref_space
+export to_physical_space
+export is_in_physical_space
+export is_in_ref_space
+export cell_field_from_function
 
 export GridTopology
 export num_cells
diff --git a/src/Geometry/Triangulations.jl b/src/Geometry/Triangulations.jl
index 78f63f354..fe5ff704f 100644
--- a/src/Geometry/Triangulations.jl
+++ b/src/Geometry/Triangulations.jl
@@ -233,7 +233,8 @@ end
     restrict(cf::CellField,trian::Triangulation)
 """
 function restrict(cf::CellField,trian::Triangulation)
-  a = get_array(cf)
+  _cf = to_ref_space(cf)
+  a = get_array(_cf)
   r = restrict(a,trian)
   _restrict_cell_field(r,trian)
 end
diff --git a/src/Inference/Inference.jl b/src/Inference/Inference.jl
index 53fb44196..c5286e58f 100644
--- a/src/Inference/Inference.jl
+++ b/src/Inference/Inference.jl
@@ -32,8 +32,8 @@ The underlying implementation uses the function [`testargs`](@ref) to generate
 some test values in order to call the function and determine the returned type.
 This mechanism does not use `Base._return_type`. One of the advantages is
 that the given function `f` is called, and thus, meaningful error messages
-will be displayed if there is any error in `f`. 
-    
+will be displayed if there is any error in `f`.
+
 """
 function return_type(f::Function,Ts...)
   args = testargs(f,Ts...)
@@ -115,11 +115,12 @@ Returns an arbitrary instance of type `T`. It defaults to `zero(T)` for
 non-array types and to an empty array for array types.
 This function is used to compute the default test arguments in
 [`testargs`](@ref).
-It can be overloaded for new types `T` if `zero(T)` does not makes sense. 
+It can be overloaded for new types `T` if `zero(T)` does not makes sense.
 """
 function testvalue end
 
 testvalue(::Type{T}) where T = zero(T)
+testvalue(v) = testvalue(typeof(v))
 
 function testvalue(::Type{T}) where T<:AbstractArray{E,N} where {E,N}
    similar(T,fill(0,N)...)
diff --git a/src/MultiField/MultiField.jl b/src/MultiField/MultiField.jl
index 8eada98ef..7d9d844bb 100644
--- a/src/MultiField/MultiField.jl
+++ b/src/MultiField/MultiField.jl
@@ -41,6 +41,7 @@ import Gridap.Geometry: restrict
 import Gridap.Fields: integrate
 import Gridap.Fields: evaluate
 import Gridap.FESpaces: TrialStyle
+import Gridap.FESpaces: RefStyle
 import Gridap.FESpaces: FECellBasisStyle
 import Gridap.FESpaces: FEFunctionStyle
 import Gridap.FESpaces: num_free_dofs
diff --git a/src/MultiField/MultiFieldCellBases.jl b/src/MultiField/MultiFieldCellBases.jl
index e8dcac096..d1b9ca326 100644
--- a/src/MultiField/MultiFieldCellBases.jl
+++ b/src/MultiField/MultiFieldCellBases.jl
@@ -1,12 +1,15 @@
 
-struct CellBasisWithFieldID{S} <: CellBasis
+struct CellBasisWithFieldID{S,R} <: CellBasis
   trial_style::Val{S}
   cell_basis::CellBasis
   field_id::Int
+  ref_style::Val{R}
   function CellBasisWithFieldID(cell_basis::CellBasis,field_id::Integer)
     trial_style = TrialStyle(cell_basis)
     S = get_val_parameter(trial_style)
-    new{S}(trial_style,cell_basis,field_id)
+    ref_style = RefStyle(cell_basis)
+    R = get_val_parameter(ref_style)
+    new{S,R}(trial_style,cell_basis,field_id,ref_style)
   end
 end
 
@@ -14,7 +17,9 @@ get_array(cb::CellBasisWithFieldID) = get_array(cb.cell_basis)
 
 get_cell_map(cb::CellBasisWithFieldID) = get_cell_map(cb.cell_basis)
 
-TrialStyle(::Type{CellBasisWithFieldID{S}}) where S = Val{S}()
+TrialStyle(::Type{CellBasisWithFieldID{S,R}}) where {S,R} = Val{S}()
+
+RefStyle(::Type{CellBasisWithFieldID{S,R}}) where {S,R} = Val{R}()
 
 function similar_object(cb::CellBasisWithFieldID,array::AbstractArray)
   cell_basis = similar_object(cb.cell_basis,array)
@@ -45,7 +50,8 @@ function  _similar_cell_basis_with_field_id_test_trial(v,a,b)
   field_id_rows = a.field_id
   field_id_cols = b.field_id
   cell_matrix_field = similar_object(a.cell_basis,b.cell_basis,v)
-  CellMatrixFieldWithFieldIds(cell_matrix_field,field_id_rows,field_id_cols)
+  @assert is_in_ref_space(a) == is_in_ref_space(b)
+  CellMatrixFieldWithFieldIds(cell_matrix_field,field_id_rows,field_id_cols,RefStyle(a))
 end
 
 function operate(op,a::CellBasisWithFieldID,b::CellBasisWithFieldID)
@@ -76,12 +82,15 @@ function _operate_cell_basis_with_field_id(op,a,b,atrial,btrial)
   _operate_cell_basis(op,a,b,TrialStyle(a),TrialStyle(b))
 end
 
-struct CellMatrixFieldWithFieldIds <: CellMatrixField
+struct CellMatrixFieldWithFieldIds{R} <: CellMatrixField
   cell_matrix_field::CellMatrixField
   field_id_rows::Int
   field_id_cols::Int
+  ref_style::Val{R}
 end
 
+RefStyle(::Type{CellMatrixFieldWithFieldIds{R}}) where R = Val{R}()
+
 get_array(a::CellMatrixFieldWithFieldIds) = get_array(a.cell_matrix_field)
 
 get_cell_map(a::CellMatrixFieldWithFieldIds) = get_cell_map(a.cell_matrix_field)
@@ -94,7 +103,7 @@ function similar_object(cf::CellMatrixFieldWithFieldIds,a::AbstractArray)
   cell_matrix_field = similar_object(cf.cell_matrix_field,a)
   field_id_rows = cf.field_id_rows
   field_id_cols = cf.field_id_cols
-  CellMatrixFieldWithFieldIds(cell_matrix_field,field_id_rows,field_id_cols)
+  CellMatrixFieldWithFieldIds(cell_matrix_field,field_id_rows,field_id_cols,RefStyle(cf))
 end
 
 function similar_object(a::CellMatrixFieldWithFieldIds,b::CellMatrixFieldWithFieldIds,v::AbstractArray)
@@ -171,24 +180,28 @@ function integrate(cb::BlockTracker,trian::Triangulation,quad::CellQuadrature)
   MultiFieldCellArray(blocks,block_ids)
 end
 
-struct MultiCellBasis{S} <: GridapType
+struct MultiCellBasis{S,R} <: GridapType
   blocks::Vector{CellBasisWithFieldID{S}}
   function MultiCellBasis(blocks::Vector{<:CellBasis})
     cb = first(blocks)
     S = is_trial(cb)
+    R = is_in_ref_space(cb)
     new_blocks = CellBasisWithFieldID{S}[]
     for (field_id, cell_basis) in enumerate(blocks)
-      @assert is_trial(cell_basis) == S "All the provided bases need to be either test or trial"
+      @assert is_trial(cell_basis) == S "All the provided bases must be of the same type (trial or test)"
+      @assert is_in_ref_space(cell_basis) == R "All the provided bases must be defined in the same space (reference or physical)"
       block = CellBasisWithFieldID(cell_basis,field_id)
       push!(new_blocks,block)
     end
-    new{S}(new_blocks)
+    new{S,R}(new_blocks)
   end
 end
 
 FECellBasisStyle(::Type{<:MultiCellBasis}) = Val{true}()
 
-TrialStyle(::Type{MultiCellBasis{S}}) where S = Val{S}()
+TrialStyle(::Type{MultiCellBasis{S,R}}) where {S,R} = Val{S}()
+
+RefStyle(::Type{MultiCellBasis{S,R}}) where {S,R} = Val{R}()
 
 num_fields(m::MultiCellBasis) = length(m.blocks)
 
@@ -224,6 +237,3 @@ end
   end
   vec
 end
-
-
-
diff --git a/src/MultiField/MultiFieldFESpaces.jl b/src/MultiField/MultiFieldFESpaces.jl
index daa4879c7..2f8dbb639 100644
--- a/src/MultiField/MultiFieldFESpaces.jl
+++ b/src/MultiField/MultiFieldFESpaces.jl
@@ -1,7 +1,7 @@
 
 abstract type MultiFieldStyle end
 
-struct ConsequtiveMultiFieldStyle <: MultiFieldStyle end
+struct ConsecutiveMultiFieldStyle <: MultiFieldStyle end
 
 struct StridedMultiFieldStyle <: MultiFieldStyle end
 
@@ -30,7 +30,7 @@ end
     MultiFieldFESpace(spaces::Vector{<:SingleFieldFESpace})
 """
 function MultiFieldFESpace(spaces::Vector{<:SingleFieldFESpace})
-  MultiFieldFESpace(spaces,ConsequtiveMultiFieldStyle())
+  MultiFieldFESpace(spaces,ConsecutiveMultiFieldStyle())
 end
 
 MultiFieldStyle(::Type{MultiFieldFESpace{S,B}}) where {S,B} = S()
@@ -129,7 +129,7 @@ function  _restrict_to_field(f,::MultiFieldStyle,free_values,field)
   @notimplemented
 end
 
-function  _restrict_to_field(f,::ConsequtiveMultiFieldStyle,free_values,field)
+function  _restrict_to_field(f,::ConsecutiveMultiFieldStyle,free_values,field)
   offsets = compute_field_offsets(f)
   U = f.spaces
   pini = offsets[field] + 1
@@ -146,7 +146,7 @@ function _get_cell_dofs(f,::MultiFieldStyle)
 end
 
 
-function _get_cell_dofs(f,::ConsequtiveMultiFieldStyle)
+function _get_cell_dofs(f,::ConsecutiveMultiFieldStyle)
   offsets = compute_field_offsets(f)
   spaces = f.spaces
   function fun(i,space)
@@ -171,13 +171,13 @@ function _sum_if_first_positive(a,b)
   end
 end
 
-# API for the ConsequtiveMultiFieldStyle
+# API for the ConsecutiveMultiFieldStyle
 
 """
     compute_field_offsets(f::MultiFieldFESpace)
 """
 function compute_field_offsets(f::MultiFieldFESpace)
-  @assert MultiFieldStyle(f) == ConsequtiveMultiFieldStyle()
+  @assert MultiFieldStyle(f) == ConsecutiveMultiFieldStyle()
   U = f.spaces
   n = length(U)
   offsets = zeros(Int,n)
diff --git a/src/Polynomials/ChangeBasis.jl b/src/Polynomials/ChangeBasis.jl
index ab87e9f93..1177ffa0b 100644
--- a/src/Polynomials/ChangeBasis.jl
+++ b/src/Polynomials/ChangeBasis.jl
@@ -30,6 +30,10 @@ function change_basis(basis,changeofbasis::AbstractMatrix)
   BasisFromChangeOfBasis(basis,changeofbasis)
 end
 
+function kernel_return_type(::typeof(change_basis),prebasis,matrix_inv)
+  typeof(change_basis(prebasis,matrix_inv))
+end
+
 struct BasisFromChangeOfBasis{B,M} <: Field
   basis::B
   change::M
diff --git a/src/Polynomials/Polynomials.jl b/src/Polynomials/Polynomials.jl
index 8ec871c0b..3fed6e78c 100644
--- a/src/Polynomials/Polynomials.jl
+++ b/src/Polynomials/Polynomials.jl
@@ -21,6 +21,7 @@ import Gridap.Fields: evaluate_gradient!
 import Gridap.Fields: gradient_cache
 import Gridap.Fields: evaluate_hessian!
 import Gridap.Fields: hessian_cache
+import Gridap.Arrays: kernel_return_type
 
 export MonomialBasis
 export QGradMonomialBasis
diff --git a/src/ReferenceFEs/Dofs.jl b/src/ReferenceFEs/Dofs.jl
index 41e27f5f3..92001e1c8 100644
--- a/src/ReferenceFEs/Dofs.jl
+++ b/src/ReferenceFEs/Dofs.jl
@@ -160,4 +160,3 @@ end
 function kernel_return_type(k::DofEval,dof,field)
   dof_return_type(dof,field)
 end
-
diff --git a/src/ReferenceFEs/LagrangianDofBases.jl b/src/ReferenceFEs/LagrangianDofBases.jl
index ba9598acc..50e254942 100644
--- a/src/ReferenceFEs/LagrangianDofBases.jl
+++ b/src/ReferenceFEs/LagrangianDofBases.jl
@@ -36,6 +36,8 @@ function LagrangianDofBasis(::Type{T},nodes::Vector{<:Point}) where T
   LagrangianDofBasis(nodes,r...)
 end
 
+get_nodes(b::LagrangianDofBasis) = b.nodes
+
 function _generate_dof_layout_node_major(::Type{<:Real},nnodes::Integer)
   ndofs = nnodes
   dof_to_comp = ones(Int,ndofs)
@@ -127,5 +129,3 @@ function _evaluate_lagr_dof!(c::AbstractMatrix,node_pdof_comp_to_val,node_and_co
   end
   r
 end
-
-
diff --git a/src/ReferenceFEs/NodalReferenceFEs.jl b/src/ReferenceFEs/NodalReferenceFEs.jl
index 9357a1b5d..9109b8b7f 100644
--- a/src/ReferenceFEs/NodalReferenceFEs.jl
+++ b/src/ReferenceFEs/NodalReferenceFEs.jl
@@ -76,7 +76,7 @@ function test_nodal_reference_fe(reffe::NodalReferenceFE)
   @test isa(get_face_nodes(reffe),Vector{Vector{Int}})
 end
 
-# Dafault API
+# Default API
 
 """
     num_nodes(reffe::NodalReferenceFE)
@@ -147,9 +147,8 @@ function get_face_nodes(reffe::NodalReferenceFE,d::Integer)
 end
 
 # Generic implementation
-
 """
-  struct GenericNodalRefFE{D,T,V} <: NodalReferenceFE{D}
+  struct GenericNodalCartesianRefFE{D,T,V} <: NodalReferenceFE{D}
     reffe::GenericRefFE{D}
     node_coordinates::Vector{Point{D,T}}
     node_and_comp_to_dof::Vector{V}
@@ -158,7 +157,7 @@ end
     face_nodes::Vector{Vector{Int}}
   end
 """
-struct GenericNodalRefFE{D,T,V} <: NodalReferenceFE{D}
+struct GenericNodalCartesianRefFE{D,T,V} <: NodalReferenceFE{D}
   reffe::GenericRefFE{D}
   node_coordinates::Vector{Point{D,T}}
   node_and_comp_to_dof::Vector{V}
@@ -169,31 +168,30 @@ end
 
 # NodalReffe
 
-get_node_coordinates(reffe::GenericNodalRefFE) = reffe.node_coordinates
-
-get_node_and_comp_to_dof(reffe::GenericNodalRefFE) = reffe.node_and_comp_to_dof
+get_node_coordinates(reffe::GenericNodalCartesianRefFE) = reffe.node_coordinates
 
-get_face_own_nodes(reffe::GenericNodalRefFE) = reffe.face_own_nodes
+get_node_and_comp_to_dof(reffe::GenericNodalCartesianRefFE) = reffe.node_and_comp_to_dof
 
-get_face_own_nodes_permutations(reffe::GenericNodalRefFE) = reffe.face_own_nodes_permutations
+get_face_own_nodes(reffe::GenericNodalCartesianRefFE) = reffe.face_own_nodes
 
-get_face_nodes(reffe::GenericNodalRefFE) = reffe.face_nodes
+get_face_own_nodes_permutations(reffe::GenericNodalCartesianRefFE) = reffe.face_own_nodes_permutations
 
-# Reffe 
+get_face_nodes(reffe::GenericNodalCartesianRefFE) = reffe.face_nodes
 
-num_dofs(reffe::GenericNodalRefFE) = reffe.reffe.ndofs
+# Reffe
 
-get_polytope(reffe::GenericNodalRefFE) = reffe.reffe.polytope
+num_dofs(reffe::GenericNodalCartesianRefFE) = reffe.reffe.ndofs
 
-get_prebasis(reffe::GenericNodalRefFE) = reffe.reffe.prebasis
+get_polytope(reffe::GenericNodalCartesianRefFE) = reffe.reffe.polytope
 
-get_dof_basis(reffe::GenericNodalRefFE) = reffe.reffe.dofs
+get_prebasis(reffe::GenericNodalCartesianRefFE) = reffe.reffe.prebasis
 
-get_face_own_dofs(reffe::GenericNodalRefFE) = reffe.reffe.face_own_dofs
+get_dof_basis(reffe::GenericNodalCartesianRefFE) = reffe.reffe.dofs
 
-get_face_own_dofs_permutations(reffe::GenericNodalRefFE) = reffe.reffe.face_own_dofs_permutations
+get_face_own_dofs(reffe::GenericNodalCartesianRefFE) = reffe.reffe.face_own_dofs
 
-get_face_dofs(reffe::GenericNodalRefFE) = reffe.reffe.face_own_dofs
+get_face_own_dofs_permutations(reffe::GenericNodalCartesianRefFE) = reffe.reffe.face_own_dofs_permutations
 
-get_shapefuns(reffe::GenericNodalRefFE) = reffe.reffe.shapefuns
+get_face_dofs(reffe::GenericNodalCartesianRefFE) = reffe.reffe.face_own_dofs
 
+get_shapefuns(reffe::GenericNodalCartesianRefFE) = reffe.reffe.shapefuns
diff --git a/src/ReferenceFEs/RaviartThomasRefFEs.jl b/src/ReferenceFEs/RaviartThomasRefFEs.jl
index 87fb7f8c8..33a9464ca 100644
--- a/src/ReferenceFEs/RaviartThomasRefFEs.jl
+++ b/src/ReferenceFEs/RaviartThomasRefFEs.jl
@@ -179,6 +179,12 @@ struct MomentBasedDofBasis{P,V} <: Dof
   face_moments::Vector{Array{V}}
   face_nodes::Vector{UnitRange{Int}}
 
+  function MomentBasedDofBasis(nodes,f_moments,f_nodes)
+    P = eltype(nodes)
+    V = eltype(eltype(f_moments))
+    new{P,V}(nodes,f_moments,f_nodes)
+  end
+
   function MomentBasedDofBasis(f_nodes,f_moments)
     P = eltype(eltype(f_nodes))
     V = eltype(eltype(f_moments))
@@ -200,6 +206,10 @@ struct MomentBasedDofBasis{P,V} <: Dof
   end
 end
 
+get_nodes(b::MomentBasedDofBasis) = b.nodes
+get_face_moments(b::MomentBasedDofBasis) = b.face_moments
+get_face_nodes_dofs(b::MomentBasedDofBasis) = b.face_nodes
+
 function num_dofs(b::MomentBasedDofBasis)
   n = 0
   for m in b.face_moments
diff --git a/src/ReferenceFEs/ReferenceFEs.jl b/src/ReferenceFEs/ReferenceFEs.jl
index 07d46ae68..d2f22b93d 100644
--- a/src/ReferenceFEs/ReferenceFEs.jl
+++ b/src/ReferenceFEs/ReferenceFEs.jl
@@ -80,12 +80,16 @@ export TET_AXIS
 export INVALID_PERM
 
 export Dof
+export get_nodes
+export get_face_moments
+export get_face_nodes_dofs
+export get_nodes
 export evaluate_dof!
 export evaluate_dof
 export dof_cache
 export dof_return_type
 export test_dof
-export evaluate_dof_array
+# export evaluate_dof_array
 
 export ReferenceFE
 export GenericRefFE
@@ -102,9 +106,7 @@ export test_reference_fe
 export num_dofs
 
 export NodalReferenceFE
-export GenericNodalRefFE
-export get_face_own_nodes
-export get_face_nodes
+export GenericNodalCartesianRefFE
 export get_face_own_nodes_permutations
 export get_own_nodes_permutations
 export get_node_coordinates
@@ -129,6 +131,10 @@ export is_Q
 export is_P
 export is_S
 
+export MomentBasedDofBasis
+export get_face_own_nodes
+export get_face_nodes
+
 export VERTEX1
 export SEG2
 export TRI3
diff --git a/src/TensorValues/Operations.jl b/src/TensorValues/Operations.jl
index 485efde75..6c7a3c799 100644
--- a/src/TensorValues/Operations.jl
+++ b/src/TensorValues/Operations.jl
@@ -209,7 +209,7 @@ end
 
 """
 """
-@generated function symmetic_part(v::TensorValue{D}) where D
+@generated function symmetric_part(v::TensorValue{D}) where D
   str = "("
   for j in 1:D
     for i in 1:D
@@ -222,7 +222,7 @@ end
 
 # Define new operations for Gridap types
 
-for op in (:symmetic_part,)
+for op in (:symmetric_part,)
   @eval begin
     function ($op)(a::GridapType)
       operate($op,a)
diff --git a/src/TensorValues/TensorValues.jl b/src/TensorValues/TensorValues.jl
index 2cac90dc3..eec716873 100644
--- a/src/TensorValues/TensorValues.jl
+++ b/src/TensorValues/TensorValues.jl
@@ -41,7 +41,7 @@ export VectorValue
 export inner, outer, meas
 #export det, inv, tr, dot, norm
 export mutable
-export symmetic_part
+export symmetric_part
 export n_components
 export change_eltype
 export diagonal_tensor
diff --git a/src/Visualization/Vtk.jl b/src/Visualization/Vtk.jl
index a843ffd54..9c44391a6 100644
--- a/src/Visualization/Vtk.jl
+++ b/src/Visualization/Vtk.jl
@@ -399,8 +399,8 @@ end
 function _prepare_pdata(trian,cellfields,samplingpoints)
   pdata = Dict()
   for (k,v) in cellfields
-    _v = convert_to_cell_field(v,get_cell_map(trian))
-    pdata[k], = _prepare_node_to_coords(evaluate_field_array(get_array(_v),samplingpoints))
+    _v = CellField(v,trian)
+    pdata[k], = _prepare_node_to_coords(evaluate(_v,samplingpoints))
   end
   pdata
 end
diff --git a/test/FESpacesTests/CellBasesTests.jl b/test/FESpacesTests/CellBasesTests.jl
index dbb96e652..580305622 100644
--- a/test/FESpacesTests/CellBasesTests.jl
+++ b/test/FESpacesTests/CellBasesTests.jl
@@ -43,8 +43,8 @@ f(x) = sin(4*pi*(x[1]-x[2]^2))+1
 uh = interpolate_everywhere(V0,f)
 uhx = collect(evaluate(uh,q))
 
-v = GenericCellBasis(Val{false}(),get_array(cell_basis),get_cell_map(cell_basis))
-u = GenericCellBasis(Val{true}(),get_array(cell_basis),get_cell_map(cell_basis))
+v = GenericCellBasis(Val{false}(),get_array(cell_basis),get_cell_map(cell_basis),Val{true}())
+u = GenericCellBasis(Val{true}(),get_array(cell_basis),get_cell_map(cell_basis),Val{true}())
 
 dv = v
 du = u
diff --git a/test/FESpacesTests/CellDofBasesTests.jl b/test/FESpacesTests/CellDofBasesTests.jl
new file mode 100644
index 000000000..44cb8976d
--- /dev/null
+++ b/test/FESpacesTests/CellDofBasesTests.jl
@@ -0,0 +1,52 @@
+module CellDofBasesTests
+
+using Test
+using Gridap.Geometry
+using Gridap.FESpaces
+using Gridap.Fields
+
+domain =(0,1,0,1)
+partition = (3,3)
+model = CartesianDiscreteModel(domain,partition)
+
+order = 1
+
+V = TestFESpace(
+  reffe = :Lagrangian,
+  conformity = :H1,
+  valuetype=Float64,
+  order = order,
+  model = model,
+  dirichlet_tags = [1,6],
+  dof_space = :physical)
+
+cell_dof_basis = get_cell_dof_basis(V)
+@test is_in_physical_space(cell_dof_basis)
+
+cell_basis = get_cell_basis(V)
+@test is_in_physical_space(cell_basis)
+
+test_cell_dof_basis(cell_dof_basis,cell_basis)
+
+# TODO not working: it should be the cell-wise identity...
+#display(evaluate(cell_dof_basis,cell_basis))
+
+V = TestFESpace(
+  reffe = :Lagrangian,
+  conformity = :H1,
+  valuetype=Float64,
+  order = order,
+  model = model,
+  dirichlet_tags = [1,6])
+
+cell_dof_basis = get_cell_dof_basis(V)
+@test is_in_ref_space(cell_dof_basis)
+
+cell_basis = get_cell_basis(V)
+@test is_in_ref_space(cell_basis)
+
+test_cell_dof_basis(cell_dof_basis,cell_basis)
+
+#display(evaluate(cell_dof_basis,cell_basis))
+
+end # module
diff --git a/test/FESpacesTests/ConformingFESpacesTests.jl b/test/FESpacesTests/ConformingFESpacesTests.jl
index 0858e38fe..c39da7b30 100644
--- a/test/FESpacesTests/ConformingFESpacesTests.jl
+++ b/test/FESpacesTests/ConformingFESpacesTests.jl
@@ -18,11 +18,11 @@ grid_topology = get_grid_topology(model)
 polytopes = get_polytopes(grid_topology)
 reffes = [LagrangianRefFE(Float64,p,order) for p in polytopes]
 
-face_labeing = get_face_labeling(model)
+face_labeling = get_face_labeling(model)
 dirichlet_tags = ["tag_1","tag_6"]
 
 cell_dofs, nfree, ndiri, dirichlet_dof_tag, dirichlet_cells = compute_conforming_cell_dofs(
-  reffes, grid_topology, face_labeing, dirichlet_tags)
+  reffes, grid_topology, face_labeling, dirichlet_tags)
 
 r = [
   [-1,1,4,5,14,15,16,17,35],[1,2,5,6,18,19,17,20,36],[2,3,6,7,21,22,20,23,37],
@@ -40,7 +40,7 @@ reffes = [LagrangianRefFE(VectorValue{2,Float64},p,order) for p in polytopes]
 dirichlet_components = [(true,true), (false,true)]
 
 cell_dofs, nfree, ndiri, dirichlet_dof_tag, dirichlet_cells = compute_conforming_cell_dofs(
-  reffes, grid_topology, face_labeing, dirichlet_tags, dirichlet_components)
+  reffes, grid_topology, face_labeling, dirichlet_tags, dirichlet_components)
 
 r = [
   [-1,1,7,9,-2,2,8,10],[1,3,9,11,2,4,10,12],[3,5,11,13,4,6,12,14],
@@ -59,7 +59,7 @@ reffes = [LagrangianRefFE(VectorValue{2,Float64},p,order) for p in polytopes]
 dirichlet_components = [(true,true), (false,true)]
 
 cell_dofs, nfree, ndiri, dirichlet_dof_tag, dirichlet_cells = compute_conforming_cell_dofs(
-  reffes, grid_topology, face_labeing, dirichlet_tags, dirichlet_components)
+  reffes, grid_topology, face_labeling, dirichlet_tags, dirichlet_components)
 
 V = GradConformingFESpace(reffes,model,dirichlet_tags)
 test_single_field_fe_space(V)
diff --git a/test/FESpacesTests/CurlConformingFESpacesTests.jl b/test/FESpacesTests/CurlConformingFESpacesTests.jl
index 9b451018e..9a976f248 100644
--- a/test/FESpacesTests/CurlConformingFESpacesTests.jl
+++ b/test/FESpacesTests/CurlConformingFESpacesTests.jl
@@ -1,5 +1,5 @@
 module CurlConformingFESpacesTests
-##
+
 using Test
 using Gridap
 using LinearAlgebra
@@ -43,6 +43,5 @@ el2 = sqrt(sum(integrate(inner(e,e),trian,quad)))
 
 # writevtk(trian,"trian",nsubcells=10,cellfields=["uh"=>uh])
 
-##
 
 end # module
diff --git a/test/FESpacesTests/ExtendedFESpacesTests.jl b/test/FESpacesTests/ExtendedFESpacesTests.jl
index 3efc978aa..ccabd1ba1 100644
--- a/test/FESpacesTests/ExtendedFESpacesTests.jl
+++ b/test/FESpacesTests/ExtendedFESpacesTests.jl
@@ -76,6 +76,7 @@ u(x) = x[1]+x[2]
 
 uh = interpolate(U,u)
 
+
 uh_in = restrict(uh,trian_in)
 
 uh_Γ = restrict(uh,trian_Γ)
diff --git a/test/FESpacesTests/FESpaceFactoriesTests.jl b/test/FESpacesTests/FESpaceFactoriesTests.jl
index 3f71acffa..3fd2f4ff1 100644
--- a/test/FESpacesTests/FESpaceFactoriesTests.jl
+++ b/test/FESpacesTests/FESpaceFactoriesTests.jl
@@ -19,7 +19,7 @@ V = FESpace(
  order=order,
  conformity=:L2)
 
-@test isa(V,UnsconstrainedFESpace)
+@test isa(V,UnconstrainedFESpace)
 
 V = FESpace(
  triangulation=get_triangulation(model),
@@ -28,7 +28,7 @@ V = FESpace(
  order=order,
  conformity=:L2)
 
-@test isa(V,UnsconstrainedFESpace)
+@test isa(V,UnconstrainedFESpace)
 
 V = FESpace(
  model=model,
@@ -37,7 +37,7 @@ V = FESpace(
  order=order,
  conformity=:H1)
 
-@test isa(V,UnsconstrainedFESpace)
+@test isa(V,UnconstrainedFESpace)
 
 V = FESpace(
  model=model,
@@ -47,7 +47,7 @@ V = FESpace(
  order=order,
  dirichlet_tags="boundary")
 
-@test isa(V,UnsconstrainedFESpace)
+@test isa(V,UnconstrainedFESpace)
 
 D = num_point_dims(model)
 
@@ -60,7 +60,7 @@ V = FESpace(
  dirichlet_tags="boundary",
  dirichlet_masks=(true,false))
 
-@test isa(V,UnsconstrainedFESpace)
+@test isa(V,UnconstrainedFESpace)
 
 V = FESpace(
  model=model,
@@ -71,7 +71,7 @@ V = FESpace(
  dirichlet_tags=[1,2],
  dirichlet_masks=[(true,false),(false,true)])
 
-@test isa(V,UnsconstrainedFESpace)
+@test isa(V,UnconstrainedFESpace)
 
 V = FESpace(
  model=model,
@@ -84,7 +84,7 @@ V = FESpace(
 #uh = FEFunction(V,rand(num_free_dofs(V)))
 #writevtk(get_triangulation(model),"trian",nsubcells=20,cellfields=["uh"=>uh])
 
-@test isa(V,UnsconstrainedFESpace)
+@test isa(V,UnconstrainedFESpace)
 
 V = FESpace(
  model=model,
diff --git a/test/FESpacesTests/PhysicalBasesTests.jl b/test/FESpacesTests/PhysicalBasesTests.jl
new file mode 100644
index 000000000..2ce8ecc47
--- /dev/null
+++ b/test/FESpacesTests/PhysicalBasesTests.jl
@@ -0,0 +1,141 @@
+module PhysicalBasesTests
+
+using Gridap
+using Gridap.ReferenceFEs
+using Gridap.Geometry
+using Gridap.Arrays
+using Gridap.Fields
+using Gridap.FESpaces
+using Gridap.Polynomials
+using Gridap.Integration
+using Test
+
+# Start with a PhysicalSpaceCellBasis
+a = 1
+b = ( a == 1 )
+
+# domain = (0,1)
+# partition = (3,)
+domain = (0,1,0,1)
+partition = (2,2)
+model = CartesianDiscreteModel(domain,partition)
+order = 1
+# order = 2
+
+trian = get_triangulation(model)
+quad = CellQuadrature(trian,order)
+q = get_coordinates(quad)
+
+polytopes = get_polytopes(model)
+cell_to_ctype = get_cell_type(model)
+
+grid = get_grid(model)
+cell_map = get_cell_map(grid)
+
+# Test against the ref approach...
+
+T = Float64
+reffes = [LagrangianRefFE(T,p,order) for p in polytopes]
+
+dof_bases = map(get_dof_basis,reffes)
+
+FEM = Gridap.FESpaces
+cell_dof_basis = FEM._cell_dof_basis_physical_space(dof_bases,cell_to_ctype,cell_map)
+cell_dof_basis = GenericCellDofBasis(Val{false}(),cell_dof_basis)
+
+prebasis =  map(get_prebasis,reffes)
+cell_prebasis = CompressedArray(prebasis,cell_to_ctype)
+cell_prebasis_new = GenericCellBasis(Val{false}(),cell_prebasis,cell_map,Val{false}())
+
+typeof(cell_prebasis)
+isa(cell_prebasis,CellBasis)
+
+# cell_matrix = evaluate(cell_dof_basis,cell_prebasis)
+
+# cell_shapefuns = _cell_shape_functions_physical_space(cell_prebasis,cell_dof_basis,cell_map)
+
+psfs, x  = Gridap.FESpaces.compute_cell_space_physical(reffes, cell_to_ctype, cell_map)
+sfs, x  = Gridap.FESpaces.compute_cell_space(reffes, cell_to_ctype, cell_map)
+
+# T = VectorValue{2,Float64}
+# reffes = [LagrangianRefFE(T,p,order) for p in polytopes]
+
+r = evaluate(sfs,q)
+rg = evaluate(gradient(sfs),q)
+rp = evaluate(psfs,q)
+rgp = evaluate(gradient(psfs),q)
+
+@test all([ rg[i] ≈ rgp[i] for i in 1:length(rg) ])
+@test all([ r[i] ≈ rp[i] for i in 1:length(rg) ])
+
+T = VectorValue{2,Float64}
+reffes = [LagrangianRefFE(T,p,order) for p in polytopes]
+
+psfs, x  = Gridap.FESpaces.compute_cell_space_physical(reffes, cell_to_ctype, cell_map)
+sfs, x  = Gridap.FESpaces.compute_cell_space(reffes, cell_to_ctype, cell_map)
+
+r = evaluate(sfs,q)
+rg = evaluate(gradient(sfs),q)
+rp = evaluate(psfs,q)
+rgp = evaluate(gradient(psfs),q)
+
+@test all([ rg[i] ≈ rgp[i] for i in 1:length(rg) ])
+@test all([ r[i] ≈ rp[i] for i in 1:length(rg) ])
+
+func(x) = x
+
+cell_field = convert_to_cell_field(func,cell_map)
+isa(cell_field,CellField)
+
+# Now RT elements
+
+T = Float64
+order = 0
+reffes = [RaviartThomasRefFE(T,p,order) for p in polytopes]
+
+psfs, dofp  = Gridap.FESpaces.compute_cell_space_physical(reffes, cell_to_ctype, cell_map)
+sfs, dof  = Gridap.FESpaces.compute_cell_space(reffes, cell_to_ctype, cell_map)
+
+r = evaluate(sfs,q)
+rg = evaluate(gradient(sfs),q)
+rp = evaluate(psfs,q)
+rgp = evaluate(gradient(psfs),q)
+
+@test all([ r[i] ≈ rp[i] for i in 1:length(rp) ])
+@test all([ rg[i] ≈ rgp[i] for i in 1:length(rg) ])
+
+end #module
+
+# # If I want new evaluation...
+# function kernel_evaluate(k::typeof{change_basis},x,cell_prebasis,cell_matrix_inv)
+#    cell_prebasis_x = evaluate_field_array(cell_prebasis,x)
+#   apply(mul,cell_prebasis_x,cell_prebasis,cell_matrix_inv)
+# end
+# function apply_gradient(k::typeof(change_basis),cell_prebasis,cell_matrix_inv)
+#    cell_prebasis_grad = gradient(cell_prebasis)
+#    apply(change_basis,cell_prebasis_grad,cell_matrix_inv)
+# end
+# ##
+# # Optimisation : evaluate_field_array for AbstractArray with FieldLike
+# # Define a new kernel that better treats the inverse
+# struct InvKernel <: Kernel end
+# function kernel_cache(k::InvKernel,mat)
+# end
+# function apply_kernel!(cache,k::InvKernel,mat)
+# end
+# function kernel_cache(k::InvKernel,mat)
+# CachedArray(copy(mat))
+# end
+# function apply_kernel!(cache,k::InvKernel,mat)
+#   setsize!(cache,size(mat))
+#   m = cache.array
+#   fill!(m,zero(m))
+#   for i:size(m,1); m[i] = 1; end
+#   ldiv!(mat,m)
+#   m
+# end
+# k = InvKernel()
+#
+# isa(cell_prebasis,CellBasis)
+#
+# change_basis(cell_prebasis[1],cell_matrix_inv[1])
diff --git a/test/FESpacesTests/PhysicalFESpacesTests.jl b/test/FESpacesTests/PhysicalFESpacesTests.jl
new file mode 100644
index 000000000..647026504
--- /dev/null
+++ b/test/FESpacesTests/PhysicalFESpacesTests.jl
@@ -0,0 +1,128 @@
+module PhysicalFESpacesTests
+
+using Gridap
+using Gridap.ReferenceFEs
+using Gridap.Geometry
+using Gridap.Arrays
+using Gridap.Fields
+using Gridap.FESpaces
+using Gridap.Polynomials
+using Test
+
+# domain = (0,1)
+# partition = (3,)
+domain = (0,1,0,1)
+partition = (2,2)
+model = CartesianDiscreteModel(domain,partition)
+order = 1
+
+trian = get_triangulation(model)
+quad = CellQuadrature(trian,order)
+
+u(x) = x[1]
+T = Float64
+
+u(x) = VectorValue(x[1],x[2])
+T = VectorValue{2,Float64}
+
+# Juno.@run TestFESpace( reffe = :Lagrangian, conformity = :H1, valuetype=T, order = order, model = model, dirichlet_tags = [1,6], dof_space = :physical)
+# Juno.@enter TestFESpace( reffe = :Lagrangian, conformity = :H1, valuetype=T, order = order, model = model, dirichlet_tags = [1,6], dof_space = :physical)
+# Juno.@enter TestFESpace( reffe = :Lagrangian, conformity = :H1, valuetype=T, order = order, model = model, dirichlet_tags = [1,6])
+
+Vp = TestFESpace(
+  reffe = :Lagrangian,
+  conformity = :H1,
+  valuetype=T,
+  order = order,
+  model = model,
+  dirichlet_tags = [1,6],
+  dof_space = :physical)
+
+V = TestFESpace(
+  reffe = :Lagrangian,
+  conformity = :H1,
+  valuetype=VectorValue{2,Float64},
+  order = order,
+  model = model,
+  dirichlet_tags = [1,6])
+
+  V.cell_basis.ref_trait
+
+# Vp = TestFESpace(
+#   reffe = :RaviartThomas,
+#   conformity = :Hdiv,
+#   order = order,
+#   model = model,
+#   dirichlet_tags = [1,6],
+#   dof_space = :physical)
+#
+# V = TestFESpace(
+#   reffe = :RaviartThomas,
+#   conformity = :Hdiv,
+#   order = order,
+#   model = model,
+#   dirichlet_tags = [1,6])
+# test_single_field_fe_space(V)
+
+Up = TrialFESpace(Vp,u)
+U = TrialFESpace(V,u)
+U.cell_basis.ref_trait
+U.space.cell_basis.ref_trait
+
+# FEM = Gridap.FESpaces
+# Juno.@enter TrialFESpace(Vp,u)
+# dirichlet_values = FEM.compute_dirichlet_values_for_tags(V,u)
+
+q = get_coordinates(quad)
+r = evaluate(U.cell_basis,q)
+rp = evaluate(Up.cell_basis,q)
+
+@test all(collect(r) .≈ collect(rp))
+
+@test Gridap.FESpaces.RefStyle(Up.space.cell_dof_basis) == Val{false}()
+@test Gridap.FESpaces.RefStyle(U.space.cell_dof_basis) == Val{true}()
+
+@test Gridap.FESpaces.RefStyle(Up.space.cell_basis) == Val{false}()
+@test Gridap.FESpaces.RefStyle(U.space.cell_basis) == Val{true}()
+
+@test Gridap.FESpaces.RefStyle(Vp.cell_basis) == Val{false}()
+@test Gridap.FESpaces.RefStyle(V.cell_basis) == Val{true}()
+
+uhf = Gridap.FESpaces.interpolate(Up,u)
+uhd = Gridap.FESpaces.interpolate_dirichlet(Up,u)
+uhdf = Gridap.FESpaces.interpolate_everywhere(Up,u)
+uhf.dirichlet_values
+uhdf.dirichlet_values
+uhd.dirichlet_values
+
+uhf_r = Gridap.FESpaces.interpolate(U,u)
+uhd_r= Gridap.FESpaces.interpolate_dirichlet(U,u)
+uhdf_r = Gridap.FESpaces.interpolate_everywhere(U,u)
+uhd_r.dirichlet_values
+
+
+@test uhf.free_values == uhf_r.free_values
+@test uhf.dirichlet_values == uhf_r.dirichlet_values
+@test uhd.dirichlet_values == uhd_r.dirichlet_values
+@test uhdf.free_values == uhdf_r.free_values
+@test uhdf.dirichlet_values == uhdf_r.dirichlet_values
+
+@test uhf.free_values == uhdf.free_values
+@test uhf.dirichlet_values == uhd.dirichlet_values
+@test uhdf.free_values == uhf.free_values
+@test uhdf.dirichlet_values == uhd.dirichlet_values
+@test uhdf.dirichlet_values == uhdf_r.dirichlet_values
+
+e = u - uhf_r
+
+el2 = sqrt(sum(integrate(inner(e,e),trian,quad)))
+
+@test el2 < 1.0e-10
+
+e = u - uhf
+
+el2 = sqrt(sum(integrate(inner(e,e),trian,quad)))
+
+@test el2 < 1.0e-10
+
+end #module
diff --git a/test/FESpacesTests/UnconstrainedFESpacesTests.jl b/test/FESpacesTests/UnconstrainedFESpacesTests.jl
index a5d41133b..a29d3ee6c 100644
--- a/test/FESpacesTests/UnconstrainedFESpacesTests.jl
+++ b/test/FESpacesTests/UnconstrainedFESpacesTests.jl
@@ -1,4 +1,4 @@
-module UnsconstrainedFESpacesTests
+module UnconstrainedFESpacesTests
 
 using Test
 using Gridap.Arrays
diff --git a/test/FESpacesTests/runtests.jl b/test/FESpacesTests/runtests.jl
index ac4b8f236..57a94e48f 100644
--- a/test/FESpacesTests/runtests.jl
+++ b/test/FESpacesTests/runtests.jl
@@ -4,6 +4,8 @@ using Test
 
 @testset "CellBases" begin include("CellBasesTests.jl") end
 
+@testset "CellDofBases" begin include("CellDofBasesTests.jl") end
+
 @testset "Law" begin include("LawTests.jl") end
 
 @testset "FEFunctions" begin include("FEFunctionsTests.jl") end
@@ -52,6 +54,10 @@ using Test
 
 @testset "ExtendedFESpaces" begin include("ExtendedFESpacesTests.jl") end
 
+@testset "PhysicalBasesTests" begin include("PhysicalBasesTests.jl") end
+
+@testset "PhysicalFESpaces" begin include("PhysicalFESpacesTests.jl") end
+
 @testset "FESpaceFactories" begin include("FESpaceFactoriesTests.jl") end
 
 @testset "StateLaws" begin include("StateLawsTests.jl") end
diff --git a/test/FieldsTests/DiffOperatorsTests.jl b/test/FieldsTests/DiffOperatorsTests.jl
index 2402aacaa..76d86ae5d 100644
--- a/test/FieldsTests/DiffOperatorsTests.jl
+++ b/test/FieldsTests/DiffOperatorsTests.jl
@@ -34,7 +34,7 @@ for f in (_f,_af)
   
   @test outer(f,∇) == transpose(∇(f))
 
-  @test ε(f) == symmetic_part(gradient(f))
+  @test ε(f) == symmetric_part(gradient(f))
 
   @test Δ(f) == ∇*∇(f)
 
diff --git a/test/FieldsTests/FunctionFieldsTests.jl b/test/FieldsTests/FunctionFieldsTests.jl
new file mode 100644
index 000000000..811595143
--- /dev/null
+++ b/test/FieldsTests/FunctionFieldsTests.jl
@@ -0,0 +1,17 @@
+using Gridap
+using Gridap.Inference
+using Gridap.Fields
+using Test
+
+u(x) = x
+
+p1 = Point(1,2)
+p2 = Point(2,1)
+p3 = Point(4,3)
+p4 = Point(6,1)
+x = [p1,p2,p3,p4]
+
+f = FunctionField(u)
+test_field(f,x,x)
+
+Fill(f,4)
diff --git a/test/GeometryTests/CellFieldsTests.jl b/test/GeometryTests/CellFieldsTests.jl
index 531937e05..7304a0d7a 100644
--- a/test/GeometryTests/CellFieldsTests.jl
+++ b/test/GeometryTests/CellFieldsTests.jl
@@ -29,6 +29,9 @@ r = [[0.0], [2.0], [0.0], [2.0]]
 g = Vector{VectorValue{2,Float64}}[[(1, 0)], [(1, 0)], [(1, 0)], [(1, 0)]]
 test_cell_field(cf1,q,r,grad=g)
 
+isa(cf1,Geometry.CellFieldLike)
+@test RefStyle(cf1) == Val{true}()
+
 cf2 = cf1 + 4
 r2 = map( (i) -> i.+4, r)
 test_cell_field(cf2,q,r2)
diff --git a/test/GridapTests/PhysicalPoissonTests.jl b/test/GridapTests/PhysicalPoissonTests.jl
new file mode 100644
index 000000000..87dad1763
--- /dev/null
+++ b/test/GridapTests/PhysicalPoissonTests.jl
@@ -0,0 +1,111 @@
+module PhysicalPoissonTests
+
+using Test
+using Gridap
+import Gridap: ∇
+
+domain = (0,1,0,1)
+partition = (4,4)
+model = CartesianDiscreteModel(domain,partition)
+order = 2
+
+const h = (domain[2]-domain[1]) / partition[1]
+const γ = 10
+
+labels = get_face_labeling(model)
+add_tag_from_tags!(labels,"dirichlet",[1,2,5])
+add_tag_from_tags!(labels,"neumann",[7,8])
+add_tag_from_tags!(labels,"nitsche",6)
+
+trian = get_triangulation(model)
+degree = order
+quad = CellQuadrature(trian,degree)
+
+ntrian = BoundaryTriangulation(model,labels,"neumann")
+ndegree = order
+nquad = CellQuadrature(ntrian,ndegree)
+const nn = get_normal_vector(ntrian)
+
+dtrian = BoundaryTriangulation(model,labels,"nitsche")
+ddegree = order
+dquad = CellQuadrature(dtrian,ddegree)
+const dn = get_normal_vector(dtrian)
+
+u_scal(x) = x[1]^2 + x[2]
+∇u_scal(x) = VectorValue( 2*x[1], one(x[2]) )
+Δu_scal(x) = 2
+f_scal(x) = - Δu_scal(x)
+∇(::typeof(u_scal)) = ∇u_scal
+
+scalar_data = Dict{Symbol,Any}()
+scalar_data[:valuetype] = Float64
+scalar_data[:u] = u_scal
+scalar_data[:f] = f_scal
+
+u_vec(x) = VectorValue( x[1]^2 + x[2], 4*x[1] - x[2]^2 )
+∇u_vec(x) = TensorValue( 2*x[1], one(x[2]), 4*one(x[1]), - 2*x[2] )
+Δu_vec(x) = VectorValue( 2, -2 )
+f_vec(x) = - Δu_vec(x)
+∇(::typeof(u_vec)) = ∇u_vec
+
+vector_data = Dict{Symbol,Any}()
+vector_data[:valuetype] = VectorValue{2,Float64}
+vector_data[:u] = u_vec
+vector_data[:f] = f_vec
+
+for data in [ vector_data, scalar_data ]
+
+  T = data[:valuetype]
+  u = data[:u]
+  f = data[:f]
+
+  V = TestFESpace(
+   model=model,
+   order=order,
+   reffe=:Lagrangian,
+   labels=labels,
+   valuetype=T,
+   dirichlet_tags="dirichlet",
+   dof_space = :physical)
+   # dirichlet_tags="dirichlet")
+  # Juno.@enter TestFESpace( model=model, order=order, reffe=:Lagrangian, labels=labels, valuetype=T, dirichlet_tags="dirichlet", dof_space = :physical)
+
+  U = TrialFESpace(V,u)
+
+  uh = interpolate(U,u)
+
+  a(u,v) = inner(∇(v),∇(u))
+  l(v) = v*f
+  t_Ω = AffineFETerm(a,l,trian,quad)
+
+  uh_Γn = restrict(uh,ntrian)
+  uh_Γd = restrict(uh,dtrian)
+
+  l_Γn(v) = v*(nn*∇(uh_Γn))
+  t_Γn = FESource(l_Γn,ntrian,nquad)
+
+  a_Γd(u,v) = (γ/h)*v*u  - v*(dn*∇(u)) - (dn*∇(v))*u
+  l_Γd(v) = (γ/h)*v*uh_Γd - (dn*∇(v))*u
+  t_Γd = AffineFETerm(a_Γd,l_Γd,dtrian,dquad)
+
+  op = AffineFEOperator(U,V,t_Ω,t_Γn,t_Γd)
+
+  uh = solve(op)
+
+  e = u - uh
+
+  l2(u) = inner(u,u)
+  sh1(u) = a(u,u)
+  h1(u) = sh1(u) + l2(u)
+
+  el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
+  eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))
+  ul2 = sqrt(sum( integrate(l2(uh),trian,quad) ))
+  uh1 = sqrt(sum( integrate(h1(uh),trian,quad) ))
+
+  @test el2/ul2 < 1.e-8
+  @test eh1/uh1 < 1.e-7
+
+end
+
+end # module
diff --git a/test/GridapTests/StokesTaylorHoodTests.jl b/test/GridapTests/StokesTaylorHoodTests.jl
index db1935a23..13fde5591 100644
--- a/test/GridapTests/StokesTaylorHoodTests.jl
+++ b/test/GridapTests/StokesTaylorHoodTests.jl
@@ -29,73 +29,80 @@ add_tag_from_tags!(labels,"neumann",[6,7,8])
 
 order = 2
 
-V = TestFESpace(
- model=model,
- order=order,
- reffe=:Lagrangian,
- labels=labels,
- valuetype=VectorValue{2,Float64},
- dirichlet_tags="dirichlet",
- conformity=:H1)
-
-Q = TestFESpace(
- model=model,
- order=order-1,
- reffe=:Lagrangian,
- valuetype=Float64,
- conformity=:H1)
-
-U = TrialFESpace(V,u)
-P = TrialFESpace(Q)
-
-Y = MultiFieldFESpace([V,Q])
-X = MultiFieldFESpace([U,P])
-
-trian = get_triangulation(model)
-degree = order
-quad = CellQuadrature(trian,degree)
-
-btrian = BoundaryTriangulation(model,labels,"neumann")
-bdegree = order
-bquad = CellQuadrature(btrian,bdegree)
-const n = get_normal_vector(btrian)
-
-function a(x,y)
-  u,p = x
-  v,q = y
-  inner(∇(v),∇(u)) - (∇*v)*p + q*(∇*u)
+ref_style = [:reference,:physical]
+
+for ref_st in ref_style
+  V = TestFESpace(
+  model=model,
+  order=order,
+  reffe=:Lagrangian,
+  labels=labels,
+  valuetype=VectorValue{2,Float64},
+  dirichlet_tags="dirichlet",
+  dof_space=ref_st,
+  conformity=:H1)
+
+  Q = TestFESpace(
+  model=model,
+  order=order-1,
+  reffe=:Lagrangian,
+  valuetype=Float64,
+  dof_space=ref_st,
+  conformity=:H1)
+
+  U = TrialFESpace(V,u)
+  P = TrialFESpace(Q)
+
+  Y = MultiFieldFESpace([V,Q])
+  X = MultiFieldFESpace([U,P])
+
+  trian = get_triangulation(model)
+  degree = order
+  quad = CellQuadrature(trian,degree)
+
+  btrian = BoundaryTriangulation(model,labels,"neumann")
+  bdegree = order
+  bquad = CellQuadrature(btrian,bdegree)
+  n = get_normal_vector(btrian)
+
+  function a(x,y)
+    u,p = x
+    v,q = y
+    inner(∇(v),∇(u)) - (∇*v)*p + q*(∇*u)
+  end
+
+  function l(y)
+    v,q = y
+    v*f + q*g
+  end
+
+  function l_Γb(y)
+    v,q = y
+    v*(n*∇u) - (n*v)*p
+  end
+
+  t_Ω = AffineFETerm(a,l,trian,quad)
+  t_Γb = FESource(l_Γb,btrian,bquad)
+
+  op = AffineFEOperator(X,Y,t_Ω,t_Γb)
+
+  uh, ph = solve(op)
+
+  eu = u - uh
+  ep = p - ph
+
+  l2(v) = v*v
+  h1(v) = v*v + inner(∇(v),∇(v))
+
+  eu_l2 = sqrt(sum(integrate(l2(eu),trian,quad)))
+  eu_h1 = sqrt(sum(integrate(h1(eu),trian,quad)))
+  ep_l2 = sqrt(sum(integrate(l2(ep),trian,quad)))
+
+  tol = 1.0e-9
+  @test eu_l2 < tol
+  @test eu_h1 < tol
+  @test ep_l2 < tol
 end
 
-function l(y)
-  v,q = y
-  v*f + q*g
-end
-
-function l_Γb(y)
-  v,q = y
-  v*(n*∇u) - (n*v)*p
-end
-
-t_Ω = AffineFETerm(a,l,trian,quad)
-t_Γb = FESource(l_Γb,btrian,bquad)
-
-op = AffineFEOperator(X,Y,t_Ω,t_Γb)
-
-uh, ph = solve(op)
-
-eu = u - uh
-ep = p - ph
-
-l2(v) = v*v
-h1(v) = v*v + inner(∇(v),∇(v))
-
-eu_l2 = sqrt(sum(integrate(l2(eu),trian,quad)))
-eu_h1 = sqrt(sum(integrate(h1(eu),trian,quad)))
-ep_l2 = sqrt(sum(integrate(l2(ep),trian,quad)))
-
-tol = 1.0e-9
-@test eu_l2 < tol
-@test eu_h1 < tol
-@test ep_l2 < tol
 
 end # module
diff --git a/test/GridapTests/runtests.jl b/test/GridapTests/runtests.jl
index 4f6421deb..c9e11e633 100644
--- a/test/GridapTests/runtests.jl
+++ b/test/GridapTests/runtests.jl
@@ -20,4 +20,6 @@ using Test
 
 @testset "IsotropicDamage" begin include("IsotropicDamageTests.jl") end
 
+@testset "PhysicalPoisson" begin include("PhysicalPoissonTests.jl") end
+
 end # module
diff --git a/test/MultiFieldTests/MultiFieldCellBasesTests.jl b/test/MultiFieldTests/MultiFieldCellBasesTests.jl
index 94b29f062..7aad5392d 100644
--- a/test/MultiFieldTests/MultiFieldCellBasesTests.jl
+++ b/test/MultiFieldTests/MultiFieldCellBasesTests.jl
@@ -22,12 +22,14 @@ quad = CellQuadrature(trian,degree)
 q = get_coordinates(quad)
 
 V = TestFESpace(model=model,reffe=:Lagrangian,order=order,conformity=:H1,valuetype=Float64)
+Vp = TestFESpace(model=model,reffe=:Lagrangian,order=order,conformity=:H1,valuetype=Float64,dof_space=:physical)
 U = TrialFESpace(V)
 
 cell_basis_v = get_cell_basis(V)
 field_id_v = 3
 v = CellBasisWithFieldID(cell_basis_v,field_id_v)
 vq = collect(evaluate(v,q))
+@test is_in_ref_space(v)
 
 field_id_a = 2
 a = CellBasisWithFieldID(cell_basis_v,field_id_a)
@@ -42,6 +44,12 @@ field_id_b = 5
 b = CellBasisWithFieldID(cell_basis_u,field_id_b)
 bq = collect(evaluate(b,q))
 
+cell_basis_vp = get_cell_basis(Vp)
+field_id_vp = 6
+vp = CellBasisWithFieldID(cell_basis_vp,field_id_vp)
+vpq = collect(evaluate(vp,q))
+@test !is_in_ref_space(vp)
+
 r = 2*v
 test_cell_basis(r,q,2*vq)
 @test isa(r,CellBasisWithFieldID)
diff --git a/test/MultiFieldTests/MultiFieldFESpacesTests.jl b/test/MultiFieldTests/MultiFieldFESpacesTests.jl
index f2edd5f79..3a37ab2e2 100644
--- a/test/MultiFieldTests/MultiFieldFESpacesTests.jl
+++ b/test/MultiFieldTests/MultiFieldFESpacesTests.jl
@@ -10,7 +10,7 @@ using Test
 using Gridap.MultiField
 using Gridap.MultiField: MultiFieldFESpace
 using Gridap.MultiField: MultiFieldCellArray
-using Gridap.MultiField: ConsequtiveMultiFieldStyle
+using Gridap.MultiField: ConsecutiveMultiFieldStyle
 
 order = 2
 
@@ -22,64 +22,69 @@ trian = get_triangulation(model)
 degree = order
 quad = CellQuadrature(trian,degree)
 
-V = TestFESpace(model=model,order=order,reffe=:Lagrangian,conformity=:H1,valuetype=Float64)
-Q = TestFESpace(model=model,order=order-1,reffe=:Lagrangian,conformity=:L2,valuetype=Float64)
+ref_style = [:reference,:physical]
 
-U = TrialFESpace(V)
-P = TrialFESpace(Q)
+for ref_st in ref_style
+  
+  V = TestFESpace(model=model,order=order,reffe=:Lagrangian,conformity=:H1,valuetype=Float64,dof_space=ref_st)
+  Q = TestFESpace(model=model,order=order-1,reffe=:Lagrangian,conformity=:L2,valuetype=Float64,dof_space=ref_st)
 
-multi_field_style = ConsequtiveMultiFieldStyle()
+  U = TrialFESpace(V)
+  P = TrialFESpace(Q)
 
-Y = MultiFieldFESpace([V,Q],multi_field_style)
-X = MultiFieldFESpace([U,P],multi_field_style)
+  multi_field_style = ConsecutiveMultiFieldStyle()
 
-@test num_free_dofs(X) == num_free_dofs(U) + num_free_dofs(P)
-@test num_free_dofs(X) == num_free_dofs(Y)
+  Y = MultiFieldFESpace([V,Q],multi_field_style)
+  X = MultiFieldFESpace([U,P],multi_field_style)
 
-free_values = rand(num_free_dofs(X))
-xh = FEFunction(X,free_values)
-test_fe_function(xh)
-@test is_a_fe_function(xh)
-uh, ph = xh
-@test is_a_fe_function(uh)
-@test is_a_fe_function(ph)
+  @test num_free_dofs(X) == num_free_dofs(U) + num_free_dofs(P)
+  @test num_free_dofs(X) == num_free_dofs(Y)
 
-dy = get_cell_basis(Y)
-@test is_test(dy)
-@test is_a_fe_cell_basis(dy)
-dv, dq = dy
-@test is_a_fe_cell_basis(dv)
-@test is_a_fe_cell_basis(dq)
+  free_values = rand(num_free_dofs(X))
+  xh = FEFunction(X,free_values)
+  test_fe_function(xh)
+  @test is_a_fe_function(xh)
+  uh, ph = xh
+  @test is_a_fe_function(uh)
+  @test is_a_fe_function(ph)
 
-dx = get_cell_basis(X)
-@test is_a_fe_cell_basis(dx)
-@test is_trial(dx)
-du, dp = dx
-@test is_a_fe_cell_basis(du)
-@test is_a_fe_cell_basis(dp)
+  dy = get_cell_basis(Y)
+  @test is_test(dy)
+  @test is_a_fe_cell_basis(dy)
+  dv, dq = dy
+  @test is_a_fe_cell_basis(dv)
+  @test is_a_fe_cell_basis(dq)
 
-cellmat = integrate(dv*du,trian,quad)
-cellvec = integrate(dv*2,trian,quad)
-cellids = get_cell_id(trian)
-cellmatvec = pair_arrays(cellmat,cellvec)
+  dx = get_cell_basis(X)
+  @test is_a_fe_cell_basis(dx)
+  @test is_trial(dx)
+  du, dp = dx
+  @test is_a_fe_cell_basis(du)
+  @test is_a_fe_cell_basis(dp)
 
-matvecdata = (cellmatvec,cellids,cellids)
-matdata = (cellmat,cellids,cellids)
-vecdata = (cellvec,cellids)
+  cellmat = integrate(dv*du,trian,quad)
+  cellvec = integrate(dv*2,trian,quad)
+  cellids = get_cell_id(trian)
+  cellmatvec = pair_arrays(cellmat,cellvec)
 
-test_fe_space(V,matvecdata,matdata,vecdata)
-test_fe_space(U,matvecdata,matdata,vecdata)
+  matvecdata = (cellmatvec,cellids,cellids)
+  matdata = (cellmat,cellids,cellids)
+  vecdata = (cellvec,cellids)
 
-#using Gridap.Visualization
-#writevtk(trian,"trian";nsubcells=30,cellfields=["uh" => uh, "ph"=> ph])
+  test_fe_space(V,matvecdata,matdata,vecdata)
+  test_fe_space(U,matvecdata,matdata,vecdata)
 
-cell_dofs = get_cell_dofs(X)
-@test isa(cell_dofs,MultiFieldCellArray)
+  #using Gridap.Visualization
+  #writevtk(trian,"trian";nsubcells=30,cellfields=["uh" => uh, "ph"=> ph])
 
-cellids = [3,5,2]
+  cell_dofs = get_cell_dofs(X)
+  @test isa(cell_dofs,MultiFieldCellArray)
 
-cell_dofs_new = reindex(cell_dofs,cellids)
-@test isa(cell_dofs_new,MultiFieldCellArray)
-@test cell_dofs_new.block_ids == cell_dofs.block_ids
+  cellids = [3,5,2]
+
+  cell_dofs_new = reindex(cell_dofs,cellids)
+  @test isa(cell_dofs_new,MultiFieldCellArray)
+  @test cell_dofs_new.block_ids == cell_dofs.block_ids
+end
 
 end # module
diff --git a/test/ReferenceFEsTests/NodalReferenceFEsTests.jl b/test/ReferenceFEsTests/NodalReferenceFEsTests.jl
index 2b40f8562..68e6057cb 100644
--- a/test/ReferenceFEsTests/NodalReferenceFEsTests.jl
+++ b/test/ReferenceFEsTests/NodalReferenceFEsTests.jl
@@ -32,7 +32,7 @@ face_own_nodes = face_own_dofs
 face_own_nodes_permutations = face_own_dofs_permutations
 face_nodes = face_dofs
 
-reffe = GenericNodalRefFE(
+reffe = GenericNodalCartesianRefFE(
   reffe,
   node_coordinates,
   node_and_comp_to_dof,
diff --git a/test/TensorValuesTests/OperationsTests.jl b/test/TensorValuesTests/OperationsTests.jl
index d4b81578f..4ac20a2f0 100644
--- a/test/TensorValuesTests/OperationsTests.jl
+++ b/test/TensorValuesTests/OperationsTests.jl
@@ -251,7 +251,7 @@ t = TensorValue(1,2,3,4,5,6,7,8,9)
 @test tr(t) == 15
 @test tr(t) == 15
 
-@test symmetic_part(t) == TensorValue(1.0, 3.0, 5.0, 3.0, 5.0, 7.0, 5.0, 7.0, 9.0)
+@test symmetric_part(t) == TensorValue(1.0, 3.0, 5.0, 3.0, 5.0, 7.0, 5.0, 7.0, 9.0)
 
 a = TensorValue(1,2,3,4)
 b = a'