From 20c9e72fcd5fe0f8f0dc78a93fcd7491b63602dc Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 30 Sep 2025 17:53:24 +0200 Subject: [PATCH 01/20] add automatic reordering --- src/general/initial_condition.jl | 13 ++++++ .../structure/total_lagrangian_sph/system.jl | 46 +++++++++++++++---- 2 files changed, 49 insertions(+), 10 deletions(-) diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index 2e7ac773b8..475f844223 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -415,3 +415,16 @@ function find_too_close_particles(coords, min_distance) return result end + +function move_particles_to_end!(ic::InitialCondition, particle_ids_to_move) + sort_key = [i in particle_ids_to_move ? 1 : 0 for i in eachparticle(ic)] + permutation = sortperm(sort_key) + + ic.coordinates .= ic.coordinates[:, permutation] + ic.velocity .= ic.velocity[:, permutation] + ic.mass .= ic.mass[permutation] + ic.density .= ic.density[permutation] + ic.pressure .= ic.pressure[permutation] + + return ic +end diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index 8d69638256..fd53f3ca50 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -1,10 +1,11 @@ @doc raw""" TotalLagrangianSPHSystem(initial_condition, - smoothing_kernel, smoothing_length, - young_modulus, poisson_ratio; - n_clamped_particles=0, boundary_model=nothing, - acceleration=ntuple(_ -> 0.0, NDIMS), - penalty_force=nothing, source_terms=nothing) + smoothing_kernel, smoothing_length, + young_modulus, poisson_ratio; + n_clamped_particles=0, clamped_particles=Int[], + acceleration=ntuple(_ -> 0.0, NDIMS), + penalty_force=nothing, viscosity=nothing, + source_terms=nothing, boundary_model=nothing) System for particles of an elastic structure. @@ -13,6 +14,16 @@ all relevant quantities and operators are measured with respect to the initial configuration (O’Connor & Rogers 2021, Belytschko et al. 2000). See [Total Lagrangian SPH](@ref tlsph) for more details on the method. +There are two approaches to specify clamped particles in the system: + +1. **Manual ordering**: Provide the number of clamped particles via `n_clamped_particles`, + ensuring that these particles are the last entries in the `InitialCondition` + (see the info box below for details). + +2. **Automatic ordering**: Pass the indices of the particles to be clamped via `clamped_particles`. + In this case, the specified particles will be internally reordered so that they become + the last particles in the `InitialCondition`. + # Arguments - `initial_condition`: Initial condition representing the system's particles. - `young_modulus`: Young's modulus. @@ -24,8 +35,10 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method. # Keyword Arguments - `n_clamped_particles`: Number of clamped particles which are fixed and not integrated - to clamp the structure. Note that the clamped particles must be the **last** - particles in the `InitialCondition`. See the info box below. + to clamp the structure. Note that the clamped particles must be the **last** + particles in the `InitialCondition`. See the info box below. +- `clamped_particles`: Indices of the clamped particles which are fixed and not integrated + to clamp the structure. - `boundary_model`: Boundary model to compute the hydrodynamic density and pressure for fluid-structure interaction (see [Boundary Models](@ref boundary_models)). - `penalty_force`: Penalty force to ensure regular particle position under large deformations @@ -40,7 +53,8 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method. See, for example, [`SourceTermDamping`](@ref). !!! note - The clamped particles must be the **last** particles in the `InitialCondition`. + If specifing the clamped particles manually (via `n_clamped_particles`), + the clamped particles must be the **last** particles in the `InitialCondition`. To do so, e.g. use the `union` function: ```jldoctest; output = false, setup = :(clamped_particles = RectangularShape(0.1, (1, 4), (0.0, 0.0), density=1.0); beam = RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0)) structure = union(beam, clamped_particles) @@ -84,11 +98,11 @@ end function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length, young_modulus, poisson_ratio; - n_clamped_particles=0, boundary_model=nothing, + n_clamped_particles=0, clamped_particles=Int[], acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), penalty_force=nothing, viscosity=nothing, - source_terms=nothing) + source_terms=nothing, boundary_model=nothing) NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) n_particles = nparticles(initial_condition) @@ -103,6 +117,18 @@ function TotalLagrangianSPHSystem(initial_condition, throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end + # Handle clamped particles: + # If only n_clamped_particles > 0: assume manual ordering (do nothing). + # If both are 0/empty: no clamped particles (do nothing). + if !isempty(clamped_particles) + # Check consistency if both methods are used + if n_clamped_particles > 0 && n_clamped_particles != length(clamped_particles) + throw(ArgumentError("`n_clamped_particles` must equal length of `clamped_particles` when both are specified")) + end + n_clamped_particles = length(clamped_particles) + move_particles_to_end!(initial_condition, clamped_particles) + end + initial_coordinates = copy(initial_condition.coordinates) current_coordinates = copy(initial_condition.coordinates) mass = copy(initial_condition.mass) From da82198a6d8436b5c95c97894264e6f37bbad3bd Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 30 Sep 2025 17:59:33 +0200 Subject: [PATCH 02/20] fix typo --- src/schemes/structure/total_lagrangian_sph/system.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index fd53f3ca50..92f311685e 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -53,7 +53,7 @@ There are two approaches to specify clamped particles in the system: See, for example, [`SourceTermDamping`](@ref). !!! note - If specifing the clamped particles manually (via `n_clamped_particles`), + If specifying the clamped particles manually (via `n_clamped_particles`), the clamped particles must be the **last** particles in the `InitialCondition`. To do so, e.g. use the `union` function: ```jldoctest; output = false, setup = :(clamped_particles = RectangularShape(0.1, (1, 4), (0.0, 0.0), density=1.0); beam = RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0)) From 51affb8d399cd3eb2b06f0f7852f043230595fd9 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 30 Sep 2025 18:02:54 +0200 Subject: [PATCH 03/20] revise docs --- src/schemes/structure/total_lagrangian_sph/system.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index 92f311685e..83490a892a 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -2,7 +2,8 @@ TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length, young_modulus, poisson_ratio; - n_clamped_particles=0, clamped_particles=Int[], + n_clamped_particles=0, + clamped_particles::Vector{Int}=Int[], acceleration=ntuple(_ -> 0.0, NDIMS), penalty_force=nothing, viscosity=nothing, source_terms=nothing, boundary_model=nothing) @@ -37,8 +38,8 @@ There are two approaches to specify clamped particles in the system: - `n_clamped_particles`: Number of clamped particles which are fixed and not integrated to clamp the structure. Note that the clamped particles must be the **last** particles in the `InitialCondition`. See the info box below. -- `clamped_particles`: Indices of the clamped particles which are fixed and not integrated - to clamp the structure. +- `clamped_particles`: A vector of indices specifying the clamped particles which are fixed + and not integrated to clamp the structure. - `boundary_model`: Boundary model to compute the hydrodynamic density and pressure for fluid-structure interaction (see [Boundary Models](@ref boundary_models)). - `penalty_force`: Penalty force to ensure regular particle position under large deformations @@ -98,7 +99,8 @@ end function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length, young_modulus, poisson_ratio; - n_clamped_particles=0, clamped_particles=Int[], + n_clamped_particles=0, + clamped_particles::Vector{Int}=Int[], acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), penalty_force=nothing, viscosity=nothing, From 270dd86606bacdbccd9512b8d162b779d272b783 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 30 Sep 2025 18:05:02 +0200 Subject: [PATCH 04/20] indentation --- .../structure/total_lagrangian_sph/system.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index 83490a892a..a3f430dac8 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -1,12 +1,12 @@ @doc raw""" TotalLagrangianSPHSystem(initial_condition, - smoothing_kernel, smoothing_length, - young_modulus, poisson_ratio; - n_clamped_particles=0, - clamped_particles::Vector{Int}=Int[], - acceleration=ntuple(_ -> 0.0, NDIMS), - penalty_force=nothing, viscosity=nothing, - source_terms=nothing, boundary_model=nothing) + smoothing_kernel, smoothing_length, + young_modulus, poisson_ratio; + n_clamped_particles=0, + clamped_particles::Vector{Int}=Int[], + acceleration=ntuple(_ -> 0.0, NDIMS), + penalty_force=nothing, viscosity=nothing, + source_terms=nothing, boundary_model=nothing) System for particles of an elastic structure. From 1ad1a2be8573fe983819f991e2c38f5307c89945 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 2 Oct 2025 14:30:57 +0200 Subject: [PATCH 05/20] add tests --- src/general/initial_condition.jl | 6 ++- test/general/initial_condition.jl | 63 +++++++++++++++++++++++++++++++ 2 files changed, 68 insertions(+), 1 deletion(-) diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index 475f844223..7d2a03e2a1 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -416,8 +416,12 @@ function find_too_close_particles(coords, min_distance) return result end -function move_particles_to_end!(ic::InitialCondition, particle_ids_to_move) +function move_particles_to_end!(ic::InitialCondition, particle_ids_to_move::Vector) + invalid_id = findfirst(i -> i <= 0 || i > nparticles(ic), particle_ids_to_move) + isnothing(invalid_id) || throw(BoundsError(ic, invalid_id)) + sort_key = [i in particle_ids_to_move ? 1 : 0 for i in eachparticle(ic)] + # determine a permutation that sorts 'sort_key' in ascending order permutation = sortperm(sort_key) ic.coordinates .= ic.coordinates[:, permutation] diff --git a/test/general/initial_condition.jl b/test/general/initial_condition.jl index 38bab06c58..e98485e1b8 100644 --- a/test/general/initial_condition.jl +++ b/test/general/initial_condition.jl @@ -333,4 +333,67 @@ @test initial_condition.coordinates ≈ expected_coords end end + @testset verbose=true "Move Particles to End" begin + @testset "Valid Particle Movement" begin + # Create a simple initial condition with known particle arrangement + coordinates = [1.0 2.0 3.0 4.0 5.0; + 1.0 2.0 3.0 4.0 5.0] + velocity = [0.1 0.2 0.3 0.4 0.5; + 0.1 0.2 0.3 0.4 0.5] + mass = [1.1, 2.2, 3.3, 4.4, 5.5] + density = [10.0, 20.0, 30.0, 40.0, 50.0] + pressure = [100.0, 200.0, 300.0, 400.0, 500.0] + + ic = InitialCondition(coordinates=coordinates, velocity=velocity, + mass=mass, density=density, pressure=pressure) + + # Move particles 2 and 4 to the end + particle_ids_to_move = [2, 4] + TrixiParticles.move_particles_to_end!(ic, particle_ids_to_move) + + # Expected order after moving: [1, 3, 5, 2, 4] + expected_coordinates = [1.0 3.0 5.0 2.0 4.0; + 1.0 3.0 5.0 2.0 4.0] + expected_velocity = [0.1 0.3 0.5 0.2 0.4; + 0.1 0.3 0.5 0.2 0.4] + expected_mass = [1.1, 3.3, 5.5, 2.2, 4.4] + expected_density = [10.0, 30.0, 50.0, 20.0, 40.0] + expected_pressure = [100.0, 300.0, 500.0, 200.0, 400.0] + + @test ic.coordinates == expected_coordinates + @test ic.velocity == expected_velocity + @test ic.mass == expected_mass + @test ic.density == expected_density + @test ic.pressure == expected_pressure + end + + @testset "Duplicate Particle IDs" begin + coordinates = [1.0 2.0 3.0 4.0; 1.0 2.0 3.0 4.0] + velocity = [0.1 0.2 0.3 0.4; 0.1 0.2 0.3 0.4] + mass = [1.1, 2.2, 3.3, 4.4] + density = [10.0, 20.0, 30.0, 40.0] + + ic = InitialCondition(coordinates=coordinates, velocity=velocity, + mass=mass, density=density) + + # Move particle 2 multiple times (should only move once) + TrixiParticles.move_particles_to_end!(ic, [2, 2, 3]) + + # Expected order: [1, 4, 2, 3] (2 and 3 moved to end) + expected_coordinates = [1.0 4.0 2.0 3.0; 1.0 4.0 2.0 3.0] + expected_mass = [1.1, 4.4, 2.2, 3.3] + + @test ic.coordinates == expected_coordinates + @test ic.mass == expected_mass + end + + @testset "Out of Bounds Particle IDs" begin + ic = rectangular_patch(0.1, (2, 4)) + + # Test with invalid particle ID + @test_throws BoundsError TrixiParticles.move_particles_to_end!(ic, [9]) + @test_throws BoundsError TrixiParticles.move_particles_to_end!(ic, [0]) + @test_throws BoundsError TrixiParticles.move_particles_to_end!(ic, [-1]) + end + end end From 37310d28d6d41d9627e0099a916ecf3b3bfd3f21 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 14 Oct 2025 10:07:00 +0200 Subject: [PATCH 06/20] add vector variant --- src/general/initial_condition.jl | 21 ++- .../structure/total_lagrangian_sph/system.jl | 2 + test/general/initial_condition.jl | 139 +++++++++++------- 3 files changed, 107 insertions(+), 55 deletions(-) diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index 7d2a03e2a1..f09883993a 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -426,9 +426,24 @@ function move_particles_to_end!(ic::InitialCondition, particle_ids_to_move::Vect ic.coordinates .= ic.coordinates[:, permutation] ic.velocity .= ic.velocity[:, permutation] - ic.mass .= ic.mass[permutation] - ic.density .= ic.density[permutation] - ic.pressure .= ic.pressure[permutation] + move_particles_to_end!(ic.mass, particle_ids_to_move) + move_particles_to_end!(ic.density, particle_ids_to_move) + move_particles_to_end!(ic.pressure, particle_ids_to_move) return ic end + +function move_particles_to_end!(a::AbstractVector, particle_ids_to_move::Vector) + invalid_id = findfirst(i -> i <= 0 || i > length(a), particle_ids_to_move) + isnothing(invalid_id) || throw(BoundsError(a, invalid_id)) + + sort_key = [i in particle_ids_to_move ? 1 : 0 for i in eachindex(a)] + # determine a permutation that sorts 'sort_key' in ascending order + permutation = sortperm(sort_key) + + a .= a[permutation] + + return a +end + +move_particles_to_end!(a::Real, particle_ids_to_move::Vector) = a diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index 4ed99b40d8..273f03df8a 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -136,6 +136,8 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing end n_clamped_particles = length(clamped_particles) move_particles_to_end!(initial_condition, clamped_particles) + move_particles_to_end!(young_modulus, clamped_particles) + move_particles_to_end!(poisson_ratio, clamped_particles) end initial_coordinates = copy(initial_condition.coordinates) diff --git a/test/general/initial_condition.jl b/test/general/initial_condition.jl index e98485e1b8..5e168fd4dd 100644 --- a/test/general/initial_condition.jl +++ b/test/general/initial_condition.jl @@ -334,66 +334,101 @@ end end @testset verbose=true "Move Particles to End" begin - @testset "Valid Particle Movement" begin - # Create a simple initial condition with known particle arrangement - coordinates = [1.0 2.0 3.0 4.0 5.0; - 1.0 2.0 3.0 4.0 5.0] - velocity = [0.1 0.2 0.3 0.4 0.5; - 0.1 0.2 0.3 0.4 0.5] - mass = [1.1, 2.2, 3.3, 4.4, 5.5] - density = [10.0, 20.0, 30.0, 40.0, 50.0] - pressure = [100.0, 200.0, 300.0, 400.0, 500.0] - - ic = InitialCondition(coordinates=coordinates, velocity=velocity, - mass=mass, density=density, pressure=pressure) - - # Move particles 2 and 4 to the end - particle_ids_to_move = [2, 4] - TrixiParticles.move_particles_to_end!(ic, particle_ids_to_move) - - # Expected order after moving: [1, 3, 5, 2, 4] - expected_coordinates = [1.0 3.0 5.0 2.0 4.0; - 1.0 3.0 5.0 2.0 4.0] - expected_velocity = [0.1 0.3 0.5 0.2 0.4; - 0.1 0.3 0.5 0.2 0.4] - expected_mass = [1.1, 3.3, 5.5, 2.2, 4.4] - expected_density = [10.0, 30.0, 50.0, 20.0, 40.0] - expected_pressure = [100.0, 300.0, 500.0, 200.0, 400.0] - - @test ic.coordinates == expected_coordinates - @test ic.velocity == expected_velocity - @test ic.mass == expected_mass - @test ic.density == expected_density - @test ic.pressure == expected_pressure - end + @testset verbose=true "InitialCondition Variant" begin + @testset verbose=true "Valid Particle Movement" begin + # Create a simple initial condition with known particle arrangement + coordinates = [1.0 2.0 3.0 4.0 5.0; + 1.0 2.0 3.0 4.0 5.0] + velocity = [0.1 0.2 0.3 0.4 0.5; + 0.1 0.2 0.3 0.4 0.5] + mass = [1.1, 2.2, 3.3, 4.4, 5.5] + density = [10.0, 20.0, 30.0, 40.0, 50.0] + pressure = [100.0, 200.0, 300.0, 400.0, 500.0] + + ic = InitialCondition(coordinates=coordinates, velocity=velocity, + mass=mass, density=density, pressure=pressure) + + # Move particles 2 and 4 to the end + particle_ids_to_move = [2, 4] + TrixiParticles.move_particles_to_end!(ic, particle_ids_to_move) + + # Expected order after moving: [1, 3, 5, 2, 4] + expected_coordinates = [1.0 3.0 5.0 2.0 4.0; + 1.0 3.0 5.0 2.0 4.0] + expected_velocity = [0.1 0.3 0.5 0.2 0.4; + 0.1 0.3 0.5 0.2 0.4] + expected_mass = [1.1, 3.3, 5.5, 2.2, 4.4] + expected_density = [10.0, 30.0, 50.0, 20.0, 40.0] + expected_pressure = [100.0, 300.0, 500.0, 200.0, 400.0] + + @test ic.coordinates == expected_coordinates + @test ic.velocity == expected_velocity + @test ic.mass == expected_mass + @test ic.density == expected_density + @test ic.pressure == expected_pressure + end + + @testset verbose=true "Duplicate Particle IDs" begin + coordinates = [1.0 2.0 3.0 4.0; 1.0 2.0 3.0 4.0] + velocity = [0.1 0.2 0.3 0.4; 0.1 0.2 0.3 0.4] + mass = [1.1, 2.2, 3.3, 4.4] + density = [10.0, 20.0, 30.0, 40.0] + + ic = InitialCondition(coordinates=coordinates, velocity=velocity, + mass=mass, density=density) - @testset "Duplicate Particle IDs" begin - coordinates = [1.0 2.0 3.0 4.0; 1.0 2.0 3.0 4.0] - velocity = [0.1 0.2 0.3 0.4; 0.1 0.2 0.3 0.4] - mass = [1.1, 2.2, 3.3, 4.4] - density = [10.0, 20.0, 30.0, 40.0] + # Move particle 2 multiple times (should only move once) + TrixiParticles.move_particles_to_end!(ic, [2, 2, 3]) - ic = InitialCondition(coordinates=coordinates, velocity=velocity, - mass=mass, density=density) + # Expected order: [1, 4, 2, 3] (2 and 3 moved to end) + expected_coordinates = [1.0 4.0 2.0 3.0; 1.0 4.0 2.0 3.0] + expected_mass = [1.1, 4.4, 2.2, 3.3] - # Move particle 2 multiple times (should only move once) - TrixiParticles.move_particles_to_end!(ic, [2, 2, 3]) + @test ic.coordinates == expected_coordinates + @test ic.mass == expected_mass + end - # Expected order: [1, 4, 2, 3] (2 and 3 moved to end) - expected_coordinates = [1.0 4.0 2.0 3.0; 1.0 4.0 2.0 3.0] - expected_mass = [1.1, 4.4, 2.2, 3.3] + @testset verbose=true "Out of Bounds Particle IDs" begin + ic = rectangular_patch(0.1, (2, 4)) - @test ic.coordinates == expected_coordinates - @test ic.mass == expected_mass + # Test with invalid particle ID + @test_throws BoundsError TrixiParticles.move_particles_to_end!(ic, [9]) + @test_throws BoundsError TrixiParticles.move_particles_to_end!(ic, [0]) + @test_throws BoundsError TrixiParticles.move_particles_to_end!(ic, [-1]) + end end - @testset "Out of Bounds Particle IDs" begin - ic = rectangular_patch(0.1, (2, 4)) + @testset verbose=true "Vector Variant" begin + @testset verbose=true "Valid Particle Movement" begin + a = [1.0, 2.0, 3.0, 4.0, 5.0] + returned = TrixiParticles.move_particles_to_end!(a, [2, 4]) + + expected = [1.0, 3.0, 5.0, 2.0, 4.0] + @test a == expected + @test returned == expected + end - # Test with invalid particle ID - @test_throws BoundsError TrixiParticles.move_particles_to_end!(ic, [9]) - @test_throws BoundsError TrixiParticles.move_particles_to_end!(ic, [0]) - @test_throws BoundsError TrixiParticles.move_particles_to_end!(ic, [-1]) + @testset verbose=true "Duplicate Particle IDs" begin + a = [1, 2, 3, 4] + TrixiParticles.move_particles_to_end!(a, [2, 2, 3]) + + expected = [1, 4, 2, 3] + @test a == expected + end + + @testset verbose=true "Out of Bounds Particle IDs" begin + a = collect(1:8) + + @test_throws BoundsError TrixiParticles.move_particles_to_end!(a, [9]) + @test_throws BoundsError TrixiParticles.move_particles_to_end!(a, [0]) + @test_throws BoundsError TrixiParticles.move_particles_to_end!(a, [-1]) + end + + @testset verbose=true "Empty IDs" begin + a = [10, 20, 30] + TrixiParticles.move_particles_to_end!(a, Int[]) + @test a == [10, 20, 30] + end end end end From d2b8fa7dd5d7a9927cb215deae1db741d72c7c5e Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 14 Oct 2025 11:45:08 +0200 Subject: [PATCH 07/20] force unique indices --- src/schemes/structure/total_lagrangian_sph/system.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index 273f03df8a..f2a277de1b 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -134,6 +134,7 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing if n_clamped_particles > 0 && n_clamped_particles != length(clamped_particles) throw(ArgumentError("`n_clamped_particles` must equal length of `clamped_particles` when both are specified")) end + unique!(clamped_particles) n_clamped_particles = length(clamped_particles) move_particles_to_end!(initial_condition, clamped_particles) move_particles_to_end!(young_modulus, clamped_particles) From 88b09e39355eeda0afbef560f79bef15b71ce530 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 22 Oct 2025 19:27:40 +0200 Subject: [PATCH 08/20] add depwarning --- .../structure/total_lagrangian_sph/system.jl | 34 +++++++++---------- 1 file changed, 16 insertions(+), 18 deletions(-) diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index f2a277de1b..e0d2017b46 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -15,16 +15,6 @@ all relevant quantities and operators are measured with respect to the initial configuration (O’Connor & Rogers 2021, Belytschko et al. 2000). See [Total Lagrangian SPH](@ref tlsph) for more details on the method. -There are two approaches to specify clamped particles in the system: - -1. **Manual ordering**: Provide the number of clamped particles via `n_clamped_particles`, - ensuring that these particles are the last entries in the `InitialCondition` - (see the info box below for details). - -2. **Automatic ordering**: Pass the indices of the particles to be clamped via `clamped_particles`. - In this case, the specified particles will be internally reordered so that they become - the last particles in the `InitialCondition`. - # Arguments - `initial_condition`: Initial condition representing the system's particles. - `young_modulus`: Young's modulus. @@ -35,9 +25,11 @@ There are two approaches to specify clamped particles in the system: See [Smoothing Kernels](@ref smoothing_kernel). # Keywords -- `n_clamped_particles`: Number of clamped particles which are fixed and not integrated +- `n_clamped_particles`(deprecated): Number of clamped particles which are fixed and not integrated to clamp the structure. Note that the clamped particles must be the **last** particles in the `InitialCondition`. See the info box below. + This keyword is deprecated and will be removed in a future release. + Prefer passing `clamped_particles` with the explicit particle indices to be clamped. - `clamped_particles`: A vector of indices specifying the clamped particles which are fixed and not integrated to clamp the structure. - `clamped_particles_motion`: Prescribed motion of the clamped particles. @@ -126,14 +118,20 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end - # Handle clamped particles: - # If only n_clamped_particles > 0: assume manual ordering (do nothing). - # If both are 0/empty: no clamped particles (do nothing). - if !isempty(clamped_particles) - # Check consistency if both methods are used - if n_clamped_particles > 0 && n_clamped_particles != length(clamped_particles) - throw(ArgumentError("`n_clamped_particles` must equal length of `clamped_particles` when both are specified")) + # Backwards compatibility: `n_clamped_particles` is deprecated. + # Emit a deprecation warning and (if the user didn't supply explicit indices) + # convert the old `n_clamped_particles` convention to `clamped_particles`. + if n_clamped_particles != 0 + Base.depwarn("keyword `n_clamped_particles` is deprecated and will be removed in a future release; " * + "pass `clamped_particles` (Vector{Int} of indices) instead.", + :n_clamped_particles) + if isempty(clamped_particles) + clamped_particles = collect((n_particles - n_clamped_particles + 1):n_particles) end + end + + # Handle clamped particles + if !isempty(clamped_particles) unique!(clamped_particles) n_clamped_particles = length(clamped_particles) move_particles_to_end!(initial_condition, clamped_particles) From ed3330442824c77f6b672c15e13aeb2d4dea6121 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 22 Oct 2025 19:37:31 +0200 Subject: [PATCH 09/20] drop `n_particles_clamped` --- .../structure/total_lagrangian_sph/system.jl | 37 ------------------- 1 file changed, 37 deletions(-) diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index e0d2017b46..fc593029c3 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -1,7 +1,6 @@ @doc raw""" TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length, young_modulus, poisson_ratio; - n_clamped_particles=0, clamped_particles::Vector{Int}=Int[], clamped_particles_motion=nothing, acceleration=ntuple(_ -> 0.0, NDIMS), @@ -25,11 +24,6 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method. See [Smoothing Kernels](@ref smoothing_kernel). # Keywords -- `n_clamped_particles`(deprecated): Number of clamped particles which are fixed and not integrated - to clamp the structure. Note that the clamped particles must be the **last** - particles in the `InitialCondition`. See the info box below. - This keyword is deprecated and will be removed in a future release. - Prefer passing `clamped_particles` with the explicit particle indices to be clamped. - `clamped_particles`: A vector of indices specifying the clamped particles which are fixed and not integrated to clamp the structure. - `clamped_particles_motion`: Prescribed motion of the clamped particles. @@ -47,24 +41,6 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method. (which are the quantities of a single particle), returning a `Tuple` or `SVector` that is to be added to the acceleration of that particle. See, for example, [`SourceTermDamping`](@ref). - -!!! note - If specifying the clamped particles manually (via `n_clamped_particles`), - the clamped particles must be the **last** particles in the `InitialCondition`. - To do so, e.g. use the `union` function: - ```jldoctest; output = false, setup = :(clamped_particles = RectangularShape(0.1, (1, 4), (0.0, 0.0), density=1.0); beam = RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0)) - structure = union(beam, clamped_particles) - - # output - ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ - │ InitialCondition{Float64} │ - │ ═════════════════════════ │ - │ #dimensions: ……………………………………………… 2 │ - │ #particles: ………………………………………………… 16 │ - │ particle spacing: ………………………………… 0.1 │ - └──────────────────────────────────────────────────────────────────────────────────────────────────┘ - ``` - where `beam` and `clamped_particles` are of type [`InitialCondition`](@ref). """ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, ARRAY3D, YM, PR, LL, LM, K, PF, V, ST, M, IM, @@ -97,7 +73,6 @@ end function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length, young_modulus, poisson_ratio; - n_clamped_particles=0, clamped_particles::Vector{Int}=Int[], clamped_particles_motion=nothing, acceleration=ntuple(_ -> 0.0, @@ -118,18 +93,6 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end - # Backwards compatibility: `n_clamped_particles` is deprecated. - # Emit a deprecation warning and (if the user didn't supply explicit indices) - # convert the old `n_clamped_particles` convention to `clamped_particles`. - if n_clamped_particles != 0 - Base.depwarn("keyword `n_clamped_particles` is deprecated and will be removed in a future release; " * - "pass `clamped_particles` (Vector{Int} of indices) instead.", - :n_clamped_particles) - if isempty(clamped_particles) - clamped_particles = collect((n_particles - n_clamped_particles + 1):n_particles) - end - end - # Handle clamped particles if !isempty(clamped_particles) unique!(clamped_particles) From c3c2d13063e09bd6a8451c8c9b6516409a7b421d Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 28 Oct 2025 09:27:07 +0100 Subject: [PATCH 10/20] fix regex --- test/examples/examples.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 91ec13ef13..6e3e58adb3 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -8,7 +8,10 @@ @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "structure", "oscillating_beam_2d.jl"), - tspan=(0.0, 0.1)) + tspan=(0.0, 0.1)) [ + r"┌ Warning: keyword `n_clamped_particles` is deprecated.*\n", + r"└ @ Core.*\n" + ] @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end From ffaf028d1ca8451b85f2a2914749e7fbaebd6ae4 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 28 Oct 2025 09:28:44 +0100 Subject: [PATCH 11/20] rm regex --- test/examples/examples.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 6e3e58adb3..91ec13ef13 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -8,10 +8,7 @@ @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "structure", "oscillating_beam_2d.jl"), - tspan=(0.0, 0.1)) [ - r"┌ Warning: keyword `n_clamped_particles` is deprecated.*\n", - r"└ @ Core.*\n" - ] + tspan=(0.0, 0.1)) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end From 76040782c02dfb0423e1c9692b7a0776e5f4bf46 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 28 Oct 2025 11:07:27 +0100 Subject: [PATCH 12/20] fix again --- test/examples/examples.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 6e3e58adb3..0551f4a020 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -10,6 +10,7 @@ "oscillating_beam_2d.jl"), tspan=(0.0, 0.1)) [ r"┌ Warning: keyword `n_clamped_particles` is deprecated.*\n", + r"│\s+caller = .*\n", # allow the caller line inside the boxed warning r"└ @ Core.*\n" ] @test sol.retcode == ReturnCode.Success From 407ef63b31a47ddfbf01058772a714eaceb0df91 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 28 Oct 2025 17:38:28 +0100 Subject: [PATCH 13/20] implement suggestions --- examples/structure/oscillating_beam_2d.jl | 2 +- src/general/initial_condition.jl | 2 +- .../structure/total_lagrangian_sph/system.jl | 18 +++++++++++------- test/examples/examples.jl | 6 +----- 4 files changed, 14 insertions(+), 14 deletions(-) diff --git a/examples/structure/oscillating_beam_2d.jl b/examples/structure/oscillating_beam_2d.jl index 10a0e2137f..a26819b798 100644 --- a/examples/structure/oscillating_beam_2d.jl +++ b/examples/structure/oscillating_beam_2d.jl @@ -61,7 +61,7 @@ smoothing_kernel = WendlandC2Kernel{2}() structure_system = TotalLagrangianSPHSystem(structure, smoothing_kernel, smoothing_length, material.E, material.nu, - n_clamped_particles=nparticles(clamped_particles), + clamped_particles=collect((nparticles(beam) + 1):nparticles(structure)), acceleration=(0.0, -gravity), penalty_force=nothing, viscosity=nothing, clamped_particles_motion=nothing) diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index f09883993a..aeab236ca3 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -421,7 +421,7 @@ function move_particles_to_end!(ic::InitialCondition, particle_ids_to_move::Vect isnothing(invalid_id) || throw(BoundsError(ic, invalid_id)) sort_key = [i in particle_ids_to_move ? 1 : 0 for i in eachparticle(ic)] - # determine a permutation that sorts 'sort_key' in ascending order + # Determine a permutation that sorts 'sort_key' in ascending order permutation = sortperm(sort_key) ic.coordinates .= ic.coordinates[:, permutation] diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index 835c6a7799..1b70e445ee 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -25,11 +25,11 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method. See [Smoothing Kernels](@ref smoothing_kernel). # Keywords -- `n_clamped_particles`(deprecated): Number of clamped particles which are fixed and not integrated +- `n_clamped_particles` (deprecated): Number of clamped particles which are fixed and not integrated to clamp the structure. Note that the clamped particles must be the **last** particles in the `InitialCondition`. See the info box below. This keyword is deprecated and will be removed in a future release. - Prefer passing `clamped_particles` with the explicit particle indices to be clamped. + Instead pass `clamped_particles` with the explicit particle indices to be clamped. - `clamped_particles`: A vector of indices specifying the clamped particles which are fixed and not integrated to clamp the structure. - `clamped_particles_motion`: Prescribed motion of the clamped particles. @@ -98,7 +98,7 @@ end function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length, young_modulus, poisson_ratio; n_clamped_particles=0, - clamped_particles::Vector{Int}=Int[], + clamped_particles::Union{Vector{Int}, AbstractUnitRange}=Int[], clamped_particles_motion=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), @@ -127,16 +127,20 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing :n_clamped_particles) if isempty(clamped_particles) clamped_particles = collect((n_particles - n_clamped_particles + 1):n_particles) + else + throw(ArgumentError("Either `n_clamped_particles` or `clamped_particles` can be specified, not both.")) end end # Handle clamped particles if !isempty(clamped_particles) - unique!(clamped_particles) + + @assert allunique(clamped_particles) "`clamped_particles` contains duplicate particle indices." + n_clamped_particles = length(clamped_particles) - move_particles_to_end!(initial_condition, clamped_particles) - move_particles_to_end!(young_modulus, clamped_particles) - move_particles_to_end!(poisson_ratio, clamped_particles) + move_particles_to_end!(initial_condition, collect(clamped_particles)) + move_particles_to_end!(young_modulus, collect(clamped_particles)) + move_particles_to_end!(poisson_ratio, collect(clamped_particles)) end initial_coordinates = copy(initial_condition.coordinates) diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 0551f4a020..91ec13ef13 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -8,11 +8,7 @@ @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "structure", "oscillating_beam_2d.jl"), - tspan=(0.0, 0.1)) [ - r"┌ Warning: keyword `n_clamped_particles` is deprecated.*\n", - r"│\s+caller = .*\n", # allow the caller line inside the boxed warning - r"└ @ Core.*\n" - ] + tspan=(0.0, 0.1)) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end From aa26aa7fa6cd4fa5b5ed88fd42cb448a23ccc69e Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 28 Oct 2025 17:40:08 +0100 Subject: [PATCH 14/20] apply formatter --- src/schemes/structure/total_lagrangian_sph/system.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index 1b70e445ee..4d237255f5 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -134,7 +134,6 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing # Handle clamped particles if !isempty(clamped_particles) - @assert allunique(clamped_particles) "`clamped_particles` contains duplicate particle indices." n_clamped_particles = length(clamped_particles) From 0e33f31311cb1d07d2e3ab0e8543575276587ff7 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 29 Oct 2025 13:52:21 +0100 Subject: [PATCH 15/20] fix examples --- examples/fsi/falling_water_column_2d.jl | 2 +- examples/structure/oscillating_beam_2d.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/fsi/falling_water_column_2d.jl b/examples/fsi/falling_water_column_2d.jl index 6fec694bfe..7f2688a0ae 100644 --- a/examples/fsi/falling_water_column_2d.jl +++ b/examples/fsi/falling_water_column_2d.jl @@ -68,7 +68,7 @@ structure_system = TotalLagrangianSPHSystem(structure, smoothing_kernel, smoothing_length, material.E, material.nu, boundary_model=boundary_model, - n_clamped_particles=nparticles(clamped_particles), + clamped_particles=(nparticles(beam) + 1):nparticles(structure), acceleration=(0.0, -gravity)) # ========================================================================================== diff --git a/examples/structure/oscillating_beam_2d.jl b/examples/structure/oscillating_beam_2d.jl index a26819b798..13cd830c48 100644 --- a/examples/structure/oscillating_beam_2d.jl +++ b/examples/structure/oscillating_beam_2d.jl @@ -61,7 +61,7 @@ smoothing_kernel = WendlandC2Kernel{2}() structure_system = TotalLagrangianSPHSystem(structure, smoothing_kernel, smoothing_length, material.E, material.nu, - clamped_particles=collect((nparticles(beam) + 1):nparticles(structure)), + clamped_particles=(nparticles(beam) + 1):nparticles(structure), acceleration=(0.0, -gravity), penalty_force=nothing, viscosity=nothing, clamped_particles_motion=nothing) From 213321693e1118b266a2e4a64d04568194c8949b Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 29 Oct 2025 16:14:45 +0100 Subject: [PATCH 16/20] fix --- examples/fsi/dam_break_gate_2d.jl | 2 +- examples/fsi/dam_break_plate_2d.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/fsi/dam_break_gate_2d.jl b/examples/fsi/dam_break_gate_2d.jl index 015232adf9..2bef6f0f50 100644 --- a/examples/fsi/dam_break_gate_2d.jl +++ b/examples/fsi/dam_break_gate_2d.jl @@ -147,7 +147,7 @@ structure_system = TotalLagrangianSPHSystem(structure, structure_smoothing_kernel, structure_smoothing_length, E, nu, boundary_model=boundary_model_structure, - n_clamped_particles=n_particles_x, + clamped_particles=(nparticles(plate) + 1):nparticles(structure), acceleration=(0.0, -gravity)) # ========================================================================================== diff --git a/examples/fsi/dam_break_plate_2d.jl b/examples/fsi/dam_break_plate_2d.jl index 15b795eabf..467c902c0c 100644 --- a/examples/fsi/dam_break_plate_2d.jl +++ b/examples/fsi/dam_break_plate_2d.jl @@ -127,7 +127,7 @@ structure_system = TotalLagrangianSPHSystem(structure, structure_smoothing_kernel, structure_smoothing_length, E, nu, boundary_model=boundary_model_structure, - n_clamped_particles=n_particles_x, + clamped_particles=(nparticles(plate) + 1):nparticles(structure), acceleration=(0.0, -gravity), penalty_force=PenaltyForceGanzenmueller(alpha=0.01)) From 7f849781be3d8da6278684b45bc338045167212c Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 30 Oct 2025 12:41:50 +0100 Subject: [PATCH 17/20] implement suggestions --- src/general/initial_condition.jl | 6 +++--- src/schemes/structure/total_lagrangian_sph/system.jl | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index aeab236ca3..c182a9c334 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -426,9 +426,9 @@ function move_particles_to_end!(ic::InitialCondition, particle_ids_to_move::Vect ic.coordinates .= ic.coordinates[:, permutation] ic.velocity .= ic.velocity[:, permutation] - move_particles_to_end!(ic.mass, particle_ids_to_move) - move_particles_to_end!(ic.density, particle_ids_to_move) - move_particles_to_end!(ic.pressure, particle_ids_to_move) + ic.mass .= ic.mass[permutation] + ic.density .= ic.density[permutation] + ic.pressure .= ic.pressure[permutation] return ic end diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index 4d237255f5..dda54c95eb 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -134,7 +134,7 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing # Handle clamped particles if !isempty(clamped_particles) - @assert allunique(clamped_particles) "`clamped_particles` contains duplicate particle indices." + @assert allunique(clamped_particles) "`clamped_particles` contains duplicate particle indices" n_clamped_particles = length(clamped_particles) move_particles_to_end!(initial_condition, collect(clamped_particles)) From 34f4c3b009688b9a48ce55a3136b8e7efaffd979 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 30 Oct 2025 14:30:05 +0100 Subject: [PATCH 18/20] rm datatype restriction --- src/schemes/structure/total_lagrangian_sph/system.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index dda54c95eb..a3a15f8fc5 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -98,7 +98,7 @@ end function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length, young_modulus, poisson_ratio; n_clamped_particles=0, - clamped_particles::Union{Vector{Int}, AbstractUnitRange}=Int[], + clamped_particles=Int[], clamped_particles_motion=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), From acd73dfef937962bf56efdd5cfd33c76255ae1a4 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 31 Oct 2025 09:54:59 +0100 Subject: [PATCH 19/20] implement suggestions --- .../structure/total_lagrangian_sph/system.jl | 39 +++++++++++-------- 1 file changed, 22 insertions(+), 17 deletions(-) diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index a3a15f8fc5..20cbce72c4 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -2,7 +2,7 @@ TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length, young_modulus, poisson_ratio; n_clamped_particles=0, - clamped_particles::Vector{Int}=Int[], + clamped_particles=Int[], clamped_particles_motion=nothing, acceleration=ntuple(_ -> 0.0, NDIMS), penalty_force=nothing, viscosity=nothing, @@ -25,12 +25,12 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method. See [Smoothing Kernels](@ref smoothing_kernel). # Keywords -- `n_clamped_particles` (deprecated): Number of clamped particles which are fixed and not integrated +- `n_clamped_particles` (deprecated): Number of clamped particles that are fixed and not integrated to clamp the structure. Note that the clamped particles must be the **last** particles in the `InitialCondition`. See the info box below. This keyword is deprecated and will be removed in a future release. Instead pass `clamped_particles` with the explicit particle indices to be clamped. -- `clamped_particles`: A vector of indices specifying the clamped particles which are fixed +- `clamped_particles`: Indices specifying the clamped particles that are fixed and not integrated to clamp the structure. - `clamped_particles_motion`: Prescribed motion of the clamped particles. If `nothing` (default), the clamped particles are fixed. @@ -137,35 +137,40 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing @assert allunique(clamped_particles) "`clamped_particles` contains duplicate particle indices" n_clamped_particles = length(clamped_particles) - move_particles_to_end!(initial_condition, collect(clamped_particles)) - move_particles_to_end!(young_modulus, collect(clamped_particles)) - move_particles_to_end!(poisson_ratio, collect(clamped_particles)) + ic_sorted_indices = deepcopy(initial_condition) + young_modulus_sorted_indices = copy(young_modulus) + poisson_ratio_sorted_indices = copy(poisson_ratio) + move_particles_to_end!(ic_sorted_indices, collect(clamped_particles)) + move_particles_to_end!(young_modulus_sorted_indices, collect(clamped_particles)) + move_particles_to_end!(poisson_ratio_sorted_indices, collect(clamped_particles)) end - initial_coordinates = copy(initial_condition.coordinates) - current_coordinates = copy(initial_condition.coordinates) - mass = copy(initial_condition.mass) - material_density = copy(initial_condition.density) + initial_coordinates = copy(ic_sorted_indices.coordinates) + current_coordinates = copy(ic_sorted_indices.coordinates) + mass = copy(ic_sorted_indices.mass) + material_density = copy(ic_sorted_indices.density) correction_matrix = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) pk1_corrected = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) deformation_grad = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) n_integrated_particles = n_particles - n_clamped_particles - lame_lambda = @. young_modulus * poisson_ratio / - ((1 + poisson_ratio) * (1 - 2 * poisson_ratio)) - lame_mu = @. (young_modulus / 2) / (1 + poisson_ratio) + lame_lambda = @. young_modulus_sorted_indices * poisson_ratio_sorted_indices / + ((1 + poisson_ratio_sorted_indices) * + (1 - 2 * poisson_ratio_sorted_indices)) + lame_mu = @. (young_modulus_sorted_indices / 2) / (1 + poisson_ratio_sorted_indices) ismoving = Ref(!isnothing(clamped_particles_motion)) - initialize_prescribed_motion!(clamped_particles_motion, initial_condition, + initialize_prescribed_motion!(clamped_particles_motion, ic_sorted_indices, n_clamped_particles) - cache = create_cache_tlsph(clamped_particles_motion, initial_condition) + cache = create_cache_tlsph(clamped_particles_motion, ic_sorted_indices) - return TotalLagrangianSPHSystem(initial_condition, initial_coordinates, + return TotalLagrangianSPHSystem(ic_sorted_indices, initial_coordinates, current_coordinates, mass, correction_matrix, pk1_corrected, deformation_grad, material_density, - n_integrated_particles, young_modulus, poisson_ratio, + n_integrated_particles, young_modulus_sorted_indices, + poisson_ratio_sorted_indices, lame_lambda, lame_mu, smoothing_kernel, smoothing_length, acceleration_, boundary_model, penalty_force, viscosity, source_terms, From 79fc6aaaa8abb88f2ee90da717b83ea0f8a9c08b Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 31 Oct 2025 09:56:11 +0100 Subject: [PATCH 20/20] adapt example files --- examples/fsi/dam_break_gate_2d.jl | 4 ++-- examples/fsi/dam_break_plate_2d.jl | 4 ++-- examples/fsi/falling_water_column_2d.jl | 2 +- examples/structure/oscillating_beam_2d.jl | 4 ++-- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/examples/fsi/dam_break_gate_2d.jl b/examples/fsi/dam_break_gate_2d.jl index 2bef6f0f50..7b02e5cc9a 100644 --- a/examples/fsi/dam_break_gate_2d.jl +++ b/examples/fsi/dam_break_gate_2d.jl @@ -96,7 +96,7 @@ clamped_particles = RectangularShape(structure_particle_spacing, (n_particles_x, 1), (plate_position, 0.0), density=structure_density, place_on_shell=true) -structure = union(plate, clamped_particles) +structure = union(clamped_particles, plate) # ========================================================================================== # ==== Fluid @@ -147,7 +147,7 @@ structure_system = TotalLagrangianSPHSystem(structure, structure_smoothing_kernel, structure_smoothing_length, E, nu, boundary_model=boundary_model_structure, - clamped_particles=(nparticles(plate) + 1):nparticles(structure), + clamped_particles=1:nparticles(clamped_particles), acceleration=(0.0, -gravity)) # ========================================================================================== diff --git a/examples/fsi/dam_break_plate_2d.jl b/examples/fsi/dam_break_plate_2d.jl index 467c902c0c..44d86c5317 100644 --- a/examples/fsi/dam_break_plate_2d.jl +++ b/examples/fsi/dam_break_plate_2d.jl @@ -70,7 +70,7 @@ clamped_particles = RectangularShape(structure_particle_spacing, (n_particles_x, 1), plate_position, density=structure_density, place_on_shell=true) -structure = union(plate, clamped_particles) +structure = union(clamped_particles, plate) # ========================================================================================== # ==== Fluid @@ -127,7 +127,7 @@ structure_system = TotalLagrangianSPHSystem(structure, structure_smoothing_kernel, structure_smoothing_length, E, nu, boundary_model=boundary_model_structure, - clamped_particles=(nparticles(plate) + 1):nparticles(structure), + clamped_particles=1:nparticles(clamped_particles), acceleration=(0.0, -gravity), penalty_force=PenaltyForceGanzenmueller(alpha=0.01)) diff --git a/examples/fsi/falling_water_column_2d.jl b/examples/fsi/falling_water_column_2d.jl index 7f2688a0ae..45898c7240 100644 --- a/examples/fsi/falling_water_column_2d.jl +++ b/examples/fsi/falling_water_column_2d.jl @@ -68,7 +68,7 @@ structure_system = TotalLagrangianSPHSystem(structure, smoothing_kernel, smoothing_length, material.E, material.nu, boundary_model=boundary_model, - clamped_particles=(nparticles(beam) + 1):nparticles(structure), + clamped_particles=1:nparticles(clamped_particles), acceleration=(0.0, -gravity)) # ========================================================================================== diff --git a/examples/structure/oscillating_beam_2d.jl b/examples/structure/oscillating_beam_2d.jl index 13cd830c48..697e8e475d 100644 --- a/examples/structure/oscillating_beam_2d.jl +++ b/examples/structure/oscillating_beam_2d.jl @@ -52,7 +52,7 @@ n_particles_per_dimension = (round(Int, elastic_beam.length / particle_spacing) beam = RectangularShape(particle_spacing, n_particles_per_dimension, (0.0, 0.0), density=material.density, place_on_shell=true) -structure = union(beam, clamped_particles) +structure = union(clamped_particles, beam) # ========================================================================================== # ==== Structure @@ -61,7 +61,7 @@ smoothing_kernel = WendlandC2Kernel{2}() structure_system = TotalLagrangianSPHSystem(structure, smoothing_kernel, smoothing_length, material.E, material.nu, - clamped_particles=(nparticles(beam) + 1):nparticles(structure), + clamped_particles=1:nparticles(clamped_particles), acceleration=(0.0, -gravity), penalty_force=nothing, viscosity=nothing, clamped_particles_motion=nothing)