diff --git a/CHANGELOG.md b/CHANGELOG.md index 10113c654e..72b1d453e0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -244,6 +244,22 @@ more discussion). ``` +- **Sparsity pattern and global matrix construction**: since there is now explicit support + for working with the sparsity pattern before instantiating a matrix the function + `create_sparsity_pattern` has been removed. To recover the old functionality that return a + sparse matrix from the DofHandler directly use `allocate_matrix` instead. + + Examples: + ```diff + # Create sparse matrix from DofHandler + - K = create_sparsity_pattern(dh) + + K = allocate_matrix(dh) + + # Create condensed sparse matrix from DofHandler + ConstraintHandler + - K = create_sparsity_pattern(dh, ch) + + K = allocate_matrix(dh, ch) + ``` + ### Added - `InterfaceValues` for computing jumps and averages over interfaces. ([#743][github-743]) @@ -317,6 +333,16 @@ more discussion). a given grid (based on its node coordinates), and returns the minimum and maximum vertices of the bounding box. ([#880][github-880]) +- Support for working with sparsity patterns has been added. This means that Ferrite exposes + the intermediate "state" between the DofHandler and the instantiated matrix as the new + struct `SparsityPattern`. This make it possible to insert custom equations or couplings in + the pattern before instantiating the matrix. The function `create_sparsity_pattern` have + been removed. The new function `allocate_matrix` is instead used to instantiate the + matrix. Refer to the documentation for more details. ([#888][github-888]) + + **To upgrade**: if you want to recover the old functionality and don't need to work with + the pattern, replace any usage of `create_sparsity_pattern` with `allocate_matrix`. + - A new function, `geometric_interpolation`, is exported, which gives the geometric interpolation for each cell type. This is equivalent to the deprecated `Ferrite.default_interpolation` function. ([#953][github-953]) @@ -330,9 +356,9 @@ more discussion). ### Changed -- `create_sparsity_pattern` now supports cross-element dof coupling by passing kwarg - `topology` along with an optional `cross_coupling` matrix that behaves similar to - the `coupling` kwarg. ([#710][github-#710]) +- It is now possible to create sparsity patterns with interface couplings, see the new + function `add_interface_entries!` and the rework of sparsity pattern construction. + ([#710][github-#710]) - The `AbstractCell` interface has been reworked. This change should not affect user code, but may in some cases be relevant for code parsing external mesh files. In particular, the @@ -1003,6 +1029,7 @@ poking into Ferrite internals: [github-835]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/835 [github-855]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/855 [github-880]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/880 +[github-888]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/888 [github-914]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/914 [github-924]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/924 [github-938]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/938 diff --git a/benchmark/helper.jl b/benchmark/helper.jl index 8602e48a3f..369681b1db 100644 --- a/benchmark/helper.jl +++ b/benchmark/helper.jl @@ -1,6 +1,7 @@ module FerriteBenchmarkHelper using Ferrite +using LinearAlgebra: Symmetric function geo_types_for_spatial_dim(spatial_dim) spatial_dim == 1 && return [Line, QuadraticLine] @@ -177,7 +178,7 @@ function _assemble_mass(dh, cellvalues, sym) Me = zeros(n_basefuncs, n_basefuncs) fe = zeros(n_basefuncs) - M = sym ? create_symmetric_sparsity_pattern(dh) : create_sparsity_pattern(dh); + M = sym ? Symmetric(allocate_matrix(dh)) : allocate_matrix(dh); f = zeros(ndofs(dh)) assembler = start_assemble(M, f); diff --git a/docs/liveserver.jl b/docs/liveserver.jl index f2c099eaa2..a589e910c1 100755 --- a/docs/liveserver.jl +++ b/docs/liveserver.jl @@ -14,12 +14,14 @@ push!(ARGS, "liveserver") # Run LiveServer.servedocs(...) import LiveServer LiveServer.servedocs(; + host = "0.0.0.0", # Documentation root where make.jl and src/ are located foldername = joinpath(repo_root, "docs"), # Extra source folder to watch for changes include_dirs = [ - # Watch the src folder so docstrings can be Revise'd + # Watch the src and ext folder so docstrings can be Revise'd joinpath(repo_root, "src"), + joinpath(repo_root, "ext"), ], skip_dirs = [ # Skip the folder where Literate.jl output is written. This is needed diff --git a/docs/make.jl b/docs/make.jl index 133333cc45..f6dae950f0 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -10,7 +10,10 @@ if liveserver @timeit dto "Revise.revise()" Revise.revise() end -using Documenter, DocumenterCitations, Ferrite, FerriteGmsh, FerriteMeshParser +using Documenter, DocumenterCitations, Ferrite, FerriteGmsh, FerriteMeshParser, SparseArrays, LinearAlgebra + +using BlockArrays +const FerriteBlockArrays = Base.get_extension(Ferrite, :FerriteBlockArrays) const is_ci = haskey(ENV, "GITHUB_ACTIONS") @@ -66,6 +69,7 @@ bibtex_plugin = CitationBibliography( "topics/reference_shapes.md", "topics/FEValues.md", "topics/degrees_of_freedom.md", + "topics/sparse_matrix.md", "topics/assembly.md", "topics/boundary_conditions.md", "topics/constraints.md", @@ -78,6 +82,7 @@ bibtex_plugin = CitationBibliography( "reference/interpolations.md", "reference/fevalues.md", "reference/dofhandler.md", + "reference/sparsity_pattern.md", "reference/assembly.md", "reference/boundary_conditions.md", "reference/grid.md", diff --git a/docs/src/devdocs/dofhandler.md b/docs/src/devdocs/dofhandler.md index 96a5cd6c04..3b4feda17d 100644 --- a/docs/src/devdocs/dofhandler.md +++ b/docs/src/devdocs/dofhandler.md @@ -25,6 +25,4 @@ Ferrite.find_field(dh::DofHandler, field_name::Symbol) Ferrite._close_subdofhandler! Ferrite._distribute_dofs_for_cell! Ferrite.permute_and_push! -Ferrite.cross_element_coupling! -Ferrite._add_cross_coupling ``` diff --git a/docs/src/literate-gallery/helmholtz.jl b/docs/src/literate-gallery/helmholtz.jl index 17afe674cd..59668d1a56 100644 --- a/docs/src/literate-gallery/helmholtz.jl +++ b/docs/src/literate-gallery/helmholtz.jl @@ -85,7 +85,7 @@ add!(dbcs, dbc) close!(dbcs) update!(dbcs, 0.0) -K = create_sparsity_pattern(dh); +K = allocate_matrix(dh); function doassemble(cellvalues::CellValues, facetvalues::FacetValues, K::SparseMatrixCSC, dh::DofHandler) diff --git a/docs/src/literate-gallery/landau.jl b/docs/src/literate-gallery/landau.jl index b84283c869..654f813a36 100644 --- a/docs/src/literate-gallery/landau.jl +++ b/docs/src/literate-gallery/landau.jl @@ -105,7 +105,7 @@ function LandauModel(α, G, gridsize, left::Vec{DIM, T}, right::Vec{DIM, T}, elp apply!(dofvector, boundaryconds) - hessian = create_sparsity_pattern(dofhandler) + hessian = allocate_matrix(dofhandler) dpc = ndofs_per_cell(dofhandler) cpc = length(grid.cells[1].nodes) caches = [ThreadCache(dpc, cpc, copy(cvP), ModelParams(α, G), elpotential) for t=1:nthreads()] @@ -193,7 +193,7 @@ function minimize!(model; kwargs...) dh = model.dofhandler dofs = model.dofs ∇f = fill(0.0, length(dofs)) - ∇²f = create_sparsity_pattern(dh) + ∇²f = allocate_matrix(dh) function g!(storage, x) ∇F!(storage, x, model) apply_zero!(storage, model.boundaryconds) diff --git a/docs/src/literate-gallery/quasi_incompressible_hyperelasticity.jl b/docs/src/literate-gallery/quasi_incompressible_hyperelasticity.jl index 6f23947b6d..96287579a0 100644 --- a/docs/src/literate-gallery/quasi_incompressible_hyperelasticity.jl +++ b/docs/src/literate-gallery/quasi_incompressible_hyperelasticity.jl @@ -296,7 +296,7 @@ function solve(interpolation_u, interpolation_p) apply!(w, dbc) ## Create the sparse matrix and residual vector - K = create_sparsity_pattern(dh) + K = allocate_matrix(dh) f = zeros(_ndofs) ## We run the simulation parameterized by a time like parameter. `Tf` denotes the final value diff --git a/docs/src/literate-gallery/topology_optimization.jl b/docs/src/literate-gallery/topology_optimization.jl index a498a2ae04..823756d13e 100644 --- a/docs/src/literate-gallery/topology_optimization.jl +++ b/docs/src/literate-gallery/topology_optimization.jl @@ -428,7 +428,7 @@ function topopt(ra,ρ,n,filename; output=:false) χ = zeros(getncells(dh.grid)) r = zeros(n_dofs) # residual - K = create_sparsity_pattern(dh) # stiffness matrix + K = allocate_matrix(dh) # stiffness matrix i_max = 300 ## maximum number of iteration steps tol = 1e-4 diff --git a/docs/src/literate-howto/threaded_assembly.jl b/docs/src/literate-howto/threaded_assembly.jl index 2f21fb33b1..3b814a380e 100644 --- a/docs/src/literate-howto/threaded_assembly.jl +++ b/docs/src/literate-howto/threaded_assembly.jl @@ -182,7 +182,7 @@ function run_assemble() ip = Lagrange{RefHexahedron,1}()^3 dh = create_dofhandler(grid, ip); - K = create_sparsity_pattern(dh); + K = allocate_matrix(dh); C = create_stiffness(Val{3}()); ## compilation doassemble(K, colors, grid, dh, C, ip); diff --git a/docs/src/literate-tutorials/computational_homogenization.jl b/docs/src/literate-tutorials/computational_homogenization.jl index 95eca9c6f5..945c02fe06 100644 --- a/docs/src/literate-tutorials/computational_homogenization.jl +++ b/docs/src/literate-tutorials/computational_homogenization.jl @@ -277,8 +277,8 @@ ch = (dirichlet = ch_dirichlet, periodic = ch_periodic); # and the constraint handler. K = ( - dirichlet = create_sparsity_pattern(dh), - periodic = create_sparsity_pattern(dh, ch.periodic), + dirichlet = allocate_matrix(dh), + periodic = allocate_matrix(dh, ch.periodic), ); # We define the fourth order elasticity tensor for the matrix material, and define the diff --git a/docs/src/literate-tutorials/dg_heat_equation.jl b/docs/src/literate-tutorials/dg_heat_equation.jl index f5c2495837..1e6bd55ceb 100644 --- a/docs/src/literate-tutorials/dg_heat_equation.jl +++ b/docs/src/literate-tutorials/dg_heat_equation.jl @@ -167,7 +167,7 @@ close!(dh); # However, when generating the sparsity pattern we need to pass the topology and the cross-element coupling matrix when we're using # discontinuous interpolations. The cross-element coupling matrix is of size [1,1] in this case as # we have only one field and one DofHandler. -K = create_sparsity_pattern(dh, topology = topology, cross_coupling = trues(1,1)); +K = allocate_matrix(dh, topology = topology, interface_coupling = trues(1,1)); # ### Boundary conditions # The Dirichlet boundary conditions are treated diff --git a/docs/src/literate-tutorials/heat_equation.jl b/docs/src/literate-tutorials/heat_equation.jl index b924ce653a..639472162e 100644 --- a/docs/src/literate-tutorials/heat_equation.jl +++ b/docs/src/literate-tutorials/heat_equation.jl @@ -71,9 +71,9 @@ add!(dh, :u, ip) 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 +# using `allocate_matrix`. This function returns a sparse matrix # with the correct entries stored. -K = create_sparsity_pattern(dh) +K = allocate_matrix(dh) # ### Boundary conditions # In Ferrite constraints like Dirichlet boundary conditions diff --git a/docs/src/literate-tutorials/hyperelasticity.jl b/docs/src/literate-tutorials/hyperelasticity.jl index 9f63a6638a..9195a7374e 100644 --- a/docs/src/literate-tutorials/hyperelasticity.jl +++ b/docs/src/literate-tutorials/hyperelasticity.jl @@ -374,14 +374,14 @@ function solve() apply!(un, dbcs) ## Create sparse matrix and residual vector - K = create_sparsity_pattern(dh) + K = allocate_matrix(dh) g = zeros(_ndofs) ## Perform Newton iterations newton_itr = -1 NEWTON_TOL = 1e-8 NEWTON_MAXITER = 30 - prog = ProgressMeter.ProgressThresh(NEWTON_TOL, "Solving:") + prog = ProgressMeter.ProgressThresh(NEWTON_TOL; desc = "Solving:") while true; newton_itr += 1 ## Construct the current guess diff --git a/docs/src/literate-tutorials/incompressible_elasticity.jl b/docs/src/literate-tutorials/incompressible_elasticity.jl index ea5b097c22..5060ff86bd 100644 --- a/docs/src/literate-tutorials/incompressible_elasticity.jl +++ b/docs/src/literate-tutorials/incompressible_elasticity.jl @@ -258,7 +258,7 @@ function solve(ν, interpolation_u, interpolation_p) cellvalues_u, cellvalues_p, facetvalues_u = create_values(interpolation_u, interpolation_p) ## Assembly and solve - K = create_sparsity_pattern(dh) + K = allocate_matrix(dh) K, f = doassemble(cellvalues_u, cellvalues_p, facetvalues_u, K, grid, dh, mp) apply!(K, f, dbc) u = K \ f diff --git a/docs/src/literate-tutorials/linear_shell.jl b/docs/src/literate-tutorials/linear_shell.jl index cbddd7b9b1..b6236b9bb1 100644 --- a/docs/src/literate-tutorials/linear_shell.jl +++ b/docs/src/literate-tutorials/linear_shell.jl @@ -84,7 +84,7 @@ data = (thickness = 1.0, C = C); #Named tuple nnodes = getnbasefunctions(ip) ndofs_shell = ndofs_per_cell(dh) -K = create_sparsity_pattern(dh) +K = allocate_matrix(dh) f = zeros(Float64, ndofs(dh)) ke = zeros(ndofs_shell, ndofs_shell) diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl index 2a58297708..28db2ea391 100644 --- a/docs/src/literate-tutorials/ns_vs_diffeq.jl +++ b/docs/src/literate-tutorials/ns_vs_diffeq.jl @@ -346,10 +346,10 @@ if IS_CI #hide end #hide Δt_save = 0.1 -M = create_sparsity_pattern(dh); +M = allocate_matrix(dh); M = assemble_mass_matrix(cellvalues_v, cellvalues_p, M, dh); -K = create_sparsity_pattern(dh); +K = allocate_matrix(dh); K = assemble_stokes_matrix(cellvalues_v, cellvalues_p, ν, K, dh); # These are our initial conditions. We start from the zero solution, because it diff --git a/docs/src/literate-tutorials/plasticity.jl b/docs/src/literate-tutorials/plasticity.jl index 6b34be0e2e..acce3ecfe2 100644 --- a/docs/src/literate-tutorials/plasticity.jl +++ b/docs/src/literate-tutorials/plasticity.jl @@ -282,7 +282,7 @@ function solve() u = zeros(n_dofs) # solution vector Δu = zeros(n_dofs) # displacement correction r = zeros(n_dofs) # residual - K = create_sparsity_pattern(dh); # tangent stiffness matrix + K = allocate_matrix(dh); # tangent stiffness matrix ## Create material states. One array for each cell, where each element is an array of material- ## states - one for each integration point diff --git a/docs/src/literate-tutorials/porous_media.jl b/docs/src/literate-tutorials/porous_media.jl index 7e8a31b510..65c211dd43 100644 --- a/docs/src/literate-tutorials/porous_media.jl +++ b/docs/src/literate-tutorials/porous_media.jl @@ -331,7 +331,7 @@ end; # Given the `DofHandler`, `ConstraintHandler`, and `CellValues`, # we can solve the problem by stepping through the time history function solve(dh, ch, domains; Δt=0.025, t_total=1.0) - K = create_sparsity_pattern(dh); + K = allocate_matrix(dh) r = zeros(ndofs(dh)) a = zeros(ndofs(dh)) a_old = copy(a) diff --git a/docs/src/literate-tutorials/reactive_surface.jl b/docs/src/literate-tutorials/reactive_surface.jl index 9eb15808aa..781a3d7751 100644 --- a/docs/src/literate-tutorials/reactive_surface.jl +++ b/docs/src/literate-tutorials/reactive_surface.jl @@ -261,8 +261,8 @@ function gray_scott_on_sphere(material::GrayScottMaterial, Δt::Real, T::Real, r close!(dh); ## We can save some memory by telling the sparsity pattern that the matrices are not coupled. - M = create_sparsity_pattern(dh; coupling=[true false;false true]) - D = create_sparsity_pattern(dh; coupling=[true false;false true]) + M = allocate_matrix(dh; coupling=[true false; false true]) + D = allocate_matrix(dh; coupling=[true false; false true]) ## Since the heat problem is linear and has no time dependent parameters, we precompute the ## decomposition of the system matrix to speed up the linear system solver. diff --git a/docs/src/literate-tutorials/stokes-flow.jl b/docs/src/literate-tutorials/stokes-flow.jl index 0971fd1b4f..da3bf26dde 100644 --- a/docs/src/literate-tutorials/stokes-flow.jl +++ b/docs/src/literate-tutorials/stokes-flow.jl @@ -505,7 +505,7 @@ function main() ch = setup_constraints(dh, fvp) ## Global tangent matrix and rhs coupling = [true true; true false] # no coupling between pressure test/trial functions - K = create_sparsity_pattern(dh, ch; coupling=coupling) + K = allocate_matrix(dh, ch; coupling=coupling) f = zeros(ndofs(dh)) ## Assemble system assemble_system!(K, f, dh, cvu, cvp) diff --git a/docs/src/literate-tutorials/transient_heat_equation.jl b/docs/src/literate-tutorials/transient_heat_equation.jl index df64dc7913..30042f25c3 100644 --- a/docs/src/literate-tutorials/transient_heat_equation.jl +++ b/docs/src/literate-tutorials/transient_heat_equation.jl @@ -80,8 +80,8 @@ close!(dh); # M_{ij} = \int_{\Omega} v_i \, u_j \ \mathrm{d}\Omega, # ``` # where $u_i$ and $v_j$ are trial and test functions, respectively. -K = create_sparsity_pattern(dh); -M = create_sparsity_pattern(dh); +K = allocate_matrix(dh); +M = allocate_matrix(dh); # We also preallocate the right hand side f = zeros(ndofs(dh)); diff --git a/docs/src/reference/assembly.md b/docs/src/reference/assembly.md index 79a2a1576e..0a2111aca6 100644 --- a/docs/src/reference/assembly.md +++ b/docs/src/reference/assembly.md @@ -9,8 +9,3 @@ start_assemble assemble! finish_assemble ``` - -```@docs -create_sparsity_pattern -create_symmetric_sparsity_pattern -``` diff --git a/docs/src/reference/sparsity_pattern.md b/docs/src/reference/sparsity_pattern.md new file mode 100644 index 0000000000..a6648896b5 --- /dev/null +++ b/docs/src/reference/sparsity_pattern.md @@ -0,0 +1,57 @@ +# Sparsity pattern and sparse matrices + +This is the reference documentation for sparsity patterns and sparse matrix instantiation. +See the topic section on [Sparsity pattern and sparse matrices](@ref topic-sparse-matrix). + +## Sparsity patterns + +### `AbstractSparsityPattern` + +The following applies to all subtypes of `AbstractSparsityPattern`: + +```@docs +Ferrite.AbstractSparsityPattern +add_sparsity_entries! +add_cell_entries! +add_interface_entries! +add_constraint_entries! +Ferrite.add_entry! +``` + +### `SparsityPattern` + +```@docs +SparsityPattern(::Int, ::Int) +allocate_matrix(::SparsityPattern) +SparsityPattern +``` + +### `BlockSparsityPattern` + +!!! note "Package extension" + This functionality is only enabled when the package + [BlockArrays.jl](https://github.com/JuliaArrays/BlockArrays.jl) is installed (`pkg> add + BlockArrays`) and loaded (`using BlockArrays`) in the session. + +```@docs +BlockSparsityPattern(::Vector{Int}) +Main.FerriteBlockArrays.BlockSparsityPattern +allocate_matrix(::Main.FerriteBlockArrays.BlockSparsityPattern) +allocate_matrix(::Type{<:BlockMatrix{T, Matrix{S}}}, sp::Main.FerriteBlockArrays.BlockSparsityPattern) where {T, S <: AbstractMatrix{T}} +``` + +## Sparse matrices + +### Creating matrix from `SparsityPattern` + +```@docs +allocate_matrix(::Type{S}, ::Ferrite.AbstractSparsityPattern) where {Tv, Ti, S <: SparseMatrixCSC{Tv, Ti}} +allocate_matrix(::Type{Symmetric{Tv, S}}, ::Ferrite.AbstractSparsityPattern) where {Tv, Ti, S <: SparseMatrixCSC{Tv, Ti}} +``` + +### Creating matrix from `DofHandler` + +```@docs +allocate_matrix(::Type{MatrixType}, ::DofHandler, args...; kwargs...) where {MatrixType} +allocate_matrix(::DofHandler, args...; kwargs...) +``` diff --git a/docs/src/topics/assembly.md b/docs/src/topics/assembly.md index 90d2e506a7..0a7d1fb8ce 100644 --- a/docs/src/topics/assembly.md +++ b/docs/src/topics/assembly.md @@ -27,32 +27,6 @@ e.g. iterative solvers or time dependent problems where the sparse matrix structure, or [Sparsity Pattern](@ref) will stay the same in every iteration/ time step. -## Sparsity Pattern - -Given a `DofHandler` we can obtain the corresponding sparse matrix by using the -[`create_sparsity_pattern`](@ref) function. This will setup a `SparseMatrixCSC` -with stored values on all the places corresponding to the degree of freedom numbering -in the `DofHandler`. This means that when we assemble into the global stiffness matrix -there is no need to change the internal representation of the sparse matrix since the -sparse structure will not change. - -Often the finite element problem is symmetric and will result in a symmetric sparse -matrix. This information is often something that the sparse solver can take advantage of. -If the solver only needs half the matrix there is no need to assemble both halves. -For this purpose there is a [`create_symmetric_sparsity_pattern`](@ref) function that -will only create the upper half of the matrix, and return a `Symmetric` wrapped -`SparseMatrixCSC`. - -Given a `DofHandler` `dh` we can obtain the (symmetric) sparsity pattern as - -```julia -K = create_sparsity_pattern(dh) -K = create_symmetric_sparsity_pattern(dh) -``` - -The returned sparse matrix will be used together with an `Assembler`, which -assembles efficiently into the matrix, without modifying the internal representation. - ## `Assembler` Assembling efficiently into the sparse matrix requires some extra workspace. @@ -96,7 +70,7 @@ In such cases it is enough to construct the global matrix `K` once. Below is some pseudo-code for how to do this for a time-dependent problem: ```julia -K = create_sparsity_pattern(dh) +K = allocate_matrix(dh) f = zeros(ndofs(dh)) for t in 1:timesteps diff --git a/docs/src/topics/constraints.md b/docs/src/topics/constraints.md index ac2d39a19f..c2baa243fa 100644 --- a/docs/src/topics/constraints.md +++ b/docs/src/topics/constraints.md @@ -39,7 +39,7 @@ Affine constraints will affect the sparsity pattern of the stiffness matrix, and the `ConstraintHandler` as an argument when creating the sparsity pattern: ```julia -K = create_sparsity_pattern(dh, ch) +K = allocate_matrix(dh, ch) ``` ### Solving linear problems diff --git a/docs/src/topics/sparse_matrix.md b/docs/src/topics/sparse_matrix.md new file mode 100644 index 0000000000..923315ed4a --- /dev/null +++ b/docs/src/topics/sparse_matrix.md @@ -0,0 +1,184 @@ +# [Sparsity pattern and sparse matrices](@id topic-sparse-matrix) + +An important property of the finite element method is that it results in *sparse matrices* +for the linear systems to be solved. On this page the topic of sparsity and sparse matrices +are discussed. + +```@contents +Pages = ["sparse_matrix.md"] +Depth = 2:2 +``` + +## Sparsity pattern + +The sparse structure of the linear system depends on many factors such as e.g. the weak +form, the discretization, and the choice of interpolation(s). In the end it boils down to +how the degrees of freedom (DoFs) *couple* with each other. The most common reason that two +DoFs couple is because they belong to the same element. Note, however, that this is not +guaranteed to result in a coupling since it depends on the specific weak form that is being +discretized, see e.g. [Increasing the sparsity](@ref). Boundary conditions and constraints +can also result in additional DoF couplings. + +If DoFs `i` and `j` couple, then the computed value in the eventual matrix will be +*structurally nonzero*[^1]. In this case the entry `(i, j)` should be included in the +sparsity pattern. Conversely, if DoFs `i` and `j` *don't* couple, then the computed value +will be *zero*. In this case the entry `(i, j)` should *not* be included in the sparsity +pattern since there is no need to allocate memory for entries that will be zero. + +The sparsity, i.e. the ratio of zero-entries to the total number of entries, is often[^2] +*very* high and taking advantage of this results in huge savings in terms of memory. For +example, in a problem with ``10^6`` DoFs there will be a matrix of size ``10^6 \times +10^6``. If all ``10^{12}`` entries of this matrix had to be stored (0% sparsity) as double +precision (`Float64`, 8 bytes) it would require 8 TB of memory. If instead the sparsity is +99.9973% (which is the case when solving the heat equation on a three dimensional hypercube +with linear Lagrange interpolation) this would be reduced to 216 MB. + +[1]: Structurally nonzero means that there is a possibility of a nonzero value even though + the computed value might become zero in the end for various reasons. + +[2]: At least for most practical problems using low order interpolations. + + +!!! details "Sparsity pattern example" + + To give an example, in this one-dimensional heat problem (see the [Heat + equation](../tutorials/heat_equation.md) tutorial for the weak form) we have 4 nodes + with 3 elements in between. For simplicitly DoF numbers and node numbers are the same + but this is not true in general since nodes and DoFs can be numbered independently (and + in fact are numbered independently in Ferrite). + + ``` + 1 ----- 2 ----- 3 ----- 4 + ``` + + Assuming we use linear Lagrange interpolation (the "hat functions") this will give the + following connections according to the weak form: + - Trial function 1 couples with test functions 1 and 2 (entries `(1, 1)` and `(1, 2)` + included in the sparsity pattern) + - Trial function 2 couples with test functions 1, 2, and 3 (entries `(2, 1)`, `(2, 2)`, + and `(2, 3)` included in the sparsity pattern) + - Trial function 3 couples with test functions 2, 3, and 4 (entries `(3, 2)`, `(3, 3)`, + and `(3, 4)` included in the sparsity pattern) + - Trial function 4 couples with test functions 3 and 4 (entries `(4, 3)` and `(4, 4)` + included in the sparsity pattern) + + The resulting sparsity pattern would look like this: + + ``` + 4×4 SparseArrays.SparseMatrixCSC{Float64, Int64} with 10 stored entries: + 0.0 0.0 ⋅ ⋅ + 0.0 0.0 0.0 ⋅ + ⋅ 0.0 0.0 0.0 + ⋅ ⋅ 0.0 0.0 + ``` + + Moreover, if the problem is solved with periodic boundary conditions, for example by + constraining the value on the right side to the value on the left side, there will be + additional couplings. In the example above, this means that DoF 4 should be equal to DoF + 1. Since DoF 4 is constrained it has to be eliminated from the system. Existing entries + that include DoF 4 are `(3, 4)`, `(4, 3)`, and `(4, 4)`. Given the simple constraint in + this case we can simply replace DoF 4 with DoF 1 in these entries and we end up with + entries `(3, 1)`, `(1, 3)`, and `(1, 1)`. This results in two new entries: `(3, 1)` and + `(1, 3)` (entry `(1, 1)` is already included). + +## Creating sparsity patterns + +Creating a sparsity pattern can be quite expensive if not done properly and therefore +Ferrite provides efficient methods and data structures for this. In general the sparsity +pattern is not known in advance and has to be created incrementally. To make this +incremental construction efficient it is necessary to use a dynamic data structure which +allow for fast insertions. + +The sparsity pattern also serves as a "matrix builder". When all entries are inserted into +the sparsity pattern the dynamic data structure is typically converted, or "compressed", +into a sparse matrix format such as e.g. the [*compressed sparse row +(CSR)*](https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_row_(CSR,_CRS_or_Yale_format)) +format or the [*compressed sparse column +(CSC)*](https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_column_(CSC_or_CCS)) +format, where the latter is the default sparse matrix type implemented in the [SparseArrays +standard library](https://github.com/JuliaSparse/SparseArrays.jl). These matrix formats +allow for fast linear algebra operations, such as factorizations and matrix-vector +multiplications, that are needed when the linear system is solved. See [Instantiating the +sparse matrix](@ref) for more details. + +In summary, a dynamic structure is more efficient when incrementally building the pattern by +inserting new entries, and a static or compressed structure is more efficient for linear +algebra operations. + +### Basic sparsity patterns construction + +Working with the sparsity pattern explicitly is in many cases not necessary. For basic +usage (e.g. when only one matrix needed, when no customization of the pattern is +required, etc) there exist convenience methods of [`allocate_matrix`](@ref) that return +the matrix directly. Most examples in this documentation don't deal with the sparsity +pattern explicitly because the basic method suffice. +See also [Instantiating the sparse matrix](@ref) for more details. + +### Custom sparsity pattern construction + +In more advanced cases there might be a need for more fine grained control of the sparsity +pattern. The following steps are typically taken when constructing a sparsity pattern in +Ferrite: + + 1. **Initialize an empty pattern:** This can be done by either using the + [`init_sparsity_pattern(dh)`](@ref) function or by using a constructor directly. + `init_sparsity_pattern` will return a default pattern type that is compatible with the + DofHandler. In some cases you might require another type of pattern (for example a + blocked pattern, see [Blocked sparsity pattern](@ref)) and in that case you can use the + constructor directly. + + 2. **Add entries to the pattern:** There are a number of functions that add entries to the + pattern: + - [`add_sparsity_entries!`](@ref) is a convenience method for performing the common + task of calling `add_cell_entries!`, `add_interface_entries!`, and + `add_constraint_entries!` after each other (see below). + - [`add_cell_entries!`](@ref) adds entries for all couplings between the DoFs within + each element. These entries correspond to assembling the standard element matrix and + is thus almost always required. + - [`add_interface_entries!`](@ref) adds entries for couplings between the DoFs in + neighboring elements. These entries are required when integrating along internal + interfaces between elements (e.g. for discontinuous Galerkin methods). + - [`add_constraint_entries!`](@ref) adds entries required from constraints and boundary + conditions in the ConstraintHandler. Note that this operation depends on existing + entries in the pattern and *must* be called as the last operation on the pattern. + - [`Ferrite.add_entry!`](@ref) adds a single entry to the pattern. This can be used if + you need to add custom entries that are not covered by the other functions. + + 3. **Instantiate the matrix:** A sparse matrix can be created from the sparsity pattern + using [`allocate_matrix`](@ref), see [Instantiating the sparse matrix](@ref) below for + more details. + +### Increasing the sparsity + +By default, when creating a sparsity pattern, it is assumed that each DoF within an element +couple with with *all* other DoFs in the element. + +!!! todo + - Discuss the `coupling` keyword argument. + - Discuss the `keep_constrained` keyword argument. + +### Blocked sparsity pattern + +!!! todo + Discuss `BlockSparsityPattern` and `BlockArrays` extension. + +## Instantiating the sparse matrix + +As mentioned above, for many simple cases there is no need to work with the sparsity pattern +directly and using methods of [`allocate_matrix`](@ref) that take the DofHandler as input is +enough, for example: + +```julia +K = allocate_matrix(dh, ch) +``` + +`allocate_matrix` is also used to instantiate a matrix from a sparsity pattern, for example: + +```julia +K = allocate_matrix(sp) +``` + +!!! note "Multiple matrices with the same pattern" + For some problems there is a need for multiple matrices with the same sparsity pattern, + for example a mass matrix and a stiffness matrix. In this case it is more efficient to + create the sparsity pattern once and then instantiate both matrices from it. diff --git a/ext/FerriteBlockArrays.jl b/ext/FerriteBlockArrays.jl index 42f46a1f00..b7bdb9ef7c 100644 --- a/ext/FerriteBlockArrays.jl +++ b/ext/FerriteBlockArrays.jl @@ -2,60 +2,61 @@ module FerriteBlockArrays using BlockArrays: BlockArray, BlockIndex, BlockMatrix, BlockVector, block, blockaxes, blockindex, blocks, findblockindex +using BlockArrays: Block, BlockArray, BlockIndex, BlockMatrix, BlockVector, block, + blockaxes, blockindex, blocks, findblockindex, undef_blocks using Ferrite: - Ferrite, CellIterator, ConstraintHandler, DofHandler, addindex!, assemble!, celldofs, - create_sparsity_pattern, dof_range, fillzero!, ndofs - -# TODO: Move into Ferrite and enable for mixed grids / subdomains -function global_dof_range(dh::DofHandler, f::Symbol) - set = Set{Int}() - frange = dof_range(dh, f) - for cc in CellIterator(dh) - union!(set, @view celldofs(cc)[frange]) - end - dofmin, dofmax = extrema(set) - r = dofmin:dofmax - if length(set) != length(r) - error("renumber by blocks you donkey") - end - return r -end + Ferrite, BlockSparsityPattern, ConstraintHandler, addindex!, allocate_matrix, assemble!, + fillzero! +using SparseArrays: SparseMatrixCSC + + +############################## +## Instantiating the matrix ## +############################## -################################### -## Creating the sparsity pattern ## -################################### - -# Note: -# Creating the full unblocked matrix and then splitting into blocks inside the BlockArray -# constructor (i.e. by `getindex(::SparseMatrixCSC, ::UnitRange, ::UnitRange)`) is -# consistently faster than creating individual blocks directly. However, the latter approach -# uses less than half of the memory (measured for a 2x2 block system and various problem -# sizes), so might be useful in the future to provide an option on what algorithm to use. - -# TODO: Could potentially extract the element type and matrix type for the individual blocks -# by allowing e.g. create_sparsity_pattern(BlockMatrix{Float32}, ...) but that is not -# even supported by regular pattern right now. -function Ferrite.create_sparsity_pattern(::Type{<:BlockMatrix}, dh, ch; kwargs...) - K = create_sparsity_pattern(dh, ch; kwargs...) - # Infer block sizes from the fields in the DofHandler - block_sizes = [length(global_dof_range(dh, f)) for f in dh.field_names] - return BlockArray(K, block_sizes, block_sizes) +# function Ferrite.allocate_matrix(::Type{B}, dh, ch, ...) where B <: BlockMatrix +# # TODO: Create BSP from the induced field blocks in dh +# end + +# Fill in missing matrix type, this allows allocate_matrix(BlockMatrix, sp) +function Ferrite.allocate_matrix(::Type{<:BlockMatrix}, sp::BlockSparsityPattern) + return allocate_matrix(BlockMatrix{Float64, Matrix{SparseMatrixCSC{Float64, Int}}}, sp) end -function Ferrite.create_sparsity_pattern(B::BlockMatrix, dh, ch; kwargs...) - if !(size(B, 1) == size(B, 2) == ndofs(dh)) - error("size of input matrix ($(size(B))) does not match number of dofs ($(ndofs(dh)))") - end - K = create_sparsity_pattern(dh, ch; kwargs...) - ax = axes(B) - for block_j in blockaxes(B, 2), block_i in blockaxes(B, 1) - range_j = ax[2][block_j] - range_i = ax[1][block_i] - B[block_i, block_j] = K[range_i, range_j] +""" + 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 type +allocate_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) +``` + +!!! note "Package extension" + This functionality is only enabled when the package + [BlockArrays.jl](https://github.com/JuliaArrays/BlockArrays.jl) is installed (`pkg> add + BlockArrays`) and loaded (`using BlockArrays`) in the session. +""" +function Ferrite.allocate_matrix(::Type{<:BlockMatrix{T, Matrix{S}}}, sp::BlockSparsityPattern) where {T, S <: AbstractMatrix{T}} + @assert isconcretetype(S) + block_sizes = sp.block_sizes + K = BlockArray(undef_blocks, S, block_sizes, block_sizes) + for j in 1:length(block_sizes), i in 1:length(block_sizes) + K[Block(i), Block(j)] = allocate_matrix(S, sp.blocks[i, j]) end - return B + return K end + ########################################### ## BlockAssembler and associated methods ## ########################################### diff --git a/ext/FerriteMetis.jl b/ext/FerriteMetis.jl index b1f9b47193..053d9a725d 100644 --- a/ext/FerriteMetis.jl +++ b/ext/FerriteMetis.jl @@ -41,7 +41,7 @@ function Ferrite.compute_renumber_permutation( if coupling !== nothing # Set sym = true since Metis.permutation requires a symmetric graph. # TODO: Perhaps just symmetrize it: coupling = coupling' .| coupling - couplings = Ferrite._coupling_to_local_dof_coupling(dh, coupling, #= sym =# true) + couplings = Ferrite._coupling_to_local_dof_coupling(dh, coupling) end # Create the CSR (CSC, but pattern is symmetric so equivalent) using diff --git a/src/Dofs/block_sparsity_pattern.jl b/src/Dofs/block_sparsity_pattern.jl new file mode 100644 index 0000000000..465f5f5223 --- /dev/null +++ b/src/Dofs/block_sparsity_pattern.jl @@ -0,0 +1,131 @@ +######################## +# BlockSparsityPattern # +######################## + +""" + struct BlockSparsityPattern <: AbstractSparsityPattern + +Data structure representing non-zero entries for an eventual *blocked* sparse matrix. + +See the constructor [`BlockSparsityPattern(::Vector{Int})`](@ref +BlockSparsityPattern(::Vector{Int})) for the user-facing documentation. + +# Struct fields + - `nrows::Int`: number of rows + - `ncols::Int`: number of column + - `block_sizes::Vector{Int}`: row and column block sizes + - `blocks::Matrix{SparsityPattern}`: matrix of size `length(block_sizes) × + length(block_sizes)` where `blocks[i, j]` is a [`SparsityPattern`](@ref) corresponding to + block `(i, j)`. + +!!! warning "Internal struct" + The specific implementation of this struct, such as struct fields, type layout and type + parameters, are internal and should not be relied upon. +""" +struct BlockSparsityPattern <: AbstractSparsityPattern + nrows::Int + ncols::Int + block_sizes::Vector{Int} + blocks::Matrix{SparsityPattern} +end + +""" + BlockSparsityPattern(block_sizes::Vector{Int}) + +Create an empty `BlockSparsityPattern` with row and column block sizes given by +`block_sizes`. + +# Examples +```julia +# Create a block sparsity pattern with block size 10 x 5 +sparsity_pattern = BlockSparsityPattern([10, 5]) +``` + +# Methods +The following methods apply to `BlockSparsityPattern` (see their respective documentation +for more details): + + - [`add_sparsity_entries!`](@ref): convenience method for calling + [`add_cell_entries!`](@ref), [`add_interface_entries!`](@ref), and + [`add_constraint_entries!`](@ref). + - [`add_cell_entries!`](@ref): add entries corresponding to DoF couplings within the cells. + - [`add_interface_entries!`](@ref): add entries corresponding to DoF couplings on the + interface between cells. + - [`add_constraint_entries!`](@ref): add entries resulting from constraints. + - [`allocate_matrix`](@ref allocate_matrix(::SparsityPattern)): instantiate a (block) + matrix from the pattern. The default matrix type is `BlockMatrix{Float64, + Matrix{SparseMatrixCSC{Float64, Int}}}`, i.e. a `BlockMatrix`, where the individual + blocks are of type `SparseMatrixCSC{Float64, Int}`. + +!!! note "Package extension" + This functionality is only enabled when the package + [BlockArrays.jl](https://github.com/JuliaArrays/BlockArrays.jl) is installed (`pkg> add + BlockArrays`) and loaded (`using BlockArrays`) in the session. +""" +function BlockSparsityPattern(blk_sizes::AbstractVector{<:Integer}) + block_sizes = collect(Int, blk_sizes) + nrows = ncols = sum(block_sizes) + nblocks = length(block_sizes) + # TODO: Maybe all of these could/should share the same PoolAllocator? + blocks = [SparsityPattern(block_sizes[i], block_sizes[j]) for i in 1:nblocks, j in 1:nblocks] + return BlockSparsityPattern(nrows, ncols, block_sizes, blocks) +end + +getnrows(bsp::BlockSparsityPattern) = bsp.nrows +getncols(bsp::BlockSparsityPattern) = bsp.ncols + +# Compute block index and local index into that block +@inline function _find_block(block_sizes::Vector{Int}, i::Int) + accumulated_block_size = 0 + block_index = 1 + while !(accumulated_block_size < i <= accumulated_block_size + block_sizes[block_index]) + accumulated_block_size += block_sizes[block_index] + block_index += 1 + end + local_index = i - accumulated_block_size + return block_index, local_index +end + +@inline function add_entry!(bsp::BlockSparsityPattern, row::Int, col::Int) + @boundscheck 1 <= row <= getnrows(bsp) && 1 <= col <= getncols(bsp) + row_block, row_local = _find_block(bsp.block_sizes, row) + col_block, col_local = _find_block(bsp.block_sizes, col) + add_entry!(bsp.blocks[row_block, col_block], row_local, col_local) + return +end + +# Helper struct to iterate over the rows. Behaves similar to +# Iterators.flatten([eachrow(bsp.blocks[row_block, col_block) for col_block in 1:nblocks]) +# but we need to add the offset to the iterated values. +struct BSPRowIterator + bsp::BlockSparsityPattern + row::Int + row_block::Int + row_local::Int + function BSPRowIterator(bsp::BlockSparsityPattern, row::Int) + @assert 1 <= row <= getnrows(bsp) + row_block, row_local = _find_block(bsp.block_sizes, row) + return new(bsp, row, row_block, row_local) + end +end + +function Base.iterate(it::BSPRowIterator, state = (1, 1)) + col_block, idx = state + bsp = it.bsp + col_block > length(bsp.block_sizes) && return nothing + block = bsp.blocks[it.row_block, col_block] + colidxs = eachrow(block, it.row_local) + if idx > length(colidxs) + # Advance the col_block and reset idx to 1 + return iterate(it, (col_block + 1, 1)) + else + # Compute global col idx and advance idx + col_local = colidxs[idx] + offset = sum((bsp.block_sizes[i] for i in 1:col_block-1); init = 0) + return offset + col_local, (col_block, idx + 1) + end +end + +# TODO: eltype of the generator do not infer; might need another auxiliary struct. +eachrow(bsp::BlockSparsityPattern) = (BSPRowIterator(bsp, row) for row in 1:getnrows(bsp)) +eachrow(bsp::BlockSparsityPattern, row::Int) = BSPRowIterator(bsp, row) diff --git a/src/Dofs/sparsity_pattern.jl b/src/Dofs/sparsity_pattern.jl index a4bbdd8a22..fd51dacfa7 100644 --- a/src/Dofs/sparsity_pattern.jl +++ b/src/Dofs/sparsity_pattern.jl @@ -1,72 +1,449 @@ +########################### +# AbstractSparsityPattern # +########################### + +""" + Ferrite.AbstractSparsityPattern + +Supertype for sparsity pattern implementations, e.g. [`SparsityPattern`](@ref) and +[`BlockSparsityPattern`](@ref). +""" +abstract type AbstractSparsityPattern end + +""" + getnrows(sp::AbstractSparsityPattern) + +Return the number of rows in the sparsity pattern `sp`. +""" +getnrows(sp::AbstractSparsityPattern) + +""" + getncols(sp::AbstractSparsityPattern) + +Return the number of columns in the sparsity pattern `sp`. +""" +getncols(sp::AbstractSparsityPattern) + +""" + add_entry!(sp::AbstractSparsityPattern, row::Int, col::Int) + +Add an entry to the sparsity pattern `sp` at row `row` and column `col`. +""" +add_entry!(sp::AbstractSparsityPattern, row::Int, col::Int) + +# This is necessary to avoid warning about not importing Base.eachrow when +# adding docstring before the definitions further down. +function eachrow end + +""" + eachrow(sp::AbstractSparsityPattern) + +Return an iterator over the rows of the sparsity pattern `sp`. +Each element of the iterator iterates indices of the stored *columns* for that row. +""" +eachrow(sp::AbstractSparsityPattern) + """ - create_sparsity_pattern(dh::DofHandler; coupling, [topology::Union{Nothing, AbstractTopology}], [cross_coupling]) + eachrow(sp::AbstractSparsityPattern, row::Int) + +Return an iterator over *column* indices in row `row` of the sparsity pattern. -Create the sparsity pattern corresponding to the degree of freedom -numbering in the [`DofHandler`](@ref). Return a `SparseMatrixCSC` -with stored values in the correct places. +Conceptually this is equivalent to [`eachrow(sp)[row]`](@ref +eachrow(::AbstractSparsityPattern)). However, the iterator `eachrow(sp)` isn't always +indexable. This method should be used when a specific row needs to be "random access"d. +""" +eachrow(sp::AbstractSparsityPattern, row::Int) -The keyword arguments `coupling` and `cross_coupling` can be used to specify how fields (or components) in the dof -handler couple to each other. `coupling` and `cross_coupling` should be square matrices of booleans with -number of rows/columns equal to the total number of fields, or total number of components, -in the DofHandler with `true` if fields are coupled and `false` if -not. By default full coupling is assumed inside the element with no coupling between elements. -If `topology` and `cross_coupling` are passed, dof of fields with discontinuous interpolations are coupled between elements according to `cross_coupling`. +################### +# SparsityPattern # +################### -See the [Sparsity Pattern](@ref) section of the manual. """ -function create_sparsity_pattern(dh::AbstractDofHandler; coupling=nothing, - topology::Union{Nothing, AbstractTopology} = nothing, cross_coupling = nothing) - return _create_sparsity_pattern(dh, nothing, false, true, coupling, topology, cross_coupling) + struct SparsityPattern <: AbstractSparsityPattern + +Data structure representing non-zero entries in the eventual sparse matrix. + +See the constructor [`SparsityPattern(::Int, ::Int)`](@ref) for the user-facing +documentation. + +# Struct fields + - `nrows::Int`: number of rows + - `ncols::Int`: number of column + - `rows::Vector{Vector{Int}}`: vector of length `nrows`, where `rows[i]` is a + *sorted* vector of column indices for non zero entries in row `i`. + +!!! warning "Internal struct" + The specific implementation of this struct, such as struct fields, type layout and type + parameters, are internal and should not be relied upon. +""" +struct SparsityPattern <: AbstractSparsityPattern + nrows::Int + ncols::Int + mempool::PoolAllocator.MemoryPool{Int} + rows::Vector{PoolAllocator.PoolVector{Int}} end """ - create_symmetric_sparsity_pattern(dh::DofHandler; coupling, topology::Union{Nothing, AbstractTopology}, cross_coupling) + SparsityPattern(nrows::Int, ncols::Int; nnz_per_row::Int = 8) + +Create an empty [`SparsityPattern`](@ref) with `nrows` rows and `ncols` columns. +`nnz_per_row` is used as a memory hint for the number of non zero entries per +row. + +`SparsityPattern` is the default sparsity pattern type for the standard DofHandler and is +therefore commonly constructed using [`init_sparsity_pattern`](@ref) instead of with this +constructor. + +# Examples +```julia +# Create a sparsity pattern for an 100 x 100 matrix, hinting at 10 entries per row +sparsity_pattern = SparsityPattern(100, 100; nnz_per_row = 10) +``` + +# Methods +The following methods apply to `SparsityPattern` (see their respective documentation for +more details): + - [`add_sparsity_entries!`](@ref): convenience method for calling + [`add_cell_entries!`](@ref), [`add_interface_entries!`](@ref), and + [`add_constraint_entries!`](@ref). + - [`add_cell_entries!`](@ref): add entries corresponding to DoF couplings within the cells. + - [`add_interface_entries!`](@ref): add entries corresponding to DoF couplings on the + interface between cells. + - [`add_constraint_entries!`](@ref): add entries resulting from constraints. + - [`allocate_matrix`](@ref allocate_matrix(::SparsityPattern)): instantiate a matrix from + the pattern. The default matrix type is `SparseMatrixCSC{Float64, Int}`. +""" +function SparsityPattern(nrows::Int, ncols::Int; nnz_per_row::Int = 8) + mempool = PoolAllocator.MemoryPool{Int}() + rows = Vector{PoolAllocator.PoolVector{Int}}(undef, nrows) + for i in 1:nrows + rows[i] = PoolAllocator.resize(PoolAllocator.malloc(mempool, nnz_per_row), 0) + end + sp = SparsityPattern(nrows, ncols, mempool, rows) + return sp +end + +function Base.show(io::IO, ::MIME"text/plain", sp::SparsityPattern) + iob = IOBuffer() + println(iob, "$(getnrows(sp))×$(getncols(sp)) $(sprint(show, typeof(sp))):") + # Collect min/max/avg entries per row + min_entries = typemax(Int) + max_entries = typemin(Int) + stored_entries = 0 + for r in eachrow(sp) + l = length(r) + stored_entries += l + min_entries = min(min_entries, l) + max_entries = max(max_entries, l) + end + # Print sparsity + sparsity_pct = round( + (getnrows(sp) * getncols(sp) - stored_entries) / (getnrows(sp) * getncols(sp)) * 100 * 1000 + ) / 1000 + println(iob, " - Sparsity: $(sparsity_pct)% ($(stored_entries) stored entries)") + # Print row stats + avg_entries = round(stored_entries / getnrows(sp) * 10) / 10 + println(iob, " - Entries per row (min, max, avg): $(min_entries), $(max_entries), $(avg_entries)") + # Compute memory estimate + @assert getnrows(sp) * sizeof(eltype(sp.rows)) == sizeof(sp.rows) + bytes_used = sizeof(sp.rows) + stored_entries * sizeof(Int) + bytes_allocated = sizeof(sp.rows) + PoolAllocator.mempool_stats(sp.mempool)[2] + print(iob, " - Memory estimate: $(Base.format_bytes(bytes_used)) used, $(Base.format_bytes(bytes_allocated)) allocated") + write(io, seekstart(iob)) + return +end + +getnrows(sp::SparsityPattern) = sp.nrows +getncols(sp::SparsityPattern) = sp.ncols + +@inline function add_entry!(sp::SparsityPattern, row::Int, col::Int) + @boundscheck (1 <= row <= getnrows(sp) && 1 <= col <= getncols(sp)) || throw(BoundsError(sp, (row, col))) + r = @inbounds sp.rows[row] + r = insert_sorted(r, col) + @inbounds sp.rows[row] = r + return +end + +@inline function insert_sorted(x::PoolAllocator.PoolVector{Int}, item::Int) + k = searchsortedfirst(x, item) + if k == length(x) + 1 || @inbounds(x[k]) != item + x = PoolAllocator.insert(x, k, item) + end + return x +end + +eachrow(sp::SparsityPattern) = sp.rows +eachrow(sp::SparsityPattern, row::Int) = sp.rows[row] -Create the symmetric sparsity pattern corresponding to the degree of freedom -numbering in the [`DofHandler`](@ref) by only considering the upper -triangle of the matrix. Return a `Symmetric{SparseMatrixCSC}`. -See the [Sparsity Pattern](@ref) section of the manual. +################################################ +## Adding entries to AbstractSparsityPatterns ## +################################################ + +""" + init_sparsity_pattern(dh::DofHandler; nnz_per_row::Int) + +Initialize an empty [`SparsityPattern`](@ref) with `ndofs(dh)` rows and `ndofs(dh)` columns. + +# Keyword arguments + - `nnz_per_row`: memory optimization hint for the number of non-zero entries per row that + will be added to the pattern. """ -function create_symmetric_sparsity_pattern(dh::AbstractDofHandler; coupling=nothing, - topology::Union{Nothing, AbstractTopology} = nothing, cross_coupling = nothing) - return Symmetric(_create_sparsity_pattern(dh, nothing, true, true, coupling, topology, cross_coupling), :U) +function init_sparsity_pattern( + dh::DofHandler; + # TODO: What is a good estimate for nnz_per_row? + nnz_per_row::Int = 2 * ndofs_per_cell(dh.subdofhandlers[1]), # FIXME + ) + sp = SparsityPattern(ndofs(dh), ndofs(dh); nnz_per_row = nnz_per_row) + return sp end """ - create_symmetric_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler; coupling, topology::Union{Nothing, AbstractTopology}, cross_coupling) + add_sparsity_entries!( + sp::AbstractSparsityPattern, + dh::DofHandler, + ch::Union{ConstraintHandler, Nothing} = nothing; + topology = nothing, + keep_constrained::Bool = true, + coupling = nothing, + interface_coupling = nothing, + ) + +Convenience method for doing the common task of calling [`add_cell_entries!`](@ref), +[`add_interface_entries!`](@ref), and [`add_constraint_entries!`](@ref), depending on what +arguments are passed: + - `add_cell_entries!` is always called + - `add_interface_entries!` is called if `topology` is provided (i.e. not `nothing`) + - `add_constraint_entries!` is called if the ConstraintHandler is provided + +For more details about arguments and keyword arguments, see the respective functions. +""" +function add_sparsity_entries!( + sp::AbstractSparsityPattern, dh::DofHandler, + ch::Union{ConstraintHandler, Nothing} = nothing; + keep_constrained::Bool = true, + coupling::Union{AbstractMatrix{Bool}, Nothing} = nothing, + interface_coupling::Union{AbstractMatrix{Bool}, Nothing} = nothing, + topology = nothing, + ) + # Argument checking + isclosed(dh) || error("the DofHandler must be closed") + if getnrows(sp) < ndofs(dh) || getncols(sp) < ndofs(dh) + error("number of rows ($(getnrows(sp))) or columns ($(getncols(sp))) in the sparsity pattern is smaller than number of dofs ($(ndofs(dh)))") + end + # Add all entries + add_diagonal_entries!(sp) + add_cell_entries!(sp, dh, ch; keep_constrained, coupling) + if topology !== nothing + add_interface_entries!(sp, dh, ch; topology, keep_constrained, interface_coupling) + end + if ch !== nothing + add_constraint_entries!(sp, ch; keep_constrained) + end + return sp +end -Create a symmetric sparsity pattern accounting for affine constraints in `ch`. See -the Affine Constraints section of the manual for further details. """ -function create_symmetric_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler; - keep_constrained::Bool=true, coupling=nothing, topology::Union{Nothing, AbstractTopology} = nothing, - cross_coupling = nothing) - return Symmetric(_create_sparsity_pattern(dh, ch, true, keep_constrained, coupling, topology, cross_coupling), :U) + add_cell_entries!( + sp::AbstractSparsityPattern, + dh::DofHandler, + ch::Union{ConstraintHandler, Nothing} = nothing; + keep_constrained::Bool = true, + coupling::Union{AbstractMatrix{Bool}, Nothing}, = nothing + ) + +Add entries to the sparsity pattern `sp` corresponding to DoF couplings within the cells as +described by the DofHandler `dh`. + +# Keyword arguments + - `keep_constrained`: whether or not entries for constrained DoFs should be kept + (`keep_constrained = true`) or eliminated (`keep_constrained = false`) from the sparsity + pattern. `keep_constrained = false` requires passing the ConstraintHandler `ch`. + - `coupling`: the coupling between fields/components within each cell. By default + (`coupling = nothing`) it is assumed that all DoFs in each cell couple with each other. +""" +function add_cell_entries!( + sp::AbstractSparsityPattern, + dh::DofHandler, ch::Union{ConstraintHandler, Nothing} = nothing; + keep_constrained::Bool = true, coupling::Union{AbstractMatrix{Bool}, Nothing} = nothing, + ) + # Expand coupling from nfields × nfields to ndofs_per_cell × ndofs_per_cell + # TODO: Perhaps this can be done in the loop over SubDofHandlers instead. + if coupling !== nothing + coupling = _coupling_to_local_dof_coupling(dh, coupling) + end + if !keep_constrained + ch === nothing && error("must pass ConstraintHandler when `keep_constrained = true`") + isclosed(ch) || error("the ConstraintHandler must be closed") + ch.dh === dh || error("the DofHandler and the ConstraintHandler's DofHandler must be the same") + end + return _add_cell_entries!(sp, dh, ch, keep_constrained, coupling) end """ - create_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler; coupling, topology::Union{Nothing, AbstractTopology} = nothing) + add_interface_entries!( + sp::SparsityPattern, dh::DofHandler, ch::Union{ConstraintHandler, Nothing}; + topology::ExclusiveTopology, keep_constrained::Bool = true, + interface_coupling::AbstractMatrix{Bool}, + ) + +Add entries to the sparsity pattern `sp` corresponding to DoF couplings on the interface +between cells as described by the DofHandler `dh`. + +# Keyword arguments + - `topology`: the topology corresponding to the grid. + - `keep_constrained`: whether or not entries for constrained DoFs should be kept + (`keep_constrained = true`) or eliminated (`keep_constrained = false`) from the sparsity + pattern. `keep_constrained = false` requires passing the ConstraintHandler `ch`. + - `interface_coupling`: the coupling between fields/components across the interface. +""" +function add_interface_entries!( + sp::SparsityPattern, dh::DofHandler, ch::Union{ConstraintHandler, Nothing} = nothing; + topology::ExclusiveTopology, keep_constrained::Bool = true, + interface_coupling::AbstractMatrix{Bool}, + ) + if !keep_constrained + ch === nothing && error("must pass ConstraintHandler when `keep_constrained = true`") + isclosed(ch) || error("the ConstraintHandler must be closed") + ch.dh === dh || error("the DofHandler and the ConstraintHandler's DofHandler must be the same") + end + return _add_interface_entries!(sp, dh, ch, topology, keep_constrained, interface_coupling) +end -Create a sparsity pattern accounting for affine constraints in `ch`. See -the Affine Constraints section of the manual for further details. """ -function create_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler; - keep_constrained::Bool=true, coupling=nothing, topology::Union{Nothing, AbstractTopology} = nothing, - cross_coupling = nothing) - return _create_sparsity_pattern(dh, ch, false, keep_constrained, coupling, topology, cross_coupling) + add_constraint_entries!( + sp::AbstractSparsityPattern, ch::ConstraintHandler; + keep_constrained::Bool = true, + ) + +Add all entries resulting from constraints in the ConstraintHandler `ch` to the sparsity +pattern. Note that, since this operation depends on existing entries in the pattern, this +function must be called as the *last* step when creating the sparsity pattern. + +# Keyword arguments + - `keep_constrained`: whether or not entries for constrained DoFs should be kept + (`keep_constrained = true`) or eliminated (`keep_constrained = false`) from the sparsity + pattern. +""" +function add_constraint_entries!( + sp::AbstractSparsityPattern, ch::ConstraintHandler; + keep_constrained::Bool = true, +) + return _add_constraint_entries!(sp, ch.dofcoefficients, ch.dofmapping, keep_constrained) end +function add_diagonal_entries!(sp::AbstractSparsityPattern) + for d in 1:min(getnrows(sp), getncols(sp)) + add_entry!(sp, d, d) + end + return sp +end + + +############################################################ +# Sparse matrix instantiation from AbstractSparsityPattern # +############################################################ + +""" + allocate_matrix(::Type{SparseMatrixCSC{Tv, Ti}}, sp::SparsityPattern) + +Allocate a sparse matrix of type `SparseMatrixCSC{Tv, Ti}` from the sparsity pattern `sp`. +""" +function allocate_matrix(::Type{S}, sp::AbstractSparsityPattern) where {Tv, Ti, S <: SparseMatrixCSC{Tv, Ti}} + return _allocate_matrix(S, sp, #=sym=# false) +end + +""" + 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. +""" +function allocate_matrix(::Type{Symmetric{Tv, S}}, sp::AbstractSparsityPattern) where {Tv, Ti, S <: SparseMatrixCSC{Tv, Ti}} + return Symmetric(_allocate_matrix(S, sp, #=sym=# true)) +end + +""" + 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}}). +""" +allocate_matrix(sp::SparsityPattern) = allocate_matrix(SparseMatrixCSC{Float64, Int}, sp) + +""" + 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`](@ref allocate_matrix(::Type{<:Any}, ::SparsityPattern)) for +supported matrix types, and to [`create_sparsity_pattern`](@ref) 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 + ```julia + sp = init_sparsity_pattern(dh) + add_sparsity_entries!(sp, dh) + K = allocate_matrix(sp) + M = allocate_matrix(sp) + ``` + instead of + ```julia + 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. +""" +function allocate_matrix(::Type{MatrixType}, dh::DofHandler, args...; kwargs...) where {MatrixType} + sp = init_sparsity_pattern(dh) + add_sparsity_entries!(sp, dh, args...; kwargs...) + return allocate_matrix(MatrixType, sp) +end + +""" + 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...)`](@ref allocate_matrix(::Type{MatrixType}, ::DofHandler, args...; +kwargs...) where {MatrixType}) -- refer to that method for details. +""" +function allocate_matrix(dh::DofHandler, args...; kwargs...) + return allocate_matrix(SparseMatrixCSC{Float64, Int}, dh, args...; kwargs...) +end + + +############################## +# Sparsity pattern internals # +############################## + # Compute a coupling matrix of size (ndofs_per_cell × ndofs_per_cell) based on the input # coupling which can be of size i) (nfields × nfields) specifying coupling between fields, # ii) (ncomponents × ncomponents) specifying coupling between components, or iii) # (ndofs_per_cell × ndofs_per_cell) specifying coupling between all local dofs, i.e. a # "template" local matrix. -function _coupling_to_local_dof_coupling(dh::DofHandler, coupling::AbstractMatrix{Bool}, sym::Bool) +function _coupling_to_local_dof_coupling(dh::DofHandler, coupling::AbstractMatrix{Bool}) sz = size(coupling, 1) sz == size(coupling, 2) || error("coupling not square") - sym && (issymmetric(coupling) || error("coupling not symmetric")) # Return one matrix per (potential) sub-domain outs = Matrix{Bool}[] @@ -102,38 +479,136 @@ function _coupling_to_local_dof_coupling(dh::DofHandler, coupling::AbstractMatri return outs end -""" - _add_cross_coupling(coupling_sdh, dof_i, dof_j, cell_field_dofs, neighbor_field_dofs, i, j, sym, keep_constrained, ch, cnt, I, J) +function _add_cell_entries!( + sp::AbstractSparsityPattern, dh::DofHandler, ch::Union{ConstraintHandler, Nothing}, + keep_constrained::Bool, coupling::Union{Vector{<:AbstractMatrix{Bool}}, Nothing}, + ) + # Add all connections between dofs for every cell while filtering based + # on a) constraints, and b) field/dof coupling. + cc = CellCache(dh) + for (sdhi, sdh) in pairs(dh.subdofhandlers) + set = BitSet(sdh.cellset) + coupling === nothing || (coupling_sdh = coupling[sdhi]) + for cell_id in set + reinit!(cc, cell_id) + for (i, row) in pairs(cc.dofs) + # a) check constraint for row + !keep_constrained && haskey(ch.dofmapping, row) && continue + # TODO: Extracting the row here and reinserting after the j-loop + # should give some nice speedup + for (j, col) in pairs(cc.dofs) + # b) check coupling between (local) dofs i and j + coupling === nothing || coupling_sdh[i, j] || continue + # a) check constraint for col + !keep_constrained && haskey(ch.dofmapping, col) && continue + # Insert col as a non zero index for this row + add_entry!(sp, row, col) + end + end + end + end + return sp +end -Helper function used to mutate `I` and `J` to add cross-element coupling. -""" -function _add_cross_coupling(coupling_sdh::Matrix{Bool}, dof_i::Int, dof_j::Int, +function _add_constraint_entries!( + sp::AbstractSparsityPattern, dofcoefficients::Vector{Union{DofCoefficients{T}, Nothing}}, + dofmapping::Dict{Int,Int}, keep_constrained::Bool, + ) where {T} + + # Return early if there are no non-trivial affine constraints + any(i -> !(i === nothing || isempty(i)), dofcoefficients) || return + + # New entries tracked separately and inserted after since it is not possible to modify + # the datastructure while looping over it. + mempool = PoolAllocator.MemoryPool{Int}() + sp′ = Dict{Int, PoolAllocator.PoolVector{Int}}() + + for (row, colidxs) in zip(1:getnrows(sp), eachrow(sp)) # pairs(eachrow(sp)) + row_coeffs = coefficients_for_dof(dofmapping, dofcoefficients, row) + if row_coeffs === nothing + # This row is _not_ constrained, check columns of this row... + !keep_constrained && haskey(dofmapping, row) && continue + for col in colidxs + col_coeffs = coefficients_for_dof(dofmapping, dofcoefficients, col) + if col_coeffs === nothing + # ... this column is _not_ constrained, done. + continue + else + # ... this column _is_ constrained, distribute to columns. + for (col′, _) in col_coeffs + r = get(sp′, row) do + PoolAllocator.resize(PoolAllocator.malloc(mempool, 8), 0) + end + r = insert_sorted(r, col′) + sp′[row] = r + end + end + end + else + # This row _is_ constrained, check columns of this row... + for col in colidxs + col_coeffs = coefficients_for_dof(dofmapping, dofcoefficients, col) + if col_coeffs === nothing + # ... this column is _not_ constrained, distribute to rows. + !keep_constrained && haskey(dofmapping, col) && continue + for (row′, _) in row_coeffs + r = get(sp′, row′) do + PoolAllocator.resize(PoolAllocator.malloc(mempool, 8), 0) + end + r = insert_sorted(r, col) + sp′[row′] = r + end + else + # ... this column _is_ constrained, double-distribute to columns/rows. + for (row′, _) in row_coeffs + !keep_constrained && haskey(dofmapping, row′) && continue + for (col′, _) in col_coeffs + !keep_constrained && haskey(dofmapping, col′) && continue + r = get(sp′, row′) do + PoolAllocator.resize(PoolAllocator.malloc(mempool, 8), 0) + end + r = insert_sorted(r, col′) + sp′[row′] = r + end + end + end + end + end + end + + # Insert new entries into the sparsity pattern + for (row, colidxs) in sp′ + # TODO: Extract row here and just insert_sorted + for col in colidxs + add_entry!(sp, row, col) + end + end + + return sp +end + +function _add_interface_entry(sp::SparsityPattern, coupling_sdh::Matrix{Bool}, dof_i::Int, dof_j::Int, cell_field_dofs::Union{Vector{Int}, SubArray}, neighbor_field_dofs::Union{Vector{Int}, SubArray}, - i::Int, j::Int, sym::Bool, keep_constrained::Bool, ch::Union{ConstraintHandler, Nothing}, cnt::Int, I::Vector{Int}, J::Vector{Int}) + i::Int, j::Int, keep_constrained::Bool, ch::Union{ConstraintHandler, Nothing}) - coupling_sdh[dof_i, dof_j] || return cnt + coupling_sdh[dof_i, dof_j] || return dofi = cell_field_dofs[i] dofj = neighbor_field_dofs[j] - sym && (dofj > dofi && return cnt) - !keep_constrained && (haskey(ch.dofmapping, dofi) || haskey(ch.dofmapping, dofj)) && return cnt - cnt += 1 - _add_or_grow(cnt, I, J, dofi, dofj) - return cnt + # sym && (dofj > dofi && return cnt) + !keep_constrained && (haskey(ch.dofmapping, dofi) || haskey(ch.dofmapping, dofj)) && return + add_entry!(sp, dofi, dofj) + return end -""" - cross_element_coupling!(dh::DofHandler, topology::ExclusiveTopology, sym::Bool, keep_constrained::Bool, couplings::Union{AbstractVector{<:AbstractMatrix{Bool}},Nothing}, cnt::Int, I::Vector{Int}, J::Vector{Int}) - -Mutates `I, J` to account for cross-element coupling by calling [`_add_cross_coupling`](@ref). -Returns the updated value of `cnt`. - -Used internally for sparsity patterns with cross-element coupling. -""" -function cross_element_coupling!(dh::DofHandler, ch::Union{ConstraintHandler, Nothing}, topology::ExclusiveTopology, sym::Bool, keep_constrained::Bool, couplings::AbstractVector{<:AbstractMatrix{Bool}}, cnt::Int, I::Vector{Int}, J::Vector{Int}) - fca = FacetCache(CellCache(dh, UpdateFlags(false, false, true)), Int[], ScalarWrapper(0)) - fcb = FacetCache(CellCache(dh, UpdateFlags(false, false, true)), Int[], ScalarWrapper(0)) - ic = InterfaceCache(fca, fcb, Int[]) - for ic in InterfaceIterator(ic, dh.grid, topology) +function _add_interface_entries!( + sp::SparsityPattern, dh::DofHandler, ch::Union{ConstraintHandler, Nothing}, + topology::ExclusiveTopology, keep_constrained::Bool, + interface_coupling::AbstractMatrix{Bool}, + ) + couplings = _coupling_to_local_dof_coupling(dh, interface_coupling) + for ic in InterfaceIterator(dh, topology) + # TODO: This looks like it can be optimized for the common case where + # the cells are in the same subdofhandler sdhs_idx = dh.cell_to_subdofhandler[cellid.([ic.a, ic.b])] sdhs = dh.subdofhandlers[sdhs_idx] for (i, sdh) in pairs(sdhs) @@ -144,7 +619,7 @@ function cross_element_coupling!(dh::DofHandler, ch::Union{ConstraintHandler, No cell_dofs = celldofs(sdh_idx == 1 ? ic.a : ic.b) cell_field_dofs = @view cell_dofs[dofrange1] for neighbor_field in sdh.field_names - sdh2 = sdhs[i==1 ? 2 : 1] + sdh2 = sdhs[i == 1 ? 2 : 1] neighbor_field ∈ sdh2.field_names || continue dofrange2 = dof_range(sdh2, neighbor_field) neighbor_dofs = celldofs(sdh_idx == 2 ? ic.a : ic.b) @@ -153,152 +628,44 @@ function cross_element_coupling!(dh::DofHandler, ch::Union{ConstraintHandler, No for (j, dof_j) in pairs(dofrange2), (i, dof_i) in pairs(dofrange1) # This line to avoid coupling the shared dof in continuous interpolations as cross-element. They're coupled in the local coupling matrix. (cell_field_dofs[i] ∈ neighbor_dofs || neighbor_field_dofs[j] ∈ cell_dofs) && continue - cnt = _add_cross_coupling(coupling_sdh, dof_i, dof_j, cell_field_dofs, neighbor_field_dofs, i, j, sym, keep_constrained, ch, cnt, I, J) - cnt = _add_cross_coupling(coupling_sdh, dof_j, dof_i, neighbor_field_dofs, cell_field_dofs, j, i, sym, keep_constrained, ch, cnt, I, J) + _add_interface_entry(sp, coupling_sdh, dof_i, dof_j, cell_field_dofs, neighbor_field_dofs, i, j, keep_constrained, ch) + _add_interface_entry(sp, coupling_sdh, dof_j, dof_i, neighbor_field_dofs, cell_field_dofs, j, i, keep_constrained, ch) end end end end end - return cnt + return sp end -function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{ConstraintHandler, Nothing}=#, sym::Bool, keep_constrained::Bool, coupling::Union{AbstractMatrix{Bool},Nothing}, - topology::Union{Nothing, AbstractTopology}, cross_coupling::Union{AbstractMatrix{Bool},Nothing}) - @assert isclosed(dh) - if !keep_constrained - @assert ch !== nothing && isclosed(ch) - end - - couplings = isnothing(coupling) ? nothing : _coupling_to_local_dof_coupling(dh, coupling, sym) - cross_couplings = isnothing(cross_coupling) ? nothing : _coupling_to_local_dof_coupling(dh, cross_coupling, sym) - - # Allocate buffers. Compute an upper bound for the buffer length and allocate it all up - # front since they will become large and expensive to re-allocate. The bound is exact - # when keeping constrained dofs (default) and if not it only over-estimates with number - # of entries eliminated by constraints. - max_buffer_length = ndofs(dh) # diagonal elements - for (sdh_idx, sdh) in pairs(dh.subdofhandlers) - set = sdh.cellset - n = ndofs_per_cell(sdh) - entries_per_cell = if coupling === nothing - sym ? div(n * (n + 1), 2) : n^2 - else - coupling_sdh = couplings[sdh_idx] - count(coupling_sdh[i, j] for i in 1:n for j in (sym ? i : 1):n) +# Internal matrix instantiation for SparseMatrixCSC and Symmetric{SparseMatrixCSC} +function _allocate_matrix(::Type{SparseMatrixCSC{Tv, Ti}}, sp::AbstractSparsityPattern, sym::Bool) where {Tv, Ti} + # 1. Setup colptr + colptr = zeros(Ti, getncols(sp) + 1) + colptr[1] = 1 + for (row, colidxs) in enumerate(eachrow(sp)) + for col in colidxs + sym && row > col && continue + colptr[col+1] += 1 end - max_buffer_length += entries_per_cell * length(set) end - I = Vector{Int}(undef, max_buffer_length) - J = Vector{Int}(undef, max_buffer_length) - global_dofs = Int[] - cnt = 0 - - for (sdh_idx, sdh) in pairs(dh.subdofhandlers) - coupling === nothing || (coupling_sdh = couplings[sdh_idx]) - # TODO: Remove BitSet construction when SubDofHandler ensures sorted collections - set = BitSet(sdh.cellset) - n = ndofs_per_cell(sdh) - resize!(global_dofs, n) - @inbounds for element_id in set - celldofs!(global_dofs, dh, element_id) - for j in eachindex(global_dofs), i in eachindex(global_dofs) - coupling === nothing || coupling_sdh[i, j] || continue - dofi = global_dofs[i] - dofj = global_dofs[j] - sym && (dofi > dofj && continue) - !keep_constrained && (haskey(ch.dofmapping, dofi) || haskey(ch.dofmapping, dofj)) && continue - cnt += 1 - I[cnt] = dofi - J[cnt] = dofj - end + cumsum!(colptr, colptr) + nnz = colptr[end] - 1 + # 2. Allocate rowval and nzval now that nnz is known + rowval = Vector{Ti}(undef, nnz) + nzval = zeros(Tv, nnz) + # 3. Populate rowval. Since SparsityPattern is row-based we need to allocate an extra + # work buffer here to keep track of the next index into rowval + nextinds = copy(colptr) + for (row, colidxs) in zip(1:getnrows(sp), eachrow(sp)) # pairs(eachrow(sp)) + for col in colidxs + sym && row > col && continue + k = nextinds[col] + rowval[k] = row + nextinds[col] = k + 1 end end - if !isnothing(topology) && !isnothing(cross_coupling) && any(cross_coupling) - cnt = cross_element_coupling!(dh, ch, topology, sym, keep_constrained, cross_couplings, cnt, I, J) - end - # Always add diagonal entries - resize!(I, cnt + ndofs(dh)) - resize!(J, cnt + ndofs(dh)) - @inbounds for d in 1:ndofs(dh) - cnt += 1 - I[cnt] = d - J[cnt] = d - end - @assert length(I) == length(J) == cnt - - K = spzeros!!(Float64, I, J, ndofs(dh), ndofs(dh)) - - # If ConstraintHandler is given, create the condensation pattern due to affine constraints - if ch !== nothing - @assert isclosed(ch) - fill!(K.nzval, 1) - _condense_sparsity_pattern!(K, ch.dofcoefficients, ch.dofmapping, keep_constrained) - fillzero!(K) - end - - return K -end - -# Similar to Ferrite._condense!(K, ch), but only add the non-zero entries to K (that arises from the condensation process) -function _condense_sparsity_pattern!(K::SparseMatrixCSC{T}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, dofmapping::Dict{Int,Int}, keep_constrained::Bool) where T - ndofs = size(K, 1) - - # Return early if there are no non-trivial affine constraints - any(i -> !(i === nothing || isempty(i)), dofcoefficients) || return - - # Adding new entries to K is extremely slow, so create a new sparsity triplet for the - # condensed sparsity pattern - N = 2 * length(dofcoefficients) # TODO: Better size estimate for additional condensed sparsity pattern. - I = Int[]; resize!(I, N) - J = Int[]; resize!(J, N) - - cnt = 0 - for col in 1:ndofs - col_coeffs = coefficients_for_dof(dofmapping, dofcoefficients, col) - if col_coeffs === nothing - !keep_constrained && haskey(dofmapping, col) && continue - for ri in nzrange(K, col) - row = K.rowval[ri] - row_coeffs = coefficients_for_dof(dofmapping, dofcoefficients, row) - row_coeffs === nothing && continue - for (d, _) in row_coeffs - cnt += 1 - _add_or_grow(cnt, I, J, d, col) - end - end - else - for ri in nzrange(K, col) - row = K.rowval[ri] - row_coeffs = coefficients_for_dof(dofmapping, dofcoefficients, row) - if row_coeffs === nothing - !keep_constrained && haskey(dofmapping, row) && continue - for (d, _) in col_coeffs - cnt += 1 - _add_or_grow(cnt, I, J, row, d) - end - else - for (d1, _) in col_coeffs - !keep_constrained && haskey(dofmapping, d1) && continue - for (d2, _) in row_coeffs - !keep_constrained && haskey(dofmapping, d2) && continue - cnt += 1 - _add_or_grow(cnt, I, J, d1, d2) - end - end - end - end - end - end - - resize!(I, cnt) - resize!(J, cnt) - - # Fill the sparse matrix with a non-zero value so that :+ operation does not remove entries with value zero. - K2 = spzeros!!(Float64, I, J, ndofs, ndofs) - fill!(K2.nzval, 1) - - K .+= K2 - - return nothing + @assert all(i -> nextinds[i] == colptr[i + 1], 1:getncols(sp)) + S = SparseMatrixCSC(getnrows(sp), getncols(sp), colptr, rowval, nzval) + return S end diff --git a/src/Ferrite.jl b/src/Ferrite.jl index 41be4a4534..ca5ef421e5 100644 --- a/src/Ferrite.jl +++ b/src/Ferrite.jl @@ -8,8 +8,7 @@ using Base: using EnumX: EnumX, @enumx using LinearAlgebra: - LinearAlgebra, Symmetric, Transpose, cholesky, det, issymmetric, norm, - pinv, tr + LinearAlgebra, Symmetric, Transpose, cholesky, det, norm, pinv, tr using NearestNeighbors: NearestNeighbors, KDTree, knn using OrderedCollections: @@ -107,6 +106,7 @@ const AbstractVecOrSet{T} = Union{AbstractSet{T}, AbstractVector{T}} const IntegerCollection = AbstractVecOrSet{<:Integer} include("utils.jl") +include("PoolAllocator.jl") # Matrix/Vector utilities include("arrayutils.jl") @@ -140,6 +140,7 @@ include("Dofs/DofHandler.jl") include("Dofs/ConstraintHandler.jl") include("Dofs/apply_analytical.jl") include("Dofs/sparsity_pattern.jl") +include("Dofs/block_sparsity_pattern.jl") include("Dofs/DofRenumbering.jl") include("iterators.jl") diff --git a/src/L2_projection.jl b/src/L2_projection.jl index ec0a15a0fb..34da061225 100644 --- a/src/L2_projection.jl +++ b/src/L2_projection.jl @@ -130,7 +130,7 @@ of the mass-matrix. function close!(proj::L2Projector) close!(proj.dh) M = _assemble_L2_matrix(proj.dh, proj.qrs_lhs) - proj.M_cholesky = cholesky(M) + proj.M_cholesky = cholesky(Symmetric(M)) return proj end @@ -144,7 +144,7 @@ end _mass_qr(ip::VectorizedInterpolation) = _mass_qr(ip.ip) function _assemble_L2_matrix(dh::DofHandler, qrs_lhs::Vector{<:QuadratureRule}) - M = create_symmetric_sparsity_pattern(dh) + M = Symmetric(allocate_matrix(dh)) assembler = start_assemble(M) for (sdh, qr_lhs) in zip(dh.subdofhandlers, qrs_lhs) ip_fun = only(sdh.field_interpolations) @@ -160,7 +160,7 @@ function _assemble_L2_matrix!(assembler, cellvalues::CellValues, sdh::SubDofHand n = getnbasefunctions(cellvalues) Me = zeros(n, n) - function symmetrize_to_lower!(K) + function symmetrize_to_lower!(K::Matrix) for i in 1:size(K, 1) for j in i+1:size(K, 1) K[j, i] = K[i, j] @@ -286,7 +286,7 @@ function _project(proj::L2Projector, qrs_rhs::Vector{<:QuadratureRule}, vars::Un # Recast to original input type make_T(vals) = T <: AbstractTensor ? T(Tuple(vals)) : vals[1] - return T[make_T(x) for x in eachrow(projected_vals)] + return T[make_T(x) for x in Base.eachrow(projected_vals)] end function assemble_proj_rhs!(f::Matrix, cellvalues::CellValues, sdh::SubDofHandler, vars::Union{AbstractVector, AbstractDict}) diff --git a/src/PoolAllocator.jl b/src/PoolAllocator.jl new file mode 100644 index 0000000000..2ac84ceb71 --- /dev/null +++ b/src/PoolAllocator.jl @@ -0,0 +1,237 @@ +module PoolAllocator + +# Checkmate LanguageServer.jl +const var"@propagate_inbounds" = Base.var"@propagate_inbounds" + +@eval macro $(Symbol("const"))(field) + if VERSION >= v"1.8.0-DEV.1148" + Expr(:const, esc(field)) + else + return esc(field) + end +end + +const PAGE_SIZE = 4 * 1024 * 1024 # 4 MiB + +# A page corresponds to a memory block of size `PAGE_SIZE` bytes. +# Allocations of arrays are views into this block. +mutable struct Page{T} + @const buf::Vector{T} # data buffer (TODO: Memory in recent Julias?) + @const blocksize::Int # blocksize for this page + @const freelist::BitVector # block is free/used + n_free::Int # number of free blocks + function Page{T}(blocksize::Int) where T + @assert isbitstype(T) + buf = Vector{T}(undef, PAGE_SIZE ÷ sizeof(T)) + n_blocks, r = divrem(length(buf), blocksize) + @assert r == 0 + return new{T}(buf, blocksize, trues(n_blocks), n_blocks) + end +end + +# Find a free block and mark it as used +function _malloc(page::Page, size::Int) + if size != page.blocksize + error("malloc: requested size does not match the blocksize of this page") + end + # Return early if the page is already full + if page.n_free == 0 + return nothing + end + # Find the first free block + blockindex = findfirst(page.freelist)::Int + if !@inbounds(page.freelist[blockindex]) + error("malloc: block already in use") + end + @inbounds page.freelist[blockindex] = false + offset = (blockindex - 1) * page.blocksize + page.n_free -= 1 + return offset +end + +# Mark a block as free +function _free(page::Page, offset::Int) + blockindex = offset ÷ page.blocksize + 1 + if @inbounds page.freelist[blockindex] + error("free: block already free'd") + end + @inbounds page.freelist[blockindex] = true + page.n_free += 1 + # TODO: If this page is completely unused it can be collected and reused. + return +end + +# A book is a collection of pages with a specific blocksize +struct Book{T} + blocksize::Int + pages::Vector{Page{T}} +end + +# Find a page with a free block of the requested size +function _malloc(book::Book{T}, size::Int) where {T} + @assert book.blocksize == size + # Check existing pages + for page in book.pages + offset = _malloc(page, size) + if offset !== nothing + return (page, offset) + end + end + # Allocate a new page + page = Page{T}(book.blocksize) + push!(book.pages, page) + # Allocate block in the new page + offset = _malloc(page, size) + @assert offset !== nothing + return (page, offset) +end + +struct MemoryPool{T} + books::Vector{Book{T}} # blocksizes 2, 4, 6, 8, ... + function MemoryPool{T}() where T + mempool = new(Book{T}[]) + return mempool + end +end + +# Free all pages by resizing all page containers to 0 +function free(mempool::MemoryPool) + for i in 1:length(mempool.books) + isassigned(mempool.books, i) || continue + resize!(mempool.books[i].pages, 0) + end + resize!(mempool.books, 0) + return +end + +function mempool_stats(mempool::MemoryPool{T}) where T + bytes_used = 0 + bytes_allocated = 0 + for bookidx in 1:length(mempool.books) + isassigned(mempool.books, bookidx) || continue + book = mempool.books[bookidx] + bytes_allocated += length(book.pages) * PAGE_SIZE + for page in book.pages + bytes_used += count(!, page.freelist) * page.blocksize * sizeof(T) + end + end + return bytes_used, bytes_allocated +end + +function Base.show(io::IO, ::MIME"text/plain", mempool::MemoryPool{T}) where T + n_books = count(i -> isassigned(mempool.books, i), 1:length(mempool.books)) + print(io, "PoolAllocator.MemoryPool{$(T)} with $(n_books) fixed size pools") + n_books == 0 && return + println(io, ":") + for idx in 1:length(mempool.books) + isassigned(mempool.books, idx) || continue + h = mempool.books[idx] + blocksize = h.blocksize + # @assert blocksize == 2^idx + npages = length(h.pages) + n_free = mapreduce(p -> p.n_free, +, h.pages; init=0) + n_tot = npages * PAGE_SIZE ÷ blocksize ÷ sizeof(T) + println(io, " - blocksize: $(blocksize), npages: $(npages), usage: $(n_tot - n_free) / $(n_tot)") + end + return +end + +function bookindex_from_blocksize(blocksize::Int) + return (8 * sizeof(Int) - leading_zeros(blocksize)) % Int +end + +function malloc(mempool::MemoryPool{T}, dims::NTuple{N, Int}) where {T, N} + @assert prod(dims) > 0 + blocksize = nextpow(2, prod(dims)) + bookidx = bookindex_from_blocksize(blocksize) + if length(mempool.books) < bookidx + resize!(mempool.books, bookidx) + end + if !isassigned(mempool.books, bookidx) + book = Book(blocksize, Page{T}[]) + mempool.books[bookidx] = book + else + book = mempool.books[bookidx] + end + page, offset = _malloc(book, blocksize) + + return PoolArray{T, N}(mempool, page, offset, dims) +end + + +# PoolArray is a view into a page that also has a reference to the MemoryPool so that it can +# be resized/reallocated. +struct PoolArray{T, N} <: AbstractArray{T, N} + mempool::MemoryPool{T} + page::Page{T} + offset::Int + size::NTuple{N, Int} +end + +const PoolVector{T} = PoolArray{T, 1} + +# Constructors +function malloc(mempool::MemoryPool, dim1::Int) + return malloc(mempool, (dim1, )) +end +function malloc(mempool::MemoryPool, dim1::Int, dim2::Int, dimx::Int...) + dims = (dim1, dim2, map(Int, dimx)..., ) + return malloc(mempool, dims) +end + +function free(x::PoolArray) + _free(x.page, x.offset) + return +end + +function realloc(x::PoolArray{T}, newsize::Int) where T + @assert newsize > length(x) # TODO: Allow shrinkage? + @assert newsize <= PAGE_SIZE ÷ sizeof(T) # TODO: Might be required + # Find the page for the block to make sure it was allocated in this mempool + # page = find_page(x.mempool, ) + # Allocate the new block + x′ = malloc(x.mempool, newsize) + # Copy the data + copyto!(x′, x) + # Free the old block and return + _free(x.page, x.offset) + return x′ +end + +# AbstractArray interface +Base.size(mv::PoolArray) = mv.size +allocated_length(mv::PoolArray) = mv.page.blocksize +Base.IndexStyle(::Type{<:PoolArray}) = IndexLinear() +@propagate_inbounds function Base.getindex(mv::PoolArray, i::Int) + @boundscheck checkbounds(mv, i) + return @inbounds mv.page.buf[mv.offset + i] +end +@propagate_inbounds function Base.setindex!(mv::PoolArray{T}, v::T, i::Int) where T + @boundscheck checkbounds(mv, i) + @inbounds mv.page.buf[mv.offset + i] = v + return mv +end + +# Utilities needed for the sparsity pattern +@inline function resize(x::PoolVector{T}, n::Int) where T + if n > allocated_length(x) + return realloc(x, n) + else + return PoolVector{T}(x.mempool, x.page, x.offset, (n, )) + end +end + +@inline function insert(x::PoolVector{T}, k::Int, item::T) where T + lx = length(x) + # Make room + x = resize(x, lx + 1) + # Shift elements after the insertion point to the back + @inbounds for i in lx:-1:k + x[i + 1] = x[i] + end + # Insert the new element + @inbounds x[k] = item + return x +end + +end # module PoolAllocator diff --git a/src/deprecations.jl b/src/deprecations.jl index a62d61fef8..bbd9d8f682 100644 --- a/src/deprecations.jl +++ b/src/deprecations.jl @@ -377,3 +377,5 @@ function default_interpolation(::Type{C}) where {C <: AbstractCell} @warn "Ferrite.default_interpolation is deprecated, use the exported `geometric_interpolation` instead" maxlog=1 return geometric_interpolation(C) end + +@deprecate create_sparsity_pattern allocate_matrix diff --git a/src/exports.jl b/src/exports.jl index a8fe8701f2..3469d8a5cf 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -116,14 +116,23 @@ export ndofs_per_cell, celldofs!, celldofs, - create_sparsity_pattern, - create_symmetric_sparsity_pattern, dof_range, renumber!, DofOrder, evaluate_at_grid_nodes, apply_analytical!, +# Sparsity pattern + # AbstractSparsityPattern, + SparsityPattern, + BlockSparsityPattern, + init_sparsity_pattern, + add_sparsity_entries!, + add_cell_entries!, + add_interface_entries!, + add_constraint_entries!, + allocate_matrix, + # Constraints ConstraintHandler, Dirichlet, diff --git a/test/PoolAllocator.jl b/test/PoolAllocator.jl new file mode 100644 index 0000000000..4c50003cb7 --- /dev/null +++ b/test/PoolAllocator.jl @@ -0,0 +1,81 @@ +using Test, Ferrite.PoolAllocator + +@testset "PoolAllocator.jl" begin + + # Basic malloc, realloc, free + mempool = PoolAllocator.MemoryPool{Int}() + x = PoolAllocator.malloc(mempool, 1024) + @test x isa PoolAllocator.PoolArray{Int} + x .= 1:1024 + x′ = PoolAllocator.realloc(x, 2048) + @test x′ isa PoolAllocator.PoolArray{Int} + @test_throws ErrorException("free: block already free'd") PoolAllocator.free(x) + @test x′[1:1024] == 1:1024 + PoolAllocator.free(x′) + @test_throws ErrorException("free: block already free'd") PoolAllocator.free(x′) + + # Internal page allocation: exhaust some pages + mempool = PoolAllocator.MemoryPool{Int}() + xs = PoolAllocator.PoolArray{Int}[] + for _ in 1:(PoolAllocator.PAGE_SIZE ÷ 512 ÷ sizeof(Int) * 2 + 1) + x = PoolAllocator.malloc(mempool, 512) + push!(xs, x) + end + @test length(mempool.books[10].pages) == 3 + @test all(!, mempool.books[10].pages[1].freelist) + @test all(!, mempool.books[10].pages[2].freelist) + @test !mempool.books[10].pages[3].freelist[1] + @test all(mempool.books[10].pages[3].freelist[2:end]) + xs′ = PoolAllocator.PoolArray{Int}[] + for x in xs + x′ = PoolAllocator.realloc(x, 1024) + @test_throws ErrorException("free: block already free'd") PoolAllocator.free(x) + push!(xs′, x′) + end + @test length(mempool.books[10].pages) == 3 # TODO + @test all(mempool.books[10].pages[1].freelist) + @test all(mempool.books[10].pages[2].freelist) + @test all(mempool.books[10].pages[3].freelist) + @test length(mempool.books[11].pages) == 5 + @test all(!, mempool.books[11].pages[1].freelist) + @test all(!, mempool.books[11].pages[2].freelist) + @test all(!, mempool.books[11].pages[3].freelist) + @test all(!, mempool.books[11].pages[4].freelist) + @test !mempool.books[11].pages[5].freelist[1] + @test all(mempool.books[11].pages[5].freelist[2:end]) + for x in xs′ + PoolAllocator.free(x) + @test_throws ErrorException("free: block already free'd") PoolAllocator.free(x) + end + PoolAllocator.free(mempool) + @test length(mempool.books) == 0 + + # Array functions + mempool = PoolAllocator.MemoryPool{Int}() + x = PoolAllocator.malloc(mempool, 8) + @test length(x) == 8 + x = PoolAllocator.resize(x, 0) + @test length(x) == 0 + x = PoolAllocator.resize(x, 16) + @test length(x) == 16 + x .= 1:16 + x = PoolAllocator.resize(x, 8) + @test x == 1:8 + x = PoolAllocator.resize(x, 8) + x = PoolAllocator.insert(x, 1, -1) + x = PoolAllocator.insert(x, length(x) + 1, -1) + x = PoolAllocator.insert(x, 2, -2) + x = PoolAllocator.insert(x, length(x), -2) + @test x == [-1; -2; 1:8; -2; -1] + + # n-d arrays + mempool = PoolAllocator.MemoryPool{Int}() + A = PoolAllocator.malloc(mempool, 64, 64) + @test size(A) == (64, 64) + B = PoolAllocator.malloc(mempool, (64, 32)) + @test size(B) == (64, 32) + + # Smoke test for `show` + show(devnull, MIME"text/plain"(), mempool) + +end diff --git a/test/blockarrays.jl b/test/blockarrays.jl index 59b1d8db9a..81b1f627c8 100644 --- a/test/blockarrays.jl +++ b/test/blockarrays.jl @@ -19,9 +19,12 @@ using Ferrite, BlockArrays, SparseArrays, Test close!(ch) update!(ch, 0) - K = create_sparsity_pattern(dh, ch) + K = allocate_matrix(dh, ch) f = zeros(axes(K, 1)) - KB = create_sparsity_pattern(BlockMatrix, dh, ch) + # TODO: allocate_matrix(BlockMatrix, ...) should work and default to field blocking + bsp = BlockSparsityPattern([2nd, 1nd]) + add_sparsity_entries!(bsp, dh, ch) + KB = allocate_matrix(BlockMatrix, bsp) @test KB isa BlockMatrix @test blocksize(KB) == (2, 2) @test size(KB[Block(1), Block(1)]) == (2nd, 2nd) @@ -66,16 +69,16 @@ using Ferrite, BlockArrays, SparseArrays, Test # Global application of BC not supported yet @test_throws ErrorException apply!(KB, fB, ch) - # Custom blocking by passing a partially initialized matrix + # Custom blocking perm = invperm([ch.free_dofs; ch.prescribed_dofs]) renumber!(dh, ch, perm) nfree = length(ch.free_dofs) npres = length(ch.prescribed_dofs) - K = create_sparsity_pattern(dh, ch) + K = allocate_matrix(dh, ch) block_sizes = [nfree, npres] - KBtmp = BlockArray(undef_blocks, SparseMatrixCSC{Float64, Int}, block_sizes, block_sizes) - KB = create_sparsity_pattern(KBtmp, dh, ch) - @test KBtmp === KB + bsp = BlockSparsityPattern(block_sizes) + add_sparsity_entries!(bsp, dh, ch) + KB = allocate_matrix(BlockMatrix, bsp) @test blocksize(KB) == (2, 2) @test size(KB[Block(1), Block(1)]) == (nfree, nfree) @test size(KB[Block(2), Block(1)]) == (npres, nfree) diff --git a/test/integration/test_simple_scalar_convergence.jl b/test/integration/test_simple_scalar_convergence.jl index 75052af5de..238efe9fac 100644 --- a/test/integration/test_simple_scalar_convergence.jl +++ b/test/integration/test_simple_scalar_convergence.jl @@ -112,7 +112,7 @@ end # Assemble and solve function solve(dh, ch, cellvalues) - K, f = assemble_global(cellvalues, create_sparsity_pattern(dh), dh); + K, f = assemble_global(cellvalues, allocate_matrix(dh), dh); apply!(K, f, ch) u = K \ f; end diff --git a/test/runtests.jl b/test/runtests.jl index 201d99eb45..50e77155e2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -41,6 +41,7 @@ include("test_interfacevalues.jl") include("test_quadrules.jl") include("test_assemble.jl") include("test_dofs.jl") +include("test_sparsity_patterns.jl") include("test_constraints.jl") include("test_grid_dofhandler_vtk.jl") include("test_vtk_export.jl") @@ -53,6 +54,7 @@ include("test_pointevaluation.jl") # include("test_notebooks.jl") include("test_apply_rhs.jl") include("test_apply_analytical.jl") +include("PoolAllocator.jl") include("test_deprecations.jl") HAS_EXTENSIONS && include("blockarrays.jl") include("test_examples.jl") diff --git a/test/test_abstractgrid.jl b/test/test_abstractgrid.jl index f0896adc77..a069e175d7 100644 --- a/test/test_abstractgrid.jl +++ b/test/test_abstractgrid.jl @@ -69,7 +69,7 @@ add!(ch, dbc) close!(ch) update!(ch, 0.0) - K = create_sparsity_pattern(dh); + K = allocate_matrix(dh); K, f = doassemble!(cellvalues, K, dh); apply!(K, f, ch) sol = K \ f diff --git a/test/test_apply_rhs.jl b/test/test_apply_rhs.jl index bcd2517d22..d21a8fc6f9 100644 --- a/test/test_apply_rhs.jl +++ b/test/test_apply_rhs.jl @@ -8,7 +8,7 @@ function test_apply_rhs() add!(dh, :u, ip) close!(dh) - K = create_sparsity_pattern(dh) + K = allocate_matrix(dh) ch = ConstraintHandler(dh) @@ -62,7 +62,7 @@ function test_apply_rhs() end K, f = doassemble!(cellvalues, K, dh) - A = create_sparsity_pattern(dh) + A = allocate_matrix(dh) A, g = doassemble!(cellvalues, A, dh) rhsdata = get_rhs_data(ch, A) diff --git a/test/test_constraints.jl b/test/test_constraints.jl index 8beafec108..b660c20998 100644 --- a/test/test_constraints.jl +++ b/test/test_constraints.jl @@ -343,9 +343,9 @@ end C, g = Ferrite.create_constraint_matrix(ch) # Assemble - K = create_sparsity_pattern(dh, ch) + K = allocate_matrix(dh, ch) f = zeros(ndofs(dh)); f[end] = 1.0 - Kl = create_sparsity_pattern(dh, ch) + Kl = allocate_matrix(dh, ch) fl = copy(f) assembler = start_assemble(Kl, fl) for cell in CellIterator(dh) @@ -410,7 +410,7 @@ end add!(ch, AffineConstraint(1, [3=>params.a], params.b)) close!(ch) - K = create_sparsity_pattern(dh, ch) + K = allocate_matrix(dh, ch) r = zeros(ndofs(dh)) a = zeros(ndofs(dh)) @@ -1293,10 +1293,10 @@ end # testset close!(ch2) update!(ch2, 0) - K1 = create_sparsity_pattern(dh, ch1) + K1 = allocate_matrix(dh, ch1) f1 = zeros(ndofs(dh)) a1 = start_assemble(K1, f1) - K2 = create_sparsity_pattern(dh, ch2) + K2 = allocate_matrix(dh, ch2) f2 = zeros(ndofs(dh)) a2 = start_assemble(K2, f2) @@ -1333,9 +1333,9 @@ end # subtestset add!(ch2, Dirichlet(:u, getfacetset(grid, "right"), (x, t) -> 2.0t)) close!(ch2) - K1 = create_sparsity_pattern(dh, ch1) + K1 = allocate_matrix(dh, ch1) f1 = zeros(ndofs(dh)) - K2 = create_sparsity_pattern(dh, ch2) + K2 = allocate_matrix(dh, ch2) f2 = zeros(ndofs(dh)) for t in (1.0, 2.0) @@ -1415,7 +1415,7 @@ end # testset for azero in (nothing, false, true) - S = create_sparsity_pattern(dh) + S = allocate_matrix(dh) f = zeros(ndofs(dh)) K_dbc_standard = copy(S) @@ -1428,7 +1428,7 @@ end # testset f_dbc_local = copy(f) assembler_dbc_local = start_assemble(K_dbc_local, f_dbc_local) - S = create_sparsity_pattern(dh, ch_ac) + S = allocate_matrix(dh, ch_ac) K_ac_standard = copy(S) f_ac_standard = copy(f) @@ -1440,7 +1440,7 @@ end # testset f_ac_local = copy(f) assembler_ac_local = start_assemble(K_ac_local, f_ac_local) - S = create_sparsity_pattern(dh, ch_p) + S = allocate_matrix(dh, ch_p) K_p_standard = copy(S) f_p_standard = copy(f) @@ -1573,8 +1573,8 @@ end # testset ch = ConstraintHandler(dh) add!(ch, Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> 0)) close!(ch) - Kfull = create_sparsity_pattern(dh, ch) - K = create_sparsity_pattern(dh, ch; keep_constrained=false) + Kfull = allocate_matrix(dh, ch) + K = allocate_matrix(dh, ch; keep_constrained=false) # Pattern tests nonzero_edges = Set( (i, j) for d in 1:getncells(grid) diff --git a/test/test_dofs.jl b/test/test_dofs.jl index 573d087304..f0417f2d5f 100644 --- a/test/test_dofs.jl +++ b/test/test_dofs.jl @@ -189,12 +189,12 @@ end @test original_dofcoefficients == ch.dofcoefficients # Integration tests - K = create_sparsity_pattern(dh, ch) + K = allocate_matrix(dh, ch) f = zeros(ndofs(dh)) a = start_assemble(K, f) dhp, _, chp = dhmdhch() renumber!(dhp, chp, perm) - Kp = create_sparsity_pattern(dhp, chp) + Kp = allocate_matrix(dhp, chp) fp = zeros(ndofs(dhp)) ap = start_assemble(Kp, fp) for cellid in 1:getncells(dh.grid) @@ -438,11 +438,17 @@ end end return false end + function is_stored(sparsity_pattern::SparsityPattern, i, j) + return findfirst(k -> k == j, sparsity_pattern.rows[i]) !== nothing + end # Full coupling (default) - K = create_sparsity_pattern(dh) + sparsity_pattern = init_sparsity_pattern(dh) + add_sparsity_entries!(sparsity_pattern, dh) + K = allocate_matrix(sparsity_pattern) @test eltype(K) == Float64 for j in 1:ndofs(dh), i in 1:ndofs(dh) + @test is_stored(sparsity_pattern, i, j) @test is_stored(K, i, j) end @@ -452,25 +458,30 @@ end true true # v true false # q ] - K = create_sparsity_pattern(dh; coupling=coupling) - Kch = create_sparsity_pattern(dh, ch; coupling=coupling) - @test K.rowval == Kch.rowval - @test K.colptr == Kch.colptr - KS = create_symmetric_sparsity_pattern(dh; coupling=coupling) - KSch = create_symmetric_sparsity_pattern(dh, ch; coupling=coupling) - @test KS.data.rowval == KSch.data.rowval - @test KS.data.colptr == KSch.data.colptr + sparsity_pattern = init_sparsity_pattern(dh) + add_sparsity_entries!(sparsity_pattern, dh; coupling=coupling) + K = allocate_matrix(sparsity_pattern) + # Kch = allocate_matrix(dh, ch; coupling=coupling) + # @test K.rowval == Kch.rowval + # @test K.colptr == Kch.colptr + # KS = create_symmetric_sparsity_pattern(dh; coupling=coupling) + # KSch = create_symmetric_sparsity_pattern(dh, ch; coupling=coupling) + # @test KS.data.rowval == KSch.data.rowval + # @test KS.data.colptr == KSch.data.colptr for j in udofs, i in Iterators.flatten((vdofs, qdofs)) + @test is_stored(sparsity_pattern, i, j) @test is_stored(K, i, j) - @test is_stored(KS, i, j) == (i <= j) + # @test is_stored(KS, i, j) == (i <= j) end for j in pdofs, i in vdofs + @test is_stored(sparsity_pattern, i, j) @test is_stored(K, i, j) - @test is_stored(KS, i, j) + # @test is_stored(KS, i, j) end for j in pdofs, i in qdofs + @test is_stored(sparsity_pattern, i, j) == (i == j) @test is_stored(K, i, j) == (i == j) - @test is_stored(KS, i, j) == (i == j) + # @test is_stored(KS, i, j) == (i == j) end # Component coupling @@ -480,37 +491,45 @@ end true false true # v2 false true true # q ] - K = create_sparsity_pattern(dh; coupling=coupling) - KS = create_symmetric_sparsity_pattern(dh; coupling=coupling) + sparsity_pattern = init_sparsity_pattern(dh) + add_sparsity_entries!(sparsity_pattern, dh; coupling=coupling) + K = allocate_matrix(sparsity_pattern) + # KS = create_symmetric_sparsity_pattern(dh; coupling=coupling) for j in u1dofs, i in vdofs + @test is_stored(sparsity_pattern, i, j) @test is_stored(K, i, j) - @test is_stored(KS, i, j) == (i <= j) + # @test is_stored(KS, i, j) == (i <= j) end for j in u1dofs, i in qdofs + @test !is_stored(sparsity_pattern, i, j) @test !is_stored(K, i, j) - @test !is_stored(KS, i, j) + # @test !is_stored(KS, i, j) end for j in u2dofs, i in Iterators.flatten((v1dofs, qdofs)) + @test is_stored(sparsity_pattern, i, j) @test is_stored(K, i, j) - @test is_stored(KS, i, j) == (i <= j) + # @test is_stored(KS, i, j) == (i <= j) end for j in u2dofs, i in v2dofs + @test is_stored(sparsity_pattern, i, j) == (i == j) @test is_stored(K, i, j) == (i == j) - @test is_stored(KS, i, j) == (i == j) + # @test is_stored(KS, i, j) == (i == j) end for j in pdofs, i in v1dofs + @test !is_stored(sparsity_pattern, i, j) @test !is_stored(K, i, j) - @test !is_stored(KS, i, j) + # @test !is_stored(KS, i, j) end for j in pdofs, i in Iterators.flatten((v2dofs, qdofs)) + @test is_stored(sparsity_pattern, i, j) @test is_stored(K, i, j) - @test is_stored(KS, i, j) == (i <= j) + # @test is_stored(KS, i, j) == (i <= j) end # Error paths - @test_throws ErrorException("coupling not square") create_sparsity_pattern(dh; coupling=[true true]) - @test_throws ErrorException("coupling not symmetric") create_symmetric_sparsity_pattern(dh; coupling=[true true; false true]) - @test_throws ErrorException("could not create coupling") create_symmetric_sparsity_pattern(dh; coupling=falses(100, 100)) + @test_throws ErrorException("coupling not square") allocate_matrix(dh; coupling=[true true]) + # @test_throws ErrorException("coupling not symmetric") create_symmetric_sparsity_pattern(dh; coupling=[true true; false true]) + # @test_throws ErrorException("could not create coupling") create_symmetric_sparsity_pattern(dh; coupling=falses(100, 100)) # Test coupling with subdomains grid = generate_grid(Quadrilateral, (1, 2)) @@ -521,28 +540,35 @@ end sdh2 = SubDofHandler(dh, Set(2)) add!(sdh2, :u, Lagrange{RefQuadrilateral,1}()^2) close!(dh) - K = create_sparsity_pattern(dh; coupling = [true true; true false]) - KS = create_symmetric_sparsity_pattern(dh; coupling = [true true; true false]) + + sparsity_pattern = init_sparsity_pattern(dh) + add_sparsity_entries!(sparsity_pattern, dh; coupling = [true true; true false]) + K = allocate_matrix(sparsity_pattern) + KS = Symmetric(allocate_matrix(dh; #= symmetric=true, =# coupling = [true true; true false])) # Subdomain 1: u and p udofs = celldofs(dh, 1)[dof_range(sdh1, :u)] pdofs = celldofs(dh, 1)[dof_range(sdh1, :p)] for j in udofs, i in Iterators.flatten((udofs, pdofs)) + @test is_stored(sparsity_pattern, i, j) @test is_stored(K, i, j) - @test is_stored(KS, i, j) == (i <= j) + # @test is_stored(KS, i, j) == (i <= j) end for j in pdofs, i in udofs + @test is_stored(sparsity_pattern, i, j) @test is_stored(K, i, j) - @test is_stored(KS, i, j) + # @test is_stored(KS, i, j) end for j in pdofs, i in pdofs + @test is_stored(sparsity_pattern, i, j) == (i == j) @test is_stored(K, i, j) == (i == j) - @test is_stored(KS, i, j) == (i == j) + # @test is_stored(KS, i, j) == (i == j) end # Subdomain 2: u udofs = celldofs(dh, 2)[dof_range(sdh2, :u)] for j in udofs, i in udofs + @test is_stored(sparsity_pattern, i, j) @test is_stored(K, i, j) - @test is_stored(KS, i, j) == (i <= j) + # @test is_stored(KS, i, j) == (i <= j) end end @@ -624,18 +650,18 @@ end end end end - function check_coupling(dh, topology, K, coupling, cross_coupling) + function check_coupling(dh, topology, K, coupling, interface_coupling) for cell_idx in eachindex(getcells(dh.grid)) sdh = dh.subdofhandlers[dh.cell_to_subdofhandler[cell_idx]] coupling_idx = [1,1] - cross_coupling_idx = [1,1] + interface_coupling_idx = [1,1] vdim = [1,1] # test inner coupling _check_dofs(K, dh, sdh, cell_idx, coupling, coupling_idx, vdim, [cell_idx], false) # test cross-element coupling neighborhood = Ferrite.getrefdim(dh.grid.cells[1]) > 1 ? topology.face_face_neighbor : topology.vertex_vertex_neighbor neighbors = neighborhood[cell_idx, :] - _check_dofs(K, dh, sdh, cell_idx, cross_coupling, cross_coupling_idx, vdim, [i[1][1] for i in neighbors[.!isempty.(neighbors)]], true) + _check_dofs(K, dh, sdh, cell_idx, interface_coupling, interface_coupling_idx, vdim, [i[1][1] for i in neighbors[.!isempty.(neighbors)]], true) end end grid = generate_grid(Quadrilateral, (2, 2)) @@ -645,16 +671,16 @@ end add!(dh, :p, DiscontinuousLagrange{RefQuadrilateral,1}()) add!(dh, :w, Lagrange{RefQuadrilateral,1}()) close!(dh) - for coupling in couplings, cross_coupling in couplings - K = create_sparsity_pattern(dh; coupling=coupling, topology = topology, cross_coupling = cross_coupling) - all(coupling) && @test K == create_sparsity_pattern(dh, topology = topology, cross_coupling = cross_coupling) - check_coupling(dh, topology, K, coupling, cross_coupling) + for coupling in couplings, interface_coupling in couplings + K = allocate_matrix(dh; coupling=coupling, topology = topology, interface_coupling = interface_coupling) + all(coupling) && @test K == allocate_matrix(dh, topology = topology, interface_coupling = interface_coupling) + check_coupling(dh, topology, K, coupling, interface_coupling) end # Error paths - @test_throws ErrorException("coupling not square") create_sparsity_pattern(dh; coupling=[true true]) - @test_throws ErrorException("coupling not symmetric") create_symmetric_sparsity_pattern(dh; coupling=[true true; false true]) - @test_throws ErrorException("could not create coupling") create_symmetric_sparsity_pattern(dh; coupling=falses(100, 100)) + @test_throws ErrorException("coupling not square") allocate_matrix(dh; coupling=[true true]) + # @test_throws ErrorException("coupling not symmetric") allocate_matrix(dh; coupling=[true true; false true]) + @test_throws ErrorException("could not create coupling") allocate_matrix(dh; coupling=falses(100, 100)) # Test coupling with subdomains # Note: `check_coupling` works for this case only because the second domain has dofs from the first domain in order. Otherwise tests like in continuous ip are required. @@ -670,10 +696,10 @@ end add!(sdh2, :u, DiscontinuousLagrange{RefQuadrilateral,1}()^2) close!(dh) - for coupling in couplings, cross_coupling in couplings - K = create_sparsity_pattern(dh; coupling=coupling, topology = topology, cross_coupling = cross_coupling) - all(coupling) && @test K == create_sparsity_pattern(dh, topology = topology, cross_coupling = cross_coupling) - check_coupling(dh, topology, K, coupling, cross_coupling) + for coupling in couplings, interface_coupling in couplings + K = allocate_matrix(dh; coupling=coupling, topology = topology, interface_coupling = interface_coupling) + all(coupling) && @test K == allocate_matrix(dh, topology = topology, interface_coupling = interface_coupling) + check_coupling(dh, topology, K, coupling, interface_coupling) end # Testing Crouzeix-Raviart coupling @@ -683,9 +709,9 @@ end add!(dh, :u, CrouzeixRaviart{RefTriangle,1}()) close!(dh) coupling = trues(3,3) - K = create_sparsity_pattern(dh; coupling=coupling, topology = topology, cross_coupling = coupling) - K_cont = create_sparsity_pattern(dh; coupling=coupling, topology = topology, cross_coupling = falses(3,3)) - K_default = create_sparsity_pattern(dh) + K = allocate_matrix(dh; coupling=coupling, topology = topology, interface_coupling = coupling) + K_cont = allocate_matrix(dh; coupling=coupling, topology = topology, interface_coupling = falses(3,3)) + K_default = allocate_matrix(dh) @test K == K_cont == K_default end diff --git a/test/test_l2_projection.jl b/test/test_l2_projection.jl index 147368af52..bdb3df58d8 100644 --- a/test/test_l2_projection.jl +++ b/test/test_l2_projection.jl @@ -424,7 +424,7 @@ function test_export(;subset::Bool) end -function test_show() +function test_show_l2() grid = generate_grid(Triangle, (2,2)) ip = Lagrange{RefTriangle, 1}() proj = L2Projector(ip, grid) @@ -496,6 +496,6 @@ end test_projection_mixedgrid() test_export(subset=false) test_export(subset=true) - test_show() + test_show_l2() test_l2proj_errorpaths() end diff --git a/test/test_mixeddofhandler.jl b/test/test_mixeddofhandler.jl index cf9ad6dbf1..ac0e8ada1f 100644 --- a/test/test_mixeddofhandler.jl +++ b/test/test_mixeddofhandler.jl @@ -335,7 +335,7 @@ function test_2_element_heat_eq() end end - K = create_sparsity_pattern(dh) + K = allocate_matrix(dh) f = zeros(ndofs(dh)); assembler = start_assemble(K, f); # Use the same assemble function since it is the same weak form for both cell-types diff --git a/test/test_pointevaluation.jl b/test/test_pointevaluation.jl index eec7184797..a52c594fae 100644 --- a/test/test_pointevaluation.jl +++ b/test/test_pointevaluation.jl @@ -127,7 +127,7 @@ function dofhandler() end function _pointeval_dofhandler2_manual_projection(dh, csv, cvv, f_s, f_v) - M = create_sparsity_pattern(dh) + M = allocate_matrix(dh) f = zeros(ndofs(dh)) asm = start_assemble(M, f) me = zeros(ndofs_per_cell(dh), ndofs_per_cell(dh)) diff --git a/test/test_sparsity_patterns.jl b/test/test_sparsity_patterns.jl new file mode 100644 index 0000000000..d9f063d192 --- /dev/null +++ b/test/test_sparsity_patterns.jl @@ -0,0 +1,332 @@ +using Ferrite, Test, SparseArrays, Random + +# Minimal implementation of a custom sparsity pattern +struct TestPattern <: Ferrite.AbstractSparsityPattern + nrowscols::Tuple{Int, Int} + data::Vector{Vector{Int}} + function TestPattern(m::Int, n::Int) + return new((m, n), Vector{Int}[Int[] for _ in 1:m]) + end +end +Ferrite.getnrows(tp::TestPattern) = tp.nrowscols[1] +Ferrite.getncols(tp::TestPattern) = tp.nrowscols[2] +function Ferrite.add_entry!(tp::TestPattern, row::Int, col::Int) + if !(1 <= row <= tp.nrowscols[1] && 1 <= col <= tp.nrowscols[2]) + error("out of bounds") + end + r = tp.data[row] + k = searchsortedfirst(r, col) + if k == lastindex(r) + 1 || r[k] != col + insert!(r, k, col) + end + return +end +Ferrite.eachrow(tp::TestPattern) = tp.data +Ferrite.eachrow(tp::TestPattern, r::Int) = tp.data[r] + +function compare_patterns(p1, px...) + @test all(p -> Ferrite.getnrows(p1) == Ferrite.getnrows(p), px) + @test all(p -> Ferrite.getncols(p1) == Ferrite.getncols(p), px) + for rs in zip(Ferrite.eachrow.((p1, px...,))...) + for cs in zip(rs...) + @test all(c -> cs[1] == c, cs) + end + end +end + +# Compare the storage of SparseMatrixCSC +function compare_matrices(A1, Ax...) + @assert A1 isa SparseMatrixCSC + @assert length(Ax) > 0 + @assert all(A -> A isa SparseMatrixCSC, Ax) + @test all(A -> size(A1) == size(A), Ax) + @test all(A -> A1.colptr == A.colptr, Ax) + @test all(A -> A1.rowval == A.rowval, Ax) + return +end + +function is_stored(dsp::SparsityPattern, i, j) + return findfirst(k -> k == j, dsp.rows[i]) !== nothing +end + +@testset "SparsityPattern" begin + + # Ferrite.add_entry! + for (m, n) in ((5, 5), (3, 5), (5, 3)) + dsp = SparsityPattern(m, n) + @test Ferrite.getnrows(dsp) == m + @test Ferrite.getncols(dsp) == n + for r in randperm(m), c in randperm(n) + @test !is_stored(dsp, r, c) + Ferrite.add_entry!(dsp, r, c) + @test is_stored(dsp, r, c) + end + A = allocate_matrix(dsp) + fill!(A.nzval, 1) + @test A == ones(m, n) + # Error paths + @test_throws BoundsError Ferrite.add_entry!(dsp, 0, 1) + @test_throws BoundsError Ferrite.add_entry!(dsp, 1, 0) + @test_throws BoundsError Ferrite.add_entry!(dsp, m+1, 1) + @test_throws BoundsError Ferrite.add_entry!(dsp, 1, n+1) + end + + function testdhch() + local grid, dh, ch + grid = generate_grid(Quadrilateral, (2, 1)) + dh = DofHandler(grid) + add!(dh, :v, Lagrange{RefQuadrilateral,1}()^2) + add!(dh, :s, Lagrange{RefQuadrilateral,1}()) + close!(dh) + ch = ConstraintHandler(dh) + add!(ch, Dirichlet(:v, getfacetset(grid, "left"), (x, t) -> 0, [2])) + add!(ch, Dirichlet(:s, getfacetset(grid, "left"), (x, t) -> 0)) + add!(ch, AffineConstraint(15, [1 => 0.5, 7 => 0.5], 0.0)) + close!(ch) + return dh, ch + end + + dh, ch = testdhch() + + # Mismatching size + + # Test show method + dsp = SparsityPattern(ndofs(dh), ndofs(dh)) + str = sprint(show, "text/plain", dsp) + @test contains(str, "$(ndofs(dh))×$(ndofs(dh))") + @test contains(str, r" - Sparsity: 100.0% \(0 stored entries\)$"m) + @test contains(str, r" - Entries per row \(min, max, avg\): 0, 0, 0.0$"m) + @test contains(str, r" - Memory estimate: .* used, .* allocated$"m) + add_sparsity_entries!(dsp, dh) + str = sprint(show, "text/plain", dsp) + @test contains(str, "$(ndofs(dh))×$(ndofs(dh))") + @test contains(str, r" - Sparsity: .*% \(252 stored entries\)$"m) + @test contains(str, r" - Entries per row \(min, max, avg\): 12, 18, 14\.0$"m) + @test contains(str, r" - Memory estimate: .* used, .* allocated$"m) + + # Test all the possible entrypoints using SparsityPattern with a DofHandler + compare_matrices( + # Reference matrix from COO representation + let I = Int[], J = Int[] + for c in CellIterator(dh) + for row in c.dofs, col in c.dofs + push!(I, row) + push!(J, col) + end + end + sparse(I, J, zeros(Float64, length(I)), ndofs(dh), ndofs(dh)) + end, + let + A = allocate_matrix(dh) + @test A isa SparseMatrixCSC{Float64, Int} + A + end, + let + A = allocate_matrix(SparseMatrixCSC{Float32, Int}, dh) + @test A isa SparseMatrixCSC{Float32, Int} + A + end, + let + dsp = init_sparsity_pattern(dh) + add_sparsity_entries!(dsp, dh) + @test dsp isa SparsityPattern + A = allocate_matrix(dsp) + @test A isa SparseMatrixCSC{Float64, Int} + A + end, + let + dsp = init_sparsity_pattern(dh) + add_sparsity_entries!(dsp, dh) + A = allocate_matrix(SparseMatrixCSC{Float32, Int32}, dsp) + @test A isa SparseMatrixCSC{Float32, Int32} + A + end, + let + dsp = SparsityPattern(ndofs(dh), ndofs(dh); nnz_per_row = 5) + allocate_matrix(add_sparsity_entries!(dsp, dh)) + end, + ) + + # Test entrypoints with a DofHandler + ConstraintHandler + compare_matrices( + let + A = allocate_matrix(dh, ch) + @test A isa SparseMatrixCSC{Float64, Int} + A + end, + let + A = allocate_matrix(SparseMatrixCSC{Float32, Int}, dh, ch) + @test A isa SparseMatrixCSC{Float32, Int} + A + end, + let + dsp = init_sparsity_pattern(dh) + add_sparsity_entries!(dsp, dh, ch) + @test dsp isa SparsityPattern + A = allocate_matrix(dsp) + @test A isa SparseMatrixCSC{Float64, Int} + A + end, + let + dsp = init_sparsity_pattern(dh) + add_sparsity_entries!(dsp, dh, ch) + A = allocate_matrix(SparseMatrixCSC{Float32, Int32}, dsp) + @test A isa SparseMatrixCSC{Float32, Int32} + A + end, + let + dsp = SparsityPattern(ndofs(dh), ndofs(dh)) + allocate_matrix(add_sparsity_entries!(dsp, dh, ch)) + end, + ) + + # Test entrypoints with a DofHandler + coupling + remove constrained + kwargs = (; coupling = [true true; false true], keep_constrained = false) + compare_matrices( + let + A = allocate_matrix(dh, ch; kwargs...) + @test A isa SparseMatrixCSC{Float64, Int} + A + end, + let + A = allocate_matrix(SparseMatrixCSC{Float32, Int}, dh, ch; kwargs...) + @test A isa SparseMatrixCSC{Float32, Int} + A + end, + let + dsp = init_sparsity_pattern(dh) + add_sparsity_entries!(dsp, dh, ch; kwargs...) + @test dsp isa SparsityPattern + A = allocate_matrix(dsp) + @test A isa SparseMatrixCSC{Float64, Int} + A + end, + let + dsp = init_sparsity_pattern(dh) + add_sparsity_entries!(dsp, dh, ch; kwargs...) + A = allocate_matrix(SparseMatrixCSC{Float32, Int32}, dsp) + @test A isa SparseMatrixCSC{Float32, Int32} + A + end, + let + dsp = SparsityPattern(ndofs(dh), ndofs(dh)) + allocate_matrix(add_sparsity_entries!(dsp, dh, ch; kwargs...)) + end, + ) + +end + +@testset "Sparsity pattern generics" begin + + # Test setup + grid = generate_grid(Hexahedron, (5, 5, 5)) + dh = DofHandler(grid) + add!(dh, :u, Lagrange{RefHexahedron, 2}()^3) + add!(dh, :p, Lagrange{RefHexahedron, 1}()) + close!(dh) + ch = ConstraintHandler(dh) + add!(ch, Dirichlet(:p, union(getfacetset.(Ref(grid), ("left", "right", "top", "bottom", "front", "back"),)...), x -> 0)) + add!(ch, PeriodicDirichlet(:u, collect_periodic_facets(grid))) + close!(ch) + + function make_patterns(dh) + nd = ndofs(dh) + tp = TestPattern(nd, nd) + sp = SparsityPattern(nd, nd) + bp = BlockSparsityPattern([nd ÷ 2, nd - nd ÷ 2]) + return tp, sp, bp + end + + # DofHandler + ps = make_patterns(dh) + for p in ps + add_sparsity_entries!(p, dh) + end + compare_patterns(ps...) + + # DofHandler + ConstraintHandler + ps = make_patterns(dh) + for p in ps + add_sparsity_entries!(p, dh, ch) + end + compare_patterns(ps...) + + # DofHandler + ConstraintHandler later + ps = make_patterns(dh) + for p in ps + add_sparsity_entries!(p, dh) + add_constraint_entries!(p, ch) + end + compare_patterns(ps...) + + # Individual pieces + ps = make_patterns(dh) + for p in ps + add_cell_entries!(p, dh) + add_constraint_entries!(p, ch) + end + compare_patterns(ps...) + + # Ignore constrained dofs + ps = make_patterns(dh) + for p in ps + add_sparsity_entries!(p, dh, ch; keep_constrained=false) + # Test that prescribed dofs only have diagonal entry + for row in ch.prescribed_dofs + r = Ferrite.eachrow(p, row) + col, state = iterate(r) + @test col == row + @test iterate(r, state) === nothing + end + end + compare_patterns(ps...) + + # Coupling + ps = make_patterns(dh) + for p in ps + add_sparsity_entries!(p, dh, ch; coupling = [true true; false true]) + end + compare_patterns(ps...) + ps = make_patterns(dh) + for p in ps + coupling = ones(Bool, 4, 4) + coupling[4, 1:3] .= false + add_sparsity_entries!(p, dh, ch; coupling = coupling) + end + compare_patterns(ps...) + ps = make_patterns(dh) + for p in ps + coupling = ones(Bool, ndofs_per_cell(dh), ndofs_per_cell(dh)) + coupling[1:2:(ndofs_per_cell(dh) ÷ 2), :] .= false + add_sparsity_entries!(p, dh, ch; coupling = coupling) + end + compare_patterns(ps...) + + # Error paths + dh_open = DofHandler(grid) + for p in (SparsityPattern(2, 2), TestPattern(2, 2), BlockSparsityPattern([1, 1])) + @test_throws ErrorException("the DofHandler must be closed") add_sparsity_entries!(p, dh_open) + end + for p in (SparsityPattern(ndofs(dh), 2), TestPattern(ndofs(dh), 2), BlockSparsityPattern([2, 2])) + @test_throws ErrorException add_sparsity_entries!(p, dh) + end + for p in (SparsityPattern(2, ndofs(dh)), TestPattern(2, ndofs(dh))) + @test_throws ErrorException add_sparsity_entries!(p, dh) + end + patterns = ( + SparsityPattern(ndofs(dh), ndofs(dh)), + TestPattern(ndofs(dh), ndofs(dh)), + BlockSparsityPattern([ndofs(dh) ÷ 2, ndofs(dh) - ndofs(dh) ÷ 2]), + ) + for p in patterns + @test_throws ErrorException add_sparsity_entries!(p, dh; keep_constrained=false) + end + ch_open = ConstraintHandler(dh) + for p in patterns + @test_throws ErrorException add_sparsity_entries!(p, dh, ch_open; keep_constrained=false) + end + ch_bad = ConstraintHandler(close!(DofHandler(grid))) + for p in patterns + @test_throws ErrorException add_sparsity_entries!(p, dh, ch_bad; keep_constrained=false) + end +end