diff --git a/CHANGELOG.md b/CHANGELOG.md index af1d16cba2..95a03f9732 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,13 +6,19 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] +### Added + - [Metis.jl](https://github.com/JuliaSparse/Metis.jl) extension for fill-reducing DoF + permutation. This uses Julias new package extension mechanism (requires Julia 1.10) to + support a new DoF renumbering order `DofOrder.Ext{Metis}()` that can be passed to + `renumber!` to renumber DoFs using the Metis.jl library. ([#393][github-393], + [#549][github-549]) ## [0.3.10] - 2022-12-11 ### Added - New functions `apply_local!` and `apply_assemble!` for applying constraints locally on the element level before assembling to the global system. ([#528][github-528]) - New functionality to renumber DoFs by fields or by components. This is useful when you - need the global matrix to be blocked. ([#545][github-545]) + need the global matrix to be blocked. ([#378][github-378], [#545][github-545]) - Functionality to renumber DoFs in DofHandler and ConstraintHandler simultaneously: `renumber!(dh::DofHandler, ch::ConstraintHandler, order)`. Previously renumbering had to be done *before* creating the ConstraintHandler since otherwise DoF numbers would be @@ -153,10 +159,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 [github-352]: https://github.com/Ferrite-FEM/Ferrite.jl/issues/352 [github-363]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/363 +[github-378]: https://github.com/Ferrite-FEM/Ferrite.jl/issues/378 [github-385]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/385 [github-386]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/386 [github-390]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/390 [github-392]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/392 +[github-393]: https://github.com/Ferrite-FEM/Ferrite.jl/issues/393 [github-401]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/401 [github-402]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/402 [github-404]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/404 @@ -215,6 +223,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 [github-544]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/544 [github-545]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/545 [github-547]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/547 +[github-549]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/549 [github-550]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/550 [Unreleased]: https://github.com/Ferrite-FEM/Ferrite.jl/compare/v0.3.10...HEAD diff --git a/Project.toml b/Project.toml index 751f8c79af..4bd79e1ead 100644 --- a/Project.toml +++ b/Project.toml @@ -12,8 +12,15 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" +[weakdeps] +Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b" + +[extensions] +MetisExt = "Metis" + [compat] EnumX = "1" +Metis = "1.3" NearestNeighbors = "0.4" Preferences = "1" Reexport = "1" @@ -28,6 +35,7 @@ FerriteGmsh = "4f95f4f8-b27c-4ae5-9a39-ea55e634e36b" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" +Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b" NBInclude = "0db19996-df87-5ea3-a455-e3a50d440464" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -37,4 +45,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [targets] -test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "NBInclude", "ProgressMeter", "Random", "SHA", "StableRNGs", "Test", "TimerOutputs"] +test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "NBInclude", "ProgressMeter", "Random", "SHA", "StableRNGs", "Test", "TimerOutputs"] diff --git a/ext/MetisExt.jl b/ext/MetisExt.jl new file mode 100644 index 0000000000..e43d7faa6c --- /dev/null +++ b/ext/MetisExt.jl @@ -0,0 +1,83 @@ +module MetisExt + +using Ferrite +using Ferrite: AbstractDofHandler +using Metis.LibMetis: idx_t +using Metis: Metis +using SparseArrays: sparse + +struct MetisOrder <: DofOrder.Ext{Metis} + coupling::Union{Matrix{Bool},Nothing} +end + +""" + DofOrder.Ext{Metis}(; coupling) + +Fill-reducing permutation order from [Metis.jl](https://github.com/JuliaSparse/Metis.jl). + +Since computing the permutation involves constructing the structural couplings between all +DoFs the field/component coupling can be provided; see [`create_sparsity_pattern`](@ref) for +details. +""" +function DofOrder.Ext{Metis}(; + coupling::Union{AbstractMatrix{Bool},Nothing}=nothing, +) + return MetisOrder(coupling) +end + +function Ferrite.compute_renumber_permutation( + dh::AbstractDofHandler, + ch::Union{ConstraintHandler,Nothing}, + order::DofOrder.Ext{Metis} +) + + # Expand the coupling matrix to size ndofs_per_cell × ndofs_per_cell + coupling = order.coupling + if coupling === nothing + n = ndofs_per_cell(dh) + entries_per_cell = n * (n - 1) + else # coupling !== nothing + # Set sym = true since Metis.permutation requires a symmetric graph. + # TODO: Perhaps just symmetrize it: coupling = coupling' .| coupling + coupling = Ferrite._coupling_to_local_dof_coupling(dh, coupling, #= sym =# true) + # Compute entries per cell, subtract diagonal elements + entries_per_cell = + count(coupling[i, j] for i in axes(coupling, 1), j in axes(coupling, 2) if i != j) + end + + # Create the CSR (CSC, but pattern is symmetric so equivalent) using + # Metis.idx_t as the integer type + L = entries_per_cell * getncells(dh.grid) + I = Vector{idx_t}(undef, L) + J = Vector{idx_t}(undef, L) + idx = 0 + @inbounds for cc in CellIterator(dh) + dofs = celldofs(cc) + for (i, dofi) in pairs(dofs), (j, dofj) in pairs(dofs) + dofi == dofj && continue # Metis doesn't want the diagonal + coupling === nothing || coupling[i, j] || continue + idx += 1 + I[idx] = dofi + J[idx] = dofj + end + end + @assert idx == L + N = ndofs(dh) + # TODO: Use spzeros! in Julia 1.10. + S = sparse(I, J, zeros(Float32, length(I)), N, N) + + # Add entries from affine constraints + if ch !== nothing + error("TODO: Use constraints.") + end + + # Construct a Metis.Graph + G = Metis.Graph(idx_t(N), S.colptr, S.rowval) + + # Compute the permutation + _, perm = Metis.permutation(G) + + return perm +end + +end # module MetisExt diff --git a/src/Dofs/DofRenumbering.jl b/src/Dofs/DofRenumbering.jl index 410dfc930c..25e0ba1356 100644 --- a/src/Dofs/DofRenumbering.jl +++ b/src/Dofs/DofRenumbering.jl @@ -25,6 +25,18 @@ module DofOrder ComponentWise(x=Int[]) = new(_check_target_blocks(x)) end + """ + DofOrder.Ext{T} + + DoF permutation order from external package `T`. Currently supported extensions: + - `DofOrder.Ext{Metis}`: Fill-reducing permutation from + [Metis.jl](https://github.com/JuliaSparse/Metis.jl). + """ + abstract type Ext{T} end + function Ext{T}(args...; kwargs...) where T + throw(ArgumentError("Unknown external order DofOrder.Ext{$T}. See documentation for `DofOrder.Ext` for details.")) + end + end # module DofOrder """ @@ -39,6 +51,8 @@ ordering `order`. `perm[i]`. - [`DofOrder.FieldWise()`](@ref) for renumbering dofs field wise. - [`DofOrder.ComponentWise()`](@ref) for renumbering dofs component wise. + - `DofOrder.Ext{T}` for "external" renumber permutations, see documentation for + `DofOrder.Ext` for details. !!! warning The dof numbering in the DofHandler and ConstraintHandler *must always be consistent*. @@ -204,3 +218,7 @@ function compute_renumber_permutation(dh::DofHandler, _, order::DofOrder.Compone perm = invperm(iperm) return perm end + +function compute_renumber_permutation(dh::AbstractDofHandler, ::Union{ConstraintHandler,Nothing}, ::DofOrder.Ext{M}) where M + error("Renumbering extension based on package $M not available.") +end diff --git a/test/runtests.jl b/test/runtests.jl index def659697e..c9a2e3eaf1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,11 @@ using Random using LinearAlgebra using SparseArrays +const HAS_EXTENSIONS = isdefined(Base, :get_extension) +if HAS_EXTENSIONS + import Metis +end + include("test_utils.jl") include("test_interpolations.jl") include("test_cellvalues.jl") diff --git a/test/test_dofs.jl b/test/test_dofs.jl index 0e005d3b44..c4204bdd0b 100644 --- a/test/test_dofs.jl +++ b/test/test_dofs.jl @@ -254,6 +254,16 @@ end for el in 1:2, r in [dof_range(dh, :v)[1:2:end], dof_range(dh, :v)[2:2:end], dof_range(dh, :s)] @test sign.(diff(celldofs(dh, el)[r])) == sign.(diff(celldofs(dho, el)[r])) end + + # Metis ordering + if HAS_EXTENSIONS + # TODO: Should probably test that the new order result in less fill-in + dh, ch = testdhch() + renumber!(dh, DofOrder.Ext{Metis}()) + @test_throws ErrorException renumber!(dh, ch, DofOrder.Ext{Metis}()) + renumber!(dh, DofOrder.Ext{Metis}(coupling=[true true; true false])) + @test_throws ErrorException renumber!(dh, ch, DofOrder.Ext{Metis}(coupling=[true true; true false])) + end end @testset "dof coupling" begin