From c5be3d433080083cc173ac47eedd128a79b9fd7f Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Thu, 12 Dec 2024 15:40:41 +1100
Subject: [PATCH] Fixed Nedelec

---
 src/ReferenceFEs/MomentBasedReferenceFEs.jl   | 34 ++++++++++++++-----
 src/ReferenceFEs/NedelecRefFEs.jl             | 10 ++++--
 .../CurlConformingFESpacesTests.jl            |  4 +--
 3 files changed, 35 insertions(+), 13 deletions(-)

diff --git a/src/ReferenceFEs/MomentBasedReferenceFEs.jl b/src/ReferenceFEs/MomentBasedReferenceFEs.jl
index 4632b3ad3..27e551928 100644
--- a/src/ReferenceFEs/MomentBasedReferenceFEs.jl
+++ b/src/ReferenceFEs/MomentBasedReferenceFEs.jl
@@ -52,7 +52,7 @@ function num_dofs(b::MomentBasedDofBasis)
   n
 end
 
-function return_cache(b::MomentBasedDofBasis{P,V},field) where {P,V}
+function return_cache(b::MomentBasedDofBasis{P,V}, field) where {P,V}
   alloc_cache(vals::AbstractVector,T,ndofs) = zeros(T,ndofs)
   alloc_cache(vals::AbstractMatrix,T,ndofs) = zeros(T,ndofs,size(vals,2))
   cf = return_cache(field,b.nodes)
@@ -60,18 +60,15 @@ function return_cache(b::MomentBasedDofBasis{P,V},field) where {P,V}
   T = typeof(dot(zero(V),zero(eltype(vals))))
   r = alloc_cache(vals,T,num_dofs(b))
   c = CachedArray(r)
-  (c, cf)
+  return c, cf
 end
 
-function evaluate!(cache,b::MomentBasedDofBasis,field)
+function evaluate!(cache, b::MomentBasedDofBasis, field::Field)
   c, cf = cache
+  setsize!(c, size(b))
   vals = evaluate!(cf,field,b.nodes)
   dofs = c.array
-  _eval_moment_dof_basis!(dofs,vals,b)
-  dofs
-end
 
-function _eval_moment_dof_basis!(dofs,vals::AbstractVector,b)
   o = 1
   z = zero(eltype(dofs))
   face_nodes = b.face_nodes
@@ -90,9 +87,16 @@ function _eval_moment_dof_basis!(dofs,vals::AbstractVector,b)
       end
     end
   end
+
+  return dofs
 end
 
-function _eval_moment_dof_basis!(dofs,vals::AbstractMatrix,b)
+function evaluate!(cache, b::MomentBasedDofBasis, field::AbstractVector{<:Field})
+  c, cf = cache
+  setsize!(c, (size(b,1),length(field)))
+  vals = evaluate!(cf,field,b.nodes)
+  dofs = c.array
+
   o = 1
   na = size(vals,2)
   z = zero(eltype(dofs))
@@ -114,6 +118,8 @@ function _eval_moment_dof_basis!(dofs,vals::AbstractMatrix,b)
       end
     end
   end
+
+  return dofs
 end
 
 # MomentBasedReferenceFE
@@ -160,11 +166,21 @@ function get_edge_tangent(m::FaceMeasure{1,Dc}) where {Dc}
   return ConstantField(t[m.face])
 end
 
+# Extends a Df-dimensional vector to a Dc-dimensional one that 
+# lives in the tangent space of the Dc-embedded Df-dimensional manifold.
+# Equivalent to transpose(∇(fmap))
 function get_extension(m::FaceMeasure{Df,Dc}) where {Df,Dc}
   @assert Df == Dc - 1
   vs = ReferenceFEs._nfaces_vertices(Float64,m.cpoly,Df)[m.face]
-  return ConstantField(TensorValue(hcat([vs[2]-vs[1]...],[vs[3]-vs[1]...])))
+  J = TensorValue(hcat([vs[2]-vs[1]...],[vs[3]-vs[1]...]))
+  return ConstantField(J/meas(transpose(J)))
 end
+# function get_extension(m::FaceMeasure{Df,Dc}) where {Df,Dc}
+#   @assert Df == Dc - 1
+#   fmap = m.fmaps[m.face]
+#   J = Broadcasting(∇)(fmap)
+#   return Operation(*)(Operation(transpose)(J),Operation(x -> 1/meas(x))(J))
+# end
 
 # TO DO: Bug in 3D, when n==2; n==4 and D==3. Also, working on this to make better and more general.
 function get_facet_measure(p::Polytope{D}, face::Int) where D
diff --git a/src/ReferenceFEs/NedelecRefFEs.jl b/src/ReferenceFEs/NedelecRefFEs.jl
index 7cdbb98e5..3faa07361 100644
--- a/src/ReferenceFEs/NedelecRefFEs.jl
+++ b/src/ReferenceFEs/NedelecRefFEs.jl
@@ -30,7 +30,7 @@ function NedelecRefFE(::Type{et},p::Polytope,order::Integer) where et
   function cmom(φ,μ,ds) # Cell moment function: σ_K(φ,μ) = ∫(φ⋅μ)dK
     Broadcasting(Operation(⋅))(φ,μ)
   end
-  function fmom(φ,μ,ds) # Face moment function: σ_F(φ,μ) = ∫((φ×n)⋅μ)dF
+  function fmom_HEX(φ,μ,ds) # Face moment function: σ_F(φ,μ) = ∫((φ×n)⋅μ)dF
     o = get_facet_orientations(ds.cpoly)[ds.face] # This is a hack to avoid a sign map
     n = o*get_facet_normal(ds)
     E = get_extension(ds)
@@ -38,6 +38,11 @@ function NedelecRefFE(::Type{et},p::Polytope,order::Integer) where et
     φn = Broadcasting(Operation(×))(n,φ)
     Broadcasting(Operation(⋅))(φn,Eμ)
   end
+  function fmom_TET(φ,μ,ds) # Face moment function: σ_F(φ,μ) = ∫((φ×n)⋅μ)dF
+    E = get_extension(ds)
+    Eμ = Broadcasting(Operation(⋅))(E,μ) # We have to extend the basis to 3D
+    Broadcasting(Operation(⋅))(φ,Eμ)
+  end
   function emom(φ,μ,ds) # Edge moment function: σ_E(φ,μ) = ∫((φ⋅t)*μ)dE
     t = get_edge_tangent(ds)
     φt = Broadcasting(Operation(⋅))(φ,t)
@@ -47,7 +52,8 @@ function NedelecRefFE(::Type{et},p::Polytope,order::Integer) where et
   moments = Tuple[
     (get_dimrange(p,1),emom,eb), # Edge moments
   ]
-  if D == 3 && order > 0
+  if (D == 3) && order > 0
+    fmom = ifelse(is_n_cube(p),fmom_HEX,fmom_TET)
     push!(moments,(get_dimrange(p,D-1),fmom,fb)) # Face moments
   end
   if (is_n_cube(p) && order > 0) || (is_simplex(p) && order > D-2)
diff --git a/test/FESpacesTests/CurlConformingFESpacesTests.jl b/test/FESpacesTests/CurlConformingFESpacesTests.jl
index bd46c3383..d5e983822 100644
--- a/test/FESpacesTests/CurlConformingFESpacesTests.jl
+++ b/test/FESpacesTests/CurlConformingFESpacesTests.jl
@@ -101,7 +101,7 @@ e = u - uh
 el2 = sqrt(sum( ∫( e⋅e )*dΩ ))
 @test el2 < 1.0e-10
 
-domain =(0,1,0,1,0,1)
+domain = (0,1,0,1,0,1)
 partition = (3,3,3)
 model = CartesianDiscreteModel(domain,partition)
 
@@ -128,7 +128,7 @@ e = u - uh
 dΩ = Measure(Ω,order)
 
 el2 = sqrt(sum( ∫( e⋅e )*dΩ ))
-@test el2 < 1.0e-10
+@test el2 < 2.0e-10
 
 
 # using Gridap.Visualization