diff --git a/docs/src/literate-howto/postprocessing.jl b/docs/src/literate-howto/postprocessing.jl
index 2eeb09b3c3..0bb4eeb912 100644
--- a/docs/src/literate-howto/postprocessing.jl
+++ b/docs/src/literate-howto/postprocessing.jl
@@ -13,12 +13,28 @@
 #
 # ## Introduction
 #
-# After running a simulation, we usually want to visualize the results in different ways.
-# The `L2Projector` and the `PointEvalHandler` build a pipeline for doing so. With the `L2Projector`,
-# integration point quantities can be projected to the nodes. The `PointEvalHandler` enables evaluation of
-# the finite element approximated function in any coordinate in the domain. Thus with the combination of both functionalities,
-# both nodal quantities and integration point quantities can be evaluated in any coordinate, allowing for example
-# cut-planes through 3D structures or cut-lines through 2D-structures.
+# After running a simulation, we usually want to postprocess and visualize the results in
+# different ways. Ferrite provides several tools to facilitate these tasks:
+#
+#  - L2 projection of (discrete) data onto FE interpolations using the `L2Projector`
+#  - Evalutation of fields (solutions, projections, etc) at arbitrary, user-defined, points
+#    using the `PointEvalHandler`
+#  - Builtin functionality for exporting data (solutions, cell data, projections, etc) to
+#    the VTK format
+#  - [Makie.jl](https://docs.makie.org/) based plotting using
+#    [FerriteViz.jl](https://ferrite-fem.github.io/FerriteViz.jl/)
+#
+# This how-to demonstrates the VTK export, the `L2Projector` for projecting discrete
+# quadrature point data onto a continuous FE interpolation, and the `PointEvalHandler` for
+# evaluating the FE solution, and the projection, along a user-defined cut line through the
+# domain.
+
+# !!! warning "Custom visualization"
+#     A common assumption is that the numbering of degrees of freedom matche the numbering
+#     of the nodes in the grid. This is *NOT* the case in Ferrite. If the available tools
+#     don't suit your needs and you decide to "roll your own" visualization you need to be
+#     aware of this and take it into account. For the specific case of evaluating the
+#     solution at the grid nodes you can use [`evaluate_at_grid_nodes`](@ref).
 #
 # This example continues from the Heat equation example, where the temperature field was
 # determined on a square domain. In this example, we first compute the heat flux in each
diff --git a/docs/src/literate-tutorials/heat_equation.jl b/docs/src/literate-tutorials/heat_equation.jl
index 9c190d4ee7..46316464eb 100644
--- a/docs/src/literate-tutorials/heat_equation.jl
+++ b/docs/src/literate-tutorials/heat_equation.jl
@@ -70,6 +70,11 @@ dh = DofHandler(grid)
 add!(dh, :u, ip)
 close!(dh);
 
+# !!! warning "Numbering of degrees of freedom"
+#     A common assumption is that the numbering of degrees of freedom follows the global
+#     numbering of the nodes in the grid. This is *NOT* the case in Ferrite. For more
+#     details, see the [Ferrite numbering rules](@ref "Ordering-of-Dofs").
+
 # Now that we have distributed all our dofs we can create our tangent matrix,
 # using `allocate_matrix`. This function returns a sparse matrix
 # with the correct entries stored.
@@ -210,6 +215,11 @@ K, f = assemble_global(cellvalues, K, dh);
 apply!(K, f, ch)
 u = K \ f;
 
+# !!! warning "Numbering of degrees of freedom"
+#     Once again, recall that numbering of degrees of freedom does *NOT* follow the global
+#     numbering of the nodes in the grid. Specifically, `u[i]` is *NOT* the temperature at
+#     node `i`.
+
 # ### Exporting to VTK
 # To visualize the result we export the grid and our field `u`
 # to a VTK-file, which can be viewed in e.g. [ParaView](https://www.paraview.org/).
diff --git a/docs/src/literate-tutorials/linear_elasticity.jl b/docs/src/literate-tutorials/linear_elasticity.jl
index c067e5d212..16b018fa08 100644
--- a/docs/src/literate-tutorials/linear_elasticity.jl
+++ b/docs/src/literate-tutorials/linear_elasticity.jl
@@ -150,6 +150,11 @@ dh = DofHandler(grid)
 add!(dh, :u, ip)
 close!(dh);
 
+# !!! warning "Numbering of degrees of freedom"
+#     A common assumption is that the numbering of degrees of freedom follows the global
+#     numbering of the nodes in the grid. This is *NOT* the case in Ferrite. For more
+#     details, see the [Ferrite numbering rules](@ref "Ordering-of-Dofs").
+
 # ### Boundary conditions
 # We set Dirichlet boundary conditions by fixing the motion normal to the bottom and left
 # boundaries. The last argument to `Dirichlet` determines which components of the field should be
@@ -315,6 +320,11 @@ assemble_external_forces!(f_ext, dh, getfacetset(grid, "top"), facetvalues, trac
 apply!(K, f_ext, ch)
 u = K \ f_ext;
 
+# !!! warning "Numbering of degrees of freedom"
+#     Once again, recall that numbering of degrees of freedom does *NOT* follow the global
+#     numbering of the nodes in the grid. Specifically, `u[2 * i - 1]` and `u[2 * i]` are
+#     *NOT* the displacements at node `i`.
+
 # ### Postprocessing
 # In this case, we want to analyze the displacements, as well as the stress field.
 # We calculate the stress in each quadrature point, and then export it in two different
diff --git a/docs/src/topics/fe_intro.md b/docs/src/topics/fe_intro.md
index 957c8b21fa..85ae6a484c 100644
--- a/docs/src/topics/fe_intro.md
+++ b/docs/src/topics/fe_intro.md
@@ -72,13 +72,18 @@ set of *elements* or *cells*. We call this geometric discretization *grid* (or *
 and denote it with $\Omega_h$. In this example the corners of the triangles are called
 *nodes*.
 
-Next we introduce the finite element approximation $u_\mathrm{h} \approx u$ as a sum of N nodal
-*shape functions*, where we denote each of these function by $\phi_i$ and the corresponding
-*nodal values* $\hat{u}_i$. Note that *shape functions* are sometimes referred to as
-*basis functions* or *trial functions*, and instead of $\phi_i$ they are sometimes denoted $N_i$.
-In this example we choose to approximate the test function in the same way. This approach is known
-as the *Galerkin finite element method*. Formally we write the evaluation of our approximations
-at a specific point $\mathbf{x}$ in our domain $\Omega$ as:
+Next, we introduce the finite element approximation $u_\mathrm{h} \approx u$ as linear combination of
+*shape functions*, $\phi_i$. The corresponding weights are usually called *degree of freedoms* (dofs), $\hat{u}_i$.
+Sometimes, the dofs are called *weights* or *nodal values*. In Ferrite, the numbering of the dofs does not correspond
+to the node numbers in the grid. While such numbering is common in basic finite element codes,
+Ferrite supports different approximations of the finite element fields and the geometry, prohibiting
+such basic numbering. For more details, see the [Ferrite numbering rules](@ref "Ordering-of-Dofs").
+
+Note that *shape functions* are sometimes referred to as *basis functions* or *trial functions*,
+and instead of $\phi_i$ they are sometimes denoted $N_i$. In this example we choose to approximate
+the test function in the same way. This approach is known as the *Bubnov-Galerkin finite element
+method*. Formally we write the evaluation of our approximations at a specific point $\mathbf{x}$
+in our domain $\Omega$ as:
 
 ```math
 u_\mathrm{h}(\mathbf{x}) = \sum_{i=1}^{\mathrm{N}} \phi_i(\mathbf{x}) \, \hat{u}_i,\qquad
diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl
index 84f59e62b7..6ada38eb90 100644
--- a/src/Dofs/DofHandler.jl
+++ b/src/Dofs/DofHandler.jl
@@ -129,6 +129,10 @@ add!(dh, :u, ip_u)
 add!(dh, :p, ip_p)
 close!(dh)
 ```
+
+!!! note "Numbering of degree of freedom"
+    Note that the numbering of degrees of freedom do *NOT* follow the global numbering of
+    nodes in the associated grid.
 """
 function DofHandler(grid::G) where {dim, G <: AbstractGrid{dim}}
     ncells = getncells(grid)