diff --git a/docs/src/systems/rigid_body.md b/docs/src/systems/rigid_body.md
index 08e50e13c1..27928499ce 100644
--- a/docs/src/systems/rigid_body.md
+++ b/docs/src/systems/rigid_body.md
@@ -17,6 +17,37 @@ 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` currently defines a normal spring-dashpot contact law with the
+parameters `normal_stiffness`, `normal_damping`, and `contact_distance`.
+
+The current implementation uses the same model for rigid-wall and rigid-rigid contact:
+
+- rigid-wall contact groups penetrating wall neighbors into a small number of contact
+  manifolds per rigid particle and applies one normal contact force per manifold,
+- rigid-rigid contact evaluates direct pairwise normal contact forces between rigid
+  particles,
+- and both paths are currently normal-only, i.e. there are no tangential/frictional
+  forces or contact-history terms yet.
+
+Here, a contact manifold is a discrete approximation of one locally smooth contact
+patch. A rigid particle touching a flat wall will usually produce one manifold,
+while corners or edges can produce several.
+
+The number of cached rigid-wall manifolds per rigid particle is controlled by the
+`RigidBodySystem(...; max_manifolds=8)` keyword argument. If more wall-contact
+patches are detected than cached manifold slots are available, the implementation
+falls back to the best-matching existing manifold for that particle.
+
+`contact_distance` defines when contact starts. If `contact_distance == 0`, the
+particle spacing of the `RigidBodySystem` is used.
+
+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/src/general/abstract_system.jl b/src/general/abstract_system.jl
index fde3b92a7b..ce3feb424c 100644
--- a/src/general/abstract_system.jl
+++ b/src/general/abstract_system.jl
@@ -179,7 +179,11 @@ 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
+
+function reset_interaction_caches!(system)
     return system
 end
 
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index 87a9070702..6f8f59197e 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -529,6 +529,16 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t)
     end
 end
 
+# Some systems accumulate pairwise interaction state outside `dv_ode`. Reset that state once
+# at the beginning of every explicitly assembled interaction pass.
+function reset_interaction_caches!(semi::Union{NamedTuple, Semidiscretization})
+    foreach_system(semi) do system
+        reset_interaction_caches!(system)
+    end
+
+    return semi
+end
+
 # The `SplitIntegrationCallback` overwrites `semi_wrap` to use a different
 # semidiscretization for wrapping arrays.
 # TODO `semi` is not used yet, but will be used when the source terms API is modified
@@ -668,6 +678,8 @@ end
 end
 
 function system_interaction!(dv_ode, v_ode, u_ode, semi)
+    reset_interaction_caches!(semi)
+
     # Call `interact!` for each pair of systems
     foreach_system(semi) do system
         foreach_system(semi) do neighbor
@@ -700,6 +712,8 @@ end
 # Function barrier to make benchmarking interactions easier.
 # One can benchmark, e.g. the fluid-fluid interaction, with:
 # dv_ode, du_ode = copy(sol.u[end]).x; v_ode, u_ode = copy(sol.u[end]).x;
+# For manual multi-pair interaction assembly, call `reset_interaction_caches!(semi)` once
+# before the first direct `interact!` call.
 # @btime TrixiParticles.interact!($dv_ode, $v_ode, $u_ode, $fluid_system, $fluid_system, $semi);
 @inline function interact!(dv_ode, v_ode, u_ode, system, neighbor, semi; timer_str="")
     dv = wrap_v(dv_ode, system, semi)
diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl
index b8216db0e6..a36cd45aa5 100644
--- a/src/io/write_vtk.jl
+++ b/src/io/write_vtk.jl
@@ -433,6 +433,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 7a9b290965..108a075ae6 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_models.jl b/src/schemes/structure/rigid_body/contact_models.jl
index a75b50261e..0a91d50f90 100644
--- a/src/schemes/structure/rigid_body/contact_models.jl
+++ b/src/schemes/structure/rigid_body/contact_models.jl
@@ -5,10 +5,10 @@ abstract type AbstractRigidContactModel end
                       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.
+Shared rigid-contact model used by the active rigid-wall and rigid-rigid contact paths.
+The current 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.
 
 !!! warning "Experimental implementation"
     This is an experimental feature and may change in future releases.
diff --git a/src/schemes/structure/rigid_body/rhs.jl b/src/schemes/structure/rigid_body/rhs.jl
index 20c685f697..48fc8d26c8 100644
--- a/src/schemes/structure/rigid_body/rhs.jl
+++ b/src/schemes/structure/rigid_body/rhs.jl
@@ -108,6 +108,10 @@ function interact!(dv, v_particle_system, u_particle_system,
 
     NDIMS = ndims(particle_system)
     ELTYPE = eltype(particle_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
 
     # 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.
@@ -128,6 +132,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)
@@ -160,7 +166,6 @@ 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
@@ -174,13 +179,24 @@ function interact!(dv, v_particle_system, u_particle_system,
                                 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
 
@@ -331,6 +347,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;
@@ -354,10 +378,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
 
@@ -366,8 +388,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 0ba34b1fce..033ea9f1da 100644
--- a/src/schemes/structure/rigid_body/system.jl
+++ b/src/schemes/structure/rigid_body/system.jl
@@ -119,7 +119,8 @@ function RigidBodySystem(initial_condition; boundary_model=nothing,
         inverse_inertia = Ref(zero(SMatrix{3, 3, ELTYPE, 9}))
     end
 
-    cache = (;
+    cache = (; 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 +150,17 @@ function create_cache_contact_manifold(::Nothing, ::Val{NDIMS}, ELTYPE,
     return (;)
 end
 
-# Allocate per-particle manifold scratch arrays for rigid-wall contact.
+# 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),
@@ -374,6 +377,14 @@ function update_boundary_interpolation!(boundary_model, system::RigidBodySystem,
     return system
 end
 
+function reset_interaction_caches!(system::RigidBodySystem)
+    set_zero!(system.force_per_particle)
+    system.cache.contact_count[] = 0
+    system.cache.max_contact_penetration[] = zero(eltype(system))
+
+    return system
+end
+
 function update_final!(system::RigidBodySystem, v, u, v_ode, u_ode, semi, t)
     system_coords = current_coordinates(u, system)
     system_velocity = current_velocity(v, system)
@@ -387,10 +398,6 @@ 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)
-
     system.center_of_mass[] = center_of_mass
     system.center_of_mass_velocity[] = center_of_mass_velocity
     system.inertia[] = rotational_kinematics.inertia
@@ -555,6 +562,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,
@@ -564,6 +573,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
 
@@ -573,6 +583,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/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 13cb8affe7..a1ea2b6de9 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
 
@@ -783,6 +787,7 @@
         TrixiParticles.update_final!(rigid_system_2, v_rigid_2, u_rigid_2,
                                      v_ode_rigid, u_ode_rigid, semi_rigid, 0.0)
 
+        TrixiParticles.reset_interaction_caches!(semi_rigid)
         TrixiParticles.interact!(dv_ode_rigid, v_ode_rigid, u_ode_rigid,
                                  rigid_system_1, rigid_system_2, semi_rigid)
         force_after_forward_1 = copy(rigid_system_1.force_per_particle)
@@ -812,6 +817,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
@@ -841,6 +850,34 @@
         @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.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
+        @test vec(rigid_system_1.force_per_particle[:, 1]) ≈ collect(expected_force)
+        @test vec(rigid_system_2.force_per_particle[:, 1]) ≈ collect(-expected_force)
+
+        TrixiParticles.update_systems_and_nhs(v_ode_rigid, u_ode_rigid, semi_rigid, 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
+        @test vec(rigid_system_1.force_per_particle[:, 1]) ≈ collect(expected_force)
+        @test vec(rigid_system_2.force_per_particle[:, 1]) ≈ collect(-expected_force)
+
+        TrixiParticles.set_zero!(dv_ode_reset)
+        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
+        @test vec(rigid_system_1.force_per_particle[:, 1]) ≈ collect(expected_force)
+        @test vec(rigid_system_2.force_per_particle[:, 1]) ≈ collect(-expected_force)
+
         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,
@@ -856,6 +893,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]
@@ -944,6 +1004,12 @@
                                                    contact_model=contact_model,
                                                    max_manifolds=0)
 
+        system_meta_data = Dict{String, Any}()
+        TrixiParticles.add_system_data!(system_meta_data, rigid_system)
+        @test system_meta_data["contact_model"]["normal_stiffness"] ≈ 2.0e4
+        @test system_meta_data["contact_model"]["normal_damping"] ≈ 20.0
+        @test system_meta_data["contact_model"]["contact_distance"] ≈ 0.1
+
         semi = Semidiscretization(rigid_system, boundary_system)
         ode = semidiscretize(semi, (0.0, 0.01))
         v_ode, u_ode = ode.u0.x
@@ -973,6 +1039,7 @@
         @test kick_dv[2, 1] > 0
         @test kick_rigid_system.resultant_force[][2] > 0
 
+        TrixiParticles.reset_interaction_caches!(semi)
         TrixiParticles.interact!(dv_ode, v_ode, u_ode, rigid_system, boundary_system, semi)
         dv = TrixiParticles.wrap_v(dv_ode, rigid_system, semi)
         v_rigid = TrixiParticles.wrap_v(v_ode, rigid_system, semi)
@@ -981,6 +1048,35 @@
                                              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
+        direct_force = copy(rigid_system.force_per_particle)
+        direct_resultant_force = rigid_system.resultant_force[]
+
+        TrixiParticles.set_zero!(dv_ode)
+        TrixiParticles.update_final!(rigid_system, v_rigid, u_rigid, v_ode, u_ode, semi,
+                                     0.0)
+        TrixiParticles.reset_interaction_caches!(semi)
+        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.cache.contact_count[] == 1
+        @test rigid_system.cache.max_contact_penetration[] ≈ 0.05
+        @test rigid_system.force_per_particle == direct_force
+        @test rigid_system.resultant_force[] ≈ direct_resultant_force
+
+        TrixiParticles.set_zero!(dv_ode)
+        TrixiParticles.update_final!(rigid_system, v_rigid, u_rigid, v_ode, u_ode, semi,
+                                     0.0)
+        TrixiParticles.reset_interaction_caches!(semi)
+        TrixiParticles.finalize_interaction!(rigid_system, dv, v_rigid, u_rigid,
+                                             dv_ode, v_ode, u_ode, semi)
+
+        @test all(iszero, dv)
+        @test iszero(rigid_system.resultant_force[])
+        @test iszero(rigid_system.resultant_torque[])
+        @test iszero(rigid_system.angular_acceleration_force[])
 
         far_rigid_ic = InitialCondition(; coordinates=reshape([0.0, 0.09], 2, 1),
                                         velocity=rigid_velocity,
@@ -1002,6 +1098,7 @@
         short_support_v_ode, short_support_u_ode = short_support_ode.u0.x
         short_support_dv_ode = zero(short_support_v_ode)
 
+        TrixiParticles.reset_interaction_caches!(short_support_semi)
         TrixiParticles.interact!(short_support_dv_ode, short_support_v_ode,
                                  short_support_u_ode, far_rigid_system,
                                  short_support_boundary, short_support_semi)