Crouzeix Raviart Finite Element

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

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. 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)
+  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)

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")

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

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)

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
+
+
+