Skip to content

Commit 8de6ad2

Browse files
KnutAMKnutAM
authored
Add storage eltype to abstract assembler types (#1290)
* Add data storage type to assembler types * Describe `AbstractAssembler` in devdocs * Add test for correct number type * Add tests for numeric types to CSR assembler * Add type-check for block assembler
1 parent f42a81b commit 8de6ad2

5 files changed

Lines changed: 103 additions & 85 deletions

File tree

‎docs/src/devdocs/assembly.md‎

Lines changed: 14 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,19 @@
11
# [Assembly](@id devdocs-assembly)
22

3-
The assembler handles the insertion of the element matrices and element vectors into the system matrix and right hand side.
3+
An assembler handles the insertion of the element matrices and element vectors into the system matrix and vector.
4+
and should *normally* (the exact interface is yet to be fully established) subtype `AbstractAssembler{T}`. Here `T` is the
5+
`eltype` of the contained system matrix and vector. This allows the user to infer the eltype when preallocating the element
6+
matrix and vector, e.g.
7+
```julia
8+
function doassemble!(assembler::Ferrite.AbstractAssembler{T}, ...) where {T}
9+
Ke = zeros(T, n, n) # n = dofs per cell
10+
fe = zeros(T, n)
11+
for cell in CellIterator(...)
12+
element_routine!(Ke, fe, cell, ...)
13+
assemble!(assembler, celldofs(cell), Ke, fe)
14+
end
15+
end
16+
```
417

518
## Custom matrix formats
619
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

‎ext/FerriteBlockArrays.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -61,7 +61,7 @@ end
6161
## BlockAssembler and associated methods ##
6262
###########################################
6363

64-
struct BlockAssembler{BM, Bv} <: Ferrite.AbstractAssembler
64+
struct BlockAssembler{Tv, BM <: BlockMatrix{Tv}, Bv <: AbstractVector{Tv}} <: Ferrite.AbstractAssembler{Tv}
6565
K::BM
6666
f::Bv
6767
blockindices::Vector{BlockIndex{1}}

‎src/assembler.jl‎

Lines changed: 17 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,14 @@
1-
abstract type AbstractAssembler end
2-
abstract type AbstractCSCAssembler <: AbstractAssembler end
3-
abstract type AbstractCSRAssembler <: AbstractAssembler end
1+
abstract type AbstractAssembler{Tv} end
2+
abstract type AbstractCSCAssembler{Tv} <: AbstractAssembler{Tv} end
3+
abstract type AbstractCSRAssembler{Tv} <: AbstractAssembler{Tv} end
44

55
"""
66
struct COOAssembler{Tv, Ti}
77
88
This assembler creates a COO (**coo**rdinate format) representation of a sparse matrix
99
during assembly and converts it into a `SparseMatrixCSC{Tv, Ti}` on finalization.
1010
"""
11-
struct COOAssembler{Tv, Ti} # <: AbstractAssembler
11+
struct COOAssembler{Tv, Ti} # <: AbstractAssembler{Tv}
1212
nrows::Int
1313
ncols::Int
1414
f::Vector{Tv}
@@ -173,7 +173,7 @@ matrix_handle, vector_handle
173173
"""
174174
Assembler for sparse matrix with CSC storage type.
175175
"""
176-
struct CSCAssembler{Tv, Ti, MT <: AbstractSparseMatrixCSC{Tv, Ti}} <: AbstractCSCAssembler
176+
struct CSCAssembler{Tv, Ti, MT <: AbstractSparseMatrixCSC{Tv, Ti}} <: AbstractCSCAssembler{Tv}
177177
K::MT
178178
f::Vector{Tv}
179179
rowpermutation::Vector{Int}
@@ -185,7 +185,7 @@ end
185185
"""
186186
Assembler for sparse matrix with CSR storage type.
187187
"""
188-
struct CSRAssembler{Tv, Ti, MT <: AbstractSparseMatrix{Tv, Ti}} <: AbstractCSRAssembler #AbstractSparseMatrixCSR does not exist
188+
struct CSRAssembler{Tv, Ti, MT <: AbstractSparseMatrix{Tv, Ti}} <: AbstractCSRAssembler{Tv} #AbstractSparseMatrixCSR does not exist
189189
K::MT
190190
f::Vector{Tv}
191191
rowpermutation::Vector{Int}
@@ -197,7 +197,7 @@ end
197197
"""
198198
Assembler for symmetric sparse matrix with CSC storage type.
199199
"""
200-
struct SymmetricCSCAssembler{Tv, Ti, MT <: Symmetric{Tv, <:AbstractSparseMatrixCSC{Tv, Ti}}} <: AbstractCSCAssembler
200+
struct SymmetricCSCAssembler{Tv, Ti, MT <: Symmetric{Tv, <:AbstractSparseMatrixCSC{Tv, Ti}}} <: AbstractCSCAssembler{Tv}
201201
K::MT
202202
f::Vector{Tv}
203203
rowpermutation::Vector{Int} # Symmetric assembly doesn't need separate row and
@@ -222,15 +222,15 @@ matrix_handle(a::SymmetricCSCAssembler) = a.K.data
222222
vector_handle(a::Union{AbstractCSCAssembler, AbstractCSRAssembler}) = a.f
223223

224224
"""
225-
start_assemble(K::AbstractSparseMatrixCSC; fillzero::Bool=true) -> CSCAssembler
226-
start_assemble(K::AbstractSparseMatrixCSC, f::Vector; fillzero::Bool=true) -> CSCAssembler
225+
start_assemble(K::AbstractSparseMatrixCSC{Tv}; fillzero::Bool=true) -> CSCAssembler{Tv}
226+
start_assemble(K::AbstractSparseMatrixCSC{Tv}, f::Vector{Tv}; fillzero::Bool=true) -> CSCAssembler{Tv}
227227
228-
Create a `CSCAssembler` from the matrix `K` and optional vector `f`.
228+
Create a `CSCAssembler{Tv}` from the matrix `K` and optional vector `f` with value type `Tv`.
229229
230-
start_assemble(K::Symmetric{AbstractSparseMatrixCSC}; fillzero::Bool=true) -> SymmetricCSCAssembler
231-
start_assemble(K::Symmetric{AbstractSparseMatrixCSC}, f::Vector=Td[]; fillzero::Bool=true) -> SymmetricCSCAssembler
230+
start_assemble(K::Symmetric{AbstractSparseMatrixCSC{Tv}}; fillzero::Bool=true) -> SymmetricCSCAssembler{Tv}
231+
start_assemble(K::Symmetric{AbstractSparseMatrixCSC{Tv}}, f::Vector=Tv[]; fillzero::Bool=true) -> SymmetricCSCAssembler{Tv}
232232
233-
Create a `SymmetricCSCAssembler` from the matrix `K` and optional vector `f`.
233+
Create a `SymmetricCSCAssembler{Tv}` from the matrix `K` and optional vector `f` with value type `Tv`.
234234
235235
`CSCAssembler` and `SymmetricCSCAssembler` allocate workspace
236236
necessary for efficient matrix assembly. To assemble the contribution from an element, use
@@ -259,16 +259,16 @@ function finish_assemble(a::Union{CSCAssembler, CSRAssembler, SymmetricCSCAssemb
259259
end
260260

261261
"""
262-
assemble!(A::AbstractAssembler, dofs::AbstractVector{Int}, Ke::AbstractMatrix)
263-
assemble!(A::AbstractAssembler, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector)
262+
assemble!(A::Ferrite.AbstractAssembler, dofs::AbstractVector{Int}, Ke::AbstractMatrix)
263+
assemble!(A::Ferrite.AbstractAssembler, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector)
264264
265265
Assemble the square element stiffness matrix `Ke` (and optional force vector `fe`) into the global
266266
stiffness (and force) in `A`, given the element degrees of freedom `dofs`.
267267
268268
This is equivalent to `K[dofs, dofs] += Ke` and `f[dofs] += fe`, where `K` is the global stiffness matrix and `f` the global force/residual vector, but more efficient.
269269
270-
assemble!(A::AbstractAssembler, rowdofs::AbstractVector{Int}, coldofs::AbstractVector{Int}, Ke::AbstractMatrix)
271-
assemble!(A::AbstractAssembler, rowdofs::AbstractVector{Int}, coldofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector)
270+
assemble!(A::Ferrite.AbstractAssembler, rowdofs::AbstractVector{Int}, coldofs::AbstractVector{Int}, Ke::AbstractMatrix)
271+
assemble!(A::Ferrite.AbstractAssembler, rowdofs::AbstractVector{Int}, coldofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector)
272272
273273
Assemble the element stiffness matrix `Ke` (and optional force vector `fe`) into the global
274274
stiffness (and force) in `A`, given the element row degrees of freedom, `rowdofs`, and element column degrees of freedom, `coldofs`.

‎test/blockarrays.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,7 @@ using Ferrite, BlockArrays, SparseArrays, Test
4343
@test iszero(K)
4444
@test iszero(f)
4545
block_assembler = start_assemble(KB, fB)
46+
@test isa(block_assembler, Ferrite.AbstractAssembler{Float64})
4647
@test iszero(KB)
4748
@test iszero(fB)
4849

‎test/test_assemble.jl‎

Lines changed: 70 additions & 66 deletions
Original file line numberDiff line numberDiff line change
@@ -58,75 +58,79 @@ import LinearAlgebra: Symmetric
5858
# CSCAssembler: assemble with different row and col dofs
5959
I = [1, 1, 4, 4, 6, 6]
6060
J = [1, 3, 1, 3, 1, 3]
61-
V = zeros(length(I))
62-
K = sparse(I, J, V)
63-
f = zeros(6)
64-
assembler = start_assemble(K, f)
65-
rdofs = [1, 4, 6]
66-
cdofs = [1, 3]
67-
Ke = rand(length(rdofs), length(cdofs))
68-
fe = rand(length(rdofs))
69-
assemble!(assembler, rdofs, cdofs, Ke, fe)
70-
assemble!(assembler, rdofs, cdofs, Ke, fe)
71-
@test_throws ArgumentError assemble!(assembler, rdofs, Ke, fe) # Not in sparsity pattern
72-
@test all(K[rdofs, cdofs] .== 2Ke)
73-
@test all(f[rdofs] .== 2fe)
74-
75-
# CSCAssembler: Assemble rectangular part in quadratic matrix
76-
K = SparseMatrixCSC(6, 6, [K.colptr..., 7, 7, 7], K.rowval, K.nzval)
77-
assembler = start_assemble(K, f)
78-
rdofs = [1, 4, 6]
79-
cdofs = [1, 3]
80-
Ke = rand(length(rdofs), length(cdofs))
81-
fe = rand(length(rdofs))
82-
assemble!(assembler, rdofs, cdofs, Ke, fe)
83-
assemble!(assembler, rdofs, cdofs, Ke, fe)
84-
@test_throws ArgumentError assemble!(assembler, rdofs, Ke, fe) # Not in sparsity pattern
85-
@test all(K[rdofs, cdofs] .== 2Ke)
86-
@test all(f[rdofs] .== 2fe)
87-
88-
# SparseMatrix assembler
89-
K = spzeros(10, 10)
90-
f = zeros(10)
91-
ke = [rand(4, 4), rand(4, 4)]
92-
fe = [rand(4), rand(4)]
93-
dofs = [[1, 5, 3, 7], [10, 8, 2, 5]]
94-
for i in 1:2
95-
K[dofs[i], dofs[i]] += ke[i]
96-
f[dofs[i]] += fe[i]
97-
end
98-
99-
Kc = copy(K)
100-
fc = copy(f)
101-
102-
assembler = start_assemble(Kc)
103-
@test all(iszero, Kc.nzval) # start_assemble zeroes
104-
for i in 1:2
105-
assemble!(assembler, dofs[i], ke[i])
106-
end
107-
@test Kc ≈ K
108-
109-
assembler = start_assemble(Kc, fc)
110-
@test all(iszero, Kc.nzval)
111-
@test all(iszero, fc)
112-
for i in 1:2
113-
assemble!(assembler, dofs[i], ke[i], fe[i])
114-
end
115-
@test Kc ≈ K
116-
@test fc ≈ f
117-
118-
# No zero filling
119-
assembler = start_assemble(Kc, fc; fillzero = false)
120-
@test Kc ≈ K
121-
@test fc ≈ f
122-
for i in 1:2
123-
assemble!(assembler, dofs[i], ke[i], fe[i])
61+
for T in (Float32, Float64)
62+
V = zeros(T, length(I))
63+
K = sparse(I, J, V)
64+
f = zeros(T, 6)
65+
assembler = start_assemble(K, f)
66+
@test isa(assembler, Ferrite.AbstractAssembler{T})
67+
rdofs = [1, 4, 6]
68+
cdofs = [1, 3]
69+
Ke = rand(T, length(rdofs), length(cdofs))
70+
fe = rand(T, length(rdofs))
71+
assemble!(assembler, rdofs, cdofs, Ke, fe)
72+
assemble!(assembler, rdofs, cdofs, Ke, fe)
73+
@test_throws ArgumentError assemble!(assembler, rdofs, Ke, fe) # Not in sparsity pattern
74+
@test all(K[rdofs, cdofs] .== 2Ke)
75+
@test all(f[rdofs] .== 2fe)
76+
77+
# CSCAssembler: Assemble rectangular part in quadratic matrix
78+
K = SparseMatrixCSC(6, 6, [K.colptr..., 7, 7, 7], K.rowval, K.nzval)
79+
assembler = start_assemble(K, f)
80+
rdofs = [1, 4, 6]
81+
cdofs = [1, 3]
82+
Ke = rand(T, length(rdofs), length(cdofs))
83+
fe = rand(T, length(rdofs))
84+
assemble!(assembler, rdofs, cdofs, Ke, fe)
85+
assemble!(assembler, rdofs, cdofs, Ke, fe)
86+
@test_throws ArgumentError assemble!(assembler, rdofs, Ke, fe) # Not in sparsity pattern
87+
@test all(K[rdofs, cdofs] .== 2Ke)
88+
@test all(f[rdofs] .== 2fe)
89+
90+
# SparseMatrix assembler
91+
K = spzeros(T, 10, 10)
92+
f = zeros(T, 10)
93+
ke = [rand(T, 4, 4), rand(T, 4, 4)]
94+
fe = [rand(T, 4), rand(T, 4)]
95+
dofs = [[1, 5, 3, 7], [10, 8, 2, 5]]
96+
for i in 1:2
97+
K[dofs[i], dofs[i]] += ke[i]
98+
f[dofs[i]] += fe[i]
99+
end
100+
101+
Kc = copy(K)
102+
fc = copy(f)
103+
104+
assembler = start_assemble(Kc)
105+
@test all(iszero, Kc.nzval) # start_assemble zeroes
106+
for i in 1:2
107+
assemble!(assembler, dofs[i], ke[i])
108+
end
109+
@test Kc ≈ K
110+
111+
assembler = start_assemble(Kc, fc)
112+
@test all(iszero, Kc.nzval)
113+
@test all(iszero, fc)
114+
for i in 1:2
115+
assemble!(assembler, dofs[i], ke[i], fe[i])
116+
end
117+
@test Kc ≈ K
118+
@test fc ≈ f
119+
120+
# No zero filling
121+
assembler = start_assemble(Kc, fc; fillzero = false)
122+
@test Kc ≈ K
123+
@test fc ≈ f
124+
for i in 1:2
125+
assemble!(assembler, dofs[i], ke[i], fe[i])
126+
end
127+
@test Kc ≈ 2K
128+
@test fc ≈ 2f
124129
end
125-
@test Kc ≈ 2K
126-
@test fc ≈ 2f
127130

128131
# Error paths
129-
assembler = start_assemble(Kc, fc)
132+
K = sparse(I, J, zeros(length(I)))
133+
assembler = start_assemble(K, zeros(size(K, 1)))
130134
@test_throws BoundsError assemble!(assembler, [11, 1, 2, 3], rand(4, 4))
131135
@test_throws BoundsError assemble!(assembler, [11, 1, 2, 3], rand(4, 4), rand(4))
132136
@test_throws BoundsError assemble!(assembler, [11, 1, 2], rand(4, 4))

0 commit comments

Comments
 (0)