Soft Body Physics Simulation

Role: Physics Programmer
Engine: Unity, C#
Project Type: Solo, Personal Project
Duration: ~3 months

About this project

I started this project during my 2nd year at Staffordshire university, not too long after completing my module on Applied Mathematics. I really enjoyed physics programming, so I decided to work on this personal project in my spare time. Its a custom soft body physics simulation built in Unity using Verlet integration, iterative constraint solving, pressure-based volume preservation, collision resolution, and CCD inverse kinematics. I designed the project as a ground-up physics experiment without relying on Unity’s built-in rigidbody joints, which taught me to debug difficult Heisenbugs, crucial maths for creating stable, optimised and accurate physics.

Technical Overview

Custom Verlet Integration

Implemented position-based integration storing current and previous positions to derive implicit velocity and enable stable constraint solving.

Iterative Constraint Solver

Distance constraints solved using weighted inverse-mass distribution across multiple solver iterations per physics step.

Surface Collision Resolution

Custom signed-distance plane collision with tangent/normal velocity decomposition for stable bounce control.

Volume Preservation

Polygon area calculated using the Shoelace formula and corrected via pressure-based force distribution.

CCD Inverse Kinematics

Cyclic Coordinate Descent solver integrated directly with point-mass system.

Modular Physics Solver

Centralized solver orchestrating integration, constraints, collision, area correction, and IK per FixedUpdate step.

Verlet Integration Core

Implemented position-based Verlet integration for stable simulation without explicit velocity storage.

Integration Step

Velocity is implicitly derived from position delta. Damping applied to control energy growth. Continuous forces accumulated before integration. 


// Verlet integration
Vector2 implicitVelocity = position - prevPosition;
implicitVelocity *= 0.98f;

prevPosition = position;
position += implicitVelocity;
position += netForce * Time.deltaTime * Time.deltaTime;

This approach improves numerical stability when combined with iterative constraint correction.

Distance & Follower Constraints

Implemented weighted distance constraints with inverse-mass distribution.

Distance Constraint Resolution


Vector2 delta = c.pointB.position - c.pointA.position;
float distance = delta.magnitude;

float difference = (distance - c.restLength) / distance;
Vector2 offset = delta * difference * c.springCoefficient;

float totalWeight = c.pointA.inverseMass + c.pointB.inverseMass;

c.pointA.position += offset * (c.pointA.inverseMass / totalWeight);
c.pointB.position -= offset * (c.pointB.inverseMass / totalWeight);

Constraints are solved iteratively per frame to improve convergence. Inverse mass enables fixed or dynamic points.

Custom Surface Collision

Implemented signed distance collision against planar surfaces with velocity decomposition for stable response. 


float signedDistance = Vector2.Dot(toPoint, normal);

if (signedDistance <= point.radius)
{
    float penetration = point.radius - signedDistance;
    point.position += normal * penetration;

    Vector2 v = point.position - point.prevPosition;

    Vector2 vNormal = normal * Vector2.Dot(v, normal);
    Vector2 vTangent = v - vNormal;

    point.prevPosition = point.position - vTangent;
}

Separating normal and tangential components allows controlled bounce and prevents instability.

Pressure-Based Volume Constraint

Maintains soft-body volume using polygon area preservation.

Area Calculation (Shoelace Formula)


for (int i = 0; i < points.Count; i++)
{
    int next = (i + 1) % points.Count;
    area += (points[i].position.x * points[next].position.y 
           - points[next].position.x * points[i].position.y);
}
area *= 0.5f;

Pressure force derived from area error and distributed using partial derivatives of the polygon area.

CCD Inverse Kinematics

Implemented Cyclic Coordinate Descent solver integrated directly into the same point-mass simulation. 


float cross = (r2.x * r1.y) - (r2.y * r1.x);
float dot = Vector2.Dot(r2, r1);
float angle = Mathf.Atan2(cross, dot);

Uses cross and dot products to compute signed angle correction. Joints above the pivot are rotated incrementally for convergence.

Physics Solver Architecture

Centralized solver executes in FixedUpdate: 


void FixedUpdate()
{
    foreach (var point in _points)
        point.VerletStep();

    for (int i = 0; i < 7; i++)
    {
        ResolveCollisions();
        ConstraintStep();
        CalculateArea();
        CCDSolver();
    }
}

Multiple solver passes improve constraint stability and reduce jitter.

Engineering Challenges Solved