diff --git a/src/FESpaces/DivConformingFESpaces.jl b/src/FESpaces/DivConformingFESpaces.jl index 85697b73b..5decd79ab 100644 --- a/src/FESpaces/DivConformingFESpaces.jl +++ b/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,7 +37,8 @@ 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}() @@ -45,7 +46,7 @@ function get_cell_dof_basis(model::DiscreteModel, 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 diff --git a/src/Polynomials/QCurlGradMonomialBases.jl b/src/Polynomials/QCurlGradMonomialBases.jl index 21f456824..17e5f4850 100644 --- a/src/Polynomials/QCurlGradMonomialBases.jl +++ b/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}) diff --git a/src/ReferenceFEs/BDMRefFEs.jl b/src/ReferenceFEs/BDMRefFEs.jl new file mode 100644 index 000000000..9e778e64f --- /dev/null +++ b/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 diff --git a/src/ReferenceFEs/NedelecRefFEs.jl b/src/ReferenceFEs/NedelecRefFEs.jl index 217b9f8ba..d9b3fa070 100644 --- a/src/ReferenceFEs/NedelecRefFEs.jl +++ b/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) diff --git a/src/ReferenceFEs/RaviartThomasRefFEs.jl b/src/ReferenceFEs/RaviartThomasRefFEs.jl index c29471652..d4388ea94 100644 --- a/src/ReferenceFEs/RaviartThomasRefFEs.jl +++ b/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) diff --git a/src/ReferenceFEs/ReferenceFEs.jl b/src/ReferenceFEs/ReferenceFEs.jl index 502b395ad..f7cd472ca 100644 --- a/src/ReferenceFEs/ReferenceFEs.jl +++ b/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") diff --git a/test/FESpacesTests/CellFieldsTests.jl b/test/FESpacesTests/CellFieldsTests.jl index abbbf12f6..ea9811c6d 100644 --- a/test/FESpacesTests/CellFieldsTests.jl +++ b/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 diff --git a/test/FESpacesTests/DivConformingFESpacesTests.jl b/test/FESpacesTests/DivConformingFESpacesTests.jl index b8c0dc8c9..3c073d975 100644 --- a/test/FESpacesTests/DivConformingFESpacesTests.jl +++ b/test/FESpacesTests/DivConformingFESpacesTests.jl @@ -47,15 +47,122 @@ function test_div_v_q_equiv(U,V,P,Q,Ω) @test norm(v1-v2) < tol end +@testset "Test Raviart-Thomas" begin + + domain =(0,1,0,1) + partition = (3,3) + model = CartesianDiscreteModel(domain,partition) + + order = 1 + + u(x) = x + + reffe = ReferenceFE(raviart_thomas,order) + + V = TestFESpace(model,reffe,dirichlet_tags = [1,6]) + test_single_field_fe_space(V) + U = TrialFESpace(V,u) + + reffe = ReferenceFE(lagrangian,Float64,order) + Q = TestFESpace(model,reffe,conformity=:L2) + P = TrialFESpace(Q) + + uh = interpolate(u,U) + + e = u - uh + + Ω = Triangulation(model) + dΩ = Measure(Ω,2*order) + + el2 = sqrt(sum( ∫( e⋅e )*dΩ )) + @test el2 < 1.0e-10 + + test_div_v_q_equiv(U,V,P,Q,Ω) + + #using Gridap.Visualization + # + #writevtk(Ω,"trian",nsubcells=10,cellfields=["uh"=>uh]) + + order = 1 + + reffe = ReferenceFE(TET,raviart_thomas,order) + + domain =(0,1,0,1,0,1) + partition = (3,3,3) + model = simplexify(CartesianDiscreteModel(domain,partition)) + + labels = get_face_labeling(model) + dir_tags = Array{Integer}(undef,0) + + V = FESpace(model,reffe,conformity=DivConformity()) + U = TrialFESpace(V,u) + + reffe = ReferenceFE(lagrangian,Float64,order) + Q = TestFESpace(model,reffe,conformity=:L2) + P = TrialFESpace(Q) + + v(x) = VectorValue(-0.5*x[1]+1.0,-0.5*x[2],-0.5*x[3]) + vh = interpolate(v,V) + + e = v - vh + + Ω = Triangulation(model) + dΩ = Measure(Ω,2*order) + + el2 = sqrt(sum( ∫( e⋅e )*dΩ )) + @test el2 < 1.0e-10 + + test_div_v_q_equiv(U,V,P,Q,Ω) + + #using Gridap.Visualization + # + #writevtk(trian,"test",order=3,cellfields=["vh"=>vh]) + + # Tests on manifold + + # Create domain + domain = (0,1,0,1,0,1) + cells = (2,2,2) + model = CartesianDiscreteModel(domain,cells) + + # Restrict model to cube surface + labels = get_face_labeling(model) + bgface_to_mask = get_face_mask(labels,"boundary",2) + Γface_to_bgface = findall(bgface_to_mask) + Dc2Dp3model = DiscreteModelPortion(DiscreteModel(Polytope{2},model),Γface_to_bgface) + + order = 0 + degree = 1 + + reffe_rt = ReferenceFE(raviart_thomas,Float64,order) + V = FESpace(Dc2Dp3model, reffe_rt ; conformity=:HDiv) + U = TrialFESpace(V,u) + reffe = ReferenceFE(lagrangian,Float64,order) + Q = TestFESpace(Dc2Dp3model,reffe,conformity=:L2) + P = TrialFESpace(Q) + uh = FEFunction(V,rand(num_free_dofs(V))) + vh = interpolate_everywhere(uh,V) + + Ω = Triangulation(Dc2Dp3model) + dΩ = Measure(Ω,2*order) + e=sqrt(sum(∫((uh-vh)⋅(uh-vh))dΩ)) + @test e < 1.0e-12 + + test_div_v_q_equiv(U,V,P,Q,Ω) + +end + +@testset "Test BDM" begin + domain =(0,1,0,1) partition = (3,3) -model = CartesianDiscreteModel(domain,partition) +model = CartesianDiscreteModel(domain,partition) |> simplexify order = 1 u(x) = x -reffe = ReferenceFE(raviart_thomas,order) +reffe = ReferenceFE(bdm,order) V = TestFESpace(model,reffe,dirichlet_tags = [1,6]) test_single_field_fe_space(V) @@ -77,13 +184,9 @@ el2 = sqrt(sum( ∫( e⋅e )*dΩ )) test_div_v_q_equiv(U,V,P,Q,Ω) -#using Gridap.Visualization -# -#writevtk(Ω,"trian",nsubcells=10,cellfields=["uh"=>uh]) - order = 1 -reffe = ReferenceFE(TET,raviart_thomas,order) +reffe = ReferenceFE(TET,bdm,order) domain =(0,1,0,1,0,1) partition = (3,3,3) @@ -110,18 +213,10 @@ dΩ = Measure(Ω,2*order) el2 = sqrt(sum( ∫( e⋅e )*dΩ )) @test el2 < 1.0e-10 -test_div_v_q_equiv(U,V,P,Q,Ω) - -#using Gridap.Visualization -# -#writevtk(trian,"test",order=3,cellfields=["vh"=>vh]) - -# Tests on manifold - # Create domain domain = (0,1,0,1,0,1) cells = (2,2,2) -model = CartesianDiscreteModel(domain,cells) +model = CartesianDiscreteModel(domain,cells) |> simplexify # Restrict model to cube surface labels = get_face_labeling(model) @@ -129,11 +224,11 @@ bgface_to_mask = get_face_mask(labels,"boundary",2) Γface_to_bgface = findall(bgface_to_mask) Dc2Dp3model = DiscreteModelPortion(DiscreteModel(Polytope{2},model),Γface_to_bgface) -order = 0 -degree = 1 +order = 1 +degree = 2 -reffe_rt = ReferenceFE(raviart_thomas,Float64,order) -V = FESpace(Dc2Dp3model, reffe_rt ; conformity=:HDiv) +reffe_bdm = ReferenceFE(bdm,Float64,order) +V = FESpace(Dc2Dp3model, reffe_bdm ; conformity=:HDiv) U = TrialFESpace(V,u) reffe = ReferenceFE(lagrangian,Float64,order) Q = TestFESpace(Dc2Dp3model,reffe,conformity=:L2) @@ -148,5 +243,6 @@ e=sqrt(sum(∫((uh-vh)⋅(uh-vh))dΩ)) test_div_v_q_equiv(U,V,P,Q,Ω) +end end # module diff --git a/test/ReferenceFEsTests/BDMRefFEsTests.jl b/test/ReferenceFEsTests/BDMRefFEsTests.jl new file mode 100644 index 000000000..f695d2c19 --- /dev/null +++ b/test/ReferenceFEsTests/BDMRefFEsTests.jl @@ -0,0 +1,101 @@ +module BDMRefFEsTest + +using Test +using Gridap.Polynomials +using Gridap.Fields +using Gridap.TensorValues +using Gridap.Fields: MockField +using Gridap.ReferenceFEs +using Gridap.Geometry +using Gridap.CellData +using Gridap.Arrays +using Gridap.Visualization + +using Gridap +using Gridap.FESpaces + +p = TRI +D = num_dims(TRI) +et = Float64 +order = 1 + +reffe = BDMRefFE(et,p,order) + +@test length(get_prebasis(reffe)) == 6 +@test get_order(get_prebasis(reffe)) == 1 +@test num_dofs(reffe) == 6 +@test Conformity(reffe) == DivConformity() + +prebasis = get_prebasis(reffe) +dof_basis = get_dof_basis(reffe) + +v = VectorValue(3.0,0.0) +field = GenericField(x->v*x[1]) + +cache = return_cache(dof_basis,field) +r = evaluate!(cache, dof_basis, field) +test_dof_array(dof_basis,field,r) + +cache = return_cache(dof_basis,prebasis) +r = evaluate!(cache, dof_basis, prebasis) +test_dof_array(dof_basis,prebasis,r) + + +p = TET +D = num_dims(TET) +et = Float64 +order = 1 + +reffe = BDMRefFE(et,p,order) +test_reference_fe(reffe) +@test length(get_prebasis(reffe)) == 12 +@test num_dofs(reffe) == 12 +@test get_order(get_prebasis(reffe)) == 1 +@test Conformity(reffe) == DivConformity() + +p = TET +D = num_dims(p) +et = Float64 +order = 2 + +reffe = BDMRefFE(et,p,order) +test_reference_fe(reffe) +@test length(get_prebasis(reffe)) == 30 +@test num_dofs(reffe) == 30 +@test get_order(get_prebasis(reffe)) == 2 +@test Conformity(reffe) == DivConformity() + +prebasis = get_prebasis(reffe) +dof_basis = get_dof_basis(reffe) + +v = VectorValue(0.0,3.0,0.0) +field = GenericField(x->v) + +cache = return_cache(dof_basis,field) +r = evaluate!(cache, dof_basis, field) +test_dof_array(dof_basis,field,r) + +cache = return_cache(dof_basis,prebasis) +r = evaluate!(cache, dof_basis, prebasis) +test_dof_array(dof_basis,prebasis,r) + +# Factory function +reffe = ReferenceFE(TET,bdm,1) +@test length(get_prebasis(reffe)) == 12 +@test get_order(get_prebasis(reffe)) == 1 +@test num_dofs(reffe) == 12 +@test Conformity(reffe) == DivConformity() + +reffe = ReferenceFE(TET,bdm,Float64,1) +@test length(get_prebasis(reffe)) == 12 +@test get_order(get_prebasis(reffe)) == 1 +@test num_dofs(reffe) == 12 +@test Conformity(reffe) == DivConformity() + +@test Conformity(reffe,:L2) == L2Conformity() +@test Conformity(reffe,:Hdiv) == DivConformity() +@test Conformity(reffe,:HDiv) == DivConformity() + +@test BDM() == bdm + +end # module diff --git a/test/ReferenceFEsTests/RaviartThomasRefFEsTests.jl b/test/ReferenceFEsTests/RaviartThomasRefFEsTests.jl index c89cc31b8..1bd40ca7b 100644 --- a/test/ReferenceFEsTests/RaviartThomasRefFEsTests.jl +++ b/test/ReferenceFEsTests/RaviartThomasRefFEsTests.jl @@ -50,6 +50,26 @@ cache = return_cache(dof_basis,prebasis) r = evaluate!(cache, dof_basis, prebasis) test_dof_array(dof_basis,prebasis,r) +### + +order = 0 +p = TRI +D = num_dims(TRI) +et = Float64 + +reffe = RaviartThomasRefFE(et,p,order) + +dofs = get_dof_basis(reffe) +nodes, nf_nodes, nf_moments = get_nodes(dofs), + get_face_nodes_dofs(dofs), + get_face_moments(dofs) +nodes +nf_nodes +nf_moments + +### + + p = TET D = num_dims(TET)