From 7b82c7ff5f3f80171c9426c68ea66cc072e766a6 Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Tue, 18 Jun 2024 23:45:10 +0200 Subject: [PATCH] Implement support for working with sparsity patterns This patch changes the sparsity pattern interface to decouple the pattern creation from matrix instantiation. This will make it easier to support many different types of matrices from the same sparstiy pattern data structure. Main additions: - New sparsity pattern datastructures `SparsityPattern` and `BlockSparsityPattern`. - New functions for adding entries: `add_cell_entries!`, `add_interface_entries!`, `add_constraint_entries!`, and `Ferrite.add_entry!`. There is also `create_sparsity_pattern!` which calls the three first methods. - New function `allocate_matrix` which instantiates a matrix of chosen type from a sparsity pattern. - `allocate_matrix` is also a convenience method that can be used when no customization of the patterns is necessary, therefore, the new `allocate_matrix` can be used as a direct replacement for the old `create_sparsity_pattern`. --- CHANGELOG.md | 33 +- benchmark/helper.jl | 3 +- docs/Manifest.toml | 4 +- docs/liveserver.jl | 4 +- docs/make.jl | 7 +- docs/src/devdocs/dofhandler.md | 2 - docs/src/literate-gallery/helmholtz.jl | 2 +- docs/src/literate-gallery/landau.jl | 4 +- .../quasi_incompressible_hyperelasticity.jl | 2 +- .../literate-gallery/topology_optimization.jl | 2 +- docs/src/literate-howto/threaded_assembly.jl | 2 +- .../computational_homogenization.jl | 4 +- .../literate-tutorials/dg_heat_equation.jl | 2 +- docs/src/literate-tutorials/heat_equation.jl | 4 +- .../src/literate-tutorials/hyperelasticity.jl | 4 +- .../incompressible_elasticity.jl | 2 +- docs/src/literate-tutorials/linear_shell.jl | 2 +- docs/src/literate-tutorials/ns_vs_diffeq.jl | 4 +- docs/src/literate-tutorials/plasticity.jl | 2 +- docs/src/literate-tutorials/porous_media.jl | 2 +- .../literate-tutorials/reactive_surface.jl | 4 +- docs/src/literate-tutorials/stokes-flow.jl | 2 +- .../transient_heat_equation.jl | 4 +- docs/src/reference/assembly.md | 5 - docs/src/reference/sparsity_pattern.md | 57 ++ docs/src/topics/assembly.md | 28 +- docs/src/topics/constraints.md | 2 +- docs/src/topics/sparse_matrix.md | 165 ++++ ext/FerriteBlockArrays.jl | 91 +- ext/FerriteMetis.jl | 2 +- src/Dofs/block_sparsity_pattern.jl | 130 +++ src/Dofs/sparsity_pattern.jl | 795 +++++++++++++----- src/Ferrite.jl | 5 +- src/HeapAllocator.jl | 315 +++++++ src/L2_projection.jl | 9 +- src/deprecations.jl | 2 + src/exports.jl | 13 +- test/HeapAllocator.jl | 111 +++ test/blockarrays.jl | 17 +- .../test_simple_scalar_convergence.jl | 2 +- test/runtests.jl | 2 + test/test_abstractgrid.jl | 2 +- test/test_apply_rhs.jl | 4 +- test/test_constraints.jl | 24 +- test/test_dofs.jl | 122 +-- test/test_l2_projection.jl | 4 +- test/test_mixeddofhandler.jl | 2 +- test/test_pointevaluation.jl | 2 +- test/test_sparsity_patterns.jl | 332 ++++++++ 49 files changed, 1941 insertions(+), 403 deletions(-) create mode 100644 docs/src/reference/sparsity_pattern.md create mode 100644 docs/src/topics/sparse_matrix.md create mode 100644 src/Dofs/block_sparsity_pattern.jl create mode 100644 src/HeapAllocator.jl create mode 100644 test/HeapAllocator.jl create mode 100644 test/test_sparsity_patterns.jl 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/Manifest.toml b/docs/Manifest.toml index 80623f405c..de55bee597 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -1,6 +1,6 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.10.3" +julia_version = "1.10.4" manifest_format = "2.0" project_hash = "8722b6a67abf1b81e2c7808ddfcfb1d737cd0040" @@ -1297,7 +1297,7 @@ version = "1.3.5+1" [[deps.OpenBLAS_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.23+2" +version = "0.3.23+4" [[deps.OpenLibm_jll]] deps = ["Artifacts", "Libdl"] 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 9db1554cb7..84633be4df 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..003a88e5f8 --- /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 +AbstractSparsityPattern +create_sparsity_pattern! +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}, ::AbstractSparsityPattern) where {Tv, Ti, S <: SparseMatrixCSC{Tv, Ti}} +allocate_matrix(::Type{Symmetric{Tv, S}}, ::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..4d205cfad8 --- /dev/null +++ b/docs/src/topics/sparse_matrix.md @@ -0,0 +1,165 @@ +# [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 *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. + + +!!! 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 datastructures 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. 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: + - [`create_sparsity_pattern!`](@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 + +!!! todo + Write text. diff --git a/ext/FerriteBlockArrays.jl b/ext/FerriteBlockArrays.jl index 90e50a25f2..d9e77c25bd 100644 --- a/ext/FerriteBlockArrays.jl +++ b/ext/FerriteBlockArrays.jl @@ -1,60 +1,59 @@ 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 using Ferrite: addindex!, fillzero! +using SparseArrays: SparseMatrixCSC -# 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 -################################### -## 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) +############################## +## Instantiating the matrix ## +############################## + +# 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 d2a82935ad..ed10563c7c 100644 --- a/ext/FerriteMetis.jl +++ b/ext/FerriteMetis.jl @@ -39,7 +39,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..65b8494efb --- /dev/null +++ b/src/Dofs/block_sparsity_pattern.jl @@ -0,0 +1,130 @@ +######################## +# 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): + + - [`create_sparsity_pattern!`](@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) + 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 + +n_rows(bsp::BlockSparsityPattern) = bsp.nrows +n_cols(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 <= n_rows(bsp) && 1 <= col <= n_cols(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 <= n_rows(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:n_rows(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..a8c0262b7b 100644 --- a/src/Dofs/sparsity_pattern.jl +++ b/src/Dofs/sparsity_pattern.jl @@ -1,83 +1,473 @@ +########################### +# AbstractSparsityPattern # +########################### + +""" + AbstractSparsityPattern + +Supertype for sparsity pattern implementations, e.g. [`SparsityPattern`](@ref) and +[`BlockSparsityPattern`](@ref). +""" +abstract type AbstractSparsityPattern end + +""" + n_rows(sp::AbstractSparsityPattern) + +Return the number of rows in the sparsity pattern `sp`. +""" +n_rows(sp::AbstractSparsityPattern) + +""" + n_cols(sp::AbstractSparsityPattern) + +Return the number of columns in the sparsity pattern `sp`. +""" +n_cols(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) + +""" + eachrow(sp::AbstractSparsityPattern, row::Int) + +Return an iterator over *column* indices in row `row` of the sparsity pattern. + +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) + + +################### +# SparsityPattern # +################### + +""" + 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 + heap::HeapAllocator.Heap + rows::Vector{HeapAllocator.HeapVector{Int}} +end + +""" + 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): + - [`create_sparsity_pattern!`](@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) + heap = HeapAllocator.Heap() + rows = Vector{HeapAllocator.HeapVector{Int}}(undef, nrows) + for i in 1:nrows + rows[i] = HeapAllocator.resize(HeapAllocator.alloc_array(heap, Int, nnz_per_row), 0) + end + sp = SparsityPattern(nrows, ncols, heap, rows) + return sp +end + +function Base.show(io::IO, ::MIME"text/plain", sp::SparsityPattern) + iob = IOBuffer() + println(iob, "$(n_rows(sp))×$(n_cols(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( + (n_rows(sp) * n_cols(sp) - stored_entries) / (n_rows(sp) * n_cols(sp)) * 100 * 1000 + ) / 1000 + println(iob, " - Sparsity: $(sparsity_pct)% ($(stored_entries) stored entries)") + # Print row stats + avg_entries = round(stored_entries / n_rows(sp) * 10) / 10 + println(iob, " - Entries per row (min, max, avg): $(min_entries), $(max_entries), $(avg_entries)") + # Compute memory estimate + @assert n_rows(sp) * sizeof(eltype(sp.rows)) == sizeof(sp.rows) + bytes_used = sizeof(sp.rows) + stored_entries * sizeof(Int) + bytes_allocated = sizeof(sp.rows) + HeapAllocator.heap_stats(sp.heap)[2] + print(iob, " - Memory estimate: $(Base.format_bytes(bytes_used)) used, $(Base.format_bytes(bytes_allocated)) allocated") + write(io, seekstart(iob)) + return +end + +n_rows(sp::SparsityPattern) = sp.nrows +n_cols(sp::SparsityPattern) = sp.ncols + +@inline function add_entry!(sp::SparsityPattern, row::Int, col::Int) + @boundscheck (1 <= row <= n_rows(sp) && 1 <= col <= n_cols(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::HeapAllocator.HeapVector{Int}, item::Int) + k = searchsortedfirst(x, item) + if k == length(x) + 1 || @inbounds(x[k]) != item + x = HeapAllocator.insert(x, k, item) + end + return x +end + +eachrow(sp::SparsityPattern) = sp.rows +eachrow(sp::SparsityPattern, row::Int) = sp.rows[row] + + +################################################ +## 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 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_sparsity_pattern!( + 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. """ - create_sparsity_pattern(dh::DofHandler; coupling, [topology::Union{Nothing, AbstractTopology}], [cross_coupling]) +function create_sparsity_pattern!( + 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 n_rows(sp) < ndofs(dh) || n_cols(sp) < ndofs(dh) + error("number of rows ($(n_rows(sp))) or columns ($(n_cols(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 + # return _create_sparsity_pattern!(sp, dh, ch, keep_constrained, coupling, + # interface_coupling, topology) +end -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. -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. +# # Arguments +# - `sp`: the sparsity pattern to add entries to +# - `dh`: the DofHandler for which to generate the pattern +# - `ch`: the ConstraintHandler. Used to add the non-zero entries resulting from +# `AffineConstraint`s (and `PeriodicDirichlet` constraints, which uses `AffineConstraint`s +# internally). See [`add_constraint_entries!`](@ref). +# # 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. +# - `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. -If `topology` and `cross_coupling` are passed, dof of fields with discontinuous interpolations are coupled between elements according to `cross_coupling`. +""" + 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 -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) + 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, 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 + +""" + 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(n_rows(sp), n_cols(sp)) + add_entry!(sp, d, d) + end + return sp +end + + +############################################################ +# Sparse matrix instantiation from AbstractSparsityPattern # +############################################################ + """ - create_symmetric_sparsity_pattern(dh::DofHandler; coupling, topology::Union{Nothing, AbstractTopology}, cross_coupling) + allocate_matrix(::Type{SparseMatrixCSC{Tv, Ti}}, sp::SparsityPattern) -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}`. +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) -See the [Sparsity Pattern](@ref) section of the manual. +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 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 allocate_matrix(::Type{Symmetric{Tv, S}}, sp::AbstractSparsityPattern) where {Tv, Ti, S <: SparseMatrixCSC{Tv, Ti}} + return Symmetric(_allocate_matrix(S, sp, #=sym=# true)) end """ - create_symmetric_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler; coupling, topology::Union{Nothing, AbstractTopology}, cross_coupling) + allocate_matrix(sp::SparsityPattern) + +Allocate a sparse matrix of type `SparseMatrixCSC{Float64, Int}` from the sparsity pattern +`sp`. -Create a symmetric sparsity pattern accounting for affine constraints in `ch`. See -the Affine Constraints section of the manual for further details. +This method is a shorthand for the equivalent +[`allocate_matrix(SparseMatrixCSC{Float64, Int}, sp)`] +(@ref allocate_matrix(::Type{S}, sp::AbstractSparsityPattern) where {Tv, Ti, S <: SparseMatrixCSC{Tv, Ti}}). """ -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) +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) +create_sparsity_pattern!(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) + create_sparsity_pattern!(sp, dh) + K = allocate_matrix(sp) + M = allocate_matrix(sp) + ``` + instead of + ```julia + K = allocate_matrix(dh) + M = allocate_matrix(dh) + ``` +""" +function allocate_matrix(::Type{MatrixType}, dh::DofHandler, args...; kwargs...) where {MatrixType} + sp = init_sparsity_pattern(dh) + create_sparsity_pattern!(sp, dh, args...; kwargs...) + return allocate_matrix(MatrixType, sp) end """ - create_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler; coupling, topology::Union{Nothing, AbstractTopology} = nothing) + allocate_matrix(dh::DofHandler, args...; kwargs...) + +Allocate a matrix of type `SparseMatrixCSC{Float64, Int}` from the DofHandler `dh`. -Create a sparsity pattern accounting for affine constraints in `ch`. See -the Affine Constraints section of the manual for further details. +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 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) +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}[] field_dims = map(fieldname -> n_components(dh, fieldname), dh.field_names) - for sdh in dh.subdofhandlers - out = zeros(Bool, ndofs_per_cell(sdh), ndofs_per_cell(sdh)) + for fh in dh.subdofhandlers + out = zeros(Bool, ndofs_per_cell(fh), ndofs_per_cell(fh)) push!(outs, out) - dof_ranges = [dof_range(sdh, f) for f in sdh.field_names] - global_idxs = [findfirst(x -> x === f, dh.field_names) for f in sdh.field_names] + dof_ranges = [dof_range(fh, f) for f in fh.field_names] + global_idxs = [findfirst(x -> x === f, dh.field_names) for f in fh.field_names] if sz == length(dh.field_names) # Coupling given by fields for (j, jrange) in pairs(dof_ranges), (i, irange) in pairs(dof_ranges) @@ -92,7 +482,7 @@ function _coupling_to_local_dof_coupling(dh::DofHandler, coupling::AbstractMatri out[I, J] = coupling[ic, jc] end end - elseif sz == ndofs_per_cell(sdh) # Coupling given by template local matrix + elseif sz == ndofs_per_cell(fh) # Coupling given by template local matrix # TODO: coupling[fieldhandler_idx] if different template per subddomain out .= coupling else @@ -102,38 +492,139 @@ 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. + heap = HeapAllocator.Heap() + sp′ = Dict{Int, HeapAllocator.HeapVector{Int}}() + + for (row, colidxs) in zip(1:n_rows(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 + HeapAllocator.resize(HeapAllocator.alloc_array(heap, Int, 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 + HeapAllocator.resize(HeapAllocator.alloc_array(heap, Int, 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 + HeapAllocator.resize(HeapAllocator.alloc_array(heap, Int, 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 + + # Release the memory in the temporary heap + HeapAllocator.free(heap) + + 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 +635,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 +644,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) - 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 +# 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, n_cols(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 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 + 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:n_rows(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 - - 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:n_cols(sp)) + S = SparseMatrixCSC(n_rows(sp), n_cols(sp), colptr, rowval, nzval) + return S end diff --git a/src/Ferrite.jl b/src/Ferrite.jl index 5e794c9c93..2e288754ce 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("HeapAllocator.jl") # Matrix/Vector utilities include("arrayutils.jl") @@ -139,6 +139,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/HeapAllocator.jl b/src/HeapAllocator.jl new file mode 100644 index 0000000000..0927b62360 --- /dev/null +++ b/src/HeapAllocator.jl @@ -0,0 +1,315 @@ +module HeapAllocator + +@eval macro $(Symbol("const"))(field) + if VERSION >= v"1.8.0-DEV.1148" + Expr(:const, esc(field)) + else + return esc(field) + end +end + +const MALLOC_PAGE_SIZE = 4 * 1024 * 1024 % UInt # 4 MiB + +# Like Ptr{T} but also stores the number of bytes allocated +struct SizedPtr{T} + ptr::Ptr{T} + size::UInt +end + +SizedPtr{T}(ptr::SizedPtr) where {T} = SizedPtr{T}(Ptr{T}(ptr.ptr), ptr.size) +Ptr{T}(ptr::SizedPtr) where {T} = Ptr{T}(ptr.ptr) + +function memcpy(dst::SizedPtr, src::SizedPtr) + @assert dst.size >= src.size + ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), dst.ptr, src.ptr, src.size) + return +end + +# A page corresponds to a larger Libc.malloc call (MALLOC_PAGE_SIZE). Each page +# is split into smaller blocks to minimize the number of Libc.malloc/Libc.free +# calls. +mutable struct Page + @const ptr::SizedPtr{UInt8} # malloc'd pointer + @const blocksize::UInt # blocksize for this page + @const freelist::BitVector # block is free/used + free::SizedPtr{UInt8} # head of the free-list + n_free::UInt # number of free blocks +end + +function Page(ptr::SizedPtr{UInt8}, blocksize::UInt) + n_blocks, r = divrem(ptr.size, blocksize) + @assert r == 0 + # The first free pointer is the first block + free = SizedPtr(ptr.ptr, blocksize) + # Set up the free list + for i in 0:(n_blocks - 1) + ptr_this = ptr.ptr + i * blocksize + ptr_next = i == n_blocks - 1 ? Ptr{UInt8}() : ptr_this + blocksize + unsafe_store!(Ptr{Ptr{UInt8}}(ptr_this), ptr_next) + end + return Page(ptr, blocksize, trues(n_blocks), free, n_blocks) +end + +function _malloc(page::Page, size::UInt) + if size != page.blocksize + error("malloc: requested size does not match the blocksize") + end + # Return early with null if the page is full + if page.n_free == 0 + @assert page.free.ptr == C_NULL + return SizedPtr{UInt8}(Ptr{UInt8}(), size) + end + # Read the pointer to be returned + ret = page.free + @assert ret.ptr != C_NULL + # Make sure the pointer is not in use + blockindex = (ret.ptr - page.ptr.ptr) ÷ page.blocksize + 1 + if !@inbounds(page.freelist[blockindex]) + error("malloc: pointer already in use") + end + @inbounds page.freelist[blockindex] = false + # Look up and store the next free pointer + page.free = SizedPtr{UInt8}(unsafe_load(Ptr{Ptr{UInt8}}(ret)), size) + page.n_free -= 1 + return ret +end + +function _free(page::Page, ptr::SizedPtr{UInt8}) + if !(ptr.size == page.blocksize && page.ptr.ptr <= ptr.ptr < page.ptr.ptr + page.ptr.size) + error("free: not allocated in this page") + end + blockindex = (ptr.ptr - page.ptr.ptr) ÷ page.blocksize + 1 + if @inbounds page.freelist[blockindex] + error("free: double free") + end + @inbounds page.freelist[blockindex] = true + # Write the current free pointer to the pointer to be freed + unsafe_store!(Ptr{Ptr{UInt8}}(ptr), Ptr{UInt8}(page.free)) + # Store the just-freed pointer and increment the availability counter + page.free = ptr + page.n_free += 1 + # TODO: If this page is completely unused it can be collected and reused. + return +end + +# Collection of pages for a specific size +struct FixedSizeHeap + blocksize::UInt + pages::Vector{Page} +end + +function _malloc(fheap::FixedSizeHeap, size::UInt) + @assert fheap.blocksize == size + # Try all existing pages + for page in fheap.pages + ptr = _malloc(page, size) + ptr.ptr == C_NULL || return ptr + end + # Allocate a new page + # TODO: Replace Libc.malloc with Memory in recent Julias? + mptr = SizedPtr{UInt8}(Libc.malloc(MALLOC_PAGE_SIZE), MALLOC_PAGE_SIZE) + page = Page(mptr, fheap.blocksize) + push!(fheap.pages, page) + # Allocate block in the new page + ptr = _malloc(page, size) + @assert ptr.ptr != C_NULL + return ptr +end + +mutable struct Heap + @const size_heaps::Vector{FixedSizeHeap} # 2, 4, 6, 8, ... + function Heap() + heap = new(FixedSizeHeap[]) + finalizer(free, heap) + return heap + end +end + +function free(heap::Heap) + for i in 1:length(heap.size_heaps) + isassigned(heap.size_heaps, i) || continue + for page in heap.size_heaps[i].pages + Libc.free(page.ptr.ptr) + end + end + resize!(heap.size_heaps, 0) + return +end + +function heap_stats(heap::Heap) + bytes_used = 0 + bytes_allocated = 0 + for heapidx in 1:length(heap.size_heaps) + isassigned(heap.size_heaps, heapidx) || continue + size_heap = heap.size_heaps[heapidx] + bytes_allocated += length(size_heap.pages) * Int(MALLOC_PAGE_SIZE) + for page in size_heap.pages + bytes_used += count(!, page.freelist) * page.blocksize + end + end + return bytes_used, bytes_allocated +end + +function Base.show(io::IO, ::MIME"text/plain", heap::Heap) + iob = IOBuffer() + println(iob, "HeapAllocator.Heap with blockheaps:") + for idx in 1:length(heap.size_heaps) + isassigned(heap.size_heaps, idx) || continue + h = heap.size_heaps[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 * MALLOC_PAGE_SIZE ÷ blocksize + println(iob, " - Blocksize: $(blocksize), npages: $(npages), usage: $(n_tot - n_free) / $(n_tot)") + end + write(io, seekstart(iob)) + return +end + +function heapindex_from_blocksize(blocksize::UInt) + return (8 * sizeof(UInt) - leading_zeros(blocksize)) % Int +end + +function malloc(heap::Heap, size::UInt) + size = max(size, sizeof(Ptr{UInt8})) + blocksize = nextpow(2, size) + fheapidx = heapindex_from_blocksize(blocksize) + if length(heap.size_heaps) < fheapidx + resize!(heap.size_heaps, fheapidx) + end + if !isassigned(heap.size_heaps, fheapidx) + fheap = FixedSizeHeap(blocksize, Page[]) + heap.size_heaps[fheapidx] = fheap + else + fheap = heap.size_heaps[fheapidx] + end + return _malloc(fheap, blocksize) +end + +malloc(heap::Heap, size::Integer) = malloc(heap, size % UInt) + +function malloc(heap::Heap, ::Type{T}, size::Integer) where T + @assert isbitstype(T) + ptr = malloc(heap, (sizeof(T) * size) % UInt) + return SizedPtr{T}(ptr) +end + + +@inline function find_page(heap::Heap, ptr::SizedPtr{UInt8}) + heapidx = heapindex_from_blocksize(ptr.size) + if !isassigned(heap.size_heaps, heapidx) + error("pointer not malloc'd in this heap") + end + size_heap = heap.size_heaps[heapidx] + # Search for the page containing the pointer + # TODO: Insert pages in sorted order and use searchsortedfirst? + # TODO: Align the pages and store the base pointer -> page index mapping in a dict? + pageidx = findfirst(p -> p.ptr.ptr <= ptr.ptr < (p.ptr.ptr + p.ptr.size), size_heap.pages) + if pageidx === nothing + error("pointer not malloc'd in this heap") + end + return size_heap.pages[pageidx] +end + +function free(heap::Heap, ptr::SizedPtr{UInt8}) + _free(find_page(heap, ptr), ptr) + return +end + +function free(heap::Heap, ptr::SizedPtr) + return free(heap, SizedPtr{UInt8}(ptr)) +end + +function realloc(heap::Heap, ptr::SizedPtr{UInt8}, newsize::UInt) + @assert newsize > ptr.size # TODO: Allow shrinkage? + # Find the page for the pointer to make sure it was allocated in this heap + page = find_page(heap, ptr) + # Allocate the new pointer + ptr′ = malloc(heap, newsize) + # Copy the data + memcpy(ptr′, ptr) + # Free the old pointer and return + _free(page, ptr) + return ptr′ +end + +# Rewrap to the same pointer type +function realloc(heap::Heap, ptr::SizedPtr{T}, newsize::Integer) where T + ptr′ = realloc(heap, SizedPtr{UInt8}(ptr), (newsize * sizeof(T)) % UInt) + return SizedPtr{T}(ptr′) +end + +# Minimal AbstractVector implementation on top of a raw pointer + a length +struct HeapArray{T, N} <: AbstractArray{T, N} + ptr::SizedPtr{T} + size::NTuple{N, Int} + heap::Heap +end + +const HeapVector{T} = HeapArray{T, 1} + +Base.size(mv::HeapArray) = mv.size +allocated_length(mv::HeapArray{T}) where T = mv.ptr.size ÷ sizeof(T) +Base.IndexStyle(::Type{<:HeapArray}) = IndexLinear() +Base.@propagate_inbounds function Base.getindex(mv::HeapArray, i::Int) + @boundscheck checkbounds(mv, i) + return unsafe_load(mv.ptr.ptr, i) +end +Base.@propagate_inbounds function Base.setindex!(mv::HeapArray{T}, v::T, i::Int) where T + @boundscheck checkbounds(mv, i) + unsafe_store!(mv.ptr.ptr, v, i) + return mv +end + +function alloc_array(heap::Heap, ::Type{T}, size::NTuple{N, Int}) where {T, N} + ptr = malloc(heap, T, prod(size)) + return HeapArray{T, N}(ptr, size, heap) +end + +# Allow vararg dims +function alloc_array(heap::Heap, ::Type{T}, size::Int) where T + return alloc_array(heap, T, (size, )) +end +function alloc_array(heap::Heap, ::Type{T}, size1::Int, size2::Int, sizes::Int...) where T + size = (size1, size2, map(Int, sizes)..., ) + return alloc_array(heap, T, size) +end + +function realloc(heap::Heap, x::HeapVector{T}, size::Int) where T + @assert heap === x.heap + ptr = realloc(heap, x.ptr, size) + return HeapVector{T}(ptr, (size, ), heap) +end + +realloc(x::HeapVector, n::Int) = realloc(x.heap, x, n) + +function free(heap::Heap, x::HeapVector) + @assert heap === x.heap + return free(heap, x.ptr) +end + +free(x::HeapVector) = free(x.heap, x) + +@inline function resize(x::HeapVector{T}, n::Int) where T + if n > allocated_length(x) + return realloc(x, n) + else + return HeapVector{T}(x.ptr, (n, ), x.heap) + end +end + +@inline function insert(x::HeapVector{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 HeapAllocator diff --git a/src/L2_projection.jl b/src/L2_projection.jl index ec0a15a0fb..e3d86fff76 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,8 @@ 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) + n = getnbasefunctions(fe_values) + M = Symmetric(allocate_matrix(create_sparsity_pattern!(init_sparsity_pattern(dh; nnz_per_row = 2 * n), dh))) assembler = start_assemble(M) for (sdh, qr_lhs) in zip(dh.subdofhandlers, qrs_lhs) ip_fun = only(sdh.field_interpolations) @@ -160,7 +161,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 +287,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/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..8b433f9733 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, + create_sparsity_pattern!, + add_cell_entries!, + add_interface_entries!, + add_constraint_entries!, + allocate_matrix, + # Constraints ConstraintHandler, Dirichlet, diff --git a/test/HeapAllocator.jl b/test/HeapAllocator.jl new file mode 100644 index 0000000000..5d55165a75 --- /dev/null +++ b/test/HeapAllocator.jl @@ -0,0 +1,111 @@ +using Test, Ferrite.HeapAllocator + +@testset "HeapAllocator.jl" begin + + # SizedPtr{T} + sz = 1024 % UInt + mptr = Ptr{UInt8}(Libc.malloc(sz)) + ptr = HeapAllocator.SizedPtr(mptr, sz) + @test Ptr{UInt8}(ptr) === ptr.ptr === mptr + @test Ptr{Int}(ptr) === Ptr{Int}(ptr.ptr) === Ptr{Int}(mptr) + @test ptr.size == sz + + # Basic malloc, realloc, free + heap = HeapAllocator.Heap() + ptr = HeapAllocator.malloc(heap, 1024) + @test ptr isa HeapAllocator.SizedPtr{UInt8} + ptr′ = HeapAllocator.realloc(heap, ptr, 2048) + @test ptr′ isa HeapAllocator.SizedPtr{UInt8} + @test_throws ErrorException("free: double free") HeapAllocator.free(heap, ptr) + HeapAllocator.free(heap, ptr′) + @test_throws ErrorException("free: double free") HeapAllocator.free(heap, ptr′) + + # Internal page allocation: exhaust some pages + heap = HeapAllocator.Heap() + ptrs = HeapAllocator.SizedPtr{UInt8}[] + for _ in 1:(HeapAllocator.MALLOC_PAGE_SIZE ÷ 512 * 2 + 1) + ptr = HeapAllocator.malloc(heap, 512) + push!(ptrs, ptr) + end + @test length(heap.size_heaps[10].pages) == 3 + @test all(!, heap.size_heaps[10].pages[1].freelist) + @test all(!, heap.size_heaps[10].pages[2].freelist) + @test !heap.size_heaps[10].pages[3].freelist[1] + @test all(heap.size_heaps[10].pages[3].freelist[2:end]) + ptrs′ = HeapAllocator.SizedPtr{UInt8}[] + for ptr in ptrs + ptr′ = HeapAllocator.realloc(heap, ptr, 1024) + @test_throws ErrorException("free: double free") HeapAllocator.free(heap, ptr) + push!(ptrs′, ptr′) + end + @test length(heap.size_heaps[10].pages) == 3 # TODO + @test all(heap.size_heaps[10].pages[1].freelist) + @test all(heap.size_heaps[10].pages[2].freelist) + @test all(heap.size_heaps[10].pages[3].freelist) + @test length(heap.size_heaps[11].pages) == 5 + @test all(!, heap.size_heaps[11].pages[1].freelist) + @test all(!, heap.size_heaps[11].pages[2].freelist) + @test all(!, heap.size_heaps[11].pages[3].freelist) + @test all(!, heap.size_heaps[11].pages[4].freelist) + @test !heap.size_heaps[11].pages[5].freelist[1] + @test all(heap.size_heaps[11].pages[5].freelist[2:end]) + for ptr in ptrs′ + HeapAllocator.free(heap, ptr) + @test_throws ErrorException("free: double free") HeapAllocator.free(heap, ptr) + end + + # malloc, realloc, free with custom type + heap = HeapAllocator.Heap() + ptr = HeapAllocator.malloc(heap, Int, 1024) + @test ptr isa HeapAllocator.SizedPtr{Int} + @test ptr.size == 1024 * sizeof(Int) + ptr′ = HeapAllocator.realloc(heap, ptr, 2048) + @test ptr′ isa HeapAllocator.SizedPtr{Int} + @test ptr′.size == 2048 * sizeof(Int) + @test_throws ErrorException("free: double free") HeapAllocator.free(heap, ptr) + HeapAllocator.free(heap, ptr′) + @test_throws ErrorException("free: double free") HeapAllocator.free(heap, ptr′) + + # malloc, realloc, free with arrays + heap = HeapAllocator.Heap() + x = HeapAllocator.alloc_array(heap, Int, 1024) + x .= 1:1024 + x′ = HeapAllocator.realloc(heap, x, 2048) + @test_throws ErrorException("free: double free") HeapAllocator.free(x) + @test_throws ErrorException("free: double free") HeapAllocator.free(heap, x) + @test_throws ErrorException("free: double free") HeapAllocator.free(heap, x.ptr) + @test x′[1:1024] == 1:1024 + + # Array functions + heap = HeapAllocator.Heap() + x = HeapAllocator.alloc_array(heap, Int, 8) + x .= 1:8 + @test length(x) == 8 + x = HeapAllocator.resize(x, 0) + @test length(x) == 0 + x = HeapAllocator.resize(x, 16) + @test length(x) == 16 + @test x[1:8] == 1:8 + x = HeapAllocator.resize(x, 8) + x = HeapAllocator.insert(x, 1, -1) + x = HeapAllocator.insert(x, length(x) + 1, -1) + x = HeapAllocator.insert(x, 2, -2) + x = HeapAllocator.insert(x, length(x), -2) + @test x == [-1; -2; 1:8; -2; -1] + + # n-d arrays + heap = HeapAllocator.Heap() + A = HeapAllocator.alloc_array(heap, Int, 64, 64) + B = HeapAllocator.alloc_array(heap, Int, (64, 32)) + + # Error paths + heap = HeapAllocator.Heap() + @test_throws ErrorException("pointer not malloc'd in this heap") HeapAllocator.free(heap, HeapAllocator.SizedPtr{UInt8}(0xbadc0ffee0ddf00d, 8 % UInt)) + HeapAllocator.malloc(heap, 8) + @test_throws ErrorException("pointer not malloc'd in this heap") HeapAllocator.free(heap, HeapAllocator.SizedPtr{UInt8}(0xbadc0ffee0ddf00d, 8 % UInt)) + + # Smoke show + show(devnull, MIME"text/plain"(), heap) + HeapAllocator.heap_stats(heap) + +end diff --git a/test/blockarrays.jl b/test/blockarrays.jl index 59b1d8db9a..4f9da0db17 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]) + create_sparsity_pattern!(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) + create_sparsity_pattern!(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..417e9463f2 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("HeapAllocator.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 39b1c38d70..d7704a5025 100644 --- a/test/test_dofs.jl +++ b/test/test_dofs.jl @@ -180,12 +180,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) @@ -429,11 +429,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) + create_sparsity_pattern!(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 @@ -443,25 +449,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) + create_sparsity_pattern!(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 @@ -471,37 +482,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) + create_sparsity_pattern!(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)) @@ -512,28 +531,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) + create_sparsity_pattern!(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 @@ -615,18 +641,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)) @@ -636,16 +662,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. @@ -661,10 +687,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 @@ -674,9 +700,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 e0a7213217..f2ebf77938 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..3b02c41873 --- /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.n_rows(tp::TestPattern) = tp.nrowscols[1] +Ferrite.n_cols(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.n_rows(p1) == Ferrite.n_rows(p), px) + @test all(p -> Ferrite.n_cols(p1) == Ferrite.n_cols(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.n_rows(dsp) == m + @test Ferrite.n_cols(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) + create_sparsity_pattern!(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) + create_sparsity_pattern!(dsp, dh) + @test dsp isa SparsityPattern + A = allocate_matrix(dsp) + @test A isa SparseMatrixCSC{Float64, Int} + A + end, + let + dsp = init_sparsity_pattern(dh) + create_sparsity_pattern!(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(create_sparsity_pattern!(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) + create_sparsity_pattern!(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) + create_sparsity_pattern!(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(create_sparsity_pattern!(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) + create_sparsity_pattern!(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) + create_sparsity_pattern!(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(create_sparsity_pattern!(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 + create_sparsity_pattern!(p, dh) + end + compare_patterns(ps...) + + # DofHandler + ConstraintHandler + ps = make_patterns(dh) + for p in ps + create_sparsity_pattern!(p, dh, ch) + end + compare_patterns(ps...) + + # DofHandler + ConstraintHandler later + ps = make_patterns(dh) + for p in ps + create_sparsity_pattern!(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 + create_sparsity_pattern!(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 + create_sparsity_pattern!(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 + create_sparsity_pattern!(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 + create_sparsity_pattern!(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") create_sparsity_pattern!(p, dh_open) + end + for p in (SparsityPattern(ndofs(dh), 2), TestPattern(ndofs(dh), 2), BlockSparsityPattern([2, 2])) + @test_throws ErrorException create_sparsity_pattern!(p, dh) + end + for p in (SparsityPattern(2, ndofs(dh)), TestPattern(2, ndofs(dh))) + @test_throws ErrorException create_sparsity_pattern!(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 create_sparsity_pattern!(p, dh; keep_constrained=false) + end + ch_open = ConstraintHandler(dh) + for p in patterns + @test_throws ErrorException create_sparsity_pattern!(p, dh, ch_open; keep_constrained=false) + end + ch_bad = ConstraintHandler(close!(DofHandler(grid))) + for p in patterns + @test_throws ErrorException create_sparsity_pattern!(p, dh, ch_bad; keep_constrained=false) + end +end