ConstraintHandler for discontinuous interpolations

docs/src/devdocs/interpolations.md

@@ -31,6 +31,8 @@ Ferrite.edgedof_interior_indices(::Interpolation)
 Ferrite.celldof_interior_indices(::Interpolation)
 Ferrite.getnbasefunctions(::Interpolation)
 Ferrite.reference_coordinates(::Interpolation)
+Ferrite.IsDiscontinuous(::Interpolation)
+Ferrite.get_continuous_interpolation(::Interpolation)
 ```
 
 for all entities which exist on that reference element. The dof functions default to having no

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,13 @@ function add!(ch::ConstraintHandler, dbc::Dirichlet)
         if interpolation isa VectorizedInterpolation
             interpolation = interpolation.ip
         end
+
+        if getorder(interpolation) == 0 
+            @warn("No dof prescribed for order 0 interpolations")
+            dbc_added = true
+            continue
+        end
+
         # 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)
@@ -845,7 +852,7 @@ function add!(ch::ConstraintHandler, dbc::Dirichlet)
             EntityType = FaceIndex
         end
         CT = getcelltype(ch.dh.grid, first(fh.cellset)) # Same celltype enforced in FieldHandler constructor
-        bcvalues = BCValues(interpolation, default_interpolation(CT), EntityType)
+        bcvalues = BCValues(interpolation , default_interpolation(CT), EntityType)
         # Recreate the Dirichlet(...) struct with the filtered set and call internal add!
         filtered_dbc = Dirichlet(dbc.field_name, filtered_set, dbc.f, components)
         _add!(
@@ -1173,7 +1180,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 +1201,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

@@ -86,6 +86,7 @@ struct InterpolationInfo
     reference_dim::Int
     adjust_during_distribution::Bool
     n_copies::Int
+    is_discontinuous::Bool
     function InterpolationInfo(interpolation::InterpolationByDim{3})
         n_copies = 1
         if interpolation isa VectorizedInterpolation
@@ -100,6 +101,7 @@ struct InterpolationInfo
             3,
             adjust_dofs_during_distribution(interpolation),
             n_copies,
+            IsDiscontinuous(interpolation)
         )
     end
     function InterpolationInfo(interpolation::InterpolationByDim{2})
@@ -116,6 +118,7 @@ struct InterpolationInfo
             2,
             adjust_dofs_during_distribution(interpolation),
             n_copies,
+            IsDiscontinuous(interpolation)
         )
     end
     function InterpolationInfo(interpolation::InterpolationByDim{1})
@@ -131,7 +134,8 @@ struct InterpolationInfo
             length(celldof_interior_indices(interpolation)),
             1,
             adjust_dofs_during_distribution(interpolation),
-            n_copies
+            n_copies,
+            IsDiscontinuous(interpolation)
         )
     end
 end
@@ -249,6 +253,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 +267,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 +292,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 +336,34 @@ 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
+
+"""
+    IsDiscontinuous(::Interpolation)
+    IsDiscontinuous(::Type{<:Interpolation})
+Checks whether the interpolation is discontinuous (i.e. `DiscontinuousLagrange`)
+"""
+IsDiscontinuous(::Interpolation) = false
+IsDiscontinuous(::Type{<:Interpolation}) = false
+
+"""
+    get_continuous_interpolation(::Interpolation)
+    get_continuous_interpolation(::Type{<:Interpolation})
+Returns the continuous version (i.e. `Lagrange`) of a discontinuous interpolation is  (i.e. `DiscontinuousLagrange`).
+Used internally for associating `cell_dofs` with faces.
+The return type follows that of the argument.
+```julia-repl
+julia> Ferrite.get_continuous_interpolation(DiscontinuousLagrange{RefQuadrilateral,1}())
+Lagrange{RefQuadrilateral, 1}()
+julia> Ferrite.get_continuous_interpolation(DiscontinuousLagrange{RefQuadrilateral,1})
+Lagrange{RefQuadrilateral, 1}
+```
+"""
+get_continuous_interpolation(ip::Interpolation) = throw(ArgumentError("Interpolation $ip is already continuous."))
+get_continuous_interpolation(ip::Type{<:Interpolation}) = throw(ArgumentError("Interpolation $ip is already continuous."))
+
 #########################
 # DiscontinuousLagrange #
 #########################
@@ -339,6 +377,12 @@ struct DiscontinuousLagrange{shape, order, unused} <: ScalarInterpolation{shape,
     end
 end
 
+IsDiscontinuous(::Type{<:DiscontinuousLagrange}) = true
+IsDiscontinuous(::DiscontinuousLagrange) = true
+
+get_continuous_interpolation(::DiscontinuousLagrange{ref_shape,order}) where {ref_shape, order} = Lagrange{ref_shape,order}()
+get_continuous_interpolation(::Type{<:DiscontinuousLagrange{ref_shape,order}}) where {ref_shape, order} = Lagrange{ref_shape,order, Nothing}
+
 getlowerorder(::DiscontinuousLagrange{shape,order}) where {shape,order} = DiscontinuousLagrange{shape,order-1}()
 
 getnbasefunctions(::DiscontinuousLagrange{shape,order}) where {shape,order} = getnbasefunctions(Lagrange{shape,order}())
@@ -347,6 +391,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) = dirichlet_facedof_indices(get_continuous_interpolation(ip))
+dirichlet_edge_indices(ip::DiscontinuousLagrange) = dirichlet_edgedof_indices(get_continuous_interpolation(ip))
+dirichlet_vertex_indices(ip::DiscontinuousLagrange) = dirichlet_vertex_indices(get_continuous_interpolation(ip))
+
 # Mirror the Lagrange element for now.
 function reference_coordinates(ip::DiscontinuousLagrange{shape, order}) where {shape, order}
     return reference_coordinates(Lagrange{shape,order}())
@@ -1296,3 +1345,9 @@ function shape_gradient_and_value(ipv::VectorizedInterpolation{vdim, shape}, ξ:
 end
 
 reference_coordinates(ip::VectorizedInterpolation) = reference_coordinates(ip.ip)
+
+IsDiscontinuous(ipv::VectorizedInterpolation) = IsDiscontinuous(ipv.ip)
+IsDiscontinuous(::Type{<:VectorizedInterpolation{vdim, refshape, order, ip}}) where {vdim, refshape, order, ip<:DiscontinuousLagrange}= IsDiscontinuous(ip)
+
+get_continuous_interpolation(ipv::VectorizedInterpolation{vdim}) where {vdim} = VectorizedInterpolation{vdim}(get_continuous_interpolation(ipv.ip))
+get_continuous_interpolation(::Type{<:VectorizedInterpolation{vdim, refshape, order, ip}}) where {vdim, refshape, order, ip<:DiscontinuousLagrange} = VectorizedInterpolation{vdim, refshape, order, get_continuous_interpolation(ip)}

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_logs (:warn,"No dof prescribed for order 0 interpolations") add!(ch0, Dirichlet(:w, Γ, (x, t) -> 0))
     ## Scalar
     dbc = Dirichlet(:s, Γ, (x, t) -> 0)
     add!(ch, dbc)

test/test_interpolations.jl

@@ -163,4 +163,21 @@ end
 @test Ferrite.reference_coordinates(DiscontinuousLagrange{RefTetrahedron,0}()) ≈ [Vec{3,Float64}((1/4,1/4,1/4))]
 @test Ferrite.reference_coordinates(DiscontinuousLagrange{RefHexahedron,0}()) ≈ [Vec{3,Float64}((0,0,0))]
 
+# Test discontinuous interpolations related functions
+d_ip = DiscontinuousLagrange{RefQuadrilateral,1}()
+d_ip_t = DiscontinuousLagrange{RefQuadrilateral,1}
+
+ip = Lagrange{RefQuadrilateral,1}()
+ip_t = Lagrange{RefQuadrilateral,1}
+
+@test Ferrite.IsDiscontinuous(ip) == false
+@test Ferrite.IsDiscontinuous(ip_t) == false
+@test Ferrite.IsDiscontinuous(d_ip) == true
+@test Ferrite.IsDiscontinuous(d_ip_t) == true
+
+@test Ferrite.get_continuous_interpolation(d_ip_t) == typeof(Ferrite.get_continuous_interpolation(d_ip))
+
+# Error paths
+@test_throws ArgumentError("Interpolation $ip is already continuous.") Ferrite.get_continuous_interpolation(ip)
+@test_throws ArgumentError("Interpolation $ip_t is already continuous.") Ferrite.get_continuous_interpolation(ip_t)
 end