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`.
Ferrite-FEM · Jun 20, 2024 · 7b82c7f · 7b82c7f
1 parent 54ed4d9
commit 7b82c7f
Show file tree

Hide file tree

Showing 49 changed files with 1,941 additions and 403 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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...)
+```