Merge branch 'main' into test_nightly

.github/workflows/SpellCheck.yml

@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/typos@v1.23.6
+        uses: crate-ci/typos@v1.24.3

README.md

@@ -12,7 +12,8 @@
 
 **TrixiParticles.jl** is a numerical simulation framework designed for particle-based numerical methods, with an emphasis on multiphysics applications, written in [Julia](https://julialang.org).
 A primary goal of the framework is to be user-friendly for engineering, science, and educational purposes. In addition to its extensible design and optimized implementation, we prioritize the user experience, including installation, pre- and postprocessing.
-Its features include:
+
+[![YouTube](https://github.com/user-attachments/assets/dc2be627-a799-4bfd-9226-2077f737c4b0)](https://www.youtube.com/watch?v=V7FWl4YumcA&t=4667s)
 
 ## Features
 - Incompressible Navier-Stokes

docs/make.jl

@@ -104,8 +104,8 @@ makedocs(sitename="TrixiParticles.jl",
              "Tutorial" => "tutorial.md",
              "Examples" => "examples.md",
              "Visualization" => "visualization.md",
-             "preprocessing" => [
-                 "Sample Geometries" => joinpath("preprocessing", "preprocessing.md"),
+             "Preprocessing" => [
+                 "Sampling of Geometries" => joinpath("preprocessing", "preprocessing.md"),
              ],
              "Components" => [
                  "Overview" => "overview.md",

docs/src/install.md

@@ -1,8 +1,8 @@
 # [Installation](@id installation)
 
 ## Setting up Julia
-If you have not yet installed Julia, please [follow the instructions for your
-operating system](https://julialang.org/downloads/platform/). TrixiParticles.jl works
+If you have not yet installed Julia, please [follow the instructions on the
+official website](https://julialang.org/downloads/). TrixiParticles.jl works
 with Julia v1.9 and newer. We recommend using the latest stable release of Julia.
 
 ## For users

docs/src/preprocessing/preprocessing.md

@@ -1,10 +1,249 @@
+# Sampling of Geometries
 
-TODO
+Generating the initial configuration of a simulation requires filling volumes (3D) or surfaces (2D) of complex geometries with particles.
+The algorithm to sample a complex geometry should be robust and fast,
+since for large problems (large numbers of particles) or complex geometries (many geometry faces),
+generating the initial configuration is not trivial and can be very expensive in terms of computational cost.
+We therefore use a [winding number](https://en.wikipedia.org/wiki/Winding_number) approach for an inside-outside segmentation of an object.
+The winding number ``w(\mathbf{p})`` is a signed integer-valued function of a point ``\mathbf{p}`` and is defined as
 
+```math
+w(\mathbf{p}) = \frac{1}{2 \pi} \sum^n_{i=1} \Theta_i.
+```
+
+Here, ``\Theta_i`` is the *signed* angle between ``\mathbf{c}_i - \mathbf{p}`` and ``\mathbf{c}_{i+1} - \mathbf{p}`` where ``\mathbf{c}_i`` and ``\mathbf{c}_{i+1}`` are two consecutive vertices on a curve.
+In 3D, we refer to the solid angle of an *oriented* triangle with respect to ``\mathbf{p}``.
+
+We provide the following methods to calculate ``w(\mathbf{p})``:
+- Hormann et al. (2001) evaluate the winding number combined with an even-odd rule, but only for 2D polygons (see [WindingNumberHormann](@ref)).
+- Naive winding: Jacobson et al. (2013) generalized the winding number so that the algorithm can be applied for both 2D and 3D geometries (see [WindingNumberJacobson](@ref)).
+- Hierarchical winding: Jacobson et al. (2013) also introduced a fast hierarchical evaluation of the winding number. For further information see the description below.
+
+## [Hierarchical Winding](@id hierarchical_winding)
+According to Jacobson et al. (2013) the winding number with respect to a polygon (2D) or triangle mesh (3D) is the sum of the winding numbers with respect to each edge (2D) or face (3D).
+We can show this with the following example in which we determine the winding number for each edge of a triangle separately and sum them up:
+
+```julia
+using TrixiParticles
+using Plots
+
+triangle = [125.0 375.0 250.0 125.0;
+            175.0 175.0 350.0 175.0]
+
+# Delete all edges but one
+edge1 = deleteat!(TrixiParticles.Polygon(triangle), [2, 3])
+edge2 = deleteat!(TrixiParticles.Polygon(triangle), [1, 3])
+edge3 = deleteat!(TrixiParticles.Polygon(triangle), [1, 2])
+
+algorithm = WindingNumberJacobson()
+
+grid = hcat(([x, y] for x in 1:500, y in 1:500)...)
+
+_, w1 = algorithm(edge1, grid; store_winding_number=true)
+_, w2 = algorithm(edge2, grid; store_winding_number=true)
+_, w3 = algorithm(edge3, grid; store_winding_number=true)
+
+w = w1 + w2 + w3
+
+heatmap(1:500, 1:500, reshape(w1, 500, 500)', color=:coolwarm, showaxis=false,
+        tickfontsize=12, size=(570, 500), margin=6 * Plots.mm)
+heatmap(1:500, 1:500, reshape(w2, 500, 500)', color=:coolwarm, showaxis=false,
+        tickfontsize=12, size=(570, 500), margin=6 * Plots.mm)
+heatmap(1:500, 1:500, reshape(w3, 500, 500)', color=:coolwarm, showaxis=false,
+        tickfontsize=12, size=(570, 500), margin=6 * Plots.mm)
+heatmap(1:500, 1:500, reshape(w, 500, 500)', color=:coolwarm, showaxis=false,
+        tickfontsize=12, size=(570, 500), margin=6 * Plots.mm, clims=(-1, 1))
+
+```
+
+```@raw html
+<figure>
+  <img src="https://github.com/user-attachments/assets/bf491b2d-740e-4136-8a7b-e321f26f86fd" alt="triangle"/>
+</figure>
+```
+
+This summation property has some interesting consequences that we can utilize for an efficient computation of the winding number.
+Let ``\mathcal{S}`` be an open surface and ``\bar{\mathcal{S}}`` an arbitrary closing surface, such that
+
+```math
+\partial \bar{\mathcal{S}} = \partial \mathcal{S}
+```
+
+and ``\mathcal{B} = \bar{\mathcal{S}} \cup \mathcal{S}`` is some closed oriented surface.
+For any query point ``\mathbf{p}`` outside of ``\mathcal{B}``, we know that
+
+```math
+w_{\mathcal{S}}(\mathbf{p}) + w_{\bar{\mathcal{S}}}(\mathbf{p}) = w_{\mathcal{B}}(\mathbf{p}) = 0.
+```
+
+This means
+
+```math
+w_{\mathcal{S}}(\mathbf{p}) = - w_{\bar{\mathcal{S}}}(\mathbf{p}),
+```
+
+regardless of how ``\bar{\mathcal{S}}`` is constructed (as long as ``\mathbf{p}`` is outside of ``\mathcal{B}``).
+
+We can use this property in the discrete case to efficiently compute the winding number of a query point
+by partitioning the polygon or mesh in a "small" part (as in consisting of a small number of edges/faces) and a "large" part.
+For the small part we just compute the winding number, and for the large part we construct a small closing and compute its winding number.
+The partitioning is based on a hierarchical construction of bounding boxes.
+
+### Bounding volume hierarchy
+
+To efficiently find a "small part" and a "large part" as mentioned above, we construct a hierarchy of bounding boxes by starting with the whole domain and recursively splitting it in two equally sized boxes.
+The resulting hierarchy is a binary tree.
+
+The algorithm by Jacobsen et al. (Algorithm 2, p. 5) traverses this binary tree recursively until we find the leaf in which the query point is located.
+The recursion stops with the following criteria:
+
+- if the bounding box ``T`` is a leaf then ``T.\mathcal{S} = \mathcal{S} \cap T``, the part of ``\mathcal{S}``
+  that lies inside ``T``, is the "small part" mentioned above, so evaluate the winding number naively as ``w(\mathbf{p}, T.\mathcal{S})``.
+- else if ``\mathbf{p}`` is outside ``T`` then ``T.\mathcal{S}`` is the "large part", so evaluate the winding number naively
+  as ``-w(\mathbf{p}, T.\bar{\mathcal{S}})``, where ``T.\bar{\mathcal{S}}`` is the closing surface of ``T.\mathcal{S}``.
+
+#### Continuous example
+
+Now consider the following continuous (not discretized to a polygon) 2D example.
+We compute the winding number of the point ``\mathbf{p}`` with respect to ``\mathcal{S}`` using the depicted hierarchy of bounding boxes.
+
+```@raw html
+<figure>
+  <img src="https://github.com/user-attachments/assets/0ca2f475-6dd5-43f9-8b0c-87a0612ecdf4" alt="continuous closing"/>
+</figure>
+```
+
+(1):
+- Recurse left: ``w_{\text{left}} = \text{\texttt{hierarchical\_winding}} (\mathbf{p}, T.\text{left})``
+- Recurse right: ``w_{\text{right}} = \text{\texttt{hierarchical\_winding}} (\mathbf{p},T.\text{right})``
+
+(2):
+- Query point ``\mathbf{p}`` is outside bounding box ``T``, so don't recurse deeper.
+- Compute ``w_{\mathcal{S}}(\mathbf{p}) = - w_{\bar{\mathcal{S}}}(\mathbf{p})`` with the closure ``T.\bar{\mathcal{S}}``, which is generally much smaller (fewer edges in the discrete version) than ``T.\mathcal{S}``:
+
+```math
+w_{\text{left}} = -\text{\texttt{naive\_winding}} (\mathbf{p}, T.\bar{\mathcal{S}})
+```
+
+(3):
+- Bounding box ``T`` is a leaf. Use open surface ``T.\mathcal{S}``:
+
+```math
+w_{\text{right}} = \text{\texttt{naive\_winding}} (\mathbf{p}, T.\mathcal{S})
+```
+
+The reconstructed surface will then look as in the following image.
+
+```@raw html
+<figure>
+  <img src="https://github.com/user-attachments/assets/920bb4f1-1336-4e77-b06d-d5b46ca0d8d5" alt="reconstructed surface"/>
+</figure>
+```
+
+We finally sum up the winding numbers
+
+```math
+w = w_{\text{left}} + w_{\text{right} } = -w_{T_{\text{left}}.\bar{\mathcal{S}}} + w_{T_{\text{right}}.\mathcal{S}}
+```
+
+#### Discrete example
+
+We will now go through the discrete version of the example above.
+
+```@raw html
+<figure>
+  <img src="https://github.com/user-attachments/assets/a9b59cc3-5421-40af-b0b0-f4c18a5a7078" alt="discrete geometry"/>
+</figure>
+```
+
+To construct the hierarchy for the discrete piecewise-linear example in (1), we have to do the following.
+
+(2):
+Each edge is distributed to the child whose box contains the edge's barycenter (red dots in (2)).
+Splitting stops when the number of a box's edges slips below a
+threshold (usually ``\approx 100`` faces in 3D, here: 6 edges).
+
+(3):
+For the closure, Jacobson et al. (2013) define *exterior vertices* (*exterior edges* in 3D)
+as boundary vertices of such a segmentation (red dots in (3)).
+To find them, we traverse around each edge (face in 3D) in order, and
+increment or decrement for each vertex (edge) a specific counter.
+
+```julia
+v1 = edge_vertices_ids[edge][1]
+v2 = edge_vertices_ids[edge][2]
+
+vertex_count[v1] += 1
+vertex_count[v2] -= 1
+```
+
+In 2D, a vertex is declared as exterior if `vertex_count(vertex) != 0`, so there is not the same amount of edges in this box going into versus out of the vertex.
+To construct the closing surface, the exterior vertices are then connected to one arbitrary
+exterior vertex using appropriately oriented line segments:
+
+```julia
+edge = vertex_count[v] > 0 ? (closing_vertex, v) : (v, closing_vertex)
+```
+
+The resulting closed surface ``T.S \cup T.\bar{S}`` then has the same number of edges going into and out of each vertex.
+
+#### Incorrect evaluation
+
+If we follow the algorithm, we know that recursion stops if
+
+- the bounding box ``T`` is a leaf or
+- the query point ``\mathbf{p}`` is outside the box.
+
+```@raw html
+<figure>
+  <img src="https://github.com/user-attachments/assets/7bae164a-8d5b-4761-9d54-9abf99fca94a" alt="incorrect evaluation"/>
+</figure>
+```
+
+(1): The query point ``\mathbf{p}`` is outside the box, so we calculate the winding number with the (red) closure of the box.
+
+(2): The query point ``\mathbf{p}`` is inside the box, so we use the (blue) edges distributed to the box.
+
+(3): In this case, it leads to an incorrect evaluation of the winding number.
+The query point is clearly inside the box, but not inside the reconstructed surface.
+This is because the property ``w_{\mathcal{S}}(\mathbf{p}) = - w_{\bar{\mathcal{S}}}(\mathbf{p})``
+only holds when ``\mathbf{p}`` is outside of ``\mathcal{B}``, which is not the case here.
+
+#### Correct evaluation
+Jacobson et al. (2013) don't mention this problem or provide a solution to it.
+We contacted the authors and found that they know about this problem and solve it
+by resizing the bounding box to fully include the closing surface
+of the neighboring box, since it doesn't matter if the boxes overlap.
+
+```@raw html
+<figure>
+  <img src="https://github.com/user-attachments/assets/097f01f4-1f37-48e4-968a-4c0970548b24" alt="correct evaluation resizing"/>
+</figure>
+```
+
+To avoid resizing, we take a different approach and calculate the closure of the bounding box differently:
+- Exclude 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 flipped orientation,
+  so that the closing is actually a closing of the exterior plus intersecting edges.
+
+```@raw html
+<figure>
+  <img src="https://github.com/user-attachments/assets/a8ff9a7e-e6d6-44d1-9a29-7debddf2803d" alt="correct evaluation intersecting" width=60%/>
+</figure>
+```
+
+The evaluation then looks as follows.
+
+```@raw html
+<figure>
+  <img src="https://github.com/user-attachments/assets/9bb2d2ad-14e8-4bd0-a9bd-3c824932affd" alt="correct evaluation intersecting 2"/>
+</figure>
+```
 
 ```@autodocs
 Modules = [TrixiParticles]
-Pages = [joinpath("preprocessing", "point_in_poly", "winding_number_horman.jl")]
+Pages = [joinpath("preprocessing", "point_in_poly", "winding_number_hormann.jl")]
 ```
 
 ```@autodocs
@@ -17,11 +256,10 @@ Modules = [TrixiParticles]
 Pages = [joinpath("preprocessing", "geometries", "io.jl")]
 ```
 
-
 ### [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://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".
+- Kai Hormann, 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)

examples/fluid/pipe_flow_2d.jl

@@ -87,7 +87,7 @@ end
 inflow = InFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), flow_direction,
                 open_boundary_layers, density=fluid_density, particle_spacing)
 
-open_boundary_in = OpenBoundarySPHSystem(inflow; sound_speed, fluid_system,
+open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system,
                                          boundary_model=BoundaryModelLastiwka(),
                                          buffer_size=n_buffer_particles,
                                          reference_density=fluid_density,
@@ -98,7 +98,7 @@ outflow = OutFlow(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2
                   flow_direction, open_boundary_layers, density=fluid_density,
                   particle_spacing)
 
-open_boundary_out = OpenBoundarySPHSystem(outflow; sound_speed, fluid_system,
+open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system,
                                           boundary_model=BoundaryModelLastiwka(),
                                           buffer_size=n_buffer_particles,
                                           reference_density=fluid_density,

src/TrixiParticles.jl

@@ -73,7 +73,7 @@ export BoundaryMovement
 export examples_dir, validation_dir, trixi_include
 export trixi2vtk
 export RectangularTank, RectangularShape, SphereShape, ComplexShape
-export WindingNumberHorman, WindingNumberJacobson
+export WindingNumberHormann, WindingNumberJacobson
 export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry
 export SourceTermDamping
 export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection,

src/general/buffer.jl

@@ -1,10 +1,11 @@
 struct SystemBuffer{V}
-    active_particle :: BitVector
+    active_particle :: Vector{Bool}
     eachparticle    :: V # Vector{Int}
     buffer_size     :: Int
 
     function SystemBuffer(active_size, buffer_size::Integer)
-        active_particle = vcat(trues(active_size), falses(buffer_size))
+        # We cannot use a `BitVector` here, as writing to a `BitVector` is not thread-safe
+        active_particle = vcat(fill(true, active_size), fill(false, buffer_size))
         eachparticle = collect(1:active_size)
 
         return new{typeof(eachparticle)}(active_particle, eachparticle, buffer_size)

src/preprocessing/geometries/triangle_mesh.jl

@@ -194,7 +194,8 @@ function unique_sorted(vertices)
     # Sort by the first entry of the vectors
     compare_first_element = (x, y) -> x[1] < y[1]
     vertices_sorted = sort!(vertices, lt=compare_first_element)
-    keep = trues(length(vertices_sorted))
+    # We cannot use a `BitVector` here, as writing to a `BitVector` is not thread-safe
+    keep = fill(true, length(vertices_sorted))
 
     PointNeighbors.@threaded vertices_sorted for i in eachindex(vertices_sorted)
         # We only sorted by the first entry, so we have to check all previous vertices

src/preprocessing/point_in_poly/point_in_poly.jl

@@ -1,3 +1,3 @@
 include("hierarchical_winding.jl")
-include("winding_number_horman.jl")
+include("winding_number_hormann.jl")
 include("winding_number_jacobson.jl")

src/preprocessing/point_in_poly/winding_number_horman.jl → src/preprocessing/point_in_poly/winding_number_hormann.jl

@@ -1,20 +1,18 @@
 """
-    WindingNumberHorman()
+    WindingNumberHormann()
 
-Algorithm for inside-outside segmentation of a complex geometry proposed by [Horman et al. (2001)](@ref references_complex_shape).
+Algorithm for inside-outside segmentation of a complex geometry proposed by [Hormann et al. (2001)](@ref references_complex_shape).
 It is only supported for 2D geometries.
-[`WindingNumberHorman`](@ref) might handle edge cases a bit better, since the winding number is an integer value.
-Also, it is faster than [`WindingNumberJacobson`](@ref) for 2D geometries with more than about 100 edges.
-
+[`WindingNumberHormann`](@ref) might handle edge cases a bit better, since the winding number is an integer value.
 
 !!! warning "Experimental Implementation"
     This is an experimental feature and may change in any future releases.
 """
-struct WindingNumberHorman end
+struct WindingNumberHormann end
 
-# Algorithm 2 from Horman et al. (2001) "The point in polygon problem for arbitrary polygons"
+# Algorithm 2 from Hormann et al. (2001) "The point in polygon problem for arbitrary polygons"
 # https://doi.org/10.1016/S0925-7721(01)00012-8
-function (point_in_poly::WindingNumberHorman)(geometry, points; store_winding_number=false)
+function (point_in_poly::WindingNumberHormann)(geometry, points; store_winding_number=false)
     (; edge_vertices) = geometry
 
     # We cannot use a `BitVector` here, as writing to a `BitVector` is not thread-safe

src/preprocessing/point_in_poly/winding_number_jacobson.jl

@@ -58,7 +58,7 @@ Algorithm for inside-outside segmentation of a complex geometry proposed by Jaco
 - `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,
-                          which gives a significant speedup.
+                          which gives a significant speedup. For further information see [Hierarchical Winding](@ref hierarchical_winding).
 - `winding_number_factor`: For leaky geometries, a factor of `0.4` will give a better inside-outside segmentation.
 
 !!! warning "Experimental Implementation"