GitHub

Assembly

An assembler handles the insertion of the element matrices and element vectors into the system matrix and vector, and should normally (the exact interface is yet to be fully established) subtype AbstractAssembler{T}. Here T is the eltype of the contained system matrix and vector. This allows the user to infer the eltype when preallocating the element matrix and vector, e.g.

function doassemble!(assembler::Ferrite.AbstractAssembler{T}, ...) where {T}    Ke = zeros(T, n, n) # n = dofs per cell    fe = zeros(T, n)    for cell in CellIterator(...)        element_routine!(Ke, fe, cell, ...)        assemble!(assembler, celldofs(cell), Ke, fe)    endend

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 and Ferrite.CSRAssembler 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

Ferrite.allocate_matrix — Function
allocate_matrix(::Type{SparseMatrixCSC{Tv, Ti}}, sp::SparsityPattern)

Allocate a sparse matrix of type SparseMatrixCSC{Tv, Ti} from the sparsity pattern sp.

source
allocate_matrix(::Type{Symmetric{Tv, SparseMatrixCSC{Tv, Ti}}}, sp::SparsityPattern)

Instantiate a sparse matrix of type Symmetric{Tv, SparseMatrixCSC{Tv, Ti}}, i.e. a LinearAlgebra.Symmetric-wrapped SparseMatrixCSC, from the sparsity pattern sp. The resulting matrix will only store entries above, and including, the diagonal.

source
allocate_matrix(sp::SparsityPattern)

Allocate a sparse matrix of type SparseMatrixCSC{Float64, Int} from the sparsity pattern sp.

This method is a shorthand for the equivalent [allocate_matrix(SparseMatrixCSC{Float64, Int}, sp)] (@ref allocate_matrix(::Type{S}, sp::Ferrite.AbstractSparsityPattern) where {Tv, Ti, S <: SparseMatrixCSC{Tv, Ti}}).

source
allocate_matrix(MatrixType, dh::DofHandler, args...; kwargs...)

Allocate a matrix of type MatrixType from the DofHandler dh.

This is a convenience method and is equivalent to:

julia sp = init_sparsity_pattern(dh) add_sparsity_entries!(sp, dh, args...; kwargs...) allocate_matrix(MatrixType, sp)`

Refer to allocate_matrix for supported matrix types, and to init_sparsity_pattern for details about supported arguments args and keyword arguments kwargs.

Note

If more than one sparse matrix is needed (e.g. a stiffness and a mass matrix) it is more efficient to explicitly create the sparsity pattern instead of using this method, i.e. use

sp = init_sparsity_pattern(dh)add_sparsity_entries!(sp, dh)K = allocate_matrix(sp)M = allocate_matrix(sp)

instead of

K = allocate_matrix(dh)M = allocate_matrix(dh)

Note that for some matrix types it is possible to copy the instantiated matrix (M = copy(K)) instead.

source
allocate_matrix(dh::DofHandler, args...; kwargs...)

Allocate a matrix of type SparseMatrixCSC{Float64, Int} from the DofHandler dh.

This method is a shorthand for the equivalent allocate_matrix(SparseMatrixCSC{Float64, Int}, dh, args...; kwargs...) – refer to that method for details.

source
allocate_matrix(::Type{BlockMatrix}, sp::BlockSparsityPattern)
allocate_matrix(::Type{BlockMatrix{T, Matrix{S}}}, sp::BlockSparsityPattern)

Instantiate a blocked sparse matrix from the blocked sparsity pattern sp.

The type of the returned matrix is a BlockMatrix with blocks of type S (defaults to SparseMatrixCSC{T, Int}).

Examples

# Create a sparse matrix with default block typeallocate_matrix(BlockMatrix, sparsity_pattern)# Create a sparse matrix with blocks of type SparseMatrixCSC{Float32, Int}allocate_matrix(BlockMatrix{Float32, Matrix{SparseMatrixCSC{Float32, Int}}}, sparsity_pattern)
Package extension

This functionality is only enabled when the package BlockArrays.jl is installed (pkg> add BlockArrays) and loaded (using BlockArrays) in the session.

source

the internal interface

Ferrite.zero_out_columns! — Function
zero_out_columns!(K::AbstractMatrix, ch::ConstraintHandler)
zero_out_columns!(K::AbstractMatrix, columns::AbstractVector{<:Integer}, mask::AbstractVector{Bool})

Set the values of all columns associated with constrained dofs to zero.

The three argument form is the one a matrix format has to dispatch; the two argument form forwards to it. columns is the sorted list of the columns to zero and mask flags the same set (mask[j] is true iff j ∈ columns). Both forms of the same information are passed because which one can be used efficiently depends on the storage: a column-compressed format walks columns, a row-compressed one scans its stored column indices against mask. The indices are local to K, so the same method serves a whole matrix and a block of a blocked matrix.

source
Ferrite.add_inhomogeneities! — Function
add_inhomogeneities!(f::AbstractVector, K::AbstractMatrix, ch::ConstraintHandler)

Compute "f -= K*inhomogeneities".

Forwards to the four argument form below, which is the one a matrix format dispatches. The only type needing its own method here is Symmetric, which has to account for the triangle that is not stored.

source
add_inhomogeneities!(f::AbstractVector, K::AbstractMatrix, columns::AbstractVector{<:Integer}, inhomogeneities::AbstractVector)

Compute "f -= K * g", where g is zero except at columns, where it takes the values inhomogeneities. The indices are local to K, so the same method serves a whole matrix and a block of a blocked matrix (f is then a view of the corresponding rows).

This is the form a matrix format dispatches. The default below is a generic SpMSpV kernel; a format that can walk its stored entries directly should replace it, as the sparse formats Ferrite ships with do.

source
Ferrite.condense_into! — Function
condense_into!(Kdst::AbstractMatrix, K::AbstractMatrix, rowoffset::Int, coloffset::Int, dofcoefficients::Vector, dofmapping::Dict)

Condense the stored entries of K by adding their affine contributions to Kdst. rowoffset/coloffset translate the indices of K into the dof numbering of Kdst; both are zero when Kdst === K and equal the block offsets when K is a block of a blocked Kdst. This is the matrix half of Ferrite._condense!; the right hand side is condensed separately and independently of the matrix, see Ferrite._condense_rhs_column!.

This is the form a matrix format dispatches in order to support affine constraints, both on its own and as a block of a blocked matrix. Note that the contributions are written through Ferrite.addindex! on Kdst, which for a blocked Kdst is not the same object as K.

source
Ferrite.addindex! — Function
addindex!(A::AbstractMatrix{T}, v::T, i::Integer, j::Integer, ::Val{atomic} = Val(false))
addindex!(b::AbstractVector{T}, v::T, i::Integer, ::Val{atomic} = Val(false))

Equivalent to A[i, j] += v but more efficient. The optional atomic input controls whether the operation should be performed atomically (i.e. concurrency-safe) or not.

A[i, j] += v is lowered to A[i, j] = A[i, j] + v which requires a double lookup of the memory location for index (i, j) – one time for the read, and one time for the write. This method avoids the double lookup.

Zeros are ignored (i.e. if iszero(v)) by returning early. If the index (i, j) is not existing in the sparsity pattern of A this method throws a SparsityError.

Fallback: A[i, j] += v.

source
Ferrite._assemble_inner! — Function
Ferrite._assemble_inner!(K, Ke, rowdofs, sortedrowdofs, rowpermutation, coldofs, sortedcoldofs, colpermutation, sym, atomic, rowoffset, coloffset)

Scatter the element matrix Ke into the (already allocated) entries of K, i.e. K[rowdofs, coldofs] += Ke. The dofs are passed both in element order (rowdofs, coldofs) and sorted ascending (sortedrowdofs, sortedcoldofs), the latter together with the permutations mapping a sorted position back to its index in Ke, so that a format storing its entries in sorted order can walk them and the element matrix in a single pass. If sym is true only the upper triangle of Ke is read. atomic is a Val{Bool} selecting whether the accumulation is concurrency safe.

rowoffset and coloffset place the matrix within a larger system; they are only used to report global indices when an entry is missing from the sparsity pattern, and are nonzero when K is used as a block of a blocked matrix.

The default implementation writes the entries one by one with Ferrite.addindex!, which is all a custom format has to provide. A format that stores its entries sorted should specialize this and walk them together with the element matrix, as CSC and CSR do.

source

and the AbstractMatrix interface for their custom matrix type. apply! itself is generic and dispatches on AbstractMatrix, so a custom format is supported as soon as the functions above are dispatched – there is no need to add an apply! method.

Each of these takes its data – an element matrix or the constraint data – explicitly, plus index offsets where relevant, rather than an assembler or a ConstraintHandler. That is deliberate: the very same methods are then used both for a matrix on its own and for a matrix used as a block of a blocked matrix, where the indices are block local and the offsets place the block in the global system. A format that implements them therefore works with the BlockArrays extension without any further work, and without that extension having to know anything about it.

Three conventions are worth calling out:

  • _assemble_inner! is the only one of these with a working default implementation: it writes the element matrix entry by entry with addindex!, so a format is assembled into as soon as it implements that. Specializing it pays off for a format that stores its entries sorted, since the element matrix and the stored entries can then be walked in a single pass instead of searching for each entry – which is what makes it worth roughly a factor of three for CSC and CSR.
  • zero_out_rows! and zero_out_columns! receive the set of indices to zero twice, once as a sorted list and once as a boolean mask. Which one can be used efficiently depends on the storage: a column-compressed format walks the listed columns directly, while a row-compressed format has to scan its stored column indices against the mask. Passing both avoids forcing every format to build the representation it does not have.
  • condense_into! writes into a destination matrix that is not necessarily the matrix it reads, which is what lets a block condense into the blocked matrix it belongs to. It only has to handle the matrix; the right-hand side is condensed separately by Ferrite._condense!, so the order in which stored entries are visited does not matter.

Finally, Ferrite._condense! itself is dispatched per format, but is a one-liner over Ferrite.condense_into! for anything that implements it:

Ferrite._condense! — Function
_condense!(K::AbstractMatrix, f::AbstractVector, dofcoefficients::Vector{Union{Nothing, DofCoefficients{T}}}, dofmapping::Dict{<:Integer, <:Integer}, sym::Bool = false)

Condenses affine constraints K := C'KC and f := C'*f in-place, assuming the sparsity pattern is correct.

source

CSC and CSR are mirror images of one another – one stores columns contiguously, the other rows – so their implementations of the above are shared, parameterised by which index is the contiguous one. Ferrite.minor_indices is the accessor that abstracts the difference:

Ferrite.minor_indices — Function
Ferrite.minor_indices(K::AbstractSparseMatrix)

For a sparse matrix that stores the entries of one index contiguously ("the major"), return the vector of the other index ("the minor") of every stored entry, parallel to nonzeros(K). This is rowvals(K) for column-compressed storage (CSC) and colvals(K) for row-compressed storage (CSR).

This is an internal helper shared by the CSC and CSR implementations of the constraint application interface (see the devdocs on assembly); it is not part of that interface. A format that does not store scalar entries in flat arrays parallel to nonzeros – a blocked format such as BSR, for instance – simply does not define it, and implements the interface functions directly.

source

This is an implementation detail of those two formats, not part of the interface. A format that does not store scalar entries in flat arrays parallel to nonzeros – a blocked format such as BSR, say – simply implements the interface functions directly and never defines it.

Custom assembler

In case the default assembler is insufficient, users can implement a custom assembler. For this, they can create a custom type and dispatch the following functions.

Ferrite.start_assemble — Function
start_assemble(K::AbstractSparseMatrixCSC{Tv}; fillzero = true, atomic = false) -> CSCAssembler{Tv}
start_assemble(K::AbstractSparseMatrixCSC{Tv}, f::Vector{Tv}; fillzero = true, atomic = false) -> CSCAssembler{Tv}

Create a CSCAssembler{Tv} from the matrix K and optional vector f with value type Tv.

start_assemble(K::Symmetric{AbstractSparseMatrixCSC{Tv}}; fillzero = true, atomic = false) -> SymmetricCSCAssembler{Tv}start_assemble(K::Symmetric{AbstractSparseMatrixCSC{Tv}}, f::Vector = Tv[]; fillzero = true, atomic = false) -> SymmetricCSCAssembler{Tv}

Create a SymmetricCSCAssembler{Tv} from the matrix K and optional vector f with value type Tv.

CSCAssembler and SymmetricCSCAssembler allocate workspace necessary for efficient matrix assembly. To assemble the contribution from an element, use assemble!.

The keyword argument fillzero can be set to false if K and f should not be zeroed out, but instead keep their current values.

The keyword argument atomic can be set to true to make the accumulation into K and f use atomic additions. This makes it safe to assemble from multiple concurrent tasks without partitioning the cells into independent sets ("grid coloring"), at the cost of some overhead and a non-deterministic result: the order in which contributions are added to a given entry depends on the task scheduling, and floating point addition is not associative. Atomic accumulation is only supported for the value types Float16, Float32, and Float64, and Complex of these (other value types throw an ArgumentError). Note that each task still needs its own assembler since the assembler contains buffers that are modified during assemble!. Note also that the value of atomic determines a type parameter of the returned assembler, so for a type stable setup the value should be a literal (or otherwise a compile time constant). See the howto on multithreaded assembly for more details.

Depending on the loaded extensions more assembly formats become available through this interface.

source
Ferrite.assemble! — Function
assemble!(a::COOAssembler, dofs, Ke)
assemble!(a::COOAssembler, dofs, Ke, fe)

Assembles the element matrix Ke and element vector fe into a.

source
assemble!(a::COOAssembler, rowdofs, coldofs, Ke)

Assembles the matrix Ke into a according to the dofs specified by rowdofs and coldofs.

source
assemble!(g, dofs, ge, atomic = Val(false))

Assembles the element residual ge into the global residual vector g.

source
assemble!(A::Ferrite.AbstractAssembler, dofs::AbstractVector{Int}, Ke::AbstractMatrix)
assemble!(A::Ferrite.AbstractAssembler, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector)

Assemble the square element stiffness matrix Ke (and optional force vector fe) into the global stiffness (and force) in A, given the element degrees of freedom dofs.

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.

assemble!(A::Ferrite.AbstractAssembler, rowdofs::AbstractVector{Int}, coldofs::AbstractVector{Int}, Ke::AbstractMatrix)assemble!(A::Ferrite.AbstractAssembler, rowdofs::AbstractVector{Int}, coldofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector)

Assemble the element stiffness matrix Ke (and optional force vector fe) into the global stiffness (and force) in A, given the element row degrees of freedom, rowdofs, and element column degrees of freedom, coldofs. This is equivalent to K[rowdofs, coldofs] += Ke and f[rowdofs] += fe, but more efficient.

source

For local elimination support the following functions might also need custom dispatches

Ferrite._condense_local! — Function
_condense_local!(
    local_matrix::AbstractMatrix, local_vector::AbstractVector,
    global_matrix #=::SparseMatrixCSC=#, global_vector #=::Vector=#,
    global_dofs::AbstractVector, dofmapping::Dict, dofcoefficients::Vector,
    isconstrained::BitVector, atomic::Val = Val(false)
)

Condensation of affine constraints on element level. If possible this function only modifies the local arrays. Constraints reaching outside of global_dofs are written directly into global_matrix/global_vector; atomic controls whether those writes use atomic additions, which is required when assembling concurrently without grid coloring.

source

Note that apply_assemble! passes the assembler's atomic flag on to Ferrite._condense_local!, so a custom assembler supporting atomic accumulation should report it through Ferrite._is_atomic in order to make the global writes of non-local constraints concurrency-safe as well.

Type definitions

Ferrite.COOAssembler — Type
struct COOAssembler{Tv, Ti}

This assembler creates a COO (coordinate format) representation of a sparse matrix during assembly and converts it into a SparseMatrixCSC{Tv, Ti} on finalization.

source

Utility functions

Ferrite.matrix_handle — Function
matrix_handle(a::AbstractAssembler)
vector_handle(a::AbstractAssembler)

Return a reference to the underlying matrix/vector of the assembler used during assembly operations.

source
Ferrite.vector_handle — Function
matrix_handle(a::AbstractAssembler)
vector_handle(a::AbstractAssembler)

Return a reference to the underlying matrix/vector of the assembler used during assembly operations.

source
Ferrite._sortdofs_for_assembly! — Function
_sortdofs_for_assembly!(permutation::Vector{Int}, sorteddofs::Vector{Int}, dofs::AbstractVector)

Sorts the dofs into a separate buffer and returns it together with a permutation vector.

source
Ferrite.sortperm2! — Function
sortperm2!(data::AbstractVector, permutation::AbstractVector)

Sort the input vector inplace and compute the corresponding permutation.

source