Merge branch 'master' into add_missing_multifield_sparsematrixassemblerconstructor

Author: fverdugo

NEWS.md

@@ -8,18 +8,21 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
   - Added missing SparseMatrixAssembler constructor for MultiFieldFESpaces. Since PR [#320](https://github.com/gridap/Gridap.jl/pull/320/).
+  - kw-argument `space` to `LagrangianRefFE` constructor in order to select the type of underlying polynomial space, i.e., `:Q`, `:S`, or `:P`. Since PR [#321](https://github.com/gridap/Gridap.jl/pull/321).
 
 ### Changed
   - The meaning of `inward/outward` has slightly changed for `SkeletonCellBasis` objects. Now, by accessing to these properties a `ReducedSkeletonCellBasis` is returned, which allows to use the result in a more flexible way (in particular, the result can be used in a similar way than the result of `jump` or `mean`). Since PR [#317](https://github.com/gridap/Gridap.jl/pull/317).
-  - Major refactoring in `ReferenceFEs` module. Since PR [#319](https://github.com/gridap/Gridap.jl/pull/319). In particular:
+  - Major refactoring in `ReferenceFEs` module. Since PR [#319](https://github.com/gridap/Gridap.jl/pull/319) and [#321](https://github.com/gridap/Gridap.jl/pull/321). In particular:
     - `NodalReferenceFE` has been replaced by a new abstract type `LagrangianRefFE`.
     - `GenericNodalCartesianRefFE` has been replaced by `GenericLagrangianRefFE`.
-    - `DiscRefFE` replaced by `PDiscRefFE`.
 
 ### Removed
   - Removals associated with the `ReferenceFEs` refactoring in PR [#319](https://github.com/gridap/Gridap.jl/pull/319):
     - Removed `QDiscRefFE` constructor. Use a standard `LagrangianRefFE` and `L2Conformity` instead.
+    - Removed `PDiscRefFE` constructor. Use `LagrangianRefFE` constructor with the kw-argument `space=:P`.
+    - Removed `CDLagrangianRefFE` constructor. Use a standard `LagrangianRefFE` and `CDConformity` instead.
     - Removed fields `face_own_dofs` and `face_own_dof_permutations` from `GenericRefFE`.
+    - Removed struct `DiscRefFE`.
 
 ### Fixed
   - Replaced `+=` by `add_entry!`. Since PR [#316](https://github.com/gridap/Gridap.jl/pull/316).

src/FESpaces/FESpaceFactories.jl

@@ -248,9 +248,7 @@ function _setup_lagrange_spaces(kwargs)
       end
     else
       if reffe == :PLagrangian
-        _reffes = [PDiscRefFE(T,p,order) for p in polytopes]
-      elseif reffe == :QLagrangian
-        _reffes = [QDiscRefFE(T,p,order) for p in polytopes]
+        _reffes = [LagrangianRefFE(T,p,order;space=:P) for p in polytopes]
       else
         @unreachable "Not possible to use a $reffe reffe on polytopes $(polytopes...)"
       end
@@ -287,74 +285,6 @@ function _setup_lagrange_spaces(kwargs)
 
 end
 
-#function _setup_lagrange_spaces(kwargs)
-#
-#  reffe = _get_kwarg(:reffe,kwargs)
-#  model = _get_kwarg(:model,kwargs)
-#  labels = _get_kwarg(:labels,kwargs,nothing)
-#  conformity = _get_kwarg(:conformity,kwargs,true)
-#  diritags = _get_kwarg(:dirichlet_tags,kwargs,Int[])
-#  dirimasks = _get_kwarg(:dirichlet_masks,kwargs,nothing)
-#  order = _get_kwarg(:order,kwargs)
-#
-#  polytopes = get_polytopes(model)
-#  trian = get_triangulation(model)
-#
-#  T = _get_kwarg(:valuetype,kwargs,nothing)
-#  if T == nothing
-#    @unreachable "valuetype is a mandatory keyword argument in FESpace constructor for Lagrangian reference FEs"
-#  end
-#
-#  if _is_reffe_lagrangian_compatible_with_polytopes(reffe,polytopes)
-#
-#    if reffe == :SLagrangian
-#      _reffes = [SerendipityRefFE(T,p,order) for p in polytopes]
-#    else
-#      _reffes = [LagrangianRefFE(T,p,order) for p in polytopes]
-#    end
-#
-#    if conformity in [false, :L2]
-#
-#      s = "Strong dirichlet conditions cannot be imposed in discontinuous spaces for the moment"
-#      @notimplementedif diritags != Int[] s
-#      @notimplementedif dirimasks != nothing s
-#
-#      return  DiscontinuousFESpace(_reffes,trian)
-#
-#    elseif conformity in [true, :default, :H1, :C0]
-#      if labels == nothing
-#        return GradConformingFESpace(_reffes,model,diritags,dirimasks)
-#      else
-#        return GradConformingFESpace(_reffes,model,labels,diritags,dirimasks)
-#      end
-#
-#    else
-#      s = "Conformity $conformity not implemented for $reffe reference FE on polytopes $(polytopes...)"
-#      @unreachable s
-#
-#    end
-#
-#  elseif reffe == :PLagrangian
-#
-#      if conformity in [false, :L2]
-#
-#        _reffes = [PDiscRefFE(T,p,order) for p in polytopes]
-#        return  DiscontinuousFESpace(_reffes,trian)
-#
-#      else
-#
-#        @unreachable "Conformity $conformity not possible for $reffe reference FE on $(polytopes...)"
-#
-#      end
-#
-#  else
-#
-#    @notimplemented "Reference element $reffe not implemented on $(polytopes...)"
-#
-#  end
-#
-#end
-
 function _get_kwarg(kwarg,kwargs)
   try
     return kwargs[kwarg]

src/ReferenceFEs/CDLagrangianRefFEs.jl

@@ -15,20 +15,12 @@ end
 
 # Constructors
 
-function LagrangianRefFE(::Type{T},p::ExtrusionPolytope{D},order::Int,cont) where {T,D}
-  CDLagrangianRefFE(T,p,order,cont)
-end
-
-function LagrangianRefFE(::Type{T},p::ExtrusionPolytope{D},orders,cont) where {T,D}
-  CDLagrangianRefFE(T,p,orders,cont)
-end
-
-function CDLagrangianRefFE(::Type{T},p::ExtrusionPolytope{D},order::Int,cont) where {T,D}
+function _CDLagrangianRefFE(::Type{T},p::ExtrusionPolytope{D},order::Int,cont) where {T,D}
   orders = tfill(order,Val{D}())
-  CDLagrangianRefFE(T,p,orders,cont)
+  _CDLagrangianRefFE(T,p,orders,cont)
 end
 
-function CDLagrangianRefFE(::Type{T},p::ExtrusionPolytope{D},orders,cont) where {T,D}
+function _CDLagrangianRefFE(::Type{T},p::ExtrusionPolytope{D},orders,cont) where {T,D}
   cond(c,o) = ( o > 0 || c == DISC )
   @assert all([cond(c,o) for (c,o) in zip(cont,orders)])
   _cd_lagrangian_ref_fe(T,p,orders,cont)

src/ReferenceFEs/CLagrangianRefFEs.jl

@@ -201,12 +201,26 @@ new polytope types increasing customization possibilities.
 - [`compute_own_nodes_permutations(p::Polytope, interior_nodes)`](@ref)
 - [`compute_lagrangian_reffaces(::Type{T},p::Polytope,orders) where T`](@ref)
 """
-function LagrangianRefFE(::Type{T},p::Polytope{D},orders) where {T,D}
-  if any(orders.==0) && !all(orders.==0)
-    cont = map(i -> i == 0 ? DISC : CONT,orders)
-    return _cd_lagrangian_ref_fe(T,p,orders,cont)
+function LagrangianRefFE(::Type{T},p::Polytope{D},orders;space::Symbol=_default_space(p)) where {T,D}
+  if space == :P && is_n_cube(p)
+    return _PDiscRefFE(T,p,orders)
+  elseif space == :S && is_n_cube(p)
+    SerendipityRefFE(T,p,orders)
   else
-    return _lagrangian_ref_fe(T,p,orders)
+    if any(orders.==0) && !all(orders.==0)
+      cont = map(i -> i == 0 ? DISC : CONT,orders)
+      return _cd_lagrangian_ref_fe(T,p,orders,cont)
+    else
+      return _lagrangian_ref_fe(T,p,orders)
+    end
+  end
+end
+
+function _default_space(p)
+  if is_n_cube(p)
+    :Q
+  else
+    :P
   end
 end
 
@@ -320,43 +334,11 @@ function _generate_face_nodes_aux(
   face_fnode_to_node
 end
 
-#function _find_own_dof_permutaions(node_perms,node_and_comp_to_dof,nfacenodeids,nfacedofsids)
-#  dof_perms = Vector{Int}[]
-#  T = eltype(node_and_comp_to_dof)
-#  ncomps = num_components(T)
-#  idof_to_dof = nfacedofsids[end]
-#  inode_to_node = nfacenodeids[end]
-#  for inode_to_pinode in node_perms
-#    ninodes = length(inode_to_pinode)
-#    nidofs = ncomps*ninodes
-#    idof_to_pidof = fill(INVALID_PERM,nidofs)
-#    for (inode,ipnode) in enumerate(inode_to_pinode)
-#      if ipnode == INVALID_PERM
-#        continue
-#      end
-#      node = inode_to_node[inode]
-#      pnode = inode_to_node[ipnode]
-#      comp_to_pdof = node_and_comp_to_dof[pnode]
-#      comp_to_dof = node_and_comp_to_dof[node]
-#      for comp in 1:ncomps
-#        dof = comp_to_dof[comp]
-#        pdof = comp_to_pdof[comp]
-#        idof = findfirst(i->i==dof,idof_to_dof)
-#        ipdof = findfirst(i->i==pdof,idof_to_dof)
-#        idof_to_pidof[idof] = ipdof
-#      end
-#    end
-#    push!(dof_perms,idof_to_pidof)
-#  end
-#  dof_perms
-#end
-
-
 # Constructors taking Int
 
-function LagrangianRefFE(::Type{T},p::Polytope{D},order::Int) where {T,D}
+function LagrangianRefFE(::Type{T},p::Polytope{D},order::Int;space::Symbol=_default_space(p)) where {T,D}
   orders = tfill(order,Val{D}())
-  LagrangianRefFE(T,p,orders)
+  LagrangianRefFE(T,p,orders;space=space)
 end
 
 function MonomialBasis(::Type{T},p::Polytope{D},order::Int) where {D,T}

src/ReferenceFEs/PDiscRefFEs.jl

@@ -1,38 +1,29 @@
-"""
-    struct PDiscRefFE{D} <: LagrangianRefFE{D}
-      # Private fields
-    end
-"""
-struct PDiscRefFE{D} <: LagrangianRefFE{D}
-  reffe::LagrangianRefFE{D}
-  polytope::Polytope{D}
+
+function _PDiscRefFE(::Type{T},p::Polytope,orders) where T
+    order = first(orders)
+    @notimplementedif any( orders .!= order ) "Anisotropic serentopity FEs not allowed"
+    _PDiscRefFE(T,p,order)
 end
 
-function PDiscRefFE(::Type{T},p::Polytope,order::Integer) where T
+function _PDiscRefFE(::Type{T},p::Polytope,order::Integer) where T
   D = num_cell_dims(p)
   extrusion = tfill(TET_AXIS,Val{D}())
   simplex = ExtrusionPolytope(extrusion)
   reffe = LagrangianRefFE(T,simplex,order)
-  PDiscRefFE{D}(reffe,p)
+  metadata = nothing
+  face_nodes = [Int[] for face in 1:num_faces(p)]
+  face_nodes[end] = collect(1:num_nodes(reffe))
+  dofs = get_dof_basis(reffe)
+  face_dofs = _generate_face_own_dofs(face_nodes,dofs.node_and_comp_to_dof)
+
+  reffe = GenericRefFE(
+    num_dofs(reffe),
+    p,
+    get_prebasis(reffe),
+    get_dof_basis(reffe),
+    L2Conformity(),
+    metadata,
+    face_dofs)
+  GenericLagrangianRefFE(reffe,face_nodes)
 end
 
-# LagrangianRefFE
-
-get_face_nodes(reffe::PDiscRefFE) = get_face_own_nodes(reffe)
-
-# Reffe
-
-num_dofs(reffe::PDiscRefFE) = num_dofs(reffe.reffe)
-
-get_polytope(reffe::PDiscRefFE) = reffe.polytope
-
-get_prebasis(reffe::PDiscRefFE) = get_prebasis(reffe.reffe)
-
-get_dof_basis(reffe::PDiscRefFE) = get_dof_basis(reffe.reffe)
-
-get_default_conformity(reffe::PDiscRefFE) = L2Conformity()
-
-get_face_dofs(reffe::PDiscRefFE) = get_face_own_dofs(reffe)
-
-get_shapefuns(reffe::PDiscRefFE) = get_shapefuns(reffe.reffe)
-

src/ReferenceFEs/ReferenceFEs.jl

@@ -154,10 +154,8 @@ export DivConformity
 export CDConformity
 
 export SerendipityRefFE
-export PDiscRefFE
 export RaviartThomasRefFE
 export NedelecRefFE
-export CDLagrangianRefFE
 
 include("Polytopes.jl")
 

test/FESpacesTests/CDLagrangianFESpacesTests.jl

@@ -16,93 +16,51 @@ orders = (2,1)
 order = 2*max(orders...)
 
 T = VectorValue{2,Float64}
-reffe = CDLagrangianRefFE(T,QUAD,orders,(DISC,CONT))
-V = FESpace(model=model,reffe=reffe)#),dirichlet_tags = [1,6])
+reffe = LagrangianRefFE(T,QUAD,orders)
+V = FESpace(model=model,reffe=reffe,conformity=CDConformity((DISC,CONT)))#),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
 
-reffe = CDLagrangianRefFE(Float64,QUAD,(2,2),(CONT,DISC))
-
-quad9 = LagrangianRefFE(T,QUAD,2)
-V = FESpace(model=model,reffe=quad9,conformity=CDConformity((CONT,DISC)))
-
+reffe = LagrangianRefFE(T,QUAD,2)
+V = FESpace(model=model,reffe=reffe,conformity=CDConformity((CONT,DISC)))
 U = TrialFESpace(V,u)
 uh = interpolate(U,u)
 e = u - uh
-
 el2 = sqrt(sum(integrate(inner(e,e),trian,quad)))
 @test el2 < 1.0e-10
 
-cdquad9 = CDLagrangianRefFE(T,QUAD,2,(CONT,DISC))
-V = FESpace(model=model,reffe=cdquad9)
-
+reffe = LagrangianRefFE(T,QUAD,(2,1))
+V = FESpace(model=model,reffe=reffe,conformity=CDConformity((DISC,CONT)))
 U = TrialFESpace(V,u)
 uh = interpolate(U,u)
 e = u - uh
-
 el2 = sqrt(sum(integrate(inner(e,e),trian,quad)))
 @test el2 < 1.0e-10
 
-cdquad9 = CDLagrangianRefFE(T,QUAD,2,(CONT,DISC))
-V = FESpace(model=model,reffe=cdquad9,conformity=CDConformity((CONT,DISC)))
-
-U = TrialFESpace(V,u)
-uh = interpolate(U,u)
-e = u - uh
-
-el2 = sqrt(sum(integrate(inner(e,e),trian,quad)))
-@test el2 < 1.0e-10
-
-cdquad9 = CDLagrangianRefFE(T,QUAD,2,(CONT,DISC))
-V = FESpace(model=model,reffe=cdquad9,conformity=CDConformity((DISC,CONT)))
-
-U = TrialFESpace(V,u)
-uh = interpolate(U,u)
-e = u - uh
-
-el2 = sqrt(sum(integrate(inner(e,e),trian,quad)))
-@test el2 < 1.0e-10
-
-cdquad9 = LagrangianRefFE(T,QUAD,2,(CONT,DISC))
-V = FESpace(model=model,reffe=cdquad9,conformity=CDConformity((DISC,CONT)))
-
-U = TrialFESpace(V,u)
-uh = interpolate(U,u)
-e = u - uh
-
-el2 = sqrt(sum(integrate(inner(e,e),trian,quad)))
-@test el2 < 1.0e-10
-
-cdquad9 = LagrangianRefFE(T,QUAD,(2,1))
-V = FESpace(model=model,reffe=cdquad9,conformity=CDConformity((DISC,CONT)))
-
+reffe = LagrangianRefFE(T,QUAD,(2,0))
+V = FESpace(model=model,reffe=reffe,conformity=CDConformity((CONT,DISC)))
+u(x) = VectorValue(x[1],0.0)
 U = TrialFESpace(V,u)
 uh = interpolate(U,u)
 e = u - uh
-
 el2 = sqrt(sum(integrate(inner(e,e),trian,quad)))
 @test el2 < 1.0e-10
 
-cdquad9 = LagrangianRefFE(T,QUAD,(2,0))
-V = FESpace(model=model,reffe=cdquad9,conformity=CDConformity((CONT,DISC)))
-
+reffe = LagrangianRefFE(T,QUAD,(2,0))
+V = FESpace(model=model,reffe=reffe,conformity=CDConformity((DISC,DISC)))
 u(x) = VectorValue(x[1],0.0)
 U = TrialFESpace(V,u)
 uh = interpolate(U,u)
 e = u - uh
-
 el2 = sqrt(sum(integrate(inner(e,e),trian,quad)))
 @test el2 < 1.0e-10
 

test/FESpacesTests/DiscontinuousFESpacesTests.jl

@@ -31,7 +31,7 @@ uh = FEFunction(V,rand(num_free_dofs(V)))
 #
 #writevtk(trian,"trian",nsubcells=40,cellfields=["fh"=>fh, "uh"=>uh])
 
-reffes = [PDiscRefFE(Float64,p,order) for p in polytopes]
+reffes = [LagrangianRefFE(Float64,p,order,space=:P) for p in polytopes]
 V = DiscontinuousFESpace(reffes,trian)
 test_single_field_fe_space(V)