ConstraintHandler support for discontinuous interpolations (#729)

docs/src/devdocs/interpolations.md

@@ -14,6 +14,7 @@ Ferrite.getorder(::Interpolation)
 Ferrite.value(::Interpolation{dim}, ::Vec{dim,T}) where {dim,T}
 Ferrite.derivative(::Interpolation{dim}, ::Vec{dim}) where {dim}
 Ferrite.boundarydof_indices
+Ferrite.dirichlet_boundarydof_indices
 ```
 
 ### Required methods to implement for all subtypes of `Interpolation` to define a new finite element
@@ -23,9 +24,12 @@ Depending on the dimension of the reference element the following functions have
 ```@docs
 Ferrite.value(::Interpolation, ::Int, ::Vec)
 Ferrite.vertexdof_indices(::Interpolation)
+Ferrite.dirichlet_vertexdof_indices(::Interpolation)
 Ferrite.facedof_indices(::Interpolation)
+Ferrite.dirichlet_facedof_indices(::Interpolation)
 Ferrite.facedof_interior_indices(::Interpolation)
 Ferrite.edgedof_indices(::Interpolation)
+Ferrite.dirichlet_edgedof_indices(::Interpolation)
 Ferrite.edgedof_interior_indices(::Interpolation)
 Ferrite.celldof_interior_indices(::Interpolation)
 Ferrite.getnbasefunctions(::Interpolation)

src/Dofs/ConstraintHandler.jl

@@ -275,7 +275,7 @@ end
 # Dirichlet on (face|edge|vertex)set
 function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcfaces::Set{Index}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, _) where {Index<:BoundaryIndex}
     local_face_dofs, local_face_dofs_offset =
-        _local_face_dofs_for_bc(interpolation, field_dim, dbc.components, offset, boundarydof_indices(eltype(bcfaces)))
+        _local_face_dofs_for_bc(interpolation, field_dim, dbc.components, offset, dirichlet_boundarydof_indices(eltype(bcfaces)))
     copy!(dbc.local_face_dofs, local_face_dofs)
     copy!(dbc.local_face_dofs_offset, local_face_dofs_offset)
 
@@ -300,7 +300,7 @@ end
 
 # Calculate which local dof index live on each face:
 # face `i` have dofs `local_face_dofs[local_face_dofs_offset[i]:local_face_dofs_offset[i+1]-1]
-function _local_face_dofs_for_bc(interpolation, field_dim, components, offset, boundaryfunc::F=facedof_indices) where F
+function _local_face_dofs_for_bc(interpolation, field_dim, components, offset, boundaryfunc::F=dirichlet_facedof_indices) where F
     @assert issorted(components)
     local_face_dofs = Int[]
     local_face_dofs_offset = Int[1]
@@ -833,6 +833,7 @@ function add!(ch::ConstraintHandler, dbc::Dirichlet)
         if interpolation isa VectorizedInterpolation
             interpolation = interpolation.ip
         end
+        getorder(interpolation) == 0 && error("No dof prescribed for order 0 interpolations")
         # Set up components to prescribe (empty input means prescribe all components)
         components = isempty(dbc.components) ? collect(Int, 1:n_comp) : dbc.components
         if !all(c -> 0 < c <= n_comp, components)
@@ -1174,7 +1175,7 @@ end
 function mirror_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange{<:Union{RefQuadrilateral,RefTriangle}}, n::Int)
     # For 2D we always permute since Ferrite defines dofs counter-clockwise
     ret = collect(1:length(local_face_dofs))
-    for (i, f) in enumerate(facedof_indices(ip))
+    for (i, f) in enumerate(dirichlet_facedof_indices(ip))
         this_offset = local_face_dofs_offset[i]
         other_offset = this_offset + n
         for d in 1:n
@@ -1195,7 +1196,7 @@ function mirror_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange
     ret = collect(1:length(local_face_dofs))
 
     # Mirror by changing from counter-clockwise to clockwise
-    for (i, f) in enumerate(facedof_indices(ip))
+    for (i, f) in enumerate(dirichlet_facedof_indices(ip))
         r = local_face_dofs_offset[i]:(local_face_dofs_offset[i+1] - 1)
         # 1. Rotate the corners
         vertex_range = r[1:(N*n)]

src/FEValues/face_values.jl

@@ -207,7 +207,7 @@ function BCValues(::Type{T}, func_interpol::Interpolation{refshape}, geom_interp
     interpolation_coords = reference_coordinates(func_interpol)
 
     qrs = QuadratureRule{refshape,T,dim}[]
-    for boundarydofs in boundarydof_indices(boundary_type)(func_interpol)
+    for boundarydofs in dirichlet_boundarydof_indices(boundary_type)(func_interpol)
         dofcoords = Vec{dim,T}[]
         for boundarydof in boundarydofs
             push!(dofcoords, interpolation_coords[boundarydof])

src/interpolations.jl

@@ -249,6 +249,20 @@ match the vertex enumeration of the corresponding geometrical cell.
 """
 vertexdof_indices(ip::Interpolation) = ntuple(_ -> (), nvertices(ip))
 
+"""
+    dirichlet_vertexdof_indices(ip::Interpolation)
+
+A tuple containing tuples of local dof indices for the respective vertex in local
+enumeration on a cell defined by [`vertices(::Cell)`](@ref). The vertex enumeration must
+match the vertex enumeration of the corresponding geometrical cell.
+Used internally in [`ConstraintHandler`](@ref) and defaults to [`vertexdof_indices(ip::Interpolation)`](@ref) for continuous interpolation.
+
+!!! note
+    The dofs appearing in the tuple must be continuous and increasing! The first dof must be
+    the 1, as vertex dofs are enumerated first.
+"""
+dirichlet_vertexdof_indices(ip::Interpolation) = vertexdof_indices(ip)
+
 """
     edgedof_indices(ip::Interpolation)
 
@@ -261,6 +275,19 @@ Here the first entries are the vertex dofs, followed by the edge interior dofs.
 """
 edgedof_indices(::Interpolation)
 
+"""
+    dirichlet_edgedof_indices(ip::Interpolation)
+
+A tuple containing tuples of local dof indices for the respective edge in local enumeration
+on a cell defined by [`edges(::Cell)`](@ref). The edge enumeration must match the edge
+enumeration of the corresponding geometrical cell.
+Used internally in [`ConstraintHandler`](@ref) and defaults to [`edgedof_indices(ip::Interpolation)`](@ref) for continuous interpolation.
+
+The dofs are guaranteed to be aligned with the local ordering of the entities on the oriented edge.
+Here the first entries are the vertex dofs, followed by the edge interior dofs.
+"""
+dirichlet_edgedof_indices(ip::Interpolation) = edgedof_indices(ip)
+
 """
     edgedof_interior_indices(ip::Interpolation)
 
@@ -284,6 +311,16 @@ the face enumeration of the corresponding geometrical cell.
 """
 facedof_indices(::Interpolation)
 
+"""
+    dirichlet_facedof_indices(ip::Interpolation)
+
+A tuple containing tuples of all local dof indices for the respective face in local
+enumeration on a cell defined by [`faces(::Cell)`](@ref). The face enumeration must match
+the face enumeration of the corresponding geometrical cell.
+Used internally in [`ConstraintHandler`](@ref) and defaults to [`facedof_indices(ip::Interpolation)`](@ref) for continuous interpolation.
+"""
+dirichlet_facedof_indices(ip::Interpolation) = facedof_indices(ip)
+
 """
     facedof_interior_indices(ip::Interpolation)
 
@@ -326,6 +363,18 @@ boundarydof_indices(::Type{FaceIndex}) = Ferrite.facedof_indices
 boundarydof_indices(::Type{EdgeIndex}) = Ferrite.edgedof_indices
 boundarydof_indices(::Type{VertexIndex}) = Ferrite.vertexdof_indices
 
+"""
+    dirichlet_boundarydof_indices(::Type{<:BoundaryIndex})
+
+Helper function to generically dispatch on the correct dof sets of a boundary entity.
+Used internally in [`ConstraintHandler`](@ref) and defaults to [`boundarydof_indices(ip::Interpolation)`](@ref) for continuous interpolation.
+"""
+dirichlet_boundarydof_indices(::Type{<:BoundaryIndex})
+
+dirichlet_boundarydof_indices(::Type{FaceIndex}) = Ferrite.dirichlet_facedof_indices
+dirichlet_boundarydof_indices(::Type{EdgeIndex}) = Ferrite.dirichlet_edgedof_indices
+dirichlet_boundarydof_indices(::Type{VertexIndex}) = Ferrite.dirichlet_vertexdof_indices
+
 #########################
 # DiscontinuousLagrange #
 #########################
@@ -339,6 +388,8 @@ struct DiscontinuousLagrange{shape, order, unused} <: ScalarInterpolation{shape,
     end
 end
 
+adjust_dofs_during_distribution(::DiscontinuousLagrange) = false
+
 getlowerorder(::DiscontinuousLagrange{shape,order}) where {shape,order} = DiscontinuousLagrange{shape,order-1}()
 
 getnbasefunctions(::DiscontinuousLagrange{shape,order}) where {shape,order} = getnbasefunctions(Lagrange{shape,order}())
@@ -347,6 +398,11 @@ getnbasefunctions(::DiscontinuousLagrange{shape,0}) where {shape} = 1
 # This just moves all dofs into the interior of the element.
 celldof_interior_indices(ip::DiscontinuousLagrange) = ntuple(i->i, getnbasefunctions(ip))
 
+# Mirror the Lagrange element for now to avoid repeating.
+dirichlet_facedof_indices(ip::DiscontinuousLagrange{shape, order}) where {shape, order} = dirichlet_facedof_indices(Lagrange{shape, order}())
+dirichlet_edgedof_indices(ip::DiscontinuousLagrange{shape, order}) where {shape, order} = dirichlet_edgedof_indices(Lagrange{shape, order}())
+dirichlet_vertexdof_indices(ip::DiscontinuousLagrange{shape, order}) where {shape, order} = dirichlet_vertexdof_indices(Lagrange{shape, order}())
+
 # Mirror the Lagrange element for now.
 function reference_coordinates(ip::DiscontinuousLagrange{shape, order}) where {shape, order}
     return reference_coordinates(Lagrange{shape,order}())

test/test_constraints.jl

@@ -7,34 +7,38 @@
     dh = DofHandler(grid)
     add!(dh, :s, Lagrange{RefTriangle,1}())
     add!(dh, :v, Lagrange{RefTriangle,1}()^2)
+    add!(dh, :z, DiscontinuousLagrange{RefTriangle,0}())
+    add!(dh, :sd, DiscontinuousLagrange{RefTriangle,1}())
+    add!(dh, :vd, DiscontinuousLagrange{RefTriangle,1}()^2)
     close!(dh)
     ch = ConstraintHandler(dh)
-
     # Dirichlet
     @test_throws ErrorException("components are empty: $(Int)[]") Dirichlet(:u, Γ, (x, t) -> 0, Int[])
     @test_throws ErrorException("components not sorted: [2, 1]") Dirichlet(:u, Γ, (x, t) -> 0, Int[2, 1])
     @test_throws ErrorException("components not unique: [2, 2]") Dirichlet(:u, Γ, (x, t) -> 0, Int[2, 2])
-    ## Scalar
-    dbc = Dirichlet(:s, Γ, (x, t) -> 0)
-    add!(ch, dbc)
-    @test ch.dbcs[end].components == [1]
-    dbc = Dirichlet(:s, Γ, (x, t) -> 0, [1])
-    add!(ch, dbc)
-    @test ch.dbcs[end].components == [1]
-    dbc = Dirichlet(:s, Γ, (x, t) -> 0, [1, 2])
-    @test_throws ErrorException("components [1, 2] not within range of field :s (1 dimension(s))") add!(ch, dbc)
-    dbc = Dirichlet(:p, Γ, (x, t) -> 0)
-    @test_throws ErrorException("No overlap between dbc::Dirichlet and fields in the ConstraintHandler's DofHandler") add!(ch, dbc)
-    #    ## Vector
-    dbc = Dirichlet(:v, Γ, (x, t) -> 0)
-    add!(ch, dbc)
-    @test ch.dbcs[end].components == [1, 2]
-    dbc = Dirichlet(:v, Γ, (x, t) -> 0, [1])
-    add!(ch, dbc)
-    @test ch.dbcs[end].components == [1]
-    dbc = Dirichlet(:v, Γ, (x, t) -> 0, [2, 3])
-    @test_throws ErrorException("components [2, 3] not within range of field :v (2 dimension(s))") add!(ch, dbc)
-
+    @test_throws ErrorException("No dof prescribed for order 0 interpolations") add!(ch, Dirichlet(:z, Γ, (x, t) -> 0))
+    for (s, v) in [(:s, :v), (:sd, :vd)] 
+        ## Scalar
+        dbc = Dirichlet(s, Γ, (x, t) -> 0)
+        add!(ch, dbc)
+        @test ch.dbcs[end].components == [1]
+        dbc = Dirichlet(s, Γ, (x, t) -> 0, [1])
+        add!(ch, dbc)
+        @test ch.dbcs[end].components == [1]
+        dbc = Dirichlet(s, Γ, (x, t) -> 0, [1, 2])
+        @test_throws ErrorException("components [1, 2] not within range of field :$s (1 dimension(s))") add!(ch, dbc)
+        dbc = Dirichlet(:p, Γ, (x, t) -> 0)
+        @test_throws ErrorException("No overlap between dbc::Dirichlet and fields in the ConstraintHandler's DofHandler") add!(ch, dbc)
+        ## Vector
+        dbc = Dirichlet(v, Γ, (x, t) -> 0)
+        add!(ch, dbc)
+        @test ch.dbcs[end].components == [1, 2]
+        dbc = Dirichlet(v, Γ, (x, t) -> 0, [1])
+        add!(ch, dbc)
+        @test ch.dbcs[end].components == [1]
+        dbc = Dirichlet(v, Γ, (x, t) -> 0, [2, 3])
+        @test_throws ErrorException("components [2, 3] not within range of field :$v (2 dimension(s))") add!(ch, dbc)
+    end
     # PeriodicDirichlet
     @test_throws ErrorException("components are empty: $(Int)[]") PeriodicDirichlet(:u, face_map, Int[])
     @test_throws ErrorException("components not sorted: [2, 1]") PeriodicDirichlet(:u, face_map, Int[2, 1])
@@ -52,7 +56,7 @@
     pdbc = PeriodicDirichlet(:s, face_map, [1, 2])
     @test_throws ErrorException("components [1, 2] not within range of field :s (1 dimension(s))") add!(ConstraintHandler(dh), pdbc)
     pdbc = PeriodicDirichlet(:p, face_map)
-    @test_throws ErrorException("Did not find field :p in DofHandler (existing fields: [:s, :v]).") add!(ConstraintHandler(dh), pdbc)
+    @test_throws ErrorException("Did not find field :p in DofHandler (existing fields: [:s, :v, :z, :sd, :vd]).") add!(ConstraintHandler(dh), pdbc)
     ## Vector
     pdbc = PeriodicDirichlet(:v, face_map)
     add!(ConstraintHandler(dh), pdbc)
@@ -187,6 +191,53 @@ end
     @test ch.prescribed_dofs == [10, 11, 14, 25, 27]
 end
 
+@testset "discontinuous ip constraints" begin
+    grid = generate_grid(Hexahedron, (1, 1, 1))
+    addedgeset!(grid, "bottom", x-> x[3] ≈ -1.0)
+    dh = DofHandler(grid)
+    add!(dh, :u, DiscontinuousLagrange{RefHexahedron,1}()^3)
+    add!(dh, :p, DiscontinuousLagrange{RefHexahedron,1}())
+    close!(dh)
+
+    face_ch = ConstraintHandler(dh)
+    face_dbc = Dirichlet(:u, getfaceset(grid, "bottom"), (x,t) -> x, [1, 2, 3])
+    add!(face_ch, face_dbc)
+    close!(face_ch)
+    update!(face_ch)
+
+    edge_ch = ConstraintHandler(dh)
+    edge_dbc = Dirichlet(:u, getedgeset(grid, "bottom"), (x,t) -> x, [1, 2, 3])
+    add!(edge_ch, edge_dbc)
+    close!(edge_ch)
+    update!(edge_ch)
+
+    @test edge_ch.prescribed_dofs == face_ch.prescribed_dofs == collect(1:12)
+    @test edge_ch.inhomogeneities == face_ch.inhomogeneities == reduce(vcat, getcoordinates(grid,1)[1:4])
+
+    # This can be merged with the continuous test or removed.
+    # Shell mesh edge bcs
+    nodes = [Node{3,Float64}(Vec(0.0,0.0,0.0)), Node{3,Float64}(Vec(1.0,0.0,0.0)), 
+             Node{3,Float64}(Vec(1.0,1.0,0.0)), Node{3,Float64}(Vec(0.0,1.0,0.0)),
+             Node{3,Float64}(Vec(2.0,0.0,0.0)), Node{3,Float64}(Vec(2.0,2.0,0.0))]
+
+    cells = [Quadrilateral((1,2,3,4)), Quadrilateral((2,5,6,3))]
+    grid = Grid(cells,nodes)
+
+    # 3d quad with 1st order 2d interpolation
+    dh = DofHandler(grid)
+    add!(dh, :u, DiscontinuousLagrange{RefQuadrilateral,2}())
+    add!(dh, :θ, DiscontinuousLagrange{RefQuadrilateral,2}())
+    close!(dh)
+
+    addedgeset!(grid, "bottom", x -> x[2] ≈ 0.0) #bottom edge
+    edge_ch = ConstraintHandler(dh)
+    edge_dbc = Dirichlet(:θ, getedgeset(grid, "bottom"), (x,t) -> (0.0,), [1])
+    add!(edge_ch, edge_dbc)
+    close!(edge_ch)
+    update!(edge_ch)
+
+    @test edge_ch.prescribed_dofs == [10, 11, 14, 28, 29, 32]
+end
 #@testset "edge bc mixed grid" begin
 #    # Test mesh
 #    # 6---7---8---9---10