Replicated Raviart-Thomas refes with new machinery

src/Fields/FieldArrays.jl

@@ -509,6 +509,22 @@ for T in (:(Point),:(AbstractArray{<:Point}))
       evaluate!(r,bm,rs...)
     end
 
+    function return_cache(k::BroadcastOpFieldArray{typeof(∘)},x::$T)
+      f, g = k.args
+      cg = return_cache(g,x)
+      gx = evaluate!(cg,g,x)
+      cf = return_cache(f,gx)
+      return cg, cf
+    end
+
+    function evaluate!(cache, k::BroadcastOpFieldArray{typeof(∘)},x::$T)
+      cg, cf = cache
+      f, g = k.args
+      gx = evaluate!(cg,g,x)
+      fgx = evaluate!(cf,f,gx)
+      return fgx
+    end
+
   end
 end
 

src/Polynomials/PCurlGradJacobiPolynomialBases.jl

@@ -51,7 +51,7 @@ end
 
 get_order(f::PCurlGradJacobiPolynomialBasis{D,T}) where {D,T} = f.order
 
-return_type(::PCurlGradJacobiPolynomialBasis{D,T}) where {D,T} = T
+return_type(::PCurlGradJacobiPolynomialBasis{D,T}) where {D,T} = VectorValue{D,T}
 
 function return_cache(f::PCurlGradJacobiPolynomialBasis{D,T},x::AbstractVector{<:Point}) where {D,T}
   @check D == length(eltype(x)) "Incorrect number of point components"

src/Polynomials/PCurlGradMonomialBases.jl

@@ -52,7 +52,7 @@ end
 
 get_order(f::PCurlGradMonomialBasis{D,T}) where {D,T} = f.order
 
-return_type(::PCurlGradMonomialBasis{D,T}) where {D,T} = T
+return_type(::PCurlGradMonomialBasis{D,T}) where {D,T} = VectorValue{D,T}
 
 function return_cache(f::PCurlGradMonomialBasis{D,T},x::AbstractVector{<:Point}) where {D,T}
   @check D == length(eltype(x)) "Incorrect number of point components"

src/Polynomials/QCurlGradMonomialBases.jl

@@ -39,8 +39,7 @@ function QCurlGradMonomialBasis{D}(::Type{T},order::Int) where {D,T}
   QCurlGradMonomialBasis(T,order,terms,perms)
 end
 
-# @santiagobadia: This is dirty, I would put here VectorValue{D,T}
-return_type(::QCurlGradMonomialBasis{D,T}) where {D,T} = T
+return_type(::QCurlGradMonomialBasis{D,T}) where {D,T} = VectorValue{D,T}
 
 function return_cache(f::QCurlGradMonomialBasis,x::AbstractVector{<:Point})
   return_cache(f.qgrad,x)

src/Polynomials/QGradJacobiPolynomialBases.jl

@@ -46,6 +46,8 @@ num_terms(f::QGradJacobiPolynomialBasis{D,T}) where {D,T} = length(f.terms)*D
 
 get_order(f::QGradJacobiPolynomialBasis) = f.order
 
+return_type(::QGradJacobiPolynomialBasis{D,T}) where {D,T} = VectorValue{D,T}
+
 function return_cache(f::QGradJacobiPolynomialBasis{D,T},x::AbstractVector{<:Point}) where {D,T}
   @check D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
@@ -266,7 +268,7 @@ function QCurlGradJacobiPolynomialBasis{D}(::Type{T},order::Int) where {D,T}
   QCurlGradJacobiPolynomialBasis(T,order,terms,perms)
 end
 
-return_type(::QCurlGradJacobiPolynomialBasis{D,T}) where {D,T} = T
+return_type(::QCurlGradJacobiPolynomialBasis{D,T}) where {D,T} = VectorValue{D,T}
 
 function return_cache(f::QCurlGradJacobiPolynomialBasis,x::AbstractVector{<:Point})
   return_cache(f.qgrad,x)

src/Polynomials/QGradMonomialBases.jl

@@ -47,6 +47,8 @@ num_terms(f::QGradMonomialBasis{D,T}) where {D,T} = length(f.terms)*D
 
 get_order(f::QGradMonomialBasis) = f.order
 
+return_type(::QGradMonomialBasis{D,T}) where {D,T} = VectorValue{D,T}
+
 function return_cache(f::QGradMonomialBasis{D,T},x::AbstractVector{<:Point}) where {D,T}
   @check D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)

src/ReferenceFEs/DivConformingReferenceFEs.jl

@@ -0,0 +1,143 @@
+
+struct DivConformity <: Conformity end
+abstract type DivConforming <: ReferenceFEName end
+
+# RaviartThomas
+
+struct RaviartThomas <: DivConforming end
+const raviart_thomas = RaviartThomas()
+
+"""
+    RaviartThomasRefFE(::Type{et},p::Polytope,order::Integer) where et
+
+The `order` argument has the following meaning: the divergence of the  functions in this basis
+is in the Q space of degree `order`.
+"""
+function RaviartThomasRefFE(
+  ::Type{T},p::Polytope{D},order::Integer
+) where {T,D}
+
+  if is_n_cube(p)
+    prebasis = QCurlGradJacobiPolynomialBasis{D}(T,order)        # Prebasis
+    cb = QGradJacobiPolynomialBasis{D}(T,order-1)                # Cell basis
+    fb = JacobiPolynomialBasis{D}(T,order,Polynomials._q_filter) # Face basis
+  elseif is_simplex(p)
+    prebasis = PCurlGradMonomialBasis{D}(et,order)               # Prebasis
+    cb = MonomialBasis{D}(T,order-1,Polynomials._p_filter)       # Cell basis
+    fb = JacobiPolynomialBasis{D}(T,order,Polynomials._p_filter) # Face basis
+  else
+    @notimplemented "H(div) Reference FE only available for cubes and simplices"
+  end
+
+  function cmom(φ,μ,ds) # Cell moment function
+    Broadcasting(Operation(⋅))(φ,μ)
+  end
+  function fmom(φ,μ,ds) # Face moment function
+    n = get_normal(ds)
+    φn = Broadcasting(Operation(⋅))(φ,n)
+    Broadcasting(Operation(*))(φn,μ)
+  end
+  moments = [
+    (get_dimrange(p,D-1),fmom,fb), # Face moments
+    (get_dimrange(p,D),cmom,cb)    # Cell moments
+  ]
+
+  return MomentBasedReferenceFE(RaviartThomas(),p,prebasis,moments,DivConformity())
+end
+
+function ReferenceFE(p::Polytope,::RaviartThomas,order;kwargs...)
+  RaviartThomasRefFE(Float64,p,order;kwargs...)
+end
+
+function ReferenceFE(p::Polytope,::RaviartThomas,::Type{T},order;kwargs...) where T
+  RaviartThomasRefFE(T,p,order;kwargs...)
+end
+
+function Conformity(reffe::GenericRefFE{RaviartThomas},sym::Symbol)
+  hdiv = (:Hdiv,:HDiv)
+  if sym == :L2
+    L2Conformity()
+  elseif sym in hdiv
+    DivConformity()
+  else
+    @unreachable """\n
+    It is not possible to use conformity = $sym on a Raviart Thomas reference FE.
+
+    Possible values of conformity for this reference fe are $((:L2, hdiv...)).
+    """
+  end
+end
+
+function get_face_own_dofs(reffe::GenericRefFE{RaviartThomas}, conf::DivConformity)
+  get_face_dofs(reffe)
+end
+
+# TODO: Please remove me
+function JacobiBasis(::Type{T},p::Polytope,orders) where T
+  compute_jacobi_basis(T,p,orders)
+end
+function JacobiBasis(::Type{T},p::Polytope{D},order::Int) where {D,T}
+  orders = tfill(order,Val{D}())
+  JacobiBasis(T,p,orders)
+end
+function compute_jacobi_basis(::Type{T},p::ExtrusionPolytope{D},orders) where {D,T}
+  extrusion = Tuple(p.extrusion)
+  terms = _monomial_terms(extrusion,orders)
+  JacobiPolynomialBasis{D}(T,orders,terms)
+end
+
+# ContraVariantPiolaMap
+
+struct ContraVariantPiolaMap <: Map end
+
+function evaluate!(
+  cache,
+  ::Broadcasting{typeof(∇)},
+  a::Fields.BroadcastOpFieldArray{ContraVariantPiolaMap}
+)
+  v, Jt, sign_flip = a.args
+  ∇v = Broadcasting(∇)(v)
+  k = ContraVariantPiolaMap()
+  Broadcasting(Operation(k))(∇v,Jt,sign_flip)
+end
+
+function lazy_map(
+  ::Broadcasting{typeof(gradient)},
+  a::LazyArray{<:Fill{Broadcasting{Operation{ContraVariantPiolaMap}}}}
+)
+  v, Jt, sign_flip = a.args
+  ∇v = lazy_map(Broadcasting(∇),v)
+  k = ContraVariantPiolaMap()
+  lazy_map(Broadcasting(Operation(k)),∇v,Jt,sign_flip)
+end
+
+function lazy_map(
+  k::ContraVariantPiolaMap,
+  cell_ref_shapefuns::AbstractArray{<:AbstractArray{<:Field}},
+  cell_map::AbstractArray{<:Field},
+  sign_flip::AbstractArray{<:AbstractArray{<:Field}}
+)
+  cell_Jt = lazy_map(∇,cell_map)
+  lazy_map(Broadcasting(Operation(k)),cell_ref_shapefuns,cell_Jt,sign_flip)
+end
+
+function evaluate!(
+  cache,::ContraVariantPiolaMap,
+  v::Number,
+  Jt::Number,
+  sign_flip::Bool
+)
+  idetJ = 1/meas(Jt)
+  ((-1)^sign_flip*v)⋅(idetJ*Jt)
+end
+
+function evaluate!(
+  cache,
+  k::ContraVariantPiolaMap,
+  v::AbstractVector{<:Field},
+  phi::Field,
+  sign_flip::AbstractVector{<:Field}
+)
+  Jt = ∇(phi)
+  Broadcasting(Operation(k))(v,Jt,sign_flip)
+end