Simulating rigid body collisions in 3d I have been reading about physics engines and I am confused on how one approaches simulating collision responses.
I read about the coefficient of restitution:
https://en.wikipedia.org/wiki/Coefficient_of_restitution
But it wasn't clear how this idea could be extended to 3d with rigid bodies. Conceptually you could make contacts non-hard and integrate a force + damping term over a number of iterations dependent upon penetration depth but that is not how engines with hard contact constraints work which is what I am interested in.
 A: $\newcommand{\b}  {\mathbf}$
Assume two particle collied in 3D
the equations are
\begin{align*}
 &m_1\,(\mathbf v_1-\mathbf u_1)=-\lambda\,\mathbf n\tag 1
 \end{align*}
\begin{align*}
  &m_2\,(\mathbf v_2-\mathbf u_2)=\lambda\,\mathbf n\tag 2
  \end{align*}
\begin{align*}
  &\left[(\mathbf v_2-\mathbf{v}_1)+\epsilon\,(\mathbf u_2-\mathbf u_1)\right]\cdot\mathbf n=0\tag 3
\end{align*}
you have 7 scalar equations for the 7 unknowns;  the  6 components of the vectors $~\mathbf v_i~$ and $~\lambda$
Adding equation (1) and (2) you obtain the conservation of the linear momentum and for $~\epsilon=1~$ the conservation of the energy .
where

*

*$\b v_1~,\b v_2~$ velocity after the collision

*$\b u_1~,\b u_2~$ velocity bevor the collision

*$m_i~$ particle masses

*$\b n~$ collision direction vector $\quad,\b n\cdot\b n=1$

*$\epsilon~$ coefficient of restitution
$\quad,\epsilon=0~$ perfectly inelastic collision
$\quad,\epsilon=1~$ perfectly elastic collision


Theory
starting with Newton equation immediately after the collision
\begin{align*}
 &m_i\,\frac{d\mathbf v'_i}{dt}= \pm\, F_c\,\mathbf n\quad\Rightarrow\\
 &m_i\,\int_{\mathbf u_i}^{\mathbf v_i}\,d\b v'_i=\pm\int  F_c\,\mathbf n\,dt=-\lambda\mathbf{n}
 \end{align*}
\begin{align*}
 &m_i\,(\mathbf v_i-\mathbf u_i)=\pm\lambda\,\mathbf n\quad i=1,2
\end{align*}
$~i=1~$ minus ,   $~i=2~$ plus
where $~ F_c~$ is the constraint force
Conservation of the energy
\begin{align*}
  &E=\frac{1}{2}\left(m_1\,(\mathbf{v}_1)^2+m_2\,(\mathbf{v}_2)^2-
  m_1\,(\mathbf{u}_1)^2-m_2\,(\mathbf{u}_2)^2\right)=0\\
  &2\,E=\left(m_1\,\left [(\mathbf{v}_1)^2- (\mathbf{u}_1)^2\right]
  +m_2\,\left[(\mathbf{v}_2)^2-
 (\mathbf{u}_2)^2\right]\right)=0\\
  &2\,E=\left(m_1\,\left [\mathbf{v}_1- \mathbf{u}_1\right]\cdot
  \left [\mathbf{v}_1+ \mathbf{u}_1\right]
  +m_2\,\left[\mathbf{v}_2-\mathbf{u}_2\right]
  \cdot \left[\mathbf{v}_2+\mathbf{u}_2\right]\right)=0\\
  &\text{with}\quad \mathbf{v}_1- \mathbf{u}_1=-\frac{\lambda}{m_1}\,\mathbf n
  \quad, \mathbf{v}_2- \mathbf{u}_2=\frac{\lambda}{m_2}\,\mathbf n\\
  &2\,E=\left(m_1\,\left [-\frac{\lambda}{m_1}\,\mathbf n\right]\cdot
  \left [\mathbf{v}_1+ \mathbf{u}_1\right]
  +m_2\,\left[\frac{\lambda}{m_2}\,\mathbf n\right]
  \cdot \left[\mathbf{v}_2+\mathbf{u}_2\right]\right)=0\quad\Rightarrow\\
  &2E=\left[(\mathbf v_2-\mathbf{v}_1)+(\mathbf u_2-\mathbf u_1)\right]\cdot\mathbf n=0
 \end{align*}
and  with the coefficient of restitution $~\epsilon~$
\begin{align*}
  &\left[(\mathbf v_2-\mathbf{v}_1)+\epsilon\,(\mathbf u_2-\mathbf u_1)\right]\cdot\mathbf n=0
 \end{align*}
thus for $~\epsilon=1~$ you obtain the conservation of the energy

Example
assume one dimensional
$$\b u_1=[u,0,0]^T~,\b u_2=[0,0,0]^T~,\b n=[1,0,0]^T$$
you obtain
$$\b v_1=\left[{\frac {u \left( m_{{1}}+m_{{2}}\epsilon  \right) }{m_{{2}}+m_{{1}}}}~,0~,0]^t\right]$$
$$\b v_2=\left[-{\frac {m_{{1}}u \left( -1+\epsilon  \right) }{m_{{2}}+m_{{1}}}}~,0~,0\right]^T$$
$$\lambda=-{\frac {m_{{2}}m_{{1}}u \left( -1+\epsilon  \right) }{m_{{2}}+m_{{1}}
}}
$$
$$2\,E={\frac {m_{{1}}{u}^{2}m_{{2}} \left( -1+{\epsilon }^{2} \right) }{m_{{
2}}+m_{{1}}}}
$$
A: In all three scenarios of two particles, a particle and a 3D body and two 3D bodies the calculation of the impulse magnitude $J$ is the same, once the reduced mass $m^\star$ of the contact is found.
$$ \boxed{ J = (1+\epsilon)\, m^\star \; v_{\rm imp}} $$
where $v_{\rm imp}$ is the relative speed of approach of the two bodies at the point of contact, $\epsilon$ is the coefficient of restitution and $m^\star$ is the reduced mass of the system.
For the three cases above, this is how to calculate the reduced mass given the contact normal direction $\boldsymbol{n}$.

*

*Two Particles of mass $m_1$ and $m_2$
$$ m^\star = \dfrac{1}{\frac{1}{m_1} + \frac{1}{m_2}} $$


*One Body and One Particle with masses $m_1$ and $m_2$ respectively, and mass moment of inertia tensor ${\bf I}_1$ for the body, as well as the position of the center of mass, relative to the contact point, denoted with the vector $\boldsymbol{d}_1$
$$ m^\star = \dfrac{1}{\frac{1}{m_1} + ( \boldsymbol{n}\times \boldsymbol{d}_1) \cdot {\bf I}_1^{-1} (\boldsymbol{n} \times \boldsymbol{d}_1) + \frac{1}{m_2}} $$


*Two Bodies with masses $m_1$ and $m_2$ respectively, and mass moment of inertia tensor ${\bf I}_1$ and ${\rm I}_2$, as well as the positions of the center of mass, relative to the contact point, denoted with the vectors $\boldsymbol{d}_1$ and $\boldsymbol{d}_2$
$$ m^\star = \dfrac{1}{\frac{1}{m_1} + ( \boldsymbol{n}\times \boldsymbol{d}_1) \cdot {\bf I}_1^{-1} (\boldsymbol{n} \times \boldsymbol{d}_1) + \frac{1}{m_2} + ( \boldsymbol{n}\times \boldsymbol{d}_2) \cdot {\bf I}_2^{-1} (\boldsymbol{n} \times \boldsymbol{d}_2)} $$
The above is the same you will find in the Collision Response Wikipedia article, as equation (5).
Note that $\times$ is the vector cross product, and $\cdot$ is the vector dot product.
It is also identical to equation (8-18) in the Physically Based Modeling, Lecture Notes II for rigid body simulations by Andrew Witkin and David Baraff. SIGCOURSE link

The exact calculation for $v_{\rm imp}$ is
$$ v_{\rm imp} =  \boldsymbol{n} \cdot ( \left(\boldsymbol{v}_1 + \boldsymbol{d}_1 \times \boldsymbol{\omega}_1\right) - \left( \boldsymbol{v}_2 + \boldsymbol{d}_2 \times \boldsymbol{\omega}_2 \right) ) $$
The above-calculated impulse $J$ obeys the law of contact
$$ v_{\rm bounce} = -\epsilon \; v_{\rm imp} $$
Here is an example implementation of the above calculation for two bodies in C#
public class Contact
{
    public Vector3 Position { get; set; }
    public Vector3 Normal { get; set; }
    public double Epsilon { get; set; }

    public double GetImpulse(RigidBody2 target, RigidBody2 contact, out double impactSpeed)
    {
        double m_1 = target.MassProperties.Mass;
        double m_2 = contact.MassProperties.Mass;
        Matrix3 I_1_inv = LinearAlgebra.Inverse(target.MassProperties.MMoi);
        Matrix3 I_2_inv = LinearAlgebra.Inverse(contact.MassProperties.MMoi);
        Vector3 d_1 = target.MassProperties.CG - Position;
        Vector3 d_2 = contact.MassProperties.CG - Position;
    
        Vector3 h_1 = LinearAlgebra.Cross(Normal, d_1);
        Vector3 h_2 = LinearAlgebra.Cross(Normal, d_2);
    
        double m_reduced = 1 / (
              1 / m_1 + LinearAlgebra.Dot(h_1, I_1_inv*h_1)
            + 1 / m_2 + LinearAlgebra.Dot(h_2, I_2_inv*h_2));
    
        impactSpeed = LinearAlgebra.Dot(Normal, 
             (contact.Velocity + LinearAlgebra.Cross(d_1, contact.Omega)) 
           - (target.Velocity) + LinearAlgebra.Cross(d_2, target.Omega));
    
        return (1 + Epsilon) * m_reduced * impactSpeed;
    }
}

