Merge branch 'main' into adami_pressure_changes

.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.22.9
+        uses: crate-ci/typos@v1.23.6

.gitignore

@@ -20,4 +20,9 @@ run
 *.vtu
 *.pvd
 *.mp4
-*.csv
+
+# ParaView states
+*.pvsm
+
+# FreeCAD models
+*.FCStd

.typos.toml

@@ -1,3 +1,6 @@
+[default.extend-identifiers]
+SPlisHSPlasH = "SPlisHSPlasH"
+
 [default.extend-words]
 ba = "ba"
 Shepard = "Shepard"

NEWS.md

@@ -2,23 +2,21 @@
 
 TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
-We aim at 3 to 4 month between major release versions and about 2 weeks between minor versions.
 
-## Version 0.3.x
+## Version 0.2.2
 
 ### Highlights
+Hotfix for threaded sampling of complex geometries.
 
-### Added
+## Version 0.2.1
 
-### Removed
-
-### Deprecated
+### Highlights
+Particle sampling of complex geometries from `.stl` and `.asc` files.
 
 ## Version 0.2.0
 
 ### Removed
-- Use of the internal neighborhood search has been removed and replaced with PointNeighbors.jl.
-
+Use of the internal neighborhood search has been removed and replaced with PointNeighbors.jl.
 
 ## Development Cycle 0.1
 

Project.toml

@@ -1,15 +1,17 @@
 name = "TrixiParticles"
 uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
 authors = ["erik.faulhaber <[email protected]>"]
-version = "0.2.1-dev"
+version = "0.2.3-dev"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 FastPow = "c0e83750-1142-43a8-81cf-6c956b72b4d1"
+FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
@@ -32,15 +34,17 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 Adapt = "3, 4"
 CSV = "0.10"
 DataFrames = "1.6"
+DelimitedFiles = "1"
 DiffEqCallbacks = "2.25, 3"
 FastPow = "0.1"
+FileIO = "1"
 ForwardDiff = "0.10"
 GPUArraysCore = "0.1"
 JSON = "0.21"
 KernelAbstractions = "0.9"
 MuladdMacro = "0.2"
 PointNeighbors = "0.4.2"
-Polyester = "0.7.5"
+Polyester = "0.7.10"
 RecipesBase = "1"
 Reexport = "1"
 SciMLBase = "1, 2"

README.md

@@ -21,6 +21,7 @@ Its features include:
 - Solid-body mechanics
   - Methods:  Total Lagrangian SPH (TLSPH), Discrete Element Method (DEM)
 - Fluid-Structure Interaction
+- Particle sampling of complex geometries from `.stl` and `.asc` files.
 - Output formats:
   - VTK
 

docs/make.jl

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

docs/src/preprocessing/preprocessing.md

@@ -0,0 +1,265 @@
+# Sampling of Geometries
+
+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})``:
+- Horman et al. (2001) evaluate the winding number combined with an even-odd rule, but only for 2D polygons (see [WindingNumberHorman](@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")]
+```
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("preprocessing", "point_in_poly", "winding_number_jacobson.jl")]
+```
+
+```@autodocs
+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".
+  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)