ConstraintHandler for discontinuous interpolations

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)
@@ -1173,7 +1174,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
@@ -1194,7 +1195,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/interpolations.jl

@@ -249,6 +249,8 @@ match the vertex enumeration of the corresponding geometrical cell.
 """
 vertexdof_indices(ip::Interpolation) = ntuple(_ -> (), nvertices(ip))
 
+dirichlet_vertexdof_indices(ip::Interpolation) = vertexdof_indices(ip)
+
 """
     edgedof_indices(ip::Interpolation)
 
@@ -261,6 +263,8 @@ Here the first entries are the vertex dofs, followed by the edge interior dofs.
 """
 edgedof_indices(::Interpolation)
 
+dirichlet_edgedof_indices(ip::Interpolation) = edgedof_indices(ip)
+
 """
     edgedof_interior_indices(ip::Interpolation)
 
@@ -284,6 +288,8 @@ the face enumeration of the corresponding geometrical cell.
 """
 facedof_indices(::Interpolation)
 
+dirichlet_facedof_indices(ip::Interpolation) = facedof_indices(ip)
+
 """
     facedof_interior_indices(ip::Interpolation)
 
@@ -326,6 +332,10 @@ 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{FaceIndex}) = Ferrite.dirichlet_facedof_indices
+dirichlet_boundarydof_indices(::Type{EdgeIndex}) = Ferrite.dirichlet_edgedof_indices
+dirichlet_boundarydof_indices(::Type{VertexIndex}) = Ferrite.dirichlet_vertexdof_indices
+
 #########################
 # DiscontinuousLagrange #
 #########################
@@ -347,6 +357,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

@@ -10,10 +10,16 @@
     close!(dh)
     ch = ConstraintHandler(dh)
 
+    dh0 = DofHandler(grid)
+    add!(dh0, :w, DiscontinuousLagrange{RefTriangle,0}())
+    close!(dh0)
+    ch0 = ConstraintHandler(dh0)
+
     # 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])
+    @test_throws ErrorException("No dof prescribed for order 0 interpolations") add!(ch0, Dirichlet(:w, Γ, (x, t) -> 0))
     ## Scalar
     dbc = Dirichlet(:s, Γ, (x, t) -> 0)
     add!(ch, dbc)