Sparsity patterns without constrained dofs

This patch adds a new keyword argument `keep_constrained::Bool=true` to
`create_sparsity_pattern` which, when set to `false`, do not create
entries in the global matrix for rows/values that correspond to
constrained degrees of freedom. Naturally this can only be used with
local condensation of constraints (cf. #528).

src/Dofs/ConstraintHandler.jl

@@ -669,7 +669,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
@@ -685,6 +685,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)
@@ -701,17 +702,22 @@ 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
                         if !_addindex_sparsematrix!(K, 0.0, row, d)
                             cnt += 1
                             _add_or_grow(cnt, I, J, row, d)
                         end
                     end
                 else
-                    for (d1, _) in col_coeffs, (d2, _) in row_coeffs
-                        if !_addindex_sparsematrix!(K, 0.0, d1, d2)
-                            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
+                            if !_addindex_sparsematrix!(K, 0.0, d1, d2)
+                                cnt += 1
+                                _add_or_grow(cnt, I, J, d1, d2)
+                            end
                         end
                     end
                 end
@@ -934,10 +940,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

src/Dofs/DofHandler.jl

@@ -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
 
 """
@@ -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)
@@ -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))
@@ -438,7 +444,7 @@ function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{Constraint
 
         V = ones(length(I))
         K = sparse(I, J, V, ndofs(dh), ndofs(dh))
-        _condense_sparsity_pattern!(K, ch.dofcoefficients, ch.dofmapping)
+        _condense_sparsity_pattern!(K, ch.dofcoefficients, ch.dofmapping, keep_constrained)
         fill!(K.nzval, 0.0)
     else
         V = zeros(length(I))

test/test_constraints.jl

@@ -1375,3 +1375,71 @@ 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
+    function is_stored(A, i, j)
+        for m in nzrange(A, j)
+            A.rowval[m] == i && return true
+        end
+        return false
+    end
+    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 is_stored(K, i, j) == is_stored(Kfull, i, j) == false
+    end
+    for (i, j) in nonzero_edges
+        if i != j && (haskey(ch.dofmapping, i) || haskey(ch.dofmapping, j))
+            @test !is_stored(K, i, j)
+        else
+            @test is_stored(K, i, j)
+        end
+        @test is_stored(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
+
+    # Error paths
+    K = spdiagm(0 => zeros(2))
+    a = start_assemble(K)
+    err = ErrorException("some row indices were not found")
+    ## Errors below diagonal
+    @test_throws err assemble!(a, [1, 2], [1.0 0.0; 3.0 4.0])
+    @test_throws err assemble!(a, [2, 1], [1.0 2.0; 0.0 4.0])
+    ## Errors above diagonal
+    @test_throws err assemble!(a, [1, 2], [1.0 2.0; 0.0 4.0])
+    @test_throws err assemble!(a, [2, 1], [1.0 0.0; 3.0 4.0])
+end # testset