Dof permutations for moment-based reffes

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.
 

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
 
 """

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
 )

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 #
 #################

src/Polynomials/BarycentricPΛBases.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

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,25 +51,25 @@ 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
 
 
   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

src/ReferenceFEs/BDMRefFEs.jl

@@ -1,7 +1,7 @@
 """
     struct BDM <: ReferenceFEName
 """
-struct BDM <: ReferenceFEName end
+struct BDM <: MomentBasedRefFEName end
 
 """
     const bdm = BDM()

src/ReferenceFEs/CrouzeixRaviartRefFEs.jl

@@ -1,7 +1,7 @@
 """
     struct CrouzeixRaviart  <: ReferenceFEName
 """
-struct CrouzeixRaviart <: ReferenceFEName end
+struct CrouzeixRaviart <: MomentBasedRefFEName end
 
 """
     const crouzeix_raviart = CrouzeixRaviart()