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).

Another necessary change included in this patch is that assembly now
ignores entries in the local matrix that are zero (in particular
rows/columns in the local matrix which have been zeroed out in the local
condensation).

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,8 +940,14 @@ function _check_cellset_dirichlet(grid::AbstractGrid, cellset::Set{Int}, nodeset
     end
 end
 
-create_symmetric_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler) = Symmetric(_create_sparsity_pattern(dh, ch, true), :U)
-create_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler) = _create_sparsity_pattern(dh, ch, false)
+function create_symmetric_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler;
+                                           keep_constrained::Bool=true)
+    return Symmetric(_create_sparsity_pattern(dh, ch, true, keep_constrained), :U)
+end
+function create_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler;
+                                 keep_constrained::Bool=true)
+    return _create_sparsity_pattern(dh, ch, false, keep_constrained)
+end
 
 struct PeriodicFacePair
     mirror::FaceIndex

src/Dofs/DofHandler.jl

@@ -333,7 +333,7 @@ with stored values in the correct places.
 
 See the [Sparsity Pattern](@ref) section of the manual.
 """
-create_sparsity_pattern(dh::AbstractDofHandler) = _create_sparsity_pattern(dh, nothing, false)
+create_sparsity_pattern(dh::AbstractDofHandler) = _create_sparsity_pattern(dh, nothing, false, true)
 
 """
     create_symmetric_sparsity_pattern(dh::DofHandler)
@@ -344,10 +344,13 @@ triangle of the matrix. Return a `Symmetric{SparseMatrixCSC}`.
 
 See the [Sparsity Pattern](@ref) section of the manual.
 """
-create_symmetric_sparsity_pattern(dh::AbstractDofHandler) = Symmetric(_create_sparsity_pattern(dh, nothing, true), :U)
+create_symmetric_sparsity_pattern(dh::AbstractDofHandler) = Symmetric(_create_sparsity_pattern(dh, nothing, true, true), :U)
 
-function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{ConstraintHandler, Nothing}=#, sym::Bool)
+function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{ConstraintHandler, Nothing}=#, sym::Bool, keep_constrained::Bool)
     @assert isclosed(dh)
+    if !keep_constrained
+        @assert ch !== nothing && isclosed(ch)
+    end
     ncells = getncells(dh.grid)
     # Compute approximate size for the buffers using the dofs in the first element
     n = ndofs_per_cell(dh)
@@ -365,6 +368,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))
@@ -394,7 +398,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))

src/assembler.jl

@@ -181,18 +181,38 @@ end
     current_col = 1
     @inbounds for Kcol in sorteddofs
         maxlookups = sym ? current_col : ld
+        Kecol = permutation[current_col]
         current_idx = 1
-        for r in nzrange(K, Kcol)
+        ri = 1
+        nzr = nzrange(K, Kcol)
+        while ri <= length(nzr) && current_idx <= maxlookups
+            r = nzr[ri]
             Kerow = permutation[current_idx]
+            val = Ke[Kerow, Kecol]
             if K.rowval[r] == dofs[Kerow]
-                Kecol = permutation[current_col]
-                K.nzval[r] += Ke[Kerow, Kecol]
+                # Match: add the value
+                if !iszero(val)
+                    K.nzval[r] += val
+                end
+                ri += 1
+                current_idx += 1
+            elseif K.rowval[r] < dofs[Kerow]
+                # No match yet: advance the global row pointer
+                ri += 1
+            elseif K.rowval[r] > dofs[Kerow]
+                # No match: no entry exist in the global matrix for this dof. This is
+                # allowed as long as the value which would have been inserted is zero.
+                if !iszero(val)
+                    error("some row indices were not found")
+                end
                 current_idx += 1
             end
-            current_idx > maxlookups && break
         end
-        if current_idx <= maxlookups
-            error("some row indices were not found")
+        # Make sure that remaining entries in this column of the local matrix is zero
+        for i in current_idx:maxlookups
+            if !iszero(Ke[permutation[i], Kecol])
+                error("some row indices were not found")
+            end
         end
         current_col += 1
     end

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