Merge branch 'main' into boundary-zone-from-stl-file

.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/accelerated_tank_2d.jl

@@ -4,63 +4,17 @@
 using TrixiParticles
 using OrdinaryDiffEq
 
-# ==========================================================================================
-# ==== Resolution
+# Resolution
 fluid_particle_spacing = 0.05
 
-# Make sure that the kernel support of fluid particles at a boundary is always fully sampled
-boundary_layers = 3
-
-# ==========================================================================================
-# ==== Experiment Setup
-tspan = (0.0, 1.0)
-
-# Boundary geometry and initial fluid particle positions
-initial_fluid_size = (1.0, 0.9)
-tank_size = (1.0, 1.0)
-
-fluid_density = 1000.0
-sound_speed = 10.0
-state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
-                                   exponent=7, clip_negative_pressure=false)
-
-tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density,
-                       acceleration=(0.0, -9.81), state_equation=state_equation,
-                       n_layers=boundary_layers, spacing_ratio=1.0)
-
 # Function for moving boundaries
 movement_function(t) = SVector(0.0, 0.5 * 9.81 * t^2)
 
 is_moving(t) = true
 
 boundary_movement = BoundaryMovement(movement_function, is_moving)
 
-# Import the setup from `hydrostatic_water_column_2d.jl`
 trixi_include(@__MODULE__,
               joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"),
               fluid_particle_spacing=fluid_particle_spacing, movement=boundary_movement,
-              acceleration=(0.0, 0.0), tank=tank, semi=nothing, ode=nothing,
-              sol=nothing) # Overwrite `sol` assignment to skip time integration
-
-# ==========================================================================================
-# ==== Simulation
-semi = Semidiscretization(fluid_system, boundary_system)
-ode = semidiscretize(semi, tspan)
-
-info_callback = InfoCallback(interval=50)
-saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
-
-callbacks = CallbackSet(info_callback, saving_callback)
-
-# Use a Runge-Kutta method with automatic (error based) time step size control.
-# Limiting of the maximum stepsize is necessary to prevent crashing.
-# When particles are approaching a wall in a uniform way, they can be advanced
-# with large time steps. Close to the wall, the stepsize has to be reduced drastically.
-# Sometimes, the method fails to do so because forces become extremely large when
-# fluid particles are very close to boundary particles, and the time integration method
-# interprets this as an instability.
-sol = solve(ode, RDPK3SpFSAL35(),
-            abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
-            reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
-            dtmax=1e-2, # Limit stepsize to prevent crashing
-            save_everystep=false, callback=callbacks);
+              tspan=(0.0, 1.0), system_acceleration=(0.0, 0.0))