From bb7e2528dbe917709bacbfba06f9d3b37c32c784 Mon Sep 17 00:00:00 2001 From: MU00232207 Date: Thu, 5 Dec 2024 13:59:04 +1100 Subject: [PATCH 1/2] CR working in 2D --- src/ReferenceFEs/CRRefFEs.jl | 59 +++++++++++++++++++++ src/ReferenceFEs/MomentBasedReferenceFEs.jl | 43 +++++++++++++++ src/ReferenceFEs/ReferenceFEs.jl | 5 ++ test/ReferenceFEsTests/BDMRefFEsTests.jl | 6 +-- 4 files changed, 110 insertions(+), 3 deletions(-) create mode 100644 src/ReferenceFEs/CRRefFEs.jl diff --git a/src/ReferenceFEs/CRRefFEs.jl b/src/ReferenceFEs/CRRefFEs.jl new file mode 100644 index 000000000..261fcd64b --- /dev/null +++ b/src/ReferenceFEs/CRRefFEs.jl @@ -0,0 +1,59 @@ +struct CR <: ReferenceFEName end +const cr = CR() + +""" +CRRefFE(::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 P space of degree `order-1`. + +""" +function CRRefFE(::Type{T},p::Polytope,order::Integer) where T + D = num_dims(p) + + if is_simplex(p) && order == 1 + prebasis = MonomialBasis{D}(T,order,Polynomials._p_filter) + fb = MonomialBasis{D-1}(T,0,Polynomials._p_filter) + else + @notimplemented "CR Reference FE only available for simplices and lowest order" + end + + function fmom(φ,μ,ds) # Face moment function : σ_F(φ,μ) = 1/|F| ( ∫((φ)*μ)dF ) + D = num_dims(ds.cpoly) + facet_measure = get_facet_measure(ds.cpoly, D-1) + facet_measure_1 = Gridap.Fields.ConstantField(1 / facet_measure[ds.face]) + φμ = Broadcasting(Operation(⋅))(φ,μ) + Broadcasting(Operation(*))(φμ,facet_measure_1) + end + + moments = Tuple[ + (get_dimrange(p,D-1),fmom,fb), # Face moments + ] + + return Gridap.ReferenceFEs.MomentBasedReferenceFE(CR(),p,prebasis,moments,L2Conformity()) +end + + +function ReferenceFE(p::Polytope,::CR, order) + CRRefFE(Float64,p,order) +end + +function ReferenceFE(p::Polytope,::CR,::Type{T}, order) where T + CRRefFE(T,p,order) +end + +function Conformity(reffe::GenericRefFE{CR},sym::Symbol) + if sym == :L2 + L2Conformity() + else + @unreachable """\n + It is not possible to use conformity = $sym on a CR reference FE. + + Possible values of conformity for this reference fe are $((:L2, hdiv...)). + """ + end +end + +function get_face_own_dofs(reffe::GenericRefFE{CR}, conf::L2Conformity) + get_face_dofs(reffe) +end diff --git a/src/ReferenceFEs/MomentBasedReferenceFEs.jl b/src/ReferenceFEs/MomentBasedReferenceFEs.jl index c032e4c89..7974016d2 100644 --- a/src/ReferenceFEs/MomentBasedReferenceFEs.jl +++ b/src/ReferenceFEs/MomentBasedReferenceFEs.jl @@ -166,6 +166,49 @@ function get_extension(m::FaceMeasure{Df,Dc}) where {Df,Dc} return ConstantField(TensorValue(hcat([vs[2]-vs[1]...],[vs[3]-vs[1]...]))) end +# TO DO: Bug in 3D, when n==2; n==4 and D==3 +function get_facet_measure(p::Polytope{D}, face::Int) where D + measures = Float64[] + facet_entities = get_face_coordinates(p) + for entity in facet_entities + n = length(entity) + if n == 1 + push!(measures, 0.0) # A point has zero measure + elseif n == 2 + # Length of an edge + p1, p2 = entity + push!(measures, norm(p2-p1)) + elseif n == 3 && D == 2 + # Perimeter of the closed polygon + n = length(entity) + perimeter = 0.0 + for i in 1:n + p1, p2 = entity[i], entity[mod1(i+1, n)] # cyclic indices + perimeter += norm([p2[1] - p1[1], p2[2] - p1[2]]) + end + push!(measures, perimeter) + elseif n == 3 && D == 3 + # Area of a simplex + p1, p2, p3 = entity + v1 = [p2[i] - p1[i] for i in 1:D] + v2 = [p3[i] - p1[i] for i in 1:D] + area = 0.5 * norm(cross(v1, v2)) + push!(measures, area) + elseif n == 4 && D == 3 + # Volume of a tetrahedron ( To do: Should be perimeter of the tetrahedron.) + p1, p2, p3, p4 = entity + v1 = [p2[i] - p1[i] for i in 1:D] + v2 = [p3[i] - p1[i] for i in 1:D] + v3 = [p4[i] - p1[i] for i in 1:D] + volume = abs(dot(v1, cross(v2, v3))) / 6 + push!(measures, volume) + end + + end + dim = get_dimranges(p)[face+1] + return measures[dim] +end + function Arrays.return_cache( σ::Function, # σ(φ,μ,ds) -> Field/Array{Field} φ::AbstractArray{<:Field}, # φ: prebasis (defined on the cell) diff --git a/src/ReferenceFEs/ReferenceFEs.jl b/src/ReferenceFEs/ReferenceFEs.jl index 44cbc6bee..4799cb396 100644 --- a/src/ReferenceFEs/ReferenceFEs.jl +++ b/src/ReferenceFEs/ReferenceFEs.jl @@ -56,6 +56,7 @@ export get_dimranges export get_dimrange export get_vertex_coordinates export get_facet_normal +export get_facet_measure export get_facet_orientations export get_edge_tangent export get_vertex_permutations @@ -182,6 +183,7 @@ export BDMRefFE export NedelecRefFE export BezierRefFE export ModalC0RefFE +export CRRefFE export Lagrangian export DivConforming @@ -197,6 +199,7 @@ export bdm export nedelec export bezier export modalC0 +export cr export Quadrature export QuadratureName @@ -250,6 +253,8 @@ include("BDMRefFEs.jl") include("NedelecRefFEs.jl") +include("CRRefFEs.jl") + include("MockDofs.jl") include("BezierRefFEs.jl") diff --git a/test/ReferenceFEsTests/BDMRefFEsTests.jl b/test/ReferenceFEsTests/BDMRefFEsTests.jl index f695d2c19..7b124ae18 100644 --- a/test/ReferenceFEsTests/BDMRefFEsTests.jl +++ b/test/ReferenceFEsTests/BDMRefFEsTests.jl @@ -1,4 +1,4 @@ -module BDMRefFEsTest +# module BDMRefFEsTest using Test using Gridap.Polynomials @@ -34,7 +34,7 @@ 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) +@enter test_dof_array(dof_basis,field,r) cache = return_cache(dof_basis,prebasis) r = evaluate!(cache, dof_basis, prebasis) @@ -98,4 +98,4 @@ reffe = ReferenceFE(TET,bdm,Float64,1) @test BDM() == bdm -end # module +# end # module From 289476aaee50609f7434f93fc779075eaad6029e Mon Sep 17 00:00:00 2001 From: MU00232207 Date: Mon, 9 Dec 2024 15:14:00 +1100 Subject: [PATCH 2/2] CRRefFEs working in 2D --- src/ReferenceFEs/MomentBasedReferenceFEs.jl | 2 +- test/ReferenceFEsTests/CRRefFEsTests.jl | 33 ++++ .../ReferenceFEsTests/CrouziexRaviartTests.jl | 148 ++++++++++++++++++ 3 files changed, 182 insertions(+), 1 deletion(-) create mode 100644 test/ReferenceFEsTests/CRRefFEsTests.jl create mode 100644 test/ReferenceFEsTests/CrouziexRaviartTests.jl diff --git a/src/ReferenceFEs/MomentBasedReferenceFEs.jl b/src/ReferenceFEs/MomentBasedReferenceFEs.jl index 7974016d2..1d6a94e02 100644 --- a/src/ReferenceFEs/MomentBasedReferenceFEs.jl +++ b/src/ReferenceFEs/MomentBasedReferenceFEs.jl @@ -166,7 +166,7 @@ function get_extension(m::FaceMeasure{Df,Dc}) where {Df,Dc} return ConstantField(TensorValue(hcat([vs[2]-vs[1]...],[vs[3]-vs[1]...]))) end -# TO DO: Bug in 3D, when n==2; n==4 and D==3 +# TO DO: Bug in 3D, when n==2; n==4 and D==3. Also, working on this to make better and more general. function get_facet_measure(p::Polytope{D}, face::Int) where D measures = Float64[] facet_entities = get_face_coordinates(p) diff --git a/test/ReferenceFEsTests/CRRefFEsTests.jl b/test/ReferenceFEsTests/CRRefFEsTests.jl new file mode 100644 index 000000000..24869c0d8 --- /dev/null +++ b/test/ReferenceFEsTests/CRRefFEsTests.jl @@ -0,0 +1,33 @@ +using Gridap +using Gridap.ReferenceFEs +using Gridap.Geometry +using Gridap.Fields +using Gridap.Arrays +using Gridap.ReferenceFEs +using Gridap.Polynomials +using Gridap.Helpers + +using FillArrays + + +p = TRI +D = num_dims(p) + +# T = Float64 +T = VectorValue{D,Float64} + + +cr_reffe = CRRefFE(T,p,1) + +cr_dofs = get_dof_basis(cr_reffe) +cr_prebasis = get_prebasis(cr_reffe) +cr_shapefuns = get_shapefuns(cr_reffe) + +M = evaluate(cr_dofs, cr_shapefuns) + +partition = (0,1,0,1) +cells = (2,2) +model = simplexify(CartesianDiscreteModel(partition, cells)) + +V = FESpace(model,cr_reffe) +get_cell_dof_ids(V) \ No newline at end of file diff --git a/test/ReferenceFEsTests/CrouziexRaviartTests.jl b/test/ReferenceFEsTests/CrouziexRaviartTests.jl new file mode 100644 index 000000000..ed867bc3a --- /dev/null +++ b/test/ReferenceFEsTests/CrouziexRaviartTests.jl @@ -0,0 +1,148 @@ +module CrouziexRaviartTests + + using Gridap + using Gridap.ReferenceFEs, Gridap.Geometry, Gridap.FESpaces, Gridap.Arrays, Gridap.TensorValues + using Gridap.Helpers + + + function solve_crScalarPoisson(partition, cells, u_exact) + + f(x) = - Δ(u_exact)(x) + + model = simplexify(CartesianDiscreteModel(partition, cells)) + + # reffe = CRRefFE(VectorValue{2,Float64},TRI,1) + reffe = CRRefFE(Float64,TRI,1) + V = FESpace(model,reffe,dirichlet_tags="boundary") + U = TrialFESpace(V,u_exact) + + Ω = Triangulation(model) + dΩ = Measure(Ω,3) + + a(u,v) = ∫( ∇(u)⋅∇(v) )*dΩ + l(v) = ∫( f*v )*dΩ + + op = AffineFEOperator(a,l,U,V) + uh = solve(op) + e = uh-u_exact + + return sqrt(sum( ∫( e⋅e )*dΩ )) + end + + + function solve_crVectorPoisson(partition, cells, u_exact) + + f(x) = - Δ(u_exact)(x) + + model = simplexify(CartesianDiscreteModel(partition, cells)) + + reffe = CRRefFE(VectorValue{2,Float64},TRI,1) + V = FESpace(model,reffe,dirichlet_tags="boundary") + U = TrialFESpace(V, u_exact) + + Ω = Triangulation(model) + dΩ = Measure(Ω,3) + + a(u,v) = ∫( ∇(u) ⊙ ∇(v) )*dΩ + l(v) = ∫( f ⋅ v )*dΩ + + op = AffineFEOperator(a,l,U,V) + uh = solve(op) + e = uh-u_exact + + return sqrt(sum( ∫( e⋅e )*dΩ )) + end + + function solve_crStokes(partition, cells, u_exact, p_exact) + f(x) = - Δ(u_exact)(x) + ∇(p_exact)(x) + model = simplexify(CartesianDiscreteModel(partition, cells)) + + reffe_u = CRRefFE(VectorValue{2,Float64},TRI,1) + reffe_p = ReferenceFE(lagrangian, Float64, 0; space=:P) + V = FESpace(model,reffe_u,dirichlet_tags="boundary") + Q = FESpace(model,reffe_p,conformity=:L2,constraint=:zeromean) + Y = MultiFieldFESpace([V,Q]) + + U = TrialFESpace(V, u_exact) + P = TrialFESpace(Q) + X = MultiFieldFESpace([U,P]) + + Ω = Triangulation(model) + dΩ = Measure(Ω,3) + + a((u,p),(v,q)) = ∫( ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u) )dΩ + l((v,q)) = ∫( v⋅f )dΩ + + op = AffineFEOperator(a,l,X,Y) + uh, ph = solve(op) + eu = uh-u_exact + ep = ph-p_exact + + return sqrt(sum( ∫( eu⋅eu )*dΩ )), sqrt(sum( ∫( ep⋅ep )*dΩ )) + end + + function conv_test_Stokes(partition,ns,u,p) + el2u = Float64[] + el2p = Float64[] + hs = Float64[] + for n in ns + l2u, l2p = solve_crStokes(partition,(n,n),u,p) + h = 1.0/n + push!(el2u,l2u) + push!(el2p,l2p) + push!(hs,h) + end + println(el2u) + println(el2p) + el2u, el2p, hs + end + + function conv_test_Poisson(partition,ns,u) + el2 = Float64[] + hs = Float64[] + for n in ns + l2 = solve_crScalarPoisson(partition,(n,n),u) + println(l2) + h = 1.0/n + push!(el2,l2) + push!(hs,h) + end + println(el2) + el2, hs + end + + function slope(hs,errors) + x = log10.(hs) + y = log10.(errors) + linreg = hcat(fill!(similar(x), 1), x) \ y + linreg[2] + end + + + partition = (0,1,0,1) + + # Stokes + u_exact(x) = VectorValue( [sin(pi*x[1])^2*sin(pi*x[2])*cos(pi*x[2]), -sin(pi*x[2])^2*sin(pi*x[1])*cos(pi*x[1])] ) + p_exact(x) = sin(2*pi*x[1])*sin(2*pi*x[2]) + + # Scalar Poisson + u_exact_p(x) = sin(2*π*x[1])*sin(2*π*x[2]) + + # Vector Poisson + # u_exact(x) = VectorValue( [sin(2*π*x[1])*sin(2*π*x[2]), x[1]*(x[1]-1)*x[2]*(x[2]-1)] ) + + ns = [4,8,16,32,64] + + el, hs = conv_test_Poisson(partition,ns,u_exact_p) + # println("Slope L2-norm u_Poisson: $(slope(hs,el))") + + elu, elp, hs = conv_test_Stokes(partition,ns,u_exact,p_exact) + println("Slope L2-norm u_Stokes: $(slope(hs,elu))") + println("Slope L2-norm p_Stokes: $(slope(hs,elp))") + + println("Slope L2-norm u_Poisson: $(slope(hs,el))") + +end # module + + +