diff --git a/docs/literate/src/tut_rigid_body_fsi.jl b/docs/literate/src/tut_rigid_body_fsi.jl
index 5f11a3d561..12482b70de 100644
--- a/docs/literate/src/tut_rigid_body_fsi.jl
+++ b/docs/literate/src/tut_rigid_body_fsi.jl
@@ -227,10 +227,15 @@ nothing # hide
 
 # ## Step 3: With contact model
 # Finally, we add a `contact_model` to handle collisions between rigid bodies and between
-# rigid bodies and the tank.
+# rigid bodies and the tank. Here we also enable the frictional rigid-wall path, so we need
+# `UpdateCallback()` to update the tangential contact history between time steps.
 
 contact_model = RigidContactModel(; normal_stiffness=2.0e5,
                                   normal_damping=150.0,
+                                  static_friction_coefficient=0.6,
+                                  kinetic_friction_coefficient=0.4,
+                                  tangential_stiffness=1.0e5,
+                                  tangential_damping=150.0,
                                   contact_distance=2.0 * structure_particle_spacing)
 nothing # hide
 
@@ -253,15 +258,15 @@ nothing # hide
 
 info_callback = InfoCallback(interval=100) # hide
 saving_callback_step3 = SolutionSavingCallback(dt=0.02, prefix="step3") # hide
-callbacks_step3 = CallbackSet(info_callback, saving_callback_step3) # hide
+callbacks_step3 = CallbackSet(info_callback, saving_callback_step3, UpdateCallback()) # hide
 nothing # hide
 
 # ```julia
 # sol_step3 = solve(ode_step3, RDPK3SpFSAL49(), save_everystep=false,
-#                   callback=callbacks,  abstol=1e-6, reltol=1e-4, dtmax=2e-3)
+#                   callback=callbacks_step3,  abstol=1e-6, reltol=1e-4, dtmax=2e-3)
 # ```
 sol_step3 = solve(ode_step3, RDPK3SpFSAL49(), abstol=1e-6, reltol=1e-4, dtmax=2e-3, # hide
-                  save_everystep=false) # hide
+                  save_everystep=false, callback=callbacks_step3) # hide
 nothing # hide
 
 # The plot now shows the full simulation. The squares collide with the tank bottom and each other.
@@ -320,15 +325,15 @@ nothing # hide
 
 info_callback = InfoCallback(interval=100) # hide
 saving_callback_step4 = SolutionSavingCallback(dt=0.02, prefix="step4") # hide
-callbacks_step4 = CallbackSet(info_callback, saving_callback_step4) # hide
+callbacks_step4 = CallbackSet(info_callback, saving_callback_step4, UpdateCallback()) # hide
 nothing # hide
 
 # ```julia
 # sol_step4 = solve(ode_step4, RDPK3SpFSAL49(), save_everystep=false,
-#                   callback=callbacks,  abstol=1e-6, reltol=1e-4, dtmax=2e-3)
+#                   callback=callbacks_step4,  abstol=1e-6, reltol=1e-4, dtmax=2e-3)
 # ```
 sol_step4 = solve(ode_step4, RDPK3SpFSAL49(), abstol=1e-6, reltol=1e-4, dtmax=2e-3, # hide
-                  save_everystep=false) # hide
+                  save_everystep=false, callback=callbacks_step4) # hide
 nothing # hide
 
 # And here is the final plot with circles instead of squares.
@@ -372,9 +377,10 @@ semi_next = Semidiscretization(fluid_system, boundary_system,
                                small_sphere_systems...)
 
 ode_step_next = semidiscretize(semi_next, tspan) #hide
+callbacks_step_next = CallbackSet(UpdateCallback()) # hide
 sol_step_next = solve(ode_step_next, RDPK3SpFSAL49(), # hide
                       abstol=1e-6, reltol=1e-4, dtmax=2e-3, # hide
-                      save_everystep=false) # hide
+                      save_everystep=false, callback=callbacks_step_next) # hide
 nothing # hide
 
 plot(sol_step_next, legend=nothing) #hide
@@ -428,8 +434,9 @@ hexagon_system = RigidBodySystem(hexagon_shape;
 semi_hexagon = Semidiscretization(fluid_system, boundary_system, hexagon_system)
 ode_step_hex = semidiscretize(semi_hexagon, (0.0, 0.4))
 
+callbacks_step_hex = CallbackSet(UpdateCallback()) # hide
 sol_step_hex = solve(ode_step_hex, RDPK3SpFSAL49(), abstol=1e-6, reltol=1e-5, dtmax=1e-3,
-                     save_everystep=false)
+                     save_everystep=false, callback=callbacks_step_hex)
 plot(sol_step_hex, legend=nothing) #hide
 plot!(dpi=200) # hide
 savefig("tut_rigid_body_fsi_hex.png"); # hide
diff --git a/docs/src/callbacks.md b/docs/src/callbacks.md
index c96ccc7dba..c900458f48 100644
--- a/docs/src/callbacks.md
+++ b/docs/src/callbacks.md
@@ -1,5 +1,10 @@
 # Callbacks
 
+[`UpdateCallback`](@ref) is required for systems that keep mutable state between time
+steps. In the current rigid-contact implementation, this applies when a
+[`RigidContactModel`](@ref) uses tangential history, i.e. the rigid-wall friction path
+needs the tangential displacement cache to be updated between steps.
+
 ```@autodocs
 Modules = [TrixiParticles]
 Pages = map(file -> joinpath("callbacks", file), readdir(joinpath("..", "src", "callbacks")))
diff --git a/docs/src/systems/rigid_body.md b/docs/src/systems/rigid_body.md
index 08e50e13c1..c3201af6bf 100644
--- a/docs/src/systems/rigid_body.md
+++ b/docs/src/systems/rigid_body.md
@@ -17,6 +17,39 @@ Rigid contact is configured through the contact model. This is separate from the
 boundary model used for fluid-structure interaction; see
 [Boundary Models](@ref boundary_models) for that part of the rigid-body setup.
 
+`RigidContactModel` defines the rigid-contact law shared by rigid-wall and rigid-rigid
+interaction. The always-active parameters are `normal_stiffness`, `normal_damping`, and
+`contact_distance`.
+
+For rigid-wall contact, the current implementation also supports tangential friction with
+the parameters `static_friction_coefficient`, `kinetic_friction_coefficient`,
+`tangential_stiffness`, `tangential_damping`, `stick_velocity_tolerance`, and
+`penetration_slop`. When the tangential spring history is active, this requires
+[`UpdateCallback`](@ref) so the tangential displacement cache is updated between time
+steps.
+
+The current implementation uses the same model for rigid-wall and rigid-rigid contact, but
+the active behavior is not identical yet:
+
+- rigid-wall contact groups penetrating wall neighbors into a small number of contact
+  manifolds per rigid particle and applies one normal-plus-tangential contact force per
+  manifold,
+- rigid-rigid contact evaluates direct pairwise normal contact forces between rigid
+  particles,
+- and rigid-rigid remains normal-only for now, even though the tangential-history key
+  design is already shared with the wall path.
+
+`contact_distance` defines when contact starts. If `contact_distance == 0`, the
+particle spacing of the `RigidBodySystem` is used when the contact model is adapted to
+the runtime system.
+
+If no `contact_model` is specified for a rigid body, rigid-wall and rigid-rigid contact
+for that system are disabled.
+
+For output and postprocessing, rigid bodies also expose the diagnostics
+`contact_count` and `max_contact_penetration`. They are available through rigid-body
+system data and VTK output.
+
 ```@autodocs
 Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "structure", "rigid_body", "contact_models.jl")]
diff --git a/examples/fsi/falling_rigid_spheres_2d.jl b/examples/fsi/falling_rigid_spheres_2d.jl
index 026e35793e..51b8a7a7f7 100644
--- a/examples/fsi/falling_rigid_spheres_2d.jl
+++ b/examples/fsi/falling_rigid_spheres_2d.jl
@@ -99,9 +99,14 @@ boundary_model_structure_2 = BoundaryModelDummyParticles(hydrodynamic_densities_
                                                          fluid_smoothing_kernel,
                                                          fluid_smoothing_length)
 
-# Basic rigid contact model used for both rigid bodies.
+# Use the frictional rigid-wall contact path, which requires `UpdateCallback()` to keep
+# the tangential contact history in sync between time steps.
 contact_model = RigidContactModel(; normal_stiffness=2.0e5,
                                   normal_damping=200.0,
+                                  static_friction_coefficient=0.6,
+                                  kinetic_friction_coefficient=0.4,
+                                  tangential_stiffness=1.0e5,
+                                  tangential_damping=150.0,
                                   contact_distance=2.0 *
                                                    structure_particle_spacing)
 
@@ -125,7 +130,7 @@ ode = semidiscretize(semi, tspan)
 info_callback = InfoCallback(interval=50)
 saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", prefix="")
 
-callbacks = CallbackSet(info_callback, saving_callback)
+callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback())
 
 # Use a Runge-Kutta method with automatic (error based) time step size control.
 sol = solve(ode, RDPK3SpFSAL49(),
diff --git a/examples/fsi/falling_rotating_rigid_squares_2d.jl b/examples/fsi/falling_rotating_rigid_squares_2d.jl
index 5eb1160ea0..61dcc284ef 100644
--- a/examples/fsi/falling_rotating_rigid_squares_2d.jl
+++ b/examples/fsi/falling_rotating_rigid_squares_2d.jl
@@ -109,13 +109,22 @@ end
 boundary_model_structure_1 = structure_boundary_model(square1)
 boundary_model_structure_2 = structure_boundary_model(square2)
 
-# Use a less dissipative wall contact for the denser square so its rebound is more visible.
+# Use frictional rigid-wall contact so the rotating squares exchange tangential impulses with
+# the tank floor as well as normal contact forces. This requires `UpdateCallback()`.
 contact_model_1 = RigidContactModel(; normal_stiffness=2.0e5,
                                     normal_damping=200.0,
+                                    static_friction_coefficient=0.6,
+                                    kinetic_friction_coefficient=0.4,
+                                    tangential_stiffness=1.0e5,
+                                    tangential_damping=180.0,
                                     contact_distance=2.0 *
                                                      structure_particle_spacing)
 contact_model_2 = RigidContactModel(; normal_stiffness=2.0e5,
                                     normal_damping=80.0,
+                                    static_friction_coefficient=0.5,
+                                    kinetic_friction_coefficient=0.3,
+                                    tangential_stiffness=8.0e4,
+                                    tangential_damping=120.0,
                                     contact_distance=2.0 *
                                                      structure_particle_spacing)
 
@@ -144,7 +153,7 @@ saving_callback = SolutionSavingCallback(dt=0.01,
                                          output_directory="out",
                                          prefix="")
 
-callbacks = CallbackSet(info_callback, saving_callback)
+callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback())
 
 # Use a Runge-Kutta method with automatic (error based) time step size control.
 # To prevent penetration of fluid particles through the rigid bodies or the boundary
diff --git a/examples/fsi/falling_rotating_rigid_squares_w_buoys_2d.jl b/examples/fsi/falling_rotating_rigid_squares_w_buoys_2d.jl
index 37efb627a9..2e647fcbf2 100644
--- a/examples/fsi/falling_rotating_rigid_squares_w_buoys_2d.jl
+++ b/examples/fsi/falling_rotating_rigid_squares_w_buoys_2d.jl
@@ -20,6 +20,10 @@ small_sphere_y = initial_fluid_size[2] + small_sphere_radius
 small_sphere_x_positions = 0.2:(3 * small_sphere_radius):1.8
 small_sphere_contact_model = RigidContactModel(; normal_stiffness=2.0e5,
                                                normal_damping=120.0,
+                                               static_friction_coefficient=0.6,
+                                               kinetic_friction_coefficient=0.4,
+                                               tangential_stiffness=1.0e5,
+                                               tangential_damping=150.0,
                                                contact_distance=2.0 *
                                                                 structure_particle_spacing)
 extra_structure_systems = [begin
diff --git a/examples/structure/sliding_rigid_squares_friction_2d.jl b/examples/structure/sliding_rigid_squares_friction_2d.jl
new file mode 100644
index 0000000000..c28f9e6566
--- /dev/null
+++ b/examples/structure/sliding_rigid_squares_friction_2d.jl
@@ -0,0 +1,99 @@
+# ==========================================================================================
+# 2D Sliding Rigid Squares with and without Wall Friction
+#
+# Two identical rigid squares slide on the same floor. The left square uses the normal-only
+# rigid contact model from PR1, while the right square uses the frictional wall-contact path
+# added in PR2. The frictional square slows down and starts rotating due to tangential wall
+# forces, whereas the normal-only square keeps sliding without spin-up.
+#
+# In ParaView, compare the trajectories and the rigid-body field data such as
+# `angular_velocity`, `contact_count`, and `max_contact_penetration`.
+# ==========================================================================================
+
+using TrixiParticles
+using OrdinaryDiffEq
+
+# ==========================================================================================
+# ==== Resolution
+particle_spacing = 0.03
+boundary_layers = 3
+
+# ==========================================================================================
+# ==== Experiment Setup
+gravity = 9.81
+tspan = (0.0, 0.8)
+
+square_side_length = 0.18
+square_density = 1000.0
+square_particles_per_side = round(Int, square_side_length / particle_spacing)
+square_bottom_y = 0.03
+
+square_frictionless = RectangularShape(particle_spacing,
+                                       (square_particles_per_side,
+                                        square_particles_per_side),
+                                       (-1.0, square_bottom_y),
+                                       density=square_density,
+                                       velocity=(1.0, 0.0))
+square_frictional = RectangularShape(particle_spacing,
+                                     (square_particles_per_side,
+                                      square_particles_per_side),
+                                     (0.55, square_bottom_y),
+                                     density=square_density,
+                                     velocity=(1.0, 0.0))
+
+# ==========================================================================================
+# ==== Wall Boundary
+floor_length = 3.0
+floor_height = 0.03
+wall_density = 1000.0
+
+floor = RectangularTank(particle_spacing, (0.0, 0.0), (floor_length, floor_height),
+                        wall_density, n_layers=boundary_layers,
+                        min_coordinates=(-1.5, 0.0),
+                        faces=(false, false, true, false))
+
+boundary_model = BoundaryModelMonaghanKajtar(10.0, 1.0, particle_spacing,
+                                             floor.boundary.mass)
+boundary_system = WallBoundarySystem(floor.boundary, boundary_model)
+
+# ==========================================================================================
+# ==== Rigid Structures
+contact_model_frictionless = RigidContactModel(; normal_stiffness=2.0e5,
+                                               normal_damping=180.0,
+                                               contact_distance=2.0 * particle_spacing)
+
+contact_model_frictional = RigidContactModel(; normal_stiffness=2.0e5,
+                                             normal_damping=180.0,
+                                             static_friction_coefficient=0.6,
+                                             kinetic_friction_coefficient=0.4,
+                                             tangential_stiffness=1.0e5,
+                                             tangential_damping=150.0,
+                                             contact_distance=2.0 * particle_spacing)
+
+structure_system_frictionless = RigidBodySystem(square_frictionless;
+                                                contact_model=contact_model_frictionless,
+                                                acceleration=(0.0, -gravity),
+                                                particle_spacing=particle_spacing,
+                                                color_value=1)
+structure_system_frictional = RigidBodySystem(square_frictional;
+                                              contact_model=contact_model_frictional,
+                                              acceleration=(0.0, -gravity),
+                                              particle_spacing=particle_spacing,
+                                              color_value=2)
+
+# ==========================================================================================
+# ==== Simulation
+semi = Semidiscretization(structure_system_frictionless, structure_system_frictional,
+                          boundary_system)
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=50)
+saving_callback = SolutionSavingCallback(dt=0.02, output_directory="out", prefix="")
+
+callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback())
+
+sol = solve(ode, RDPK3SpFSAL49(),
+            abstol=1e-6,
+            reltol=1e-4,
+            dtmax=1e-3,
+            save_everystep=false, callback=callbacks);
diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl
index 42b85fd669..c85fb3119d 100644
--- a/src/callbacks/post_process.jl
+++ b/src/callbacks/post_process.jl
@@ -159,7 +159,8 @@ function (pp::PostprocessCallback)(integrator)
         # still have the values from the last stage of the previous step if not updated here.
         @trixi_timeit timer() "update systems and nhs" begin
             # Don't create sub-timers here to avoid cluttering the timer output
-            @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t)
+            @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t;
+                                                     reset_interaction_caches=false)
         end
 
         # Transfer to CPU if data is on the GPU. Do nothing if already on CPU.
diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl
index ab76dfacc2..cd5df79952 100644
--- a/src/callbacks/update.jl
+++ b/src/callbacks/update.jl
@@ -96,6 +96,10 @@ function (update_callback!::UpdateCallback)(integrator)
             update_particle_packing(system, v_ode, u_ode, semi, integrator)
         end
 
+        foreach_system(semi) do system
+            update_rigid_contact_eachstep!(system, v_ode, u_ode, semi, t, integrator)
+        end
+
         # This is only used by the particle packing system and should be removed in the future
         foreach_system(semi) do system
             update_transport_velocity!(system, v_ode, semi, integrator)
diff --git a/src/general/abstract_system.jl b/src/general/abstract_system.jl
index 10ccb6261f..0f9fa8517f 100644
--- a/src/general/abstract_system.jl
+++ b/src/general/abstract_system.jl
@@ -153,7 +153,7 @@ function update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t)
     return system
 end
 
-function update_final!(system, v, u, v_ode, u_ode, semi, t)
+function update_final!(system, v, u, v_ode, u_ode, semi, t; kwargs...)
     return system
 end
 
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index cd407858e4..07792a9456 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -491,7 +491,7 @@ end
 
 # Update the systems and neighborhood searches (NHS) for a simulation
 # before calling `interact!` to compute forces.
-function update_systems_and_nhs(v_ode, u_ode, semi, t)
+function update_systems_and_nhs(v_ode, u_ode, semi, t; reset_interaction_caches=true)
     # First update step before updating the NHS
     # (for example for writing the current coordinates in the TLSPH system)
     foreach_system(semi) do system
@@ -539,7 +539,8 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t)
         v = wrap_v(v_ode, system, semi)
         u = wrap_u(u_ode, system, semi)
 
-        update_final!(system, v, u, v_ode, u_ode, semi, t)
+        update_final!(system, v, u, v_ode, u_ode, semi, t;
+                      reset_interaction_caches)
     end
 end
 
diff --git a/src/io/io.jl b/src/io/io.jl
index e164098c4a..e5c52f29c7 100644
--- a/src/io/io.jl
+++ b/src/io/io.jl
@@ -128,6 +128,7 @@ function add_system_data!(system_data, system::RigidBodySystem)
     system_data["adhesion_coefficient"] = system.adhesion_coefficient
     system_data["color"] = system.cache.color
     add_system_data!(system_data, system.boundary_model)
+    add_system_data!(system_data, system.contact_model)
 end
 
 function add_system_data!(system_data, system::WallBoundarySystem)
@@ -219,6 +220,20 @@ function add_system_data!(system_data, contact_model::LinearContactModel)
     system_data["contact_model"]["normal_stiffness"] = contact_model.normal_stiffness
 end
 
+function add_system_data!(system_data, contact_model::RigidContactModel)
+    system_data["contact_model"] = Dict{String, Any}()
+    system_data["contact_model"]["model"] = type2string(contact_model)
+    system_data["contact_model"]["normal_stiffness"] = contact_model.normal_stiffness
+    system_data["contact_model"]["normal_damping"] = contact_model.normal_damping
+    system_data["contact_model"]["static_friction_coefficient"] = contact_model.static_friction_coefficient
+    system_data["contact_model"]["kinetic_friction_coefficient"] = contact_model.kinetic_friction_coefficient
+    system_data["contact_model"]["tangential_stiffness"] = contact_model.tangential_stiffness
+    system_data["contact_model"]["tangential_damping"] = contact_model.tangential_damping
+    system_data["contact_model"]["contact_distance"] = contact_model.contact_distance
+    system_data["contact_model"]["stick_velocity_tolerance"] = contact_model.stick_velocity_tolerance
+    system_data["contact_model"]["penetration_slop"] = contact_model.penetration_slop
+end
+
 function add_system_data!(system_data, state_equation::StateEquationCole)
     system_data["state_equation"] = Dict{String, Any}()
     system_data["state_equation"]["model"] = type2string(state_equation)
diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl
index 56eb4ae12d..d6b1a9a6e1 100644
--- a/src/io/write_vtk.jl
+++ b/src/io/write_vtk.jl
@@ -68,7 +68,8 @@ function trixi2vtk(dvdu_ode, vu_ode, semi, t; iter=nothing, output_directory="ou
     @trixi_timeit timer() "update systems" begin
         v_ode, u_ode = vu_ode.x
         # Don't create sub-timers here to avoid cluttering the timer output
-        @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t)
+        @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t;
+                                                 reset_interaction_caches=false)
     end
 
     filenames = system_names(systems)
@@ -438,6 +439,8 @@ function write2vtk!(vtk, v, u, t, system::RigidBodySystem)
     vtk["resultant_torque"] = [system.resultant_torque[]]
     vtk["angular_acceleration_force"] = [system.angular_acceleration_force[]]
     vtk["gyroscopic_acceleration"] = [system.gyroscopic_acceleration[]]
+    vtk["contact_count"] = [system.cache.contact_count[]]
+    vtk["max_contact_penetration"] = [system.cache.max_contact_penetration[]]
 
     write2vtk!(vtk, v, u, t, system.boundary_model, system)
 end
diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl
index 33f80f43a4..413c0d6a3d 100644
--- a/src/schemes/fluid/entropically_damped_sph/system.jl
+++ b/src/schemes/fluid/entropically_damped_sph/system.jl
@@ -306,7 +306,8 @@ function update_pressure!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_od
     compute_surface_delta_function!(system, system.surface_tension, semi)
 end
 
-function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t)
+function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t;
+                       kwargs...)
     (; surface_tension) = system
 
     # Surface normal of neighbor and boundary needs to have been calculated already
diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl
index f84849f096..514fea7009 100644
--- a/src/schemes/fluid/weakly_compressible_sph/system.jl
+++ b/src/schemes/fluid/weakly_compressible_sph/system.jl
@@ -338,7 +338,8 @@ function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_od
     return system
 end
 
-function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t)
+function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t;
+                       kwargs...)
     (; surface_tension) = system
 
     # Surface normal of neighbor and boundary needs to have been calculated already
diff --git a/src/schemes/structure/rigid_body/contact.jl b/src/schemes/structure/rigid_body/contact.jl
index 658be79c04..62b6211bd4 100644
--- a/src/schemes/structure/rigid_body/contact.jl
+++ b/src/schemes/structure/rigid_body/contact.jl
@@ -1 +1,2 @@
 include("contact_models.jl")
+include("contact_forces.jl")
diff --git a/src/schemes/structure/rigid_body/contact_forces.jl b/src/schemes/structure/rigid_body/contact_forces.jl
new file mode 100644
index 0000000000..8270b640a3
--- /dev/null
+++ b/src/schemes/structure/rigid_body/contact_forces.jl
@@ -0,0 +1,246 @@
+@inline function requires_update_callback(contact_model::RigidContactModel)
+    return contact_model.tangential_stiffness > 0 &&
+           contact_model.static_friction_coefficient > 0
+end
+
+function create_cache_contact_history(contact_model::RigidContactModel, ::Val{NDIMS},
+                                      ::Type{ELTYPE}) where {NDIMS, ELTYPE}
+    contact_tangential_displacement = requires_update_callback(contact_model) ?
+                                      Dict{RigidContactKey, SVector{NDIMS, ELTYPE}}() :
+                                      nothing
+
+    return (; contact_tangential_displacement)
+end
+
+@inline function requires_update_callback(system::RigidBodySystem)
+    return !isnothing(system.contact_model) &&
+           requires_update_callback(system.contact_model)
+end
+
+@inline function normal_friction_reference_force(contact_model::RigidContactModel,
+                                                 penetration, normal_velocity)
+    elastic_force = contact_model.normal_stiffness * penetration
+    damping_force = -contact_model.normal_damping * normal_velocity
+
+    return max(elastic_force + damping_force, zero(elastic_force))
+end
+
+function tangential_contact_force(contact_model::RigidContactModel,
+                                  tangential_displacement,
+                                  tangential_velocity,
+                                  normal_force_friction_reference)
+    force_trial = -contact_model.tangential_stiffness * tangential_displacement -
+                  contact_model.tangential_damping * tangential_velocity
+
+    trial_norm = norm(force_trial)
+    static_limit = contact_model.static_friction_coefficient *
+                   normal_force_friction_reference
+    if trial_norm <= static_limit
+        return force_trial
+    end
+
+    kinetic_limit = contact_model.kinetic_friction_coefficient *
+                    normal_force_friction_reference
+    kinetic_limit <= eps(eltype(tangential_velocity)) && return zero(force_trial)
+
+    tangential_speed = norm(tangential_velocity)
+    regularization_velocity = max(contact_model.stick_velocity_tolerance,
+                                  sqrt(eps(eltype(tangential_velocity))))
+
+    if tangential_speed > eps(eltype(tangential_velocity))
+        speed_factor = tanh(tangential_speed / regularization_velocity)
+        return -kinetic_limit * speed_factor * tangential_velocity / tangential_speed
+    end
+
+    if trial_norm > eps(eltype(tangential_velocity))
+        return -kinetic_limit * force_trial / trial_norm
+    end
+
+    return zero(force_trial)
+end
+
+update_rigid_contact_eachstep!(system, v_ode, u_ode, semi, t, integrator) = system
+
+function update_rigid_contact_eachstep!(system::RigidBodySystem{<:Any, <:Any, NDIMS},
+                                        v_ode, u_ode, semi, t,
+                                        integrator) where {NDIMS}
+    requires_update_callback(system) || return system
+
+    v_system = wrap_v(v_ode, system, semi)
+    u_system = wrap_u(u_ode, system, semi)
+    active_contact_keys = Set{RigidContactKey}()
+
+    foreach_system(semi) do neighbor_system
+        neighbor_system === system && return
+        update_contact_history_pair!(system, neighbor_system, v_system, u_system,
+                                     v_ode, u_ode, semi, integrator.dt,
+                                     active_contact_keys)
+    end
+
+    contact_map = system.cache.contact_tangential_displacement
+    for key in collect(keys(contact_map))
+        key in active_contact_keys && continue
+        delete!(contact_map, key)
+    end
+
+    return system
+end
+
+update_contact_history_pair!(system, neighbor_system, v_system, u_system, v_ode, u_ode,
+                             semi, dt, active_contact_keys) = active_contact_keys
+
+function update_contact_history_pair!(system::RigidBodySystem{<:Any, <:Any, NDIMS},
+                                      neighbor_system::WallBoundarySystem,
+                                      v_system, u_system,
+                                      v_ode, u_ode,
+                                      semi, dt,
+                                      active_contact_keys) where {NDIMS}
+    contact_model = system.contact_model
+    isnothing(contact_model) && return active_contact_keys
+
+    set_zero!(system.cache.contact_manifold_count)
+    set_zero!(system.cache.contact_manifold_weight_sum)
+    set_zero!(system.cache.contact_manifold_penetration_sum)
+    set_zero!(system.cache.contact_manifold_normal_sum)
+    set_zero!(system.cache.contact_manifold_wall_velocity_sum)
+
+    v_neighbor = wrap_v(v_ode, neighbor_system, semi)
+    u_neighbor = wrap_u(u_ode, neighbor_system, semi)
+    system_coords = current_coordinates(u_system, system)
+    neighbor_coords = current_coordinates(u_neighbor, neighbor_system)
+
+    foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, semi;
+                           points=each_integrated_particle(system),
+                           parallelization_backend=SerialBackend()) do particle, neighbor,
+                                                                       pos_diff, distance
+        accumulate_wall_contact_pair!(system, v_neighbor, neighbor_system,
+                                      particle, neighbor, pos_diff, distance,
+                                      contact_model)
+    end
+
+    neighbor_system_index = system_indices(neighbor_system, semi)
+    ELTYPE = eltype(system)
+    zero_tangential = zero(SVector{NDIMS, ELTYPE})
+
+    for particle in each_integrated_particle(system)
+        n_manifolds = system.cache.contact_manifold_count[particle]
+        n_manifolds == 0 && continue
+
+        particle_velocity = current_velocity(v_system, system, particle)
+
+        for manifold_index in 1:n_manifolds
+            weight_sum = system.cache.contact_manifold_weight_sum[manifold_index, particle]
+            weight_sum <= eps(ELTYPE) && continue
+
+            normal = extract_svector(system.cache.contact_manifold_normal_sum, Val(NDIMS),
+                                     manifold_index, particle) / weight_sum
+            normal_norm = norm(normal)
+            normal_norm <= eps(ELTYPE) && continue
+            normal /= normal_norm
+
+            wall_velocity = extract_svector(system.cache.contact_manifold_wall_velocity_sum,
+                                            Val(NDIMS), manifold_index, particle) /
+                            weight_sum
+            penetration_effective = system.cache.contact_manifold_penetration_sum[manifold_index,
+                                                                                  particle] /
+                                    weight_sum
+            relative_velocity = particle_velocity - wall_velocity
+            normal_velocity = dot(relative_velocity, normal)
+            tangential_velocity = relative_velocity - normal_velocity * normal
+
+            contact_key = wall_contact_key(neighbor_system_index, particle, manifold_index)
+            push!(active_contact_keys, contact_key)
+            update_contact_tangential_history!(system, contact_key, tangential_velocity,
+                                               normal, penetration_effective,
+                                               normal_velocity, dt, contact_model,
+                                               zero_tangential)
+        end
+    end
+
+    return active_contact_keys
+end
+
+function update_contact_history_pair!(system::RigidBodySystem{<:Any, <:Any, NDIMS},
+                                      neighbor_system::RigidBodySystem,
+                                      v_system, u_system,
+                                      v_ode, u_ode,
+                                      semi, dt,
+                                      active_contact_keys) where {NDIMS}
+    contact_model = system.contact_model
+    neighbor_contact_model = neighbor_system.contact_model
+    if isnothing(contact_model) || isnothing(neighbor_contact_model)
+        return active_contact_keys
+    end
+
+    v_neighbor = wrap_v(v_ode, neighbor_system, semi)
+    u_neighbor = wrap_u(u_ode, neighbor_system, semi)
+    system_coords = current_coordinates(u_system, system)
+    neighbor_coords = current_coordinates(u_neighbor, neighbor_system)
+
+    neighbor_system_index = system_indices(neighbor_system, semi)
+    ELTYPE = eltype(system)
+    zero_tangential = zero(SVector{NDIMS, ELTYPE})
+
+    foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, semi;
+                           points=each_integrated_particle(system),
+                           parallelization_backend=SerialBackend()) do particle, neighbor,
+                                                                       pos_diff, distance
+        distance <= eps(ELTYPE) && return
+
+        penetration = max(contact_model.contact_distance,
+                          neighbor_contact_model.contact_distance) - distance
+        penetration_effective = penetration - contact_model.penetration_slop
+        penetration_effective <= 0 && return
+
+        normal = pos_diff / distance
+        particle_velocity = current_velocity(v_system, system, particle)
+        neighbor_velocity = current_velocity(v_neighbor, neighbor_system, neighbor)
+        relative_velocity = particle_velocity - neighbor_velocity
+        normal_velocity = dot(relative_velocity, normal)
+        tangential_velocity = relative_velocity - normal_velocity * normal
+
+        contact_key = rigid_rigid_contact_key(neighbor_system_index, particle, neighbor)
+        push!(active_contact_keys, contact_key)
+        update_contact_tangential_history!(system, contact_key, tangential_velocity,
+                                           normal, penetration_effective,
+                                           normal_velocity, dt, contact_model,
+                                           zero_tangential)
+    end
+
+    return active_contact_keys
+end
+
+function update_contact_tangential_history!(system::RigidBodySystem, contact_key,
+                                            tangential_velocity, normal,
+                                            penetration_effective, normal_velocity, dt,
+                                            contact_model::RigidContactModel,
+                                            zero_tangential)
+    contact_map = system.cache.contact_tangential_displacement
+    isnothing(contact_map) && return contact_map
+
+    dt_ = isfinite(dt) && dt > 0 ? convert(eltype(system), dt) : zero(eltype(system))
+    tangential_displacement = get(contact_map, contact_key, zero_tangential)
+    tangential_displacement += dt_ * tangential_velocity
+    tangential_displacement -= dot(tangential_displacement, normal) * normal
+
+    if contact_model.tangential_stiffness > eps(eltype(system))
+        normal_force_reference = normal_friction_reference_force(contact_model,
+                                                                 penetration_effective,
+                                                                 normal_velocity)
+        max_displacement = contact_model.static_friction_coefficient *
+                           normal_force_reference /
+                           contact_model.tangential_stiffness
+        displacement_norm = norm(tangential_displacement)
+
+        if displacement_norm > max_displacement &&
+           displacement_norm > eps(eltype(system))
+            tangential_displacement *= max_displacement / displacement_norm
+        end
+    else
+        tangential_displacement = zero_tangential
+    end
+
+    contact_map[contact_key] = tangential_displacement
+
+    return contact_map
+end
diff --git a/src/schemes/structure/rigid_body/contact_models.jl b/src/schemes/structure/rigid_body/contact_models.jl
index c32fe38b80..48bc1226e5 100644
--- a/src/schemes/structure/rigid_body/contact_models.jl
+++ b/src/schemes/structure/rigid_body/contact_models.jl
@@ -1,14 +1,65 @@
 abstract type AbstractRigidContactModel end
 
+@enum RigidContactKind::UInt8 begin
+    WallContact = 1
+    RigidRigidContact = 2
+end
+
+"""
+    RigidContactKey(neighbor_system_index, local_particle, contact_slot, contact_kind)
+
+Shared tangential-history key for rigid contact.
+
+`contact_slot` stores the wall-manifold index for rigid-wall contact and the neighbor
+particle index for rigid-rigid contact.
+"""
+struct RigidContactKey
+    neighbor_system_index::Int
+    local_particle::Int
+    contact_slot::Int
+    contact_kind::RigidContactKind
+end
+
+@inline wall_contact_key(neighbor_system_index, local_particle, manifold_index) =
+    RigidContactKey(neighbor_system_index, local_particle, manifold_index, WallContact)
+
+@inline rigid_rigid_contact_key(neighbor_system_index, local_particle, neighbor_particle) =
+    RigidContactKey(neighbor_system_index, local_particle, neighbor_particle,
+                    RigidRigidContact)
+
+@inline function Base.:(==)(lhs::RigidContactKey, rhs::RigidContactKey)
+    return lhs.neighbor_system_index == rhs.neighbor_system_index &&
+           lhs.local_particle == rhs.local_particle &&
+           lhs.contact_slot == rhs.contact_slot &&
+           lhs.contact_kind == rhs.contact_kind
+end
+
+@inline Base.isequal(lhs::RigidContactKey, rhs::RigidContactKey) = lhs == rhs
+
+@inline function Base.hash(key::RigidContactKey, h::UInt)
+    return hash((key.neighbor_system_index, key.local_particle,
+                 key.contact_slot, key.contact_kind), h)
+end
+
 """
     RigidContactModel(; normal_stiffness,
                       normal_damping=0.0,
-                      contact_distance=0.0)
-
-Basic rigid contact model stored on a rigid body.
-It is currently used for both rigid-wall and rigid-rigid contact.
-The contact force consists of a linear normal spring-dashpot contribution only.
-If `contact_distance == 0`, the particle spacing of the `RigidBodySystem` will be used as contact distance.
+                      static_friction_coefficient=nothing,
+                      kinetic_friction_coefficient=nothing,
+                      tangential_stiffness=nothing,
+                      tangential_damping=nothing,
+                      contact_distance=0.0,
+                      stick_velocity_tolerance=nothing,
+                      penetration_slop=nothing)
+
+Shared rigid-contact model used by the active rigid-wall and rigid-rigid contact paths.
+Rigid-wall contact currently combines the linear normal spring-dashpot law with
+history-driven tangential friction updated through `UpdateCallback`, while rigid-rigid
+contact remains normal-only but shares the same runtime model and history-key design.
+
+If `contact_distance == 0`, the particle spacing of the `RigidBodySystem` will be used
+as contact distance when the model is adapted via
+`copy_contact_model(model, particle_spacing, ELTYPE)`.
 
 !!! warning "Experimental implementation"
     This is an experimental feature and may change in future releases.
@@ -16,28 +67,86 @@ If `contact_distance == 0`, the particle spacing of the `RigidBodySystem` will b
 struct RigidContactModel{ELTYPE <: Real} <: AbstractRigidContactModel
     normal_stiffness::ELTYPE
     normal_damping::ELTYPE
+    static_friction_coefficient::ELTYPE
+    kinetic_friction_coefficient::ELTYPE
+    tangential_stiffness::ELTYPE
+    tangential_damping::ELTYPE
     contact_distance::ELTYPE
+    stick_velocity_tolerance::ELTYPE
+    penetration_slop::ELTYPE
 end
 
 function RigidContactModel(; normal_stiffness,
                            normal_damping=0.0,
-                           contact_distance=0.0)
+                           static_friction_coefficient=nothing,
+                           kinetic_friction_coefficient=nothing,
+                           tangential_stiffness=nothing,
+                           tangential_damping=nothing,
+                           contact_distance=0.0,
+                           stick_velocity_tolerance=nothing,
+                           penetration_slop=nothing)
+    tangential_mode = !isnothing(static_friction_coefficient) ||
+                      !isnothing(kinetic_friction_coefficient) ||
+                      !isnothing(tangential_stiffness) ||
+                      !isnothing(tangential_damping)
+
+    static_friction_coefficient = something(static_friction_coefficient,
+                                            tangential_mode ? 0.5 : 0.0)
+    kinetic_friction_coefficient = something(kinetic_friction_coefficient,
+                                             tangential_mode ? 0.4 : 0.0)
+    tangential_stiffness = something(tangential_stiffness, 0.0)
+    tangential_damping = something(tangential_damping, 0.0)
+    stick_velocity_tolerance = something(stick_velocity_tolerance, 1.0e-6)
+    penetration_slop = something(penetration_slop, 0.0)
     ELTYPE = promote_type(typeof(normal_stiffness),
                           typeof(normal_damping),
-                          typeof(contact_distance))
+                          typeof(static_friction_coefficient),
+                          typeof(kinetic_friction_coefficient),
+                          typeof(tangential_stiffness),
+                          typeof(tangential_damping),
+                          typeof(contact_distance),
+                          typeof(stick_velocity_tolerance),
+                          typeof(penetration_slop))
 
     normal_stiffness_ = convert(ELTYPE, normal_stiffness)
     normal_damping_ = convert(ELTYPE, normal_damping)
+    static_friction_coefficient_ = convert(ELTYPE, static_friction_coefficient)
+    kinetic_friction_coefficient_ = convert(ELTYPE, kinetic_friction_coefficient)
+    tangential_stiffness_ = convert(ELTYPE, tangential_stiffness)
+    tangential_damping_ = convert(ELTYPE, tangential_damping)
     contact_distance_ = convert(ELTYPE, contact_distance)
+    stick_velocity_tolerance_ = convert(ELTYPE, stick_velocity_tolerance)
+    penetration_slop_ = convert(ELTYPE, penetration_slop)
 
     normal_stiffness_ > 0 ||
         throw(ArgumentError("`normal_stiffness` must be positive"))
     normal_damping_ >= 0 ||
         throw(ArgumentError("`normal_damping` must be non-negative"))
+    static_friction_coefficient_ >= 0 ||
+        throw(ArgumentError("`static_friction_coefficient` must be non-negative"))
+    kinetic_friction_coefficient_ >= 0 ||
+        throw(ArgumentError("`kinetic_friction_coefficient` must be non-negative"))
+    kinetic_friction_coefficient_ <= static_friction_coefficient_ ||
+        throw(ArgumentError("`kinetic_friction_coefficient` must be <= `static_friction_coefficient`"))
+    tangential_stiffness_ >= 0 ||
+        throw(ArgumentError("`tangential_stiffness` must be non-negative"))
+    tangential_damping_ >= 0 ||
+        throw(ArgumentError("`tangential_damping` must be non-negative"))
     contact_distance_ >= 0 ||
         throw(ArgumentError("`contact_distance` must be non-negative"))
-
-    return RigidContactModel(normal_stiffness_, normal_damping_, contact_distance_)
+    stick_velocity_tolerance_ >= 0 ||
+        throw(ArgumentError("`stick_velocity_tolerance` must be non-negative"))
+    penetration_slop_ >= 0 ||
+        throw(ArgumentError("`penetration_slop` must be non-negative"))
+
+    return RigidContactModel(normal_stiffness_, normal_damping_,
+                             static_friction_coefficient_,
+                             kinetic_friction_coefficient_,
+                             tangential_stiffness_,
+                             tangential_damping_,
+                             contact_distance_,
+                             stick_velocity_tolerance_,
+                             penetration_slop_)
 end
 
 function copy_contact_model(model::RigidContactModel, particle_spacing,
@@ -52,7 +161,18 @@ function copy_contact_model(model::RigidContactModel, particle_spacing,
 
     return RigidContactModel(; normal_stiffness=convert(ELTYPE, model.normal_stiffness),
                              normal_damping=convert(ELTYPE, model.normal_damping),
-                             contact_distance)
+                             static_friction_coefficient=convert(ELTYPE,
+                                                                 model.static_friction_coefficient),
+                             kinetic_friction_coefficient=convert(ELTYPE,
+                                                                  model.kinetic_friction_coefficient),
+                             tangential_stiffness=convert(ELTYPE,
+                                                          model.tangential_stiffness),
+                             tangential_damping=convert(ELTYPE,
+                                                        model.tangential_damping),
+                             contact_distance,
+                             stick_velocity_tolerance=convert(ELTYPE,
+                                                              model.stick_velocity_tolerance),
+                             penetration_slop=convert(ELTYPE, model.penetration_slop))
 end
 
 function contact_time_step(system::RigidBodySystem, semi)
@@ -108,6 +228,12 @@ function Base.show(io::IO, model::RigidContactModel)
     print(io, "RigidContactModel(")
     print(io, "normal_stiffness=", model.normal_stiffness)
     print(io, ", normal_damping=", model.normal_damping)
+    print(io, ", static_friction_coefficient=", model.static_friction_coefficient)
+    print(io, ", kinetic_friction_coefficient=", model.kinetic_friction_coefficient)
+    print(io, ", tangential_stiffness=", model.tangential_stiffness)
+    print(io, ", tangential_damping=", model.tangential_damping)
     print(io, ", contact_distance=", model.contact_distance)
+    print(io, ", stick_velocity_tolerance=", model.stick_velocity_tolerance)
+    print(io, ", penetration_slop=", model.penetration_slop)
     print(io, ")")
 end
diff --git a/src/schemes/structure/rigid_body/rhs.jl b/src/schemes/structure/rigid_body/rhs.jl
index c2bdc800b2..ce77accde4 100644
--- a/src/schemes/structure/rigid_body/rhs.jl
+++ b/src/schemes/structure/rigid_body/rhs.jl
@@ -205,6 +205,13 @@ function interact!(dv, v_particle_system, u_particle_system,
 
     NDIMS = ndims(particle_system)
     ELTYPE = eltype(particle_system)
+    zero_tangential = zero(SVector{NDIMS, ELTYPE})
+    contact_map = particle_system.cache.contact_tangential_displacement
+    neighbor_system_index = system_indices(neighbor_system, semi)
+    set_zero!(particle_system.cache.contact_count_per_particle)
+    set_zero!(particle_system.cache.max_contact_penetration_per_particle)
+    contact_count_per_particle = particle_system.cache.contact_count_per_particle
+    max_contact_penetration_per_particle = particle_system.cache.max_contact_penetration_per_particle
 
     # First gather all penetrating wall neighbors into a small set of contact manifolds per
     # rigid particle. This avoids applying one noisy contact force per wall particle.
@@ -225,6 +232,8 @@ function interact!(dv, v_particle_system, u_particle_system,
     # Apply one force contribution per manifold using the averaged normal, penetration, and
     # wall velocity stored in the cache.
     @threaded semi for particle in each_integrated_particle(particle_system)
+        local_contact_count = 0
+        local_max_contact_penetration = zero(ELTYPE)
         n_manifolds = particle_system.cache.contact_manifold_count[particle]
         if n_manifolds > 0
             v_particle = current_velocity(v_particle_system, particle_system, particle)
@@ -256,28 +265,48 @@ function interact!(dv, v_particle_system, u_particle_system,
 
                         relative_velocity = v_particle - v_boundary
                         normal_velocity = dot(relative_velocity, normal)
-
-                        # Only the normal spring-dashpot part is kept in the basic collision.
-                        elastic_force = contact_model.normal_stiffness *
-                                        penetration_effective
-                        damping_force = -contact_model.normal_damping * normal_velocity
-                        normal_force_magnitude = max(elastic_force + damping_force,
-                                                     zero(ELTYPE))
+                        tangential_velocity = relative_velocity - normal_velocity * normal
+                        normal_force_magnitude = normal_friction_reference_force(contact_model,
+                                                                                penetration_effective,
+                                                                                normal_velocity)
 
                         if normal_force_magnitude > 0
-                            interaction_force = normal_force_magnitude * normal
+                            tangential_displacement = isnothing(contact_map) ?
+                                                      zero_tangential :
+                                                      get(contact_map,
+                                                          wall_contact_key(neighbor_system_index,
+                                                                           particle,
+                                                                           manifold_index),
+                                                          zero_tangential)
+                            tangential_force = tangential_contact_force(contact_model,
+                                                                        tangential_displacement,
+                                                                        tangential_velocity,
+                                                                        normal_force_magnitude)
+                            interaction_force = normal_force_magnitude * normal +
+                                                tangential_force
 
                             for dim in eachindex(interaction_force)
                                 particle_system.force_per_particle[dim,
                                                                    particle] += interaction_force[dim]
                             end
+
+                            local_contact_count += 1
+                            local_max_contact_penetration = max(local_max_contact_penetration,
+                                                                penetration_effective)
                         end
                     end
                 end
             end
         end
+
+        contact_count_per_particle[particle] = local_contact_count
+        max_contact_penetration_per_particle[particle] = local_max_contact_penetration
     end
 
+    particle_system.cache.contact_count[] += sum(contact_count_per_particle)
+    particle_system.cache.max_contact_penetration[] = max(particle_system.cache.max_contact_penetration[],
+                                                          maximum(max_contact_penetration_per_particle))
+
     return dv
 end
 
@@ -297,7 +326,8 @@ end
     distance <= eps(ELTYPE) && return particle_system
 
     penetration = contact_model.contact_distance - distance
-    penetration <= 0 && return particle_system
+    penetration_effective = penetration - contact_model.penetration_slop
+    penetration_effective <= 0 && return particle_system
 
     normal = pos_diff / distance
     wall_velocity = current_velocity(v_neighbor_system, neighbor_system, neighbor)
@@ -323,7 +353,8 @@ end
                                                       normal,
                                                       normal_merge_cos)
     accumulate_contact_manifold_sums!(particle_system.cache, particle, manifold_index,
-                                      contact_weight, normal, wall_velocity, penetration)
+                                      contact_weight, normal, wall_velocity,
+                                      penetration_effective)
 
     return particle_system
 end
@@ -430,6 +461,14 @@ function interact!(dv, v_particle_system, u_particle_system,
     ELTYPE = eltype(particle_system)
     system_coords = current_coordinates(u_particle_system, particle_system)
     neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
+    set_zero!(particle_system.cache.contact_count_per_particle)
+    set_zero!(particle_system.cache.max_contact_penetration_per_particle)
+    contact_count_per_particle = particle_system.cache.contact_count_per_particle
+    max_contact_penetration_per_particle = particle_system.cache.max_contact_penetration_per_particle
+    pair_normal_stiffness = (contact_model.normal_stiffness +
+                             neighbor_contact_model.normal_stiffness) / 2
+    pair_normal_damping = (contact_model.normal_damping +
+                           neighbor_contact_model.normal_damping) / 2
 
     foreach_point_neighbor(particle_system, neighbor_system, system_coords, neighbor_coords,
                            semi;
@@ -453,10 +492,8 @@ function interact!(dv, v_particle_system, u_particle_system,
         relative_velocity = particle_velocity - neighbor_velocity
         normal_velocity = dot(relative_velocity, normal)
 
-        elastic_force = (contact_model.normal_stiffness +
-                         neighbor_contact_model.normal_stiffness) / 2 * penetration
-        damping_force = -(contact_model.normal_damping +
-                          neighbor_contact_model.normal_damping) / 2 * normal_velocity
+        elastic_force = pair_normal_stiffness * penetration
+        damping_force = -pair_normal_damping * normal_velocity
         normal_force_magnitude = max(elastic_force + damping_force, zero(ELTYPE))
         normal_force_magnitude <= 0 && return dv
 
@@ -465,8 +502,18 @@ function interact!(dv, v_particle_system, u_particle_system,
         for dim in 1:ndims(particle_system)
             particle_system.force_per_particle[dim, particle] += interaction_force[dim]
         end
+
+        # `foreach_point_neighbor` parallelizes the outer loop over `points`, i.e. particles.
+        # This makes these per-particle reductions race-free under the regular backends.
+        contact_count_per_particle[particle] += 1
+        max_contact_penetration_per_particle[particle] = max(max_contact_penetration_per_particle[particle],
+                                                             penetration)
     end
 
+    particle_system.cache.contact_count[] += sum(contact_count_per_particle)
+    particle_system.cache.max_contact_penetration[] = max(particle_system.cache.max_contact_penetration[],
+                                                          maximum(max_contact_penetration_per_particle))
+
     return dv
 end
 
diff --git a/src/schemes/structure/rigid_body/system.jl b/src/schemes/structure/rigid_body/system.jl
index cdccaaa6d2..072bd346fa 100644
--- a/src/schemes/structure/rigid_body/system.jl
+++ b/src/schemes/structure/rigid_body/system.jl
@@ -119,7 +119,9 @@ function RigidBodySystem(initial_condition; boundary_model=nothing,
         inverse_inertia = Ref(zero(SMatrix{3, 3, ELTYPE, 9}))
     end
 
-    cache = (;
+    cache = (; create_cache_contact_history(contact_model_, Val(NDIMS), ELTYPE)...,
+             contact_count=Ref(0),
+             max_contact_penetration=Ref(zero(ELTYPE)),
              create_cache_contact_manifold(contact_model_, Val(NDIMS), ELTYPE,
                                            nparticles(initial_condition),
                                            max_manifolds_)...,
@@ -149,15 +151,22 @@ function create_cache_contact_manifold(::Nothing, ::Val{NDIMS}, ELTYPE,
     return (;)
 end
 
-# Allocate per-particle manifold scratch arrays for rigid-wall contact.
+function create_cache_contact_history(contact_model, ::Val{NDIMS},
+                                      ::Type{ELTYPE}) where {NDIMS, ELTYPE}
+    return (; contact_tangential_displacement=nothing)
+end
+
+# Allocate per-particle scratch arrays for rigid contact.
 #
-# The cache shape is `[dimension, manifold, particle]` for vector-valued sums and
+# The manifold cache shape is `[dimension, manifold, particle]` for vector-valued sums and
 # `[manifold, particle]` for scalar sums. It is rebuilt for each rigid-wall system pair in
-# the RHS and therefore acts purely as transient manifold assembly storage; the persistent
-# collision effect is the force accumulated in `force_per_particle`.
+# the RHS and therefore acts purely as transient manifold assembly storage. The per-particle
+# diagnostic scratch is reused for both rigid-wall and rigid-rigid contact reductions.
 function create_cache_contact_manifold(contact_model, ::Val{NDIMS}, ELTYPE,
                                        n_particles, max_manifolds) where {NDIMS}
-    return (; contact_manifold_count=zeros(Int, n_particles),
+    return (; contact_count_per_particle=zeros(Int, n_particles),
+            max_contact_penetration_per_particle=zeros(ELTYPE, n_particles),
+            contact_manifold_count=zeros(Int, n_particles),
             contact_manifold_weight_sum=zeros(ELTYPE, max_manifolds, n_particles),
             contact_manifold_penetration_sum=zeros(ELTYPE, max_manifolds, n_particles),
             contact_manifold_normal_sum=zeros(ELTYPE, NDIMS, max_manifolds, n_particles),
@@ -371,7 +380,8 @@ function update_boundary_interpolation!(boundary_model, system::RigidBodySystem,
     return system
 end
 
-function update_final!(system::RigidBodySystem, v, u, v_ode, u_ode, semi, t)
+function update_final!(system::RigidBodySystem, v, u, v_ode, u_ode, semi, t;
+                       reset_interaction_caches=true)
     system_coords = current_coordinates(u, system)
     system_velocity = current_velocity(v, system)
     center_of_mass,
@@ -384,9 +394,12 @@ function update_final!(system::RigidBodySystem, v, u, v_ode, u_ode, semi, t)
                                                         center_of_mass_velocity;
                                                         relative_coordinates=system.relative_coordinates)
 
-    # Reset interaction caches before RHS assembly so pairwise rigid-fluid forces can
-    # accumulate from scratch and non-RHS update paths do not expose stale resultants.
-    set_zero!(system.force_per_particle)
+    if reset_interaction_caches
+        # Reset pairwise interaction accumulation before the next RHS assembly.
+        set_zero!(system.force_per_particle)
+        system.cache.contact_count[] = 0
+        system.cache.max_contact_penetration[] = zero(eltype(system))
+    end
 
     system.center_of_mass[] = center_of_mass
     system.center_of_mass_velocity[] = center_of_mass_velocity
@@ -540,6 +553,8 @@ function system_data(system::RigidBodySystem, dv_ode, du_ode, v_ode, u_ode, semi
     resultant_torque = system.resultant_torque[]
     angular_acceleration_force = system.angular_acceleration_force[]
     gyroscopic_acceleration = system.gyroscopic_acceleration[]
+    contact_count = system.cache.contact_count[]
+    max_contact_penetration = system.cache.max_contact_penetration[]
     relative_coordinates = system.relative_coordinates
 
     return (; coordinates, velocity, mass=system.mass,
@@ -549,6 +564,7 @@ function system_data(system::RigidBodySystem, dv_ode, du_ode, v_ode, u_ode, semi
             angular_velocity,
             resultant_force, resultant_torque,
             angular_acceleration_force, gyroscopic_acceleration,
+            contact_count, max_contact_penetration,
             density, pressure, acceleration)
 end
 
@@ -558,6 +574,7 @@ function available_data(::RigidBodySystem)
             :center_of_mass, :center_of_mass_velocity,
             :angular_velocity, :resultant_force, :resultant_torque,
             :angular_acceleration_force, :gyroscopic_acceleration,
+            :contact_count, :max_contact_penetration,
             :density, :pressure, :acceleration)
 end
 
diff --git a/test/examples/examples.jl b/test/examples/examples.jl
index 44856073af..f4d9c2b53a 100644
--- a/test/examples/examples.jl
+++ b/test/examples/examples.jl
@@ -74,6 +74,14 @@
             @test sol.retcode == ReturnCode.Success
             @test count_rhs_allocations(sol, semi) == 0
         end
+
+        @trixi_testset "structure/sliding_rigid_squares_friction_2d.jl" begin
+            @trixi_test_nowarn trixi_include(@__MODULE__,
+                                             joinpath(examples_dir(), "structure",
+                                                      "sliding_rigid_squares_friction_2d.jl"),
+                                             tspan=(0.0, 0.3))
+            @test sol.retcode == ReturnCode.Success
+        end
     end
 
     @testset verbose=true "FSI" begin
diff --git a/test/io/read_vtk.jl b/test/io/read_vtk.jl
index 264fb97c9e..adf779ce76 100644
--- a/test/io/read_vtk.jl
+++ b/test/io/read_vtk.jl
@@ -199,6 +199,10 @@
                   rigid_system.angular_acceleration_force[]
             @test only(Array(TrixiParticles.ReadVTK.get_data(field_data["gyroscopic_acceleration"]))) ==
                   rigid_system.gyroscopic_acceleration[]
+            @test only(Array(TrixiParticles.ReadVTK.get_data(field_data["contact_count"]))) ==
+                  rigid_system.cache.contact_count[]
+            @test only(Array(TrixiParticles.ReadVTK.get_data(field_data["max_contact_penetration"]))) ==
+                  rigid_system.cache.max_contact_penetration[]
         end
     end
 end
diff --git a/test/systems/rigid_system.jl b/test/systems/rigid_system.jl
index cae3b70a05..fd19eb71cb 100644
--- a/test/systems/rigid_system.jl
+++ b/test/systems/rigid_system.jl
@@ -344,7 +344,11 @@
         @test data.resultant_torque == 0.0
         @test data.angular_acceleration_force == 0.0
         @test data.gyroscopic_acceleration == 0.0
+        @test data.contact_count == 0
+        @test data.max_contact_penetration == 0.0
         @test data.relative_coordinates == rigid_system.relative_coordinates
+        @test :contact_count in fields
+        @test :max_contact_penetration in fields
         @test !(:local_coordinates in fields)
     end
 
@@ -815,6 +819,10 @@
 
         @test vec(force_after_forward_1[:, 1]) ≈ collect(expected_force)
         @test vec(rigid_system_2.force_per_particle[:, 1]) ≈ collect(-expected_force)
+        @test rigid_system_1.cache.contact_count[] == 1
+        @test rigid_system_2.cache.contact_count[] == 1
+        @test rigid_system_1.cache.max_contact_penetration[] ≈ pair_penetration
+        @test rigid_system_2.cache.max_contact_penetration[] ≈ pair_penetration
 
         @test TrixiParticles.compact_support(rigid_system_1, rigid_system_2) ≈
               pair_contact_distance
@@ -837,6 +845,30 @@
         @test TrixiParticles.calculate_dt(zero_velocity_ode, u_ode_rigid, 0.25,
                                           semi_rigid) ≈ 0.25 * pair_contact_dt
 
+        dv_ode_reset = zero(v_ode_rigid)
+        TrixiParticles.update_final!(rigid_system_1, v_rigid_1, u_rigid_1,
+                                     v_ode_rigid, u_ode_rigid, semi_rigid, 0.0)
+        TrixiParticles.update_final!(rigid_system_2, v_rigid_2, u_rigid_2,
+                                     v_ode_rigid, u_ode_rigid, semi_rigid, 0.0)
+        TrixiParticles.system_interaction!(dv_ode_reset, v_ode_rigid, u_ode_rigid,
+                                           semi_rigid)
+        @test rigid_system_1.cache.contact_count[] == 1
+        @test rigid_system_2.cache.contact_count[] == 1
+        @test rigid_system_1.cache.max_contact_penetration[] ≈ pair_penetration
+        @test rigid_system_2.cache.max_contact_penetration[] ≈ pair_penetration
+
+        TrixiParticles.set_zero!(dv_ode_reset)
+        TrixiParticles.update_final!(rigid_system_1, v_rigid_1, u_rigid_1,
+                                     v_ode_rigid, u_ode_rigid, semi_rigid, 0.0)
+        TrixiParticles.update_final!(rigid_system_2, v_rigid_2, u_rigid_2,
+                                     v_ode_rigid, u_ode_rigid, semi_rigid, 0.0)
+        TrixiParticles.system_interaction!(dv_ode_reset, v_ode_rigid, u_ode_rigid,
+                                           semi_rigid)
+        @test rigid_system_1.cache.contact_count[] == 1
+        @test rigid_system_2.cache.contact_count[] == 1
+        @test rigid_system_1.cache.max_contact_penetration[] ≈ pair_penetration
+        @test rigid_system_2.cache.max_contact_penetration[] ≈ pair_penetration
+
         dv_rigid_1 = TrixiParticles.wrap_v(dv_ode_rigid, rigid_system_1, semi_rigid)
         dv_rigid_2 = TrixiParticles.wrap_v(dv_ode_rigid, rigid_system_2, semi_rigid)
         TrixiParticles.finalize_interaction!(rigid_system_1, dv_rigid_1, v_rigid_1,
@@ -852,6 +884,29 @@
         @test dv_rigid_2[1, 1] ≈ -expected_force[1] / rigid_mass_2[1]
         @test dv_rigid_1[2, 1] ≈ 0.0
         @test dv_rigid_2[2, 1] ≈ 0.0
+
+        mktempdir() do tmp_dir
+            du_ode_rigid = zero(u_ode_rigid)
+            dvdu_ode_rigid = (; x=(dv_ode_rigid, du_ode_rigid))
+            vu_ode_rigid = (; x=(v_ode_rigid, u_ode_rigid))
+            trixi2vtk(dvdu_ode_rigid, vu_ode_rigid, semi_rigid, 0.0;
+                      output_directory=tmp_dir, iter=1)
+
+            contact_filename = TrixiParticles.system_names(semi_rigid.systems)[1]
+            vtk_contact = TrixiParticles.ReadVTK.VTKFile(joinpath(tmp_dir,
+                                                                  "$(contact_filename)_1.vtu"))
+            point_data_contact = TrixiParticles.ReadVTK.get_point_data(vtk_contact)
+
+            @test only(Array(TrixiParticles.ReadVTK.get_data(point_data_contact["contact_count"]))) ==
+                  rigid_system_1.cache.contact_count[]
+            @test only(Array(TrixiParticles.ReadVTK.get_data(point_data_contact["contact_count"]))) >
+                  0
+            @test only(Array(TrixiParticles.ReadVTK.get_data(point_data_contact["max_contact_penetration"]))) ≈
+                  rigid_system_1.cache.max_contact_penetration[]
+            @test only(Array(TrixiParticles.ReadVTK.get_data(point_data_contact["max_contact_penetration"]))) >
+                  0
+        end
+
         rigid_coordinates = reshape([0.0, 0.05], 2, 1)
         rigid_velocity = reshape([0.0, -1.0], 2, 1)
         rigid_mass = [1.0]
@@ -886,7 +941,34 @@
         @test runtime_model isa RigidContactModel
         @test runtime_model.normal_stiffness ≈ 2.0e4
         @test runtime_model.normal_damping ≈ 20.0
+        @test runtime_model.static_friction_coefficient ≈ 0.0
+        @test runtime_model.kinetic_friction_coefficient ≈ 0.0
+        @test runtime_model.tangential_stiffness ≈ 0.0
+        @test runtime_model.tangential_damping ≈ 0.0
         @test runtime_model.contact_distance ≈ 0.1
+        @test runtime_model.stick_velocity_tolerance ≈ 1.0e-6
+        @test runtime_model.penetration_slop ≈ 0.0
+
+        advanced_contact_model = RigidContactModel(; normal_stiffness=5.0,
+                                                   normal_damping=1.5,
+                                                   static_friction_coefficient=0.6,
+                                                   kinetic_friction_coefficient=0.4,
+                                                   tangential_stiffness=9.0,
+                                                   tangential_damping=2.5,
+                                                   contact_distance=0.0,
+                                                   stick_velocity_tolerance=1.0e-5,
+                                                   penetration_slop=0.01)
+        advanced_runtime_model = TrixiParticles.copy_contact_model(advanced_contact_model,
+                                                                   0.125, Float32)
+        @test advanced_runtime_model.normal_stiffness ≈ Float32(5.0)
+        @test advanced_runtime_model.normal_damping ≈ Float32(1.5)
+        @test advanced_runtime_model.static_friction_coefficient ≈ Float32(0.6)
+        @test advanced_runtime_model.kinetic_friction_coefficient ≈ Float32(0.4)
+        @test advanced_runtime_model.tangential_stiffness ≈ Float32(9.0)
+        @test advanced_runtime_model.tangential_damping ≈ Float32(2.5)
+        @test advanced_runtime_model.contact_distance ≈ Float32(0.125)
+        @test advanced_runtime_model.stick_velocity_tolerance ≈ Float32(1.0e-5)
+        @test advanced_runtime_model.penetration_slop ≈ Float32(0.01)
 
         spacing_scaled_model = RigidContactModel(; normal_stiffness=5.0)
         spacing_scaled_runtime = TrixiParticles.copy_contact_model(spacing_scaled_model,
@@ -899,10 +981,26 @@
                                                      normal_damping=-1.0)
         @test_throws ArgumentError RigidContactModel(; normal_stiffness=1.0,
                                                      contact_distance=-1.0)
+        @test_throws ArgumentError RigidContactModel(; normal_stiffness=1.0,
+                                                     static_friction_coefficient=-0.1)
+        @test_throws ArgumentError RigidContactModel(; normal_stiffness=1.0,
+                                                     static_friction_coefficient=0.3,
+                                                     kinetic_friction_coefficient=0.4)
+        @test_throws ArgumentError RigidContactModel(; normal_stiffness=1.0,
+                                                     tangential_stiffness=-1.0)
+        @test_throws ArgumentError RigidContactModel(; normal_stiffness=1.0,
+                                                     tangential_damping=-1.0)
+        @test_throws ArgumentError RigidContactModel(; normal_stiffness=1.0,
+                                                     stick_velocity_tolerance=-1.0)
+        @test_throws ArgumentError RigidContactModel(; normal_stiffness=1.0,
+                                                     penetration_slop=-1.0)
 
         rigid_system = RigidBodySystem(rigid_ic;
                                        acceleration=(0.0, 0.0),
                                        contact_model=contact_model)
+        rigid_system_advanced = RigidBodySystem(rigid_ic;
+                                                acceleration=(0.0, 0.0),
+                                                contact_model=advanced_contact_model)
         rigid_system_with_boundary = RigidBodySystem(rigid_ic;
                                                      acceleration=(0.0, 0.0),
                                                      boundary_model=boundary_model,
@@ -919,6 +1017,10 @@
         @test rigid_system.contact_model.normal_stiffness ≈ contact_model.normal_stiffness
         @test rigid_system.contact_model.normal_damping ≈ contact_model.normal_damping
         @test rigid_system.contact_model.contact_distance ≈ contact_model.contact_distance
+        @test !TrixiParticles.requires_update_callback(rigid_system)
+        @test isnothing(rigid_system.cache.contact_tangential_displacement)
+        @test TrixiParticles.requires_update_callback(rigid_system_advanced)
+        @test rigid_system_advanced.cache.contact_tangential_displacement isa Dict
         @test size(rigid_system_custom_manifolds.cache.contact_manifold_weight_sum, 1) == 3
         @test TrixiParticles.compact_support(rigid_system, boundary_system) ≈
               contact_model.contact_distance
@@ -931,6 +1033,18 @@
                                                    contact_model=contact_model,
                                                    max_manifolds=0)
 
+        system_meta_data = Dict{String, Any}()
+        TrixiParticles.add_system_data!(system_meta_data, rigid_system_advanced)
+        @test system_meta_data["contact_model"]["normal_stiffness"] ≈ 5.0
+        @test system_meta_data["contact_model"]["normal_damping"] ≈ 1.5
+        @test system_meta_data["contact_model"]["static_friction_coefficient"] ≈ 0.6
+        @test system_meta_data["contact_model"]["kinetic_friction_coefficient"] ≈ 0.4
+        @test system_meta_data["contact_model"]["tangential_stiffness"] ≈ 9.0
+        @test system_meta_data["contact_model"]["tangential_damping"] ≈ 2.5
+        @test system_meta_data["contact_model"]["contact_distance"] ≈ 0.1
+        @test system_meta_data["contact_model"]["stick_velocity_tolerance"] ≈ 1.0e-5
+        @test system_meta_data["contact_model"]["penetration_slop"] ≈ 0.01
+
         semi = Semidiscretization(rigid_system, boundary_system)
         ode = semidiscretize(semi, (0.0, 0.01))
         v_ode, u_ode = ode.u0.x
@@ -944,5 +1058,190 @@
                                              dv_ode, v_ode, u_ode, semi)
 
         @test dv[2, 1] > 0
+        @test rigid_system.cache.contact_count[] == 1
+        @test rigid_system.cache.max_contact_penetration[] ≈ 0.05
+    end
+
+    @trixi_testset "Rigid Contact History" begin
+        rigid_coordinates = reshape([0.0, 0.05], 2, 1)
+        rigid_velocity = reshape([1.0, -1.0], 2, 1)
+        rigid_mass = [1.0]
+        rigid_density = [1000.0]
+        rigid_ic = InitialCondition(; coordinates=rigid_coordinates,
+                                    velocity=rigid_velocity,
+                                    mass=rigid_mass,
+                                    density=rigid_density,
+                                    particle_spacing=0.1)
+
+        boundary_coordinates = reshape([0.0, 0.0], 2, 1)
+        boundary_mass = [1.0]
+        boundary_density = [1000.0]
+        boundary_ic = InitialCondition(; coordinates=boundary_coordinates,
+                                       mass=boundary_mass,
+                                       density=boundary_density,
+                                       particle_spacing=0.1)
+
+        smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+        smoothing_length = 0.15
+        boundary_model = BoundaryModelDummyParticles(boundary_density, boundary_mass,
+                                                     SummationDensity(),
+                                                     smoothing_kernel,
+                                                     smoothing_length)
+        boundary_system = WallBoundarySystem(boundary_ic, boundary_model)
+
+        history_model = RigidContactModel(; normal_stiffness=2.0e4,
+                                          normal_damping=20.0,
+                                          static_friction_coefficient=0.6,
+                                          kinetic_friction_coefficient=0.4,
+                                          tangential_stiffness=1.0e4,
+                                          tangential_damping=5.0,
+                                          contact_distance=0.1,
+                                          stick_velocity_tolerance=1.0e-6)
+        rigid_system = RigidBodySystem(rigid_ic;
+                                       acceleration=(0.0, 0.0),
+                                       contact_model=history_model)
+
+        @test TrixiParticles.requires_update_callback(rigid_system)
+        @test rigid_system.cache.contact_tangential_displacement isa Dict
+
+        semi = Semidiscretization(rigid_system, boundary_system)
+        ode = semidiscretize(semi, (0.0, 0.01))
+        v_ode, u_ode = ode.u0.x
+        dv_ode = zero(v_ode)
+        update_error = try
+            TrixiParticles.kick!(dv_ode, v_ode, u_ode, semi, 0.0)
+            nothing
+        catch err
+            err
+        end
+        @test update_error isa ArgumentError
+        @test occursin("`UpdateCallback` is required for `RigidBodySystem`",
+                       sprint(showerror, update_error))
+
+        TrixiParticles.update_rigid_contact_eachstep!(rigid_system, v_ode, u_ode, semi, 0.0,
+                                                      (; dt=1.0e-3))
+
+        contact_map = rigid_system.cache.contact_tangential_displacement
+        @test length(contact_map) == 1
+        contact_key = first(keys(contact_map))
+        tangential_displacement = contact_map[contact_key]
+        @test contact_key.contact_kind == TrixiParticles.WallContact
+        @test contact_key.local_particle == 1
+        @test norm(tangential_displacement) > 0
+
+        v_rigid = TrixiParticles.wrap_v(v_ode, rigid_system, semi)
+        u_rigid = TrixiParticles.wrap_u(u_ode, rigid_system, semi)
+        dv = TrixiParticles.wrap_v(dv_ode, rigid_system, semi)
+        TrixiParticles.update_final!(rigid_system, v_rigid, u_rigid, v_ode, u_ode, semi,
+                                     0.0)
+        TrixiParticles.interact!(dv_ode, v_ode, u_ode, rigid_system, boundary_system, semi)
+        TrixiParticles.finalize_interaction!(rigid_system, dv, v_rigid, u_rigid,
+                                             dv_ode, v_ode, u_ode, semi)
+
+        @test rigid_system.force_per_particle[1, 1] < 0.0
+        @test rigid_system.force_per_particle[2, 1] > 0.0
+        @test dv[1, 1] < 0.0
+        @test dv[2, 1] > 0.0
+
+        u_rigid[2, 1] = 0.2
+        TrixiParticles.update_rigid_contact_eachstep!(rigid_system, v_ode, u_ode, semi, 0.0,
+                                                      (; dt=1.0e-3))
+        @test isempty(rigid_system.cache.contact_tangential_displacement)
+
+        rigid_coordinates_1 = reshape([0.0, 0.0], 2, 1)
+        rigid_coordinates_2 = reshape([0.08, 0.0], 2, 1)
+        rigid_velocity_1 = reshape([1.0, 0.5], 2, 1)
+        rigid_velocity_2 = reshape([-0.5, -0.25], 2, 1)
+        rigid_ic_1 = InitialCondition(; coordinates=rigid_coordinates_1,
+                                      velocity=rigid_velocity_1,
+                                      mass=[2.0],
+                                      density=rigid_density,
+                                      particle_spacing=0.1)
+        rigid_ic_2 = InitialCondition(; coordinates=rigid_coordinates_2,
+                                      velocity=rigid_velocity_2,
+                                      mass=rigid_mass,
+                                      density=rigid_density,
+                                      particle_spacing=0.1)
+
+        rigid_contact_model_1 = RigidContactModel(; normal_stiffness=20.0,
+                                                  normal_damping=4.0,
+                                                  static_friction_coefficient=0.6,
+                                                  kinetic_friction_coefficient=0.4,
+                                                  tangential_stiffness=10.0,
+                                                  tangential_damping=2.0,
+                                                  contact_distance=0.1)
+        rigid_contact_model_2 = RigidContactModel(; normal_stiffness=30.0,
+                                                  normal_damping=8.0,
+                                                  static_friction_coefficient=0.5,
+                                                  kinetic_friction_coefficient=0.3,
+                                                  tangential_stiffness=8.0,
+                                                  tangential_damping=1.0,
+                                                  contact_distance=0.12)
+
+        rigid_system_1 = RigidBodySystem(rigid_ic_1;
+                                         acceleration=(0.0, 0.0),
+                                         contact_model=rigid_contact_model_1)
+        rigid_system_2 = RigidBodySystem(rigid_ic_2;
+                                         acceleration=(0.0, 0.0),
+                                         contact_model=rigid_contact_model_2)
+
+        semi_rigid = Semidiscretization(rigid_system_1, rigid_system_2)
+        ode_rigid = semidiscretize(semi_rigid, (0.0, 0.01))
+        v_ode_rigid, u_ode_rigid = ode_rigid.u0.x
+        dv_ode_rigid = zero(v_ode_rigid)
+        TrixiParticles.update_rigid_contact_eachstep!(rigid_system_1, v_ode_rigid,
+                                                      u_ode_rigid, semi_rigid, 0.0,
+                                                      (; dt=1.0e-3))
+        TrixiParticles.update_rigid_contact_eachstep!(rigid_system_2, v_ode_rigid,
+                                                      u_ode_rigid, semi_rigid, 0.0,
+                                                      (; dt=1.0e-3))
+
+        rigid_key_1 = first(keys(rigid_system_1.cache.contact_tangential_displacement))
+        rigid_key_2 = first(keys(rigid_system_2.cache.contact_tangential_displacement))
+        @test rigid_key_1.contact_kind == TrixiParticles.RigidRigidContact
+        @test rigid_key_2.contact_kind == TrixiParticles.RigidRigidContact
+
+        v_rigid_1 = TrixiParticles.wrap_v(v_ode_rigid, rigid_system_1, semi_rigid)
+        u_rigid_1 = TrixiParticles.wrap_u(u_ode_rigid, rigid_system_1, semi_rigid)
+        v_rigid_2 = TrixiParticles.wrap_v(v_ode_rigid, rigid_system_2, semi_rigid)
+        u_rigid_2 = TrixiParticles.wrap_u(u_ode_rigid, rigid_system_2, semi_rigid)
+        TrixiParticles.update_final!(rigid_system_1, v_rigid_1, u_rigid_1,
+                                     v_ode_rigid, u_ode_rigid, semi_rigid, 0.0)
+        TrixiParticles.update_final!(rigid_system_2, v_rigid_2, u_rigid_2,
+                                     v_ode_rigid, u_ode_rigid, semi_rigid, 0.0)
+        TrixiParticles.interact!(dv_ode_rigid, v_ode_rigid, u_ode_rigid,
+                                 rigid_system_1, rigid_system_2, semi_rigid)
+        TrixiParticles.interact!(dv_ode_rigid, v_ode_rigid, u_ode_rigid,
+                                 rigid_system_2, rigid_system_1, semi_rigid)
+
+        pair_contact_distance = max(rigid_contact_model_1.contact_distance,
+                                    rigid_contact_model_2.contact_distance)
+        pair_normal_stiffness = (rigid_contact_model_1.normal_stiffness +
+                                 rigid_contact_model_2.normal_stiffness) / 2
+        pair_normal_damping = (rigid_contact_model_1.normal_damping +
+                               rigid_contact_model_2.normal_damping) / 2
+        pair_penetration = pair_contact_distance - 0.08
+        normal_velocity = -1.5
+        expected_force_magnitude = pair_normal_stiffness * pair_penetration -
+                                   pair_normal_damping * normal_velocity
+
+        @test rigid_system_1.force_per_particle[1, 1] ≈ -expected_force_magnitude
+        @test rigid_system_1.force_per_particle[2, 1] ≈ 0.0
+        @test rigid_system_2.force_per_particle[1, 1] ≈ expected_force_magnitude
+        @test rigid_system_2.force_per_particle[2, 1] ≈ 0.0
+        @test rigid_system_1.cache.contact_count[] == 1
+        @test rigid_system_2.cache.contact_count[] == 1
+        @test rigid_system_1.cache.max_contact_penetration[] ≈ pair_penetration
+        @test rigid_system_2.cache.max_contact_penetration[] ≈ pair_penetration
+
+        u_rigid_2[1, 1] = 0.5
+        TrixiParticles.update_rigid_contact_eachstep!(rigid_system_1, v_ode_rigid,
+                                                      u_ode_rigid, semi_rigid, 0.0,
+                                                      (; dt=1.0e-3))
+        TrixiParticles.update_rigid_contact_eachstep!(rigid_system_2, v_ode_rigid,
+                                                      u_ode_rigid, semi_rigid, 0.0,
+                                                      (; dt=1.0e-3))
+        @test isempty(rigid_system_1.cache.contact_tangential_displacement)
+        @test isempty(rigid_system_2.cache.contact_tangential_displacement)
     end
 end