diff --git a/docs/literate/src/tut_packing.jl b/docs/literate/src/tut_packing.jl index d5b9610283..d4b55990ed 100644 --- a/docs/literate/src/tut_packing.jl +++ b/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) diff --git a/docs/src/preprocessing/preprocessing.md b/docs/src/preprocessing/preprocessing.md index 0518574b73..e08b46c7c7 100644 --- a/docs/src/preprocessing/preprocessing.md +++ b/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. diff --git a/examples/preprocessing/packing_2d.jl b/examples/preprocessing/packing_2d.jl index 02951323ef..cba529cc5d 100644 --- a/examples/preprocessing/packing_2d.jl +++ b/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) # ========================================================================================== diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index bc72e155b3..ef2fe87cd4 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/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 diff --git a/src/setups/complex_shape.jl b/src/setups/complex_shape.jl index 6a78b412e2..481f80a74c 100644 --- a/src/setups/complex_shape.jl +++ b/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 diff --git a/test/setups/complex_shape.jl b/test/setups/complex_shape.jl index 2917661225..c2422e23cd 100644 --- a/test/setups/complex_shape.jl +++ b/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 = [ diff --git a/test/systems/packing_system.jl b/test/systems/packing_system.jl index 1eb02e9c44..f7b5e31cf4 100644 --- a/test/systems/packing_system.jl +++ b/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)