Skip to content

Commit

Permalink
Ported Nedelec and BDM
Browse files Browse the repository at this point in the history
  • Loading branch information
@JordiManyer
JordiManyer committed Nov 16, 2024
1 parent e0b4a77 commit dfdb7e4
Show file tree
Hide file tree
Showing 7 changed files with 631 additions and 494 deletions.
167 changes: 20 additions & 147 deletions src/ReferenceFEs/BDMRefFEs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,39 +10,31 @@ is in the P space of degree `order-1`.
"""
function BDMRefFE(::Type{et},p::Polytope,order::Integer) where et

D = num_dims(p)

vet = VectorValue{num_dims(p),et}
vet = VectorValue{D,et}

if is_simplex(p)
prebasis = MonomialBasis(vet,p,order)
fb = MonomialBasis{D-1}(et,order,Polynomials._p_filter)
cb = Polynomials.NedelecPrebasisOnSimplex{D}(order-2)
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)
function cmom(φ,μ,ds) # Cell moment function: σ_K(φ,μ) = ∫(φ·μ)dK
Broadcasting(Operation(⋅))(φ,μ)
end
function fmom(φ,μ,ds) # Face moment function : σ_F(φ,μ) = ∫((φ·n)*μ)dF
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
]

reffe
return MomentBasedReferenceFE(BDM(),p,prebasis,moments,DivConformity())
end

function ReferenceFE(p::Polytope,::BDM, order)
Expand All @@ -65,128 +57,9 @@ function Conformity(reffe::GenericRefFE{BDM},sym::Symbol)
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
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
@check 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
function get_face_own_dofs(reffe::GenericRefFE{BDM}, conf::DivConformity)
get_face_dofs(reffe)
end
83 changes: 0 additions & 83 deletions src/ReferenceFEs/Dofs.jl
View file
Original file line number Diff line number Diff line change
@@ -1,74 +1,5 @@
abstract type Dof <: Map end

# """
# abstract type Dof <: Map

# Abstract type representing a degree of freedom (DOF), a basis of DOFs, and related objects.
# These different cases are distinguished by the return type obtained when evaluating the `Dof`
# object on a `Field` object. See function [`evaluate_dof!`](@ref) for more details.

# The following functions needs to be overloaded

# - [`dof_cache`](@ref)
# - [`evaluate_dof!`](@ref)

# The following functions can be overloaded optionally

# - [`dof_return_type`](@ref)

# The interface is tested with

# - [`test_dof`](@ref)

# In most of the cases it is not strictly needed that types that implement this interface
# inherit from `Dof`. However, we recommend to inherit from `Dof`, when possible.


# """
# abstract type Dof <: Map end

# """
# return_cache(dof,field)

# Returns the cache needed to call `evaluate_dof!(cache,dof,field)`
# """
# function return_cache(dof::Dof,field)
# @abstractmethod
# end

# """
# evaluate_dof!(cache,dof,field)

# Evaluates the dof `dof` with the field `field`. It can return either an scalar value or
# an array of scalar values depending the case. The `cache` object is computed with function
# [`dof_cache`](@ref).

# When a mathematical dof is evaluated on a physical field, a scalar number is returned. If either
# the `Dof` object is a basis of DOFs, or the `Field` object is a basis of fields,
# or both objects are bases, then the returned object is an array of scalar numbers. The first
# dimensions in the resulting array are for the `Dof` object and the last ones for the `Field`
# object. E.g, a basis of `nd` DOFs evaluated at physical field returns a vector of `nd` entries.
# A basis of `nd` DOFs evaluated at a basis of `nf` fields returns a matrix of size `(nd,nf)`.
# """
# function evaluate!(cache,dof::Dof,field)
# @abstractmethod
# end

# """
# dof_return_type(dof,field)

# Returns the type for the value obtained with evaluating `dof` with `field`.

# It defaults to

# typeof(evaluate_dof(dof,field))
# """
# function return_type(dof::Dof,field)
# typeof(evaluate(dof,field))
# end

# Testers

"""
test_dof(dof,field,v;cmp::Function=(==))
Expand All @@ -93,17 +24,3 @@ function _test_dof(dof,field,v,cmp)
@test cmp(r,v)
@test typeof(r) == return_type(dof,field)
end

#struct DofEval <: Map end
#
#function return_cache(k::DofEval,dof,field)
# return_cache(dof,field)
#end
#
#function evaluate!(cache,k::DofEval,dof,field)
# evaluate!(cache,dof,field)
#end
#
#function return_type(k::DofEval,dof,field)
# return_type(dof,field)
#end
2 changes: 1 addition & 1 deletion src/ReferenceFEs/MomentBasedReferenceFEs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,7 +213,7 @@ end

"""
A moment is given by a triplet (f,σ,μ) where
- f is vector of ids if faces Fk
- f is vector of ids of faces Fk
- σ is a function σ(φ,μ,ds) that returns a Field-like object to be integrated over each Fk
- μ is a polynomials basis on Fk
Expand Down
Loading

0 comments on commit dfdb7e4

Please sign in to comment.