Add topic section on local and global dof numbering

docs/Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.11.0"
 manifest_format = "2.0"
-project_hash = "4e52f4aa4cee9f66ec4f633f0ae538fbd227ac5e"
+project_hash = "4c773ce49d51766622c20f611ef57f359d1127e8"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "eea5d80188827b35333801ef97a40c2ed653b081"
@@ -185,6 +185,11 @@ git-tree-sha1 = "9e2a6b69137e6969bab0152632dcb3bc108c8bdd"
 uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
 version = "1.0.8+1"
 
+[[deps.CEnum]]
+git-tree-sha1 = "389ad5c84de1ae7cf0e28e381131c98ea87d54fc"
+uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
+version = "0.5.0"
+
 [[deps.CPUSummary]]
 deps = ["CpuId", "IfElse", "PrecompileTools", "Static"]
 git-tree-sha1 = "5a97e67919535d6841172016c9530fd69494e5ec"
@@ -414,10 +419,10 @@ uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
 version = "1.15.1"
 
 [[deps.DifferentiationInterface]]
-deps = ["ADTypes", "Compat", "LinearAlgebra", "PackageExtensionCompat"]
-git-tree-sha1 = "75d1716eb46e1b77304f7d79ec9ac6a4e6185d02"
+deps = ["ADTypes", "LinearAlgebra"]
+git-tree-sha1 = "42c9752ea96e117edd384b1bf063aff4138271db"
 uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
-version = "0.6.7"
+version = "0.6.8"
 
     [deps.DifferentiationInterface.extensions]
     DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore"
@@ -578,18 +583,15 @@ version = "2.0.4"
 
 [[deps.Ferrite]]
 deps = ["EnumX", "ForwardDiff", "LinearAlgebra", "NearestNeighbors", "OrderedCollections", "Preferences", "Reexport", "SparseArrays", "StaticArrays", "Tensors", "WriteVTK"]
-path = ".."
+git-tree-sha1 = "1ad6e7eab1803998ad16360b1f86c7723d361cb0"
 uuid = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
 version = "1.0.0"
+weakdeps = ["BlockArrays", "Metis"]
 
     [deps.Ferrite.extensions]
     FerriteBlockArrays = "BlockArrays"
     FerriteMetis = "Metis"
 
-    [deps.Ferrite.weakdeps]
-    BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-    Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
-
 [[deps.FerriteGmsh]]
 deps = ["Ferrite", "Gmsh"]
 git-tree-sha1 = "9669f21d4ddc68ffca0d5ea12d3ac6b438b9af06"
@@ -1299,6 +1301,22 @@ git-tree-sha1 = "c13304c81eec1ed3af7fc20e75fb6b26092a1102"
 uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
 version = "0.3.2"
 
+[[deps.Metis]]
+deps = ["CEnum", "LinearAlgebra", "METIS_jll", "SparseArrays"]
+git-tree-sha1 = "54aca4fd53d39dcd2c3f1bef367b6921e8178628"
+uuid = "2679e427-3c69-5b7f-982b-ece356f1e94b"
+version = "1.5.0"
+
+    [deps.Metis.extensions]
+    MetisGraphs = "Graphs"
+    MetisLightGraphs = "LightGraphs"
+    MetisSimpleWeightedGraphs = ["SimpleWeightedGraphs", "Graphs"]
+
+    [deps.Metis.weakdeps]
+    Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
+    LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
+    SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622"
+
 [[deps.MicrosoftMPI_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "f12a29c4400ba812841c6ace3f4efbb6dbb3ba01"
@@ -2081,10 +2099,14 @@ version = "2.22.0"
     Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [[deps.SparseMatrixColorings]]
-deps = ["ADTypes", "Compat", "DataStructures", "DocStringExtensions", "LinearAlgebra", "Random", "SparseArrays"]
-git-tree-sha1 = "ad17e2069015839e58a1f9275608b405fae1f28e"
+deps = ["ADTypes", "DataStructures", "DocStringExtensions", "LinearAlgebra", "Random", "SparseArrays"]
+git-tree-sha1 = "ccc32032d8f6790ad198c99fb8ef20d8100a0de4"
 uuid = "0a514795-09f3-496d-8182-132a7b665d35"
-version = "0.4.6"
+version = "0.4.7"
+weakdeps = ["Colors"]
+
+    [deps.SparseMatrixColorings.extensions]
+    SparseMatrixColoringsColorsExt = "Colors"
 
 [[deps.Sparspak]]
 deps = ["Libdl", "LinearAlgebra", "Logging", "OffsetArrays", "Printf", "SparseArrays", "Test"]

docs/Project.toml

@@ -13,6 +13,7 @@ LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
+Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
 OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"

docs/src/topics/degrees_of_freedom.md

@@ -21,14 +21,17 @@ dh = DofHandler(grid)
 
 Before we can distribute the dofs we need to specify fields. A field is simply the unknown
 function(s) we are solving for. To add a field we need a name (a `Symbol`) and the the
-interpolation describing the shape functions for the field. Here we add a scalar field `:p`,
-interpolated using linear (degree 1) shape functions on a triangle, and a vector field `:u`,
+interpolation describing the shape functions for the field. Here we add a scalar field `:u`,
+interpolated using linear (degree 1) shape functions on a triangle, and a vector field `:v`,
 also interpolated with linear shape functions on a triangle, but raised to the power 2 to
 indicate that it is a vector field with 2 components (for a 2D problem).
 
 ```@example dofs
-add!(dh, :p, Lagrange{RefTriangle, 1}())
-add!(dh, :u, Lagrange{RefTriangle, 1}()^2)
+interpolation_u = Lagrange{RefTriangle, 1}()
+interpolation_v = Lagrange{RefTriangle, 1}() ^ 2
+
+add!(dh, :u, interpolation_u)
+add!(dh, :v, interpolation_v)
 # hide
 ```
 
@@ -40,9 +43,188 @@ dofs for the fields we added.
 close!(dh)
 ```
 
-## Ordering of Dofs
 
-!!! todo
-    Describe dof ordering within elements (vertices -> edges -> faces ->
-    volumes) and `dof_range`.
-    Describe (global) dof renumbering
+## Local DoF indices
+
+Locally on each element the DoFs are ordered by field, in the same order they are were added
+to the DofHandler. Within each field the DoFs follow the order of the interpolation.
+Concretely this means that the local matrix is a block matrix (and the local vector a block
+vector).
+
+For the example DofHandler defined above, with the two fields `:u` and `:v`, the local
+system would look something like:
+
+```math
+\begin{bmatrix}
+    K^e_{uu} & K^e_{uv} \\
+    K^e_{vu} & K^e_{vv}
+\end{bmatrix}
+\begin{bmatrix}
+    u^e \\
+    v^e
+\end{bmatrix} =
+\begin{bmatrix}
+    f^e_u \\
+    f^e_v
+\end{bmatrix}
+```
+
+where
+
+ - ``K^e_{uu}`` couples test and trial functions for the `:u` field
+ - ``K^e_{uv}`` couples test functions for the `:u` field with trial functions for the `:v` field
+ - ``K^e_{vu}`` couples test functions for the `:v` field with trial functions for the `:u` field
+ - ``K^e_{vv}`` couples test and trial functions for the `:v` field
+
+We can query the local index ranges for the dofs of the two fields with the
+[`dof_range`](@ref) function:
+
+```@example dofs
+u_range = dof_range(dh, :u)
+v_range = dof_range(dh, :v)
+
+u_range, v_range
+```
+
+i.e. the local indices for the `:u` field are `1:3` and for the `:v` field `4:9`. This
+matches directly with the number of dofs for each field: 3 for `:u` and 6 for `:v`. The
+ranges are used when assembling the blocks of the matrix, i.e. if `Ke` is the local matrix,
+then `Ke[u_range, u_range]` corresponds to the``K_{uu}``, `Ke[u_range, v_range]` to
+``K_{uv}``, etc. See for example [Tutorial 8: Stokes flow](../tutorials/stokes-flow.md) for
+how the ranges are used.
+
+
+## Global DoF indices
+
+The global ordering of dofs, however, does *not* follow any specific order by default. In
+particular, they are *not* ordered by field, so while the local system is a block system,
+the global system is not. For all intents and purposes, the default global dof ordering
+should be considered an implementation detail.
+
+!!! warning "DoF numbering is decoupled from the node numbering"
+    A common pitfall for new Ferrite is to assume that the numbering of DoFs follows the
+    numbering of the nodes in the grid. This is *not* the case. While it would be possible
+    to align these two numberings in some special cases (single isoparametric scalar field)
+    it is not possible in general.
+
+We can query the global dofs for a cell with the `celldofs` function. For example, for cell
+#45 the global dofs are:
+
+```@example dofs
+global_dofs = celldofs(dh, 45)
+```
+
+which looks more or less random. We can also look at the mapping between local and global
+dofs for field `:u`:
+
+```@example dofs
+u_range .=> global_dofs[u_range]
+```
+
+and for field `:v`:
+
+```@example dofs
+v_range .=> global_dofs[v_range]
+```
+
+which makes it clear that `:u`-dofs and `:v`-dofs are interleaved in the global numbering.
+
+
+## Renumbering global DoF indices
+
+The global DoF order determines the sparsity pattern of the global matrix. In some cases,
+mostly depending on the linear solve strategy, it can be beneficial to reorder the global
+DoFs. For example:
+
+ - For a direct solver the sparsity pattern determines the
+   [fill-in](https://en.wikipedia.org/wiki/Sparse_matrix#Reducing_fill-in). An optimized
+   order can reduce the fill-in and thus the memory consumption and the computational cost.
+   (For a iterative solver the order doesn't matter as much since the order doesn't affect
+   number of operations in a matrix-vector product.)
+ - Using block-based solvers is possible for multi-field problems. In this case it is
+   simpler if the global system is also ordered by field and/or components similarly to the
+   local system.
+
+It is important to note that renumbering the global dofs does *not* affect the local order.
+The local system is *always* blocked by fields as described above.
+
+The default ordering (which, again, should be considered a black box) results in the
+following sparsity pattern:
+
+```@example dofs
+allocate_matrix(dh)
+nothing # hide
+```
+
+!!! details "Default sparsity pattern"
+    ```@example dofs
+    allocate_matrix(dh) # hide
+    ```
+
+This looks pretty good (all entries are concentrated around the diagonal) but it is
+important to note that this is just a happy coincidence because i) we used the built-in grid
+generator (which numbers neighboring cells consecutively) and because ii) Ferrite by default
+distributes global dofs cell by cell.
+
+### Renumbering for a global block system
+
+In order to obtain a global block system it is possible to renumber by fields, or even by
+component, using the [`renumber!`](@ref) function with the [`DofOrder.FieldWise`](@ref) and
+[`DofOrder.ComponentWise`](@ref) orders, respectively.
+
+For example, renumbering by fields gives the global block system
+```math
+\begin{bmatrix}
+    K_{uu} & K_{uv} \\
+    K_{vu} & K_{vv}
+\end{bmatrix}
+```
+
+```@example dofs
+renumber!(dh, DofOrder.FieldWise())
+allocate_matrix(dh)
+nothing # hide
+```
+
+!!! details "Sparsity pattern after field-wise renumbering"
+    ```@example dofs
+    allocate_matrix(dh) # hide
+    ```
+
+Similarly, global blocking can be done by components, and fields and/or component can be
+permuted in the global system. See[`DofOrder.FieldWise`](@ref) and
+[`DofOrder.ComponentWise`](@ref) for more details.
+
+### Renumbering to reduce fill-in
+
+Ferrite can be used together with the [Metis.jl](https://github.com/JuliaSparse/Metis.jl)
+package to optimize the global DoF ordering w.r.t. reduced fill-in:
+
+```@example dofs
+using Metis
+renumber!(dh, DofOrder.Ext{Metis}())
+allocate_matrix(dh)
+nothing # hide
+```
+
+!!! details "Sparsity pattern after Metis renumbering"
+    ```@example dofs
+    allocate_matrix(dh) # hide
+    ```
+
+The [`DofOrder.FieldWise`](@ref) and [`DofOrder.ComponentWise`](@ref) orders preserve
+internal ordering within the fields/components so they can be combined with Metis to reduce
+fill-in for the individual blocks, for example:
+
+```@example dofs
+using Metis
+renumber!(dh, DofOrder.Ext{Metis}())
+renumber!(dh, DofOrder.FieldWise())
+allocate_matrix(dh)
+nothing # hide
+```
+
+!!! details "Sparsity pattern after Metis and field-wise renumbering"
+    ```@example dofs
+    allocate_matrix(dh) # hide
+    ```