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.

Sign up for GitHub

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

Jump to bottom

Sparsity patterns without constrained dofs #539

Merged
merged 1 commit into from
Dec 20, 2022
Merged

Sparsity patterns without constrained dofs #539

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 18 additions & 8 deletions src/Dofs/ConstraintHandler.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -660,7 +660,7 @@ end
end

# Similar to Ferrite._condense!(K, ch), but only add the non-zero entries to K (that arises from the condensation process)
function _condense_sparsity_pattern!(K::SparseMatrixCSC{T}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, dofmapping::Dict{Int,Int}) where T
function _condense_sparsity_pattern!(K::SparseMatrixCSC{T}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, dofmapping::Dict{Int,Int}, keep_constrained::Bool) where T
ndofs = size(K, 1)

# Return early if there are no non-trivial affine constraints
Expand All @@ -676,6 +676,7 @@ function _condense_sparsity_pattern!(K::SparseMatrixCSC{T}, dofcoefficients::Vec
for col in 1:ndofs
col_coeffs = coefficients_for_dof(dofmapping, dofcoefficients, col)
if col_coeffs === nothing
!keep_constrained && haskey(dofmapping, col) && continue
for ri in nzrange(K, col)
row = K.rowval[ri]
row_coeffs = coefficients_for_dof(dofmapping, dofcoefficients, row)
Expand All @@ -690,14 +691,19 @@ function _condense_sparsity_pattern!(K::SparseMatrixCSC{T}, dofcoefficients::Vec
row = K.rowval[ri]
row_coeffs = coefficients_for_dof(dofmapping, dofcoefficients, row)
if row_coeffs === nothing
!keep_constrained && haskey(dofmapping, row) && continue
for (d, _) in col_coeffs
cnt += 1
_add_or_grow(cnt, I, J, row, d)
end
else
for (d1, _) in col_coeffs, (d2, _) in row_coeffs
cnt += 1
_add_or_grow(cnt, I, J, d1, d2)
for (d1, _) in col_coeffs
!keep_constrained && haskey(dofmapping, d1) && continue
for (d2, _) in row_coeffs
!keep_constrained && haskey(dofmapping, d2) && continue
cnt += 1
_add_or_grow(cnt, I, J, d1, d2)
end
end
end
end
Expand Down Expand Up @@ -899,10 +905,14 @@ function _check_cellset_dirichlet(grid::AbstractGrid, cellset::Set{Int}, nodeset
end
end

create_symmetric_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler; coupling=nothing) =
Symmetric(_create_sparsity_pattern(dh, ch, true, coupling), :U)
create_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler; coupling=nothing) =
_create_sparsity_pattern(dh, ch, false, coupling)
function create_symmetric_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler;
keep_constrained::Bool=true, coupling=nothing)
return Symmetric(_create_sparsity_pattern(dh, ch, true, keep_constrained, coupling), :U)
end
function create_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler;
keep_constrained::Bool=true, coupling=nothing)
return _create_sparsity_pattern(dh, ch, false, keep_constrained, coupling)
end

struct PeriodicFacePair
mirror::FaceIndex
Expand Down
14 changes: 10 additions & 4 deletions src/Dofs/DofHandler.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -371,7 +371,7 @@ not. By default full coupling is assumed.
See the [Sparsity Pattern](@ref) section of the manual.
"""
function create_sparsity_pattern(dh::AbstractDofHandler; coupling=nothing)
return _create_sparsity_pattern(dh, nothing, false, coupling)
return _create_sparsity_pattern(dh, nothing, false, true, coupling)
end

"""
Expand All @@ -383,10 +383,15 @@ triangle of the matrix. Return a `Symmetric{SparseMatrixCSC}`.

See the [Sparsity Pattern](@ref) section of the manual.
"""
create_symmetric_sparsity_pattern(dh::AbstractDofHandler; coupling=nothing) = Symmetric(_create_sparsity_pattern(dh, nothing, true, coupling), :U)
function create_symmetric_sparsity_pattern(dh::AbstractDofHandler; coupling=nothing)
return Symmetric(_create_sparsity_pattern(dh, nothing, true, true, coupling), :U)
end

function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{ConstraintHandler, Nothing}=#, sym::Bool, coupling::Union{AbstractMatrix{Bool}, Nothing})
function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{ConstraintHandler, Nothing}=#, sym::Bool, keep_constrained::Bool, coupling::Union{AbstractMatrix{Bool},Nothing})
@assert isclosed(dh)
if !keep_constrained
@assert ch !== nothing && isclosed(ch)
end
ncells = getncells(dh.grid)
if coupling !== nothing
# Extend coupling to be of size (ndofs_per_cell × ndofs_per_cell)
Expand All @@ -409,6 +414,7 @@ function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{Constraint
dofi = global_dofs[i]
dofj = global_dofs[j]
sym && (dofi > dofj && continue)
!keep_constrained && (haskey(ch.dofmapping, dofi) || haskey(ch.dofmapping, dofj)) && continue
cnt += 1
if cnt > length(J)
resize!(I, trunc(Int, length(I) * 1.5))
Expand Down Expand Up @@ -438,7 +444,7 @@ function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{Constraint
if ch !== nothing
@assert isclosed(ch)
fill!(K.nzval, 1)
_condense_sparsity_pattern!(K, ch.dofcoefficients, ch.dofmapping)
_condense_sparsity_pattern!(K, ch.dofcoefficients, ch.dofmapping, keep_constrained)
fillzero!(K)
end

Expand Down
51 changes: 51 additions & 0 deletions test/test_constraints.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1375,3 +1375,54 @@ end # testset
end
end
end # testset

@testset "Sparsity pattern without constrained dofs" begin
grid = generate_grid(Triangle, (5, 5))
dh = DofHandler(grid)
push!(dh, :u, 1)
close!(dh)
ch = ConstraintHandler(dh)
add!(ch, Dirichlet(:u, getfaceset(grid, "left"), (x, t) -> 0))
close!(ch)
Kfull = create_sparsity_pattern(dh, ch)
K = create_sparsity_pattern(dh, ch; keep_constrained=false)
# Pattern tests
nonzero_edges = Set(
(i, j) for d in 1:getncells(grid)
for (i, j) in Iterators.product(celldofs(dh, d), celldofs(dh, d))
)
zero_edges = setdiff(Set(Iterators.product(1:ndofs(dh), 1:ndofs(dh))), nonzero_edges)
for (i, j) in zero_edges
@test Base.isstored(K, i, j) == Base.isstored(Kfull, i, j) == false
end
for (i, j) in nonzero_edges
if i != j && (haskey(ch.dofmapping, i) || haskey(ch.dofmapping, j))
@test !Base.isstored(K, i, j)
else
@test Base.isstored(K, i, j)
end
@test Base.isstored(Kfull, i, j)
end
# Assembly
afull = start_assemble(Kfull)
a = start_assemble(K)
Kc = copy(K)
ac = start_assemble(Kc, zeros(ndofs(dh)))
for cell in CellIterator(dh)
ke = rand(ndofs_per_cell(dh), ndofs_per_cell(dh))
assemble!(afull, celldofs(cell), ke)
kec = copy(ke)
apply_local!(kec, rand(ndofs_per_cell(dh)), celldofs(cell), ch)
assemble!(a, celldofs(cell), kec)
apply_assemble!(ac, ch, celldofs(cell), ke, rand(ndofs_per_cell(dh)))
end
apply!(Kfull, ch)
fd = free_dofs(ch)
pd = ch.prescribed_dofs
@test K ≈ Kc
@test K[fd, fd] ≈ Kc[fd, fd] ≈ Kfull[fd, fd]
for d in ch.prescribed_dofs
K[d, d] = Kc[d, d] = Kfull[d, d] = 1
end
@test K ≈ Kc ≈ Kfull
end # testset