Automatic ordering of clamped particles (#931)

* add automatic reordering

* fix typo

* revise docs

* indentation

* add tests

* add vector variant

* force unique indices

* add depwarning

* fix regex

* fix again

* implement suggestions

* apply formatter

* fix examples

* fix

* implement suggestions

* rm datatype restriction

* implement suggestions

* adapt example files

* fix

---------

Co-authored-by: LasNikas <[email protected]>

Author: LasNikas

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,
-                                            n_clamped_particles=n_particles_x,
+                                            clamped_particles=1:nparticles(clamped_particles),
                                             acceleration=(0.0, -gravity))
 
 # ==========================================================================================

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,
-                                            n_clamped_particles=n_particles_x,
+                                            clamped_particles=1:nparticles(clamped_particles),
                                             acceleration=(0.0, -gravity),
                                             penalty_force=PenaltyForceGanzenmueller(alpha=0.01))
 

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=1:nparticles(clamped_particles),
                                             acceleration=(0.0, -gravity))
 
 # ==========================================================================================

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,
-                                            n_clamped_particles=nparticles(clamped_particles),
+                                            clamped_particles=1:nparticles(clamped_particles),
                                             acceleration=(0.0, -gravity),
                                             penalty_force=nothing, viscosity=nothing,
                                             clamped_particles_motion=nothing)

src/general/initial_condition.jl

@@ -415,3 +415,35 @@ function find_too_close_particles(coords, min_distance)
 
     return result
 end
+
+function move_particles_to_end!(ic::InitialCondition, particle_ids_to_move)
+    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]
+    ic.velocity .= ic.velocity[:, permutation]
+    ic.mass .= ic.mass[permutation]
+    ic.density .= ic.density[permutation]
+    ic.pressure .= ic.pressure[permutation]
+
+    return ic
+end
+
+function move_particles_to_end!(a::AbstractVector, particle_ids_to_move)
+    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) = a

src/schemes/structure/total_lagrangian_sph/system.jl

@@ -1,10 +1,12 @@
 @doc raw"""
     TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length,
                              young_modulus, poisson_ratio;
-                             n_clamped_particles=0, clamped_particles_motion=nothing,
-                             boundary_model=nothing,
+                             n_clamped_particles=0,
+                             clamped_particles=Int[],
+                             clamped_particles_motion=nothing,
                              acceleration=ntuple(_ -> 0.0, NDIMS),
-                             penalty_force=nothing, source_terms=nothing)
+                             penalty_force=nothing, viscosity=nothing,
+                             source_terms=nothing, boundary_model=nothing)
 
 System for particles of an elastic structure.
 
@@ -23,9 +25,13 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method.
                         See [Smoothing Kernels](@ref smoothing_kernel).
 
 # Keywords
-- `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.
+- `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`: 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.
                     See [`PrescribedMotion`](@ref) for details.
@@ -43,7 +49,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 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)
@@ -90,12 +97,13 @@ end
 
 function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length,
                                   young_modulus, poisson_ratio;
-                                  n_clamped_particles=0, clamped_particles_motion=nothing,
-                                  boundary_model=nothing,
+                                  n_clamped_particles=0,
+                                  clamped_particles=Int[],
+                                  clamped_particles_motion=nothing,
                                   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)
@@ -110,30 +118,63 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing
         throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem"))
     end
 
-    initial_coordinates = copy(initial_condition.coordinates)
-    current_coordinates = copy(initial_condition.coordinates)
-    mass = copy(initial_condition.mass)
-    material_density = copy(initial_condition.density)
+    # 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)
+        else
+            throw(ArgumentError("Either `n_clamped_particles` or `clamped_particles` can be specified, not both."))
+        end
+    end
+
+    # Handle clamped particles
+    if !isempty(clamped_particles)
+        @assert allunique(clamped_particles) "`clamped_particles` contains duplicate particle indices"
+
+        n_clamped_particles = length(clamped_particles)
+        initial_condition_sorted = deepcopy(initial_condition)
+        young_modulus_sorted = copy(young_modulus)
+        poisson_ratio_sorted = copy(poisson_ratio)
+        move_particles_to_end!(initial_condition_sorted, clamped_particles)
+        move_particles_to_end!(young_modulus_sorted, clamped_particles)
+        move_particles_to_end!(poisson_ratio_sorted, clamped_particles)
+    else
+        initial_condition_sorted = initial_condition
+        young_modulus_sorted = young_modulus
+        poisson_ratio_sorted = poisson_ratio
+    end
+
+    initial_coordinates = copy(initial_condition_sorted.coordinates)
+    current_coordinates = copy(initial_condition_sorted.coordinates)
+    mass = copy(initial_condition_sorted.mass)
+    material_density = copy(initial_condition_sorted.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 * poisson_ratio_sorted /
+                     ((1 + poisson_ratio_sorted) *
+                      (1 - 2 * poisson_ratio_sorted))
+    lame_mu = @. (young_modulus_sorted / 2) / (1 + poisson_ratio_sorted)
 
     ismoving = Ref(!isnothing(clamped_particles_motion))
-    initialize_prescribed_motion!(clamped_particles_motion, initial_condition,
+    initialize_prescribed_motion!(clamped_particles_motion, initial_condition_sorted,
                                   n_clamped_particles)
 
-    cache = create_cache_tlsph(clamped_particles_motion, initial_condition)
+    cache = create_cache_tlsph(clamped_particles_motion, initial_condition_sorted)
 
-    return TotalLagrangianSPHSystem(initial_condition, initial_coordinates,
+    return TotalLagrangianSPHSystem(initial_condition_sorted, 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,
+                                    poisson_ratio_sorted,
                                     lame_lambda, lame_mu, smoothing_kernel,
                                     smoothing_length, acceleration_, boundary_model,
                                     penalty_force, viscosity, source_terms,

test/general/initial_condition.jl

@@ -333,4 +333,102 @@
             @test initial_condition.coordinates ≈ expected_coords
         end
     end
+    @testset verbose=true "Move Particles to End" begin
+        @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)
+
+                # 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 verbose=true "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
+
+        @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
+
+            @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