Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix dirichlet check cellcheck #594

Merged
merged 11 commits into from
Feb 8, 2023
43 changes: 22 additions & 21 deletions src/Dofs/ConstraintHandler.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand All @@ -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

"""
Expand Down
20 changes: 11 additions & 9 deletions test/test_constraints.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand Down