Restitution 1.0 produces higher bounce every time.

[Wc4ever]

When I set everything to 0.0 but 1 rigid body to restitution 1.0 every next bounce goes higher. Is it intended behaviour? Usually, it is bouncing on the same height in other physics engines.

git patch with example 0001-restitution-bug.patch

[jrouwe]

Hello,

This is a known issue and has been discussed before, see: https://github.com/godot-jolt/godot-jolt/issues/374#issuecomment-1555939027

It is something I should look into though (at least confirm that other physics engines don't have this behavior / don't do the opposite of slowly losing energy).

[Wc4ever]

Hi, thanks for answer. I could post info of how it works in other engines if you wish.

Here is test from Bullet. Every damping and friction set to 0. Looks like bullet calculates restitution as multiplication so both objects has to be 1.0 otherwise it dont bounce at all.

https://github.com/user-attachments/assets/a77816e9-bf8e-4a81-b4f2-ea640982fd31

Open Dynamics Engine surface bounce 1.0

https://github.com/user-attachments/assets/7bfe22de-275a-4f6c-a946-d5a2a9b35db9

Unity build-in engine Gosh looks like it grows in Unity also(unexpected). Material settings Bounciness 1.0 Combine type Maximum.

https://github.com/user-attachments/assets/0c643329-8a24-4ede-a158-98fa86779385

[Wc4ever]

Flax Engine (Physx 5) Wierdest one so far. Set both material Combine type Maximum and restitution 1.0 and it grows every bounce almost the same as in Jolt. But if I set Combine type Multiply and all materials to 1.0 seems it slows down but eventualy grows up again and then slows down.

https://github.com/user-attachments/assets/e56fa9a3-c09d-48fd-93ef-17a5c9f8ff70

So seems I was wrong about "other engines" but at the same time rebouns at the same height makes more sense to me.

At this point if you belive there is no bug I think it could be just declarated as expected behaviour and closed.

[jrouwe]

I'll leave the bug open for now. It's definitively not physically correct. I'll see if there's a way to make this better.

Mainly for future me, I've simplified the simulation loop to:

// Code to simulate a sphere of radius 2 falling from 10 units high
Vec3 pos = Vec3(0, 10, 0);
Vec3 vel = Vec3::sZero();
Vec3 g(0, -9.81f, 0);
float dt = 1.0f / 60.0f;
float maxy = 0;
for (int i = 0; i < 1000; ++i)
{
	// Apply gravity
	// Equivalent to PhysicsSystem::JobApplyGravity
	vel += g * dt;

	// If we're penetrating the ground, we have to reverse the velocity (corresponds to restitution 1).
	// Equivalent to PhysicsSystem::JobSolveVelocityConstraints
	float penetration = 2.0f - pos.GetY();
	if (penetration > 0 && vel.GetY() < 0.0f)
		vel = -vel;

	// Update position
	// Equivalent to PhysicsSystem::JobIntegrateVelocity
	pos += vel * dt;

	// Record new height records (note we only trace when we're going down again to avoid tracing multiple lines for the same bounce)
	if (vel.GetY() <= 0.0f && pos.GetY() > maxy)
	{
		maxy = pos.GetY();
		Trace("Simulated new high: %f", maxy);
	}
}

which produces:

Simulated new high: 9.997275
Simulated new high: 10.209825
Simulated new high: 10.425099
Simulated new high: 10.643100
Simulated new high: 10.863826
Simulated new high: 11.087276
Simulated new high: 11.313452

When simulating the same setup in Jolt we get the exact same results:

New high point: 10.000000
New high point: 10.209825
New high point: 10.425100
New high point: 10.643101
New high point: 10.863826
New high point: 11.087276
New high point: 11.313451

Implementation: SimpleTest.zip

As said in the other ticket: When the collision has been detected, the gravity has been applied already so the bounce velocity is overestimated.

[wiltosoft]

https://github.com/user-attachments/assets/7b9c515f-9490-4a32-9d19-ebed71b9a34b

Could this be related to experiencing infinite bounce? This is with friction/restitution 0.8/0.8 on the ball and 1.0/0.9 on the plane.

[jrouwe]

If I adjust the code above with your parameters:

// Code to simulate a sphere of radius 2 falling from 10 units high
Vec3 pos = Vec3(0, 10, 0);
Vec3 vel = Vec3::sZero();
Vec3 g(0, -9.81f, 0);
float dt = 1.0f / 60.0f;
float restitution_floor = 0.9f;
float restitution_sphere = 0.8f;
float restitution_combined = max(restitution_floor, restitution_sphere);
for (int i = 0; i < 1000; ++i)
{
	bool positive_before = vel.GetY() >= 0.0f;

	// Apply gravity
	// Equivalent to PhysicsSystem::JobApplyGravity
	vel += g * dt;

	// If we're penetrating the ground, we apply restitution
	// Equivalent to PhysicsSystem::JobSolveVelocityConstraints
	float penetration = 2.0f - pos.GetY();
	if (penetration > 0 && vel.GetY() < 0.0f)
		vel = -restitution_combined * vel;

	// If velocity changes from positive to negative, we trace the height
	bool positive_after = vel.GetY() >= 0.0f;
	if (positive_before && !positive_after)
		Trace("Bounce height: %f", pos.GetY());

	// Update position
	// Equivalent to PhysicsSystem::JobIntegrateVelocity
	pos += vel * dt;
}

Then no:

Bounce height: 10.000000
Bounce height: 8.627144
Bounce height: 7.512838
Bounce height: 6.609490
Bounce height: 5.857011
Bounce height: 5.247907
Bounce height: 4.739192
Bounce height: 4.324055
Bounce height: 3.972433
Bounce height: 3.679149

But as you can see from the code above, it depends on more than just the restitution (gravity and delta time).

[wiltosoft]

Thank you for your reply @jrouwe :) Appreciated.

I took your code and changed the iterations to 4000 and I see the same effect, I think? I can provide the physics world code I am using for the above video too. I may well be doing something daft tho'!

Edit: FWIW, I think this pretty much sorts it for both situations. I'm not yet familiar enough with your code to suggest a patch tho'!

    // Code to simulate a sphere of radius 2 falling from 10 units high
    Vec3 pos = Vec3(0, 10, 0);
    Vec3 vel = Vec3::sZero();
    Vec3 g(0, -9.81f, 0);
    float dt = 1.0f / 60.0f;
    float restitution_floor = 0.9f;
    float restitution_sphere = 0.8f;
    float restitution_combined = max(restitution_floor, restitution_sphere);
    for (int i = 0; i < 10000; ++i)
    {
        bool positive_before = vel.GetY() >= 0.0f;

        // Apply gravity
        // Equivalent to PhysicsSystem::JobApplyGravity
        vel += g * dt;

        // If we're penetrating the ground, we have to reverse the velocity (corresponds to restitution 1).
        // Equivalent to PhysicsSystem::JobSolveVelocityConstraints
        float penetration = 2.0f - pos.GetY();
        if (penetration > 0 && vel.GetY() < 0.0f)
            vel = -restitution_combined * (vel - g * dt);


        // If velocity changes from positive to negative, we trace the height
        bool positive_after = vel.GetY() >= 0.0f;
        if (positive_before && !positive_after)
            Trace("Iteration %5d - Bounce height: %f", i, pos.GetY());

        // Update position
        // Equivalent to PhysicsSystem::JobIntegrateVelocity
        pos += vel * dt;
    }

[jrouwe]

You're right, I didn't simulate it for long enough. It is indeed the same thing. At some point the velocity that we gain in one time step due to gravity equals the velocity we lose during the bounce and we end up in a steady state. Your solution is similar to what I had in mind. And yes, the Jolt code is a bit more involved than this simple example, so it needs a bit more thought.

[mutexmiles]

I might make this into a PR if you think it's viable, but this is my attempt at (at least approximately) fixing this issue:

(I am user wiltosoft FYI!)

Index: Jolt/Physics/Body/MotionProperties.h
IDEA additional info:
Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP
<+>UTF-8
===================================================================
diff --git a/Jolt/Physics/Body/MotionProperties.h b/Jolt/Physics/Body/MotionProperties.h
--- a/Jolt/Physics/Body/MotionProperties.h	(revision 0a7d2501b33e1251101392d2be403c551ce5f963)
+++ b/Jolt/Physics/Body/MotionProperties.h	(date 1731011352112)
@@ -42,6 +42,9 @@
 	/// Get world space linear velocity of the center of mass
 	inline Vec3				GetLinearVelocity() const										{ JPH_ASSERT(BodyAccess::sCheckRights(BodyAccess::sVelocityAccess(), BodyAccess::EAccess::Read)); return mLinearVelocity; }
 
+    /// Get world space linear velocity of the center of mass
+    inline Vec3				GetLinearVelocityGravityDelta() const										{ JPH_ASSERT(BodyAccess::sCheckRights(BodyAccess::sVelocityAccess(), BodyAccess::EAccess::Read)); return mGravityVelocityDelta; }
+
 	/// Set world space linear velocity of the center of mass
 	void					SetLinearVelocity(Vec3Arg inLinearVelocity)						{ JPH_ASSERT(BodyAccess::sCheckRights(BodyAccess::sVelocityAccess(), BodyAccess::EAccess::ReadWrite)); JPH_ASSERT(inLinearVelocity.Length() <= mMaxLinearVelocity); mLinearVelocity = LockTranslation(inLinearVelocity); }
 
@@ -235,6 +238,7 @@
 	// 1st cache line
 	// 16 byte aligned
 	Vec3					mLinearVelocity { Vec3::sZero() };								///< World space linear velocity of the center of mass (m/s)
+	Vec3                    mGravityVelocityDelta  { Vec3::sZero() };
 	Vec3					mAngularVelocity { Vec3::sZero() };								///< World space angular velocity (rad/s)
 	Vec3					mInvInertiaDiagonal;											///< Diagonal of inverse inertia matrix: D
 	Quat					mInertiaRotation;												///< Rotation (R) that takes inverse inertia diagonal to local space: Ibody^-1 = R * D * R^-1
Index: Jolt/Physics/Body/MotionProperties.inl
IDEA additional info:
Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP
<+>UTF-8
===================================================================
diff --git a/Jolt/Physics/Body/MotionProperties.inl b/Jolt/Physics/Body/MotionProperties.inl
--- a/Jolt/Physics/Body/MotionProperties.inl	(revision 0a7d2501b33e1251101392d2be403c551ce5f963)
+++ b/Jolt/Physics/Body/MotionProperties.inl	(date 1731011352141)
@@ -123,8 +123,11 @@
 	JPH_ASSERT(mCachedBodyType == EBodyType::RigidBody);
 	JPH_ASSERT(mCachedMotionType == EMotionType::Dynamic);
 
+    // Gravity factor (stored for reference in collision checks)
+    mGravityVelocityDelta = inDeltaTime * (mGravityFactor * inGravity);
+
 	// Update linear velocity
-	mLinearVelocity = LockTranslation(mLinearVelocity + inDeltaTime * (mGravityFactor * inGravity + mInvMass * GetAccumulatedForce()));
+	mLinearVelocity = LockTranslation(mLinearVelocity + mGravityVelocityDelta + inDeltaTime * mInvMass * GetAccumulatedForce());
 
 	// Update angular velocity
 	mAngularVelocity += inDeltaTime * MultiplyWorldSpaceInverseInertiaByVector(inBodyRotation, GetAccumulatedTorque());
Index: Jolt/Physics/Constraints/ContactConstraintManager.cpp
IDEA additional info:
Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP
<+>UTF-8
===================================================================
diff --git a/Jolt/Physics/Constraints/ContactConstraintManager.cpp b/Jolt/Physics/Constraints/ContactConstraintManager.cpp
--- a/Jolt/Physics/Constraints/ContactConstraintManager.cpp	(revision 0a7d2501b33e1251101392d2be403c551ce5f963)
+++ b/Jolt/Physics/Constraints/ContactConstraintManager.cpp	(date 1731011352134)
@@ -59,16 +59,27 @@
 
 	// Calculate velocity of collision points
 	Vec3 relative_velocity;
-	if constexpr (Type1 != EMotionType::Static && Type2 != EMotionType::Static)
-		relative_velocity = inBody2.GetMotionPropertiesUnchecked()->GetPointVelocityCOM(r2) - inBody1.GetMotionPropertiesUnchecked()->GetPointVelocityCOM(r1);
-	else if constexpr (Type1 != EMotionType::Static)
-		relative_velocity = -inBody1.GetMotionPropertiesUnchecked()->GetPointVelocityCOM(r1);
-	else if constexpr (Type2 != EMotionType::Static)
-		relative_velocity = inBody2.GetMotionPropertiesUnchecked()->GetPointVelocityCOM(r2);
+    // Gravity delta is also retrieved. There are some issues here, such as possible differences
+    // in gravity factor between the two objects which ought to be taken into account, but this
+    // provides a cheap approximation that will be good enough in most situations
+    float normal_gravity_velocity;
+	if constexpr (Type1 != EMotionType::Static && Type2 != EMotionType::Static) {
+        relative_velocity = inBody2.GetMotionPropertiesUnchecked()->GetPointVelocityCOM(r2) - inBody1.GetMotionPropertiesUnchecked()->GetPointVelocityCOM(r1);
+        normal_gravity_velocity = (inBody2.GetMotionPropertiesUnchecked()->GetLinearVelocityGravityDelta() - inBody1.GetMotionPropertiesUnchecked()->GetLinearVelocityGravityDelta()).Dot(inWorldSpaceNormal);
+    }
+	else if constexpr (Type1 != EMotionType::Static) {
+        relative_velocity = -inBody1.GetMotionPropertiesUnchecked()->GetPointVelocityCOM(r1);
+        normal_gravity_velocity = inBody1.GetMotionPropertiesUnchecked()->GetLinearVelocityGravityDelta().Dot(inWorldSpaceNormal);
+    }
+	else if constexpr (Type2 != EMotionType::Static) {
+        relative_velocity = inBody2.GetMotionPropertiesUnchecked()->GetPointVelocityCOM(r2);
+        normal_gravity_velocity = inBody2.GetMotionPropertiesUnchecked()->GetLinearVelocityGravityDelta().Dot(inWorldSpaceNormal);
+    }
 	else
 	{
 		JPH_ASSERT(false); // Static vs static makes no sense
 		relative_velocity = Vec3::sZero();
+        normal_gravity_velocity = 0;
 	}
 	float normal_velocity = relative_velocity.Dot(inWorldSpaceNormal);
 
@@ -93,7 +104,7 @@
 		// position rather than from a position where it is touching the other object. This causes the object
 		// to appear to move faster for 1 frame (the opposite of time stealing).
 		if (normal_velocity < -speculative_contact_velocity_bias)
-			normal_velocity_bias = inSettings.mCombinedRestitution * normal_velocity;
+			normal_velocity_bias = inSettings.mCombinedRestitution * (normal_velocity - normal_gravity_velocity);
 		else
 			// In this case we have predicted that we don't hit the other object, but if we do (due to other constraints changing velocities)
 			// the speculative contact will prevent penetration but will not apply restitution leading to another artifact.

[jrouwe]

Thanks for looking into this! I think this is indeed more or less the fix. I don't think it is needed to actually store mGravityVelocityDelta though (I'm trying to keep the per body memory usage down).

There are some issues here, such as possible differences in gravity factor between the two objects which ought to be taken into account

What do you mean with this? I think gravity factor is taken into account here.

[mutexmiles]

So yes, it is taken into account but effectively averaged and applied to both objects equally in normal_velocity_bias? Perhaps I am mistaken, but I was thinking you'd need to apply a slightly different bias to each object for it to be accurate.

For the gravity delta, yes, I was unsure how to best access the information needed (I have only lightly nosied around your code so far!) to calculate it so this seemed simpler.

[mutexmiles]

Erm... @jrouwe so I updated my project and the bug is back :(

I think the problem is here:

Index: Jolt/Physics/Constraints/ContactConstraintManager.cpp
IDEA additional info:
Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP
<+>UTF-8
===================================================================
diff --git a/Jolt/Physics/Constraints/ContactConstraintManager.cpp b/Jolt/Physics/Constraints/ContactConstraintManager.cpp
--- a/Jolt/Physics/Constraints/ContactConstraintManager.cpp	(revision 88de579d67d75ae91d12964fd8e919317da5eca4)
+++ b/Jolt/Physics/Constraints/ContactConstraintManager.cpp	(date 1731175336373)
@@ -71,7 +71,7 @@
 		const MotionProperties *mp1 = inBody1.GetMotionPropertiesUnchecked();
 		const MotionProperties *mp2 = inBody2.GetMotionPropertiesUnchecked();
 		relative_velocity = mp2->GetPointVelocityCOM(r2) - mp1->GetPointVelocityCOM(r1);
-		gravity_dt_dot_normal = inGravityDeltaTimeDotNormal * (mp2->GetGravityFactor() - mp1->GetGravityFactor());
+		gravity_dt_dot_normal = inGravityDeltaTimeDotNormal * (mp2->GetGravityFactor() * mp1->GetGravityFactor());
 	}
 	else if constexpr (Type1 != EMotionType::Static)
 	{

[jrouwe]

I don't think that is correct.

In your original code you had (replacing inBodyX.GetMotionPropertiesUnchecked() with mpX):

normal_gravity_velocity = (mp2->GetLinearVelocityGravityDelta() - mp1->GetLinearVelocityGravityDelta()).Dot(inWorldSpaceNormal);

combining this with the implementation of GetLinearVelocityGravityDelta:

mGravityVelocityDelta = inDeltaTime * (mGravityFactor * inGravity);

we get:

normal_gravity_velocity = (inDeltaTime * (mp2->GetGravityFactor() * inGravity) - inDeltaTime * (mp1->GetGravityFactor() * inGravity)).Dot(inWorldSpaceNormal);

which equals:

normal_gravity_velocity = (inDeltaTime * inGravity).Dot(inWorldSpaceNormal) * (mp2->GetGravityFactor() - mp1->GetGravityFactor());

which equals my code:

gravity_dt_dot_normal = inGravityDeltaTimeDotNormal * (mp2->GetGravityFactor() - mp1->GetGravityFactor());

I think this also makes sense intuitively: When both bodies have a gravity factor of 1 they are both receiving the same gravity delta, which means that their relative velocity did not change during the time step, so for the restitution calculations we don't need to apply any correction.

mp2->GetGravityFactor() - mp1->GetGravityFactor() = 0, so no correction will be applied. mp2->GetGravityFactor() * mp1->GetGravityFactor() = 1, which would incorrectly apply a correction.

What is your exact repro case? If it is a body bouncing on the floor then you should not even be hitting this code path as this is for the dynamic vs dynamic case.

[mutexmiles]

Yes! You are correct :) I was hasty.

This all makes sense to me and you guessed my problem in your last paragraph! Apologies for wasting your time!