Skip to content

Commit

Permalink
Metis package extension for fill-reducing renumbering
Browse files Browse the repository at this point in the history
This patch adds new DoF renumbering possibilities using "package
extensions" which is available from Julia 1.10. These extensions hook
into the new `DofOrder.Ext{T}` type, where `T` is the external package
module providing the new ordering capability.

Currently only the fill-reducing renumbering algorithm from Metis.jl is
provided, ie.. `DofOrder.Ext{Metis}()`.

Closes #393.
  • Loading branch information
fredrikekre committed Dec 11, 2022
1 parent 9ae40b3 commit a49d32b
Show file tree
Hide file tree
Showing 6 changed files with 136 additions and 2 deletions.
11 changes: 10 additions & 1 deletion CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
10 changes: 9 additions & 1 deletion Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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"
Expand All @@ -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"]
83 changes: 83 additions & 0 deletions ext/MetisExt.jl
Original file line number Diff line number Diff line change
@@ -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
18 changes: 18 additions & 0 deletions src/Dofs/DofRenumbering.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

"""
Expand All @@ -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*.
Expand Down Expand Up @@ -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
5 changes: 5 additions & 0 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
11 changes: 11 additions & 0 deletions test/test_dofs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -254,6 +254,17 @@ 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
@test false
# 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
Expand Down

0 comments on commit a49d32b

Please sign in to comment.