From 7edd47a96886decc49fe5e1c69205ccbcec22e5f Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Sat, 23 Dec 2023 21:17:48 +0100 Subject: [PATCH 01/19] Some refactoring and trying to get the extension working. --- Project.toml | 3 ++ ext/FerriteSparseMatrixCSR.jl | 47 +++++++++++++++++++++++++++++++ src/assembler.jl | 52 ++++++++++++++--------------------- 3 files changed, 71 insertions(+), 31 deletions(-) create mode 100644 ext/FerriteSparseMatrixCSR.jl diff --git a/Project.toml b/Project.toml index c18ee1634a..34f3a5107b 100644 --- a/Project.toml +++ b/Project.toml @@ -16,15 +16,18 @@ 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" EnumX = "1" Metis = "1.3" +SparseMatricesCSR = "0.6" NearestNeighbors = "0.4" Preferences = "1" Reexport = "1" diff --git a/ext/FerriteSparseMatrixCSR.jl b/ext/FerriteSparseMatrixCSR.jl new file mode 100644 index 0000000000..7c8b28c033 --- /dev/null +++ b/ext/FerriteSparseMatrixCSR.jl @@ -0,0 +1,47 @@ +module FerriteSparseMatrixCSR + +using Ferrite, SparseMatricesCSR +import Base: @propagate_inbounds + +@propagate_inbounds function Ferrite._assemble_inner!(K::SparseMatrixCSR, Ke::AbstractMatrix, dofs::AbstractVector, sorteddofs::AbstractVector, permutation::AbstractVector, sym::Bool) + current_row = 1 + @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[ri] + 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) || _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 ri:maxlookups + if !iszero(Ke[Kerow, permutation[i]]) + _missing_sparsity_pattern_error(Krow, sorteddofs[i]) + end + end + current_row += 1 + end +end + +end diff --git a/src/assembler.jl b/src/assembler.jl index 573f4437dc..2aa908c828 100644 --- a/src/assembler.jl +++ b/src/assembler.jl @@ -89,7 +89,7 @@ Assembles the element residual `ge` into the global residual vector `g`. end end -abstract type AbstractSparseAssembler end +abstract type AbstractSparseAssembler{MT <: AbstractSparseMatrix} end """ matrix_handle(a::AbstractSparseAssembler) @@ -99,20 +99,14 @@ Return a reference to the underlying matrix/vector of the assembler. """ matrix_handle, vector_handle -struct AssemblerSparsityPattern{Tv,Ti} <: AbstractSparseAssembler - K::SparseMatrixCSC{Tv,Ti} - f::Vector{Tv} - permutation::Vector{Int} - sorteddofs::Vector{Int} -end -struct AssemblerSymmetricSparsityPattern{Tv,Ti} <: AbstractSparseAssembler - K::Symmetric{Tv,SparseMatrixCSC{Tv,Ti}} +struct AssemblerSparsityPattern{Tv, MT <: AbstractSparseMatrix{Tv}} <: AbstractSparseAssembler{MT} + K::MT f::Vector{Tv} permutation::Vector{Int} sorteddofs::Vector{Int} end -function Base.show(io::IO, ::MIME"text/plain", a::Union{AssemblerSparsityPattern,AssemblerSymmetricSparsityPattern}) +function Base.show(io::IO, ::MIME"text/plain", a::Union{AssemblerSparsityPattern}) print(io, typeof(a), " for assembling into:\n - ") summary(io, a.K) if !isempty(a.f) @@ -121,38 +115,28 @@ function Base.show(io::IO, ::MIME"text/plain", a::Union{AssemblerSparsityPattern end end -matrix_handle(a::AssemblerSparsityPattern) = a.K -matrix_handle(a::AssemblerSymmetricSparsityPattern) = a.K.data -vector_handle(a::Union{AssemblerSparsityPattern, AssemblerSymmetricSparsityPattern}) = a.f +matrix_handle(a::AssemblerSparsityPattern{<:Any, <:SparseMatrixCSC}) = a.K +matrix_handle(a::AssemblerSparsityPattern{<:Any, <:Symmetric}) = a.K.data +vector_handle(a::AssemblerSparsityPattern) = a.f """ - start_assemble(K::SparseMatrixCSC; fillzero::Bool=true) -> AssemblerSparsityPattern - start_assemble(K::SparseMatrixCSC, f::Vector; fillzero::Bool=true) -> AssemblerSparsityPattern + start_assemble(K::AbstractSparseMatrix; fillzero::Bool=true) -> AssemblerSparsityPattern + start_assemble(K::AbstractSparseMatrix, f::Vector; fillzero::Bool=true) -> AssemblerSparsityPattern Create a `AssemblerSparsityPattern` from the matrix `K` and optional vector `f`. - start_assemble(K::Symmetric{SparseMatrixCSC}; fillzero::Bool=true) -> AssemblerSymmetricSparsityPattern - start_assemble(K::Symmetric{SparseMatrixCSC}, f::Vector=Td[]; fillzero::Bool=true) -> AssemblerSymmetricSparsityPattern - -Create a `AssemblerSymmetricSparsityPattern` from the matrix `K` and optional vector `f`. - -`AssemblerSparsityPattern` and `AssemblerSymmetricSparsityPattern` allocate workspace -necessary for efficient matrix assembly. To assemble the contribution from an element, use -[`assemble!`](@ref). +`AssemblerSparsityPattern` allocate workspace necessary for efficient matrix assembly. To assemble + the contribution from an element, use [`assemble!`](@ref). The keyword argument `fillzero` can be set to `false` if `K` and `f` should not be zeroed out, but instead keep their current values. """ -start_assemble(K::Union{SparseMatrixCSC, Symmetric{<:Any,SparseMatrixCSC}}, f::Vector; fillzero::Bool) +start_assemble(K::AbstractSparseMatrix, f::Vector; fillzero::Bool) -function start_assemble(K::SparseMatrixCSC{T}, f::Vector=T[]; fillzero::Bool=true) where {T} +function start_assemble(K::Union{AbstractSparseMatrix{T}, Symmetric{T,<: AbstractSparseMatrix{T}}}, f::Vector=T[]; fillzero::Bool=true) where {T} fillzero && (fillzero!(K); fillzero!(f)) return AssemblerSparsityPattern(K, f, Int[], Int[]) end -function start_assemble(K::Symmetric{T,<:SparseMatrixCSC}, f::Vector=T[]; fillzero::Bool=true) where T - fillzero && (fillzero!(K); fillzero!(f)) - return AssemblerSymmetricSparsityPattern(K, f, Int[], Int[]) -end """ assemble!(A::AbstractSparseAssembler, dofs::AbstractVector{Int}, Ke::AbstractMatrix) @@ -172,13 +156,14 @@ end @propagate_inbounds function assemble!(A::AbstractSparseAssembler, dofs::AbstractVector{Int}, fe::AbstractVector, Ke::AbstractMatrix) assemble!(A, dofs, Ke, fe) end -@propagate_inbounds function assemble!(A::AssemblerSparsityPattern, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector) +@propagate_inbounds function assemble!(A::AssemblerSparsityPattern{<:Any, <:AbstractSparseMatrix}, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector) _assemble!(A, dofs, Ke, fe, false) end -@propagate_inbounds function assemble!(A::AssemblerSymmetricSparsityPattern, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector) +@propagate_inbounds function assemble!(A::AssemblerSparsityPattern{<:Any, <:Symmetric{<:Any,<:AbstractSparseMatrix}}, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector) _assemble!(A, dofs, Ke, fe, true) end +# Main entry point for the CPU @propagate_inbounds function _assemble!(A::AbstractSparseAssembler, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector, sym::Bool) ld = length(dofs) @boundscheck checkbounds(Ke, keys(dofs), keys(dofs)) @@ -197,7 +182,12 @@ end copyto!(sorteddofs, dofs) sortperm2!(sorteddofs, permutation) + _assemble_inner!(K, Ke, dofs, sorteddofs, permutation, sym) +end + +@propagate_inbounds function _assemble_inner!(K::SparseMatrixCSC, Ke::AbstractMatrix, dofs::AbstractVector, sorteddofs::AbstractVector, permutation::AbstractVector, sym::Bool) current_col = 1 + ld = length(dofs) @inbounds for Kcol in sorteddofs maxlookups = sym ? current_col : ld Kecol = permutation[current_col] From 7c8e7b1cd94e99fd3d15f3db3bab8ddc84342f1a Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Sat, 23 Dec 2023 22:47:06 +0100 Subject: [PATCH 02/19] Almost. Imhomogeneities are broken. --- docs/src/literate-tutorials/heat_equation.jl | 11 +-- ext/FerriteSparseMatrixCSR.jl | 57 +++++++++++- src/Dofs/ConstraintHandler.jl | 94 +++++++++++--------- src/arrayutils.jl | 4 +- src/assembler.jl | 4 +- 5 files changed, 115 insertions(+), 55 deletions(-) diff --git a/docs/src/literate-tutorials/heat_equation.jl b/docs/src/literate-tutorials/heat_equation.jl index 4a82dac7b7..9afac34209 100644 --- a/docs/src/literate-tutorials/heat_equation.jl +++ b/docs/src/literate-tutorials/heat_equation.jl @@ -41,7 +41,7 @@ #md # The full program, without comments, can be found in the next [section](@ref heat_equation-plain-program). # # First we load Ferrite, and some other packages we need -using Ferrite, SparseArrays +using Ferrite, SparseArrays, SparseMatricesCSR # We start by generating a simple grid with 20x20 quadrilateral elements # using `generate_grid`. The generator defaults to the unit square, # so we don't need to specify the corners of the domain. @@ -73,7 +73,8 @@ close!(dh); # Now that we have distributed all our dofs we can create our tangent matrix, # using `create_sparsity_pattern`. This function returns a sparse matrix # with the correct entries stored. -K = create_sparsity_pattern(dh) +K = SparseMatrixCSR(transpose(create_sparsity_pattern(dh))) +#K = create_sparsity_pattern(dh) # ### Boundary conditions # In Ferrite constraints like Dirichlet boundary conditions @@ -98,7 +99,7 @@ ch = ConstraintHandler(dh); # this case do not depend on time we define our function as $f(\textbf{x}) = 0$, i.e. # no matter what $\textbf{x}$ we return $0$. When we have # specified our constraint we `add!` it to `ch`. -dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0) +dbc = Dirichlet(:u, ∂Ω, (x, t) -> 1.0) add!(ch, dbc); # Finally we also need to `close!` our constraint handler. When we call `close!` @@ -177,7 +178,7 @@ end # versions. However, through the code we use `f` and `u` instead to reflect the strong # connection between the weak form and the Ferrite implementation. -function assemble_global(cellvalues::CellValues, K::SparseMatrixCSC, dh::DofHandler) +function assemble_global(cellvalues::CellValues, K::AbstractSparseMatrix, dh::DofHandler) ## Allocate the element stiffness matrix and element force vector n_basefuncs = getnbasefunctions(cellvalues) Ke = zeros(n_basefuncs, n_basefuncs) @@ -219,7 +220,7 @@ end ## test the result #src using Test #src -@test norm(u) ≈ 3.307743912641305 #src +@test norm(u) ≈ 23.74835617035439 #src #md # ## [Plain program](@id heat_equation-plain-program) #md # diff --git a/ext/FerriteSparseMatrixCSR.jl b/ext/FerriteSparseMatrixCSR.jl index 7c8b28c033..c396e54354 100644 --- a/ext/FerriteSparseMatrixCSR.jl +++ b/ext/FerriteSparseMatrixCSR.jl @@ -1,10 +1,12 @@ module FerriteSparseMatrixCSR -using Ferrite, SparseMatricesCSR +using Ferrite, SparseArrays, SparseMatricesCSR import Base: @propagate_inbounds +import Ferrite: Symmetric @propagate_inbounds function Ferrite._assemble_inner!(K::SparseMatrixCSR, Ke::AbstractMatrix, dofs::AbstractVector, sorteddofs::AbstractVector, permutation::AbstractVector, sym::Bool) current_row = 1 + ld = length(dofs) @inbounds for Krow in sorteddofs maxlookups = sym ? current_row : ld Kerow = permutation[current_row] @@ -14,7 +16,7 @@ import Base: @propagate_inbounds while Ci <= length(nzr) && ci <= maxlookups C = nzr[Ci] Kcol = K.colval[C] - Kecol = permutation[ri] + Kecol = permutation[ci] val = Ke[Kerow, Kecol] if Kcol == dofs[Kecol] # Match: add the value (if non-zero) and advance the pointers @@ -35,7 +37,7 @@ import Base: @propagate_inbounds end end # Make sure that remaining entries in this column of the local matrix are all zero - for i in ri:maxlookups + for i in ci:maxlookups if !iszero(Ke[Kerow, permutation[i]]) _missing_sparsity_pattern_error(Krow, sorteddofs[i]) end @@ -44,4 +46,53 @@ import Base: @propagate_inbounds end end +function Ferrite.add_inhomogeneities!(K::SparseMatrixCSR, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector, dofmapping, sym::Bool) + @inbounds for i in 1:length(inhomogeneities) + d = prescribed_dofs[i] + for j in nzrange(K, d) + c = K.colval[j] + sym && c > d && break # don't look below diagonal + if (i = get(dofmapping, c, 0); i != 0) + f[d] -= inhomogeneities[i] * K.nzval[j] + end + end + end + if sym + error("Symmetric inhomogeneities not supported for SparseMatrixCSR") + # In the symmetric case, for a constrained dof `d`, we handle the contribution + # from `K[1:d, d]` in the loop above, but we are still missing the contribution + # from `K[(d+1):size(K,1), d]`. These values are not stored, but since the + # matrix is symmetric we can instead use `K[d, (d+1):size(K,1)]`. Looping over + # rows is slow, so loop over all columns again, and check if the row is a + # constrained row. + # @inbounds for col in 1:size(K, 2) + # for ri in nzrange(K, col) + # row = K.rowval[ri] + # row >= col && break + # if (i = get(dofmapping, row, 0); i != 0) + # f[col] -= inhomogeneities[i] * K.nzval[ri] + # end + # end + # end + end +end + +function Ferrite.zero_out_rows!(K::SparseMatrixCSR, ch::ConstraintHandler) # can be removed in 0.7 with #24711 merged + @debug @assert issorted(ch.prescribed_dofs) + for row in ch.prescribed_dofs + r = nzrange(K, row) + K.nzval[r] .= 0.0 + end +end + +function Ferrite.zero_out_columns!(K::SparseMatrixCSR, ch::ConstraintHandler) + colval = K.colval + nzval = K.nzval + @inbounds for i in eachindex(colval, nzval) + if haskey(ch.dofmapping, colval[i]) + nzval[i] = 0 + end + end +end + end diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index 57c48909ac..29e742c77b 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -524,7 +524,7 @@ end """ - apply!(K::SparseMatrixCSC, rhs::AbstractVector, ch::ConstraintHandler) + apply!(K::AbstractSparseMatrix, rhs::AbstractVector, ch::ConstraintHandler) Adjust the matrix `K` and right hand side `rhs` to account for the Dirichlet boundary conditions specified in `ch` such that `K \\ rhs` gives the expected solution. @@ -616,11 +616,11 @@ function _apply_v(v::AbstractVector, ch::ConstraintHandler, apply_zero::Bool) return v end -function apply!(K::Union{SparseMatrixCSC,Symmetric}, ch::ConstraintHandler) +function apply!(K::Union{AbstractSparseMatrix,Symmetric}, ch::ConstraintHandler) apply!(K, eltype(K)[], ch, true) end -function apply_zero!(K::Union{SparseMatrixCSC,Symmetric}, f::AbstractVector, ch::ConstraintHandler) +function apply_zero!(K::Union{AbstractSparseMatrix,Symmetric}, f::AbstractVector, ch::ConstraintHandler) apply!(K, f, ch, true) end @@ -629,7 +629,7 @@ end const APPLY_TRANSPOSE = ApplyStrategy.Transpose const APPLY_INPLACE = ApplyStrategy.Inplace -function apply!(KK::Union{SparseMatrixCSC,Symmetric}, f::AbstractVector, ch::ConstraintHandler, applyzero::Bool=false; +function apply!(KK::Union{AbstractSparseMatrix,Symmetric{<:Any,AbstractSparseMatrix}}, f::AbstractVector, ch::ConstraintHandler, applyzero::Bool=false; strategy::ApplyStrategy.T=ApplyStrategy.Transpose) @assert isclosed(ch) sym = isa(KK, Symmetric) @@ -638,48 +638,18 @@ function apply!(KK::Union{SparseMatrixCSC,Symmetric}, f::AbstractVector, ch::Con @boundscheck checkbounds(K, ch.prescribed_dofs, ch.prescribed_dofs) @boundscheck length(f) == 0 || checkbounds(f, ch.prescribed_dofs) - m = meandiag(K) # Use the mean of the diagonal here to not ruin things for iterative solver - # Add inhomogeneities to f: (f - K * ch.inhomogeneities) - if !applyzero - @inbounds for i in 1:length(ch.inhomogeneities) - d = ch.prescribed_dofs[i] - v = ch.inhomogeneities[i] - if v != 0 - for j in nzrange(K, d) - r = K.rowval[j] - sym && r > d && break # don't look below diagonal - f[r] -= v * K.nzval[j] - end - end - end - if sym - # In the symmetric case, for a constrained dof `d`, we handle the contribution - # from `K[1:d, d]` in the loop above, but we are still missing the contribution - # from `K[(d+1):size(K,1), d]`. These values are not stored, but since the - # matrix is symmetric we can instead use `K[d, (d+1):size(K,1)]`. Looping over - # rows is slow, so loop over all columns again, and check if the row is a - # constrained row. - @inbounds for col in 1:size(K, 2) - for ri in nzrange(K, col) - row = K.rowval[ri] - row >= col && break - if (i = get(ch.dofmapping, row, 0); i != 0) - f[col] -= ch.inhomogeneities[i] * K.nzval[ri] - end - end - end - end - end + !applyzero && add_inhomogeneities!(K, f, ch.inhomogeneities, ch.prescribed_dofs, ch.dofmapping, sym) # Condense K (C' * K * C) and f (C' * f) _condense!(K, f, ch.dofcoefficients, ch.dofmapping, sym) # Remove constrained dofs from the matrix - zero_out_columns!(K, ch.prescribed_dofs) - zero_out_rows!(K, ch.dofmapping) + zero_out_columns!(K, ch) + zero_out_rows!(K, ch) # Add meandiag to constraint dofs + m = meandiag(K) # Use the mean of the diagonal here to not ruin things for iterative solver @inbounds for i in 1:length(ch.inhomogeneities) d = ch.prescribed_dofs[i] K[d, d] = m @@ -690,6 +660,37 @@ function apply!(KK::Union{SparseMatrixCSC,Symmetric}, f::AbstractVector, ch::Con end end +function add_inhomogeneities!(K::SparseMatrixCSC, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector, dofmapping, sym::Bool) + @inbounds for i in 1:length(inhomogeneities) + d = prescribed_dofs[i] + v = inhomogeneities[i] + if v != 0 + for j in nzrange(K, d) + r = K.rowval[j] + sym && r > d && break # don't look below diagonal + f[r] -= v * K.nzval[j] + end + end + end + if sym + # In the symmetric case, for a constrained dof `d`, we handle the contribution + # from `K[1:d, d]` in the loop above, but we are still missing the contribution + # from `K[(d+1):size(K,1), d]`. These values are not stored, but since the + # matrix is symmetric we can instead use `K[d, (d+1):size(K,1)]`. Looping over + # rows is slow, so loop over all columns again, and check if the row is a + # constrained row. + @inbounds for col in 1:size(K, 2) + for ri in nzrange(K, col) + row = K.rowval[ri] + row >= col && break + if (i = get(dofmapping, row, 0); i != 0) + f[col] -= inhomogeneities[i] * K.nzval[ri] + end + end + end + end +end + # Fetch dof coefficients for a dof prescribed by an affine constraint. Return nothing if the # dof is not prescribed, or prescribed by DBC. @inline function coefficients_for_dof(dofmapping, dofcoeffs, dof) @@ -698,6 +699,13 @@ end return dofcoeffs[idx] end + +function _condense!(K::AbstractSparseMatrix, f::AbstractVector, dofcoefficients::Vector{Union{Nothing, DofCoefficients{T}}}, dofmapping::Dict{Int,Int}, sym::Bool=false) where T + # Return early if there are no non-trivial affine constraints + any(i -> !(i === nothing || isempty(i)), dofcoefficients) || return + error("condensation of ::$(typeof(K)) matrix not supported") +end + # Condenses K and f: C'*K*C, C'*f, in-place assuming the sparsity pattern is correct function _condense!(K::SparseMatrixCSC, f::AbstractVector, dofcoefficients::Vector{Union{Nothing, DofCoefficients{T}}}, dofmapping::Dict{Int,Int}, sym::Bool=false) where T @@ -808,19 +816,19 @@ function create_constraint_matrix(ch::ConstraintHandler{dh,T}) where {dh,T} end # columns need to be stored entries, this is not checked -function zero_out_columns!(K, dofs::Vector{Int}) # can be removed in 0.7 with #24711 merged - @debug @assert issorted(dofs) - for col in dofs +function zero_out_columns!(K::SparseArrays.AbstractSparseMatrixCSC, ch::ConstraintHandler) # can be removed in 0.7 with #24711 merged + @debug @assert issorted(ch.prescribed_dofs) + for col in ch.prescribed_dofs r = nzrange(K, col) K.nzval[r] .= 0.0 end end -function zero_out_rows!(K, dofmapping) +function zero_out_rows!(K::SparseArrays.AbstractSparseMatrixCSC, ch::ConstraintHandler) rowval = K.rowval nzval = K.nzval @inbounds for i in eachindex(rowval, nzval) - if haskey(dofmapping, rowval[i]) + if haskey(ch.dofmapping, rowval[i]) nzval[i] = 0 end end diff --git a/src/arrayutils.jl b/src/arrayutils.jl index 2566ee338d..0e74e002c8 100644 --- a/src/arrayutils.jl +++ b/src/arrayutils.jl @@ -77,11 +77,11 @@ function addindex!(A::SparseMatrixCSC{Tv}, v::Tv, i::Int, j::Int) where Tv end end -function fillzero!(A::SparseMatrixCSC{T}) where T +function fillzero!(A::AbstractSparseMatrix{T}) where T fill!(nonzeros(A), zero(T)) return A end -function fillzero!(A::Symmetric{T,<:SparseMatrixCSC}) where T +function fillzero!(A::Symmetric{T,<:AbstractSparseMatrix}) where T fillzero!(A.data) return A end diff --git a/src/assembler.jl b/src/assembler.jl index 2aa908c828..99c1f65910 100644 --- a/src/assembler.jl +++ b/src/assembler.jl @@ -115,7 +115,7 @@ function Base.show(io::IO, ::MIME"text/plain", a::Union{AssemblerSparsityPattern end end -matrix_handle(a::AssemblerSparsityPattern{<:Any, <:SparseMatrixCSC}) = a.K +matrix_handle(a::AssemblerSparsityPattern{<:Any, <:AbstractSparseMatrix}) = a.K matrix_handle(a::AssemblerSparsityPattern{<:Any, <:Symmetric}) = a.K.data vector_handle(a::AssemblerSparsityPattern) = a.f @@ -165,7 +165,6 @@ end # Main entry point for the CPU @propagate_inbounds function _assemble!(A::AbstractSparseAssembler, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector, sym::Bool) - ld = length(dofs) @boundscheck checkbounds(Ke, keys(dofs), keys(dofs)) if length(fe) != 0 @boundscheck checkbounds(fe, keys(dofs)) @@ -177,6 +176,7 @@ end permutation = A.permutation sorteddofs = A.sorteddofs @boundscheck checkbounds(K, dofs, dofs) + ld = length(dofs) resize!(permutation, ld) resize!(sorteddofs, ld) copyto!(sorteddofs, dofs) From a29dbb0e9680f7a08cec80ecffcc8d12670d0090 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Sat, 23 Dec 2023 22:53:50 +0100 Subject: [PATCH 03/19] revert faulty refactoring. --- src/Dofs/ConstraintHandler.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index 29e742c77b..171fe76ab3 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -638,6 +638,8 @@ function apply!(KK::Union{AbstractSparseMatrix,Symmetric{<:Any,AbstractSparseMat @boundscheck checkbounds(K, ch.prescribed_dofs, ch.prescribed_dofs) @boundscheck length(f) == 0 || checkbounds(f, ch.prescribed_dofs) + m = meandiag(K) # Use the mean of the diagonal here to not ruin things for iterative solver + # Add inhomogeneities to f: (f - K * ch.inhomogeneities) !applyzero && add_inhomogeneities!(K, f, ch.inhomogeneities, ch.prescribed_dofs, ch.dofmapping, sym) @@ -649,7 +651,6 @@ function apply!(KK::Union{AbstractSparseMatrix,Symmetric{<:Any,AbstractSparseMat zero_out_rows!(K, ch) # Add meandiag to constraint dofs - m = meandiag(K) # Use the mean of the diagonal here to not ruin things for iterative solver @inbounds for i in 1:length(ch.inhomogeneities) d = ch.prescribed_dofs[i] K[d, d] = m From d8589eee05f9ab8e9639419cbcb7f8b6f2f441e4 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Wed, 27 Dec 2023 18:21:46 +0100 Subject: [PATCH 04/19] Refactor SpMSpV kernels and add tests --- Project.toml | 3 ++- ext/FerriteSparseMatrixCSR.jl | 32 -------------------------- src/Dofs/ConstraintHandler.jl | 16 +++++++++---- test/runtests.jl | 1 + test/test_assembler_extensions.jl | 37 +++++++++++++++++++++++++++++++ test/test_constraints.jl | 1 + 6 files changed, 53 insertions(+), 37 deletions(-) create mode 100644 test/test_assembler_extensions.jl diff --git a/Project.toml b/Project.toml index 34f3a5107b..ad71316b9d 100644 --- a/Project.toml +++ b/Project.toml @@ -49,8 +49,9 @@ NBInclude = "0db19996-df87-5ea3-a455-e3a50d440464" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" +SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [targets] -test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "JET", "Metis", "NBInclude", "ProgressMeter", "Random", "SHA", "Test", "TimerOutputs", "Logging"] +test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "JET", "Metis", "NBInclude", "ProgressMeter", "Random", "SHA", "SparseMatricesCSR", "Test", "TimerOutputs", "Logging"] diff --git a/ext/FerriteSparseMatrixCSR.jl b/ext/FerriteSparseMatrixCSR.jl index c396e54354..836945d5bb 100644 --- a/ext/FerriteSparseMatrixCSR.jl +++ b/ext/FerriteSparseMatrixCSR.jl @@ -2,7 +2,6 @@ module FerriteSparseMatrixCSR using Ferrite, SparseArrays, SparseMatricesCSR import Base: @propagate_inbounds -import Ferrite: Symmetric @propagate_inbounds function Ferrite._assemble_inner!(K::SparseMatrixCSR, Ke::AbstractMatrix, dofs::AbstractVector, sorteddofs::AbstractVector, permutation::AbstractVector, sym::Bool) current_row = 1 @@ -46,37 +45,6 @@ import Ferrite: Symmetric end end -function Ferrite.add_inhomogeneities!(K::SparseMatrixCSR, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector, dofmapping, sym::Bool) - @inbounds for i in 1:length(inhomogeneities) - d = prescribed_dofs[i] - for j in nzrange(K, d) - c = K.colval[j] - sym && c > d && break # don't look below diagonal - if (i = get(dofmapping, c, 0); i != 0) - f[d] -= inhomogeneities[i] * K.nzval[j] - end - end - end - if sym - error("Symmetric inhomogeneities not supported for SparseMatrixCSR") - # In the symmetric case, for a constrained dof `d`, we handle the contribution - # from `K[1:d, d]` in the loop above, but we are still missing the contribution - # from `K[(d+1):size(K,1), d]`. These values are not stored, but since the - # matrix is symmetric we can instead use `K[d, (d+1):size(K,1)]`. Looping over - # rows is slow, so loop over all columns again, and check if the row is a - # constrained row. - # @inbounds for col in 1:size(K, 2) - # for ri in nzrange(K, col) - # row = K.rowval[ri] - # row >= col && break - # if (i = get(dofmapping, row, 0); i != 0) - # f[col] -= inhomogeneities[i] * K.nzval[ri] - # end - # end - # end - end -end - function Ferrite.zero_out_rows!(K::SparseMatrixCSR, ch::ConstraintHandler) # can be removed in 0.7 with #24711 merged @debug @assert issorted(ch.prescribed_dofs) for row in ch.prescribed_dofs diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index 171fe76ab3..41e0eb9bba 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -629,7 +629,7 @@ end const APPLY_TRANSPOSE = ApplyStrategy.Transpose const APPLY_INPLACE = ApplyStrategy.Inplace -function apply!(KK::Union{AbstractSparseMatrix,Symmetric{<:Any,AbstractSparseMatrix}}, f::AbstractVector, ch::ConstraintHandler, applyzero::Bool=false; +function apply!(KK::Union{AbstractSparseMatrix,Symmetric{<:Any,<:AbstractSparseMatrix}}, f::AbstractVector, ch::ConstraintHandler, applyzero::Bool=false; strategy::ApplyStrategy.T=ApplyStrategy.Transpose) @assert isclosed(ch) sym = isa(KK, Symmetric) @@ -641,9 +641,9 @@ function apply!(KK::Union{AbstractSparseMatrix,Symmetric{<:Any,AbstractSparseMat m = meandiag(K) # Use the mean of the diagonal here to not ruin things for iterative solver # Add inhomogeneities to f: (f - K * ch.inhomogeneities) - !applyzero && add_inhomogeneities!(K, f, ch.inhomogeneities, ch.prescribed_dofs, ch.dofmapping, sym) + !applyzero && add_inhomogeneities!(KK, f, ch.inhomogeneities, ch.prescribed_dofs, ch.dofmapping) - # Condense K (C' * K * C) and f (C' * f) + # Condense K := (C' * K * C) and f := (C' * f) _condense!(K, f, ch.dofcoefficients, ch.dofmapping, sym) # Remove constrained dofs from the matrix @@ -661,7 +661,15 @@ function apply!(KK::Union{AbstractSparseMatrix,Symmetric{<:Any,AbstractSparseMat end end -function add_inhomogeneities!(K::SparseMatrixCSC, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector, dofmapping, sym::Bool) +# Generic version to compute "f -= K*inhomogeneities" +function add_inhomogeneities!(K, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping) + f .-= K*sparsevec(prescribed_dofs, inhomogeneities, size(K,2)) +end + +# Optimized version for SparseMatrixCSC +add_inhomogeneities!(K::SparseMatrixCSC, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping) = add_inhomogeneities_csc!(K, f, inhomogeneities, prescribed_dofs, dofmapping, false) +add_inhomogeneities!(K::Symmetric{<:Any,<:SparseMatrixCSC}, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping) = add_inhomogeneities_csc!(K, f, inhomogeneities, prescribed_dofs, dofmapping, true) +function add_inhomogeneities_csc!(K::SparseMatrixCSC, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping, sym::Bool) @inbounds for i in 1:length(inhomogeneities) d = prescribed_dofs[i] v = inhomogeneities[i] diff --git a/test/runtests.jl b/test/runtests.jl index 72f7d71e53..37f310ca6b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -50,6 +50,7 @@ include("test_apply_rhs.jl") include("test_apply_analytical.jl") include("test_deprecations.jl") HAS_EXTENSIONS && include("blockarrays.jl") +HAS_EXTENSIONS && include("text_assembler_extensions.jl") include("test_examples.jl") @test all(x -> isdefined(Ferrite, x), names(Ferrite)) # Test that all exported symbols are defined diff --git a/test/test_assembler_extensions.jl b/test/test_assembler_extensions.jl new file mode 100644 index 0000000000..2f296c4244 --- /dev/null +++ b/test/test_assembler_extensions.jl @@ -0,0 +1,37 @@ +import SparseMatricesCSR: SparseMatrixCSR +using SparseArrays, LinearAlgebra + +@testset "apply!(::SparseMatrixCSR,...)" begin + # Specifically this test that values below the diagonal of K2::Symmetric aren't touched + # and that the missing values are instead taken from above the diagonal. + grid = generate_grid(Line, (2,)) + dh = DofHandler(grid) + add!(dh, :u, Lagrange{RefLine,1}()) + close!(dh) + ch = ConstraintHandler(dh) + add!(ch, Dirichlet(:u, getfaceset(grid, "left"), x -> 1)) + close!(ch) + K0 = sparse(rand(3, 3)) + K0 = K0'*K0 + K1 = SparseMatrixCSR(transpose(copy(K0))) + K2 = Symmetric(SparseMatrixCSR(transpose(copy(K0)))) + @test K0 == K1 + @test K1 == K2 + sol = [1.0, 2.0, 3.0] + f0 = K0*sol + f1 = K1*sol + f2 = K2*sol + apply!(K0, f0, ch) + apply!(K1, f1, ch) + apply!(K2, f2, ch) + @test K0 == K1 + @test K1 == K2 + @test f0 == f1 + @test f1 == f2 + @test K1\f1 ≈ K0\f0 + # Error for affine constraints + ch = ConstraintHandler(dh) + add!(ch, AffineConstraint(1, [3 => 1.0], 1.0)) + close!(ch) + @test_throws ErrorException("condensation of ::SparseMatrixCSR{1, Float64, Int64} matrix not supported") apply!(K2, f2, ch) +end diff --git a/test/test_constraints.jl b/test/test_constraints.jl index 5b511632e5..b2cbb649ae 100644 --- a/test/test_constraints.jl +++ b/test/test_constraints.jl @@ -430,6 +430,7 @@ end # Linear solution fill!(a, NaN) # (not used) doassemble!(K, r, dh, a, params) + KT = SparseMatrixCSR(transpose(copy(K))) f = zero(r); f[end] = params.f apply!(K, f, ch) a_linear = K\f From 26c0f5830406981d29f8b0285ae67b89007fe080 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Wed, 27 Dec 2023 18:25:15 +0100 Subject: [PATCH 05/19] Revert heat example. --- docs/src/literate-tutorials/heat_equation.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/docs/src/literate-tutorials/heat_equation.jl b/docs/src/literate-tutorials/heat_equation.jl index 9afac34209..4a82dac7b7 100644 --- a/docs/src/literate-tutorials/heat_equation.jl +++ b/docs/src/literate-tutorials/heat_equation.jl @@ -41,7 +41,7 @@ #md # The full program, without comments, can be found in the next [section](@ref heat_equation-plain-program). # # First we load Ferrite, and some other packages we need -using Ferrite, SparseArrays, SparseMatricesCSR +using Ferrite, SparseArrays # We start by generating a simple grid with 20x20 quadrilateral elements # using `generate_grid`. The generator defaults to the unit square, # so we don't need to specify the corners of the domain. @@ -73,8 +73,7 @@ close!(dh); # Now that we have distributed all our dofs we can create our tangent matrix, # using `create_sparsity_pattern`. This function returns a sparse matrix # with the correct entries stored. -K = SparseMatrixCSR(transpose(create_sparsity_pattern(dh))) -#K = create_sparsity_pattern(dh) +K = create_sparsity_pattern(dh) # ### Boundary conditions # In Ferrite constraints like Dirichlet boundary conditions @@ -99,7 +98,7 @@ ch = ConstraintHandler(dh); # this case do not depend on time we define our function as $f(\textbf{x}) = 0$, i.e. # no matter what $\textbf{x}$ we return $0$. When we have # specified our constraint we `add!` it to `ch`. -dbc = Dirichlet(:u, ∂Ω, (x, t) -> 1.0) +dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0) add!(ch, dbc); # Finally we also need to `close!` our constraint handler. When we call `close!` @@ -178,7 +177,7 @@ end # versions. However, through the code we use `f` and `u` instead to reflect the strong # connection between the weak form and the Ferrite implementation. -function assemble_global(cellvalues::CellValues, K::AbstractSparseMatrix, dh::DofHandler) +function assemble_global(cellvalues::CellValues, K::SparseMatrixCSC, dh::DofHandler) ## Allocate the element stiffness matrix and element force vector n_basefuncs = getnbasefunctions(cellvalues) Ke = zeros(n_basefuncs, n_basefuncs) @@ -220,7 +219,7 @@ end ## test the result #src using Test #src -@test norm(u) ≈ 23.74835617035439 #src +@test norm(u) ≈ 3.307743912641305 #src #md # ## [Plain program](@id heat_equation-plain-program) #md # From 056f5c39debc6f51f60e2a39e4dbda3471436fe8 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Wed, 27 Dec 2023 18:47:24 +0100 Subject: [PATCH 06/19] Devdocs. --- docs/src/devdocs/assembler.md | 33 +++++++++++++++++++++++++++++++++ docs/src/devdocs/index.md | 2 +- 2 files changed, 34 insertions(+), 1 deletion(-) create mode 100644 docs/src/devdocs/assembler.md diff --git a/docs/src/devdocs/assembler.md b/docs/src/devdocs/assembler.md new file mode 100644 index 0000000000..43c5a98e71 --- /dev/null +++ b/docs/src/devdocs/assembler.md @@ -0,0 +1,33 @@ +# [Assembler](@id devdocs-assembler) + +The assembler handles the insertion of the element matrices and element vectors into the system matrix and right hand side. 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 assembler [`AssemblerSparsityPattern`](@ref) should be able to handle most cases in practice. To support a custom format users have to dispatch the following functions: + +```@docs +Ferrite._assemble_inner! +Ferrite.zero_out_rows! +Ferrite.zero_out_columns! +``` + +and the `AbstractSparseMatrix` interface for their custom matrix type. + +## 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 +Ferrite.matrix_handle +Ferrite.vector_handle +start_assemble! +finish_assemble! +assemble! +``` + +For local elimination support the following functions might also need custom dispatches + +```@docs +Ferrite._condense_local! +``` + +## Custom Matrix Type Sparsity Pattern + +TODO `create_sparsity_pattern` diff --git a/docs/src/devdocs/index.md b/docs/src/devdocs/index.md index 9c16b83d7c..f292cdb55d 100644 --- a/docs/src/devdocs/index.md +++ b/docs/src/devdocs/index.md @@ -5,5 +5,5 @@ developing the library. ```@contents Depth = 1 -Pages = ["reference_cells.md", "interpolations.md", "elements.md", "FEValues.md", "dofhandler.md", "performance.md"] +Pages = ["reference_cells.md", "interpolations.md", "elements.md", "FEValues.md", "dofhandler.md", "assembler.md", "performance.md"] ``` From f021d016d1e55e734e8b29a8d27f91e7ff206e56 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Wed, 27 Dec 2023 18:49:42 +0100 Subject: [PATCH 07/19] Symmetric assembler. --- src/assembler.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/assembler.jl b/src/assembler.jl index 99c1f65910..c8f2be7bac 100644 --- a/src/assembler.jl +++ b/src/assembler.jl @@ -99,7 +99,7 @@ Return a reference to the underlying matrix/vector of the assembler. """ matrix_handle, vector_handle -struct AssemblerSparsityPattern{Tv, MT <: AbstractSparseMatrix{Tv}} <: AbstractSparseAssembler{MT} +struct AssemblerSparsityPattern{Tv, MT <: Union{AbstractSparseMatrix{Tv}, Symmetric{Tv,<:AbstractSparseMatrix{Tv}}} <: AbstractSparseAssembler{MT} K::MT f::Vector{Tv} permutation::Vector{Int} From 46f6428d0c861657c9d7f32da3419d70ec578392 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Thu, 28 Dec 2023 00:51:00 +0100 Subject: [PATCH 08/19] Derp --- src/assembler.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/assembler.jl b/src/assembler.jl index c8f2be7bac..dbeb1918e8 100644 --- a/src/assembler.jl +++ b/src/assembler.jl @@ -99,7 +99,7 @@ Return a reference to the underlying matrix/vector of the assembler. """ matrix_handle, vector_handle -struct AssemblerSparsityPattern{Tv, MT <: Union{AbstractSparseMatrix{Tv}, Symmetric{Tv,<:AbstractSparseMatrix{Tv}}} <: AbstractSparseAssembler{MT} +struct AssemblerSparsityPattern{Tv, MT <: Union{AbstractSparseMatrix{Tv}, Symmetric{Tv,<:AbstractSparseMatrix{Tv}}}} <: AbstractSparseAssembler{MT} K::MT f::Vector{Tv} permutation::Vector{Int} From 4a47a3f0a6064ef1a7f1149c361835ecb0cadc20 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Thu, 28 Dec 2023 01:29:16 +0100 Subject: [PATCH 09/19] Revert constraints test battery. --- test/test_constraints.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_constraints.jl b/test/test_constraints.jl index b2cbb649ae..5b511632e5 100644 --- a/test/test_constraints.jl +++ b/test/test_constraints.jl @@ -430,7 +430,6 @@ end # Linear solution fill!(a, NaN) # (not used) doassemble!(K, r, dh, a, params) - KT = SparseMatrixCSR(transpose(copy(K))) f = zero(r); f[end] = params.f apply!(K, f, ch) a_linear = K\f From d227ef7b75768d95ac1a07b13d30f8f67f94c413 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Thu, 28 Dec 2023 01:30:09 +0100 Subject: [PATCH 10/19] Add sparsity pattern generator for CSR. --- ext/FerriteSparseMatrixCSR.jl | 7 +++++++ src/assembler.jl | 4 ++-- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/ext/FerriteSparseMatrixCSR.jl b/ext/FerriteSparseMatrixCSR.jl index 836945d5bb..ceecdc8ab6 100644 --- a/ext/FerriteSparseMatrixCSR.jl +++ b/ext/FerriteSparseMatrixCSR.jl @@ -63,4 +63,11 @@ function Ferrite.zero_out_columns!(K::SparseMatrixCSR, ch::ConstraintHandler) end end +function Ferrite.create_sparsity_pattern(::Type{<:SparseMatrixCSR}, dh, ch; kwargs...) + # create SparseMatrixCSC + K = create_sparsity_pattern(dh, ch; kwargs...) + # transform to SparseMatrixCSR + return SparseMatrixCSR(transpose(sparse(transpose(K)))) +end + end diff --git a/src/assembler.jl b/src/assembler.jl index dbeb1918e8..c6af883706 100644 --- a/src/assembler.jl +++ b/src/assembler.jl @@ -89,7 +89,7 @@ Assembles the element residual `ge` into the global residual vector `g`. end end -abstract type AbstractSparseAssembler{MT <: AbstractSparseMatrix} end +abstract type AbstractSparseAssembler end """ matrix_handle(a::AbstractSparseAssembler) @@ -99,7 +99,7 @@ Return a reference to the underlying matrix/vector of the assembler. """ matrix_handle, vector_handle -struct AssemblerSparsityPattern{Tv, MT <: Union{AbstractSparseMatrix{Tv}, Symmetric{Tv,<:AbstractSparseMatrix{Tv}}}} <: AbstractSparseAssembler{MT} +struct AssemblerSparsityPattern{Tv, MT <: Union{AbstractSparseMatrix{Tv}, Symmetric{Tv,<:AbstractSparseMatrix{Tv}}}} <: AbstractSparseAssembler K::MT f::Vector{Tv} permutation::Vector{Int} From 0c3efc81e297baec7161410e7fe89292ec0b687b Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Thu, 28 Dec 2023 01:52:47 +0100 Subject: [PATCH 11/19] Derp --- src/Dofs/ConstraintHandler.jl | 2 +- test/runtests.jl | 2 +- test/test_assembler_extensions.jl | 1 - 3 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index 41e0eb9bba..f5d2d8182d 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -668,7 +668,7 @@ end # Optimized version for SparseMatrixCSC add_inhomogeneities!(K::SparseMatrixCSC, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping) = add_inhomogeneities_csc!(K, f, inhomogeneities, prescribed_dofs, dofmapping, false) -add_inhomogeneities!(K::Symmetric{<:Any,<:SparseMatrixCSC}, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping) = add_inhomogeneities_csc!(K, f, inhomogeneities, prescribed_dofs, dofmapping, true) +add_inhomogeneities!(K::Symmetric{<:Any,<:SparseMatrixCSC}, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping) = add_inhomogeneities_csc!(K.data, f, inhomogeneities, prescribed_dofs, dofmapping, true) function add_inhomogeneities_csc!(K::SparseMatrixCSC, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping, sym::Bool) @inbounds for i in 1:length(inhomogeneities) d = prescribed_dofs[i] diff --git a/test/runtests.jl b/test/runtests.jl index 37f310ca6b..5b5aa8bd06 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -50,7 +50,7 @@ include("test_apply_rhs.jl") include("test_apply_analytical.jl") include("test_deprecations.jl") HAS_EXTENSIONS && include("blockarrays.jl") -HAS_EXTENSIONS && include("text_assembler_extensions.jl") +HAS_EXTENSIONS && include("test_assembler_extensions.jl") include("test_examples.jl") @test all(x -> isdefined(Ferrite, x), names(Ferrite)) # Test that all exported symbols are defined diff --git a/test/test_assembler_extensions.jl b/test/test_assembler_extensions.jl index 2f296c4244..21fb446396 100644 --- a/test/test_assembler_extensions.jl +++ b/test/test_assembler_extensions.jl @@ -28,7 +28,6 @@ using SparseArrays, LinearAlgebra @test K1 == K2 @test f0 == f1 @test f1 == f2 - @test K1\f1 ≈ K0\f0 # Error for affine constraints ch = ConstraintHandler(dh) add!(ch, AffineConstraint(1, [3 => 1.0], 1.0)) From a2160e0945796cf09205e93c9390c87f3dd2fc8c Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Thu, 28 Dec 2023 02:47:04 +0100 Subject: [PATCH 12/19] Convenience function for sparsity pattern. --- ext/FerriteSparseMatrixCSR.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/ext/FerriteSparseMatrixCSR.jl b/ext/FerriteSparseMatrixCSR.jl index ceecdc8ab6..f2bfbcfe20 100644 --- a/ext/FerriteSparseMatrixCSR.jl +++ b/ext/FerriteSparseMatrixCSR.jl @@ -70,4 +70,11 @@ function Ferrite.create_sparsity_pattern(::Type{<:SparseMatrixCSR}, dh, ch; kwar return SparseMatrixCSR(transpose(sparse(transpose(K)))) end +function Ferrite.create_sparsity_pattern(::Type{<:SparseMatrixCSR}, dh; kwargs...) + # create SparseMatrixCSC + K = create_sparsity_pattern(dh; kwargs...) + # transform to SparseMatrixCSR + return SparseMatrixCSR(transpose(sparse(transpose(K)))) +end + end From 8102774acabd2bbdf1b9c1f05e62659bd4556d1f Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Thu, 28 Dec 2023 04:10:59 +0100 Subject: [PATCH 13/19] Test coverage for sparsity pattern, assembly and Dirichlet elimination. --- test/test_assembler_extensions.jl | 45 ++++++++++++++++++++++++++++++- 1 file changed, 44 insertions(+), 1 deletion(-) diff --git a/test/test_assembler_extensions.jl b/test/test_assembler_extensions.jl index 21fb446396..328a9ff3e5 100644 --- a/test/test_assembler_extensions.jl +++ b/test/test_assembler_extensions.jl @@ -1,6 +1,8 @@ -import SparseMatricesCSR: SparseMatrixCSR +import SparseMatricesCSR: SparseMatrixCSR, sparsecsr using SparseArrays, LinearAlgebra +@testset "SparseMatricesCSR extension" begin + @testset "apply!(::SparseMatrixCSR,...)" begin # Specifically this test that values below the diagonal of K2::Symmetric aren't touched # and that the missing values are instead taken from above the diagonal. @@ -34,3 +36,44 @@ using SparseArrays, LinearAlgebra close!(ch) @test_throws ErrorException("condensation of ::SparseMatrixCSR{1, Float64, Int64} matrix not supported") apply!(K2, f2, ch) end + +@testset "assembly integration" begin + grid = generate_grid(Line, (2,)) + dh = DofHandler(grid) + add!(dh, :u, Lagrange{RefLine, 1}()) + close!(dh) + + K = create_sparsity_pattern(SparseMatrixCSR, dh) + I = [1,1,2,2,2,3,3] + J = [1,2,1,2,3,2,3] + V = zeros(7) + @test K == sparsecsr(I,J,V) + f = zeros(3) + + ch = ConstraintHandler(dh) + add!(ch, Dirichlet(:u, getfaceset(grid, "left"), (x, t) -> 1)) + close!(ch) + @test K == create_sparsity_pattern(SparseMatrixCSR, dh, ch) + + assembler = start_assemble(K, f) + ke = [-1.0 1.0; 2.0 -1.0] + fe = [1.0,2.0] + assemble!(assembler, [1,2], ke,fe) + assemble!(assembler, [3,2], ke,fe) + + I = [1,1,2,2,2,3,3] + J = [1,2,1,2,3,2,3] + V = [-1.0,1.0,2.0,-2.0,2.0,1.0,-1.0] + @test K ≈ sparsecsr(I,J,V) + @test f ≈ [1.0,4.0,1.0] + + apply!(K,f,ch) + + I = [1,1,2,2,2,3,3] + J = [1,2,1,2,3,2,3] + V = [4/3,0.0,0.0,-2.0,2.0,1.0,-1.0] + @test K ≈ sparsecsr(I,J,V) + @test f ≈ [4/3,2.0,1.0] +end + +end From cfabf75e2e427db0a3f8d5f93f2c0604ad4fae0c Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Thu, 28 Dec 2023 13:57:41 +0100 Subject: [PATCH 14/19] Devdocs for add_inhomogeneities --- docs/src/devdocs/assembler.md | 6 +++++- src/Dofs/ConstraintHandler.jl | 7 ++++++- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/docs/src/devdocs/assembler.md b/docs/src/devdocs/assembler.md index 43c5a98e71..28b9df8ade 100644 --- a/docs/src/devdocs/assembler.md +++ b/docs/src/devdocs/assembler.md @@ -8,7 +8,11 @@ Ferrite.zero_out_rows! Ferrite.zero_out_columns! ``` -and the `AbstractSparseMatrix` interface for their custom matrix type. +and the `AbstractSparseMatrix` interface for their custom matrix type. Optional dispatches to speed up operations might be + +```@docs +Ferrite.add_inhomogeneities! +``` ## Custom Assembler diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index f5d2d8182d..57ef59261d 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -661,7 +661,12 @@ function apply!(KK::Union{AbstractSparseMatrix,Symmetric{<:Any,<:AbstractSparseM end end -# Generic version to compute "f -= K*inhomogeneities" +""" + add_inhomogeneities!(K, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping::Dict{Int,Int}) + +Compute "f -= K*inhomogeneities". +By default this is a generic version via SpMSpV kernel. +""" function add_inhomogeneities!(K, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping) f .-= K*sparsevec(prescribed_dofs, inhomogeneities, size(K,2)) end From 0328f2487340199800512f5fa77952ac84da4f09 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Thu, 28 Dec 2023 15:12:13 +0100 Subject: [PATCH 15/19] Add Knut's suggestion. --- ext/FerriteSparseMatrixCSR.jl | 4 ++-- test/test_assembler_extensions.jl | 25 +++++++++++++++++++++++++ 2 files changed, 27 insertions(+), 2 deletions(-) diff --git a/ext/FerriteSparseMatrixCSR.jl b/ext/FerriteSparseMatrixCSR.jl index f2bfbcfe20..df99c55074 100644 --- a/ext/FerriteSparseMatrixCSR.jl +++ b/ext/FerriteSparseMatrixCSR.jl @@ -30,7 +30,7 @@ import Base: @propagate_inbounds 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) || _missing_sparsity_pattern_error(Krow, Kcol) + iszero(val) || Ferrite._missing_sparsity_pattern_error(Krow, Kcol) # Advance the local matrix row pointer ci += 1 end @@ -38,7 +38,7 @@ import Base: @propagate_inbounds # 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]]) - _missing_sparsity_pattern_error(Krow, sorteddofs[i]) + Ferrite._missing_sparsity_pattern_error(Krow, sorteddofs[i]) end end current_row += 1 diff --git a/test/test_assembler_extensions.jl b/test/test_assembler_extensions.jl index 328a9ff3e5..80372ffeef 100644 --- a/test/test_assembler_extensions.jl +++ b/test/test_assembler_extensions.jl @@ -74,6 +74,31 @@ end V = [4/3,0.0,0.0,-2.0,2.0,1.0,-1.0] @test K ≈ sparsecsr(I,J,V) @test f ≈ [4/3,2.0,1.0] + + + grid = generate_grid(Quadrilateral, (2,2)) + ip = Lagrange{RefQuadrilateral,1}() + dh = DofHandler(grid) + add!(dh, :u, ip) + add!(dh, :v, ip) + close!(dh) + + Ke_zeros = zeros(ndofs_per_cell(dh,1), ndofs_per_cell(dh,1)) + Ke_rand = rand(ndofs_per_cell(dh,1), ndofs_per_cell(dh,1)) + dofs = celldofs(dh,1) + + for c1 ∈ [true, false], c2 ∈ [true, false], c3 ∈ [true, false], c4 ∈ [true, false] + coupling = [c1; c2;; c3; c4] + K = create_sparsity_pattern(SparseMatrixCSR, dh; coupling) + a = start_assemble(K) + assemble!(a, dofs, Ke_zeros) + if all(coupling) + assemble!(a, dofs, Ke_rand) + @test Ke_rand ≈ K[dofs,dofs] + else + @test_throws ErrorException assemble!(a, dofs, Ke_rand) + end + end end end From ed365f1deda852563b476d0513800ae90000a638 Mon Sep 17 00:00:00 2001 From: termi-official Date: Tue, 2 Jul 2024 16:19:49 +0200 Subject: [PATCH 16/19] Trim ws --- ext/FerriteSparseMatrixCSR.jl | 4 ++-- src/assembler.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/ext/FerriteSparseMatrixCSR.jl b/ext/FerriteSparseMatrixCSR.jl index df99c55074..07ad01b4b6 100644 --- a/ext/FerriteSparseMatrixCSR.jl +++ b/ext/FerriteSparseMatrixCSR.jl @@ -67,14 +67,14 @@ function Ferrite.create_sparsity_pattern(::Type{<:SparseMatrixCSR}, dh, ch; kwar # create SparseMatrixCSC K = create_sparsity_pattern(dh, ch; kwargs...) # transform to SparseMatrixCSR - return SparseMatrixCSR(transpose(sparse(transpose(K)))) + return SparseMatrixCSR(transpose(sparse(transpose(K)))) end function Ferrite.create_sparsity_pattern(::Type{<:SparseMatrixCSR}, dh; kwargs...) # create SparseMatrixCSC K = create_sparsity_pattern(dh; kwargs...) # transform to SparseMatrixCSR - return SparseMatrixCSR(transpose(sparse(transpose(K)))) + return SparseMatrixCSR(transpose(sparse(transpose(K)))) end end diff --git a/src/assembler.jl b/src/assembler.jl index a7d5473380..ef3f5a3f8a 100644 --- a/src/assembler.jl +++ b/src/assembler.jl @@ -125,7 +125,7 @@ vector_handle(a::AssemblerSparsityPattern) = a.f Create a `AssemblerSparsityPattern` from the matrix `K` and optional vector `f`. -`AssemblerSparsityPattern` allocate workspace necessary for efficient matrix assembly. To assemble +`AssemblerSparsityPattern` allocate workspace necessary for efficient matrix assembly. To assemble the contribution from an element, use [`assemble!`](@ref). The keyword argument `fillzero` can be set to `false` if `K` and `f` should not be zeroed From 29888f6c2d649747ac5ea1dfa461f88c07f9a8cd Mon Sep 17 00:00:00 2001 From: termi-official Date: Tue, 2 Jul 2024 16:20:53 +0200 Subject: [PATCH 17/19] Trim ws again --- test/test_assembler_extensions.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_assembler_extensions.jl b/test/test_assembler_extensions.jl index 80372ffeef..e8b183ead6 100644 --- a/test/test_assembler_extensions.jl +++ b/test/test_assembler_extensions.jl @@ -82,11 +82,11 @@ end add!(dh, :u, ip) add!(dh, :v, ip) close!(dh) - + Ke_zeros = zeros(ndofs_per_cell(dh,1), ndofs_per_cell(dh,1)) Ke_rand = rand(ndofs_per_cell(dh,1), ndofs_per_cell(dh,1)) dofs = celldofs(dh,1) - + for c1 ∈ [true, false], c2 ∈ [true, false], c3 ∈ [true, false], c4 ∈ [true, false] coupling = [c1; c2;; c3; c4] K = create_sparsity_pattern(SparseMatrixCSR, dh; coupling) From 60bf45fd004db393be69e555a1203a1b7b1263fc Mon Sep 17 00:00:00 2001 From: termi-official Date: Tue, 2 Jul 2024 16:47:47 +0200 Subject: [PATCH 18/19] Update to new API --- ext/FerriteSparseMatrixCSR.jl | 19 +++++++++---------- src/Ferrite.jl | 2 +- test/test_assembler_extensions.jl | 10 +++++----- 3 files changed, 15 insertions(+), 16 deletions(-) diff --git a/ext/FerriteSparseMatrixCSR.jl b/ext/FerriteSparseMatrixCSR.jl index 07ad01b4b6..07d24605ae 100644 --- a/ext/FerriteSparseMatrixCSR.jl +++ b/ext/FerriteSparseMatrixCSR.jl @@ -1,6 +1,7 @@ module FerriteSparseMatrixCSR using Ferrite, SparseArrays, SparseMatricesCSR +import Ferrite: AbstractSparsityPattern import Base: @propagate_inbounds @propagate_inbounds function Ferrite._assemble_inner!(K::SparseMatrixCSR, Ke::AbstractMatrix, dofs::AbstractVector, sorteddofs::AbstractVector, permutation::AbstractVector, sym::Bool) @@ -63,18 +64,16 @@ function Ferrite.zero_out_columns!(K::SparseMatrixCSR, ch::ConstraintHandler) end end -function Ferrite.create_sparsity_pattern(::Type{<:SparseMatrixCSR}, dh, ch; kwargs...) - # create SparseMatrixCSC - K = create_sparsity_pattern(dh, ch; kwargs...) - # transform to SparseMatrixCSR - return SparseMatrixCSR(transpose(sparse(transpose(K)))) + +function Ferrite.allocate_matrix(::Type{SparseMatrixCSR}, sp::AbstractSparsityPattern) + return allocate_matrix(SparseMatrixCSR{1,Float64,Int}, sp) end -function Ferrite.create_sparsity_pattern(::Type{<:SparseMatrixCSR}, dh; kwargs...) - # create SparseMatrixCSC - K = create_sparsity_pattern(dh; kwargs...) - # transform to SparseMatrixCSR - return SparseMatrixCSR(transpose(sparse(transpose(K)))) +function Ferrite.allocate_matrix(::Type{MatrixType}, sp::AbstractSparsityPattern) where {Bi, Tv, Ti, MatrixType <: SparseMatrixCSR{Bi, Tv, Ti}} + # Allocate CSC first ... + K = allocate_matrix(SparseMatrixCSC{Tv, Ti}, sp) + # ... and transform to SparseMatrixCSR + return SparseMatrixCSR{Bi}(transpose(sparse(transpose(K)))) end end diff --git a/src/Ferrite.jl b/src/Ferrite.jl index ca5ef421e5..86fb1c6c94 100644 --- a/src/Ferrite.jl +++ b/src/Ferrite.jl @@ -14,7 +14,7 @@ using NearestNeighbors: using OrderedCollections: OrderedSet using SparseArrays: - SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals, sparse + AbstractSparseMatrix, SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals, sparse, sparsevec using StaticArrays: StaticArrays, MArray, MMatrix, SArray, SMatrix, SVector using WriteVTK: diff --git a/test/test_assembler_extensions.jl b/test/test_assembler_extensions.jl index e8b183ead6..6ab20e00e0 100644 --- a/test/test_assembler_extensions.jl +++ b/test/test_assembler_extensions.jl @@ -11,7 +11,7 @@ using SparseArrays, LinearAlgebra add!(dh, :u, Lagrange{RefLine,1}()) close!(dh) ch = ConstraintHandler(dh) - add!(ch, Dirichlet(:u, getfaceset(grid, "left"), x -> 1)) + add!(ch, Dirichlet(:u, getfacetset(grid, "left"), x -> 1)) close!(ch) K0 = sparse(rand(3, 3)) K0 = K0'*K0 @@ -43,7 +43,7 @@ end add!(dh, :u, Lagrange{RefLine, 1}()) close!(dh) - K = create_sparsity_pattern(SparseMatrixCSR, dh) + K = allocate_matrix(SparseMatrixCSR, dh) I = [1,1,2,2,2,3,3] J = [1,2,1,2,3,2,3] V = zeros(7) @@ -51,9 +51,9 @@ end f = zeros(3) ch = ConstraintHandler(dh) - add!(ch, Dirichlet(:u, getfaceset(grid, "left"), (x, t) -> 1)) + add!(ch, Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> 1)) close!(ch) - @test K == create_sparsity_pattern(SparseMatrixCSR, dh, ch) + @test K == allocate_matrix(SparseMatrixCSR, dh, ch) assembler = start_assemble(K, f) ke = [-1.0 1.0; 2.0 -1.0] @@ -89,7 +89,7 @@ end for c1 ∈ [true, false], c2 ∈ [true, false], c3 ∈ [true, false], c4 ∈ [true, false] coupling = [c1; c2;; c3; c4] - K = create_sparsity_pattern(SparseMatrixCSR, dh; coupling) + K = allocate_matrix(SparseMatrixCSR, dh; coupling) a = start_assemble(K) assemble!(a, dofs, Ke_zeros) if all(coupling) From 0b5a861c9a41bfcc232a9c56f9c4a249516084af Mon Sep 17 00:00:00 2001 From: termi-official Date: Tue, 2 Jul 2024 17:02:24 +0200 Subject: [PATCH 19/19] [skip ci] Comments. --- test/test_assembler_extensions.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/test/test_assembler_extensions.jl b/test/test_assembler_extensions.jl index 6ab20e00e0..feb2e2568c 100644 --- a/test/test_assembler_extensions.jl +++ b/test/test_assembler_extensions.jl @@ -38,11 +38,13 @@ using SparseArrays, LinearAlgebra end @testset "assembly integration" begin + # Setup simple problem grid = generate_grid(Line, (2,)) dh = DofHandler(grid) add!(dh, :u, Lagrange{RefLine, 1}()) close!(dh) + # Check if the matrix is correctly allocated K = allocate_matrix(SparseMatrixCSR, dh) I = [1,1,2,2,2,3,3] J = [1,2,1,2,3,2,3] @@ -50,32 +52,33 @@ end @test K == sparsecsr(I,J,V) f = zeros(3) + # Check that incuding the ch doesnot mess up the pattern ch = ConstraintHandler(dh) add!(ch, Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> 1)) close!(ch) @test K == allocate_matrix(SparseMatrixCSR, dh, ch) + # Check if assembly works assembler = start_assemble(K, f) ke = [-1.0 1.0; 2.0 -1.0] fe = [1.0,2.0] assemble!(assembler, [1,2], ke,fe) assemble!(assembler, [3,2], ke,fe) - I = [1,1,2,2,2,3,3] J = [1,2,1,2,3,2,3] V = [-1.0,1.0,2.0,-2.0,2.0,1.0,-1.0] @test K ≈ sparsecsr(I,J,V) @test f ≈ [1.0,4.0,1.0] + # Check if constraint handler integration works apply!(K,f,ch) - I = [1,1,2,2,2,3,3] J = [1,2,1,2,3,2,3] V = [4/3,0.0,0.0,-2.0,2.0,1.0,-1.0] @test K ≈ sparsecsr(I,J,V) @test f ≈ [4/3,2.0,1.0] - + # Check if coupling works grid = generate_grid(Quadrilateral, (2,2)) ip = Lagrange{RefQuadrilateral,1}() dh = DofHandler(grid)