diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup.jl
new file mode 100644
index 0000000000..c93660e9bc
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup.jl
@@ -0,0 +1,253 @@
+#src # Behind the scenes of a simulation setup
+
+# This tutorial will guide you through a simple Trixi.jl setup ("elixir"), giving an overview of
+# what happens in the background during the initialization of a simulation. While the setup
+# described herein does not cover all details, it involves relatively stable parts of Trixi.jl that
+# are unlikely to undergo significant changes in the near future. The goal is to clarify some of
+# the more fundamental, *technical* concepts that are applicable to a variety of
+# (also more complex) configurations.
+
+# Trixi.jl follows the [method of lines](http://www.scholarpedia.org/article/Method_of_lines) concept for solving partial differential equations (PDEs).
+# Firstly, the PDEs are reduced to a (potentially huge) system of
+# ordinary differential equations (ODEs) by discretizing the spatial derivatives. Subsequently,
+# these generated ODEs may be solved with methods available in OrdinaryDiffEq.jl or those specifically
+# implemented in Trixi.jl. The following steps elucidate the process of transitioning from PDEs to
+# ODEs within the framework of Trixi.jl.
+
+# ## Basic setup
+
+# Import essential libraries and specify an equation.
+
+using Trixi, OrdinaryDiffEq
+equations = LinearScalarAdvectionEquation2D((-0.2, 0.7))
+
+# Generate a spatial discretization using a [`TreeMesh`](@ref) with a pre-coarsened set of cells.
+
+coordinates_min = (-2.0, -2.0)
+coordinates_max = (2.0, 2.0)
+
+coarsening_patches = ((type = "box", coordinates_min = [0.0, -2.0],
+                       coordinates_max = [2.0, 0.0]),)
+
+mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 2,
+                n_cells_max = 30_000,
+                coarsening_patches = coarsening_patches)
+
+# The created `TreeMesh` looks like the following:
+
+# ![TreeMesh_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/d5ef76ee-8246-4730-a692-b472c06063a3)
+
+# Instantiate a [`DGSEM`](@ref) solver with a user-specified polynomial degree. The solver
+# will define `polydeg + 1` [Gauss-Lobatto nodes](https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss%E2%80%93Lobatto_rules) and their associated weights within
+# the reference interval ``[-1, 1]`` in each spatial direction. These nodes will be subsequently
+# used to approximate solutions on each leaf cell of the `TreeMesh`.
+
+solver = DGSEM(polydeg = 3)
+
+# Gauss-Lobatto nodes with `polydeg = 3`:
+
+# ![Gauss-Lobatto_nodes_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/1d894611-801e-4f75-bff0-d77ca1c672e5)
+
+# ## Overview of the [`SemidiscretizationHyperbolic`](@ref) type
+
+# At this stage, all necessary components for configuring the spatial discretization are in place.
+# The remaining task is to combine these components into a single structure that will be used
+# throughout the entire simulation process. This is where [`SemidiscretizationHyperbolic`](@ref)
+# comes into play.
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+# The constructor for the `SemidiscretizationHyperbolic` object calls numerous sub-functions to
+# perform the necessary initialization steps. A brief description of the key sub-functions is
+# provided below.
+
+
+# - `init_elements(leaf_cell_ids, mesh, equations, dg.basis, RealT, uEltype)`
+
+#   The fundamental elements for approximating the solution are the leaf
+#   cells. The solution is constructed as a polynomial of the degree specified in the `DGSEM`
+#   solver in each spatial direction on each leaf cell. This polynomial approximation is evaluated 
+#   at the Gauss-Lobatto nodes mentioned earlier. The `init_elements` function extracts
+#   these leaf cells from the `TreeMesh`, assigns them the label "elements", records their
+#   coordinates, and maps the Gauss-Lobatto nodes from the 1D interval ``[-1, 1]`` onto each coordinate axis
+#   of every element.
+
+
+#   ![elements_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/9f486670-b579-4e42-8697-439540c8bbb4)
+
+#   The visualization of elements with nodes shown here includes spaces between elements, which do
+#   not exist in reality. This spacing is included only for illustrative purposes to underscore the
+#   separation of elements and the independent projection of nodes onto each element.
+
+
+# - `init_interfaces(leaf_cell_ids, mesh, elements)`
+
+#   At this point, the elements with nodes have been defined; however, they lack the necessary
+#   communication functionality. This is crucial because the local solution polynomials on the
+#   elements are not independent of each other. Furthermore, nodes on the boundary of adjacent
+#   elements share the same spatial location, which requires a method to combine this into a
+#   meaningful solution.
+#   Here [Riemann solvers](https://en.wikipedia.org/wiki/Riemann_solver#Approximate_solvers)
+#   come into play which can handle the principal ambiguity of a multi-valued solution at the
+#   same spatial location.
+
+#   As demonstrated earlier, the elements can have varying sizes. Let us initially consider
+#   neighbors with equal size. For these elements, the `init_interfaces` function generates
+#   interfaces that store information about adjacent elements, their relative positions, and
+#   allocate containers for sharing solution data between neighbors during the solution process. 
+
+#   In our visualization, these interfaces would conceptually resemble tubes connecting the
+#   corresponding elements.
+
+#   ![interfaces_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/bc3b6b02-afbc-4371-aaf7-c7bdc5a6c540)
+
+
+# - `init_mortars(leaf_cell_ids, mesh, elements, dg.mortar)`
+
+#   Returning to the consideration of different sizes among adjacent elements, within the
+#   `TreeMesh`, adjacent leaf cells can vary in side length by a maximum factor of two. This
+#   implies that a large element has one neighbor of
+#   equal size with a connection through an interface, or two neighbors at half the size,
+#   requiring a connection through so called "mortars". In 3D, a large element would have
+#   four small neighbor elements.
+
+#   Mortars store information about the connected elements, their relative positions, and allocate
+#   containers for storing the solutions along the boundaries between these elements.
+
+#   Due to the differing sizes of adjacent elements, it is not feasible to directly map boundary
+#   nodes of adjacent elements. Therefore, the concept of mortars employs a mass-conserving
+#   interpolation function to map boundary nodes from a larger element to a smaller one.
+
+#   In our visualization, mortars are represented as branched tubes.
+
+#   ![mortars_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/43a95a60-3a31-4b1f-8724-14049e7a0481)
+
+
+# - `init_boundaries(leaf_cell_ids, mesh, elements)`
+
+#   In order to apply boundary conditions, it is necessary to identify the locations of the
+#   boundaries. Therefore, we initialize a "boundaries" object, which records the elements that
+#   contain boundaries, specifies which side of an element is a boundary, stores the coordinates
+#   of boundary nodes, and allocates containers for managing solutions at these boundaries.
+
+#   In our visualization, boundaries and their corresponding nodes are highlighted with green,
+#   semi-transparent lines.
+
+#   ![boundaries_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/21996b20-4a22-4dfb-b16a-e2c22c2f29fe)
+
+# All the structures mentioned earlier are collected as a cache of type `NamedTuple`. Subsequently,
+# an object of type `SemidiscretizationHyperbolic` is initialized using this cache, initial and
+# boundary conditions, equations, mesh and solver.
+
+# In conclusion, the primary purpose of a `SemidiscretizationHyperbolic` is to collect equations,
+# the geometric representation of the domain, and approximation instructions, creating specialized
+# structures to interconnect these components in a manner that enables their utilization for
+# the numerical solution of partial differential equations (PDEs).
+
+# As evident from the earlier description of `SemidiscretizationHyperbolic`, it comprises numerous
+# functions called subsequently. Without delving into details, the structure of the primary calls
+# are illustrated as follows:
+
+# ![SemidiscretizationHyperbolic_structure](https://github.com/trixi-framework/Trixi.jl/assets/119304909/8bf59422-0537-4d7a-9f13-d9b2253c19d7)
+
+# ## Overview of the [`semidiscretize`](@ref) function
+
+# At this stage, we have defined the equations and configured the domain's discretization. The
+# final step before solving is to select a suitable time span and apply the corresponding initial
+# conditions, which are already stored in the initialized `SemidiscretizationHyperbolic` object.
+
+# The purpose of the [`semidiscretize`](@ref) function is to wrap the semidiscretization as an
+# `ODEProblem` within the specified time interval. During this procedure the approximate solution
+#  is created at the given initial time via the specified `initial_condition` function from the 
+#  `SemidiscretizationHyperbolic` object. This `ODEProblem` can be subsequently passed to the
+# `solve` function from the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) package
+# or to [`Trixi.solve`](@ref).
+
+ode = semidiscretize(semi, (0.0, 1.0));
+
+# The `semidiscretize` function involves a deep tree of subsequent calls, with the primary ones
+# explained below.
+
+
+# - `allocate_coefficients(mesh, equations, solver, cache)`
+
+#   To apply initial conditions, a data structure ("container") needs to be generated to store the
+#   initial values of the target variables for each node within each element.
+
+#   Since only one-dimensional `Array`s are `resize!`able in Julia, we use `Vector`s as an internal
+#   storage for the target variables and `resize!` them whenever needed, e.g. to change the number
+#   of elements. Then, during the solving process the same memory is reused by `unsafe_wrap`ping
+#   multi-dimensional `Array`s around the internal storage.
+
+# - `wrap_array(u_ode, semi)`
+
+#   As previously noted, `u_ode` is constructed as a 1D vector to ensure compatibility with
+#   OrdinaryDiffEq.jl. However, for internal use within Trixi.jl, identifying which part of the
+#   vector relates to specific variables, elements, or nodes can be challenging.
+
+#   This is why the `u_ode` vector is wrapped by the `wrap_array` function using `unsafe_wrap`
+#   to form a multidimensional array `u`. In this array, the first dimension corresponds to
+#   variables, followed by N dimensions corresponding to nodes for each of N space dimensions.
+#   The last dimension corresponds to the elements. 
+#   Consequently, navigation within this multidimensional array becomes noticeably easier.
+
+#   "Wrapping" in this context involves the creation of a reference to the same storage location
+#   but with an alternative structural representation. This approach enables the use of both
+#   instances `u` and `u_ode` as needed, so that changes are simultaneously reflected in both.
+#   This is possible because, from a storage perspective, they share the same stored data, while
+#   access to this data is provided in different ways.
+
+
+# - `compute_coefficients!(u, initial_conditions, t, mesh::DG, equations, solver, cache)`
+
+#   Now the variable `u`, intended to store solutions, has been allocated and wrapped, it is time
+#   to apply the initial conditions. The `compute_coefficients!` function calculates the initial
+#   conditions for each variable at every node within each element and properly stores them in the
+#   `u` array.
+
+# At this stage, the `semidiscretize` function has all the necessary components to initialize and
+# return an `ODEProblem` object, which will be used by the `solve` function to compute the
+# solution.
+
+# In summary, the internal workings of `semidiscretize` with brief descriptions can be presented
+# as follows.
+
+# ![semidiscretize_structure](https://github.com/trixi-framework/Trixi.jl/assets/119304909/491eddc4-aadb-4e29-8c76-a7c821d0674e)
+
+# ## Functions `solve` and `rhs!`
+
+# Once the `ODEProblem` object is initialized, the `solve` function and one of the ODE solvers from
+# the OrdinaryDiffEq.jl package can be utilized to compute an approximated solution using the
+# instructions contained in the `ODEProblem` object.
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), dt = 0.01,
+            save_everystep = false);
+
+# Since the `solve` function and the ODE solver have no knowledge
+# of a particular spatial discretization, it is necessary to define a
+# "right-hand-side function", `rhs!`, within Trixi.jl.
+
+# Trixi.jl includes a set of `rhs!` functions designed to compute `du`, i.e.,
+# ``\frac{\partial u}{\partial t}`` according to the structure
+# of the setup. These `rhs!` functions calculate interface, mortars, and boundary fluxes, in
+# addition to surface and volume integrals, in order to construct the `du` vector. This `du` vector
+# is then used by the time integration method to obtain the solution at the subsequent time step.
+# The `rhs!` function is called by time integration methods in each iteration of the solve loop
+# within OrdinaryDiffEq.jl, with arguments `du`, `u`, `semidiscretization`, and the current time.
+
+# Trixi.jl uses a two-levels approach for `rhs!` functions. The first level is limited to a
+# single function for each `semidiscretization` type, and its role is to redirect data to the
+# target `rhs!` for specific solver and mesh types. This target `rhs!` function is responsible
+# for calculating `du`.
+
+# Path from the `solve` function call to the appropriate `rhs!` function call:
+
+# ![rhs_structure](https://github.com/trixi-framework/Trixi.jl/assets/119304909/dbea9a0e-25a4-4afa-855e-01f1ad619982)
+
+# Computed solution:
+
+using Plots
+plot(sol)
+pd = PlotData2D(sol)
+plot!(getmesh(pd))
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/Project.toml b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/Project.toml
new file mode 100644
index 0000000000..43aec5b7f5
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/Project.toml
@@ -0,0 +1,2 @@
+[deps]
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/README.md b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/README.md
new file mode 100644
index 0000000000..011b5c7586
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/README.md
@@ -0,0 +1,15 @@
+# Plots for the tutorial "Behind the scenes of a simulation setup"
+
+To create all the images for the tutorial, execute the following command from the directory of this `README.md`:
+```julia
+pkg> activate .
+julia> include.(readdir("src"; join=true))
+```
+To create all images from a different directory, substitute `"src"` with the path to the `src` 
+folder. The resulting images will be generated in your current directory as PNG files.
+
+To generate a specific image, run the following command while replacing `"path/to/src"` and `"file_name"` with the appropriate values:
+```julia
+pkg> activate .
+julia> include(joinpath("path/to/src", "file_name"))
+```
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/SemidiscretizationHyperbolic_structure_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/SemidiscretizationHyperbolic_structure_figure.jl
new file mode 100644
index 0000000000..cae7b19d47
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/SemidiscretizationHyperbolic_structure_figure.jl
@@ -0,0 +1,64 @@
+using Plots
+plot(Shape([(-2.3,4.5), (2.35,4.5), (2.35,2.5), (-2.3,2.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2, size=(800,600), showaxis=false, grid=false, xlim=(-2.4,2.8), ylim=(-25,5.5))
+annotate!(2.3, 3.5, ("SemidiscretizationHyperbolic(mesh, equations, initial_conditions, solver; source_terms,
+boundary_conditions, RealT, uEltype, initial_cache)          ", 10, :black, :right))
+annotate!(-2.3, 1.5, ("creates and returns SemidiscretizationHyperbolic object, initialized using a mesh, equations,
+initial_conditions, boundary_conditions, source_terms, solver and cache", 9, :black, :left))
+plot!([-1.2,-1.2],[0.6,-2],arrow=true,color=:black,linewidth=2,label="")
+plot!([-1.2,-1.4],[0.6,-2],arrow=true,color=:black,linewidth=2,label="")
+plot!([-1.2,-1.],[0.6,-2],arrow=true,color=:black,linewidth=2,label="")
+annotate!(-1, -0.7, ("specialized for mesh
+and solver types", 9, :black, :left))
+plot!([1.25,1.25],[0.6,-2],arrow=true,color=:black,linewidth=2,label="")
+plot!([1.25,1.05],[0.6,-2],arrow=true,color=:black,linewidth=2,label="")
+plot!([1.25,1.45],[0.6,-2],arrow=true,color=:black,linewidth=2,label="")
+annotate!(1.48, -0.7, ("specialized for mesh
+and boundary_conditions
+types", 9, :black, :left))
+
+plot!(Shape([(-2.3,-2), (-0.1,-2), (-0.1,-4), (-2.3,-4)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.2, -3, ("create_cache(mesh::TreeMesh, equations,
+                   solver::Dg, RealT, uEltype)", 10, :black, :center))
+plot!([-2.22,-2.22],[-4,-22],arrow=false,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-0.05,-2), (2.6,-2), (2.6,-4), (-0.05,-4)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(1.27, -3, ("digest_boundary_conditions(boundary_conditions,
+                                             mesh, solver, cache)", 10, :black, :center))
+annotate!(2.6, -5, ("if necessary, converts passed boundary_conditions
+ into a suitable form for processing by Trixi.jl", 9, :black, :right))
+
+plot!(Shape([(-2,-6), (-0.55,-6), (-0.55,-7.1), (-2,-7.1)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.95, -6.5, ("local_leaf_cells(mesh.tree)", 10, :black, :left))
+annotate!(-2, -7.5, ("returns cells for which an element needs to be created (i.e. all leaf cells)", 9, :black, :left))
+plot!([-2.22,-2],[-6.5,-6.5],arrow=true,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-2,-9), (1.73,-9), (1.73,-10.1), (-2,-10.1)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.95, -9.5, ("init_elements(leaf_cell_ids, mesh, equations, dg.basis, RealT, uEltype)", 10, :black, :left))
+annotate!(-2, -10.5, ("creates and initializes elements, projects Gauss-Lobatto basis onto each of them", 9, :black, :left))
+plot!([-2.22,-2],[-9.5,-9.5],arrow=true,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-2,-12), (0.4,-12), (0.4,-13.1), (-2,-13.1)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.95, -12.5, ("init_interfaces(leaf_cell_ids, mesh, elements)", 10, :black, :left))
+annotate!(-2, -13.5, ("creates and initializes interfaces between each pair of adjacent elements of the same size", 9, :black, :left))
+plot!([-2.22,-2],[-12.5,-12.5],arrow=true,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-2,-15), (0.5,-15), (0.5,-16.1), (-2,-16.1)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.95, -15.5, ("init_boundaries(leaf_cell_ids, mesh, elements)", 10, :black, :left))
+annotate!(-2, -17, ("creates and initializes boundaries, remembers each boundary element, as well as the coordinates of
+each boundary node", 9, :black, :left))
+plot!([-2.22,-2],[-15.5,-15.5],arrow=true,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-1.6,-18), (1.3,-18), (1.3,-19.1), (-1.6,-19.1)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.55, -18.5, ("init_mortars(leaf_cell_ids, mesh, elements, dg.mortar)", 10, :black, :left))
+annotate!(-1.6, -20, ("creates and initializes mortars (type of interfaces) between each triple of adjacent coarsened
+and corresponding small elements", 9, :black, :left))
+plot!([-2.22,-1.6],[-18.5,-18.5],arrow=true,color=:black,linewidth=2,label="")
+annotate!(-2.15, -19, ("2D and 3D", 8, :black, :left))
+
+plot!(Shape([(-2,-21), (1.5,-21), (1.5,-23.1), (-2,-23.1)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.95, -22, ("create_cache(mesh, equations, dg.volume_integral, dg, uEltype)
+for 2D and 3D create_cache(mesh, equations, dg.mortar, uEltype)", 10, :black, :left))
+annotate!(-2, -23.5, ("add specialized parts of the cache required to compute the volume integral, etc.", 9, :black, :left))
+plot!([-2.22,-2],[-22,-22],arrow=true,color=:black,linewidth=2,label="")
+
+savefig("./SemidiscretizationHyperbolic")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_boundary_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_boundary_figure.jl
new file mode 100644
index 0000000000..14475d2133
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_boundary_figure.jl
@@ -0,0 +1,190 @@
+using Plots
+
+function min(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  min=coordinates[1][i]
+  for j in coordinates
+    if min>j[i]
+      min=j[i]
+    end
+  end
+  return min
+end
+
+function max(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  max=coordinates[1][i]
+  for j in coordinates
+    if max<j[i]
+      max=j[i]
+    end
+  end
+  return max
+end
+
+function nodes_project_four(coordinates::Vector{Tuple{Float64, Float64}})
+  nodes_1d=[-1.0, -0.4472135954999579, 0.4472135954999579, 1.0]
+  nodes_2d=Array{Float64,2}(undef, 16,2)
+  for i in range(1,length(nodes_1d))
+    for j in range(1,length(nodes_1d))
+      nodes_2d[(i-1)*4+j,:]=[nodes_1d[i], nodes_1d[j]]
+    end
+  end
+  k_x = (max(coordinates,1)-min(coordinates,1))/2
+  k_y = (max(coordinates,2)-min(coordinates,2))/2
+  b_x = min(coordinates,1)
+  b_y = min(coordinates,2)
+  for i in range(1,Int(length(nodes_2d)/2))
+    nodes_2d[i,1]=k_x*(nodes_2d[i,1]+1)+b_x
+    nodes_2d[i,2]=k_y*(nodes_2d[i,2]+1)+b_y
+  end
+  return nodes_2d
+end
+
+
+import Base.+
++(a::Tuple,b::Tuple) = a .+ b
+
+# level 1
+plot(Shape([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)]), linecolor="black", fillcolor="white", label="elements",legend_position=:right, linewidth=2, size=(800,600), showaxis=false, grid=false, x_lims=(-3,5.5), ylims=(-3.5,3))
+
+nodes_2d=nodes_project_four([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = "Gauss-Lobatto nodes", markersize = 5)
+
+# mortars
+plot!([-0.25,0.125], [-1,-1], linecolor="orange", fillcolor="white", label="mortars", linewidth=3)
+plot!([-0.25,0.125], [-2.5,-2.5], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([0.125,0.125], [-1,-2.5], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([0.125,0.5], [-1.75,-1.75], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([0.75,0.75], [0,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([2.25,2.25], [0,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([2.25,0.75], [-0.325,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([1.5,1.5], [-0.75,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+
+plot!([2.5,2.5], [-0.75,-2.75], linecolor="green", label="boundaries with nodes", linewidth=15, opacity=0.3)
+plot!([.5,2.5], [-2.75,-2.75], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# level 2
+
+# upper right
+# 1
+plot!(Shape([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([1.25,1.75], [2,2], linecolor="blue", fillcolor="white", label="interfaces", linewidth=3)
+plot!([2.25,2.25], [1.5,1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+plot!([2.75,2.75], [2.5,1.5], linecolor="green", label=false, linewidth=15, opacity=0.3)
+plot!([1.75,2.75], [2.5,2.5], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# 2
+plot!(Shape([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-0.25,0.25], [2,2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([0.75,0.75], [1.5,1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+plot!([.25,1.25], [2.5,2.5], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# 3
+plot!(Shape([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([1.25,1.75], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.25,0.25], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 4
+plot!(Shape([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([2.75,2.75], [1,0], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# upper left
+# 1
+plot!(Shape([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.25,-1.75], [2,2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.75,-0.75], [1,1.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+plot!([-1.25,-0.25], [2.5,2.5], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# 2
+plot!(Shape([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-2.25,-2.25], [1,1.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+plot!([-1.75,-2.75], [2.5,2.5], linecolor="green", label=false, linewidth=15, opacity=0.3)
+plot!([-2.75,-2.75], [1.5,2.5], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# 3
+plot!(Shape([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-2.25,-2.25], [-0.5,0], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+plot!([-2.75,-2.75], [0,1], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# 4
+plot!(Shape([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-0.75,-0.75], [-0.5,0], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# lower left
+# 1
+plot!(Shape([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [-1,-1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.75,-0.75], [-1.5,-2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+
+# 2
+plot!(Shape([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-2.25,-2.25], [-1.5,-2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+plot!([-2.75,-2.75], [-0.5,-1.5], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# 3
+plot!(Shape([(-2,-2),(-2,-1),(-1,-1),(-1,-2)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2.,-2.),(-2.,-1.),(-1.,-1.),(-1.,-2.)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [-2.5,-2.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+plot!([-2.75,-2.75], [-2,-3], linecolor="green", label=false, linewidth=15, opacity=0.3)
+plot!([-2.75,-1.75], [-3,-3], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# 4
+plot!(Shape([(-1,-2),(-1,-1),(-0,-1),(-0,-2)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-1.,-2.),(-1.,-1.),(-0.,-1.),(-0.,-2.)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-0.25,-1.25], [-3,-3], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+savefig("./boundary")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_elements_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_elements_figure.jl
new file mode 100644
index 0000000000..6a98467a4e
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_elements_figure.jl
@@ -0,0 +1,117 @@
+using Plots
+
+function min(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  min=coordinates[1][i]
+  for j in coordinates
+    if min>j[i]
+      min=j[i]
+    end
+  end
+  return min
+end
+
+function max(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  max=coordinates[1][i]
+  for j in coordinates
+    if max<j[i]
+      max=j[i]
+    end
+  end
+  return max
+end
+
+function nodes_project_four(coordinates::Vector{Tuple{Float64, Float64}})
+  nodes_1d=[-1.0, -0.4472135954999579, 0.4472135954999579, 1.0]
+  nodes_2d=Array{Float64,2}(undef, 16,2)
+  for i in range(1,length(nodes_1d))
+    for j in range(1,length(nodes_1d))
+      nodes_2d[(i-1)*4+j,:]=[nodes_1d[i], nodes_1d[j]]
+    end
+  end
+  k_x = (max(coordinates,1)-min(coordinates,1))/2
+  k_y = (max(coordinates,2)-min(coordinates,2))/2
+  b_x = min(coordinates,1)
+  b_y = min(coordinates,2)
+  for i in range(1,Int(length(nodes_2d)/2))
+    nodes_2d[i,1]=k_x*(nodes_2d[i,1]+1)+b_x
+    nodes_2d[i,2]=k_y*(nodes_2d[i,2]+1)+b_y
+  end
+  return nodes_2d
+end
+
+
+import Base.+
++(a::Tuple,b::Tuple) = a .+ b
+
+# level 1
+plot(Shape([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)]), linecolor="black", fillcolor="white", label="elements",legend_position=:right, linewidth=2, size=(800,600), showaxis=false, grid=false, x_lims=(-3,5.5), ylims=(-3.5,3))
+
+nodes_2d=nodes_project_four([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = "Gauss-Lobatto nodes", markersize = 5)
+
+# level 2
+
+# upper right
+plot!(Shape([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+# upper left
+plot!(Shape([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+# lower left
+plot!(Shape([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(-2,-2),(-2,-1),(-1,-1),(-1,-2)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2.,-2.),(-2.,-1.),(-1.,-1.),(-1.,-2.)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(-1,-2),(-1,-1),(-0,-1),(-0,-2)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-1.,-2.),(-1.,-1.),(-0.,-1.),(-0.,-2.)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+savefig("./elements")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_interfaces_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_interfaces_figure.jl
new file mode 100644
index 0000000000..9d5e40722f
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_interfaces_figure.jl
@@ -0,0 +1,157 @@
+using Plots
+
+function min(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  min=coordinates[1][i]
+  for j in coordinates
+    if min>j[i]
+      min=j[i]
+    end
+  end
+  return min
+end
+
+function max(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  max=coordinates[1][i]
+  for j in coordinates
+    if max<j[i]
+      max=j[i]
+    end
+  end
+  return max
+end
+
+function nodes_project_four(coordinates::Vector{Tuple{Float64, Float64}})
+  nodes_1d=[-1.0, -0.4472135954999579, 0.4472135954999579, 1.0]
+  nodes_2d=Array{Float64,2}(undef, 16,2)
+  for i in range(1,length(nodes_1d))
+    for j in range(1,length(nodes_1d))
+      nodes_2d[(i-1)*4+j,:]=[nodes_1d[i], nodes_1d[j]]
+    end
+  end
+  k_x = (max(coordinates,1)-min(coordinates,1))/2
+  k_y = (max(coordinates,2)-min(coordinates,2))/2
+  b_x = min(coordinates,1)
+  b_y = min(coordinates,2)
+  for i in range(1,Int(length(nodes_2d)/2))
+    nodes_2d[i,1]=k_x*(nodes_2d[i,1]+1)+b_x
+    nodes_2d[i,2]=k_y*(nodes_2d[i,2]+1)+b_y
+  end
+  return nodes_2d
+end
+
+
+import Base.+
++(a::Tuple,b::Tuple) = a .+ b
+
+# level 1
+plot(Shape([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)]), linecolor="black", fillcolor="white", label="elements",legend_position=:right, linewidth=2, size=(800,600), showaxis=false, grid=false, x_lims=(-3,5.5), ylims=(-3.5,3))
+
+nodes_2d=nodes_project_four([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = "Gauss-Lobatto nodes", markersize = 5)
+
+# level 2
+
+# upper right
+# 1
+plot!(Shape([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([1.25,1.75], [2,2], linecolor="blue", fillcolor="white", label="interfaces", linewidth=3)
+plot!([2.25,2.25], [1.5,1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 2
+plot!(Shape([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-0.25,0.25], [2,2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([0.75,0.75], [1.5,1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 3
+plot!(Shape([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([1.25,1.75], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.25,0.25], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 4
+plot!(Shape([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+# upper left
+# 1
+plot!(Shape([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.25,-1.75], [2,2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.75,-0.75], [1,1.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 2
+plot!(Shape([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-2.25,-2.25], [1,1.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+
+# 3
+plot!(Shape([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-2.25,-2.25], [-0.5,0], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 4
+plot!(Shape([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-0.75,-0.75], [-0.5,0], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# lower left
+# 1
+plot!(Shape([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [-1,-1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.75,-0.75], [-1.5,-2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+
+# 2
+plot!(Shape([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-2.25,-2.25], [-1.5,-2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 3
+plot!(Shape([(-2,-2),(-2,-1),(-1,-1),(-1,-2)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2.,-2.),(-2.,-1.),(-1.,-1.),(-1.,-2.)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [-2.5,-2.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 4
+plot!(Shape([(-1,-2),(-1,-1),(-0,-1),(-0,-2)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-1.,-2.),(-1.,-1.),(-0.,-1.),(-0.,-2.)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+savefig("./interfaces")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_mortars_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_mortars_figure.jl
new file mode 100644
index 0000000000..f48b410c42
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_mortars_figure.jl
@@ -0,0 +1,166 @@
+using Plots
+
+function min(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  min=coordinates[1][i]
+  for j in coordinates
+    if min>j[i]
+      min=j[i]
+    end
+  end
+  return min
+end
+
+function max(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  max=coordinates[1][i]
+  for j in coordinates
+    if max<j[i]
+      max=j[i]
+    end
+  end
+  return max
+end
+
+function nodes_project_four(coordinates::Vector{Tuple{Float64, Float64}})
+  nodes_1d=[-1.0, -0.4472135954999579, 0.4472135954999579, 1.0]
+  nodes_2d=Array{Float64,2}(undef, 16,2)
+  for i in range(1,length(nodes_1d))
+    for j in range(1,length(nodes_1d))
+      nodes_2d[(i-1)*4+j,:]=[nodes_1d[i], nodes_1d[j]]
+    end
+  end
+  k_x = (max(coordinates,1)-min(coordinates,1))/2
+  k_y = (max(coordinates,2)-min(coordinates,2))/2
+  b_x = min(coordinates,1)
+  b_y = min(coordinates,2)
+  for i in range(1,Int(length(nodes_2d)/2))
+    nodes_2d[i,1]=k_x*(nodes_2d[i,1]+1)+b_x
+    nodes_2d[i,2]=k_y*(nodes_2d[i,2]+1)+b_y
+  end
+  return nodes_2d
+end
+
+
+import Base.+
++(a::Tuple,b::Tuple) = a .+ b
+
+# level 1
+plot(Shape([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)]), linecolor="black", fillcolor="white", label="elements",legend_position=:right, linewidth=2, size=(800,600), showaxis=false, grid=false, x_lims=(-3,5.5), ylims=(-3.5,3))
+
+nodes_2d=nodes_project_four([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = "Gauss-Lobatto nodes", markersize = 5)
+
+plot!([-0.25,0.125], [-1,-1], linecolor="orange", fillcolor="white", label="mortars", linewidth=3)
+plot!([-0.25,0.125], [-2.5,-2.5], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([0.125,0.125], [-1,-2.5], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([0.125,0.5], [-1.75,-1.75], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([0.75,0.75], [0,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([2.25,2.25], [0,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([2.25,0.75], [-0.325,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([1.5,1.5], [-0.75,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+
+# level 2
+
+# upper right
+# 1
+plot!(Shape([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([1.25,1.75], [2,2], linecolor="blue", fillcolor="white", label="interfaces", linewidth=3)
+plot!([2.25,2.25], [1.5,1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 2
+plot!(Shape([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-0.25,0.25], [2,2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([0.75,0.75], [1.5,1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 3
+plot!(Shape([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([1.25,1.75], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.25,0.25], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 4
+plot!(Shape([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+# upper left
+# 1
+plot!(Shape([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.25,-1.75], [2,2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.75,-0.75], [1,1.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 2
+plot!(Shape([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-2.25,-2.25], [1,1.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+
+# 3
+plot!(Shape([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-2.25,-2.25], [-0.5,0], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 4
+plot!(Shape([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-0.75,-0.75], [-0.5,0], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# lower left
+# 1
+plot!(Shape([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [-1,-1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.75,-0.75], [-1.5,-2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+
+# 2
+plot!(Shape([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-2.25,-2.25], [-1.5,-2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 3
+plot!(Shape([(-2,-2),(-2,-1),(-1,-1),(-1,-2)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2.,-2.),(-2.,-1.),(-1.,-1.),(-1.,-2.)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [-2.5,-2.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 4
+plot!(Shape([(-1,-2),(-1,-1),(-0,-1),(-0,-2)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-1.,-2.),(-1.,-1.),(-0.,-1.),(-0.,-2.)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+savefig("./mortars")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_nodes_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_nodes_figure.jl
new file mode 100644
index 0000000000..b983563f47
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_nodes_figure.jl
@@ -0,0 +1,6 @@
+using Plots
+
+plot([-1,1],[0,0], linecolor="black", label="reference interval", size=(600, 200), ylim=(-0.5,0.5), yaxis=false, grid=false)
+scatter!([-1.0, -0.4472135954999579, 0.4472135954999579, 1.0],[0,0,0,0], color="red", label="Gauss-Lobatto nodes")
+
+savefig("./nodes")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_treemesh_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_treemesh_figure.jl
new file mode 100644
index 0000000000..8d2a8173a1
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_treemesh_figure.jl
@@ -0,0 +1,26 @@
+using Plots
+# level 0
+plot(Shape([(2,2),(2,-2),(-2,-2),(-2,2)]), linecolor="black", fillcolor="white", label="0th level, root", linewidth=2, size=(600,600))
+# level 1
+plot!(Shape([(1.97,1.97),(1.97,0),(0,0),(0,1.97)]), linecolor="red", fillcolor="white", label="1st level", linewidth=2)
+plot!(Shape([(-1.97,1.97),(0,1.97),(0,0),(-1.97,0)]), linecolor="red", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(-1.97,0),(0,0),(0,-1.97),(-1.97,-1.97)]), linecolor="red", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(1.97,0),(1.97,-1.97),(0,-1.97),(0,0)]), linecolor="red", fillcolor="white", label=false, linewidth=2)
+# level 2
+# upper right
+plot!(Shape([(1.94,1.94),(1.94,1),(1,1),(1,1.94)]), linecolor="blue", fillcolor="white", label="2nd level", linewidth=2)
+plot!(Shape([(0.03,1.94),(1,1.94),(1,1),(0.03,1)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(0.03,0.03),(0.03,1),(1,1),(1,0.03)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(1.94,0.03),(1.94,1),(1,1),(1,0.03)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+# upper left
+plot!(Shape([(-0.03,1.94),(-0.03,1),(-1,1),(-1,1.94)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(-1.94,1.94),(-1,1.94),(-1,1),(-1.94,1)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(-1.94,0.03),(-1.94,1),(-1,1),(-1,0.03)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(-0.03,0.03),(-1,0.03),(-1,1),(-0.03,1)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+# lower left
+plot!(Shape([(-0.03,-0.03),(-0.03,-1),(-1,-1),(-1,-0.03)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(-1.94,-1),(-1.94,-0.03),(-1,-0.03),(-1,-1)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(-1.94,-1.94),(-1.94,-1),(-1,-1),(-1,-1.94)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(-1,-1.94),(-1,-1),(-0.03,-1),(-0.03,-1.94)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+
+savefig("./treemesh")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/rhs_structure_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/rhs_structure_figure.jl
new file mode 100644
index 0000000000..19083dad68
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/rhs_structure_figure.jl
@@ -0,0 +1,43 @@
+using Plots
+plot(Shape([(-0.4,0.5), (0.4,0.5), (0.4,-0.5), (-0.4,-0.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2, size=(800,600), showaxis=false, grid=false, xlim=(-2,2), ylim=(-17.5,3))
+annotate!(0, 0, ("solve(...) call", 12, :black, :center))
+plot!([0,0],[-0.5,-2.5],arrow=true,color=:black,linewidth=2,label="")
+plot!(Shape([(-1.2,1), (1.2,1), (1.2,-1), (-1.2,-1)]), linecolor=:green, fillcolor=:transparent, label=false,linewidth=1, linestyle=:dash)
+annotate!(0.8, 0, ("Trixi.jl setup", 12, :green, :center))
+
+plot!(Shape([(-1.25,-2.5), (1.25,-2.5), (1.25,-3.5), (-1.25,-3.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0, -3, ("DiffEqBase.__solve(prob, alg, args...; kwargs...)", 12, :black, :center))
+plot!([0,0],[-3.5,-5.5],arrow=true,color=:black,linewidth=2,label="")
+plot!(Shape([(-1.9,-2), (1.9,-2), (1.9,-10), (-1.9,-10)]), linecolor=:red, fillcolor=:transparent, label=false,linewidth=1, linestyle=:dash)
+annotate!(1.4, -6, ("OrdinaryDiffEq.jl", 12, :red, :center))
+
+plot!(Shape([(-0.85,-5.5), (0.85,-5.5), (0.85,-6.5), (-0.85,-6.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0, -6, ("DiffEqBase.solve!(integrator)", 12, :black, :center))
+plot!([0,0],[-6.5,-8.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,-0.5],[-6.5,-8.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,0.5],[-6.5,-8.5],arrow=true,color=:black,linewidth=2,label="")
+annotate!(0.8, -7.5, ("Specialized for time \nintegration method", 9, :black, :center))
+
+plot!(Shape([(-1.2,-8.5), (1.2,-8.5), (1.2,-9.5), (-1.2,-9.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0, -9, ("perform_step!(integrator, integrator.cache)", 12, :black, :center))
+plot!([0,0],[-9.5,-11.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,-0.5],[-9.5,-11.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,0.5],[-9.5,-11.5],arrow=true,color=:black,linewidth=2,label="")
+annotate!(0.8, -10.5, ("Specialized for \nsemidiscretization", 9, :black, :center))
+
+plot!(Shape([(-0.95,-11.5), (0.95,-11.5), (0.95,-12.5), (-0.95,-12.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0, -12, ("Trixi.rhs!(du_ode, u_ode, semi, t)", 12, :black, :center))
+plot!([0,0],[-12.5,-14.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,-0.5],[-12.5,-14.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,0.5],[-12.5,-14.5],arrow=true,color=:black,linewidth=2,label="")
+annotate!(0.8, -13.5, ("Specialized for \nmesh and dg", 9, :black, :center))
+plot!(Shape([(-1.7,-11), (1.7,-11), (1.7,-17), (-1.7,-17)]), linecolor=:blue, fillcolor=:transparent, label=false,linewidth=1, linestyle=:dash)
+annotate!(1.5, -14, ("Trixi.jl", 12, :blue, :center))
+
+plot!(Shape([(-1.4,-14.5), (1.4,-14.5), (1.4,-16.5), (-1.4,-16.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0, -15.5, ("Trixi.rhs!(du, u, t, mesh, equations, initial_condition, \nboundary_conditions, source_terms, dg, cache)", 12, :black, :center))
+
+
+
+plot!([-2,2], [2,2], linecolor="white", fillcolor="white", label=false,linewidth=2)
+savefig("./rhs!")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/semidiscretize_structure_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/semidiscretize_structure_figure.jl
new file mode 100644
index 0000000000..52b0d573d7
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/semidiscretize_structure_figure.jl
@@ -0,0 +1,51 @@
+using Plots
+plot(Shape([(-1.2,4), (1.2,4), (1.2,3), (-1.2,3)]), linecolor="black", fillcolor="white", label=false,linewidth=2, size=(800,600), showaxis=false, grid=false, xlim=(-2.6,2.6), ylim=(-21.5,4.5))
+annotate!(0, 3.5, ("semidiscretize(semi, tspan; reset_threads)", 10, :black, :center))
+annotate!(0, 2.5, ("creates and returns an ODEProblem object, initialized using rhs!, u0_ode, tspan and semi", 9, :black, :center))
+plot!([0,0],[2,0.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,-0.2],[2,0.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,0.2],[2,0.5],arrow=true,color=:black,linewidth=2,label="")
+annotate!(0.3,1.25, ("specialized for semi type", 9, :black, :left))
+
+plot!(Shape([(-1.6,0.5), (1.6,0.5), (1.6,-0.5), (-1.6,-0.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0, 0, ("compute_coefficients(t, semi::SemidiscretizationHyperbolic)", 10, :black, :center))
+plot!([0,0],[-0.5,-2],arrow=true,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-1.3,-2), (1.3,-2), (1.3,-3), (-1.3,-3)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0, -2.5, ("compute_coefficients(initial_conditions, t, semi)", 10, :black, :center))
+plot!([0,0.2],[-3,-4.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,-0.87],[-3,-9.5],arrow=true,color=:black,linewidth=2,label="")
+
+plot!(Shape([(0.1,-4.5), (2.5,-4.5), (2.5,-6.5), (0.1,-6.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0.15, -5.5, ("allocate_coefficients(mesh, equations, solver,
+                                 cache)", 10, :black, :left))
+annotate!(0.1, -8, ("initializes u_ode as a 1D zero-vector with a length
+depending on the number of variables, nodes,
+dimensions of mesh and elements, this vector is used
+by OrdinaryDiffEq.jl to keep a solution", 9, :black, :left))
+
+plot!(Shape([(-2.5,-9.5), (0.05,-9.5), (0.05,-11.5), (-2.5,-11.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-2.45, -10.5, ("compute_coefficients!(u_ode, initial_conditions,
+                                    t, semi)", 10, :black, :left))
+plot!([-2.4,-2.4],[-11.5,-18.5],arrow=false,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-2.2,-13), (-0.8,-13), (-0.8,-14), (-2.2,-14)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.5, -13.5, ("wrap_array(u_ode, semi)", 10, :black, :center))
+annotate!(-2.2, -15.2, ("to simplify processing, it wraps the 1D array u_ode created in the allocate_coefficients function in
+a multidimensional one, where dimensions are responsible for variables, nodes on each axis and
+elements", 9, :black, :left))
+plot!([-2.4,-2.2],[-13.5,-13.5],arrow=true,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-1.1,-17.5), (1.7,-17.5), (1.7,-19.5), (-1.1,-19.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.05, -18.5, ("compute_coefficients!(u, initial_conditions, t, mesh,
+                                    equations, solver, cache)", 10, :black, :left))
+annotate!(-1.1, -20.5, ("applies an initial conditions to each node of each element for each variable,
+saves in the wrapped u_ode", 9, :black, :left))
+plot!([-2.4,-1.1],[-18.5,-18.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([-2.4,-1.1],[-18.5,-17.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([-2.4,-1.1],[-18.5,-19.5],arrow=true,color=:black,linewidth=2,label="")
+annotate!(-2.4, -20.2, ("specialized for solver
+type and dimensionality
+of mesh", 9, :black, :left))
+
+savefig("./semidiscretize")
\ No newline at end of file
diff --git a/docs/literate/src/files/index.jl b/docs/literate/src/files/index.jl
index 1fc025d84d..6d5800c3b6 100644
--- a/docs/literate/src/files/index.jl
+++ b/docs/literate/src/files/index.jl
@@ -21,20 +21,29 @@
 # fundamental structure of a Trixi.jl setup, the visualization of results, and the development
 # process for Trixi.jl.
 
-# ### [2 Introduction to DG methods](@ref scalar_linear_advection_1d)
+# ### [2 Behind the scenes of a simulation setup](@ref behind_the_scenes_simulation_setup)
+#-
+# This tutorial will guide you through a simple Trixi.jl setup ("elixir"), giving an overview of
+# what happens in the background during the initialization of a simulation. While the setup
+# described herein does not cover all details, it involves relatively stable parts of Trixi.jl that
+# are unlikely to undergo significant changes in the near future. The goal is to clarify some of
+# the more fundamental, *technical* concepts that are applicable to a variety of
+# (also more complex) configurations.s
+
+# ### [3 Introduction to DG methods](@ref scalar_linear_advection_1d)
 #-
 # This tutorial gives an introduction to discontinuous Galerkin (DG) methods with the example of the
 # scalar linear advection equation in 1D. Starting with some theoretical explanations, we first implement
 # a raw version of a discontinuous Galerkin spectral element method (DGSEM). Then, we will show how
 # to use features of Trixi.jl to achieve the same result.
 
-# ### [3 DGSEM with flux differencing](@ref DGSEM_FluxDiff)
+# ### [4 DGSEM with flux differencing](@ref DGSEM_FluxDiff)
 #-
 # To improve stability often the flux differencing formulation of the DGSEM (split form) is used.
 # We want to present the idea and formulation on a basic 1D level. Then, we show how this formulation
 # can be implemented in Trixi.jl and analyse entropy conservation for two different flux combinations.
 
-# ### [4 Shock capturing with flux differencing and stage limiter](@ref shock_capturing)
+# ### [5 Shock capturing with flux differencing and stage limiter](@ref shock_capturing)
 #-
 # Using the flux differencing formulation, a simple procedure to capture shocks is a hybrid blending
 # of a high-order DG method and a low-order subcell finite volume (FV) method. We present the idea on a
@@ -42,20 +51,20 @@
 # explained and added to an exemplary simulation of the Sedov blast wave with the 2D compressible Euler
 # equations.
 
-# ### [5 Non-periodic boundary conditions](@ref non_periodic_boundaries)
+# ### [6 Non-periodic boundary conditions](@ref non_periodic_boundaries)
 #-
 # Thus far, all examples used periodic boundaries. In Trixi.jl, you can also set up a simulation with
 # non-periodic boundaries. This tutorial presents the implementation of the classical Dirichlet
 # boundary condition with a following example. Then, other non-periodic boundaries are mentioned.
 
-# ### [6 DG schemes via `DGMulti` solver](@ref DGMulti_1)
+# ### [7 DG schemes via `DGMulti` solver](@ref DGMulti_1)
 #-
 # This tutorial is about the more general DG solver [`DGMulti`](@ref), introduced [here](@ref DGMulti).
 # We are showing some examples for this solver, for instance with discretization nodes by Gauss or
 # triangular elements. Moreover, we present a simple way to include pre-defined triangulate meshes for
 # non-Cartesian domains using the package [StartUpDG.jl](https://github.com/jlchan/StartUpDG.jl).
 
-# ### [7 Other SBP schemes (FD, CGSEM) via `DGMulti` solver](@ref DGMulti_2)
+# ### [8 Other SBP schemes (FD, CGSEM) via `DGMulti` solver](@ref DGMulti_2)
 #-
 # Supplementary to the previous tutorial about DG schemes via the `DGMulti` solver we now present
 # the possibility for `DGMulti` to use other SBP schemes via the package
@@ -63,7 +72,7 @@
 # For instance, we show how to set up a finite differences (FD) scheme and a continuous Galerkin
 # (CGSEM) method.
 
-# ### [8 Upwind FD SBP schemes](@ref upwind_fdsbp)
+# ### [9 Upwind FD SBP schemes](@ref upwind_fdsbp)
 #-
 # General SBP schemes can not only be used via the [`DGMulti`](@ref) solver but
 # also with a general `DG` solver. In particular, upwind finite difference SBP
@@ -71,42 +80,42 @@
 # schemes in the `DGMulti` framework, the interface is based on the package
 # [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl).
 
-# ### [9 Adding a new scalar conservation law](@ref adding_new_scalar_equations)
+# ### [10 Adding a new scalar conservation law](@ref adding_new_scalar_equations)
 #-
 # This tutorial explains how to add a new physics model using the example of the cubic conservation
 # law. First, we define the equation using a `struct` `CubicEquation` and the physical flux. Then,
 # the corresponding standard setup in Trixi.jl (`mesh`, `solver`, `semi` and `ode`) is implemented
 # and the ODE problem is solved by OrdinaryDiffEq's `solve` method.
 
-# ### [10 Adding a non-conservative equation](@ref adding_nonconservative_equation)
+# ### [11 Adding a non-conservative equation](@ref adding_nonconservative_equation)
 #-
 # In this part, another physics model is implemented, the nonconservative linear advection equation.
 # We run two different simulations with different levels of refinement and compare the resulting errors.
 
-# ### [11 Parabolic terms](@ref parabolic_terms)
+# ### [12 Parabolic terms](@ref parabolic_terms)
 #-
 # This tutorial describes how parabolic terms are implemented in Trixi.jl, e.g.,
 # to solve the advection-diffusion equation.
 
-# ### [12 Adding new parabolic terms](@ref adding_new_parabolic_terms)
+# ### [13 Adding new parabolic terms](@ref adding_new_parabolic_terms)
 #-
 # This tutorial describes how new parabolic terms can be implemented using Trixi.jl.
 
-# ### [13 Adaptive mesh refinement](@ref adaptive_mesh_refinement)
+# ### [14 Adaptive mesh refinement](@ref adaptive_mesh_refinement)
 #-
 # Adaptive mesh refinement (AMR) helps to increase the accuracy in sensitive or turbolent regions while
 # not wasting resources for less interesting parts of the domain. This leads to much more efficient
 # simulations. This tutorial presents the implementation strategy of AMR in Trixi.jl, including the use of
 # different indicators and controllers.
 
-# ### [14 Structured mesh with curvilinear mapping](@ref structured_mesh_mapping)
+# ### [15 Structured mesh with curvilinear mapping](@ref structured_mesh_mapping)
 #-
 # In this tutorial, the use of Trixi.jl's structured curved mesh type [`StructuredMesh`](@ref) is explained.
 # We present the two basic option to initialize such a mesh. First, the curved domain boundaries
 # of a circular cylinder are set by explicit boundary functions. Then, a fully curved mesh is
 # defined by passing the transformation mapping.
 
-# ### [15 Unstructured meshes with HOHQMesh.jl](@ref hohqmesh_tutorial)
+# ### [16 Unstructured meshes with HOHQMesh.jl](@ref hohqmesh_tutorial)
 #-
 # The purpose of this tutorial is to demonstrate how to use the [`UnstructuredMesh2D`](@ref)
 # functionality of Trixi.jl. This begins by running and visualizing an available unstructured
@@ -115,26 +124,26 @@
 # software in the Trixi.jl ecosystem, and then run a simulation using Trixi.jl on said mesh.
 # In the end, the tutorial briefly explains how to simulate an example using AMR via `P4estMesh`.
 
-# ### [16 P4est mesh from gmsh](@ref p4est_from_gmsh)
+# ### [17 P4est mesh from gmsh](@ref p4est_from_gmsh)
 #-
 # This tutorial describes how to obtain a [`P4estMesh`](@ref) from an existing mesh generated
 # by [`gmsh`](https://gmsh.info/) or any other meshing software that can export to the Abaqus
 # input `.inp` format. The tutorial demonstrates how edges/faces can be associated with boundary conditions based on the physical nodesets.
 
-# ### [17 Explicit time stepping](@ref time_stepping)
+# ### [18 Explicit time stepping](@ref time_stepping)
 #-
 # This tutorial is about time integration using [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl).
 # It explains how to use their algorithms and presents two types of time step choices - with error-based
 # and CFL-based adaptive step size control.
 
-# ### [18 Differentiable programming](@ref differentiable_programming)
+# ### [19 Differentiable programming](@ref differentiable_programming)
 #-
 # This part deals with some basic differentiable programming topics. For example, a Jacobian, its
 # eigenvalues and a curve of total energy (through the simulation) are calculated and plotted for
 # a few semidiscretizations. Moreover, we calculate an example for propagating errors with Measurement.jl
 # at the end.
 
-# ### [19 Custom semidiscretization](@ref custom_semidiscretization)
+# ### [20 Custom semidiscretization](@ref custom_semidiscretization)
 #-
 # This tutorial describes the [semidiscretiations](@ref overview-semidiscretizations) of Trixi.jl
 # and explains how to extend them for custom tasks.
diff --git a/docs/make.jl b/docs/make.jl
index 946b803b71..8427c4049b 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -55,6 +55,7 @@ files = [
         "Create first setup" => ("first_steps", "create_first_setup.jl"),
         "Changing Trixi.jl itself" => ("first_steps", "changing_trixi.jl"),
     ],
+    "Behind the scenes of a simulation setup" => "behind_the_scenes_simulation_setup.jl",
     # Topic: DG semidiscretizations
     "Introduction to DG methods" => "scalar_linear_advection_1d.jl",
     "DGSEM with flux differencing" => "DGSEM_FluxDiff.jl",
@@ -76,7 +77,7 @@ files = [
     # Topic: other stuff
     "Explicit time stepping" => "time_stepping.jl",
     "Differentiable programming" => "differentiable_programming.jl",
-    "Custom semidiscretizations" => "custom_semidiscretization.jl"
+    "Custom semidiscretizations" => "custom_semidiscretization.jl",
     ]
 tutorials = create_tutorials(files)