PBD Notes

PBA

Engine

  • number of constraint iterations
  • number of substeps
  • Type: Pbd or Xpbd

Substep

  • External Forces override
  • Update particle velocity:
  • Compute predicted position
  • Generate Collision Constraints override
  • Reset Lagrange multipliers (XPBD-only)
  • Solve constraints
  • Compute new velocity
  • Update position
  • Update Velocities (Damping) override
  • Clear some instant constraints (like collision etc.)
  • Advance current time

Particle

  • : current frame position
  • : velocity
  • : predicted position
  • : forces
  • : mass
  • : 1/mass

Constraint

Abstract Class

Type: Bilateral or Unilateral

  • particles

  • stiffness (PBD)

  • compliance (XPBD)

  • (XPBD)

  • Lagrange multiplier

  • Calculate the constraint function value C(x)

  • Calculate the derivative of the constraint function grad C(x)

  • Project associated particles to the constraint manifold

Project Particle

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
// constraint function value
const double c = calculateConstraintValue();

// if it is a unilateral constraint and is satisfied
if(type == Unilateral && c >= 0.0)
    return;
    
// derivative of the constraint function
// N: number of particles involved in this constraint
Matrix<N * 3, 1> grad_c = calculateGrad();
if(grad_c.norm() < 1e-12)
    return;

//----PBD----//
// Calculate s
// inv_M : Matrix<N * 3, 1>
// inv_M[i * 3 + 0/1/2] = particles[i].w
const double s = - c / (grad_c.transpose() * inv_M.asDiagonal() * grad_c);
Matrix<N * 3, 1> delta_x = stiffness * s * inv_M.asDiagonal() * grad_c;

// Update predicted positions
for(unsigned int i = 0; i < N; ++i)
    particles[i].p += delta_x[i * 3]; // delta_x[i*3] is a 3D vector
    
//----XPBD----//
// TODO

Different Type of Constraints

Some Derivation:

where is a matrix with property

When comes to the gradient of a normalized cross product: Assume

Substituting back:

Distance Constraint

  • , particles involved
  • stiffness/compliance
  • d: rest length

C(x)

1
return (p_0->p - p1->p).norm() - d;

Gradient

Derivation: 

1
2
3
4
5
6
7
8
9
10
const Vector3d& x_0 = p_0->p;
const Vector3d& x_1 = p_1->p;
const Vector3d dir = x_0 - x_1; // from p_1 to p_0
const double dist = dir.norm();

// choose a random direction when two particles are overlaped (degenerated case)
const Vector3d n = (dist < 1e-24) ? RandomNormalizedVector3() : (1.0 / dist) * dir;

grad_c[0] = n; // p_0 should move along n
grad_c[1] = -n; // p_1 should move against n

Bending Constraint

  • : vertices on a pair of adjacent triangles
  • stiffness/compliance
  • dihedral angle

C(x)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
const Vector3 p10 = p1 - p0;
const Vector3 p20 = p2 - p0;
const Vector3 p30 = p3 - p0;

const Vector3 n_0 = p10.cross(p20).normalized();
const Vector3 n_1 = p10.cross(p30).normalized();

const double crt_dihedral_angle = acos(clamp(n_0.dot(n_1), -1.0, 1.0));

// special case handling
assert(n_0.norm > 0.0);
assert(n_1.norm > 0.0);
assert(!isnan(crt_dihedral_angle));

return crt_dihedral_angle - dihedral_angle;

Gradient