Mathematics behind Rigidbody Simulation

• Sphere Collision:
  1. Sphere-Plane Intersection: First, we make a vector from a point on the plane to the center of the sphere. Then, we get the dot product of this vector with the normal of the plane. If the result is less than or equal to the radius of the sphere, it means the sphere is colliding with the plane (either by just touching it or by getting already on the negative side of the plane)
  2. Sphere-Sphere Intersection: given spheres \( s_1, s_2 \) with radius \( R \) at \( o_1 \) and \( r \) at \( o_2 \) respectively. We calculate the squared distance between centers.\( vec_s = p_1 - p_2 \) and radius sum \( rs = R+r \). Sphere intersect if squared distance between centers is less than squared sum of radius: $$vec_s \cdot vec_s < rs * rs$$
  
• Box-Plane and Box-Box: Separating Axis Algorithm
  1. AABB-Plane Collision: To test if an AABB and plane intersect, we first have to project each vertex of the AABB onto the plane's normal. This leaves us with all the vertices of the AABB on a line. We then check the vertex that is furthest from the plane. If the vertex diagonally opposite that one is on the other side of the plane, we have an intersection. Here is the process:

     
     given plane p and AABB box b
     
     c =compute center of AABB box;
     e = get the positive extend ;
     r = compute the projection interval radius of b onto L(t) = b.c + t * p.n
     	r = e.x * p.n.x + e.y * p.n.y + e.z * p.n.z
     s = compute distance of box center from plane
     Intersection occurs when distance s falls within [-r, r]
     
     										

     
  2. Box-Box Collision: Separating Axis Theorem(SAT) project all of the points of both objects onto a line normal to the plane and observe that they form two separate non-overlapping intervals.
     If we limit ourselves to convex polyhedra, we can identify a finite number of possible axes to test. Specifically, For a box testing with arbitrary rotation, only 3 unique planes per object and only 3 unique edges. If each box has a 3x3 orientation matrix \(M_n\) whose column vectors are \(M_n = [a_n, b_n, c_n]\) then we would want to test all of the flowing axes: $$\begin{matrix} 𝐚1 & 𝐛1 & 𝐜1\\ 𝐚2 & 𝐛2 & 𝐜2\\ 𝐚1 \times 𝐚2 & 𝐚1 \times 𝐛2 & 𝐚1 \times 𝐜2 \\ 𝐛1 \times 𝐚2 & 𝐛1 \times 𝐛2 & 𝐛1 \times 𝐜2 \\ 𝐜1 \times 𝐚2 & 𝐜1 \times 𝐛2 & 𝐜1 \times 𝐜2 \end{matrix}$$
  
• Convex Polygon-Plane: Gilbert-Johnson-Keerthi (GJK) and Expanding Polytope Algorithm(EPA)
  1. The Gilbert-Johnson-Keerthi (GJK) algorithm is a popular method for computing the distance between two arbitrary convex objects. One nice feature is that it can be adapted to handle intersections between different object classes (such as spheres vs. convex polyhedral) with minimal effort. The GJK algorithm relies heavily on a concept called the Minkowski Sum. The Minkowski Sum conceptually is very easy to understand. Let’s say you have two shapes, the Minkowski Sum of those shapes is all the points in shape1 added to all the points in shape2: $$ A+B = \{a+b | a \in A, b \in B \}$$ The difference operation in the Minkowski sum is $$A-B = \{a-b | a \in A, b \in B \}$$ If two shapes are overlapping the Minkowski Difference will contain the origin as the image indicates:
     
     Therefore, we want to build a polygon that attempts to enclose the origin, if it's enclosed, we found the intersection. The polygon is called Simplex.
     To build the simplex, we need to find the support point, which is the farthest point(so that contains a maximum area) in a direction that is inside the Minkowski difference polygon. The k-simplex is the convex hull of k+1 affinely independent points in a k-dimensional space, so we get We do the iteration of different searching directions check of collision(by determing whether the simplex contains the origin). For example, here is the skeleton code for 2d GJK iterative alogirhtm

     Vector d = // choose a search direction
     // get the first Minkowski Difference point
     Simplex.add(support(A, B, d));
     // negate d for the next point
     d.negate();
     // start looping
     while (true) {
       // add a new point to the simplex because we haven't terminated yet
       Simplex.add(support(A, B, d));
       // make sure that the last point we added actually passed the origin
       if (Simplex.getLast().dot(d) <= 0) {
         // if the point added last was not past the origin in the direction of d
         // then the Minkowski Sum cannot possibly contain the origin since
         // the last point added is on the edge of the Minkowski Difference
         return false;
       } else {
         // otherwise we need to determine if the origin is in
         // the current simplex
         if (Simplex.contains(ORIGIN)) {
           // if it does then we know there is a collision
           return true;
         } else {
           // otherwise we cannot be certain so find the edge who is
           // closest to the origin and use its normal (in the direction
           // of the origin) as the new d and continue the loop
           d = getDirection(Simplex);
         }
       }
     }
     
     										

     for occasions of 3D, things are more complicated, because there are extra cases.
  2. The Expanding Polytope Algorithm(EPA) Determining Penetration Depth
     After GJK, we get a simplex that contains the origin and EPA helps to get the MTV(minimum translation vector) vector along which we can translate an intersecting shape to separate it from the other one, specificlly, the length of MTV is the penetration depth. Since the simplex contains the origin, and the point on the boundary of the Simplex that is closest to the origin is MTV. EPA finds the point by expanding inside the simplex-successively subdividing its edges with new vertices.
     The process is as following:
     1. First, we find the edge of the simplex closest to the origin
     2. Compute the support point in the Minkowski difference in a direction that is normal to the edge (i.e. perpendicular to it and pointing outside the polygon).
     3. Then we split this edge by adding this support point to it.
     4. We repeat these steps until the length and direction of the vector doesn’t change much. Once the algorithm converges, we have the MTV and the penetration depth.
     
  

Impulse-based approach operates on the velocity of the bodies and not on the force or acceleration. This means the solver calculates and applies a direct change in the linear and angular velocities of the bodies, instead of calculating and applying corrective forces and relying on integration to then change the velocities. With impulse-based dynamics, the goal is to find the impulses that result in velocities that solve the constraints. The general sequence of a simulation step using impulse-based dynamics is somewhat different from that of force-based engines:

• Compute all external forces.
• Apply the forces and determine the resulting velocities
• Calculate the constraint velocities based on the behavior functions.
• Apply the constraint velocities and simulate the resulting motion.

First let's define the impulse of force \(J_F\) as $$J_F = \int_{t_0} ^{t_1} F dt$$ During a small time slice \( \langle t_0, t_1 \rangle \). The impulse of force corresponds to the difference of linear moments $$\Delta p = p(t_1) - p(t_0) = J_F$$ while the impulsive torque of a force \(F\) applied in point \(r\) in the world space is defined as $$J_{\tau} = (r - r_c) \times J_F$$ Like the impulse of force changes the linear momentum, the impulsive torque changes the angular momentum $$\Delta b = b(t_1) - b(t_0) = J_{tau}$$ Next comes to how can we get impulse from collision process. As shown in the following figure, Object A and B collides at point P with relative position \(r^{AP}\) and \(r^{BP}\) and relative velocity as \(v ^{AB}\). The normal vector for the collision can be defined as $$v^{AB} \cdot n = (v ^{AP} - v^{BP}) \cdot n $$ We can express the impulse with a single scalar \(j\) times the normal, giving us \( jn \). According to Newton's Third law, the opposite force is applied to the other entity. The collision response equation for object can be derived by $$v_2^{A} = v_1^{A} + \frac{j}{M^{A}} n$$ $$v_2^{B} = v_1^{B} - \frac{j}{M^{B}} n$$ where M is the center of mass, so we solve the scalar j as $$j = \frac{-(1+e)v_1^{AB} \cdot n}{n \cdot n (\frac{1}{M^{A}} + \frac{1}{M^B})}$$