From 8ef4f726d59bce1bebbd517f560c84f0c903abc8 Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Sat, 5 Oct 2024 13:52:17 +0200 Subject: [PATCH] Add topic section on local and global dof numbering --- docs/src/topics/degrees_of_freedom.md | 200 ++++++++++++++++++++++++-- 1 file changed, 191 insertions(+), 9 deletions(-) diff --git a/docs/src/topics/degrees_of_freedom.md b/docs/src/topics/degrees_of_freedom.md index 144b792570..e1ad6de738 100644 --- a/docs/src/topics/degrees_of_freedom.md +++ b/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 + ```