Merge branch 'main' into oriented-bbox

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/general/interpolation.jl

@@ -300,12 +300,7 @@ function interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_syst
         throw(ArgumentError("`interpolate_plane_3d` requires a 3D simulation"))
     end
 
-    points_coords, resolution_ = sample_plane((point1, point2, point3), resolution)
-
-    if !isapprox(resolution, resolution_, rtol=5e-2)
-        @info "The desired plane size is not a multiple of the resolution $resolution." *
-              "\nNew resolution is set to $resolution_."
-    end
+    points_coords = sample_face_with_fixed_size((point1, point2, point3), resolution)
 
     # Interpolate using the generated points
     results = interpolate_points(points_coords, semi, ref_system, v_ode, u_ode;
@@ -766,3 +761,49 @@ function divide_by_shepard_coefficient!(field::AbstractArray{T, 3}, shepard_coef
 
     return field
 end
+
+function sample_face_with_fixed_size(face_points::NTuple{3}, resolution)
+    # Verify that points are in 3D space
+    if any(length.(face_points) .!= 3)
+        throw(ArgumentError("all points must be 3D coordinates"))
+    end
+
+    point1_ = SVector{3}(face_points[1])
+    point2_ = SVector{3}(face_points[2])
+    point3_ = SVector{3}(face_points[3])
+
+    # Vectors defining the edges of the parallelogram
+    edge1 = point2_ - point1_
+    edge2 = point3_ - point1_
+
+    # Check if the points are collinear
+    if isapprox(norm(cross(edge1, edge2)), 0; atol=eps())
+        throw(ArgumentError("the vectors `AB` and `AC` must not be collinear"))
+    end
+
+    # Determine the number of points along each edge
+    num_points_edge1 = ceil(Int, norm(edge1) / resolution)
+    num_points_edge2 = ceil(Int, norm(edge2) / resolution)
+
+    coords = zeros(3, (num_points_edge1 + 1) * (num_points_edge2 + 1))
+
+    index = 1
+    for i in 0:num_points_edge1
+        for j in 0:num_points_edge2
+            point_on_plane = point1_ + (i / num_points_edge1) * edge1 +
+                             (j / num_points_edge2) * edge2
+            coords[:, index] = point_on_plane
+            index += 1
+        end
+    end
+
+    resolution_new = min(norm(edge1 / num_points_edge1),
+                         norm(edge2 / num_points_edge2))
+
+    if !isapprox(resolution, resolution_new, rtol=sqrt(eps(resolution)))
+        @info "The desired face size is not a multiple of the resolution $resolution." *
+              "\nNew resolution is set to $resolution_new."
+    end
+
+    return coords
+end

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,