Enable assembly of rect blocks into nxn CSR/CSC matrix, as suggested by Knut.

Author: termi-official

src/assembler.jl

@@ -245,13 +245,7 @@ start_assemble(K::Union{AbstractSparseMatrixCSC, Symmetric{<:Any, <:AbstractSpar
 
 function start_assemble(K::AbstractSparseMatrixCSC{T}, f::Vector = T[]; fillzero::Bool = true, maxcelldofs_hint::Int = 0) where {T}
     fillzero && (fillzero!(K); fillzero!(f))
-    if size(K, 1) == size(K, 2)
-        permutation = zeros(Int, maxcelldofs_hint)
-        sorteddofs = zeros(Int, maxcelldofs_hint)
-        return CSCAssembler(K, f, permutation, permutation, sorteddofs, sorteddofs)
-    else
-        return CSCAssembler(K, f, zeros(Int, maxcelldofs_hint), zeros(Int, maxcelldofs_hint), zeros(Int, maxcelldofs_hint), zeros(Int, maxcelldofs_hint))
-    end
+    return CSCAssembler(K, f, zeros(Int, maxcelldofs_hint), zeros(Int, maxcelldofs_hint), zeros(Int, maxcelldofs_hint), zeros(Int, maxcelldofs_hint))
 end
 function start_assemble(K::Symmetric{T, <:SparseMatrixCSC}, f::Vector = T[]; fillzero::Bool = true, maxcelldofs_hint::Int = 0) where {T}
     fillzero && (fillzero!(K); fillzero!(f))
@@ -322,7 +316,7 @@ end
     # a specific order, which might not be the sorted order. Hence we sort them.
     # Note that we are not allowed to mutate `dofs` in the process.
     sortedcoldofs, colpermutation = _sortdofs_for_assembly!(A.colpermutation, A.sortedcoldofs, coldofs)
-    sortedrowdofs, rowpermutation = if A.sortedcoldofs !== A.sortedrowdofs
+    sortedrowdofs, rowpermutation = if rowdofs !== coldofs
         _sortdofs_for_assembly!(A.rowpermutation, A.sortedrowdofs, rowdofs)
     else
         sortedcoldofs, colpermutation

test/test_assemble.jl

@@ -68,7 +68,20 @@ import LinearAlgebra: Symmetric
     fe = rand(length(rdofs))
     assemble!(assembler, rdofs, cdofs, Ke, fe)
     assemble!(assembler, rdofs, cdofs, Ke, fe)
-    @test_throws ArgumentError assemble!(assembler, rdofs, Ke, fe)
+    @test_throws ArgumentError assemble!(assembler, rdofs, Ke, fe) # Not in sparsity pattern
+    @test all(K[rdofs, cdofs] .== 2Ke)
+    @test all(f[rdofs] .== 2fe)
+
+    # CSCAssembler: Assemble rectangular part in quadratic matrix
+    K = SparseMatrixCSC(6, 6, [K.colptr..., 7, 7, 7], K.rowval, K.nzval)
+    assembler = start_assemble(K, f)
+    rdofs = [1, 4, 6]
+    cdofs = [1, 3]
+    Ke = rand(length(rdofs), length(cdofs))
+    fe = rand(length(rdofs))
+    assemble!(assembler, rdofs, cdofs, Ke, fe)
+    assemble!(assembler, rdofs, cdofs, Ke, fe)
+    @test_throws ArgumentError assemble!(assembler, rdofs, Ke, fe) # Not in sparsity pattern
     @test all(K[rdofs, cdofs] .== 2Ke)
     @test all(f[rdofs] .== 2fe)
 

test/test_assembler_extensions.jl

@@ -85,16 +85,35 @@ using SparseArrays, LinearAlgebra
         @test K ≈ sparsecsr(I, J, V)
         @test f ≈ [4 / 3, 2.0, 1.0]
 
+        # CSRAssembler: assemble with different row and col dofs
         I = [1, 1, 4, 4, 6, 6]
         J = [1, 3, 1, 3, 1, 3]
         V = zeros(length(I))
         K = sparsecsr(I, J, V)
-        assembler = start_assemble(K)
+        f = zeros(6)
+        assembler = start_assemble(K, f)
+        rdofs = [1, 4, 6]
+        cdofs = [1, 3]
+        Ke = rand(length(rdofs), length(cdofs))
+        fe = rand(length(rdofs))
+        assemble!(assembler, rdofs, cdofs, Ke, fe)
+        assemble!(assembler, rdofs, cdofs, Ke, fe)
+        @test_throws ArgumentError assemble!(assembler, rdofs, Ke, fe) # Not in sparsity pattern
+        @test all(K[rdofs, cdofs] .== 2Ke)
+        @test all(f[rdofs] .== 2fe)
+
+        # CSRAssembler: Assemble rectangular part in quadratic matrix
+        K = SparseMatrixCSR(6, 6, [K.colptr..., 7, 7, 7], K.rowval, K.nzval)
+        assembler = start_assemble(K, f)
         rdofs = [1, 4, 6]
         cdofs = [1, 3]
         Ke = rand(length(rdofs), length(cdofs))
-        assemble!(assembler, rdofs, cdofs, Ke)
-        @test (K[rdofs, cdofs] .== Ke) |> all
+        fe = rand(length(rdofs))
+        assemble!(assembler, rdofs, cdofs, Ke, fe)
+        assemble!(assembler, rdofs, cdofs, Ke, fe)
+        @test_throws ArgumentError assemble!(assembler, rdofs, Ke, fe) # Not in sparsity pattern
+        @test all(K[rdofs, cdofs] .== 2Ke)
+        @test all(f[rdofs] .== 2fe)
 
         # Check if coupling works
         grid = generate_grid(Quadrilateral, (2, 2))