diff --git a/src/Adaptivity/FineToCoarseReferenceFEs.jl b/src/Adaptivity/FineToCoarseReferenceFEs.jl
index e921d479e..1899d38a0 100644
--- a/src/Adaptivity/FineToCoarseReferenceFEs.jl
+++ b/src/Adaptivity/FineToCoarseReferenceFEs.jl
@@ -45,8 +45,8 @@ function Arrays.evaluate!(cache,s::FineToCoarseDofBasis{T,<:LagrangianDofBasis},
   vals  = evaluate!(cf,field,b.nodes,s.child_ids)
   ndofs = length(b.dof_to_node)
   T2    = eltype(vals)
-  ncomps = num_components(T2)
-  @check ncomps == num_components(eltype(b.node_and_comp_to_dof)) """\n
+  ncomps = num_indep_components(T2) # use num_indep_components ?
+  @check ncomps == num_indep_components(eltype(b.node_and_comp_to_dof)) """\n
   Unable to evaluate LagrangianDofBasis. The number of components of the
   given Field does not match with the LagrangianDofBasis.
 
@@ -80,8 +80,8 @@ end
 
 
 """
-  Wrapper for a ReferenceFE which is specialised for 
-  efficiently evaluating FineToCoarseFields. 
+  Wrapper for a ReferenceFE which is specialised for
+  efficiently evaluating FineToCoarseFields.
 """
 struct FineToCoarseRefFE{T,D,A} <: ReferenceFE{D}
   reffe     :: T
@@ -132,4 +132,4 @@ function FESpaces.TestFESpace(model::DiscreteModel,rrules::AbstractVector{<:Refi
   basis, reffe_args, reffe_kwargs = reffe
   reffes = lazy_map(rr -> ReferenceFE(get_polytope(rr),rr,basis,reffe_args...;reffe_kwargs...),rrules)
   return TestFESpace(model,reffes;kwargs...)
-end
\ No newline at end of file
+end
diff --git a/src/Exports.jl b/src/Exports.jl
index 782625bba..27dbe467d 100644
--- a/src/Exports.jl
+++ b/src/Exports.jl
@@ -41,6 +41,7 @@ using Gridap.Arrays: ∑; export ∑
 @publish TensorValues outer
 @publish TensorValues diagonal_tensor
 @publish TensorValues num_components
+@publish TensorValues num_indep_components
 using Gridap.TensorValues: ⊙; export ⊙
 using Gridap.TensorValues: ⊗; export ⊗
 
diff --git a/src/FESpaces/CLagrangianFESpaces.jl b/src/FESpaces/CLagrangianFESpaces.jl
index c3ad07e1a..f90e58691 100644
--- a/src/FESpaces/CLagrangianFESpaces.jl
+++ b/src/FESpaces/CLagrangianFESpaces.jl
@@ -119,7 +119,7 @@ end
 # Helpers
 
 _default_mask(::Type) = true
-_default_mask(::Type{T}) where T <: MultiValue = ntuple(i->true,Val{length(T)}())
+_default_mask(::Type{T}) where T <: MultiValue = ntuple(i->true,Val{num_indep_components(T)}())
 
 _dof_type(::Type{T}) where T = T
 _dof_type(::Type{T}) where T<:MultiValue = eltype(T)
@@ -193,7 +193,7 @@ function _generate_node_to_dof_glue_component_major(
   z::MultiValue,node_to_tag,tag_to_masks)
   nfree_dofs = 0
   ndiri_dofs = 0
-  ncomps = length(z)
+  ncomps = num_indep_components(z)
   @check length(testitem(tag_to_masks)) == ncomps
   for (node,tag) in enumerate(node_to_tag)
     if tag == UNSET
@@ -218,7 +218,7 @@ function _generate_node_to_dof_glue_component_major(
   node_and_comp_to_dof = zeros(T,nnodes)
   nfree_dofs = 0
   ndiri_dofs = 0
-  m = zero(Mutable(T))
+  m = zeros(Int32, ncomps)
   for (node,tag) in enumerate(node_to_tag)
     if tag == UNSET
       for comp in 1:ncomps
@@ -245,7 +245,7 @@ function _generate_node_to_dof_glue_component_major(
         end
       end
     end
-    node_and_comp_to_dof[node] = m
+    node_and_comp_to_dof[node] = T(m...)
   end
   glue = NodeToDofGlue(
    free_dof_to_node,
@@ -301,7 +301,7 @@ function _generate_cell_dofs_clagrangian(
   cell_to_ctype,
   node_and_comp_to_dof)
 
-  ncomps = num_components(z)
+  ncomps = num_indep_components(z)
 
   ctype_to_lnode_to_comp_to_ldof = map(get_node_and_comp_to_dof,ctype_to_reffe)
   ctype_to_num_ldofs = map(num_dofs,ctype_to_reffe)
@@ -353,8 +353,8 @@ function _fill_cell_dofs_clagrangian!(
     p = cell_to_dofs.ptrs[cell]-1
     for (lnode, node) in enumerate(nodes)
       for comp in 1:ncomps
-        ldof = lnode_and_comp_to_ldof[lnode][comp]
-        dof = node_and_comp_to_dof[node][comp]
+        ldof = indep_comp_getindex(lnode_and_comp_to_ldof[lnode], comp)
+        dof =  indep_comp_getindex(node_and_comp_to_dof[node], comp)
         cell_to_dofs.data[p+ldof] = dof
       end
     end
diff --git a/src/FESpaces/ConstantFESpaces.jl b/src/FESpaces/ConstantFESpaces.jl
index d938cc9e0..f7115f46c 100644
--- a/src/FESpaces/ConstantFESpaces.jl
+++ b/src/FESpaces/ConstantFESpaces.jl
@@ -32,7 +32,7 @@ struct ConstantFESpace{V,T,A,B,C} <: SingleFieldFESpace
     cell_dof_basis_array = lazy_map(get_dof_basis,cell_reffe)
     cell_dof_basis = CellDof(cell_dof_basis_array,Triangulation(model),ReferenceDomain())
 
-    cell_dof_ids = Fill(Int32(1):Int32(num_components(field_type)),num_cells(model))
+    cell_dof_ids = Fill(Int32(1):Int32(num_indep_components(field_type)),num_cells(model))
     A = typeof(cell_basis)
     B = typeof(cell_dof_basis)
     C = typeof(cell_dof_ids)
diff --git a/src/Polynomials/JacobiPolynomialBases.jl b/src/Polynomials/JacobiPolynomialBases.jl
index 2b615d9c3..896fc9c2d 100644
--- a/src/Polynomials/JacobiPolynomialBases.jl
+++ b/src/Polynomials/JacobiPolynomialBases.jl
@@ -9,7 +9,7 @@ struct JacobiPolynomialBasis{D,T} <: AbstractVector{JacobiPolynomial}
   end
 end
 
-@inline Base.size(a::JacobiPolynomialBasis{D,T}) where {D,T} = (length(a.terms)*num_components(T),)
+@inline Base.size(a::JacobiPolynomialBasis{D,T}) where {D,T} = (length(a.terms)*num_indep_components(T),)
 @inline Base.getindex(a::JacobiPolynomialBasis,i::Integer) = JacobiPolynomial()
 @inline Base.IndexStyle(::JacobiPolynomialBasis) = IndexLinear()
 
@@ -49,7 +49,7 @@ return_type(::JacobiPolynomialBasis{D,T}) where {D,T} = T
 function return_cache(f::JacobiPolynomialBasis{D,T},x::AbstractVector{<:Point}) where {D,T}
   @check D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   r = CachedArray(zeros(T,(np,ndof)))
   v = CachedArray(zeros(T,(ndof,)))
@@ -60,7 +60,7 @@ end
 function evaluate!(cache,f::JacobiPolynomialBasis{D,T},x::AbstractVector{<:Point}) where {D,T}
   r, v, c = cache
   np = length(x)
-  ndof = length(f.terms)*num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -82,7 +82,7 @@ function return_cache(
   f = fg.fa
   @assert D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(V)
+  ndof = length(f)
   xi = testitem(x)
   T = gradient_type(V,xi)
   n = 1 + _maximum(f.orders)
@@ -101,7 +101,7 @@ function evaluate!(
   f = fg.fa
   r, v, c, g = cache
   np = length(x)
-  ndof = length(f.terms) * num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -124,7 +124,7 @@ function return_cache(
   f = fg.fa
   @assert D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(V)
+  ndof = length(f)
   xi = testitem(x)
   T = gradient_type(gradient_type(V,xi),xi)
   n = 1 + _maximum(f.orders)
@@ -144,7 +144,7 @@ function evaluate!(
   f = fg.fa
   r, v, c, g, h = cache
   np = length(x)
-  ndof = length(f.terms) * num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -164,7 +164,7 @@ end
 # Optimizing evaluation at a single point
 
 function return_cache(f::JacobiPolynomialBasis{D,T},x::Point) where {D,T}
-  ndof = length(f.terms)*num_components(T)
+  ndof = length(f)
   r = CachedArray(zeros(T,(ndof,)))
   xs = [x]
   cf = return_cache(f,xs)
diff --git a/src/Polynomials/ModalC0Bases.jl b/src/Polynomials/ModalC0Bases.jl
index 9c3a6b2e9..834b742cc 100644
--- a/src/Polynomials/ModalC0Bases.jl
+++ b/src/Polynomials/ModalC0Bases.jl
@@ -11,11 +11,12 @@ struct ModalC0Basis{D,T,V} <: AbstractVector{ModalC0BasisFunction}
     terms::Vector{CartesianIndex{D}},
     a::Vector{Point{D,V}},
     b::Vector{Point{D,V}}) where {D,T,V}
+
     new{D,T,V}(orders,terms,a,b)
   end
 end
 
-@inline Base.size(a::ModalC0Basis{D,T,V}) where {D,T,V} = (length(a.terms)*num_components(T),)
+@inline Base.size(a::ModalC0Basis{D,T,V}) where {D,T,V} = (length(a.terms)*num_indep_components(T),)
 @inline Base.getindex(a::ModalC0Basis,i::Integer) = ModalC0BasisFunction()
 @inline Base.IndexStyle(::ModalC0Basis) = IndexLinear()
 
@@ -101,7 +102,7 @@ return_type(::ModalC0Basis{D,T,V}) where {D,T,V} = T
 function return_cache(f::ModalC0Basis{D,T,V},x::AbstractVector{<:Point}) where {D,T,V}
   @assert D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   r = CachedArray(zeros(T,(np,ndof)))
   v = CachedArray(zeros(T,(ndof,)))
@@ -112,7 +113,7 @@ end
 function evaluate!(cache,f::ModalC0Basis{D,T,V},x::AbstractVector{<:Point}) where {D,T,V}
   r, v, c = cache
   np = length(x)
-  ndof = length(f.terms)*num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -134,7 +135,7 @@ function return_cache(
   f = fg.fa
   @assert D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(V)
+  ndof = length(f)
   xi = testitem(x)
   T = gradient_type(V,xi)
   n = 1 + _maximum(f.orders)
@@ -153,7 +154,7 @@ function evaluate!(
   f = fg.fa
   r, v, c, g = cache
   np = length(x)
-  ndof = length(f.terms) * num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -176,7 +177,7 @@ function return_cache(
   f = fg.fa
   @assert D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(V)
+  ndof = length(f)
   xi = testitem(x)
   T = gradient_type(gradient_type(V,xi),xi)
   n = 1 + _maximum(f.orders)
@@ -196,7 +197,7 @@ function evaluate!(
   f = fg.fa
   r, v, c, g, h = cache
   np = length(x)
-  ndof = length(f.terms) * num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -391,16 +392,17 @@ function _evaluate_nd_mc0!(
 end
 
 @inline function _set_value_mc0!(v::AbstractVector{V},s::T,k,l) where {V,T}
-  m = zero(Mutable(V))
+  ncomp = num_indep_components(V)
+  m = zeros(T,ncomp)
   z = zero(T)
-  js = eachindex(m)
+  js = 1:ncomp
   for j in js
     for i in js
       @inbounds m[i] = z
     end
     @inbounds m[j] = s
     i = k+l*(j-1)
-    @inbounds v[i] = m
+    @inbounds v[i] = V(m...)
   end
   k+1
 end
@@ -461,8 +463,12 @@ end
   k+1
 end
 
+# Indexing and m definition should be fixed if G contains symmetries, that is
+# if the code is  optimized for symmetric tensor V valued FESpaces
+# (if gradient_type(V) returned a symmetric higher order tensor type G)
 @inline function _set_gradient_mc0!(
   v::AbstractVector{G},s,k,l,::Type{V}) where {V,G}
+  @notimplementedif num_indep_components(G) != num_components(G) "Not implemented for symmetric Jacobian or Hessian"
 
   T = eltype(s)
   m = zero(Mutable(G))
diff --git a/src/Polynomials/MonomialBases.jl b/src/Polynomials/MonomialBases.jl
index 4f5238983..261ffb809 100644
--- a/src/Polynomials/MonomialBases.jl
+++ b/src/Polynomials/MonomialBases.jl
@@ -18,7 +18,7 @@ struct MonomialBasis{D,T} <: AbstractVector{Monomial}
   end
 end
 
-Base.size(a::MonomialBasis{D,T}) where {D,T} = (length(a.terms)*num_components(T),)
+Base.size(a::MonomialBasis{D,T}) where {D,T} = (length(a.terms)*num_indep_components(T),)
 # @santiagobadia : Not sure we want to create the monomial machinery
 Base.getindex(a::MonomialBasis,i::Integer) = Monomial()
 Base.IndexStyle(::MonomialBasis) = IndexLinear()
@@ -122,7 +122,7 @@ function return_cache(f::MonomialBasis{D,T},x::AbstractVector{<:Point}) where {D
   zxi = zero(eltype(eltype(x)))
   Tp = typeof( zT*zxi*zxi + zT*zxi*zxi  )
   np = length(x)
-  ndof = length(f.terms)*num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   r = CachedArray(zeros(Tp,(np,ndof)))
   v = CachedArray(zeros(Tp,(ndof,)))
@@ -133,7 +133,7 @@ end
 function evaluate!(cache,f::MonomialBasis{D,T},x::AbstractVector{<:Point}) where {D,T}
   r, v, c = cache
   np = length(x)
-  ndof = length(f.terms)*num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -157,7 +157,7 @@ function _return_cache(
   f = fg.fa
   @check D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(V)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   r = CachedArray(zeros(T,(np,ndof)))
   v = CachedArray(zeros(T,(ndof,)))
@@ -173,7 +173,7 @@ function _return_cache(
   TisbitsType::Val{false}) where {D,V,T}
 
   cache = _return_cache(fg,x,T,Val{true}())
-  z = CachedArray(zeros(eltype(T),D)) 
+  z = CachedArray(zeros(eltype(T),D))
   (cache...,z)
 end
 
@@ -197,7 +197,7 @@ function _evaluate!(
   r, v, c, g = cache
   z = zero(Mutable(VectorValue{D,eltype(T)}))
   np = length(x)
-  ndof = length(f.terms) * num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -222,7 +222,7 @@ function _evaluate!(
   f = fg.fa
   r, v, c, g, z = cache
   np = length(x)
-  ndof = length(f.terms) * num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -255,7 +255,7 @@ function return_cache(
   f = fg.fa
   @check D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(V)
+  ndof = length(f)
   xi = testitem(x)
   T = gradient_type(gradient_type(V,xi),xi)
   n = 1 + _maximum(f.orders)
@@ -275,7 +275,7 @@ function evaluate!(
   f = fg.fa
   r, v, c, g, h = cache
   np = length(x)
-  ndof = length(f.terms) * num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -406,15 +406,16 @@ function _evaluate_nd!(
 end
 
 function _set_value!(v::AbstractVector{V},s::T,k) where {V,T}
-  m = zero(Mutable(V))
+  ncomp = num_indep_components(V)
+  m = zeros(T,ncomp)
   z = zero(T)
-  js = eachindex(m)
+  js = 1:ncomp
   for j in js
     for i in js
       @inbounds m[i] = z
     end
     m[j] = s
-    v[k] = m
+    v[k] = V(m...)
     k += 1
   end
   k
@@ -440,7 +441,7 @@ function _gradient_nd!(
     _evaluate_1d!(c,x,orders[d],d)
     _gradient_1d!(g,x,orders[d],d)
   end
-  
+
   o = one(T)
   k = 1
 
diff --git a/src/ReferenceFEs/LagrangianDofBases.jl b/src/ReferenceFEs/LagrangianDofBases.jl
index 8573d222c..07b79e94e 100644
--- a/src/ReferenceFEs/LagrangianDofBases.jl
+++ b/src/ReferenceFEs/LagrangianDofBases.jl
@@ -68,13 +68,13 @@ end
 
 # Node major implementation
 function _generate_dof_layout_node_major(::Type{T},nnodes::Integer) where T<:MultiValue
-  ncomps = num_components(T)
   V = change_eltype(T,Int)
+  ncomps = num_indep_components(T)
   ndofs = ncomps*nnodes
   dof_to_comp = zeros(Int,ndofs)
   dof_to_node = zeros(Int,ndofs)
-  node_and_comp_to_dof = zeros(V,nnodes)
-  m = zero(Mutable(V))
+  node_and_comp_to_dof = Vector{V}(undef,nnodes)
+  m = zeros(Int,ncomps)
   for node in 1:nnodes
     for comp in 1:ncomps
       o = nnodes*(comp-1)
@@ -83,7 +83,7 @@ function _generate_dof_layout_node_major(::Type{T},nnodes::Integer) where T<:Mul
       dof_to_node[dof] = node
       m[comp] = dof
     end
-    node_and_comp_to_dof[node] = m
+    node_and_comp_to_dof[node] = V(m...)
   end
   (dof_to_node, dof_to_comp, node_and_comp_to_dof)
 end
@@ -113,8 +113,8 @@ function evaluate!(cache,b::LagrangianDofBasis,field)
   vals = evaluate!(cf,field,b.nodes)
   ndofs = length(b.dof_to_node)
   T = eltype(vals)
-  ncomps = num_components(T)
-  @check ncomps == num_components(eltype(b.node_and_comp_to_dof)) """\n
+  ncomps = num_indep_components(T)
+  @check ncomps == num_indep_components(eltype(b.node_and_comp_to_dof)) """\n
   Unable to evaluate LagrangianDofBasis. The number of components of the
   given Field does not match with the LagrangianDofBasis.
 
@@ -135,8 +135,8 @@ function _evaluate_lagr_dof!(c::AbstractVector,node_comp_to_val,node_and_comp_to
     comp_to_dof = node_and_comp_to_dof[node]
     comp_to_val = node_comp_to_val[node]
     for comp in 1:ncomps
-      dof = comp_to_dof[comp]
-      val = comp_to_val[comp]
+      dof = indep_comp_getindex(comp_to_dof,comp)
+      val = indep_comp_getindex(comp_to_val,comp)
       r[dof] = val
     end
   end
@@ -152,8 +152,8 @@ function _evaluate_lagr_dof!(c::AbstractMatrix,node_pdof_comp_to_val,node_and_co
     for pdof in 1:npdofs
       comp_to_val = node_pdof_comp_to_val[node,pdof]
       for comp in 1:ncomps
-        dof = comp_to_dof[comp]
-        val = comp_to_val[comp]
+        dof = indep_comp_getindex(comp_to_dof,comp)
+        val = indep_comp_getindex(comp_to_val,comp)
         r[dof,pdof] = val
       end
     end
diff --git a/src/ReferenceFEs/LagrangianRefFEs.jl b/src/ReferenceFEs/LagrangianRefFEs.jl
index 2a2cb10b5..af14e2cc8 100644
--- a/src/ReferenceFEs/LagrangianRefFEs.jl
+++ b/src/ReferenceFEs/LagrangianRefFEs.jl
@@ -222,7 +222,7 @@ end
 function _generate_face_own_dofs(face_own_nodes, node_and_comp_to_dof)
   faces = 1:length(face_own_nodes)
   T = eltype(node_and_comp_to_dof)
-  comps = 1:num_components(T)
+  comps = 1:num_indep_components(T)
   face_own_dofs = [Int[] for i in faces]
   for face in faces
     nodes = face_own_nodes[face]
@@ -254,7 +254,7 @@ function  _generate_face_own_dofs_permutations(
   face_own_nodes_permutations, node_and_comp_to_dof, face_own_nodes, face_own_dofs)
 
   T = eltype(node_and_comp_to_dof)
-  ncomps = num_components(T)
+  ncomps = num_indep_components(T)
 
   face_own_dofs_permutations = Vector{Vector{Int}}[]
   for  (face, pindex_to_inode_to_pinode) in enumerate(face_own_nodes_permutations)
diff --git a/src/TensorValues/MultiValueTypes.jl b/src/TensorValues/MultiValueTypes.jl
index c038f5649..c59571246 100644
--- a/src/TensorValues/MultiValueTypes.jl
+++ b/src/TensorValues/MultiValueTypes.jl
@@ -31,7 +31,10 @@ num_components(::Type{<:Number}) = 1
 num_components(::Number) = num_components(Number)
 num_components(T::Type{<:MultiValue}) = @unreachable "$T type is too abstract to count its components, provide a (parametric) concrete type"
 
-"Number of independant components, that is num_component(::Number) minus the number of components determined by symetries or constraints"
+"""
+Number of independant components, that is `num_components(::Type{T})` minus the
+number of components determined from others by symmetries or constraints.
+"""
 num_indep_components(::Type{T}) where T<:Number = num_components(T)
 num_indep_components(::T) where T<:Number = num_indep_components(T)
 
@@ -40,6 +43,31 @@ function n_components(a)
   error(msg)
 end
 
+# This should probably not be exported, as (accessing) the data field of
+# MultiValue is not a public api
 function data_index(::Type{<:MultiValue},i...)
   @abstractmethod
 end
+
+"""
+    indep_comp_getindex(a::Number,i)
+
+Get the ith independent component of `a`. It only differs from `getindex(a,i)`
+when the components of `a` are linked, see [`num_indep_components`](@ref), and
+`i` should be in `1:num_indep_components(a)`.
+"""
+function indep_comp_getindex(a::Number,i)
+  @check 1 <= i <= num_indep_components(Number)
+  a[i]
+end
+
+function indep_comp_getindex(a::T,i) where {T<:MultiValue}
+  @check 1 <= i <= num_indep_components(T)
+  _get_data(a,i)
+end
+
+# abstraction of Multivalue data access in case subtypes of MultiValue don't
+# store its data in a data field
+function _get_data(a::MultiValue,i)
+  a.data[i]
+end
diff --git a/src/TensorValues/TensorValues.jl b/src/TensorValues/TensorValues.jl
index 1ac4ee2d8..7a9fa4bf3 100644
--- a/src/TensorValues/TensorValues.jl
+++ b/src/TensorValues/TensorValues.jl
@@ -63,6 +63,7 @@ export ⋅¹
 export ⋅²
 export double_contraction
 export data_index
+export indep_comp_getindex
 
 import Base: show
 import Base: zero, one
diff --git a/test/ReferenceFEsTests/CLagrangianRefFEsTests.jl b/test/ReferenceFEsTests/CLagrangianRefFEsTests.jl
index 02855bb38..1723b85ab 100644
--- a/test/ReferenceFEsTests/CLagrangianRefFEsTests.jl
+++ b/test/ReferenceFEsTests/CLagrangianRefFEsTests.jl
@@ -35,6 +35,20 @@ dofs = LagrangianDofBasis(VectorValue{3,Float64},TET,1)
 @test dofs.nodes == Point{3,Float64}[(0,0,0), (1,0,0), (0,1,0), (0,0,1)]
 @test dofs.node_and_comp_to_dof == VectorValue{3,Int}[(1,5,9), (2,6,10), (3,7,11), (4,8,12)]
 
+dofs = LagrangianDofBasis(TensorValue{2,2,Float64},TET,1)
+@test dofs.nodes == Point{3,Float64}[(0,0,0), (1,0,0), (0,1,0), (0,0,1)]
+@test dofs.node_and_comp_to_dof == TensorValue{2,2,Int}[(1,5,9,13), (2,6,10,14), (3,7,11,15), (4,8,12,16)]
+
+dofs = LagrangianDofBasis(SymTensorValue{2,Float64},TET,1)
+@test dofs.nodes == Point{3,Float64}[(0,0,0), (1,0,0), (0,1,0), (0,0,1)]
+@test dofs.node_and_comp_to_dof == SymTensorValue{2,Int}[(1,5,9), (2,6,10), (3,7,11), (4,8,12)]
+
+# For SymTracelessTensorValue, the last index in data should not be accessed,
+# as it is minus the sum of the D-1 diagonal value, so is not free/independent
+dofs = LagrangianDofBasis(SymTracelessTensorValue{2,Float64},TET,1)
+@test dofs.nodes == Point{3,Float64}[(0,0,0), (1,0,0), (0,1,0), (0,0,1)]
+@test dofs.node_and_comp_to_dof == SymTracelessTensorValue{2,Int}[(1,5), (2,6), (3,7), (4,8)]
+
 dofs = LagrangianDofBasis(Float64,WEDGE,(2,2,2))
 r = Point{3,Float64}[
   (0.0,0.0,0.0),(1.0,0.0,0.0),(0.0,1.0,0.0),
@@ -48,6 +62,15 @@ r = Point{3,Float64}[
 dofs = LagrangianDofBasis(VectorValue{2,Int},VERTEX,())
 @test dofs.node_and_comp_to_dof == VectorValue{2,Int}[(1,2)]
 
+dofs = LagrangianDofBasis(TensorValue{2,2,Int},VERTEX,())
+@test dofs.node_and_comp_to_dof == TensorValue{2,2,Int}[(1,2,3,4)]
+
+dofs = LagrangianDofBasis(SymTensorValue{2,Int},VERTEX,())
+@test dofs.node_and_comp_to_dof == SymTensorValue{2,Int}[(1,2,3)]
+
+dofs = LagrangianDofBasis(SymTracelessTensorValue{2,Int},VERTEX,())
+@test dofs.node_and_comp_to_dof == SymTracelessTensorValue{2,Int}[(1,2)]
+
 b = MonomialBasis(VectorValue{2,Int},VERTEX,())
 @test length(b) == 2
 @test evaluate(b,Point{0,Int}[(),()]) == VectorValue{2,Int}[(1, 0) (0, 1); (1, 0) (0, 1)]
diff --git a/test/ReferenceFEsTests/LagrangianDofBasesTests.jl b/test/ReferenceFEsTests/LagrangianDofBasesTests.jl
index 82b41d107..55c071467 100644
--- a/test/ReferenceFEsTests/LagrangianDofBasesTests.jl
+++ b/test/ReferenceFEsTests/LagrangianDofBasesTests.jl
@@ -52,4 +52,75 @@ dbb = [
 
 test_dof_array(db,b,dbb)
 
+T = TensorValue{2,2,Float64}
+db = LagrangianDofBasis(T,x)
+@test db.nodes === x
+@test db.node_and_comp_to_dof == TensorValue{2,2,Int}[(1,5,9,13), (2,6,10,14), (3,7,11,15), (4,8,12,16)]
+@test db.dof_to_node == [1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4]
+@test db.dof_to_comp == [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4]
+
+v = TensorValue(1,2,3,4)
+f = GenericField(x->v*x[1])
+dbf = [0, 1, 0, 1, 0, 2, 0, 2, 0, 3, 0, 3, 0, 4, 0, 4]
+
+test_dof_array(db,f,dbf)
+
+ndof = 16
+b = fill(f,ndof)
+bx = evaluate(b,x)
+dbb = [
+  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
+  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2;
+  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3;
+  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4;]
+
+test_dof_array(db,b,dbb)
+
+
+T = SymTensorValue{2,Float64}
+db = LagrangianDofBasis(T,x)
+@test db.nodes === x
+@test db.node_and_comp_to_dof == SymTensorValue{2,Int}[(1,5,9), (2,6,10), (3,7,11), (4,8,12)]
+@test db.dof_to_node == [1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4]
+@test db.dof_to_comp == [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3]
+
+v = SymTensorValue(1,2,3)
+f = GenericField(x->v*x[1])
+dbf = [0, 1, 0, 1, 0, 2, 0, 2, 0, 3, 0, 3]
+
+test_dof_array(db,f,dbf)
+
+ndof = 12
+b = fill(f,ndof)
+bx = evaluate(b,x)
+dbb = [
+  0 0 0 0 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1 1 1 1 1; 0 0 0 0 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1 1 1 1 1;
+  0 0 0 0 0 0 0 0 0 0 0 0; 2 2 2 2 2 2 2 2 2 2 2 2; 0 0 0 0 0 0 0 0 0 0 0 0; 2 2 2 2 2 2 2 2 2 2 2 2;
+  0 0 0 0 0 0 0 0 0 0 0 0; 3 3 3 3 3 3 3 3 3 3 3 3; 0 0 0 0 0 0 0 0 0 0 0 0; 3 3 3 3 3 3 3 3 3 3 3 3;]
+
+test_dof_array(db,b,dbb)
+
+
+T = SymTracelessTensorValue{2,Float64}
+db = LagrangianDofBasis(T,x)
+@test db.nodes === x
+@test db.node_and_comp_to_dof == SymTracelessTensorValue{2,Int}[(1,5), (2,6), (3,7), (4,8)]
+@test db.dof_to_node == [1, 2, 3, 4, 1, 2, 3, 4]
+@test db.dof_to_comp == [1, 1, 1, 1, 2, 2, 2, 2]
+
+v = SymTracelessTensorValue(1,2)
+f = GenericField(x->v*x[1])
+dbf = [0, 1, 0, 1, 0, 2, 0, 2]
+
+test_dof_array(db,f,dbf)
+
+ndof = 8
+b = fill(f,ndof)
+bx = evaluate(b,x)
+dbb = [
+  0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1; 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1;
+  0 0 0 0 0 0 0 0; 2 2 2 2 2 2 2 2; 0 0 0 0 0 0 0 0; 2 2 2 2 2 2 2 2;]
+
+test_dof_array(db,b,dbb)
+
 end # module