Switch from using XPBD to classic PBD with substepping (#979)

Docs/Architecture.md

@@ -738,7 +738,7 @@ These jobs will do broadphase checks for all of the soft bodies. A thread picks
 
 ## Soft Body Simulate {#soft-body-simulate}
 
-These jobs will do the actual simulation of the soft bodies. They first collide batches of soft body vertices with the rigid bodies found during the Collide job (multiple threads can work on a single soft body) and then perform the simulation using XPBD (also partially distributing a single soft body on multiple threads).
+These jobs will do the actual simulation of the soft bodies. They first collide batches of soft body vertices with the rigid bodies found during the Collide job (multiple threads can work on a single soft body) and then perform the simulation using PBD (also partially distributing a single soft body on multiple threads).
 
 ## Soft Body Finalize {#soft-body-finalize}
 

Jolt/Physics/SoftBody/SoftBodyMotionProperties.cpp

@@ -246,12 +246,10 @@ void SoftBodyMotionProperties::IntegratePositions(const SoftBodyUpdateContext &i
 		}
 }
 
-void SoftBodyMotionProperties::ApplyVolumeConstraints(const SoftBodyUpdateContext &inContext)
+void SoftBodyMotionProperties::ApplyVolumeConstraints()
 {
 	JPH_PROFILE_FUNCTION();
 
-	float inv_dt_sq = 1.0f / Square(inContext.mSubStepDeltaTime);
-
 	// Satisfy volume constraints
 	for (const Volume &v : mSettings->mVolumeConstraints)
 	{
@@ -284,7 +282,7 @@ void SoftBodyMotionProperties::ApplyVolumeConstraints(const SoftBodyUpdateContex
 		JPH_ASSERT(w1 > 0.0f || w2 > 0.0f || w3 > 0.0f || w4 > 0.0f);
 
 		// Apply correction
-		float lambda = -c / (w1 * d1c.LengthSq() + w2 * d2c.LengthSq() + w3 * d3c.LengthSq() + w4 * d4c.LengthSq() + v.mCompliance * inv_dt_sq);
+		float lambda = -v.mStiffness * c / (w1 * d1c.LengthSq() + w2 * d2c.LengthSq() + w3 * d3c.LengthSq() + w4 * d4c.LengthSq());
 		v1.mPosition += lambda * w1 * d1c;
 		v2.mPosition += lambda * w2 * d2c;
 		v3.mPosition += lambda * w3 * d3c;
@@ -346,12 +344,10 @@ void SoftBodyMotionProperties::ApplySkinConstraints([[maybe_unused]] const SoftB
 	}
 }
 
-void SoftBodyMotionProperties::ApplyEdgeConstraints(const SoftBodyUpdateContext &inContext, uint inStartIndex, uint inEndIndex)
+void SoftBodyMotionProperties::ApplyEdgeConstraints(uint inStartIndex, uint inEndIndex)
 {
 	JPH_PROFILE_FUNCTION();
 
-	float inv_dt_sq = 1.0f / Square(inContext.mSubStepDeltaTime);
-
 	// Satisfy edge constraints
 	const Array<Edge> &edge_constraints = mSettings->mEdgeConstraints;
 	for (uint i = inStartIndex; i < inEndIndex; ++i)
@@ -367,7 +363,7 @@ void SoftBodyMotionProperties::ApplyEdgeConstraints(const SoftBodyUpdateContext
 		if (length > 0.0f)
 		{
 			// Apply correction
-			Vec3 correction = delta * (length - e.mRestLength) / (length * (v0.mInvMass + v1.mInvMass + e.mCompliance * inv_dt_sq));
+			Vec3 correction = e.mStiffness * delta * (length - e.mRestLength) / (length * (v0.mInvMass + v1.mInvMass));
 			v0.mPosition += v0.mInvMass * correction;
 			v1.mPosition -= v1.mInvMass * correction;
 		}
@@ -407,7 +403,7 @@ void SoftBodyMotionProperties::ApplyCollisionConstraintsAndUpdateVelocities(cons
 			// Remember previous velocity for restitution calculations
 			Vec3 prev_v = v.mVelocity;
 
-			// XPBD velocity update
+			// Verlet velocity update
 			v.mVelocity = (v.mPosition - v.mPreviousPosition) / dt;
 
 			// Satisfy collision constraint
@@ -634,7 +630,7 @@ void SoftBodyMotionProperties::StartNextIteration(const SoftBodyUpdateContext &i
 
 	IntegratePositions(ioContext);
 
-	ApplyVolumeConstraints(ioContext);
+	ApplyVolumeConstraints();
 }
 
 SoftBodyMotionProperties::EStatus SoftBodyMotionProperties::ParallelDetermineCollisionPlanes(SoftBodyUpdateContext &ioContext)
@@ -690,7 +686,7 @@ SoftBodyMotionProperties::EStatus SoftBodyMotionProperties::ParallelApplyEdgeCon
 					uint num_edges_to_process = non_parallel_group? edge_group_size : min(SoftBodyUpdateContext::cEdgeConstraintBatch, edge_group_size - edge_start_idx);
 					if (edge_group > 0)
 						edge_start_idx += mSettings->mEdgeGroupEndIndices[edge_group - 1];
-					ApplyEdgeConstraints(ioContext, edge_start_idx, edge_start_idx + num_edges_to_process);
+					ApplyEdgeConstraints(edge_start_idx, edge_start_idx + num_edges_to_process);
 
 					// Test if we're at the end of this group
 					uint edge_constraints_processed = ioContext.mNumEdgeConstraintsProcessed.fetch_add(num_edges_to_process, memory_order_relaxed) + num_edges_to_process;

Jolt/Physics/SoftBody/SoftBodyMotionProperties.h

@@ -27,8 +27,7 @@ class DebugRenderer;
 
 /// This class contains the runtime information of a soft body.
 //
-// Based on: XPBD, Extended Position Based Dynamics, Matthias Muller, Ten Minute Physics
-// See: https://matthias-research.github.io/pages/tenMinutePhysics/09-xpbd.pdf
+// Using a classic PBD simulation with sub stepping instead of running multiple Gauss-Seidel iterations before stepping.
 class JPH_EXPORT SoftBodyMotionProperties : public MotionProperties
 {
 public:
@@ -194,18 +193,18 @@ private:
 	void								IntegratePositions(const SoftBodyUpdateContext &inContext);
 
 	/// Enforce all volume constraints
-	void								ApplyVolumeConstraints(const SoftBodyUpdateContext &inContext);
+	void								ApplyVolumeConstraints();
 
 	/// Enforce all skin constraints
 	void								ApplySkinConstraints(const SoftBodyUpdateContext &inContext);
 
 	/// Enforce all edge constraints
-	void								ApplyEdgeConstraints(const SoftBodyUpdateContext &inContext, uint inStartIndex, uint inEndIndex);
+	void								ApplyEdgeConstraints(uint inStartIndex, uint inEndIndex);
 
 	/// Enforce all LRA constraints
 	void								ApplyLRAConstraints();
 
-	/// Enforce all collision constraints & update all velocities according the XPBD algorithm
+	/// Enforce all collision constraints & update all velocities
 	void								ApplyCollisionConstraintsAndUpdateVelocities(const SoftBodyUpdateContext &inContext);
 
 	/// Update the state of the soft body (position, velocity, bounds)

Jolt/Physics/SoftBody/SoftBodySharedSettings.cpp

@@ -32,14 +32,14 @@ JPH_IMPLEMENT_SERIALIZABLE_NON_VIRTUAL(SoftBodySharedSettings::Edge)
 {
 	JPH_ADD_ATTRIBUTE(SoftBodySharedSettings::Edge, mVertex)
 	JPH_ADD_ATTRIBUTE(SoftBodySharedSettings::Edge, mRestLength)
-	JPH_ADD_ATTRIBUTE(SoftBodySharedSettings::Edge, mCompliance)
+	JPH_ADD_ATTRIBUTE(SoftBodySharedSettings::Edge, mStiffness)
 }
 
 JPH_IMPLEMENT_SERIALIZABLE_NON_VIRTUAL(SoftBodySharedSettings::Volume)
 {
 	JPH_ADD_ATTRIBUTE(SoftBodySharedSettings::Volume, mVertex)
 	JPH_ADD_ATTRIBUTE(SoftBodySharedSettings::Volume, mSixRestVolume)
-	JPH_ADD_ATTRIBUTE(SoftBodySharedSettings::Volume, mCompliance)
+	JPH_ADD_ATTRIBUTE(SoftBodySharedSettings::Volume, mStiffness)
 }
 
 JPH_IMPLEMENT_SERIALIZABLE_NON_VIRTUAL(SoftBodySharedSettings::InvBind)

Jolt/Physics/SoftBody/SoftBodySharedSettings.h

@@ -101,11 +101,11 @@ public:
 
 		/// Constructor
 						Edge() = default;
-						Edge(uint32 inVertex1, uint32 inVertex2, float inCompliance = 0.0f) : mVertex { inVertex1, inVertex2 }, mCompliance(inCompliance) { }
+						Edge(uint32 inVertex1, uint32 inVertex2, float inStiffness = 0.9f) : mVertex { inVertex1, inVertex2 }, mStiffness(inStiffness) { }
 
 		uint32			mVertex[2];									///< Indices of the vertices that form the edge
 		float			mRestLength = 1.0f;							///< Rest length of the spring
-		float			mCompliance = 0.0f;							///< Inverse of the stiffness of the spring
+		float			mStiffness = 0.9f;							///< The stiffness of the spring. This value is normalized between 0 and 1. 1 means that the spring is completely stiff. Note that stiffness depends on the time step and the number of iterations, the smaller the time step or the higher the number of iterations, the stiffer the edge.
 	};
 
 	/// Volume constraint, keeps the volume of a tetrahedron constant
@@ -115,11 +115,11 @@ public:
 
 		/// Constructor
 						Volume() = default;
-						Volume(uint32 inVertex1, uint32 inVertex2, uint32 inVertex3, uint32 inVertex4, float inCompliance = 0.0f) : mVertex { inVertex1, inVertex2, inVertex3, inVertex4 }, mCompliance(inCompliance) { }
+						Volume(uint32 inVertex1, uint32 inVertex2, uint32 inVertex3, uint32 inVertex4, float inStiffness = 0.9f) : mVertex { inVertex1, inVertex2, inVertex3, inVertex4 }, mStiffness(inStiffness) { }
 
-		uint32			mVertex[4];									///< Indices of the vertices that form the tetrhedron
+		uint32			mVertex[4];									///< Indices of the vertices that form the tetrahedron
 		float			mSixRestVolume = 1.0f;						///< 6 times the rest volume of the tetrahedron
-		float			mCompliance = 0.0f;							///< Inverse of the stiffness of the constraint
+		float			mStiffness = 0.9f;							///< The stiffness of the constraint. This value is normalized between 0 and 1. 1 means that the volume is completely stiff. Note that stiffness depends on the time step and the number of iterations, the smaller the time step or the higher the number of iterations, the stiffer the volume.
 	};
 
 	/// An inverse bind matrix take a skinned vertex from its bind pose into joint local space

Samples/Tests/SoftBody/SoftBodyLRAConstraintTest.cpp

@@ -22,7 +22,7 @@ void SoftBodyLRAConstraintTest::Initialize()
 	auto inv_mass = [](uint, uint inZ) { return inZ == 0? 0.0f : 1.0f; };
 	Ref<SoftBodySharedSettings> cloth_settings = SoftBodyCreator::CreateCloth(cNumVerticesX, cNumVerticesZ, cVertexSpacing, inv_mass);
 	for (SoftBodySharedSettings::Edge &e : cloth_settings->mEdgeConstraints)
-		e.mCompliance = 1.0e-3f; // Soften the edges a bit so that the effect of the LRA constraints is more visible
+		e.mStiffness = 0.1f; // Soften the edges a bit so that the effect of the LRA constraints is more visible
 	SoftBodyCreationSettings cloth(cloth_settings, RVec3(-10.0f, 25.0f, 0), Quat::sIdentity(), Layers::MOVING);
 	mBodyInterface->CreateAndAddSoftBody(cloth, EActivation::Activate);
 

Samples/Tests/SoftBody/SoftBodySkinnedConstraintTest.cpp

@@ -65,7 +65,7 @@ void SoftBodySkinnedConstraintTest::Initialize()
 
 	// Make edges soft
 	for (SoftBodySharedSettings::Edge &e : settings->mEdgeConstraints)
-		e.mCompliance = 1.0e-3f;
+		e.mStiffness = 0.1f;
 
 	// Create inverse bind matrices by moving the bind pose to the center of mass space for the body
 	Array<Mat44> bind_pose = GetWorldSpacePose(0.0f);

Samples/Utils/SoftBodyCreator.cpp

@@ -41,7 +41,6 @@ Ref<SoftBodySharedSettings> CreateCloth(uint inGridSizeX, uint inGridSizeZ, floa
 		for (uint x = 0; x < inGridSizeX; ++x)
 		{
 			SoftBodySharedSettings::Edge e;
-			e.mCompliance = 0.00001f;
 			e.mVertex[0] = vertex_index(x, z);
 			if (x < inGridSizeX - 1)
 			{
@@ -279,7 +278,6 @@ Ref<SoftBodySharedSettings> CreateSphere(float inRadius, uint inNumTheta, uint i
 		for (uint theta = 0; theta < inNumTheta - 1; ++theta)
 		{
 			SoftBodySharedSettings::Edge e;
-			e.mCompliance = 0.0001f;
 			e.mVertex[0] = vertex_index(theta, phi);
 
 			e.mVertex[1] = vertex_index(theta + 1, phi);