diff --git a/examples/fsi/dam_break_gate_2d.jl b/examples/fsi/dam_break_gate_2d.jl index 015232adf9..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, - n_clamped_particles=n_particles_x, + 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 15b795eabf..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, - n_clamped_particles=n_particles_x, + 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 6fec694bfe..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, - n_clamped_particles=nparticles(clamped_particles), + 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 10a0e2137f..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, - 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) diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index 2e7ac773b8..c182a9c334 100644 --- a/src/general/initial_condition.jl +++ b/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::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] + 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::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 4bcf567874..3c3ca51754 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, clamped_particles_motion=nothing, - boundary_model=nothing, + 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 +24,8 @@ 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. +- `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. @@ -41,23 +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 - 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, @@ -90,12 +73,12 @@ end function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length, young_modulus, poisson_ratio; - n_clamped_particles=0, clamped_particles_motion=nothing, - boundary_model=nothing, + 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 +93,45 @@ 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) + # Handle clamped particles + if !isempty(clamped_particles) + @assert allunique(clamped_particles) "`clamped_particles` contains duplicate particle indices" + + n_clamped_particles = length(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(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, diff --git a/test/general/initial_condition.jl b/test/general/initial_condition.jl index 38bab06c58..5e168fd4dd 100644 --- a/test/general/initial_condition.jl +++ b/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