diff --git a/docs/src/modules/Polynomials.md b/docs/src/modules/Polynomials.md
index 0070a304a..48f8e22d0 100644
--- a/docs/src/modules/Polynomials.md
+++ b/docs/src/modules/Polynomials.md
@@ -215,7 +215,7 @@ these bases are of maximum degree K+1 and the divergence of their elements are i
 Those bases are a generalization of the scalar Bernstein bases to the spaces for
 the two principal finite element families forming a de Rham complex on simplices.
 They are respectively implemented by [`BarycentricPmΛBasis`](@ref) and
-[`BarycentricPmΛBasis`](@ref). Their definition with references, and implementation
+[`BarycentricPΛBasis`](@ref). Their definition with references, and implementation
 details are provided in [this](@ref "Bernstein-basis-generalization-for-``\mathcal{P}Λ``-spaces")
 developer note.
 
diff --git a/src/Adaptivity/MacroFEs.jl b/src/Adaptivity/MacroFEs.jl
index d0522263e..1c5f2c340 100644
--- a/src/Adaptivity/MacroFEs.jl
+++ b/src/Adaptivity/MacroFEs.jl
@@ -404,12 +404,12 @@ function ReferenceFEs.get_face_own_dofs_permutations(reffe::GenericRefFE{<:Macro
 end
 
 function ReferenceFEs.get_face_own_dofs(reffe::GenericRefFE{<:MacroRefFEName}, ::L2Conformity)
-  return ReferenceFEs._get_face_own_dofs_l2(reffe)
+  return ReferenceFEs.l2_face_own_dofs(reffe)
 end
 
 function ReferenceFEs.get_face_own_dofs_permutations(reffe::GenericRefFE{<:MacroRefFEName}, ::L2Conformity)
-  face_own_dofs = ReferenceFEs.get_face_own_dofs(reffe,L2Conformity())
-  return ReferenceFEs._trivial_face_own_dofs_permutations(face_own_dofs)
+  face_own_dofs = ReferenceFEs.get_face_own_dofs(reffe, L2Conformity())
+  return ReferenceFEs.l2_face_own_dofs_permutations(face_own_dofs)
 end
 
 """
diff --git a/src/FESpaces/ConformingFESpaces.jl b/src/FESpaces/ConformingFESpaces.jl
index 34f1bc0c9..0e5f0edf4 100644
--- a/src/FESpaces/ConformingFESpaces.jl
+++ b/src/FESpaces/ConformingFESpaces.jl
@@ -37,7 +37,7 @@ function GenericCellConformity(
   cell_lface_own_ldofs :: AbstractVector{<:AbstractVector{<:AbstractVector{<:Integer}}},
   cell_d_num_dfaces :: AbstractVector{<:AbstractVector{<:Integer}},
 )
-  cell_lface_pindex_pdofs = lazy_map(ReferenceFEs._trivial_face_own_dofs_permutations, cell_lface_own_ldofs)
+  cell_lface_pindex_pdofs = lazy_map(ReferenceFEs.l2_face_own_dofs_permutations, cell_lface_own_ldofs)
   cell_num_dofs = lazy_map(x -> sum(length, x), cell_lface_own_ldofs)
   cell_ldof_comp = lazy_map(n -> fill(0,n), cell_num_dofs)
   GenericCellConformity(
@@ -210,10 +210,10 @@ CellConformity(cell_fe::CellFE) = cell_fe.cell_conformity
 Geometry.num_cells(cell_fe::CellFE) = num_cells(CellConformity(cell_fe))
 Geometry.get_cell_type(cell_fe::CellFE) = get_cell_type(CellConformity(cell_fe))
 
-# This constructor allows for provided shapefuns and dofs, which 
+# This constructor allows for provided shapefuns and dofs, which
 # is necessary for GridapDistributed
 function CellFE(
-  model::DiscreteModel, cell_reffe::AbstractArray{<:ReferenceFE}, 
+  model::DiscreteModel, cell_reffe::AbstractArray{<:ReferenceFE},
   cell_shapefuns::AbstractArray, cell_dof_basis::AbstractArray,
   conformity::Conformity
 )
diff --git a/src/FESpaces/Pullbacks.jl b/src/FESpaces/Pullbacks.jl
index 357129804..f9b4f0ce9 100644
--- a/src/FESpaces/Pullbacks.jl
+++ b/src/FESpaces/Pullbacks.jl
@@ -91,14 +91,14 @@ function compute_cell_bases_changes(
 )
   D = num_cell_dims(model)
   poly = only(get_polytopes(model))
-  if (D==2) || is_simplex(poly)
-    # For these cases, we do not need to aply a sign flip
-    return nothing
-  elseif (D==3) && is_n_cube(poly)
-    change = get_sign_flip(model, cell_reffe)
-    return (change,change)
+  change = get_edge_flip(model, cell_reffe)
+  if (D==3) && is_n_cube(poly)
+    facet_change = get_sign_flip(model, cell_reffe)
+    change = lazy_map(*, facet_change, change)
+  elseif !((D==2) || is_simplex(poly))
+    @notimplemented
   end
-  @notimplemented
+  return (change,change)
 end
 
 using Gridap.ReferenceFEs: DoubleContraVariantPiolaMap
@@ -207,6 +207,72 @@ function compute_facet_owners(model::DiscreteModel{Dc}, select_nbor=maximum) whe
   return owners
 end
 
+###############
+# EdgeSignMap #
+###############
+
+function get_edge_flip(model::DiscreteModel, cell_reffe, sign_map = EdgeSignMap(model))
+  get_edge_own_dofs(reffe) = view(get_face_own_dofs(reffe),get_dimrange(get_polytope(reffe),1))
+  cell_edge_own_dofs = lazy_map(get_edge_own_dofs, cell_reffe)
+  cell_ids = IdentityVector(Int32(num_cells(model)))
+  return lazy_map(sign_map, cell_reffe, cell_edge_own_dofs, cell_ids)
+end
+
+"""
+    struct EdgeSignMap <: Map
+      ...
+    end
+
+The `EdgeSignMap` computes the signs to apply to the mapped reference edge
+tangents, for each edge of a physical cell.
+
+`get_edge_tangent` returns a fixed reference-polytope tangent direction
+(`v2-v1` in the cell's local edge-vertex order) for each local edge. When a
+cell's local edge-vertex order is reversed relative to the canonical order
+used to build `face_own_dofs_permutations` (`pindex == 2`), the physical
+tangent direction is reversed too. This contributes a uniform `-1` factor on
+all dofs of that edge, on top of the index permutation already accounted for
+by `face_own_dofs_permutations`.
+"""
+struct EdgeSignMap{T} <: Map
+  model::T
+end
+
+function return_value(k::EdgeSignMap,reffe,edge_own_dofs,cell)
+  Diagonal(fill(one(Float64), num_dofs(reffe)))
+end
+
+function return_cache(k::EdgeSignMap,reffe,edge_own_dofs,cell)
+  model = k.model
+  topo = get_grid_topology(model)
+
+  cell_edge_pindex = get_cell_permutations(topo, 1)
+  cell_edge_pindex_cache = array_cache(cell_edge_pindex)
+
+  return cell_edge_pindex, cell_edge_pindex_cache, CachedVector(Float64)
+end
+
+function evaluate!(cache,k::EdgeSignMap,reffe,edge_own_dofs,cell)
+  cell_edge_pindex,cell_edge_pindex_cache,dof_sign_cache = cache
+
+  setsize!(dof_sign_cache, (num_dofs(reffe),))
+  dof_sign = dof_sign_cache.array
+
+  o = one(eltype(dof_sign))
+  fill!(dof_sign, o)
+
+  pindices = getindex!(cell_edge_pindex_cache,cell_edge_pindex,cell)
+  for (ledge,pindex) in enumerate(pindices)
+    if pindex == 2
+      for dof in edge_own_dofs[ledge]
+        dof_sign[dof] = -o
+      end
+    end
+  end
+
+  return Diagonal(dof_sign)
+end
+
 #################
 # DOFScalingMap #
 #################
diff --git "a/src/Polynomials/BarycentricP\316\233Bases.jl" "b/src/Polynomials/BarycentricP\316\233Bases.jl"
index 9151cfe5e..a8f7e4ac4 100644
--- "a/src/Polynomials/BarycentricP\316\233Bases.jl"
+++ "b/src/Polynomials/BarycentricP\316\233Bases.jl"
@@ -571,6 +571,15 @@ function PmΛ_bubbles(r,k,D)
   bubbles
 end
 
+"""
+    _compute_PmΛ_basis_coefficients!(m,::Val{k},D,b,vertices,indices) where k
+
+The forms
+
+    dλᵢ₁ ⋀ ...  ⋀ dλᵢₖ
+
+for (i1, i2, ..., ik) in `_sorted_combinations(D+1,k)`, stored in place in `m`.
+"""
 function _compute_PmΛ_basis_coefficients!(m,::Val{k},D,b,vertices,indices) where k
   V = eltype(m)
   M = transpose(b.cart_to_bary_matrix[:,2:end])
@@ -781,7 +790,7 @@ function _compute_PΛ_basis_form_coefficient!(Ψ,r,k,::Val{D},b,vertices,indices
   Vk = Val(k)
   V = eltype(Ψ)
 
-  iszero(r) && return _order_0_Ψ!(Ψ)
+  iszero(r) && return _order_0_Ψ!(Ψ,b,k,Val(D),vertices,indices)
 
   T = eltype(V)
   α_prec = ntuple(_->-1, N)
@@ -812,24 +821,48 @@ end
   end
 end
 
-function _order_0_Ψ!(Ψ)
+"""
+    _order_0_Ψ!(Ψ,b,k,::Val(D))
+
+We have to chose a basis of the constant `D`-dimensional `k`-forms. There are
+`binomial(D,D-k)` such forms (the number of components).
+
+The natural choice is to choose `D` forms ``dλ_2``, ... ,``dλ_(D+1)`` and
+build the cannonical `k`-form basis generated by these forms:
+
+    dλᵢ₁ ⋀ ...  ⋀ dλᵢₖ
+
+for all increasing indices. We chose the ``dλᵢ`` forms starting at ``i=2``
+because their coordinates are δ_(i,i-1) in the reference simplex, and
+
+    dλ₁  = - Σ_(i≥2) dλᵢ
+"""
+function _order_0_Ψ!(Ψ,b,k,::Val{D},vertices,indices) where D
+  N = D+1
   V = eltype(Ψ)
-  T = eltype(V)
-  for w in eachindex(Ψ)
-    Ψ[w] = ntuple( i -> T(i==w), length(V))
+  m = zero(MVector{binomial(N,k),V})
+  _compute_PmΛ_basis_coefficients!(m,Val(k),D,b,vertices,indices)
+
+  w = 1
+  @inbounds for (i, I) in enumerate(_sorted_combinations(N,k))
+    if 1 ∉ I
+      Ψ[w] = m[i]
+      w += 1
+    end
   end
+
   nothing
 end
 
 function _compute_PΛ_basis_form_coefficient!(
-  Ψ,r,k,::Val,b,vertices::Nothing,indices)
+  Ψ,r,k,::Val{D},b,vertices::Nothing,indices) where D
 
   Vk = Val(k)
   V = eltype(Ψ)
   T = eltype(V)
   Ψw = Mutable(V)(undef)
 
-  iszero(r) && return _order_0_Ψ!(Ψ)
+  iszero(r) && return _order_0_Ψ!(Ψ,b,k,Val(D),vertices,indices)
 
   @inbounds for (F, bubble_functions) in indices.bubbles
     for (w, α, _, J) in bubble_functions
diff --git a/src/Polynomials/ExteriorCalculusBases.jl b/src/Polynomials/ExteriorCalculusBases.jl
index fd13ea4fe..234bbc5db 100644
--- a/src/Polynomials/ExteriorCalculusBases.jl
+++ b/src/Polynomials/ExteriorCalculusBases.jl
@@ -38,9 +38,9 @@ function FEEC_poly_basis(::Val{D},::Type{T},r,k,F::Symbol,PT=_default_poly_type(
   @assert PT <: Polynomial
 
   # these call FEEC_space_definition_checks internally
-  if !cart_prod
-    F == :P⁻ && PT == Bernstein && return BarycentricPmΛBasis(Val(D),T,r,k,vertices; rotate_90, DG_calc)
-    F == :P  && PT == Bernstein && return BarycentricPΛBasis( Val(D),T,r,k,vertices; rotate_90, DG_calc)
+  if !cart_prod && k ∉ (0,D) && PT == Bernstein
+    F == :P⁻ && return BarycentricPmΛBasis(Val(D),T,r,k,vertices; rotate_90, DG_calc)
+    F == :P  && return BarycentricPΛBasis( Val(D),T,r,k,vertices; rotate_90, DG_calc)
   end
 
 
@@ -51,13 +51,13 @@ function FEEC_poly_basis(::Val{D},::Type{T},r,k,F::Symbol,PT=_default_poly_type(
     # Scalar H1 conforming functions
     @notimplementedif r < 0
     if     F == :P⁻ || F == :P # Lagrange, 𝓟r space
-      PT == Bernstein && return BernsteinBasisOnSimplex{D}(T,r,vertices) # only if cart_prod
+      PT == Bernstein && return BernsteinBasisOnSimplex{D}(T,r,vertices)
       _ensure_hierarchical(PT)
       CartProdPolyBasis(PT,Val(D),T,r,_p_filter)
     elseif F == :Q⁻            # Lagrange, 𝓠r space
       CartProdPolyBasis(PT,Val(D),T,r,_q_filter)
     elseif F == :S             # Lagrange, 𝓢r space
-      PT==ModalC0 || _ensure_hierarchical(PT)
+      PT==ModalC0 || _ensure_hierarchical(PT) # because ModalC0 is hierarchical from order 1, which we need
       CartProdPolyBasis(PT,Val(D),T,r,_ser_filter)
     end
 
@@ -65,11 +65,11 @@ function FEEC_poly_basis(::Val{D},::Type{T},r,k,F::Symbol,PT=_default_poly_type(
   elseif k == D
     # Scalar L2 conforming densities
     if     F == :P⁻ # Lagrange, 𝓟r₋1 space
-      PT == Bernstein && return BernsteinBasisOnSimplex{D}(T,r-1,vertices) # only if cart_prod
+      PT == Bernstein && return BernsteinBasisOnSimplex{D}(T,r-1,vertices)
       _ensure_hierarchical(PT)
       CartProdPolyBasis(PT,Val(D),T,r-1,_p_filter)
     elseif F == :P  # Lagrange, 𝓟r space
-      PT == Bernstein && return BernsteinBasisOnSimplex{D}(T,r,vertices) # only if cart_prod
+      PT == Bernstein && return BernsteinBasisOnSimplex{D}(T,r,vertices)
       _ensure_hierarchical(PT)
       CartProdPolyBasis(PT,Val(D),T,r,_p_filter)
     elseif F == :Q⁻ # Lagrange, 𝓠r₋1 space
diff --git a/src/ReferenceFEs/BDMRefFEs.jl b/src/ReferenceFEs/BDMRefFEs.jl
index e046bfd3c..92f02f9e8 100644
--- a/src/ReferenceFEs/BDMRefFEs.jl
+++ b/src/ReferenceFEs/BDMRefFEs.jl
@@ -1,7 +1,7 @@
 """
     struct BDM <: ReferenceFEName
 """
-struct BDM <: ReferenceFEName end
+struct BDM <: MomentBasedRefFEName end
 
 """
     const bdm = BDM()
diff --git a/src/ReferenceFEs/CrouzeixRaviartRefFEs.jl b/src/ReferenceFEs/CrouzeixRaviartRefFEs.jl
index 5f2ab7b88..940603d0c 100644
--- a/src/ReferenceFEs/CrouzeixRaviartRefFEs.jl
+++ b/src/ReferenceFEs/CrouzeixRaviartRefFEs.jl
@@ -1,7 +1,7 @@
 """
     struct CrouzeixRaviart  <: ReferenceFEName
 """
-struct CrouzeixRaviart <: ReferenceFEName end
+struct CrouzeixRaviart <: MomentBasedRefFEName end
 
 """
     const crouzeix_raviart = CrouzeixRaviart()
diff --git a/src/ReferenceFEs/ExtrusionPolytopes.jl b/src/ReferenceFEs/ExtrusionPolytopes.jl
index b060ed5a0..f1f4f196b 100644
--- a/src/ReferenceFEs/ExtrusionPolytopes.jl
+++ b/src/ReferenceFEs/ExtrusionPolytopes.jl
@@ -22,6 +22,27 @@ struct DFace{D} <: GridapType
   nf_nfs::Vector{Vector{Int}}
   nf_dimranges::Vector{Vector{UnitRange{Int}}}
   nf_dims::Vector{Vector{Int}}
+
+  function DFace{D}(
+    extrusion::Point{D,Int},
+    nfaces::Vector{NFace{D}},
+    dimranges::Vector{UnitRange{Int}},
+    dims::Vector{Int},
+    nf_nfs::Vector{Vector{Int}},
+    nf_dimranges::Vector{Vector{UnitRange{Int}}},
+    nf_dims::Vector{Vector{Int}}
+  ) where D
+
+    extrusion = if isone(D)
+      Point{D}(HEX_AXIS)
+    elseif D>1 && extrusion[1] != extrusion[2]
+      Point{D}(extrusion[2], extrusion[2:end]...)
+    else
+      extrusion
+    end
+
+    new{D}(extrusion, nfaces, dimranges, dims, nf_nfs, nf_dimranges, nf_dims)
+  end
 end
 
 """
@@ -30,9 +51,15 @@ end
       # + private fields
     end
 
-Concrete type for polytopes that can be represented with an "extrusion" tuple.
-The underlying extrusion is available in the field `extrusion`. Instances of
-this type can be obtained with the constructors
+Concrete type for polytopes that can be represented with an "extrusion" tuple,
+containing [`HEX_AXIS`](@ref) and [`TET_AXIS`](@ref) constants.  See e.g.
+[`QUAD`](@ref), [`TRI`](@ref), [`PYRAMID`](@ref).
+
+The underlying extrusion is available via [`get_extrusion`](@ref). Since the
+first extrusion has no meaning, it overwritten to `HEX_AXIS` for 1D `ExtrusionPolytope`
+or to the second one for D≥2.
+
+Instances of this type can be obtained with the constructors
 
 - [`Polytope(extrusion::Int...)`](@ref)
 - [`ExtrusionPolytope(extrusion::Int...)`](@ref)
@@ -50,10 +77,11 @@ struct ExtrusionPolytope{D} <: Polytope{D}
   vertex_perms::Vector{Vector{Int}}
   face_vertex_perms::Vector{Vector{Vector{Int}}}
   m_n_to_mface_to_nface::Matrix{Vector{Vector{Int}}}
+
   function ExtrusionPolytope(dface::DFace{D}) where D
     vertex_coords = _vertices_coordinates(Float64, dface)
     face_normals, face_orientations = _face_normals(Float64, dface)
-    vertex_perms = _precompute_vertex_perms_if_possible(dface)
+    vertex_perms = _compute_vertex_perms(dface)
     face_vertex_perms = _admissible_face_vertex_permutations(dface)
     m_n_to_mface_to_nface = _precompute_m_n_to_mface_to_nface(dface)
     new{D}(
@@ -162,18 +190,7 @@ function Polytope{D}(p::ExtrusionPolytope{D}, Dfaceid::Integer) where D
 end
 
 function (==)(a::ExtrusionPolytope{D}, b::ExtrusionPolytope{D}) where D
-  #The first axis is irrelevant here
-  ea = Point(Tuple(a.extrusion)[2:end])
-  eb = Point(Tuple(b.extrusion)[2:end])
-  ea == eb
-end
-
-function (==)(a::ExtrusionPolytope{1}, b::ExtrusionPolytope{1})
-  true
-end
-
-function (==)(a::ExtrusionPolytope{0}, b::ExtrusionPolytope{0})
-  true
+  a.extrusion == b.extrusion
 end
 
 function get_vertex_coordinates(p::ExtrusionPolytope)
@@ -206,28 +223,12 @@ function get_face_vertex_permutations(p::ExtrusionPolytope)
   p.face_vertex_perms
 end
 
-function is_simplex(p::ExtrusionPolytope)
-  all(Tuple(p.extrusion) .== TET_AXIS)
+function is_simplex(p::ExtrusionPolytope{D}) where D
+   D<2 || all(Tuple(p.extrusion) .== TET_AXIS)
 end
 
-function is_n_cube(p::ExtrusionPolytope)
-  all(Tuple(p.extrusion) .== HEX_AXIS)
-end
-
-function is_simplex(p::ExtrusionPolytope{0})
-  true
-end
-
-function is_n_cube(p::ExtrusionPolytope{0})
-  true
-end
-
-function is_simplex(p::ExtrusionPolytope{1})
-  true
-end
-
-function is_n_cube(p::ExtrusionPolytope{1})
-  true
+function is_n_cube(p::ExtrusionPolytope{D}) where D
+  D<2 || all(Tuple(p.extrusion) .== HEX_AXIS)
 end
 
 function simplexify(p::ExtrusionPolytope{0}; kwargs...)
@@ -777,7 +778,7 @@ function _admissible_face_vertex_permutations_fill!(perms, p::DFace{D}, ::Val{d}
   faceids = p.dimranges[d+1]
   for faceid in faceids
     f = DFace{d}(p, faceid)
-    f_perms = _precompute_vertex_perms_if_possible(f)
+    f_perms = _compute_vertex_perms(f)
     push!(perms, f_perms)
   end
   _admissible_face_vertex_permutations_fill!(perms, p, Val{d + 1}())
@@ -785,113 +786,71 @@ function _admissible_face_vertex_permutations_fill!(perms, p::DFace{D}, ::Val{d}
 end
 
 function _admissible_face_vertex_permutations_fill!(perms, p::DFace{D}, ::Val{D}) where D
-  f_perms = _precompute_vertex_perms_if_possible(p)
+  f_perms = _compute_vertex_perms(p)
   push!(perms, f_perms)
   nothing
 end
 
-function _precompute_vertex_perms_if_possible(p::DFace{D}) where D
-  if D in (0, 1)
-    perms = _admissible_permutations_simplex(p)
-  elseif D == 2
-    perms = _admissible_permutations(p)
-  else
-    perms = _admissible_permutations_identity(p)
-  end
-  perms
-end
-
 # It generates all the admissible permutations of nodes that lead to an
 # admissible polytope
-function _admissible_permutations(p::DFace{D}) where D
-  if D > 3
-    @warn "Computing permutations for a polytope of dim > 3 is overkill"
-  end
-  if D in (0, 1) || all(map(i -> i == TET_AXIS, Tuple(p.extrusion)[2:end]))
-    perms = _admissible_permutations_simplex(p)
-  elseif all(map(i -> i == HEX_AXIS, Tuple(p.extrusion)[2:end]))
-    perms = _admissible_permutations_n_cube(p)
+function _compute_vertex_perms(p::DFace{D}) where D
+  perms = if D<2 || all(==(TET_AXIS), p.extrusion)
+    _admissible_permutations_simplex(p)
+  elseif all(==(HEX_AXIS), p.extrusion)
+    _admissible_permutations_n_cube(p)
+
+  elseif Tuple(p.extrusion) == (HEX_AXIS, HEX_AXIS, TET_AXIS) # PYRAMID
+    [[1, 2, 3, 4, 5], [1, 3, 2, 4, 5], [2, 1, 4, 3, 5], [2, 4, 1, 3, 5],
+     [3, 1, 4, 2, 5], [3, 4, 1, 2, 5], [4, 2, 3, 1, 5], [4, 3, 2, 1, 5]]
+
+  elseif Tuple(p.extrusion) == (TET_AXIS, TET_AXIS, HEX_AXIS) # WEDGE
+    [[1, 2, 3, 4, 5, 6], [4, 5, 6, 1, 2, 3],
+     [1, 3, 2, 4, 6, 5], [4, 6, 5, 1, 3, 2],
+     [2, 1, 3, 5, 4, 6], [5, 4, 6, 2, 1, 3],
+     [2, 3, 1, 5, 6, 4], [5, 6, 4, 2, 3, 1],
+     [3, 1, 2, 6, 4, 5], [6, 4, 5, 3, 1, 2],
+     [3, 2, 1, 6, 5, 4], [6, 5, 4, 3, 2, 1]]
+
   else
-    @notimplemented "admissible vertex permutations only implemented for simplices and n-cubes"
+    @warn "Vertex permutations for polytpe with extrusion $(p.extrusion) is not implemented, returning empty vector."
+    Vector{Vector{Int}}()
   end
-  perms
-end
 
-function _admissible_permutations_identity(p::DFace{D}) where D
-  vs = first(_dimfrom_fs_dimto_fs(p, D, 0))
-  num_vs = length(vs)
-  l = [i for i = 1:num_vs]
-  perms = [l,]
   perms
 end
 
 function _admissible_permutations_simplex(p::DFace{D}) where D
   vs = first(_dimfrom_fs_dimto_fs(p, D, 0))
   num_vs = length(vs)
-  l = [i for i = 1:num_vs]
-  perms = Combinatorics.permutations(l)
+  l = collect(1:num_vs)
+  perms = permutations(l)
   collect(perms)
 end
 
 function _admissible_permutations_n_cube(p::DFace{D}) where D
-  vs = first(_dimfrom_fs_dimto_fs(p, D, 0))
-  num_vs = length(vs)
-  l = [i for i = 1:num_vs]
-  perms = Combinatorics.permutations(l)
-  vertices = p.nfaces[first(p.dimranges)]
-  permuted_vertices = similar(vertices)
-  admissible_perms = Vector{Int}[]
-  grads = _setup_aux_grads(vertices)
-  vol = -1
-  for perm in perms
-    for (j, cj) in enumerate(perm)
-      permuted_vertices[j] = vertices[cj]
-    end
-    jac = _setup_aux_jacobian(grads, permuted_vertices)
-    vol_i = convert(Int, abs(det(jac)))
-    if vol < 0
-      vol = vol_i
+  isometry_perms = Vector{Vector{Int}}(undef,  factorial(D)*2^D )
+
+  n_vertices = 2^D
+  size_li = tfill(2, Val(D))
+
+  li_ref = collect(LinearIndices(size_li))
+  li = copy(li_ref)
+  noflip, flip = 1:2, 2:-1:1
+
+  p_id = 1
+  # Iterate through all 2^D possible compination of plane symmetries
+  for flips in Iterators.product(tfill((noflip,flip), Val(D))...)
+    # Iterate through all D! rotations of the frame (x, y, z,...)
+    @inbounds for axes_perm in permutations(1:D)
+      permutedims!(li, li_ref, axes_perm)
+      perm = copy(li[flips...])
+      perm = reshape(perm, (n_vertices,))
+      isometry_perms[p_id] = perm
+      p_id += 1
     end
-    if vol_i == vol
-      push!(admissible_perms, perm)
-    end
-  end
-  admissible_perms
-end
-
-function _setup_aux_grads(vertices::Vector{NFace{D}}) where D
-  grads = zeros(Point{D,Int}, length(vertices))
-  m = zero(Mutable(Point{D,Int}))
-  for (i, vertex) in enumerate(vertices)
-    x = vertex.anchor
-    for di in 1:D
-      v = 1
-      for k in 1:D
-        if k == di && x[k] == 0
-          v *= -1
-        end
-      end
-      m[di] = v
-    end
-    grads[i] = m
   end
-  grads
-end
 
-function _setup_aux_jacobian(grads, permuted_vertices::Vector{NFace{D}}) where D
-  p0 = zero(Point{D,Int})
-  m = zero(Mutable(outer(p0, p0)))
-  for (i, pvertex) in enumerate(permuted_vertices)
-    x = pvertex.anchor
-    g = grads[i]
-    for di in 1:D
-      v = g[di]
-      for dj in 1:D
-        m[di, dj] += x[dj] * v
-      end
-    end
-  end
-  TensorValue(m)
+  isometry_perms
 end
 
 # Some particular cases
diff --git a/src/ReferenceFEs/GeometricDecompositions.jl b/src/ReferenceFEs/GeometricDecompositions.jl
index b6041a2fa..30e10a30a 100644
--- a/src/ReferenceFEs/GeometricDecompositions.jl
+++ b/src/ReferenceFEs/GeometricDecompositions.jl
@@ -37,6 +37,17 @@ function _l2_conforming_own_funs(shapefuns,p)
   r
 end
 
+function get_face_own_funs_permutations(shapefuns, p::Polytope, conf::Conformity)
+  # TODO This WIP method is not correctly implemented `get_face_own_funs_permutations(b,p, $conf)`
+  face_own_fun = get_face_own_funs(shapefuns, p, conf)
+  l2_face_own_dofs_permutations(face_own_fun)
+end
+
+function get_face_own_funs_permutations(shapefuns, p::Polytope, ::L2Conformity)
+  l2_face_own_fun = _l2_conforming_own_funs(shapefuns,p)
+  l2_face_own_dofs_permutations(l2_face_own_fun )
+end
+
 """
     apply_face_signflip(shapefun, p::Polytope, conf::Conformity)
 
diff --git a/src/ReferenceFEs/LagrangianRefFEs.jl b/src/ReferenceFEs/LagrangianRefFEs.jl
index 39a8eea51..4514fba42 100644
--- a/src/ReferenceFEs/LagrangianRefFEs.jl
+++ b/src/ReferenceFEs/LagrangianRefFEs.jl
@@ -64,19 +64,21 @@ end
 """
     get_face_own_nodes_permutations(reffe::LagrangianRefFE[, conf::Conformity][, d::Integer])
 
-Like [`get_face_own_nodes_permutations`](@ref), but the indices are that of the
+Like [`get_face_own_dofs_permutations`](@ref), but the indices are that of the
 nodes instead of that of the DoFs. They are different for vector/tensor-valued
 elements, for which several DoFs are placed at the same node (one per
-inpedendent component).
+independent component).
 """
-function get_face_own_nodes_permutations(reffe::LagrangianRefFE,conf::Conformity)
-  face_own_nodes = get_face_own_nodes(reffe,conf)
-  _trivial_face_own_dofs_permutations(face_own_nodes)
+function get_face_own_nodes_permutations(reffe::LagrangianRefFE, conf::Conformity)
+  # TODO this is invalid for CDConformity, this method should in theory be
+  # restricted to L2Conformity, and another implemented for CDConformity.
+  face_own_nodes = get_face_own_nodes(reffe, conf)
+  l2_face_own_dofs_permutations(face_own_nodes)
 end
 
 function get_face_own_nodes_permutations(reffe::LagrangianRefFE)
   conf = Conformity(reffe)
-  get_face_own_nodes_permutations(reffe,conf)
+  get_face_own_nodes_permutations(reffe, conf)
 end
 
 """
@@ -320,9 +322,9 @@ function  _generate_face_own_dofs_permutations(
   face_own_dofs_permutations
 end
 
-function get_face_own_dofs_permutations(reffe::LagrangianRefFE,conf::L2Conformity)
+function get_face_own_dofs_permutations(reffe::LagrangianRefFE, conf::L2Conformity)
   face_own_dofs = get_face_own_dofs(reffe,conf)
-  _trivial_face_own_dofs_permutations(face_own_dofs)
+  l2_face_own_dofs_permutations(face_own_dofs)
 end
 
 """
diff --git a/src/ReferenceFEs/ModalScalarRefFEs.jl b/src/ReferenceFEs/ModalScalarRefFEs.jl
index 1b4fdd1f0..f2fc4d5e5 100644
--- a/src/ReferenceFEs/ModalScalarRefFEs.jl
+++ b/src/ReferenceFEs/ModalScalarRefFEs.jl
@@ -5,7 +5,7 @@
 
 where `Name` is either `Lagrangian` or `Serendipity`.
 """
-struct ModalScalar{F} <: ReferenceFEName
+struct ModalScalar{F} <: MomentBasedRefFEName
   ModalScalar{Lagrangian}() = new{Lagrangian}()
   ModalScalar(::Lagrangian) = new{Lagrangian}()
   ModalScalar{Serendipity}() = new{Serendipity}()
@@ -21,6 +21,12 @@ Singletons of the [`ModalScalar`](@ref) reference FE name.
 const modal_lagrangian = ModalScalar(lagrangian)
 const modal_serendipity  = ModalScalar(serendipity)
 
+function _modal_scalar_default_PT(p,F)
+  is_n_cube(p) && F==:S && return Polynomials.ModalC0
+  (is_n_cube(p) || is_simplex(p)) && return Bernstein
+  Monomial
+end
+
 """
     ModalScalarRefFE(::Type{T}, p::Polytope{D}, order::Integer; F::Symbol, kwargs...)
 
@@ -31,13 +37,12 @@ This is a variant of `lagrangian`/`serendipity` elements that is more accurate f
 
 The `kwargs` are [`change_dof`](@ref "`change_dof` keyword argument"),
 [`poly_type`](@ref "`poly_type` keyword argument") and
-[`mom_poly_type`](@ref "`mom_poly_type` keyword argument").
+[`mom_poly_type`](@ref "`mom_poly_type` keyword argument"). `poly_type` defaults to `ModalC0` for `F=:S`.
 
-For `F=:S`, `mom_poly_type` is changed for `Legendre` when `ModalC0` (the
-default) is given because the moment basis need be hierarchical.
+For `F=:S`, `mom_poly_type` is changed for `Legendre` if the given one is not hierarchical.
 """
 function ModalScalarRefFE(::Type{T}, p::Polytope{D}, r::Integer; F::Symbol,
-  change_dof=true, poly_type=_mom_reffe_default_PT(p), mom_poly_type=poly_type) where {T,D}
+  change_dof=true, poly_type=_modal_scalar_default_PT(p,F), mom_poly_type=poly_type) where {T,D}
 
   PT, MPT = poly_type, mom_poly_type
   cart_prod = T <: MultiValue
@@ -65,7 +70,7 @@ function ModalScalarRefFE(::Type{T}, p::Polytope{D}, r::Integer; F::Symbol,
                  : nothing) for d in 0:D ]                # Q⁻ᵨΛᵈ(□ᵈ), ρ = r-1
     elseif F==:S
       prebasis = FEEC_poly_basis(Val(D),T,r,   0,:S,PT; cart_prod)   # SᵣΛ⁰(□ᴰ)
-      MPT = MPT == Polynomials.ModalC0 ? Legendre : MPT
+      MPT = isHierarchical(MPT) ? MPT : Legendre
       mb = [ (r-2d >= 0 ?
                  FEEC_poly_basis(Val(d),T,r-2d,d,:P,MPT; cart_prod)
                  : nothing) for d in 0:D ]                # PᵨΛᵈ(□ᵈ), ρ = r-2*d
diff --git a/src/ReferenceFEs/MomentBasedReferenceFEs.jl b/src/ReferenceFEs/MomentBasedReferenceFEs.jl
index f20e1ea69..664a95db0 100644
--- a/src/ReferenceFEs/MomentBasedReferenceFEs.jl
+++ b/src/ReferenceFEs/MomentBasedReferenceFEs.jl
@@ -1,4 +1,130 @@
 
+"""
+    abstract type MomentBasedRefFEName <: ReferenceFEName end
+
+Supertype for all reference FE names whose DoF basis is a [`MomentBasedDofBasis`](@ref).
+"""
+abstract type MomentBasedRefFEName <: ReferenceFEName end
+
+# get_basis_permutations
+
+"""
+    get_basis_permutations(poly::Polytope, basis) -> Vector{Vector{Int}}
+
+For each admissible vertex permutation `σ` of `poly` — in the same order as
+`get_vertex_permutations(poly)` — return the induced permutation of basis
+function indices, so that `basis[π[i]] ∘ Mσ = basis[i]` where `Mσ` is the
+affine self-map defined by `σ`.
+"""
+function get_basis_permutations(poly::Polytope, basis)
+  # TODO: This is wrong but will not throw an error, which is bad.
+  # We should decide what to do with these.
+  _trivial_basis_permutation(poly, basis)
+end
+
+function _trivial_basis_permutation(poly::Polytope, basis)
+  nfpids = length(get_vertex_permutations(poly))
+  ndofs  = length(basis)
+  [collect(1:ndofs) for _ in 1:nfpids]
+end
+
+function get_basis_permutations(
+  poly::Polytope{D}, b::BernsteinBasisOnSimplex{D}
+) where D
+  @assert is_simplex(poly) && num_dims(poly) == D
+  K     = get_order(b)
+  verts = get_vertex_coordinates(poly)
+  nodes = map(bernstein_terms(K, D)) do α
+    sum(α[i]/K * verts[i] for i in eachindex(α))
+  end
+  compute_own_nodes_permutations(poly, nodes)
+end
+
+function get_basis_permutations(
+  poly::Polytope{D}, b::CartProdPolyBasis{D,V,Bernstein}
+) where {D,V}
+  @assert is_n_cube(poly) && num_dims(poly) == D
+  K     = get_order(b)
+  ncomp = num_indep_components(V)
+  nodes = map(ci -> Point(ntuple(d -> (ci[d]-1)/K, D)), b.terms)
+  scalar_perms = compute_own_nodes_permutations(poly, nodes)
+  isone(ncomp) && return scalar_perms
+  n_terms = length(b.terms)
+  map(scalar_perms) do term_perm
+    dof_perm = Vector{Int}(undef, n_terms * ncomp)
+    for (i, j) in enumerate(term_perm)
+      for c in 1:ncomp
+        dof_perm[(i-1)*ncomp + c] = (j-1)*ncomp + c
+      end
+    end
+    dof_perm
+  end
+end
+
+function get_basis_permutations(
+  poly::Polytope{D}, b::BarycentricPΛBasis{D}
+) where D
+  @assert is_simplex(poly) && num_dims(poly) == D
+  b.k == 0 || return _trivial_basis_permutation(poly, b)
+  K     = get_order(b)
+  verts = get_vertex_coordinates(poly)
+  nodes = Vector{eltype(verts)}(undef, length(b))
+  for (_, funs) in get_bubbles(b)
+    for (w, α, rest...) in funs
+      nodes[w] = sum(α[i]/K * verts[i] for i in eachindex(α))
+    end
+  end
+  compute_own_nodes_permutations(poly, nodes)
+end
+
+# Compute get_face_own_dofs_permutations data for a moment-based element.
+# result[f][fpid] is the permutation of face-local DOF indices for face f
+# under face-permutation id fpid. Face-local DOFs are ordered by moment index.
+function moment_face_own_dofs_permutations(
+  face_to_nfpids,
+  face_to_moments,
+  moment_to_fpid_to_perm
+)
+  map(1:length(face_to_nfpids)) do f
+    nfpids     = face_to_nfpids[f]
+    moments_on = face_to_moments[f]
+    isempty(moments_on) && return [Int[] for _ in 1:nfpids]
+    block_sizes = [length(first(moment_to_fpid_to_perm[m])) for m in moments_on]
+    n_own_dofs  = sum(block_sizes)
+    map(1:nfpids) do fpid
+      perm   = Vector{Int}(undef, n_own_dofs)
+      offset = 0
+      for (k, m) in enumerate(moments_on)
+        bs    = block_sizes[k]
+        lperm = moment_to_fpid_to_perm[m][fpid]
+        for j in 1:bs
+          perm[offset + j] = offset + lperm[j]
+        end
+        offset += bs
+      end
+      perm
+    end
+  end
+end
+
+# Build face_own_dofs_permutations from the raw moments triplets and the
+# precomputed reffaces/face_types tables (already available in the constructor).
+function moment_face_own_dofs_permutations(p::Polytope, moments, reffaces, face_types)
+  nfaces    = num_faces(p)
+  n_moments = length(moments)
+  face_to_nfpids         = map(length, get_face_vertex_permutations(p))
+  face_to_moments        = [Int[] for _ in 1:nfaces]
+  moment_to_fpid_to_perm = Vector{Vector{Vector{Int}}}(undef, n_moments)
+  for (m, (faces, σ, μ)) in enumerate(moments)
+    for f in faces
+      push!(face_to_moments[f], m)
+    end
+    fp = reffaces[face_types[first(faces)]]
+    moment_to_fpid_to_perm[m] = get_basis_permutations(fp, μ)
+  end
+  moment_face_own_dofs_permutations(face_to_nfpids, face_to_moments, moment_to_fpid_to_perm)
+end
+
 # MomentBasedDofBasis
 
 struct Moment <: Dof end
@@ -17,14 +143,17 @@ struct MomentBasedDofBasis{P,V,O,N} <: AbstractVector{Moment}
   face_moments::Vector{Array{V,N}}
   face_nodes::Vector{UnitRange{Int}}
   face_own_moms::Vector{Vector{Int}}
+  face_own_moms_permutations::Vector{Vector{Vector{Int}}}
   operator::O
 
-  function MomentBasedDofBasis(nodes,f_moments,f_nodes,f_own_moms,operator=nothing)
+  function MomentBasedDofBasis(
+    nodes, f_moments, f_nodes, f_own_moms, f_own_dofs_permutations, operator = nothing
+  )
     P = eltype(nodes)
     V = eltype(eltype(f_moments))
     N = ndims(first(f_moments))
     O = typeof(operator)
-    new{P,V,O,N}(nodes,f_moments,f_nodes,f_own_moms,operator)
+    new{P,V,O,N}(nodes, f_moments, f_nodes, f_own_moms, f_own_dofs_permutations, operator)
   end
 
   # Unused and untested
@@ -76,6 +205,8 @@ the face the moment was defined on.
 """
 get_face_own_moments(b::MomentBasedDofBasis) = b.face_own_moms
 
+get_face_own_moments_permutations(b::MomentBasedDofBasis) = b.face_own_moms_permutations
+
 """
     get_face_nodes_dofs(b::MomentBasedDofBasis)
 
@@ -178,42 +309,44 @@ function evaluate!(cache, b::MomentBasedDofBasis, field::AbstractVector{<:Field}
   return dofs
 end
 
-"""
-    MomentBasedDofBasis( p::Polytope, prebasis::AbstractVector{<:Field}, moments,
-        [, face_own_dofs], operator=nothing
-    )
-
-Creates a basis of DoFs defined by moments on faces of `p`.
-
-`moments` is a vector of moment descriptors, each one is given by a triplet
-(f,σ,μ) where
-  - f is collection of ids of faces Fₖ of `p`, that index `get_faces(p)`,
-  - σ is a function σ(φ,μ,ds) **linear** in φ and μ that takes two Field-vectors φ and μ and a `FaceMeasure` ds and returns a Field-like object to be integrated over each face Fₖ,
-  - μ is a polynomials basis on Fₖ.
-
-The moment DoFs are thus defined by φ -> ∫_Fₖ σ(φ,μᵢ,ds)dFₖ,  ∀ σ,k,i.
-In the final basis, DoFs are ordered by moment, then by face, then by "test" polynomial.
-
-All the faces in a moment must be of the same type (have same reference face).
-
-If an `operator` function -- e.g. `∇` -- is given, it is applied to `φ` (with
-respect to `p`'s coordinates) before being passed to `σ`. The moment becomes
-`φ -> σ(∇φ,μ,ds)`.
-
-If `face_own_dofs` is given, it defines the moment ownership to faces.
-"""
-function MomentBasedDofBasis(
-    p::Polytope{D},
-    prebasis::AbstractVector{<:Field},
-    moments::AbstractVector{<:Tuple},
-    face_own_dofs::AbstractVector,
-    operator=nothing,
-  ) where D
-
-  dofs = MomentBasedDofBasis(p, prebasis, moments, operator)
-  # replace moment ownership
-  MomentBasedDofBasis(dofs.nodes, dofs.face_moments, dofs.face_nodes, face_own_dofs, operator)
-end
+# TODO: Is this ever used or needed?
+# """
+#     MomentBasedDofBasis( p::Polytope, prebasis::AbstractVector{<:Field}, moments,
+#         [, face_own_dofs], operator=nothing
+#     )
+#
+# Creates a basis of DoFs defined by moments on faces of `p`.
+#
+# `moments` is a vector of moment descriptors, each one is given by a triplet
+# (f,σ,μ) where
+#   - f is collection of ids of faces Fₖ of `p`, that index `get_faces(p)`,
+#   - σ is a function σ(φ,μ,ds) **linear** in φ and μ that takes two Field-vectors φ and μ and a `FaceMeasure` ds and returns a Field-like object to be integrated over each face Fₖ,
+#   - μ is a polynomials basis on Fₖ.
+#
+# The moment DoFs are thus defined by φ -> ∫_Fₖ σ(φ,μᵢ,ds)dFₖ,  ∀ σ,k,i.
+# In the final basis, DoFs are ordered by moment, then by face, then by "test" polynomial.
+#
+# All the faces in a moment must be of the same type (have same reference face).
+#
+# If an `operator` function -- e.g. `∇` -- is given, it is applied to `φ` (with
+# respect to `p`'s coordinates) before being passed to `σ`. The moment becomes
+# `φ -> σ(∇φ,μ,ds)`.
+#
+# If `face_own_dofs` is given, it defines the moment ownership to faces.
+# """
+# function MomentBasedDofBasis(
+#     p::Polytope{D},
+#     prebasis::AbstractVector{<:Field},
+#     moments::AbstractVector{<:Tuple},
+#     face_own_dofs::AbstractVector,
+#     operator=nothing,
+#   ) where D
+#   dofs = MomentBasedDofBasis(p, prebasis, moments, operator)
+#   MomentBasedDofBasis(
+#     dofs.nodes, dofs.face_moments, dofs.face_nodes, face_own_dofs,
+#     dofs.face_own_dofs_permutations, operator
+#   )
+# end
 
 function MomentBasedDofBasis(
     p::Polytope{D},
@@ -298,7 +431,11 @@ function MomentBasedDofBasis(
     end
   end
 
-  MomentBasedDofBasis(nodes, face_moments, face_nodes, face_own_moms, operator)
+  face_own_moms_permutations = moment_face_own_dofs_permutations(
+    p, moments, reffaces, face_types
+  )
+
+  MomentBasedDofBasis(nodes, face_moments, face_nodes, face_own_moms, face_own_moms_permutations, operator)
 end
 
 # Unused and untested
@@ -512,24 +649,25 @@ function MomentBasedReferenceFE(
     shapefuns = apply_face_signflip(shapefuns, p, conformity)
     dofs = compute_dofs(predofs, shapefuns)
     face_own_dofs = get_face_own_funs(prebasis, p, conformity)
+    face_own_dofs_permutations = get_face_own_funs_permutations(prebasis, p, conformity)
 
     return GenericRefFE{typeof(name)}(
-      n_dofs, p, predofs, conformity, metadata, face_own_dofs, shapefuns, dofs
+      n_dofs, p, predofs, conformity, metadata, face_own_dofs, shapefuns, dofs, face_own_dofs_permutations
     )
   end # else, standard prebasis inversion
 
-  if isnothing(face_own_dofs)
-    face_own_dofs = get_face_own_moments(dof_basis)
-  end
+  @check isnothing(face_own_dofs) # TODO document this is forbiden for change_dof=false
+  face_own_dofs = get_face_own_moments(dof_basis)
+  face_own_dofs_permutations = get_face_own_moments_permutations(dof_basis)
   GenericRefFE{typeof(name)}(
-    n_dofs, p, prebasis, dof_basis, conformity, metadata, face_own_dofs
+    n_dofs, p, prebasis, dof_basis, conformity, metadata, face_own_dofs, face_own_dofs_permutations
   )
 end
 
 # Default polynomial type for moment based reference FEs
 function _mom_reffe_default_PT(p)
   is_simplex(p) && return Bernstein
-  is_n_cube(p) && return Polynomials.ModalC0
+  is_n_cube(p) && return Bernstein
   Monomial
 end
 
diff --git a/src/ReferenceFEs/NedelecRefFEs.jl b/src/ReferenceFEs/NedelecRefFEs.jl
index 8d5c71453..3ec453709 100644
--- a/src/ReferenceFEs/NedelecRefFEs.jl
+++ b/src/ReferenceFEs/NedelecRefFEs.jl
@@ -1,7 +1,7 @@
 """
-    struct Nedelec{kind} <: ReferenceFEName
+    struct Nedelec{kind} <: MomentBasedRefFEName
 """
-struct Nedelec{kind} <: ReferenceFEName
+struct Nedelec{kind} <: MomentBasedRefFEName
   Nedelec{1}() = new{1}()
   Nedelec{2}() = new{2}()
 end
diff --git a/src/ReferenceFEs/Polytopes.jl b/src/ReferenceFEs/Polytopes.jl
index 9c85e52ed..ef803df58 100644
--- a/src/ReferenceFEs/Polytopes.jl
+++ b/src/ReferenceFEs/Polytopes.jl
@@ -204,38 +204,38 @@ end
     get_vertex_permutations(p::Polytope) -> Vector{Vector{Int}}
 
 Returns a vector of vectors containing all admissible permutations of the
-vertices of `p`. An admissible permutation is one such that, if the vertices of
-the polytope are re-labeled according to this permutation, the resulting
-polytope preserves the shape of the original one.
+vertices of `p`. A vertex permutation is "admissible" if it is the image of the
+polytope's vertices by a diffeomorphism of the polytope to itself, like a typical
+geometrical mapping.
 
-These permutation are used to iddentify every possile way a geometrical map may
+These permutation are used to e.g. iddentify every possile way a geometrical map may
 permute the vertices of the boundary faces of a reference polytope (of dimension
-≤ 3) into the physical one. Only the iddentity permutation `[1, 2, …, N]` is
-returned for 3D polytopes with `N` vertices. Indeed, a (3≥D)-dimensional
-physical polytope `𝓟` is only mapped by one physical map, but its faces are
-mapped by the physical map of all adjascent elements.
+≤ 3) into the physical one.
 
 # Examples
 
 ```jldoctest
 using Gridap.ReferenceFEs
 
-perms = get_vertex_permutations(SEGMENT)
+perms = get_vertex_permutations(SEGMENT);
 println(perms)
 
 # output
 Array{Int,1}[[1, 2], [2, 1]]
 
-perms = get_vertex_permutations(TET)
-println(perms)
+perms = get_vertex_permutations(QUAD);
+println(length(perms))
 
 # output
-Array{Int,1}[[1, 2], [2, 1]]
+8
 ```
 The first admissible permutation for a segment is `[1,2]`,i.e., the identity.
 The second one is `[2,1]`, i.e., the first vertex is relabeled as `2` and the
 second vertex is relabeled as `1`.
 
+The QUAD has only 8 admissible vertex permutations although there are 24 possible
+permutations of its vertices, but 16 of them "pinch" the square into an "hourglass",
+creating a point where the map is not differentiable.
 """
 function get_vertex_permutations(p::Polytope)
   @abstractmethod
diff --git a/src/ReferenceFEs/RaviartThomasRefFEs.jl b/src/ReferenceFEs/RaviartThomasRefFEs.jl
index 7acc47d81..9f744ca70 100644
--- a/src/ReferenceFEs/RaviartThomasRefFEs.jl
+++ b/src/ReferenceFEs/RaviartThomasRefFEs.jl
@@ -4,7 +4,7 @@
 """
     struct RaviartThomas <: ReferenceFEName
 """
-struct RaviartThomas <: ReferenceFEName end
+struct RaviartThomas <: MomentBasedRefFEName end
 
 """
     const raviart_thomas = RaviartThomas()
diff --git a/src/ReferenceFEs/ReferenceFEInterfaces.jl b/src/ReferenceFEs/ReferenceFEInterfaces.jl
index 29da76a2a..579be5816 100644
--- a/src/ReferenceFEs/ReferenceFEInterfaces.jl
+++ b/src/ReferenceFEs/ReferenceFEInterfaces.jl
@@ -283,7 +283,7 @@ function get_face_own_dofs(reffe::ReferenceFE, conf::Nothing)
   get_face_own_dofs(reffe)
 end
 
-function _get_face_own_dofs_l2(reffe::ReferenceFE)
+function l2_face_own_dofs(reffe::ReferenceFE)
   p = get_polytope(reffe)
   r = [Int[] for i in 1:num_faces(p)]
   r[end] = collect(1:num_dofs(reffe))
@@ -311,13 +311,21 @@ anisotropic order). For these permutations the corresponding vector in
 """
 function get_face_own_dofs_permutations(reffe::ReferenceFE, conf::Conformity)
   face_own_dofs = get_face_own_dofs(reffe, conf)
-  _trivial_face_own_dofs_permutations(face_own_dofs)
+  l2_face_own_dofs_permutations(face_own_dofs)
+  #empty_face_own_dofs_permutations(face_own_dofs)
 end
 
-function _trivial_face_own_dofs_permutations(face_own_dofs)
+function l2_face_own_dofs_permutations(face_own_dofs)
   [[collect(Int, 1:length(dofs)),] for dofs in face_own_dofs]
 end
 
+"""
+    empty_face_own_dofs_permutations(face_own_dofs) = [Vector{Int}[] for _ in face_own_dofs]
+"""
+function empty_face_own_dofs_permutations(face_own_dofs)
+  [Vector{Int}[] for _ in face_own_dofs]
+end
+
 function get_face_own_dofs_permutations(reffe::ReferenceFE)
   conf = Conformity(reffe)
   get_face_own_dofs_permutations(reffe, conf)
@@ -610,9 +618,8 @@ This type is a *materialization* of the `ReferenceFE` interface. That is, it is
 This type is useful to build reference FEs from the underlying ingredients without
 the need to create a new type.
 
-Note that some fields in this `struct` are type unstable deliberately in order to simplify the
-type signature. Don't access them in computationally expensive functions,
-instead extract the required fields before and pass them to the computationally expensive function.
+Note that some fields in this `struct` are abstractly typed deliberately in
+order to simplify the type signature.
 """
 @ahe struct GenericRefFE{T,D} <: ReferenceFE{D}
   ndofs::Int
@@ -622,6 +629,7 @@ instead extract the required fields before and pass them to the computationally
   conformity::Conformity
   metadata
   face_own_dofs::Vector{Vector{Int}}
+  face_own_dofs_permutations::Vector{Vector{Vector{Int}}}
   shapefuns::AbstractVector{<:Field}
 
   @doc """
@@ -633,10 +641,13 @@ instead extract the required fields before and pass them to the computationally
           conformity::Conformity,
           metadata,
           face_own_dofs::Vector{Vector{Int}},
-          shapefuns::AbstractVector{<:Field}=compute_shapefuns(dofs,prebasis)
-        ) where {T,D}
+         [shapefuns::AbstractVector{<:Field} = compute_shapefuns(dofs,prebasis), ]
+         [face_own_dofs_permutations::Vector{Vector{Vector{Int}}} = default, ]
+        )
 
   Constructor using a DoF basis and function space pre-basis.
+
+  The default `face_own_dofs_permutations` is [`l2_face_own_dofs_permutations(face_own_dofs)`](@ref).
   """
   function GenericRefFE{T}(
     ndofs::Int,
@@ -646,7 +657,9 @@ instead extract the required fields before and pass them to the computationally
     conformity::Conformity,
     metadata,
     face_own_dofs::Vector{Vector{Int}},
-    shapefuns::AbstractVector{<:Field}=compute_shapefuns(dofs, prebasis)) where {T,D}
+    shapefuns::AbstractVector{<:Field}=compute_shapefuns(dofs, prebasis),
+    face_own_dofs_permutations::Vector{Vector{Vector{Int}}}=l2_face_own_dofs_permutations(face_own_dofs),
+  ) where {T,D}
 
     new{T,D}(
       ndofs,
@@ -656,8 +669,11 @@ instead extract the required fields before and pass them to the computationally
       conformity,
       metadata,
       face_own_dofs,
-      shapefuns)
+      face_own_dofs_permutations,
+      shapefuns,
+    )
   end
+
   @doc """
         GenericRefFE{T}(
           ndofs::Int,
@@ -668,10 +684,13 @@ instead extract the required fields before and pass them to the computationally
           metadata,
           face_own_dofs::Vector{Vector{Int}},
           shapefuns::AbstractVector{<:Field},
-          dofs::AbstractVector{<:Dof}=compute_dofs(predofs,shapefuns)
-        ) where {T,D}
+         [dofs::AbstractVector{<:Dof}=compute_dofs(predofs,shapefuns), ]
+         [face_own_dofs_permutations::Vector{Vector{Vector{Int}}} = default, ]
+        )
 
   Constructor using shape function basis and a DoF pre-basis.
+
+  The default `face_own_dofs_permutations` is [`l2_face_own_dofs_permutations(face_own_dofs)`](@ref).
   """
   function GenericRefFE{T}(
     ndofs::Int,
@@ -681,7 +700,9 @@ instead extract the required fields before and pass them to the computationally
     metadata,
     face_own_dofs::Vector{Vector{Int}},
     shapefuns::AbstractVector{<:Field},
-    dofs::AbstractVector{<:Dof}=compute_dofs(predofs, shapefuns)) where {T,D}
+    dofs::AbstractVector{<:Dof}=compute_dofs(predofs, shapefuns),
+    face_own_dofs_permutations::Vector{Vector{Vector{Int}}}=l2_face_own_dofs_permutations(face_own_dofs),
+  ) where {T,D}
 
     new{T,D}(
       ndofs,
@@ -691,14 +712,40 @@ instead extract the required fields before and pass them to the computationally
       conformity,
       metadata,
       face_own_dofs,
+      face_own_dofs_permutations,
       # Trick to be able to eval dofs af shapefuns in physical space
       # cf /test/FESpacesTests/PhysicalFESpacesTests.jl
       #linear_combination(Eye{Int}(ndofs), shapefuns)
-      linear_combination(Diagonal(ones(ndofs)), shapefuns)
+      linear_combination(Diagonal(ones(ndofs)), shapefuns),
     )
   end
 end
 
+function GenericRefFE{T}(
+  ndofs, polytope, prebasis, dofs, conformity::Conformity, metadata, face_own_dofs,
+  face_own_dofs_permutations::Vector{Vector{Vector{Int}}},
+) where T
+
+  shapefuns = compute_shapefuns(dofs, prebasis)
+  GenericRefFE{T}(
+    ndofs, polytope, prebasis, dofs, conformity, metadata, face_own_dofs,
+    shapefuns, face_own_dofs_permutations,
+  )
+end
+
+function GenericRefFE{T}(
+  ndofs, polytope, predofs, conformity::Conformity, metadata, face_own_dofs, shapefuns,
+  face_own_dofs_permutations::Vector{Vector{Vector{Int}}},
+) where T
+
+  dofs = compute_dofs(predofs, shapefuns)
+  GenericRefFE{T}(
+    ndofs, polytope, predofs, conformity, metadata, face_own_dofs, shapefuns,
+    dofs, face_own_dofs_permutations,
+  )
+end
+
+
 get_name(::Type{<:GenericRefFE{Name}}) where Name = Name()
 
 num_dofs(reffe::GenericRefFE) = reffe.ndofs
@@ -714,7 +761,7 @@ Conformity(reffe::GenericRefFE) = reffe.conformity
 get_face_own_dofs(reffe::GenericRefFE) = reffe.face_own_dofs
 function get_face_own_dofs(reffe::GenericRefFE, conf::Conformity)
   conf == Conformity(reffe) && return reffe.face_own_dofs
-  conf isa L2Conformity && return _get_face_own_dofs_l2(reffe)
+  conf isa L2Conformity && return l2_face_own_dofs(reffe)
   @unreachable """\n
   It is not possible to use conformity $conf on this reference FE.
   """
@@ -726,6 +773,15 @@ function get_face_dofs(reffe::GenericRefFE)
   face_own_data_to_face_data(poly, face_own_dofs)
 end
 
+get_face_own_dofs_permutations(reffe::GenericRefFE) = reffe.face_own_dofs_permutations
+function get_face_own_dofs_permutations(reffe::GenericRefFE, conf::Conformity)
+  conf == Conformity(reffe) && return reffe.face_own_dofs_permutations
+  conf isa L2Conformity && return l2_face_own_dofs_permutations(get_face_own_dofs(reffe, conf))
+  @unreachable """\n
+  It is not possible to use conformity $conf on this reference FE.
+  """
+end
+
 get_shapefuns(reffe::GenericRefFE) = reffe.shapefuns
 
 get_metadata(reffe::GenericRefFE) = reffe.metadata
diff --git a/src/ReferenceFEs/ReferenceFEs.jl b/src/ReferenceFEs/ReferenceFEs.jl
index e01b5b734..9642a3c70 100644
--- a/src/ReferenceFEs/ReferenceFEs.jl
+++ b/src/ReferenceFEs/ReferenceFEs.jl
@@ -12,7 +12,7 @@ using Test
 using DocStringExtensions
 using LinearAlgebra
 using StaticArrays
-using Combinatorics
+using Combinatorics: permutations, combinations
 using FillArrays
 using AutoHashEquals: @auto_hash_equals as @ahe
 using ..Gridap
@@ -121,6 +121,7 @@ export test_dof_array
 
 export ReferenceFE
 export ReferenceFEName
+export MomentBasedRefFEName
 export GenericRefFE
 export get_name
 export get_polytope
diff --git a/test/ReferenceFEsTests/CLagrangianRefFEsTests.jl b/test/ReferenceFEsTests/CLagrangianRefFEsTests.jl
index ad9c474da..abb501b72 100644
--- a/test/ReferenceFEsTests/CLagrangianRefFEsTests.jl
+++ b/test/ReferenceFEsTests/CLagrangianRefFEsTests.jl
@@ -116,10 +116,10 @@ reffe = LagrangianRefFE(VectorValue{2,Float64},SEGMENT,orders)
 
 orders = (2,3)
 reffe = LagrangianRefFE(VectorValue{2,Float64},QUAD,orders)
-@test get_own_nodes_permutations(reffe) ==[[1, 2], [0, 0], [1, 2], [0, 0], [0, 0], [2, 1], [0, 0], [2, 1]]
+@test get_own_nodes_permutations(reffe) == [[1, 2], [0, 0], [1, 2], [0, 0], [2, 1], [0, 0], [2, 1], [0, 0]]
 @test get_own_dofs_permutations(reffe) == [
   [1, 2, 3, 4], [0, 0, 0, 0], [1, 2, 3, 4], [0, 0, 0, 0],
-  [0, 0, 0, 0], [2, 1, 4, 3], [0, 0, 0, 0], [2, 1, 4, 3]]
+  [2, 1, 4, 3], [0, 0, 0, 0], [2, 1, 4, 3], [0, 0, 0, 0]]
 
 reffe = LagrangianRefFE(VectorValue{2,Float64},QUAD,2)
 @test get_node_and_comp_to_dof(reffe) == VectorValue{2,Int}[
diff --git a/test/ReferenceFEsTests/ExtrusionPolytopesTests.jl b/test/ReferenceFEsTests/ExtrusionPolytopesTests.jl
index 99d304bb9..775ee3bd8 100644
--- a/test/ReferenceFEsTests/ExtrusionPolytopesTests.jl
+++ b/test/ReferenceFEsTests/ExtrusionPolytopesTests.jl
@@ -33,7 +33,7 @@ r = Point{2,Float64}[(1, 0), (1, 0), (0, 1), (0, 1)]
 r = [
   [1, 2, 3, 4], [1, 3, 2, 4], [2, 1, 4, 3], [2, 4, 1, 3],
   [3, 1, 4, 2], [3, 4, 1, 2], [4, 2, 3, 1], [4, 3, 2, 1]]
-@test get_vertex_permutations(p) == r
+@test Set(get_vertex_permutations(p)) == Set(r)
 
 r = Point{2,Float64}[(0, -1), (0, 1), (-1, 0), (1, 0)]
 @test get_facet_normal(p) == r
@@ -72,10 +72,47 @@ test_polytope(TET,optional=true)
 test_polytope(HEX,optional=true)
 
 perm = get_vertex_permutations(TRI)
-@test perm == [[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1]]
+@test Set(perm) == Set([[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1]])
 
 perm = get_vertex_permutations(QUAD)
-@test isa(perm,Vector{Vector{Int}})
+@test Set(perm) == Set([[1, 2, 3, 4],
+                        [1, 3, 2, 4],
+                        [2, 1, 4, 3],
+                        [2, 4, 1, 3],
+                        [3, 1, 4, 2],
+                        [3, 4, 1, 2],
+                        [4, 2, 3, 1],
+                        [4, 3, 2, 1]])
+
+perm = get_vertex_permutations(PYRAMID)
+@test Set(perm) == Set([[1, 2, 3, 4, 5],
+                        [1, 3, 2, 4, 5],
+                        [2, 1, 4, 3, 5],
+                        [2, 4, 1, 3, 5],
+                        [3, 1, 4, 2, 5],
+                        [3, 4, 1, 2, 5],
+                        [4, 2, 3, 1, 5],
+                        [4, 3, 2, 1, 5]])
+
+#[1, 2, 3] ⊗ {[1, 2], [2, 1]}
+#[1, 3, 2]
+#[2, 1, 3]
+#[2, 3, 1]
+#[3, 1, 2]
+#[3, 2, 1]
+perm = get_vertex_permutations(WEDGE)
+@test Set(perm) == Set([[1, 2, 3, 4, 5, 6],
+                        [1, 3, 2, 4, 6, 5],
+                        [2, 1, 3, 5, 4, 6],
+                        [2, 3, 1, 5, 6, 4],
+                        [3, 1, 2, 6, 4, 5],
+                        [3, 2, 1, 6, 5, 4],
+                        [4, 5, 6, 1, 2, 3],
+                        [4, 6, 5, 1, 3, 2],
+                        [5, 4, 6, 2, 1, 3],
+                        [5, 6, 4, 2, 3, 1],
+                        [6, 4, 5, 3, 1, 2],
+                        [6, 5, 4, 3, 2, 1]])
 
 @test num_facets(SEGMENT) == 2
 @test num_facets(TRI) == 3
@@ -105,20 +142,32 @@ iface_to_ftype = get_face_type(WEDGE,d)
 @test length(reffaces) == 1
 @test iface_to_ftype == [1, 1, 1, 1, 1, 1, 1, 1, 1]
 
-@test get_face_vertex_permutations(QUAD,1) == [[[1,2],[2,1]],[[1,2],[2,1]],[[1,2],[2,1]],[[1,2],[2,1]]]
-@test get_face_vertex_permutations(QUAD,2) == [[[1,2,3,4],[1,3,2,4],[2,1,4,3],[2,4,1,3],[3,1,4,2],[3,4,1,2],[4,2,3,1],[4,3,2,1]]]
+expected = [[[1,2],[2,1]],[[1,2],[2,1]],[[1,2],[2,1]],[[1,2],[2,1]]]
+got = get_face_vertex_permutations(QUAD,1)
+@test all(Set(a) == Set(b) for (a,b) in zip(got, expected))
 
-r = [
+expected = [[[1,2,3,4],[1,3,2,4],[2,1,4,3],[2,4,1,3],[3,1,4,2],[3,4,1,2],[4,2,3,1],[4,3,2,1]]]
+got = get_face_vertex_permutations(QUAD,2)
+@test all(Set(a) == Set(b) for (a,b) in zip(got, expected))
+
+expected = [
   [[1]],[[1]],[[1]],[[1]],
   [[1,2],[2,1]],[[1,2],[2,1]],[[1,2],[2,1]],[[1,2],[2,1]],
   [[1,2,3,4],[1,3,2,4],[2,1,4,3],[2,4,1,3],[3,1,4,2],[3,4,1,2],[4,2,3,1],[4,3,2,1]]]
-@test get_face_vertex_permutations(QUAD) == r
+got = get_face_vertex_permutations(QUAD)
+@test all(Set(a) == Set(b) for (a,b) in zip(got, expected))
+
+expected = [[[1,2],[2,1]],[[1,2],[2,1]],[[1,2],[2,1]]]
+got = get_face_vertex_permutations(TRI,1)
+@test all(Set(a) == Set(b) for (a,b) in zip(got, expected))
 
-r = [[[1,2],[2,1]],[[1,2],[2,1]],[[1,2],[2,1]]]
-@test get_face_vertex_permutations(TRI,1) == r
+expected = [[[1, 2, 3, 4, 5, 6, 7, 8], [1, 2, 5, 6, 3, 4, 7, 8], [1, 3, 2, 4, 5, 7, 6, 8], [1, 3, 5, 7, 2, 4, 6, 8], [1, 5, 2, 6, 3, 7, 4, 8], [1, 5, 3, 7, 2, 6, 4, 8], [2, 1, 4, 3, 6, 5, 8, 7], [2, 1, 6, 5, 4, 3, 8, 7], [2, 4, 1, 3, 6, 8, 5, 7], [2, 4, 6, 8, 1, 3, 5, 7], [2, 6, 1, 5, 4, 8, 3, 7], [2, 6, 4, 8, 1, 5, 3, 7], [3, 1, 4, 2, 7, 5, 8, 6], [3, 1, 7, 5, 4, 2, 8, 6], [3, 4, 1, 2, 7, 8, 5, 6], [3, 4, 7, 8, 1, 2, 5, 6], [3, 7, 1, 5, 4, 8, 2, 6], [3, 7, 4, 8, 1, 5, 2, 6], [4, 2, 3, 1, 8, 6, 7, 5], [4, 2, 8, 6, 3, 1, 7, 5], [4, 3, 2, 1, 8, 7, 6, 5], [4, 3, 8, 7, 2, 1, 6, 5], [4, 8, 2, 6, 3, 7, 1, 5], [4, 8, 3, 7, 2, 6, 1, 5], [5, 1, 6, 2, 7, 3, 8, 4], [5, 1, 7, 3, 6, 2, 8, 4], [5, 6, 1, 2, 7, 8, 3, 4], [5, 6, 7, 8, 1, 2, 3, 4], [5, 7, 1, 3, 6, 8, 2, 4], [5, 7, 6, 8, 1, 3, 2, 4], [6, 2, 5, 1, 8, 4, 7, 3], [6, 2, 8, 4, 5, 1, 7, 3], [6, 5, 2, 1, 8, 7, 4, 3], [6, 5, 8, 7, 2, 1, 4, 3], [6, 8, 2, 4, 5, 7, 1, 3], [6, 8, 5, 7, 2, 4, 1, 3], [7, 3, 5, 1, 8, 4, 6, 2], [7, 3, 8, 4, 5, 1, 6, 2], [7, 5, 3, 1, 8, 6, 4, 2], [7, 5, 8, 6, 3, 1, 4, 2], [7, 8, 3, 4, 5, 6, 1, 2], [7, 8, 5, 6, 3, 4, 1, 2], [8, 4, 6, 2, 7, 3, 5, 1], [8, 4, 7, 3, 6, 2, 5, 1], [8, 6, 4, 2, 7, 5, 3, 1], [8, 6, 7, 5, 4, 2, 3, 1], [8, 7, 4, 3, 6, 5, 2, 1], [8, 7, 6, 5, 4, 3, 2, 1]]]
+got = get_face_vertex_permutations(HEX,3)
+@test all(Set(a) == Set(b) for (a,b) in zip(got, expected))
 
-@test get_face_vertex_permutations(HEX,3) == [[[1, 2, 3, 4, 5, 6, 7, 8]]]
-@test get_face_vertex_permutations(TET,3) == [[[1, 2, 3, 4]]]
+expected = [[[1, 2, 3, 4], [1, 2, 4, 3], [1, 3, 2, 4], [1, 3, 4, 2], [1, 4, 2, 3], [1, 4, 3, 2], [2, 1, 3, 4], [2, 1, 4, 3], [2, 3, 1, 4], [2, 3, 4, 1], [2, 4, 1, 3], [2, 4, 3, 1], [3, 1, 2, 4], [3, 1, 4, 2], [3, 2, 1, 4], [3, 2, 4, 1], [3, 4, 1, 2], [3, 4, 2, 1], [4, 1, 2, 3], [4, 1, 3, 2], [4, 2, 1, 3], [4, 2, 3, 1], [4, 3, 1, 2], [4, 3, 2, 1]]]
+got = get_face_vertex_permutations(TET,3)
+@test all(Set(a) == Set(b) for (a,b) in zip(got, expected))
 
 v,p = simplexify(QUAD)
 @test p == TRI
@@ -179,21 +228,13 @@ for dim in 0:5
     @test is_simplex(p)
 end
 
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+@test SEGMENT == ExtrusionPolytope(HEX_AXIS) == ExtrusionPolytope(TET_AXIS)
+@test TRI == ExtrusionPolytope(HEX_AXIS, TET_AXIS)
+@test QUAD == ExtrusionPolytope(TET_AXIS, HEX_AXIS)
+@test TET == ExtrusionPolytope(HEX_AXIS, TET_AXIS, TET_AXIS)
+@test HEX == ExtrusionPolytope(TET_AXIS, HEX_AXIS, HEX_AXIS)
+@test PYRAMID == ExtrusionPolytope(TET_AXIS, HEX_AXIS, TET_AXIS)
+@test WEDGE == ExtrusionPolytope(HEX_AXIS, TET_AXIS, HEX_AXIS)
 
 #using Gridap.Fields
 #
diff --git a/test/ReferenceFEsTests/ModalScalarRefFEsTests.jl b/test/ReferenceFEsTests/ModalScalarRefFEsTests.jl
index bee6cd4f7..984ed682e 100644
--- a/test/ReferenceFEsTests/ModalScalarRefFEsTests.jl
+++ b/test/ReferenceFEsTests/ModalScalarRefFEsTests.jl
@@ -14,8 +14,9 @@ reffe = ModalScalarRefFE(Float64,p,2; F)
 @test reffe == ReferenceFE(p,F,2,0)
 @test reffe == ReferenceFE(p,F,2,0; nodal)
 
-face_own_dofs = Vector{Int}[[1],[2],[3],[4],[5],[6],[]]
-face_dofs = Vector{Int}[[1],[2],[3],[1,2,4],[1,3,5],[2,3,6],[1,2,3,4,5,6]]
+# corresponds to ownership of BernsteinBasisOnSimplex
+face_own_dofs = Vector{Int}[[1], [4], [6], [2], [3], [5], []]
+face_dofs = Vector{Int}[[1], [4], [6], [1, 4, 2], [1, 6, 3], [4, 6, 5], [1, 4, 6, 2, 3, 5]]
 @test get_face_own_dofs(reffe) == face_own_dofs
 @test get_face_dofs(reffe) == face_dofs
 
@@ -54,7 +55,7 @@ test_reference_fe(reffe)
 @test Conformity(reffe,:H1) == GradConformity()
 @test Conformity(reffe,:L2) == L2Conformity()
 
-@test_throws "hierarchical" reffe = ModalScalarRefFE(Float64,HEX,order; F=:S, poly_type=Bernstein)
+@test_throws "hierarchical" ModalScalarRefFE(Float64,HEX,order; F=:S, poly_type=Bernstein)
 
 reffe = ModalScalarRefFE(Float64,HEX,order; F=:S, poly_type)
 test_reference_fe(reffe)
diff --git a/test/ReferenceFEsTests/ReferenceFEInterfacesTests.jl b/test/ReferenceFEsTests/ReferenceFEInterfacesTests.jl
index 3b5ef1b3e..b79f17152 100644
--- a/test/ReferenceFEsTests/ReferenceFEInterfacesTests.jl
+++ b/test/ReferenceFEsTests/ReferenceFEInterfacesTests.jl
@@ -1,9 +1,12 @@
 module ReferenceFEInterfacesTests
 
 using Test
+using Gridap.Arrays
+using Gridap.TensorValues
 using Gridap.Fields
 using Gridap.Polynomials
 using Gridap.ReferenceFEs
+using LinearAlgebra
 
 # Abstract API
 
@@ -91,4 +94,91 @@ dofs = get_dof_basis(reffe)
 
 @test evaluate(dofs,shapefuns) == [1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0]
 
+function test_reffe_permutations(reffe::ReferenceFE)
+  conf = Conformity(reffe)
+  p = get_polytope(reffe)
+  φ = get_shapefuns(reffe)
+  face_own_dofs = get_face_own_dofs(reffe)
+  face_own_dofs_perms = get_face_own_dofs_permutations(reffe)
+  face_dofs = get_face_dofs(reffe)
+
+  # construct a geometrical map onto p
+  p_reffe = LagrangianRefFE(p)
+  map_basis = get_shapefuns(p_reffe)
+  vertex_coords = copy(get_vertex_coordinates(p))
+  Φ = linear_combination(vertex_coords, map_basis)
+  Φ_inv = inverse_map(Φ)
+  Jt = Broadcasting(∇)(Φ)
+
+  # construct the pushed basis
+  pushforward = Pushforward(get_name(reffe), conf)
+  φ_pushed = if pushforward isa IdentityPiolaMap
+    φ
+  else
+    Broadcasting(Operation(pushforward))(φ, Jt)
+  end
+
+  # choose some points x of p at which we will compare the shape-functions values
+  r = get_order(reffe)
+  x = get_nodes(get_dof_basis(LagrangianRefFE(Float64, p, r+2)))
+  @assert length(x) > length(φ) "more evaluation points are needed to compare the bases"
+
+  _flatten_rows(φx::Array{<:MultiValue}) = reinterpret(Float64,φx)
+  _flatten_rows(φx::Array{<:Number}) = φx
+
+  φx = _flatten_rows(evaluate(φ, x))
+  cφ_pushed = return_cache(φ_pushed, x)
+  cΦ        = return_cache(Φ,        x)
+  cΦ_inv    = return_cache(Φ_inv,    x)
+  #cJ = return_cache(J, x)
+
+  n, m = length(x), length(φ)
+  D = num_indep_components(value_type(φ))
+  M = n*D
+  @assert M == size(φx, 1)
+  A = cholesky(φx'*φx)
+  C = zeros(Float64, m, m)
+
+#poly_perm = last(get_vertex_permutations(p))
+  for poly_perm in get_vertex_permutations(p)
+    # map from p to itself "permuted"
+    permute!(vertex_coords, poly_perm) # this changes Φ and Φ_inv
+    pushed_xq = evaluate!(cΦ_inv, Φ_inv, x);
+    @assert pushed_xq == evaluate!(cΦ_inv, inverse_map(Φ), x)
+    φpx = _flatten_rows(evaluate!(cφ_pushed, φ_pushed, pushed_xq));
+
+    # check that φ_pushed spans the same space as φ
+    B = φx'*φpx
+    ldiv!(C, A, B)
+    is_pushed_basis_same_space = φpx ≈ φx*C
+    is_change_signed_perm = C'*C ≈ I
+    @show poly_perm
+    @show is_pushed_basis_same_space #@test
+    @show is_change_signed_perm      #@test
+    println()
+
+    for (face, (own_dof, perms)) in enumerate(zip(face_own_dofs, face_own_dofs_perms))
+      @show face, (own_dof, perms)
+      # iddentify the permutation of face
+
+      for dof in own_dof
+      end
+    end
+
+    # reset Φ and Φ_inv
+    invpermute!(vertex_coords, poly_perm)
+  end
+  println()
+
+end
+
+reffe = ReferenceFE(TRI, :P, 3, 0)
+test_reffe_permutations(reffe)
+
+reffe = ReferenceFE(TRI, :P, 2, 1)
+test_reffe_permutations(reffe)
+
+reffe = ReferenceFE(TRI, :P⁻, 2, 1)
+test_reffe_permutations(reffe)
+
 end # module
diff --git a/test/darcy_mwe.jl b/test/darcy_mwe.jl
new file mode 100644
index 000000000..6ac1e2ec4
--- /dev/null
+++ b/test/darcy_mwe.jl
@@ -0,0 +1,44 @@
+using Gridap
+using Gridap.Geometry
+using Gridap.FESpaces
+
+D = 2
+order = 2
+domain = ifelse(D==2, (0,1,0,1), (0,1,0,1,0,1))
+partition = Tuple(fill(3, D))
+model = simplexify(CartesianDiscreteModel(domain,partition); positive=true)
+
+Ω  = Triangulation(model)
+dΩ = Measure(Ω, 2*(order+1))
+
+# u_exact: degree-2, div-free, in [P2]³ ⊂ RT2.
+#   Derived from stream potential A=(0,0,x²y): curl(A) = (x², -2xy, 0), div=0.
+# p_exact: degree-2, in P2.
+# Body force g = u_exact + ∇p_exact makes (u_exact,p_exact) the exact Darcy solution.
+u_exact(x) = ifelse(D==2,VectorValue(x[1]^2, -2*x[1]*x[2]), VectorValue(x[1]^2, -2*x[1]*x[2], 0.0))
+p_tilde(x) = x[1]^2 + x[2]^2
+p_mean = sum(∫(p_tilde)dΩ) / sum(∫(1)dΩ)
+p_exact(x) = p_tilde(x) - p_mean
+
+f(x) = u_exact(x) - ∇(p_exact)(x)
+
+reffe_u = ReferenceFE(raviart_thomas, Float64, order; change_dof=false)
+reffe_p = ReferenceFE(lagrangian, Float64, order)
+
+V = FESpace(model, reffe_u; conformity=:HDiv, dirichlet_tags="boundary")
+Q = FESpace(model, reffe_p; conformity=:L2, constraint=:zeromean)
+U = TrialFESpace(V, u_exact)
+X = MultiFieldFESpace([U, Q])
+Y = MultiFieldFESpace([V, Q])
+
+a((u,p),(v,q)) = ∫( u⋅v + p*(∇⋅v) + (∇⋅u)*q )dΩ
+l((v,q))       = ∫( f⋅v )dΩ
+
+op = AffineFEOperator(a, l, X, Y)
+xh = solve(op)
+uh, ph = xh
+
+eu = uh - u_exact
+ep = ph - p_exact
+println("‖u - u_exact‖_L2 = ", sqrt(sum(∫(eu⋅eu)dΩ)))
+println("‖p - p_exact‖_L2 = ", sqrt(sum(∫(ep*ep)dΩ)))
diff --git a/test/view_permutations.jl b/test/view_permutations.jl
new file mode 100644
index 000000000..d682ce27f
--- /dev/null
+++ b/test/view_permutations.jl
@@ -0,0 +1,361 @@
+using Printf
+using Gridap
+using Gridap.ReferenceFEs
+using Gridap.Polynomials
+using Gridap.Polynomials: BernsteinBasisOnSimplex, BarycentricPΛBasis,
+                          BarycentricPmΛBasis, CartProdPolyBasis,
+                          CompWiseTensorPolyBasis, ModalC0
+using Gridap.Fields: evaluate
+using LinearAlgebra: norm, dot
+
+# ─────────────────────────────────────────────────────────────────────────────
+# Core function
+# ─────────────────────────────────────────────────────────────────────────────
+
+"""
+    evaluate_at_permutations(poly::Polytope, basis; quad_degree=nothing)
+        -> (vertex_perms, E_ref, E_perms)
+
+Evaluates `basis` at a fixed set of quadrature points on `poly`, then at the
+corresponding permuted points for each admissible vertex permutation of `poly`.
+
+Returns:
+- `vertex_perms`: the list from `get_vertex_permutations(poly)`
+- `E_ref`:   `(nq × nbasis)` matrix, `E_ref[q, i]  = φᵢ(xq[q])`
+- `E_perms`: one `(nq × nbasis)` matrix per permutation σ,
+             `E_perms[k][q, i] = φᵢ(Mσ(xq[q]))` where Mσ is the affine
+             self-map of `poly` induced by σ.
+
+If a basis function is equivariant, then column `i` of `E_perms[k]` equals
+column `π(i)` of `E_ref` (up to a scalar for signed/scaled cases).
+"""
+function evaluate_at_permutations(poly::Polytope, basis, conf=GradConformity(); quad_degree=nothing)
+  deg = isnothing(quad_degree) ? 2*get_order(basis) + 2 : quad_degree
+  quad = Quadrature(poly, deg)
+  xq   = get_coordinates(quad)
+
+  E_ref = evaluate(basis, xq)
+  verts  = get_vertex_coordinates(poly)
+  vertex_perms = get_vertex_permutations(poly)
+
+  Push = conf isa GradConformity ? nothing :
+         conf isa CurlConformity ? CoVariantPiolaMap :
+         conf isa DivConformity ? ContraVariantPiolaMap : nothing
+
+  geomap_basis = get_shapefuns(LagrangianRefFE(Float64,poly,1))
+
+  E_perms = map(vertex_perms) do σ
+    pverts = verts[invperm(σ)]
+    geo_map = linear_combination(pverts, geomap_basis)
+    pushed_basis = if isnothing(Push)
+      basis
+    else
+      basis
+      Jt = Broadcasting(∇)(geo_map)
+      lazy_map(Push(), Fill(basis,1), Fill(Jt,1))[1]
+    end
+    pushed_xq = evaluate(inverse_map(geo_map), xq)
+    evaluate(pushed_basis, pushed_xq)
+  end
+
+  vertex_perms, E_ref, E_perms
+end
+
+# ─────────────────────────────────────────────────────────────────────────────
+# Analysis helpers
+# ─────────────────────────────────────────────────────────────────────────────
+
+"""
+    analyse_column(col_σ, E_ref; tol=1e-10)
+
+For a single column `col_σ` of a permuted evaluation matrix, find which column
+`j` of `E_ref` it is most proportional to, the proportionality constant `c`,
+and the residual `‖col_σ − c·E_ref[:,j]‖`.
+
+Returns `(j, c, residual)`.
+"""
+# Works for both scalar (Float64) and vector-valued (VectorValue) columns.
+_col_dot(a, b) = sum(ai ⊙ bi for (ai, bi) in zip(a, b))
+_col_nrm(a)   = sqrt(real(_col_dot(a, a)))
+
+# Flatten a column to a plain Float64 vector (needed for the 2-col least squares).
+function _flatten(col)
+  T = eltype(col)
+  T <: Real && return Float64.(col)
+  # VectorValue or similar MultiValue: unpack each entry's components
+  Float64[xi for v in col for xi in Tuple(v)]
+end
+
+function analyse_column(col_σ, E_ref; tol=1e-10)
+  n = size(E_ref, 2)
+  best_j, best_c, best_res = 0, NaN, Inf
+  for j in 1:n
+    col_j = E_ref[:, j]
+    nrm   = _col_nrm(col_j)
+    if nrm < tol
+      if _col_nrm(col_σ) < tol
+        best_j, best_c, best_res = j, 1.0, 0.0
+        break
+      end
+      continue
+    end
+    c   = _col_dot(col_σ, col_j) / nrm^2
+    res = _col_nrm(col_σ .- c .* col_j)
+    if res < best_res
+      best_j, best_c, best_res = j, real(c), res
+    end
+  end
+  best_j, best_c, best_res
+end
+
+"""
+    analyse_column_2(col_σ, E_ref; tol=1e-10)
+
+Like `analyse_column` but searches over all pairs of columns in `E_ref`,
+solving the 2-column least-squares problem
+
+    min_{c1,c2} ‖col_σ − c1·col_ref[:,j1] − c2·col_ref[:,j2]‖
+
+Returns `(j1, j2, c1, c2, residual)` for the best-fitting pair.
+Works for both scalar and vector-valued (VectorValue) entries by flattening
+to plain Float64 vectors before solving.
+"""
+function analyse_column_2(col_σ, E_ref; tol=1e-10)
+  n  = size(E_ref, 2)
+  b  = _flatten(col_σ)
+  cols = [_flatten(E_ref[:, j]) for j in 1:n]
+
+  best_j1, best_j2, best_c1, best_c2, best_res = 0, 0, NaN, NaN, Inf
+  for j1 in 1:n, j2 in j1+1:n
+    A   = hcat(cols[j1], cols[j2])
+    x   = A \ b          # least-squares via QR
+    res = norm(b .- A * x)
+    if res < best_res
+      best_j1, best_j2, best_c1, best_c2, best_res = j1, j2, x[1], x[2], res
+    end
+  end
+  best_j1, best_j2, best_c1, best_c2, best_res
+end
+
+"""
+    classify_permutation(E_ref, E_σ; tol=1e-8) -> Symbol
+
+Classify the relationship between `E_ref` and `E_σ`:
+- `:pure`   — each column of E_σ equals a column of E_ref (pure permutation)
+- `:signed` — each column equals ±1 times a column of E_ref (signed permutation)
+- `:scaled` — each column is proportional to some column of E_ref (scaled permutation)
+- `:none`   — some column is not proportional to any column of E_ref
+"""
+function classify_permutation(E_ref, E_σ; tol=1e-8)
+  n = size(E_ref, 2)
+  all_pure, all_signed = true, true
+  scale = _col_nrm(E_ref[:, 1])
+  for i in 1:n
+    _, c, res = analyse_column(E_σ[:, i], E_ref; tol)
+    res > tol * (scale + tol)  && return :none
+    abs(abs(c) - 1.0) > tol   && (all_pure = false; all_signed = false)
+    abs(c - 1.0) > tol        && (all_pure = false)  # c = -1 is signed, not pure
+  end
+  all_pure ? :pure : all_signed ? :signed : :scaled
+end
+
+# ─────────────────────────────────────────────────────────────────────────────
+# Pretty-printing
+# ─────────────────────────────────────────────────────────────────────────────
+
+"""
+    print_permutation_analysis(poly, basis, label=""; tol=1e-8, show_matrices=false)
+
+Run `evaluate_at_permutations` and pretty-print a table showing, for each
+vertex permutation σ, how each basis function transforms:
+
+    φᵢ(Mσ(x)) ≈ c · φⱼ(x)
+
+If `show_matrices=true`, also prints the full `E_ref` and each `E_σ` matrix
+(rows = quadrature points, columns = basis functions) so you can inspect the
+raw evaluations directly.
+"""
+function print_permutation_analysis(poly::Polytope, basis, label="", conf=GradConformity();
+                                    tol=1e-8, show_matrices=false)
+
+  vertex_perms, E_ref, E_perms = evaluate_at_permutations(poly, basis, conf)
+  n = size(E_ref, 2)
+
+  lbl = isempty(label) ? "" : " [$label]"
+  println("━"^72)
+  println("  Polytope: $poly on $(typeof(basis).name.name) with $conf")
+  println("  n_basis=$(n), n_perms=$(length(vertex_perms))$lbl")
+  println("━"^72)
+
+  if show_matrices
+    println("  E_ref  (rows=quad pts, cols=basis funs):")
+    _print_matrix(E_ref)
+    println()
+  end
+
+  for (σ, E_σ) in zip(vertex_perms, E_perms)
+    kind = classify_permutation(E_ref, E_σ; tol)
+    kind_str = Dict(:pure=>"\e[32mpure perm\e[0m", :signed=>"\e[33msigned perm\e[0m",
+                    :scaled=>"\e[35mscaled perm\e[0m", :none=>"\e[31mNOT a perm\e[0m")[kind]
+    println("  σ = $σ  →  $kind_str")
+
+    if kind === :pure
+      perm = [analyse_column(E_σ[:, i], E_ref; tol)[1] for i in 1:n]
+      println("       π = $perm")
+    elseif kind === :signed
+      for i in 1:n
+        j, c, _ = analyse_column(E_σ[:, i], E_ref; tol)
+        sgn = c > 0 ? "+" : "-"
+        println("     φ[$i](Mσ·) = $sgn φ[$j]")
+      end
+    else
+      for i in 1:n
+        j, c, res = analyse_column(E_σ[:, i], E_ref; tol)
+        scale = _col_nrm(E_ref[:, 1])
+        bad   = res > tol * (scale + tol)
+
+        if !bad
+          # Single-column fit is fine
+          c_str = @sprintf("%.4f", c)
+          r_str = res < 1e-12 ? "≈0" : @sprintf("%.2e", res)
+          println("     φ[$i](Mσ·) ≈ $c_str · φ[$j]   (res=$r_str)")
+        else
+          # Try 2-column combination
+          j1, j2, c1, c2, res2 = analyse_column_2(E_σ[:, i], E_ref; tol)
+          good2 = n > 1 && res2 < tol * (scale + tol)
+          if good2
+            c1s = @sprintf("%.4f", c1)
+            c2s = @sprintf("%.4f", c2)
+            r2s = res2 < 1e-12 ? "≈0" : @sprintf("%.2e", res2)
+            println("     φ[$i](Mσ·) ≈ $c1s·φ[$j1] + $c2s·φ[$j2]   (res=$r2s)")
+          else
+            # Neither fit worked — show best 1-col with flag
+            c_str = @sprintf("%.4f", c)
+            r_str = @sprintf("%.2e", res)
+            println("  !  φ[$i](Mσ·) ≈ $c_str · φ[$j]   (1-col res=$r_str, 2-col res=$(@sprintf("%.2e",res2)))")
+          end
+        end
+      end
+    end
+
+    if show_matrices
+      println("  E_σ  (σ=$σ):")
+      _print_matrix(E_σ)
+    end
+    println()
+  end
+end
+
+# Format a single evaluation value for matrix display.
+_fmt(x::Real)                   = @sprintf("%8.4f", x)
+_fmt(x::VectorValue)            = "[" * join(_fmt.(Tuple(x)), ", ") * "]"
+_fmt(x)                         = string(x)
+
+function _print_matrix(E::AbstractMatrix; indent="    ")
+  nq, nb = size(E)
+  for q in 1:nq
+    row = join((_fmt(E[q, i]) for i in 1:nb), "  ")
+    println(indent, "q$q: ", row)
+  end
+end
+
+# ─────────────────────────────────────────────────────────────────────────────
+# Examples
+# ─────────────────────────────────────────────────────────────────────────────
+
+println("\n" * "═"^72)
+println("  BASIS PERMUTATION ANALYSIS")
+println("  Confirming pure / signed / scaled / none cases")
+println("═"^72 * "\n")
+
+# ── 1. BernsteinBasisOnSimplex on SEGMENT (pure perm) ─────────────────────
+print_permutation_analysis(SEGMENT,
+  BernsteinBasisOnSimplex{1}(Float64, 2),
+  "BernsteinBasisOnSimplex{1} order 2  [expect: pure]")
+
+# ── 2. BarycentricPΛBasis k=0 on SEGMENT (pure perm) ─────────────────────
+print_permutation_analysis(SEGMENT,
+  BarycentricPΛBasis{1}(Float64, 2, 0),
+  "BarycentricPΛBasis{1} k=0 order 2  [expect: pure]")
+
+# ── 3. BarycentricPmΛBasis k=0 on SEGMENT, order 2 (pure perm up to r=2) ─
+print_permutation_analysis(SEGMENT,
+  BarycentricPmΛBasis{1}(Float64, 2, 0),
+  "BarycentricPmΛBasis{1} k=0 order 2  [expect: pure for r≤2]")
+
+# ── 4. BarycentricPmΛBasis k=0 on SEGMENT, order 3 (scaled, c=2) ─────────
+print_permutation_analysis(SEGMENT,
+  BarycentricPmΛBasis{1}(Float64, 3, 0),
+  "BarycentricPmΛBasis{1} k=0 order 3  [expect: scaled c=2]")
+
+# ── 5. BarycentricPΛBasis k=1 on SEGMENT  ─────────────────────────────────
+#    (RT TRI cell moment  /  ND1 TRI cell moment)
+print_permutation_analysis(SEGMENT,
+  BarycentricPΛBasis{1}(Float64, 1, 1),
+  "BarycentricPΛBasis{1} k=1 order 1  [expect: signed]")
+
+# ── 6. BarycentricPΛBasis k=0 on TRI (RT TET facet, ND1 TET edge basis) ──
+print_permutation_analysis(TRI,
+  BarycentricPΛBasis{2}(Float64, 2, 0),
+  "BarycentricPΛBasis{2} k=0 order 2  [expect: pure]")
+
+# ── 7. BarycentricPΛBasis k=1 on TRI (ND1 TET triangular face moments) ───
+print_permutation_analysis(TRI,
+  BarycentricPΛBasis{2}(Float64, 1, 1),
+  "BarycentricPΛBasis{2} k=1 order 1  [expect: signed]", CurlConformity())
+
+# ── 8. CartProdPolyBasis{1,ModalC0} order 2 on SEGMENT (RT/ND QUAD edge) ──
+print_permutation_analysis(SEGMENT,
+  CartProdPolyBasis(ModalC0, Val(1), Float64, 2),
+  "CartProdPolyBasis{1,ModalC0} order 2  [expect: signed]")
+
+# ── 9. CartProdPolyBasis{1,Bernstein} order 2 on SEGMENT (pure perm) ──────
+print_permutation_analysis(SEGMENT,
+  CartProdPolyBasis(Bernstein, Val(1), Float64, 2),
+  "CartProdPolyBasis{1,Bernstein} order 2  [expect: pure]")
+
+# ── 10. CartProdPolyBasis{2,Bernstein} on QUAD (pure perm) ────────────────
+print_permutation_analysis(QUAD,
+  CartProdPolyBasis(Bernstein, Val(2), Float64, 2),
+  "CartProdPolyBasis{2,Bernstein} Q₂ on QUAD  [expect: pure]")
+
+bQ2Λ1_Q2 = get_shapefuns(ReferenceFE(QUAD, nedelec, 1; change_dof=false))
+print_permutation_analysis(QUAD,
+  bQ2Λ1_Q2,
+  "FEEC_poly_basis(Val(2), Float64, 2, 1, :Q⁻, Bernstein) Q₂ on QUAD  [expect: signed perm]", CurlConformity())
+print_permutation_analysis(QUAD,
+  bQ2Λ1_Q2,
+  "FEEC_poly_basis(Val(2), Float64, 2, 1, :Q⁻, Bernstein) Q₂ on QUAD  [expect: signed perm]", DivConformity())
+
+bQ2Λ1_Q2 = get_shapefuns(ReferenceFE(QUAD, nedelec, 0; change_dof=true))
+print_permutation_analysis(QUAD,
+  bQ2Λ1_Q2,
+  "FEEC_poly_basis(Val(2), Float64, 2, 1, :Q⁻, Bernstein) Q₂ on QUAD  [expect: signed perm]", CurlConformity())
+
+
+# Nédélec 1 TRI
+r = 4 # first non permutables are r = 2, 4, 5, ...  (I believe that changing the TRI bubble function could lead to permutable exept r = 3k-1)
+bPm2Λ1_T2 = get_shapefuns(ReferenceFE(TRI, nedelec, r-1; change_dof=false))
+print_permutation_analysis(TRI,
+  bPm2Λ1_T2,
+  "nedelec1 order $(r-1) Bernstein TRI  [expect: signed perm]", CurlConformity())
+
+r = 1 # non permutable if  r >= 2
+bPm2Λ1_T2 = FEEC_poly_basis(Val(2), Float64, r, 1, :P⁻, Bernstein)
+print_permutation_analysis(TRI,
+  bPm2Λ1_T2,
+  "FEEC_poly_basis(Val(2), Float64, $r, 1, :P⁻, Bernstein) Q₂ on TRI  [expect: signed perm]", CurlConformity())
+
+# Nédélec 2 TRI
+r = 2 # first non permutables are r = 3, 4, ...
+bP2Λ1_T2 = get_shapefuns(ReferenceFE(TRI, nedelec2, r; change_dof=false))
+print_permutation_analysis(TRI,
+  bP2Λ1_T2,
+  "nedelec2 order $(r) Bernstein TRI  [expect: signed perm]", CurlConformity())
+
+r = 2 # non permutable if  r = 0, 3, 4, ...
+bP2Λ1_T2 = FEEC_poly_basis(Val(2), Float64, r, 1, :P, Bernstein)
+print_permutation_analysis(TRI,
+  bP2Λ1_T2,
+  "FEEC_poly_basis(Val(2), Float64, $r, 1, :P, Bernstein) Q₂ on TRI  [expect: signed perm]", CurlConformity())
+