hoho

components.h

// data structures:
struct Particle {
    vec3 x;        // position
    vec3 v;        // velocity
    vec3 f;        // force accumulator
    vec3 f_ext;    // external forces (gravity, wind, etc.)
    float mass;
    int idx;       // index into global system (idx, idx+1, idx+2 for x,y,z)
};

struct Spring {
    Particle* p_i;  // endpoint i
    Particle* p_j;  // endpoint j
    float k;        // stiffness
    float L0;       // rest length
};

explicit:

// step:
for (s in springs)
    f = compute_force(s.p_i, s.p_j)
    s.p_i.f += f, s.p_j.f += -f
for (p in particles)
    p.f += p.f_ext                  // collect external forces like gravity
    p.x += dt * p.v_old             // use the old v from the previous step
    p.v_new = p.v_old + dt * p.mass.inv * p.f
    p.v_old = p.v_new

// Since we discretized time into steps (with given dt), we have the concept of
// "old" (previous step) and "new" (current step) values. We use subscripts or
// variable names to distinguish them.

NOTE

Since we discretized time into steps (with given dt), we have the concept of "old" (previous step) and "new" (current step) values. We use subscripts or variable names to distinguish them.

components.h

export default {
  data () {
    return {
      p.v_new = p.v_old + dt * p.mass.inv * p.f 
      M * v = M * v_old + dt * (f_spring + f_ext) 
    }
  }
}

explicit (symplectic):

// step:
for (s in springs)
    f = compute_force(s.p_i, s.p_j)
s.p_i.f += f, s.p_j.f += -f
for (p in particles)
    p.f += p.f_ext
p.v_new += p.v_old + dt * p.mass.inv * p.f // update velocity first
    p.x += dt * p.v_new                // we use v_new (not v_old) to update x!
    p.v_old = p.v_new

symplectic with system matrix:

// step:
n_p = particles.length
n_s = springs.length

M = zeros(n_p * 3, n_p * 3) // our particle has 3-dof in 3d space
b, v, f = zeros(n_p * 3) // col vecs matching the size of the system dof

for (s in springs)
    f = compute_force(s.p_i, s.p_j)
    s.p_i.f += f, s.p_j.f += -f
for ([p, idx] in particles)
    p.f += p.f_ext
    b[idx:idx+2] = p.mass * p.v + dt * p.f // collect rhs b, momentum
    M[idx:idx+2, idx:idx+2] = diag(p.mass)

// we want to find the new velocities (v) that satisfy the equation:
//
// M * v = M * v_old + dt * (f_spring + f_ext)                             eqn.1
//
// lhs: M * v               (`v` is unknown - what we solve for)
// rhs: M * v_old + dt * f  (all known - from current state)
// 
// new momentum (M * v) = old momentum (M * v_old) + impulse (dt * f) 
//
// Note: momentum (m*v) and impulse (dt*f) are different physical quantities,
// but they have the same units (kg * m/s), so they can be added together.
// Impulse represents the *change* in momentum.
// this is the discrete form of Newton's 2nd law btw.

// `v_old` we already know, `v` we do not know this one, so we need to solve the
// equations above.

// if we further decorate eqn.1 to the canonical form (Ax=b), we get:
//
// M * v = b                                                               eqn.2
//
// A * x = b
// |   |   |
// M   v   (M * v_old + dt * f)
//
// where A is the system matrix, x is the unknown, and b is the right-hand side.

// Fortunately, our system matrix M is diagonal, so solving `M \ b` is just some
// element-wise division `v[i] = b[i] / m[i]` - an O(n) operation, Thank god!
// Later on, we'll encounter much nastier matrices - each with there own ominous
// nickname - that will trap us down god know what rabbit holes..

v = M \ b
for ([p, idx] in particles)
    p.v += v[idx:idx+2]
    p.x += dt * p.v