Metis extension for renumbering.

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

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"]

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

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

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")

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