Implement PrescribedMotion for clamped TLSPH particles (#896)

* Prepare for moving clamped particles

* Implement prescribed motion for clamped particles

* Fix

* Add test

* Add docs

* Reformat

* Fix tests

* Fix `current_velocity`

* Fix

* Fix `system_data`

* Fix GPU tests

* Implement suggestions

---------

Co-authored-by: Niklas Neher <[email protected]>
Co-authored-by: Sven Berger <[email protected]>

Author: 3 people

examples/structure/oscillating_beam_2d.jl

@@ -63,7 +63,8 @@ structure_system = TotalLagrangianSPHSystem(structure, smoothing_kernel, smoothi
                                             material.E, material.nu,
                                             n_clamped_particles=nparticles(clamped_particles),
                                             acceleration=(0.0, -gravity),
-                                            penalty_force=nothing, viscosity=nothing)
+                                            penalty_force=nothing, viscosity=nothing,
+                                            clamped_particles_motion=nothing)
 
 # ==========================================================================================
 # ==== Simulation

src/general/semidiscretization.jl

@@ -402,10 +402,14 @@ end
 function initialize_neighborhood_searches!(semi)
     foreach_system(semi) do system
         foreach_system(semi) do neighbor
+            # TODO Initialize after adapting to the GPU.
+            # Currently, this cannot use `semi.parallelization_backend`
+            # because data is still on the CPU.
             PointNeighbors.initialize!(get_neighborhood_search(system, neighbor, semi),
                                        initial_coordinates(system),
                                        initial_coordinates(neighbor),
-                                       eachindex_y=each_active_particle(neighbor))
+                                       eachindex_y=each_active_particle(neighbor),
+                                       parallelization_backend=PolyesterBackend())
         end
     end
 

src/io/write_vtk.jl

@@ -381,6 +381,9 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem)
                            initial_coords(system, particle)
                            for particle in eachparticle(system)]
 
+    vtk["is_clamped"] = vcat(fill(0, system.n_integrated_particles),
+                             fill(1, nparticles(system) - system.n_integrated_particles))
+
     vtk["lame_lambda"] = system.lame_lambda
     vtk["lame_mu"] = system.lame_mu
     vtk["young_modulus"] = system.young_modulus

src/schemes/boundary/prescribed_motion.jl

@@ -10,8 +10,10 @@
                        are moving or not. Its boolean return value determines
                        if the neighborhood search will be updated.
 
-# Keywords
-- `moving_particles`: Indices of moving particles. Default is each particle in the system.
+# Keyword Arguments
+- `moving_particles`: Indices of moving particles. Default is each particle in the system
+                      for the [`WallBoundarySystem`](@ref) and all clamped particles
+                      for the [`TotalLagrangianSPHSystem`](@ref).
 
 # Examples
 ```jldoctest; output = false
@@ -48,7 +50,9 @@ function PrescribedMotion(movement_function, is_moving; moving_particles=nothing
     return PrescribedMotion(movement_function, is_moving, moving_particles)
 end
 
-function initialize!(prescribed_motion::PrescribedMotion, initial_condition)
+function initialize_prescribed_motion!(prescribed_motion::PrescribedMotion,
+                                       initial_condition,
+                                       n_clamped_particles=nparticles(initial_condition))
     # Test `movement_function` return type
     pos = extract_svector(initial_condition.coordinates,
                           Val(size(initial_condition.coordinates, 1)), 1)
@@ -57,23 +61,34 @@ function initialize!(prescribed_motion::PrescribedMotion, initial_condition)
               "Returning regular `Vector`s causes allocations and significant performance overhead."
     end
 
-    # Empty `moving_particles` means all particles are moving
+    # Empty `moving_particles` means all clamped particles are moving.
+    # For boundaries, all particles are considered clamped.
     if isempty(prescribed_motion.moving_particles)
         # Default is an empty vector, since the number of particles is not known when
         # instantiating `PrescribedMotion`.
-        resize!(prescribed_motion.moving_particles, nparticles(initial_condition))
-        prescribed_motion.moving_particles .= collect(1:nparticles(initial_condition))
+        resize!(prescribed_motion.moving_particles, n_clamped_particles)
+
+        # Clamped particles for TLSPH are the last `n_clamped_particles` in the system
+        first_particle = nparticles(initial_condition) - n_clamped_particles + 1
+        clamped_particles = first_particle:nparticles(initial_condition)
+        prescribed_motion.moving_particles .= collect(clamped_particles)
     end
+
+    return prescribed_motion
+end
+
+function initialize_prescribed_motion!(::Nothing, initial_condition,
+                                       n_clamped_particles=nparticles(initial_condition))
+    return nothing
 end
 
-function (prescribed_motion::PrescribedMotion)(system, t, semi)
-    (; coordinates, cache) = system
+function (prescribed_motion::PrescribedMotion)(coordinates, velocity, acceleration,
+                                               ismoving, system, semi, t)
     (; movement_function, is_moving, moving_particles) = prescribed_motion
-    (; acceleration, velocity) = cache
 
-    system.ismoving[] = is_moving(t)
+    ismoving[] = is_moving(t)
 
-    is_moving(t) || return system
+    is_moving(t) || return nothing
 
     @threaded semi for particle in moving_particles
         pos_original = initial_coords(system, particle)
@@ -90,13 +105,7 @@ function (prescribed_motion::PrescribedMotion)(system, t, semi)
         end
     end
 
-    return system
-end
-
-function (prescribed_motion::Nothing)(system::AbstractSystem, t, semi)
-    system.ismoving[] = false
-
-    return system
+    return nothing
 end
 
 """

src/schemes/boundary/wall_boundary/system.jl

@@ -42,12 +42,11 @@ function WallBoundarySystem(initial_condition, model; prescribed_motion=nothing,
     coordinates = copy(initial_condition.coordinates)
 
     ismoving = Ref(!isnothing(prescribed_motion))
-    initialize!(prescribed_motion, initial_condition)
+    initialize_prescribed_motion!(prescribed_motion, initial_condition)
 
     cache = create_cache_boundary(prescribed_motion, initial_condition)
     cache = (cache..., color=Int(color_value))
 
-    # Because of dispatches boundary model needs to be first!
     return WallBoundarySystem(initial_condition, coordinates, model, prescribed_motion,
                               ismoving, adhesion_coefficient, cache)
 end
@@ -105,7 +104,7 @@ end
     return current_velocity(v, system, system.prescribed_motion, particle)
 end
 
-@inline function current_velocity(v, system, movement, particle)
+@inline function current_velocity(v, system, prescribed_motion, particle)
     (; cache, ismoving) = system
 
     if ismoving[]
@@ -115,7 +114,7 @@ end
     return zero(SVector{ndims(system), eltype(system)})
 end
 
-@inline function current_velocity(v, system, movement::Nothing, particle)
+@inline function current_velocity(v, system, prescribed_motion::Nothing, particle)
     return zero(SVector{ndims(system), eltype(system)})
 end
 
@@ -127,7 +126,8 @@ end
     return current_acceleration(system, system.prescribed_motion, particle)
 end
 
-@inline function current_acceleration(system, ::PrescribedMotion, particle)
+@inline function current_acceleration(system::WallBoundarySystem, ::PrescribedMotion,
+                                      particle)
     (; cache, ismoving) = system
 
     if ismoving[]
@@ -137,7 +137,7 @@ end
     return zero(SVector{ndims(system), eltype(system)})
 end
 
-@inline function current_acceleration(system, ::Nothing, particle)
+@inline function current_acceleration(system::WallBoundarySystem, ::Nothing, particle)
     return zero(SVector{ndims(system), eltype(system)})
 end
 
@@ -177,7 +177,21 @@ end
 function update_positions!(system::WallBoundarySystem, v, u, v_ode, u_ode, semi, t)
     (; prescribed_motion) = system
 
-    prescribed_motion(system, t, semi)
+    apply_prescribed_motion!(system, prescribed_motion, semi, t)
+end
+
+function apply_prescribed_motion!(system::WallBoundarySystem,
+                                  prescribed_motion::PrescribedMotion, semi, t)
+    (; ismoving, coordinates, cache) = system
+    (; acceleration, velocity) = cache
+
+    prescribed_motion(coordinates, velocity, acceleration, ismoving, system, semi, t)
+
+    return system
+end
+
+function apply_prescribed_motion!(system::WallBoundarySystem, ::Nothing, semi, t)
+    return system
 end
 
 function update_quantities!(system::WallBoundarySystem, v, u, v_ode, u_ode, semi, t)