Merge pull request #299 from gridap/CDLagrangianRefFEs

Cd lagrangian ref f es

Author: santiagobadia

src/ReferenceFEs/CDLagrangianRefFEs.jl

@@ -0,0 +1,67 @@
+const CONT = 0
+const DISC = 1
+
+function CDLagrangianRefFE(::Type{T},p::Polytope{D},order::Int,cont) where {T,D}
+  orders = tfill(order,Val{D}())
+  CDLagrangianRefFE(T,p,orders,cont)
+end
+
+function CDLagrangianRefFE(::Type{T},p::Polytope{D},orders,cont) where {T,D}
+
+  lrfe = LagrangianRefFE(T,p,orders)
+  grfe = lrfe.data
+
+  face_own_nodes = _compute_cd_face_own_nodes(p,orders,cont)
+
+  face_own_dofs = _generate_face_own_dofs(face_own_nodes, grfe.dofs.node_and_comp_to_dof)
+
+  reffaces = compute_lagrangian_reffaces(T,p,orders)
+
+  _reffaces = vcat(reffaces...)
+
+  ndofs = length(grfe.dofs.dof_to_node)
+
+  face_own_dofs_permutations = _trivial_face_own_dofs_permutations(face_own_dofs)
+
+  GenericRefFE(
+      grfe.ndofs,
+      grfe.polytope,
+      grfe.prebasis,
+      grfe.dofs,
+      face_own_dofs,
+      face_own_dofs_permutations,
+      grfe.face_dofs, 
+      grfe.shapefuns)
+
+end
+
+function _compute_cd_face_own_nodes(p::ExtrusionPolytope{D},orders::NTuple{D,<:Int},cont::NTuple{D,<:Int}) where D
+  nodes_ijk = CartesianIndices(Tuple(orders.+1))
+  nodes_i = LinearIndices(nodes_ijk)
+  face_owned_nodes = Vector{Int64}[]
+  for nf in p.dface.nfaces
+    anc = nf.anchor
+    ext = nf.extrusion
+    fns = UnitRange{Int64}[]
+    for (i,(e,a,c,o)) in enumerate(zip(ext,anc,cont,orders))
+      if e==0 && c == DISC
+        push!(fns,1:-1)
+        break
+      elseif e==0 && c == CONT
+        push!(fns,a*o+1:a*o+1)
+      elseif e==1 && c == DISC
+        push!(fns,1:o+1)
+      elseif e==1 && c == CONT
+        push!(fns,2:o)
+      end
+    end
+    n_ijk = CartesianIndices(Tuple(fns))
+    if length(n_ijk) == 0
+      n_i = Int64[]
+    else
+      n_i = Int64[nodes_i[n_ijk]...]
+    end
+    push!(face_owned_nodes,n_i)
+  end
+  face_owned_nodes
+end

src/ReferenceFEs/ReferenceFEs.jl

@@ -142,12 +142,16 @@ export QUAD4
 export TET4
 export HEX8
 
+export CONT
+export DISC
+
 export SerendipityRefFE
 export DiscRefFE
 export PDiscRefFE
 export QDiscRefFE
 export RaviartThomasRefFE
 export NedelecRefFE
+export CDLagrangianRefFE
 
 include("Polytopes.jl")
 
@@ -173,4 +177,6 @@ include("RaviartThomasRefFEs.jl")
 
 include("NedelecRefFEs.jl")
 
+include("CDLagrangianRefFEs.jl")
+
 end # module

test/FESpacesTests/CDLagrangianFESpacesTests.jl

@@ -0,0 +1,38 @@
+module CDLagrangianFESpacesTests
+
+using Test
+using Gridap
+using Gridap.Geometry
+using Gridap.ReferenceFEs
+using Gridap.FESpaces
+
+domain =(0,1,0,1)
+partition = (3,3)
+model = CartesianDiscreteModel(domain,partition)
+trian = get_triangulation(model)
+
+order = 1
+
+T = VectorValue{2,Float64}
+reffe = CDLagrangianRefFE(T,QUAD,order,(DISC,CONT))
+V = FESpace(model=model,reffe=reffe)#),dirichlet_tags = [1,6])
+test_single_field_fe_space(V)
+
+u(x) = x
+U = TrialFESpace(V,u)
+
+uh = interpolate(U,u)
+
+e = u - uh
+
+trian = Triangulation(model)
+quad = CellQuadrature(trian,order)
+
+el2 = sqrt(sum(integrate(inner(e,e),trian,quad)))
+@test el2 < 1.0e-10
+
+#using Gridap.Visualization
+#
+#writevtk(trian,"trian",nsubcells=10,cellfields=["uh"=>uh])
+
+end # module

test/FESpacesTests/runtests.jl

@@ -34,6 +34,8 @@ using Test
 
 @testset "DiscontinuousFESpaces" begin include("DiscontinuousFESpacesTests.jl") end
 
+@testset "CDLagrangianFESpaces" begin include("CDLagrangianFESpacesTests.jl") end
+
 @testset "FETerms" begin include("FETermsTests.jl") end
 
 @testset "CellKernels" begin include("CellKernelsTests.jl") end

test/ReferenceFEsTests/CDLagrangianRefFEsTests.jl

@@ -0,0 +1,18 @@
+module CDLagrangianRefFEsTests
+
+using Gridap
+using Gridap.TensorValues
+using Gridap.ReferenceFEs
+
+T = Float64
+reffe = CDLagrangianRefFE(T,SEGMENT,(2,),(DISC,))
+test_reference_fe(reffe)
+
+reffe = CDLagrangianRefFE(T,QUAD,(2,2),(CONT,DISC))
+test_reference_fe(reffe)
+
+T = VectorValue{3,Float64}
+reffe = CDLagrangianRefFE(T,HEX,(2,2,2),(CONT,CONT,DISC))
+test_reference_fe(reffe)
+
+end # module

test/ReferenceFEsTests/runtests.jl

@@ -26,4 +26,6 @@ using Test
 
 @testset "NedelecRefFEs" begin include("NedelecRefFEsTests.jl") end
 
+@testset "CDLagrangianRefFEs" begin include("CDLagrangianRefFEsTests.jl") end
+
 end # module