Modularize assembly internals and add SparseMatrixCSR extension

CHANGELOG.md

@@ -5,6 +5,14 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [Next] - xxxx-xx-xx
+
+### Added
+- Support for directly assembling to `SparseMatrixCSR` (from `SparseMatricesCSR.jl`). ([#864])
+
+### Documentation updates
+- Extended assembly docs with information on how to support direct assembly into new matrix types. ([#864])
+
 ## [v1.1.0] - 2025-05-01
 
 ### Added

Project.toml

@@ -18,17 +18,20 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 [weakdeps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
+SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
 
 [extensions]
 FerriteBlockArrays = "BlockArrays"
 FerriteMetis = "Metis"
+FerriteSparseMatrixCSR = "SparseMatricesCSR"
 
 [compat]
 BlockArrays = "0.16, 1"
 EnumX = "1"
 HCubature = "1.7.0"
 ForwardDiff = "0.10, 1"
 Metis = "1.3"
+SparseMatricesCSR = "0.6"
 NearestNeighbors = "0.4"
 OrderedCollections = "1"
 Preferences = "1"
@@ -53,9 +56,10 @@ ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
 TaskLocalValues = "ed4db957-447d-4319-bfb6-7fa9ae7ecf34"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [targets]
-test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "NBInclude", "OhMyThreads", "ProgressMeter", "Random", "SHA", "TaskLocalValues", "Test", "TimerOutputs", "Logging", "HCubature"]
+test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "NBInclude", "OhMyThreads", "ProgressMeter", "Random", "SHA", "SparseMatricesCSR", "TaskLocalValues", "Test", "TimerOutputs", "Logging", "HCubature"]

docs/src/devdocs/assembly.md

@@ -1,10 +1,48 @@
 # [Assembly](@id devdocs-assembly)
 
+The assembler handles the insertion of the element matrices and element vectors into the system matrix and right hand side.
+
+## Custom matrix formats
+While the CSC and CSR formats are the most common sparse matrix formats in practice, users might want to have optimized custom matrix formats for their specific use-case. The default assemblers [`Ferrite.CSCAssembler`](@ref) and [`Ferrite.CSRAssembler`](@ref) should be able to handle most cases in practice. To support a custom format users have to dispatch the following functions on their matrix type. There is the public interface
+
+```@docs; canonical=false
+Ferrite.allocate_matrix
+```
+
+the internal interface
+```@docs
+Ferrite.zero_out_rows!
+Ferrite.zero_out_columns!
+Ferrite._condense!
+```
+
+and the `AbstractSparseMatrix` interface for their custom matrix type. Optional dispatches to speed up operations might be
+
+```@docs
+Ferrite.add_inhomogeneities!
+```
+
+## Custom Assembler
+
+In case the default assembler is insufficient, users can implement a custom assemblers. For this, they can create a custom type and dispatch the following functions.
+
+```@docs; canonical=false
+start_assemble
+assemble!
+```
+
+For local elimination support the following functions might also need custom dispatches
+
+```@docs
+Ferrite._condense_local!
+```
+
 ## Type definitions
 
 ```@docs
 Ferrite.COOAssembler
 Ferrite.CSCAssembler
+Ferrite.CSRAssembler
 Ferrite.SymmetricCSCAssembler
 ```
 

ext/FerriteSparseMatrixCSR.jl

@@ -0,0 +1,110 @@
+module FerriteSparseMatrixCSR
+
+using Ferrite, SparseArrays, SparseMatricesCSR
+import Ferrite: AbstractSparsityPattern, CSRAssembler
+import Base: @propagate_inbounds
+
+# Could be generalized if https://github.com/JuliaSparse/SparseArrays.jl/pull/546 is merged
+function Ferrite.start_assemble(K::SparseMatrixCSR{<:Any, T}, f::Vector = T[]; fillzero::Bool = true, maxcelldofs_hint::Int = 0) where {T}
+    fillzero && (Ferrite.fillzero!(K); Ferrite.fillzero!(f))
+    return CSRAssembler(K, f, zeros(Int, maxcelldofs_hint), zeros(Int, maxcelldofs_hint))
+end
+
+@propagate_inbounds function Ferrite._assemble_inner!(K::SparseMatrixCSR, Ke::AbstractMatrix, dofs::AbstractVector, sorteddofs::AbstractVector, permutation::AbstractVector, sym::Bool)
+    current_row = 1
+    ld = length(dofs)
+    return @inbounds for Krow in sorteddofs
+        maxlookups = sym ? current_row : ld
+        Kerow = permutation[current_row]
+        ci = 1 # col index pointer for the local matrix
+        Ci = 1 # col index pointer for the global matrix
+        nzr = nzrange(K, Krow)
+        while Ci <= length(nzr) && ci <= maxlookups
+            C = nzr[Ci]
+            Kcol = K.colval[C]
+            Kecol = permutation[ci]
+            val = Ke[Kerow, Kecol]
+            if Kcol == dofs[Kecol]
+                # Match: add the value (if non-zero) and advance the pointers
+                if !iszero(val)
+                    K.nzval[C] += val
+                end
+                ci += 1
+                Ci += 1
+            elseif Kcol < dofs[Kecol]
+                # No match yet: advance the global matrix row pointer
+                Ci += 1
+            else # Kcol > dofs[Kecol]
+                # 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) || Ferrite._missing_sparsity_pattern_error(Krow, Kcol)
+                # Advance the local matrix row pointer
+                ci += 1
+            end
+        end
+        # Make sure that remaining entries in this column of the local matrix are all zero
+        for i in ci:maxlookups
+            if !iszero(Ke[Kerow, permutation[i]])
+                Ferrite._missing_sparsity_pattern_error(Krow, sorteddofs[i])
+            end
+        end
+        current_row += 1
+    end
+end
+
+function Ferrite.zero_out_rows!(K::SparseMatrixCSR, ch::ConstraintHandler)
+    @debug @assert issorted(ch.prescribed_dofs)
+    for row in ch.prescribed_dofs
+        r = nzrange(K, row)
+        K.nzval[r] .= 0.0
+    end
+    return
+end
+
+function Ferrite.zero_out_columns!(K::SparseMatrixCSR, ch::ConstraintHandler)
+    colval = K.colval
+    nzval = K.nzval
+    return @inbounds for i in eachindex(colval, nzval)
+        if haskey(ch.dofmapping, colval[i])
+            nzval[i] = 0
+        end
+    end
+end
+
+function Ferrite.allocate_matrix(::Type{SparseMatrixCSR}, sp::AbstractSparsityPattern)
+    return Ferrite.allocate_matrix(SparseMatrixCSR{1, Float64, Int64}, sp)
+end
+
+function Ferrite.allocate_matrix(::Type{SparseMatrixCSR{1, Tv, Ti}}, sp::AbstractSparsityPattern) where {Tv, Ti}
+    return _allocate_matrix(SparseMatrixCSR{1, Tv, Ti}, sp, false)
+end
+
+function _allocate_matrix(::Type{SparseMatrixCSR{1, Tv, Ti}}, sp::AbstractSparsityPattern, sym::Bool) where {Tv, Ti}
+    # 1. Setup rowptr
+    rowptr = zeros(Ti, Ferrite.getnrows(sp) + 1)
+    rowptr[1] = 1
+    for (row, colidxs) in enumerate(Ferrite.eachrow(sp))
+        for col in colidxs
+            sym && row > col && continue
+            rowptr[row + 1] += 1
+        end
+    end
+    cumsum!(rowptr, rowptr)
+    nnz = rowptr[end] - 1
+    # 2. Allocate colval and nzval now that nnz is known
+    colval = Vector{Ti}(undef, nnz)
+    nzval = zeros(Tv, nnz)
+    # 3. Populate colval.
+    k = 1
+    for (row, colidxs) in zip(1:Ferrite.getnrows(sp), Ferrite.eachrow(sp)) # pairs(eachrow(sp))
+        for col in colidxs
+            sym && row > col && continue
+            colval[k] = col
+            k += 1
+        end
+    end
+    S = SparseMatrixCSR{1}(Ferrite.getnrows(sp), Ferrite.getncols(sp), rowptr, colval, nzval)
+    return S
+end
+
+end