Fix Dirichlet cellset check, fixes #593 (#594)

src/Dofs/ConstraintHandler.jl

@@ -915,32 +915,25 @@ end
 function add!(ch::ConstraintHandler{<:MixedDofHandler}, dbc::Dirichlet)
     dbc_added = false
     for fh in ch.dh.fieldhandlers
-        if overlaps(fh, dbc)
-            # Recreating the `dbc` will create a copy of `dbc.faces`. 
-            # If `dbc` have dofs not in `fh`, these will be removed from `dbc`, 
-            # thus the need to copy `dbc.faces`.
-            # In this case, add! will warn, unless warn_check_cellset=false
+        if dbc.field_name in getfieldnames(fh) && _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_check_cellset=false)
+            add!(ch, fh, dbc_, warn_not_in_cellset=false)
             dbc_added = true
         end
     end
     dbc_added || @warn("No overlap between dbc::Dirichlet and fields in the ConstraintHandler's MixedDofHandler")
     return ch
 end
 
-function overlaps(fh::FieldHandler, dbc::Dirichlet)
-    dbc.field_name in getfieldnames(fh) || return false # Must contain the correct field
-    for (cellid, _) in dbc.faces
-        cellid in fh.cellset && return true
+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
-    return false
-end
-
-function add!(ch::ConstraintHandler, fh::FieldHandler, dbc::Dirichlet; warn_check_cellset=true)
-    warn_check_cellset && _check_cellset_dirichlet(ch.dh.grid, fh.cellset, dbc.faces)
 
     celltype = getcelltype(ch.dh.grid, first(fh.cellset)) #Assume same celltype of all cells in fh.cellset
 
@@ -966,15 +959,20 @@ function add!(ch::ConstraintHandler, fh::FieldHandler, dbc::Dirichlet; warn_chec
     return ch
 end
 
-function _check_cellset_dirichlet(::AbstractGrid, cellset::Set{Int}, faceset::Set{<:BoundaryIndex})
+# 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)
-            @warn("You are trying to add a constraint to a face that is not in the cellset of the fieldhandler. The face will be skipped.")
+        if cellid in cellset
+            all || return true
+        else
+            all && return false
         end
     end
+    return all # if not returned by now and all==false, then no `cellid`s where in cellset
 end
 
-function _check_cellset_dirichlet(grid::AbstractGrid, cellset::Set{Int}, nodeset::Set{Int})
+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
@@ -983,10 +981,13 @@ function _check_cellset_dirichlet(grid::AbstractGrid, cellset::Set{Int}, nodeset
     end
 
     for nodeid in nodeset
-        if !(nodeid ∈ nodes)
-            @warn("You are trying to add a constraint to a node that is not in the cellset of the fieldhandler. The node will be skipped.")
+        if nodeid ∈ nodes
+            all || return true 
+        else
+            all && return false
         end
     end
+    return all # if not returned by now and all==false, then no `cellid`s where in cellset
 end
 
 """

test/test_constraints.jl

@@ -227,10 +227,10 @@ end
     field_vQ = Field(:v, ip_quad, 1)
 
     # Order important for test to ensure consistent dof ordering
-    push!(dh, FieldHandler([field_uQ, field_vQ], getcellset(grid, "uandvQ")))
-    push!(dh, FieldHandler([field_uT, field_vT], getcellset(grid, "uandvT")))
-    push!(dh, FieldHandler([field_uT], getcellset(grid, "onlyuT")))
-    push!(dh, FieldHandler([field_uQ], getcellset(grid, "onlyuQ")))
+    add!(dh, FieldHandler([field_uQ, field_vQ], getcellset(grid, "uandvQ")))
+    add!(dh, FieldHandler([field_uT, field_vT], getcellset(grid, "uandvT")))
+    add!(dh, FieldHandler([field_uT], getcellset(grid, "onlyuT")))
+    add!(dh, FieldHandler([field_uQ], getcellset(grid, "onlyuQ")))
     close!(dh)
 
     # Add constraints 
@@ -240,23 +240,25 @@ end
     dB_u = Dirichlet(:u, getfaceset(grid, "B"), (x,t) -> 3.0)  # Note, overwrites dA_u on node 3 
     dB_v = Dirichlet(:v, getfaceset(grid, "B"), (x,t) -> 4.0)  # :v not on cells with "B"-faces
     dC_v = Dirichlet(:v, getfaceset(grid, "C"), (x,t) -> 5.0)  # :v not on cells with "C"-faces
+    dN_u = Dirichlet(:u, Set(10), (x,t) -> 6.0)                # Add on node 10
 
     @test_logs min_level=Logging.Warn add!(ch, dA_u)    # No warning should be issued
     @test_logs min_level=Logging.Warn add!(ch, dA_v)    # No warning should be issued
     @test_logs min_level=Logging.Warn add!(ch, dB_u)    # No warning should be issued
     @test_logs (:warn,) add!(ch, dB_v)  # Warn about :v not in cells connected with dB_v's faceset
     @test_logs (:warn,) add!(ch, dC_v)  # Warn about :v not in cells connected with dC_v's faceset
+    @test_logs min_level=Logging.Warn add!(ch, dN_u)    # No warning should be issued (add to node)
     close!(ch)
 
     # The full bottom part of the mesh has been prescribed
-    @test sort(ch.prescribed_dofs) == sort([nd[i] for nd in nodedofs[1:5] for i in 1:length(nd)])
+    @test sort(ch.prescribed_dofs) == sort(push!([nd[i] for nd in nodedofs[1:5] for i in 1:length(nd)], 15))
 
     # Test that the correct dofs have been prescribed
     update!(ch, 0.0)
-    #                 nodes       N1,  N2,  N1,  N2,  N3,  N4,  N5
-    #                 field       :u,  :u,  :v,  :v,  :u,  :u,  :u
-    #                   dof        1,   2,   5,   6,  11,  12,  14
-    @test ch.inhomogeneities == [1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 3.0]
+    #                 nodes       N1,  N2,  N1,  N2,  N3,  N4,  N5   N10
+    #                 field       :u,  :u,  :v,  :v,  :u,  :u,  :u   :u
+    #                   dof        1,   2,   5,   6,  11,  12,  14   15
+    @test ch.inhomogeneities == [1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 3.0, 6.0]
     # Note that dB_u overwrite dA_u @ N3, hence the value 3.0 there
 end