Rework code for add!(::ConstraintHandler, ::Dirichlet) (#672)

This patch cleans up the code for adding `Dirichlet` constraints to the
`ConstraintHandler` a bit after the DofHandler merge. In particular, it
deprecates the method
```
add!(::ConstraintHandler, ::FieldHandler, ::Dirichlet)
```
in favor of just
```
add!(::ConstraintHandler, ::Dirichlet)
```
which was added in #427. There is no need to be able to specify the
constrained set by both the set passed to `Dirichlet` and by the set
given to the `FieldHandler`. Looking at the tests modified in this
patch, it seems this was never the intention anyway. This patch also
removes the copy of the set in the `Dirichlet` constructor in favor of
explicit partitioning of the set over the `FieldHandler`s.

src/Dofs/ConstraintHandler.jl

@@ -40,7 +40,7 @@ struct Dirichlet # <: Constraint
     local_face_dofs_offset::Vector{Int}
 end
 function Dirichlet(field_name::Symbol, faces::Set, f::Function, components=nothing)
-    return Dirichlet(f, copy(faces), field_name, __to_components(components), Int[], Int[])
+    return Dirichlet(f, faces, field_name, __to_components(components), Int[], Int[])
 end
 
 # components=nothing is default and means that all components should be constrained
@@ -272,7 +272,8 @@ function add_prescribed_dof!(ch::ConstraintHandler, constrained_dof::Int, inhomo
     return ch
 end
 
-function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcfaces::Set{Index}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::Set{Int}=Set{Int}(1:getncells(ch.dh.grid))) where {Index<:BoundaryIndex}
+# 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)))
     copy!(dbc.local_face_dofs, local_face_dofs)
@@ -282,10 +283,6 @@ function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcfaces::Set{Index}, inter
     constrained_dofs = Int[]
     cc = CellCache(ch.dh, UpdateFlags(; nodes=false, coords=false, dofs=true))
     for (cellidx, faceidx) in bcfaces
-        if cellidx ∉ cellset
-            delete!(dbc.faces, Index(cellidx, faceidx))
-            continue # skip faces that are not part of the cellset
-        end
         reinit!(cc, cellidx)
         r = local_face_dofs_offset[faceidx]:(local_face_dofs_offset[faceidx+1]-1)
         append!(constrained_dofs, cc.dofs[local_face_dofs[r]]) # TODO: for-loop over r and simply push! to ch.prescribed_dofs
@@ -318,6 +315,7 @@ function _local_face_dofs_for_bc(interpolation, field_dim, components, offset, b
     return local_face_dofs, local_face_dofs_offset
 end
 
+# Dirichlet on nodeset
 function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcnodes::Set{Int}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::Set{Int}=Set{Int}(1:getncells(ch.dh.grid)))
     if interpolation !== default_interpolation(typeof(ch.dh.grid.cells[first(cellset)]))
         @warn("adding constraint to nodeset is not recommended for sub/super-parametric approximations.")
@@ -446,11 +444,11 @@ function _update!(inhomogeneities::Vector{Float64}, f::Function, boundary_entiti
 end
 
 # for nodes
-function _update!(inhomogeneities::Vector{Float64}, f::Function, nodes::Set{Int}, field::Symbol, nodeidxs::Vector{Int}, globaldofs::Vector{Int},
+function _update!(inhomogeneities::Vector{Float64}, f::Function, ::Set{Int}, field::Symbol, nodeidxs::Vector{Int}, globaldofs::Vector{Int},
                   components::Vector{Int}, dh::AbstractDofHandler, facevalues::BCValues,
                   dofmapping::Dict{Int,Int}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, time::Real) where T
     counter = 1
-    for (idx, nodenumber) in enumerate(nodeidxs)
+    for nodenumber in nodeidxs
         x = dh.grid.nodes[nodenumber].x
         bc_value = f(x, time)
         @assert length(bc_value) == length(components)
@@ -818,81 +816,63 @@ end
 Add a `Dirichlet` boundary condition to the `ConstraintHandler`.
 """
 function add!(ch::ConstraintHandler, dbc::Dirichlet)
+    # Duplicate the Dirichlet constraint for every FieldHandler
     dbc_added = false
     for fh in ch.dh.fieldhandlers
-        if !isnothing(_find_field(fh, dbc.field_name)) && _in_cellset(ch.dh.grid, fh.cellset, dbc.faces; all=false)
-            # Dofs in `dbc` not in `fh` will be removed, hence `dbc.faces` must be copied.
-            # Recreating the `dbc` will create a copy of `dbc.faces`.
-            # In this case, add! will warn, unless `warn_not_in_cellset=false`
-            dbc_ = Dirichlet(dbc.field_name, dbc.faces, dbc.f, 
-                isempty(dbc.components) ? nothing : dbc.components) 
-                # Check for empty already done when user created `dbc`
-            add!(ch, fh, dbc_, warn_not_in_cellset=false)
-            dbc_added = true
+        # Skip if the constrained field does not live on this sub domain
+        dbc.field_name in fh.field_names || continue
+        # Compute the intersection between dbc.set and the cellset of this
+        # FieldHandler and skip if the set is empty
+        filtered_set = filter_dbc_set(ch.dh.grid, fh.cellset, dbc.faces)
+        isempty(filtered_set) && continue
+        # Fetch information about the field on this FieldHandler
+        field_idx = find_field(fh, dbc.field_name)
+        interpolation = getfieldinterpolation(fh, field_idx)
+        field_dim = getfielddim(fh, field_idx)
+        # Set up components to prescribe (empty input means prescribe all components)
+        components = isempty(dbc.components) ? collect(Int, 1:field_dim) : dbc.components
+        if !all(c -> 0 < c <= field_dim, components)
+            error("components $(components) not within range of field :$(dbc.field_name) ($(field_dim) dimension(s))")
+        end
+        # Create BCValues for coordinate evalutation at dof-locations
+        EntityType = eltype(dbc.faces) # (Face|Edge|Vertex)Index
+        if EntityType <: Integer
+            # BCValues are just dummy for nodesets so set to FaceIndex
+            EntityType = FaceIndex
         end
+        CT = getcelltype(ch.dh.grid, first(fh.cellset)) # Same celltype enforced in FieldHandler constructor
+        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!(
+            ch, filtered_dbc, filtered_dbc.faces, interpolation, field_dim,
+            field_offset(fh, dbc.field_name), bcvalues, fh.cellset,
+        )
+        dbc_added = true
     end
     dbc_added || error("No overlap between dbc::Dirichlet and fields in the ConstraintHandler's DofHandler")
     return ch
 end
 
-function add!(ch::ConstraintHandler, fh::FieldHandler, dbc::Dirichlet; warn_not_in_cellset=true)
-    if warn_not_in_cellset && !(_in_cellset(ch.dh.grid, fh.cellset, dbc.faces; all=true))
-        @warn("You are trying to add a constraint a face/edge/node that is not in the cellset of the fieldhandler. This location will be skipped")
-    end
-
-    celltype = getcelltype(ch.dh.grid, first(fh.cellset)) #Assume same celltype of all cells in fh.cellset
-
-    # Extract stuff for the field
-    field_idx = find_field(fh, dbc.field_name)
-    interpolation = getfieldinterpolation(fh, field_idx)
-    field_dim = getfielddim(fh, field_idx)
-
-    if !all(c -> 0 < c <= field_dim, dbc.components)
-        error("components $(dbc.components) not within range of field :$(dbc.field_name) ($(field_dim) dimension(s))")
-    end
-
-    # Empty components means constrain them all
-    isempty(dbc.components) && append!(dbc.components, 1:field_dim)
-
-    if eltype(dbc.faces)==Int #Special case when dbc.faces is a nodeset
-        bcvalue = BCValues(interpolation, default_interpolation(celltype), FaceIndex) #Not used by node bcs, but still have to pass it as an argument
-    else
-        bcvalue = BCValues(interpolation, default_interpolation(celltype), eltype(dbc.faces))
-    end
-
-    Ferrite._add!(ch, dbc, dbc.faces, interpolation, field_dim, field_offset(fh, dbc.field_name), bcvalue, fh.cellset)
-    return ch
-end
-
-# If all==true, return true only if all items in faceset/nodeset are in the cellset
-# If all==false, return true if some items in faceset/nodeset are in the cellset
-function _in_cellset(::AbstractGrid, cellset::Set{Int}, faceset::Set{<:BoundaryIndex}; all=true)
-    for (cellid,faceid) in faceset
-        if cellid in cellset
-            all || return true
-        else
-            all && return false
-        end
+# Return the intersection of the FieldHandler set and the Dirichlet BC set
+function filter_dbc_set(::AbstractGrid, fhset::AbstractSet{Int}, dbcset::AbstractSet{<:BoundaryIndex})
+    ret = empty(dbcset)::typeof(dbcset)
+    for x in dbcset
+        cellid, _ = x
+        cellid in fhset && push!(ret, x)
     end
-    return all # if not returned by now and all==false, then no `cellid`s where in cellset
+    return ret
 end
-
-function _in_cellset(grid::AbstractGrid, cellset::Set{Int}, nodeset::Set{Int}; all=true)
-    nodes = Set{Int}()
-    for cellid in cellset
-        for nodeid in grid.cells[cellid].nodes
-            nodeid ∈ nodes || push!(nodes, nodeid)
-        end
+function filter_dbc_set(grid::AbstractGrid, fhset::AbstractSet{Int}, dbcset::AbstractSet{Int})
+    ret = empty(dbcset)
+    nodes_in_fhset = Set{Int}()
+    for cc in CellIterator(grid, fhset, UpdateFlags(; nodes=true, coords=false))
+        union!(nodes_in_fhset, cc.nodes)
     end
-
-    for nodeid in nodeset
-        if nodeid ∈ nodes
-            all || return true 
-        else
-            all && return false
-        end
+    for nodeid in dbcset
+        nodeid in nodes_in_fhset && push!(ret, nodeid)
     end
-    return all # if not returned by now and all==false, then no `cellid`s where in cellset
+    return ret
 end
 
 struct PeriodicFacePair

src/deprecations.jl

@@ -8,3 +8,5 @@ import Base: push!
 @deprecate faces(ip::Interpolation) facedof_indices(ip) false
 @deprecate edges(ip::Interpolation) edgedof_indices(ip) false
 @deprecate nfields(dh::AbstractDofHandler) length(getfieldnames(dh)) false
+
+@deprecate add!(ch::ConstraintHandler, fh::FieldHandler, dbc::Dirichlet) add!(ch, dbc)

test/test_constraints.jl

@@ -111,9 +111,8 @@ end
    close!(dh)
 
    ch = ConstraintHandler(dh)
-   add!(ch, dh.fieldhandlers[1], Dirichlet(:u, getnodeset(mesh, "bottom"), (x,t)->1.0, 1))
-   add!(ch, dh.fieldhandlers[2], Dirichlet(:u, getnodeset(mesh, "bottom"), (x,t)->1.0, 1))
-   add!(ch, dh.fieldhandlers[1], Dirichlet(:c, getnodeset(mesh, "bottom"), (x,t)->2.0, 1))
+   add!(ch, Dirichlet(:u, getnodeset(mesh, "bottom"), (x,t)->1.0, 1))
+   add!(ch, Dirichlet(:c, getnodeset(mesh, "bottom"), (x,t)->2.0, 1))
    close!(ch)
    update!(ch)
 
@@ -136,9 +135,8 @@ end
    close!(dh)
 
    ch = ConstraintHandler(dh)
-   add!(ch, dh.fieldhandlers[1], Dirichlet(:u, getfaceset(mesh, "bottom"), (x,t)->1.0, 1))
-   add!(ch, dh.fieldhandlers[2], Dirichlet(:u, getfaceset(mesh, "bottom"), (x,t)->1.0, 1))
-   add!(ch, dh.fieldhandlers[2], Dirichlet(:c, getfaceset(mesh, "bottom"), (x,t)->2.0, 1))
+   add!(ch, Dirichlet(:u, getfaceset(mesh, "bottom"), (x,t)->1.0, 1))
+   add!(ch, Dirichlet(:c, getfaceset(mesh, "bottom"), (x,t)->2.0, 1))
    close!(ch)
    update!(ch)
 

test/test_mixeddofhandler.jl

@@ -307,8 +307,8 @@ function test_2_element_heat_eq()
     ∂Ω2 = getfaceset(grid, "right")
     dbc1 = Dirichlet(:u, ∂Ω1, (x, t) -> 0)
     dbc2 = Dirichlet(:u, ∂Ω2, (x, t) -> 0)
-    add!(ch, dh.fieldhandlers[1], dbc1);
-    add!(ch, dh.fieldhandlers[2], dbc2);
+    add!(ch, dbc1);
+    add!(ch, dbc2);
     close!(ch)
 
     function doassemble(cellset, cellvalues, assembler, dh)