We implemented a rigid boby simulator which can model the interactions of non-deforming bodies in physically accurate ways. We used the clothsim project to initially create our gui and code structure. We've used the Forward Euler approach to simulate motion and the Discrete Element Method (Mishra 2003) to model collision forces. We also implemented a spatial map to reduce the time necessary to find collisions between particles. By the end, we could render complex scenes with many rigid bodies colliding with each other.
We used some python scripts to create scene files for complex scenes quickly and accurately. These files gave the locations of particles and rigid bodies, the initial velocity and angular velocity of rigid bodies, and the positions and normals of planes.
struct Particle {
// static values
double mass;
RigidBody* parent;
// dynamic values
Vector3D position;
Vector3D forces;
}
struct RigidBody {
// Rigid Body components
vector particles;
// Constants
double mass;
Matrix3x3 inertia_tensor;
// State variables
Vector3D position; // at the center of mass
Quaternion rotation;
Vector3D velocity;
Vector3D angular_velocity;
// Computed quantities
Vector3D force;
Vector3D torque;
};
The struct for particles kept track of their position for rendering purposes and the forces enacted on it, but its velocity was derived from the rigid body velocity, angular velocity, and the particle's location with respect to rigid body's center of mass. While initializing the rigid body, we also computed the mass, center of mass, and the inertia tensor. Let $S$ be the set of all particles. $$m = \sum_{i \in S} m_i$$ $$\vec{p} = \frac{1}{m}\sum_{in \in S} m_i \vec{p}_i$$ $$I = \sum_{i \in S} m_i \begin{bmatrix} r_{i, y}^2 + r_{i, z}^2 & -r_{i, x} r_{i, y} & -r_{i, x} r_{i, z} \\ -r_{i, y} r_{i, x} & r_{i, x}^2 + r_{i, z}^2 & -r_{i, y} r_{i, z} \\ -r_{i, z} r_{i, x} & -r_{i, z} r_{i, y} & r_{i, x}^2 + r_{i, y}^2 \end{bmatrix}$$ Where \(m_i, \vec{p}_i, r_i\) represent the mass, position, and the displacement between the position and the center of mass of particle \(i\), and \(m, \vec{p}, I\) represent the mass, position (of the center of mass) and the inertia tensor of the rigid body. These equations were given from David Baraff's note.
To perform rigid body to rigid body collision, we first take the location of each particle in the system and hash them, very similarly to how they are hashed in clothsim. Using these hashes, we build a spatial map of all the particles in our system, with each particle belonging to some 3D box volume, partitioning the space that all particles occupy. To find the particles that could possible collide with a given particle, we can compute the hash o \(\vec{p}_i\) and then check only those particles in the same bin as the particle in question. This greatly reduces the search space of potentially colliding particles, making this detection step have linear run-time (vs. quadratic run-time comparing all particles with all other particles) on average.
When two particles were within \(2 R\) of each other, where \(R\) is the radius of a particle, we said that these particles were in a collision. Note, we made each particle of the same radius, a choice which made this step slightly less complex and reduced the space needed to represent a particle. Upon a collision of particles \(i, j\), we applied the Discrete Element Method from Mishra's paper, and enacted three forces on particle \(i\): spring, dashpot, and shear. The spring force repels two particles that are too close to each other, the dashpot reduces the velocity of the colliding particles, and the shear force acts as kinetic friction.
Spring: \(\vec{F}_s = k (2R - |\vec{p}_i - \vec{p}_j|) \frac{\vec{p}_i - \vec{p}_j}{|\vec{p}_i - \vec{p}_j|} \)
Dashpot: \(\vec{F}_d = \eta (\vec{v}_j - \vec{v}_i)\)
Shear: \(\vec{F}_t = k_t (\vec{v}_j - \vec{v}_i - \text{proj}_{\vec{p}_j - \vec{p}_i}(\vec{v}_j - \vec{v}_i)) \)
We looped over every particle of every rigid body and added the necessary forces if a given particle was in collision with another particle.
After body-to-body collision, we then enacted body-to-plane collision using very similar physics. However, in order to compute the relative positions \( \vec{p}_i - \vec{p}_j \) and relative velocities \( \vec{v}_j - \vec{v}_i \), we had to compute the projection of the particle velocity onto the normal of the plane. Instead of \( 2R \) in the spring force, we also just used \( R \). These forces were also enacted on the particles in question.
After we find the force applied to each particle, we allow the rigid bodies to move for one time step. We must first compute the total force and torque on the system.
$$\vec{F} = \sum_{i \in S} F_i$$
$$\vec{\tau} = \sum_{i \in S} (\vec{p}_i - \vec{p}) \times \vec{F}_i$$
After we find the total force and torque on the system, we update the velocity and angular velocity using the Forward Euler Method.
$$\vec{v} + \frac{1}{m} \vec{F} dt \rightarrow \vec{v}$$
$$\vec{\omega} + R I^{-1} R^\intercal \vec{\tau} dt \rightarrow \vec{\omega}$$
Where \( R \) is the rotation matrix which can be derived from the rotation quaternion using library functions.
Finally, we update the position of each particle by translating and rotating it. The rotation is calculated by using euler angles and representing the rotation \( \vec{\omega} dt \) with a quaternion. We can then rotate an arbitrary displacement from a center of mass by a quaternion using a library function (call this \( f \)):
$$\vec{p}_i + \vec{v} dt + f(\vec{p}_i - \vec{p}, \vec{\omega}dt) \rightarrow \vec{p}_i$$
Finally, we update the rigid body position and rotation using
$$\vec{p} + \vec{v} dt \rightarrow \vec{p}$$
$$r_{\vec{\omega} dt} * r \rightarrow r$$
Where \( r \) is the quaternion representing rotation and \( r_{\vec{\omega} dt} \) is the quaternion which represents a rotation by \(\vec{\omega} dt\). Many of these equations, especially for rotation, came from David Baraff's note and from GPU Gems 3.
We rendered each particle as a sphere and the plane as simply a plane, and used the Blinn-Phong shading model and code from clothsim in order to render them.
One of our earliest problems was in the rendering of the plane. One can see in our milestone presentation that the plane, while physically existing in the right location and interacting with rigid bodies, was drawn through a particle. To resolve this, we actually found which particle the plane rendered at, and then just made a "ghost particle" of negligible size and mass at the location we wanted the plane to render. We also made sure to not simulate this particle or let it interact with the other rigid bodies in the system
We also found that the equation for the inertia tensor written in Baraff gave us a non-invertible matrix when the particles were co-linear. This was because when the particles were co-linear, spin along the axis of co-linearity caused the spheres themselves to spin about their centers, which required the moment of inertia of a sphere and lots of edge cases. Instead, we added a regularization term to the tensor in order to make the tensor invertible at the cost of loss of spin physicality for simpler rigid bodies.
One of the biggest things we learned throughout this project was the challenge of representing and modeling physical objects using a simplified discrete-time system. There was a lot of debate about what each struct needed to support itself, what was completely unnecessary, and the order in which we should perform operations. It was also challenging setting up our Cmake files to work for our new system (though less challenging than if we had started without any clothsim guidance).
We also learned about the balance between physical accuracy and speed of implementation/performance. The Forward Euler Method is very fast and easy, but it leads to numerical instability, which we observed when our collision constants had certain values. We also didn't have any forces that mimicked static friction, since the force acting against sliding friction was strictly velocity dependent.
We also learned that simulation is very difficult to run quickly on a CPU and that it would've been greatly benefited our ability to render more complex scenes within a reasonable time frame if we had enough time to work on a GPU version of the simulation.
Ivan completed the codebase setup, got the first particle and rigid body to form on-screen by modifying classes, GLShaders, and how the code turns JSON file inputs into rigid bodies. Implemented first iteration of collisions with static bodies (i.e. planes), though this was later replaced by a more physically accurate collision scheme. Created all input scenes (JSON files), the Python scripts used to generate said scenes, and all demos and videos for both testing and final presentation.
Jaymo Kang ported the initial code from Clothsim, set up the structure of the project, implemented collision physics and unconstrained movement, reviewed code, and created most of the slide decks and the final write-up.
Jessie Yang primarily worked on collision detection, a GUI for tuning collision parameters, and various codebase refactors to adjust rigid body APIs and properties. The bulk of her work was implementing collision detection and improving its performance; she also contributed some camera changes, a few JSON scenes, changes to JSON parsing, and color randomization of rigid bodies. She and Jaymo also worked on debugging collision physics. She also edited the milestone and final videos.
All team members collaborated to review pull requests, write and review write-ups, and compile slides.