Merge branch 'main' into open-boundaries-nhs-fix-patch-2

Author: LasNikas

docs/src/preprocessing/preprocessing.md

@@ -1,4 +1,4 @@
-# Sampling of Geometries
+# [Sampling of Geometries](@id 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,
@@ -256,9 +256,57 @@ Modules = [TrixiParticles]
 Pages = [joinpath("preprocessing", "geometries", "io.jl")]
 ```
 
-# Particle Packing
+# [Particle Packing](@id particle_packing)
+To obtain a body-fitted and isotropic particle distribution, an initial configuration (see [Sampling of Geometries](@ref sampling_of_geometries)) is first generated. This configuration is then packed using a [`ParticlePackingSystem`](@ref).
+The preprocessing pipeline consists of the following steps:
+
+- Load geometry: Fig. 1, [`load_geometry`](@ref).
+- Compute the signed distance field (SDF): Fig. 2, [`SignedDistanceField`](@ref).
+- Generate boundary particles: Fig. 3, [`sample_boundary`](@ref).
+- Initial sampling of the interior particles with inside-outside segmentation: Fig. 4, [`ComplexShape`](@ref).
+- Pack the initial configuration of interior and boundary particles (Fig. 5): Fig. 6, [`ParticlePackingSystem`](@ref).
+
+The input data can either be a 3D triangulated surface mesh represented in STL format or a 2D polygonal traversal of the geometry (see [`load_geometry`](@ref)).
+The second step involves generating the SDF (see [`SignedDistanceField`](@ref)), which is necessary for the final packing step as it requires a surface detection.
+The SDF is illustrated in Fig. 2, where the distances to the surface of the geometry are visualized as a color map.
+As shown, the SDF is computed only within a narrow band around the geometry’s surface, enabling  a face-based neighborhood search (NHS) to be used exclusively during this step.
+In the third step, the initial configuration of the boundary particles is generated (orange particles in Fig. 3).
+Boundary particles are created by copying the positions of SDF points located outside the geometry but within a predefined boundary thickness (see [`sample_boundary`](@ref)).
+In the fourth step, the initial configuration of the interior particles (green particles in Fig. 4) is generated using the hierarchical winding number approach (see [Hierarchical Winding](@ref hierarchical_winding)).
+After steps **1** through **4**, the initial configuration of both interior and boundary particles is obtained, as illustrated in Fig. 5.
+The interface of the geometry surface is not well resolved with the initial particle configuration.
+Thus, in the final step, a packing algorithm by Zhu et al. [Zhu2021](@cite) is applied utilizing the SDF to simultaneously optimize the positions of both interior and boundary particles,
+yielding an isotropic distribution while accurately preserving the geometry surface, as illustrated in Fig. 6.
+
+```@raw html
+<div style="display: flex; gap: 16px; flex-wrap: wrap;">
+  <figure style="margin: 0; text-align: center;">
+    <img src="https://github.com/user-attachments/assets/7fe9d1f1-1633-4377-8b97-a4d1778aee07" alt="geometry" style="max-width: 200px;">
+    <figcaption>(1) Geometry representation</figcaption>
+  </figure>
+  <figure style="margin: 0; text-align: center;">
+    <img src="https://github.com/user-attachments/assets/2b79188c-3148-49f1-8337-718721851bf5" alt="sdf" style="max-width: 200px;">
+    <figcaption>(2) Signed distances to the surface</figcaption>
+  </figure>
+  <figure style="margin: 0; text-align: center;">
+    <img src="https://github.com/user-attachments/assets/1501718f-d1f5-4f14-b1bc-2a2e581db476" alt="boundary" style="max-width: 200px;">
+    <figcaption>(3) Boundary particles</figcaption>
+  </figure>
+  <figure style="margin: 0; text-align: center;">
+    <img src="https://github.com/user-attachments/assets/f7376b15-324a-4da1-bb59-db01c7bd6620" alt="interior" style="max-width: 200px;">
+    <figcaption>(4) Interior particles</figcaption>
+  </figure>
+  <figure style="margin: 0; text-align: center;">
+    <img src="https://github.com/user-attachments/assets/4be889d6-e70a-4c5e-bef2-0071ea4d898c" alt="initial_config" style="max-width: 200px;">
+    <figcaption>(5) Initial configuration</figcaption>
+  </figure>
+  <figure style="margin: 0; text-align: center;">
+    <img src="https://github.com/user-attachments/assets/0f7aba29-3cf7-4ec1-8c95-841e72fe620d" alt="packed_config" style="max-width: 200px;">
+    <figcaption>(6) Packed configuration</figcaption>
+  </figure>
+</div>
+```
 
-TODO
 
 ```@autodocs
 Modules = [TrixiParticles]

docs/src/refs.bib

@@ -735,3 +735,17 @@ @book{Poling2001
   publisher = {McGraw-Hill},
   address = {New York}
 }
+
+@article{Zhu2021,
+  author    = {Zhu, Yujie and Zhang, Chi and Yu, Yongchuan and Hu, Xiangyu},
+  journal   = {Journal of Hydrodynamics},
+  title     = {A {CAD}-compatible body-fitted particle generator for arbitrarily complex geometry and its application to wave-structure interaction},
+  year      = {2021},
+  issn      = {1878-0342},
+  month     = apr,
+  number    = {2},
+  pages     = {195--206},
+  volume    = {33},
+  doi       = {10.1007/s42241-021-0031-y},
+  publisher = {Springer Science and Business Media LLC}
+}

src/preprocessing/particle_packing/system.jl

@@ -1,15 +1,15 @@
 """
     ParticlePackingSystem(shape::InitialCondition;
-                          signed_distance_field::SignedDistanceField,
+                          signed_distance_field::Union{SignedDistanceField, Nothing},
                           smoothing_kernel=SchoenbergQuinticSplineKernel{ndims(shape)}(),
                           smoothing_length=shape.particle_spacing,
                           smoothing_length_interpolation=smoothing_length,
                           is_boundary=false, boundary_compress_factor=1,
                           neighborhood_search=GridNeighborhoodSearch{ndims(shape)}(),
-                          background_pressure, tlsph=true, fixed_system=false)
+                          background_pressure, tlsph=false, fixed_system=false)
 
 System to generate body-fitted particles for complex shapes.
-For more information on the methods, see description below.
+For more information on the methods, see [particle packing](@ref particle_packing).
 
 # Arguments
 - `shape`: [`InitialCondition`](@ref) to be packed.
@@ -23,7 +23,7 @@ For more information on the methods, see description below.
                            as for fluids. When `tlsph=true`, particles will be placed
                            on the boundary of the shape.
 - `is_boundary`:           When `shape` is inside the geometry that was used to create
-                           `signed_distance_field, set `is_boundary=false`.
+                           `signed_distance_field`, set `is_boundary=false`.
                            Otherwise (`shape` is the sampled boundary), set `is_boundary=true`.
                            The thickness of the boundary is specified by creating
                            `signed_distance_field` with:
@@ -40,12 +40,19 @@ For more information on the methods, see description below.
                            the signed distance from a particle to a face is required.
                            The precalculated signed distances will be interpolated
                            to each particle during the packing procedure.
+                           Set `signed_distance_field=nothing` when packing with a fixed system
+                           (see `fixed_system` description above).
 - `smoothing_kernel`:      Smoothing kernel to be used for this system.
                            See [Smoothing Kernels](@ref smoothing_kernel).
 - `smoothing_length`:      Smoothing length to be used for the gradient estimation.
                            See [Smoothing Kernels](@ref smoothing_kernel).
 - `smoothing_length_interpolation`: Smoothing length to be used for interpolating the `SignedDistanceField` information.
                                     The default is `smoothing_length_interpolation = smoothing_length`.
+- `boundary_compress_factor`: Factor to compress the boundary particles by reducing the boundary thickness by a factor of `boundary_compress_factor`.
+                              The default value is `1`, which means no compression.
+                              Compression can be useful for highly convex geometries,
+                              where the boundary volume increases significantly while the mass of the boundary particles remains constant.
+                              Recommended values are `0.8` or `0.9`.
 """
 struct ParticlePackingSystem{S, F, NDIMS, ELTYPE <: Real, PR, C, AV,
                              IC, M, D, K, N, SD, UCU} <: FluidSystem{NDIMS}
@@ -95,8 +102,7 @@ struct ParticlePackingSystem{S, F, NDIMS, ELTYPE <: Real, PR, C, AV,
 end
 
 function ParticlePackingSystem(shape::InitialCondition;
-                               signed_distance_field::Union{SignedDistanceField,
-                                                            Nothing},
+                               signed_distance_field::Union{SignedDistanceField, Nothing},
                                smoothing_kernel=SchoenbergQuinticSplineKernel{ndims(shape)}(),
                                smoothing_length=shape.particle_spacing,
                                smoothing_length_interpolation=smoothing_length,

src/schemes/fluid/entropically_damped_sph/system.jl

@@ -306,6 +306,10 @@ function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode,
     compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
     compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
     update_average_pressure!(system, system.transport_velocity, v_ode, u_ode, semi)
+
+    # Check that TVF is only used together with `UpdateCallback`
+    check_tvf_configuration(system, system.transport_velocity, v, u, v_ode, u_ode, semi, t;
+                            update_from_callback)
 end
 
 function update_average_pressure!(system, ::Nothing, v_ode, u_ode, semi)

src/schemes/fluid/transport_velocity.jl

@@ -165,19 +165,15 @@ function update_callback_used!(system, transport_velocity)
     system.cache.update_callback_used[] = true
 end
 
-function update_final!(system::FluidSystem, v, u, v_ode, u_ode, semi, t;
-                       update_from_callback=false)
-    update_final!(system, system.transport_velocity,
-                  v, u, v_ode, u_ode, semi, t; update_from_callback)
-end
-
-function update_final!(system::FluidSystem, ::Nothing,
-                       v, u, v_ode, u_ode, semi, t; update_from_callback=false)
+function check_tvf_configuration(system::FluidSystem, ::Nothing,
+                                 v, u, v_ode, u_ode, semi, t; update_from_callback=false)
     return system
 end
 
-function update_final!(system::FluidSystem, ::TransportVelocityAdami,
-                       v, u, v_ode, u_ode, semi, t; update_from_callback=false)
+function check_tvf_configuration(system::FluidSystem, ::TransportVelocityAdami,
+                                 v, u, v_ode, u_ode, semi, t; update_from_callback=false)
+    # The `UpdateCallback` will set `system.cache.update_callback_used[]` to `true`
+    # in the initialization. However, other callbacks might update the system first, hence `update_from_callback`.
     if !update_from_callback && !(system.cache.update_callback_used[])
         throw(ArgumentError("`UpdateCallback` is required when using `TransportVelocityAdami`"))
     end

src/schemes/fluid/weakly_compressible_sph/system.jl

@@ -304,6 +304,10 @@ function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode,
     # Surface normal of neighbor and boundary needs to have been calculated already
     compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
     compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
+
+    # Check that TVF is only used together with `UpdateCallback`
+    check_tvf_configuration(system, system.transport_velocity, v, u, v_ode, u_ode, semi, t;
+                            update_from_callback)
 end
 
 function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode,