Definitely not the way I would simulate a car crash but I would try to decrease the streching stiffness after detecting such a crash in order to make the car more soft for a couple of seconds, then would restore the stiffness value to bring back the car to the undeformed configuration.
It seems that the streching stiffness in your code is 0.005, but it could be a relaxation factor, so I don't know. Anyway your position-based solver loop would look something like this:
void b3Cloth::Step(float32 h, u32 iterations)
{
if (h == 0.0f)
{
return;
}
// Damping
float32 d = exp(-h * m_kd);
// Integrate using semi-implicit Euler
for (u32 i = 0; i < m_pCount; ++i)
{
b3Particle* p = m_ps + i;
p->v += h * p->im * m_gravity;
p->v *= d;
p->p0 = p->p;
p->p += h * p->v;
}
// Solve position constraints
for (u32 i = 0; i < iterations; ++i)
{
SolveC1();
}
// Estimate current velocity
float32 inv_h = 1.0f / h;
for (u32 i = 0; i < m_pCount; ++i)
{
b3Particle* p = m_ps + i;
p->v = inv_h * (p->p - p->p0);
}
// Solve velocity constraints (e.g. friction, etc)...
}
void b3Cloth::SolveC1()
{
for (u32 i = 0; i < m_c1Count; ++i)
{
b3C1* c = m_c1s + i;
b3Particle* p1 = m_ps + c->i1;
b3Particle* p2 = m_ps + c->i2;
float32 m1 = p1->im;
float32 m2 = p2->im;
float32 mass = m1 + m2;
if (mass == 0.0f)
{
continue;
}
mass = 1.0f / mass;
b3Vec3 J2 = p2->p - p1->p;
float32 L = b3Length(J2);
if (L > B3_EPSILON)
{
J2 /= L;
}
b3Vec3 J1 = -J2;
// m_k1 is the streching stiffness in [0, 1]
float32 C = L - c->L;
float32 impulse = -m_k1 * mass * C;
p1->p += (m1 * impulse) * J1;
p2->p += (m2 * impulse) * J2;
}
}
Hope it helps.