diff --git a/docs/src/preprocessing.md b/docs/src/preprocessing.md
index f50aaeaa6..990ddecbc 100644
--- a/docs/src/preprocessing.md
+++ b/docs/src/preprocessing.md
@@ -10,7 +10,7 @@ Pages = map(file -> joinpath("preprocessing", file), readdir(joinpath("..", "src
 ### [References](@id references_complex_shape)
 - Alec Jacobson, Ladislav Kavan, and Olga Sorkine-Hornung "Robust inside-outside segmentation using generalized winding numbers".
   In: ACM Transactions on Graphics, 32.4 (2013), pages 1--12.
-  [doi: 10.1145/2461912.2461916](https://doi.org/10.1145/2461912.2461916)
+  [doi: 10.1145/2461912.2461916](https://igl.ethz.ch/projects/winding-number/robust-inside-outside-segmentation-using-generalized-winding-numbers-siggraph-2013-jacobson-et-al.pdf)
 - Kai Horman, Alexander Agathos "The point in polygon problem for arbitrary polygons".
   In: Computational Geometry, 20.3 (2001), pages 131--144.
   [doi: 10.1016/s0925-7721(01)00012-8](https://doi.org/10.1016/S0925-7721(01)00012-8)
diff --git a/examples/preprocessing/complex_shape_2d.jl b/examples/preprocessing/complex_shape_2d.jl
index 7e271f999..e673285c7 100644
--- a/examples/preprocessing/complex_shape_2d.jl
+++ b/examples/preprocessing/complex_shape_2d.jl
@@ -2,10 +2,10 @@ using TrixiParticles
 
 particle_spacing = 0.05
 
-file = "hexagon"
-filename = joinpath("examples", "preprocessing", file * ".asc")
+filename = "hexagon"
+file = joinpath("examples", "preprocessing", "data", filename * ".asc")
 
-geometry = load_geometry(filename)
+geometry = load_geometry(file)
 
 trixi2vtk(geometry)
 
@@ -13,7 +13,7 @@ point_in_geometry_algorithm = WindingNumberJacobson(; geometry,
                                                     #winding_number_factor=0.4,
                                                     hierarchical_winding=true)
 
-# Returns `InitialCondition`.
+# Returns `InitialCondition`
 shape_sampled = ComplexShape(geometry; particle_spacing, density=1.0,
                              store_winding_number=true,
                              point_in_geometry_algorithm)
diff --git a/examples/preprocessing/complex_shape_3d.jl b/examples/preprocessing/complex_shape_3d.jl
index da0b3467e..a0fa1e1ab 100644
--- a/examples/preprocessing/complex_shape_3d.jl
+++ b/examples/preprocessing/complex_shape_3d.jl
@@ -2,20 +2,16 @@ using TrixiParticles
 
 particle_spacing = 0.05
 
-file = "sphere"
-filename = joinpath("examples", "preprocessing", file * ".stl")
+filename = "sphere"
+file = joinpath("examples", "preprocessing", "data", filename * ".stl")
 
-# The following triangle mesh is corrupt.
-# For more robustness, use `winding_number_factor=0.4`.
-# filename = joinpath("examples", "preprocessing", "drive_gear.stl")
-
-geometry = load_geometry(filename)
+geometry = load_geometry(file)
 
 point_in_geometry_algorithm = WindingNumberJacobson(; geometry,
                                                     # winding_number_factor=0.4,
                                                     hierarchical_winding=true)
 
-# Returns `InitialCondition`.
+# Returns `InitialCondition`
 shape_sampled = ComplexShape(geometry; particle_spacing, density=1.0,
                              point_in_geometry_algorithm)
 
diff --git a/examples/preprocessing/bar.stl b/examples/preprocessing/data/bar.stl
similarity index 100%
rename from examples/preprocessing/bar.stl
rename to examples/preprocessing/data/bar.stl
diff --git a/examples/preprocessing/circle.asc b/examples/preprocessing/data/circle.asc
similarity index 100%
rename from examples/preprocessing/circle.asc
rename to examples/preprocessing/data/circle.asc
diff --git a/examples/preprocessing/data/gear.stl b/examples/preprocessing/data/gear.stl
new file mode 100644
index 000000000..ca5e1b1a1
Binary files /dev/null and b/examples/preprocessing/data/gear.stl differ
diff --git a/examples/preprocessing/hexagon.asc b/examples/preprocessing/data/hexagon.asc
similarity index 100%
rename from examples/preprocessing/hexagon.asc
rename to examples/preprocessing/data/hexagon.asc
diff --git a/examples/preprocessing/inverted_open_curve.asc b/examples/preprocessing/data/inverted_open_curve.asc
similarity index 100%
rename from examples/preprocessing/inverted_open_curve.asc
rename to examples/preprocessing/data/inverted_open_curve.asc
diff --git a/examples/preprocessing/pipe.stl b/examples/preprocessing/data/pipe.stl
similarity index 100%
rename from examples/preprocessing/pipe.stl
rename to examples/preprocessing/data/pipe.stl
diff --git a/examples/preprocessing/sphere.stl b/examples/preprocessing/data/sphere.stl
similarity index 100%
rename from examples/preprocessing/sphere.stl
rename to examples/preprocessing/data/sphere.stl
diff --git a/examples/preprocessing/drive_gear.stl b/examples/preprocessing/drive_gear.stl
deleted file mode 100644
index d50fd98cd..000000000
Binary files a/examples/preprocessing/drive_gear.stl and /dev/null differ
diff --git a/examples/preprocessing/gear.stl b/examples/preprocessing/gear.stl
deleted file mode 100644
index 4f9da8936..000000000
Binary files a/examples/preprocessing/gear.stl and /dev/null differ
diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
index 1fb998751..637a8bc08 100644
--- a/src/TrixiParticles.jl
+++ b/src/TrixiParticles.jl
@@ -6,10 +6,10 @@ using Adapt: Adapt
 using CSV: CSV
 using Dates
 using DataFrames: DataFrame
-using DelimitedFiles
+using DelimitedFiles: DelimitedFiles
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect, PresetTimeCallback
 using FastPow: @fastpow
-using FileIO: FileIO, File, query, @format_str, Stream, stream
+using FileIO: FileIO
 using ForwardDiff: ForwardDiff
 using GPUArraysCore: AbstractGPUArray
 using JSON: JSON
diff --git a/src/preprocessing/geometries/io.jl b/src/preprocessing/geometries/io.jl
index 4f680f96b..6f30e584b 100644
--- a/src/preprocessing/geometries/io.jl
+++ b/src/preprocessing/geometries/io.jl
@@ -18,7 +18,7 @@ function load_geometry(filename; element_type=Float64)
     if file_extension == ".asc"
         geometry = load_ascii(filename; ELTYPE, skipstart=1)
     elseif file_extension == ".stl"
-        geometry = load(query(filename); ELTYPE)
+        geometry = load(FileIO.query(filename); ELTYPE)
     else
         throw(ArgumentError("Only `.stl` and `.asc` files are supported (yet)."))
     end
@@ -37,17 +37,17 @@ end
 
 # FileIO.jl docs:
 # https://juliaio.github.io/FileIO.jl/stable/implementing/#All-at-once-I/O:-implementing-load-and-save
-function load(fn::File{format"STL_BINARY"}; element_types...)
+function load(fn::FileIO.File{FileIO.format"STL_BINARY"}; element_types...)
     open(fn) do s
         FileIO.skipmagic(s) # skip over the magic bytes
         load(s; element_types...)
     end
 end
 
-function load(fs::Stream{format"STL_BINARY"}; ELTYPE=Float64)
+function load(fs::FileIO.Stream{FileIO.format"STL_BINARY"}; ELTYPE=Float64)
     # Binary STL
     # https://en.wikipedia.org/wiki/STL_%28file_format%29#Binary_STL
-    io = stream(fs)
+    io = FileIO.stream(fs)
     read(io, 80) # Throw out 80 bytes header
     n_faces = read(io, UInt32)
 
diff --git a/src/preprocessing/point_in_poly/hierarchical_winding.jl b/src/preprocessing/point_in_poly/hierarchical_winding.jl
index 9f300691f..50eaa76da 100644
--- a/src/preprocessing/point_in_poly/hierarchical_winding.jl
+++ b/src/preprocessing/point_in_poly/hierarchical_winding.jl
@@ -55,6 +55,7 @@ struct HierarchicalWinding{BB}
     bounding_box::BB
 
     function HierarchicalWinding(geometry)
+        # Note that overlapping bounding boxes are perfectly fine
         min_corner = geometry.min_corner .- sqrt(eps())
         max_corner = geometry.max_corner .+ sqrt(eps())
 
@@ -134,6 +135,11 @@ function determine_closure!(closing_faces, min_corner, max_corner, mesh::Triangl
             end
 
         else # Face is intersecting the boundaries
+
+            # Remove intersecting faces in the calculation of the exterior edges.
+            # This way, all exterior edges are inside the bounding box, and so will be the closing surface.
+            # The intersecting faces are later added with opposite orientation,
+            # so that the closing is actually a closing of the exterior plus intersecting faces.
             push!(intersecting_faces, face)
         end
     end
@@ -164,6 +170,7 @@ function determine_closure!(closing_faces, min_corner, max_corner, mesh::Triangl
         end
     end
 
+    # See comment above as to why intersecting faces are treated seperately
     for face in intersecting_faces
         v1 = face_vertices_ids[face][1]
         v2 = face_vertices_ids[face][2]
@@ -189,28 +196,27 @@ function determine_closure!(closing_edges, min_corner, max_corner, polygon::Poly
     for edge in edges
         if is_in_box(edge_vertices[edge][1], min_corner, max_corner) &&
            is_in_box(edge_vertices[edge][2], min_corner, max_corner)
+            # Edge is completely inside the bounding box
+
             v1 = edge_vertices_ids[edge][1]
             v2 = edge_vertices_ids[edge][2]
 
             vertex_count[v1] += 1
             vertex_count[v2] -= 1
-        else
+        else # Edge is intersecting the boundaries
+
+            # Remove intersecting edges in the calculation of the exterior vertices.
+            # This way, all exterior vertices are inside the bounding box,
+            # and so will be the closing surface.
+            # The intersecting edges are later added with opposite orientation,
+            # so that the closing is actually a closing of the exterior plus intersecting edges.
             push!(intersecting_edges, edge)
         end
     end
 
     exterior_vertices = sort(findall(!iszero, vertex_count))
-
     resize!(closing_edges, 0)
 
-    for edge in intersecting_edges
-        v1 = edge_vertices_ids[edge][1]
-        v2 = edge_vertices_ids[edge][2]
-
-        # Flip order of vertices
-        push!(closing_edges, (v2, v1))
-    end
-
     if !isempty(exterior_vertices)
         closing_vertex = first(exterior_vertices)
     end
@@ -226,6 +232,15 @@ function determine_closure!(closing_edges, min_corner, max_corner, polygon::Poly
         push!(closing_edges, edge)
     end
 
+    # See comment above as to why intersecting edges are treated seperately
+    for edge in intersecting_edges
+        v1 = edge_vertices_ids[edge][1]
+        v2 = edge_vertices_ids[edge][2]
+
+        # Flip order of vertices
+        push!(closing_edges, (v2, v1))
+    end
+
     return closing_edges
 end
 
diff --git a/src/preprocessing/point_in_poly/winding_number_jacobson.jl b/src/preprocessing/point_in_poly/winding_number_jacobson.jl
index f11743543..793fe3467 100644
--- a/src/preprocessing/point_in_poly/winding_number_jacobson.jl
+++ b/src/preprocessing/point_in_poly/winding_number_jacobson.jl
@@ -51,9 +51,8 @@ Algorithm for inside-outside segmentation of a complex geometry proposed by Jaco
 # Keywords
 - `geometry`: Complex geometry returned by [`load_geometry`](@ref) and is only required when using
               `hierarchical_winding=true`.
-- `hierarchical_winding`: If set to `true`, an optimized hierarchical approach will be used
+- `hierarchical_winding`: If set to `true`, an optimized hierarchical approach will be used,
                           which gives a significant speedup.
-                          It is only supported for 3D geometries yet.
 - `winding_number_factor`: For leaky geometries, a factor of `0.4` will give a better inside-outside segmentation.
 
 !!! warning "Experimental Implementation"
@@ -104,14 +103,20 @@ end
 
 # The following functions distinguish between actual triangles (edges)
 # and reconstructed triangles (edges) in the hierarchical winding approach.
+# When we reconstruct the triangles then we pass the index of the vertices. Otherwise,
+# we pass the coordinates of the vertices.
 
-# `face` holds the coordinates of each vertex
+# This method is used, when `naive_winding` is called with `(winding::NaiveWinding)`.
+# `face` holds the coordinates of each vertex.
 @inline face_vertex(mesh, face, index) = face[index]
 
 # `edge` holds the coordinates of each vertex
 @inline edge_vertex(mesh, edge, index) = edge[index]
 
-# `face` holds the index of each vertex
+
+# This method is used, when `naive_winding` is called with `(winding::HierarchicalWinding)`
+# and the query point is outside the bounding box. That is, we use the closure of the box.
+# `face` holds the index of each vertex.
 @inline function face_vertex(mesh, face::NTuple{3, Int}, index)
     v_id = face[index]
 
@@ -126,7 +131,9 @@ end
 end
 
 
-# `face` is the index of the face
+# This method is used, when `naive_winding` is called with `(winding::HierarchicalWinding)`
+# and the bounding box is a leaf.
+# `face` is the index of the face.
 @inline function face_vertex(mesh, face::Int, index)
     (; face_vertices) = mesh
 
diff --git a/src/setups/complex_shape.jl b/src/setups/complex_shape.jl
index 300c089ec..23fe3cdc0 100644
--- a/src/setups/complex_shape.jl
+++ b/src/setups/complex_shape.jl
@@ -7,13 +7,13 @@
                  grid_offset::Real=0.0, max_nparticles=10^7,
                  pad_initial_particle_grid=2particle_spacing)
 
-Complex geometry filled with particles. Returns an [`InitialCondition`](@ref).
+Sample a complex geometry with particles. Returns an [`InitialCondition`](@ref).
 Note that an initial particle grid is generated inside the bounding box of the geometry.
-An `point_in_geometry_algorithm` checks if particles are inside the geometry or not.
+A `point_in_geometry_algorithm` checks if particles are inside the geometry or not.
 For more information about the method see [`WindingNumberJacobson`](@ref) or [`WindingNumberHorman`](@ref).
 
 # Arguments
-- `geometry`: Complex geometry returned by [`load_geometry`](@ref).
+- `geometry`: Geometry returned by [`load_geometry`](@ref).
 
 # Keywords
 - `particle_spacing`:   Spacing between the particles.
@@ -28,7 +28,6 @@ For more information about the method see [`WindingNumberJacobson`](@ref) or [`W
 - `pressure`:           Scalar to set the pressure of all particles to this value.
                         This is only used by the [`EntropicallyDampedSPHSystem`](@ref) and
                         will be overwritten when using an initial pressure function in the system.
-                        Cannot be used together with hydrostatic pressure gradient.
 - `point_in_geometry_algorithm`: Algorithm for sampling the complex geometry with particles.
                               It basically checks whether a particle is inside an object or not.
                               For more information see [`WindingNumberJacobson`](@ref) or [`WindingNumberHorman`](@ref)