Skip to content

Commit

Permalink
Merge pull request #823 from gridap/implementing_bdm
Browse files Browse the repository at this point in the history
Implementing bdm
  • Loading branch information
santiagobadia authored Aug 26, 2022
2 parents 0dd61fd + 0af90b3 commit f4b4b72
Show file tree
Hide file tree
Showing 10 changed files with 479 additions and 34 deletions.
16 changes: 9 additions & 7 deletions src/FESpaces/DivConformingFESpaces.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand All @@ -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)
Expand Down Expand Up @@ -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),
Expand All @@ -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
Expand Down
1 change: 1 addition & 0 deletions src/Polynomials/QCurlGradMonomialBases.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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})
Expand Down
192 changes: 192 additions & 0 deletions src/ReferenceFEs/BDMRefFEs.jl
Original file line number Diff line number Diff line change
@@ -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
8 changes: 4 additions & 4 deletions src/ReferenceFEs/NedelecRefFEs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
8 changes: 5 additions & 3 deletions src/ReferenceFEs/RaviartThomasRefFEs.jl
Original file line number Diff line number Diff line change
@@ -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()

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
6 changes: 6 additions & 0 deletions src/ReferenceFEs/ReferenceFEs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -235,6 +239,8 @@ include("StrangQuadratures.jl")

include("RaviartThomasRefFEs.jl")

include("BDMRefFEs.jl")

include("NedelecRefFEs.jl")

include("MockDofs.jl")
Expand Down
25 changes: 25 additions & 0 deletions test/FESpacesTests/CellFieldsTests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Loading

0 comments on commit f4b4b72

Please sign in to comment.