Do not assemble zero contributions from local to global matrix

With this patch, entries from the local matrix that are zero are ignored
and never added to the global matrix. More importantly, it is now
allowed to have missing stored entries in the global matrix if the
contribution that would have been added there is zero.

A concrete example where this is useful: Consider a two-field problem
with unknowns (u, p) and test functions (v, q) where there is no
coupling between, for example, p and q. In the global block matrix K =
[A B; C D], D will be zero, and there is therefore no need to allocate
entries in K corresponding to all couplings between p- and q-dofs. (Note
though that currently it is not possible to create such a sparsity
pattern because create_sparsity_pattern assumes full coupling between
fields.) However, when assembling the local matrix k = [a b; c d] it is
much easier to use a dense matrix, and simply initialize d to zero.
After this patch it is perfectly fine to assemble the matrix k into K,
even if the global entries for d are missing.

src/assembler.jl

@@ -82,7 +82,10 @@ Assembles the element residual `ge` into the global residual vector `g`.
 @propagate_inbounds function assemble!(g::AbstractVector{T}, dofs::AbstractVector{Int}, ge::AbstractVector{T}) where {T}
     @boundscheck checkbounds(g, dofs)
     @inbounds for i in 1:length(dofs)
-        g[dofs[i]] += ge[i]
+        val = ge[i]
+        if !iszero(val)
+            g[dofs[i]] += val
+        end
     end
 end
 
@@ -181,18 +184,38 @@ end
     current_col = 1
     @inbounds for Kcol in sorteddofs
         maxlookups = sym ? current_col : ld
-        current_idx = 1
-        for r in nzrange(K, Kcol)
-            Kerow = permutation[current_idx]
-            if K.rowval[r] == dofs[Kerow]
-                Kecol = permutation[current_col]
-                K.nzval[r] += Ke[Kerow, Kecol]
-                current_idx += 1
+        Kecol = permutation[current_col]
+        ri = 1 # row index pointer for the local matrix
+        Ri = 1 # row index pointer for the global matrix
+        nzr = nzrange(K, Kcol)
+        while Ri <= length(nzr) && ri <= maxlookups
+            R = nzr[Ri]
+            Krow = K.rowval[R]
+            Kerow = permutation[ri]
+            val = Ke[Kerow, Kecol]
+            if Krow == dofs[Kerow]
+                # Match: add the value (if non-zero) and advance the pointers
+                if !iszero(val)
+                    K.nzval[R] += val
+                end
+                ri += 1
+                Ri += 1
+            elseif Krow < dofs[Kerow]
+                # No match yet: advance the global matrix row pointer
+                Ri += 1
+            else # Krow > dofs[Kerow]
+                # No match: no entry exist in the global matrix for this row. This is
+                # allowed as long as the value which would have been inserted is zero.
+                iszero(val) || error("some row indices were not found")
+                # Advance the local matrix row pointer
+                ri += 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 are all zero
+        for i in ri:maxlookups
+            if !iszero(Ke[permutation[i], Kecol])
+                error("some row indices were not found")
+            end
         end
         current_col += 1
     end

test/test_assemble.jl

@@ -75,3 +75,53 @@
     @test_throws AssertionError assemble!(assembler, [11, 1, 2], rand(4, 4))
     @test_throws AssertionError assemble!(assembler, [11, 1, 2, 3], rand(4, 4), rand(3))
 end
+
+struct IgnoreMeIfZero
+    x::Float64
+end
+Base.iszero(x::IgnoreMeIfZero) = iszero(x.x)
+function Base.:+(y::Float64, x::IgnoreMeIfZero)
+    @test !iszero(x.x)
+    return y + x.x
+end
+
+@testset "assemble! ignoring zeros" begin
+    store_dofs    = [1, 5, 2, 8]
+    assemble_dofs = [1, 5, 4, 8]
+    I = repeat(store_dofs; outer=4)
+    J = repeat(store_dofs; inner=4)
+    V = zeros(length(I))
+    K = sparse(I, J, V)
+    D = zeros(size(K))
+
+    # Standard assembler
+    a = start_assemble(K)
+    ke = rand(4,4); ke[3, :] .= 0; ke[:, 3] .= 0; ke[2,2] = 0
+    assemble!(a, assemble_dofs, IgnoreMeIfZero.(ke))
+    D[assemble_dofs, assemble_dofs] += ke
+    @test K == D
+
+    # Symmetric assembler
+    S = Symmetric(K)
+    assembler = start_assemble(S)
+    fill!(D, 0)
+    kes = [(ke[i, j] + ke[j, i]) / 2 for i in 1:4, j in 1:4]
+    D[assemble_dofs, assemble_dofs] += kes
+    kes[2, 1] = 42 # To check we don't touch elements below diagonal
+    assemble!(assembler, assemble_dofs, IgnoreMeIfZero.(kes))
+    @test S == D
+
+    # Error paths
+    K = spdiagm(0 => zeros(2))
+    a = start_assemble(K)
+    as = start_assemble(Symmetric(K))
+    errr = ErrorException("some row indices were not found")
+    ## Errors below diagonal
+    @test_throws errr assemble!(a, [1, 2], [1.0 0.0; 3.0 4.0])
+    @test_throws errr assemble!(a, [2, 1], [1.0 2.0; 0.0 4.0])
+    ## Errors above diagonal
+    @test_throws errr assemble!(a, [1, 2], [1.0 2.0; 0.0 4.0])
+    @test_throws errr assemble!(as, [1, 2], [1.0 2.0; 0.0 4.0])
+    @test_throws errr assemble!(a, [2, 1], [1.0 0.0; 3.0 4.0])
+    @test_throws errr assemble!(as, [2, 1], [1.0 0.0; 3.0 4.0])
+end