Fix complex boundary sampling bounds

docs/literate/src/tut_packing.jl

@@ -75,8 +75,8 @@ plot!(right_margin=5Plots.mm) #hide
 # ## Creating an initial configuration of boundary particles
 
 # To create the initial configuration of the boundary particles,
-# we use the sampled points of the SDF whose signed distance lies between 0
-# and `boundary_thickness`.
+# we use the sampled points of the SDF whose signed distance lies between the
+# geometry offset implied by `place_on_shell` and `boundary_thickness`.
 # Here, we need to specify the `density` of the boundary particles.
 # As an example, we choose `1.0` for all particles.
 # This gives us an [`InitialCondition`](@ref InitialCondition) for the boundary particles.
@@ -211,7 +211,7 @@ plot!(geometry, seriestype=:path, color=:black, label=nothing, linewidth=2)
 boundary_system = ParticlePackingSystem(boundary_sampled; is_boundary=true,
                                         smoothing_kernel, smoothing_length,
                                         boundary_compress_factor=0.7, signed_distance_field,
-                                        background_pressure)
+                                        boundary_thickness, background_pressure)
 
 # We can now couple the boundary system with the interior system:
 semi = Semidiscretization(packing_system, boundary_system)

docs/src/preprocessing/preprocessing.md

@@ -325,7 +325,7 @@ The second step involves generating the SDF (see [`SignedDistanceField`](@ref)),
 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)).
+Boundary particles are created by copying the positions of SDF points located outside the geometry, starting at the offset implied by `place_on_shell` and ending at 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.

examples/preprocessing/packing_2d.jl

@@ -71,7 +71,8 @@ packing_system = ParticlePackingSystem(shape_sampled; smoothing_length,
 
 boundary_system = ParticlePackingSystem(boundary_sampled; smoothing_length,
                                         is_boundary=true, signed_distance_field,
-                                        place_on_shell, boundary_compress_factor=0.8,
+                                        place_on_shell, boundary_thickness,
+                                        boundary_compress_factor=0.8,
                                         background_pressure)
 
 # ==========================================================================================

src/preprocessing/particle_packing/system.jl

@@ -5,6 +5,7 @@
                           smoothing_length=shape.particle_spacing,
                           smoothing_length_interpolation=smoothing_length,
                           is_boundary=false, boundary_compress_factor=1,
+                          boundary_thickness=nothing,
                           neighborhood_search=GridNeighborhoodSearch{ndims(shape)}(),
                           background_pressure, place_on_shell=false, fixed_system=false)
 
@@ -26,10 +27,6 @@ For more information on the methods, see [particle packing](@ref particle_packin
 - `is_boundary`:           When `shape` is inside the geometry that was used to create
                            `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:
-                              - `use_for_boundary_packing=true`
-                              - `max_signed_distance=boundary_thickness`
                            See [`SignedDistanceField`](@ref).
 - `fixed_system`:          When set to `true`, the system remains static, meaning particles
                            will not move and the `InitialCondition` will stay unchanged.
@@ -54,6 +51,10 @@ For more information on the methods, see [particle packing](@ref particle_packin
                               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`.
+- `boundary_thickness`: Thickness of the sampled boundary when `is_boundary=true`.
+                        By default, this is `signed_distance_field.max_signed_distance`.
+                        If [`sample_boundary`](@ref) used a smaller `boundary_thickness`
+                        than the `SignedDistanceField`, pass the same value here.
 """
 struct ParticlePackingSystem{S, F, NDIMS, ELTYPE <: Real, PR, C, AV,
                              IC, M, D, K, N, SD} <: AbstractFluidSystem{NDIMS}
@@ -106,6 +107,7 @@ function ParticlePackingSystem(shape::InitialCondition;
                                smoothing_length=shape.particle_spacing,
                                smoothing_length_interpolation=smoothing_length,
                                is_boundary=false, boundary_compress_factor=1,
+                               boundary_thickness=nothing,
                                neighborhood_search=GridNeighborhoodSearch{ndims(shape)}(),
                                background_pressure, place_on_shell=false,
                                fixed_system=false)
@@ -147,10 +149,30 @@ function ParticlePackingSystem(shape::InitialCondition;
     # Its value is negative if the particle is inside the geometry.
     # Otherwise (if outside), the value is positive.
     if is_boundary
-        offset = place_on_shell ? shape.particle_spacing : shape.particle_spacing / 2
+        if isnothing(signed_distance_field)
+            fixed_system ||
+                throw(ArgumentError("`signed_distance_field` is required when `is_boundary=true`"))
+
+            shift_length = zero(ELTYPE)
+        else
+            boundary_thickness_ = isnothing(boundary_thickness) ?
+                                  signed_distance_field.max_signed_distance :
+                                  convert(ELTYPE, boundary_thickness)
+
+            if boundary_thickness_ > signed_distance_field.max_signed_distance
+                throw(ArgumentError("`boundary_thickness` is greater than " *
+                                    "`max_signed_distance` of `SignedDistanceField`."))
+            end
+
+            if boundary_thickness_ < zero(boundary_thickness_)
+                throw(ArgumentError("`boundary_thickness` must be non-negative"))
+            end
 
-        shift_length = -boundary_compress_factor *
-                       signed_distance_field.max_signed_distance - offset
+            offset = place_on_shell ? shape.particle_spacing : shape.particle_spacing / 2
+
+            shift_length = -boundary_compress_factor *
+                           boundary_thickness_ - offset
+        end
     else
         shift_length = place_on_shell ? zero(ELTYPE) : shape.particle_spacing / 2
     end

src/setups/complex_shape.jl

@@ -91,8 +91,10 @@ of the geometry.
 - `boundary_density`: Density of each boundary particle.
 - `place_on_shell`:     When `place_on_shell=true`, boundary particles will be placed
                         one particle spacing from the surface of the geometry.
-                        Otherwise when `place_on_shell=true` (simulating fluid particles),
+                        Otherwise when `place_on_shell=false` (simulating fluid particles),
                         boundary particles will be placed half particle spacing away from the surface.
+                        Thus, `boundary_thickness` must be at least one particle spacing
+                        for `place_on_shell=true` and half a particle spacing otherwise.
 
 
 # Examples
@@ -111,7 +113,7 @@ boundary_sampled = sample_boundary(signed_distance_field; boundary_density=1.0,
 │ InitialCondition                                                                                 │
 │ ════════════════                                                                                 │
 │ #dimensions: ……………………………………………… 2                                                                │
-│ #particles: ………………………………………………… 889                                                              │
+│ #particles: ………………………………………………… 677                                                              │
 │ particle spacing: ………………………………… 0.03                                                             │
 │ eltype: …………………………………………………………… Float64                                                          │
 │ coordinate eltype: ……………………………… Float64                                                          │
@@ -133,10 +135,22 @@ function sample_boundary(signed_distance_field;
     end
 
     # Only keep the required part of the signed distance field
-    distance_to_boundary = zero(particle_spacing)
-    keep_indices = (distance_to_boundary .< distances .<= max_signed_distance)
+    distance_to_boundary = place_on_shell ? particle_spacing : particle_spacing / 2
+    if boundary_thickness < distance_to_boundary
+        throw(ArgumentError("`boundary_thickness` must be at least " *
+                            "`particle_spacing` for `place_on_shell=true` and " *
+                            "half `particle_spacing` for `place_on_shell=false`."))
+    end
+
+    keep_indices = (distance_to_boundary .<= distances .<= boundary_thickness)
+    boundary_positions = positions[keep_indices]
+
+    if isempty(boundary_positions)
+        throw(ArgumentError("No boundary particles were sampled. Increase " *
+                            "`boundary_thickness` or generate a denser `SignedDistanceField`."))
+    end
 
-    boundary_coordinates = stack(positions[keep_indices])
+    boundary_coordinates = stack(boundary_positions)
     return InitialCondition(; coordinates=boundary_coordinates, density=boundary_density,
                             particle_spacing)
 end

test/setups/complex_shape.jl

@@ -2,6 +2,44 @@
     data_dir = pkgdir(TrixiParticles, "examples", "preprocessing", "data")
     validation_dir = pkgdir(TrixiParticles, "test", "preprocessing", "data")
 
+    @testset verbose=true "Sample Boundary" begin
+        particle_spacing = 0.1
+        positions = [
+            SVector(0.0, 0.0),
+            SVector(0.1, 0.0),
+            SVector(0.2, 0.0),
+            SVector(0.3, 0.0),
+            SVector(0.4, 0.0)
+        ]
+        distances = [0.01, 0.05, 0.1, 0.2, 0.3]
+
+        signed_distance_field = (; positions, distances, particle_spacing,
+                                 boundary_packing=true, max_signed_distance=0.3)
+
+        boundary = sample_boundary(signed_distance_field; boundary_density=1.0,
+                                   boundary_thickness=0.2, place_on_shell=false)
+        @test boundary.coordinates ≈ stack(positions[2:4])
+
+        boundary = sample_boundary(signed_distance_field; boundary_density=1.0,
+                                   boundary_thickness=0.2, place_on_shell=true)
+        @test boundary.coordinates ≈ stack(positions[3:4])
+
+        @test_throws ArgumentError sample_boundary(signed_distance_field;
+                                                   boundary_density=1.0,
+                                                   boundary_thickness=0.04,
+                                                   place_on_shell=false)
+
+        too_thin_sdf = (; positions, distances, particle_spacing,
+                        boundary_packing=true, max_signed_distance=0.1)
+        @test_throws ArgumentError sample_boundary(too_thin_sdf; boundary_density=1.0,
+                                                   boundary_thickness=0.2)
+
+        not_boundary_sdf = (; positions, distances, particle_spacing,
+                            boundary_packing=false, max_signed_distance=0.3)
+        @test_throws ArgumentError sample_boundary(not_boundary_sdf; boundary_density=1.0,
+                                                   boundary_thickness=0.2)
+    end
+
     @testset verbose=true "2D" begin
         @testset verbose=true "Shifted Rectangle" begin
             algorithms = [

test/systems/packing_system.jl

@@ -41,6 +41,22 @@
         └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""
         @test repr("text/plain", system) == show_box
 
+        signed_distance_field = SignedDistanceField(geometry, 0.1;
+                                                    use_for_boundary_packing=true,
+                                                    max_signed_distance=0.3)
+        boundary_sampled = sample_boundary(signed_distance_field; boundary_density=1.0,
+                                           boundary_thickness=0.2,
+                                           place_on_shell=false)
+        system = ParticlePackingSystem(boundary_sampled; signed_distance_field,
+                                       background_pressure=1.0, is_boundary=true,
+                                       boundary_thickness=0.2)
+        @test system.shift_length == -0.25
+        @test_throws ArgumentError ParticlePackingSystem(boundary_sampled;
+                                                         signed_distance_field,
+                                                         background_pressure=1.0,
+                                                         is_boundary=true,
+                                                         boundary_thickness=0.4)
+
         system = ParticlePackingSystem(initial_condition,
                                        signed_distance_field=nothing,
                                        background_pressure=1.0)