Implementing bdm

src/FESpaces/DivConformingFESpaces.jl

@@ -27,7 +27,7 @@
 struct TransformRTDofBasis{Dc,Dp} <: Map end ;
 
 function get_cell_dof_basis(model::DiscreteModel,
-                            cell_reffe::AbstractArray{<:GenericRefFE{RaviartThomas}},
+                            cell_reffe::AbstractArray{<:GenericRefFE{<:DivConforming}},
                             ::DivConformity,
                             sign_flip=get_sign_flip(model, cell_reffe))
     cell_map  = get_cell_map(Triangulation(model))
@@ -37,15 +37,16 @@ function get_cell_dof_basis(model::DiscreteModel,
     Jtx       = lazy_map(evaluate,Jt,x)
     reffe     = cell_reffe[1]
     Dc        = num_dims(reffe)
-    et        = return_type(get_prebasis(reffe))
+    # @santiagobadia: A hack here, for RT returns Float64 and for BDM VectorValue{Float64}
+    et        = eltype(return_type(get_prebasis(reffe)))
     pt        = Point{Dc,et}
     Dp        = first(size(return_type(phi,zero(pt))))
     k         = TransformRTDofBasis{Dc,Dp}()
     lazy_map(k,cell_reffe,Jtx,sign_flip)
 end
 
 function get_cell_shapefuns(model::DiscreteModel,
-                            cell_reffe::AbstractArray{<:GenericRefFE{RaviartThomas}},
+                            cell_reffe::AbstractArray{<:GenericRefFE{<:DivConforming}},
                             ::DivConformity,
                             sign_flip=get_sign_flip(model, cell_reffe))
     cell_reffe_shapefuns=lazy_map(get_shapefuns,cell_reffe)
@@ -116,20 +117,21 @@ function evaluate!(cache,k::SignFlipMap,reffe,cell_id)
 end
 
 function get_sign_flip(model::DiscreteModel,
-                       cell_reffe::AbstractArray{<:GenericRefFE{RaviartThomas}})
+                       cell_reffe::AbstractArray{<:GenericRefFE{<:DivConforming}})
     lazy_map(SignFlipMap(model),
             cell_reffe,
             IdentityVector(Int32(num_cells(model))))
 end
 
 function return_cache(::TransformRTDofBasis{Dc,Dp},
-                      reffe::GenericRefFE{RaviartThomas},
+                      reffe::GenericRefFE{<:DivConforming},
                       Jtx,
                       ::AbstractVector{Bool}) where {Dc,Dp}
   p = get_polytope(reffe)
   prebasis = get_prebasis(reffe)
   order = get_order(prebasis)
-  et = return_type(prebasis)
+  # @santiagobadia: Hack as above
+  et = eltype(return_type(prebasis))
   dofs = get_dof_basis(reffe)
   nodes, nf_nodes, nf_moments =  get_nodes(dofs),
                                  get_face_nodes_dofs(dofs),
@@ -143,7 +145,7 @@ end
 
 function evaluate!(cache,
                    ::TransformRTDofBasis,
-                   reffe::GenericRefFE{RaviartThomas},
+                   reffe::GenericRefFE{<:DivConforming},
                    Jt_q,
                    sign_flip::AbstractVector{Bool})
   nodes, nf_nodes, nf_moments, face_moments = cache

src/Polynomials/QCurlGradMonomialBases.jl

@@ -39,6 +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
 
 function return_cache(f::QCurlGradMonomialBasis,x::AbstractVector{<:Point})

src/ReferenceFEs/BDMRefFEs.jl

@@ -0,0 +1,192 @@
+struct BDM <: DivConforming end
+
+const bdm = BDM()
+
+"""
+BDMRefFE(::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 BDMRefFE(::Type{et},p::Polytope,order::Integer) where et
+
+  D = num_dims(p)
+
+  vet = VectorValue{num_dims(p),et}
+
+  if is_simplex(p)
+    prebasis = MonomialBasis(vet,p,order)
+  else
+    @notimplemented "BDM Reference FE only available for simplices"
+  end
+
+  nf_nodes, nf_moments = _BDM_nodes_and_moments(et,p,order,GenericField(identity))
+
+  face_own_dofs = _face_own_dofs_from_moments(nf_moments)
+
+  face_dofs = face_own_dofs
+
+  dof_basis = MomentBasedDofBasis(nf_nodes, nf_moments)
+
+  ndofs = num_dofs(dof_basis)
+
+  metadata = nothing
+
+  reffe = GenericRefFE{BDM}(
+  ndofs,
+  p,
+  prebasis,
+  dof_basis,
+  DivConformity(),
+  metadata,
+  face_dofs)
+
+  reffe
+end
+
+function ReferenceFE(p::Polytope,::BDM, order)
+  BDMRefFE(Float64,p,order)
+end
+
+function ReferenceFE(p::Polytope,::BDM,::Type{T}, order) where T
+  BDMRefFE(T,p,order)
+end
+
+function Conformity(reffe::GenericRefFE{BDM},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 BDM reference FE.
+
+    Possible values of conformity for this reference fe are $((:L2, hdiv...)).
+      """
+    end
+  end
+
+  function get_face_own_dofs(reffe::GenericRefFE{BDM}, conf::DivConformity)
+    get_face_dofs(reffe)
+  end
+
+  function _BDM_nodes_and_moments(::Type{et}, p::Polytope, order::Integer, phi::Field) where et
+
+    D = num_dims(p)
+    ft = VectorValue{D,et}
+    pt = Point{D,et}
+
+    nf_nodes = [ zeros(pt,0) for face in 1:num_faces(p)]
+    nf_moments = [ zeros(ft,0,0) for face in 1:num_faces(p)]
+
+    fcips, fmoments = _BDM_face_values(p,et,order,phi)
+    frange = get_dimrange(p,D-1)
+    nf_nodes[frange] = fcips
+    nf_moments[frange] = fmoments
+
+    if (order > 1)
+      ccips, cmoments = _BDM_cell_values(p,et,order,phi)
+      crange = get_dimrange(p,D)
+      nf_nodes[crange] = ccips
+      nf_moments[crange] = cmoments
+    end
+
+    nf_nodes, nf_moments
+  end
+
+  function _BDM_face_moments(p, fshfs, c_fips, fcips, fwips,phi)
+    nc = length(c_fips)
+    cfshfs = fill(fshfs, nc)
+    cvals = lazy_map(evaluate,cfshfs,c_fips)
+    cvals = [fwips[i].*cvals[i] for i in 1:nc]
+    fns = get_facet_normal(p)
+
+    # Must express the normal in terms of the real/reference system of
+    # coordinates (depending if phi≡I or phi is a mapping, resp.)
+    # Hence, J = transpose(grad(phi))
+
+    Jt = fill(∇(phi),nc)
+    Jt_inv = lazy_map(Operation(pinvJt),Jt)
+    det_Jt = lazy_map(Operation(meas),Jt)
+    change = lazy_map(*,det_Jt,Jt_inv)
+    change_ips = lazy_map(evaluate,change,fcips)
+
+    cvals = [ _broadcast(typeof(n),n,J.*b) for (n,b,J) in zip(fns,cvals,change_ips)]
+
+    return cvals
+  end
+
+  # It provides for every face the nodes and the moments arrays
+  function _BDM_face_values(p,et,order,phi)
+
+    # Reference facet
+    @assert is_simplex(p) "We are assuming that all n-faces of the same n-dim are the same."
+    fp = Polytope{num_dims(p)-1}(p,1)
+
+    # geomap from ref face to polytope faces
+    fgeomap = _ref_face_to_faces_geomap(p,fp)
+
+    # Nodes are integration points (for exact integration)
+    # Thus, we define the integration points in the reference
+    # face polytope (fips and wips). Next, we consider the
+    # n-face-wise arrays of nodes in fp (constant cell array c_fips)
+    # the one of the points in the polytope after applying the geopmap
+    # (fcips), and the weights for these nodes (fwips, a constant cell array)
+    # Nodes (fcips)
+    degree = (order)*2
+    fquad = Quadrature(fp,degree)
+    fips = get_coordinates(fquad)
+    wips = get_weights(fquad)
+
+    c_fips, fcips, fwips = _nfaces_evaluation_points_weights(p, fgeomap, fips, wips)
+
+    # Moments (fmoments)
+    # The BDM prebasis is expressed in terms of shape function
+    fshfs = MonomialBasis(et,fp,order)
+
+    # Face moments, i.e., M(Fi)_{ab} = q_RF^a(xgp_RFi^b) w_Fi^b n_Fi ⋅ ()
+    fmoments = _BDM_face_moments(p, fshfs, c_fips, fcips, fwips, phi)
+
+    return fcips, fmoments
+
+  end
+
+  function _BDM_cell_moments(p, cbasis, ccips, cwips)
+    # Interior DOFs-related basis evaluated at interior integration points
+    ishfs_iips = evaluate(cbasis,ccips)
+    return cwips.⋅ishfs_iips
+  end
+
+  # It provides for every cell the nodes and the moments arrays
+  function _BDM_cell_values(p,et,order,phi)
+    # Compute integration points at interior
+    degree = 2*(order)
+    iquad = Quadrature(p,degree)
+    ccips = get_coordinates(iquad)
+    cwips = get_weights(iquad)
+
+    # Cell moments, i.e., M(C)_{ab} = q_C^a(xgp_C^b) w_C^b ⋅ ()
+    if is_simplex(p)
+      T = VectorValue{num_dims(p),et}
+      # cbasis = GradMonomialBasis{num_dims(p)}(T,order-1)
+      cbasis = Polynomials.NedelecPrebasisOnSimplex{num_dims(p)}(order-2)
+    else
+      @notimplemented
+    end
+    cell_moments = _BDM_cell_moments(p, cbasis, ccips, cwips )
+
+    # Must scale weights using phi map to get the correct integrals
+    # scaling = meas(grad(phi))
+    Jt = ∇(phi)
+    Jt_inv = pinvJt(Jt)
+    det_Jt = meas(Jt)
+    change = det_Jt*Jt_inv
+    change_ips = evaluate(change,ccips)
+
+    cmoments = change_ips.⋅cell_moments
+
+    return [ccips], [cmoments]
+
+  end

src/ReferenceFEs/NedelecRefFEs.jl

@@ -126,7 +126,7 @@ function _Nedelec_edge_values(p,et,order)
   egeomap = _ref_face_to_faces_geomap(p,ep)
 
   # Compute integration points at all polynomial edges
-  degree = (order+1)*2
+  degree = (order)*2
   equad = Quadrature(ep,degree)
   cips = get_coordinates(equad)
   wips = get_weights(equad)
@@ -163,7 +163,7 @@ function _Nedelec_face_values(p,et,order)
   fgeomap = _ref_face_to_faces_geomap(p,fp)
 
   # Compute integration points at all polynomial edges
-  degree = (order+1)*2
+  degree = (order)*2
   fquad = Quadrature(fp,degree)
   fips = get_coordinates(fquad)
   wips = get_weights(fquad)
@@ -209,7 +209,7 @@ function _Nedelec_face_values_simplex(p,et,order)
   fgeomap = _ref_face_to_faces_geomap(p,fp)
 
   # Compute integration points at all polynomial edges
-  degree = (order+1)*2
+  degree = (order)*2
   fquad = Quadrature(fp,degree)
   fips = get_coordinates(fquad)
   wips = get_weights(fquad)
@@ -256,7 +256,7 @@ end
 function _Nedelec_cell_values(p,et,order)
 
   # Compute integration points at interior
-  degree = 2*(order+1)
+  degree = 2*(order)
   iquad = Quadrature(p,degree)
   ccips = get_coordinates(iquad)
   cwips = get_weights(iquad)

src/ReferenceFEs/RaviartThomasRefFEs.jl

@@ -1,6 +1,8 @@
 struct DivConformity <: Conformity end
 
-struct RaviartThomas <: ReferenceFEName end
+abstract type DivConforming <: ReferenceFEName end
+
+struct RaviartThomas <: DivConforming end
 
 const raviart_thomas = RaviartThomas()
 
@@ -176,7 +178,7 @@ function _RT_face_values(p,et,order,phi)
   # the one of the points in the polytope after applying the geopmap
   # (fcips), and the weights for these nodes (fwips, a constant cell array)
   # Nodes (fcips)
-  degree = (order+1)*2
+  degree = (order)*2
   fquad = Quadrature(fp,degree)
   fips = get_coordinates(fquad)
   wips = get_weights(fquad)
@@ -205,7 +207,7 @@ _p_filter(e,order) = (sum(e) <= order)
 # It provides for every cell the nodes and the moments arrays
 function _RT_cell_values(p,et,order,phi)
   # Compute integration points at interior
-  degree = 2*(order+1)
+  degree = 2*(order)
   iquad = Quadrature(p,degree)
   ccips = get_coordinates(iquad)
   cwips = get_weights(iquad)

src/ReferenceFEs/ReferenceFEs.jl

@@ -175,18 +175,22 @@ export CDConformity
 
 export SerendipityRefFE
 export RaviartThomasRefFE
+export BDMRefFE
 export NedelecRefFE
 export BezierRefFE
 export ModalC0RefFE
 
 export Lagrangian
+export DivConforming
 export RaviartThomas
+export BDM
 export Nedelec
 export Bezier
 export ModalC0
 
 export lagrangian
 export raviart_thomas
+export bdm
 export nedelec
 export bezier
 export modalC0
@@ -235,6 +239,8 @@ include("StrangQuadratures.jl")
 
 include("RaviartThomasRefFEs.jl")
 
+include("BDMRefFEs.jl")
+
 include("NedelecRefFEs.jl")
 
 include("MockDofs.jl")

test/FESpacesTests/CellFieldsTests.jl

@@ -178,5 +178,30 @@ end
   end
 end
 
+p = TRI
+D = num_dims(TRI)
+et = Float64
+source_model = CartesianDiscreteModel((0,1,0,1),(2,2)) |> simplexify
+
+
+@testset "Test interpolation BDM" begin
+  # BDM Space -> BDM Space
+
+  f(x) = VectorValue([x[1], x[2]])
+  reffe = BDMRefFE(et, p, 1)
+  V₁ = FESpace(source_model, reffe, conformity=:HDiv)
+  fh = interpolate_everywhere(f, V₁);
+  # Target RT Space
+  reffe = RaviartThomasRefFE(et, p, 1)
+  model = CartesianDiscreteModel((0,1,0,1),(40,40)) |> simplexify
+  V₂ = FESpace(model, reffe, conformity=:HDiv)
+
+  ifh = Interpolable(fh)
+  gh = interpolate_everywhere(ifh, V₂)
+  pts = [VectorValue(rand(2)) for i=1:10]
+  for pt in pts
+    @test gh(pt) ≈ fh(pt)
+  end
+end
 
 end # module