Assignment 2

Table of Contents

In this assignment you will build a rigid body simulator for rectangular boxes in 2D. The simulation includes the motion of the boxes and their interactions via collision with the environment and with one another, and it includes a basic but functional contact solver that can successfully resolve many simultaneous collisions and stacks of objects up to several boxes high. You’ll implement time integration using Symplectic Euler (which is pretty simple) and collision response by solving for impulses using constraints on objects’ relative velocities (known as a velocity-level linear complimentarity problem or LCP).

The collision detection code is provided. In addition to running the contact detection algorithm we discussed in lecture, it also handles the slightly messy details of finding a suitable pair of contact points in the common case when rectangles collide along parallel edges. This makes it a lot easier to get nice behavior for stationary stacks of objects.

This assignment closely follows the lectures and the slides and notes that go with them, so you’ll want to refer to those materials for many of the details.

1 Installation

Pull the latest A2 starter code from: https://github.coecis.cornell.edu/cs5643/assignments-sp25. You should be able to see 7 files in total:

pa2_main.py is the main file you need to run for this assignment. To debug and visualize your code, run python pa2_main.py. For this assignment, we will stick to arch=ti.cpu to make sure that the starter code is runnable on most machines. A GUI with 5 buttons corresponding to five test cases of your code will show up after you run the command. By clicking on different buttons, you can switch between different tests, and test out your code in different contexts. Note that the first time you click on a particular test, the initialization of the new test scene will take some time. Please wait for a moment if you observe a significant frame rate drop after test switching.

2 Code Structure

The code roughly follows the MVC model, where scene.py stores the states of the rigid bodies, namely their positions, linear velocities and angular velocities, and collision.py stores all the detected collisions. The visualization code is in GUI.py.

The control code of our rigid body simulation has three main parts: collision detection, collision response and numerical integration, which are found in collision.py, response.py and scene.py respectively. The focus of this programming assignment is the collision response and numerical integration of a rigid body system. Thus, you only need to implement response.py and scene.py for this assignment, and collision.py contains a complete implementation of collision detection.

2.1 Construction of Simulation Tests

Now, let’s look at the main file pa2_main.py to understand how the three main classes Scene, Collision, and CollisionResponse are connected and initialized to create a rigid simulation instance.

inits = [Init.BW_C, Init.BC_C, Init.BB_C, Init.STACK, Init.RANDOM]
rigidSims = {init:SimGUI(RigidSim(init, dt, β, Cr, μ)) for init in inits}
gui = A2GUI(width, rigidSims, Init.BW_C)

The inits is a list of enum values that represent different initial states of our simulations. We iterate through this list to initialize instances of RigidSim, each of which represents a rigid body simulation test. Specifically, each RigidSim contains three components: RigiSim.scene, RigidSim.collision , and RigidSim.response, which store information related to a scene, the collisions, and the collision response respectively. Then, we instantiate a SimGUI class for each RigidSim so that we know how to draw each rigid body simulation. Eventually, this dictionary of SimGUIs is fed to A2GUI, which is the main GUI that controls the window size and how to switch between visualization of different tests. The Init.BW_C argument here is being used to inform A2GUI about the test it should show first.

2.2 Main Loop

We move on to the rest of the code in pa2_main.py to understand that how the three classes Scene, Collision, and CollisionResponse interact to update the rigid body states.

\[ \begin{align*} \mathrm{while }\quad&\mathrm{ gui.is\_window\_running:}\\ &\color{magenta} \mathrm{gui.currsim.scene.update\_vel()}\\ &\color{blue} \mathrm{gui.currsim.collision.clearCollision()}\\ &\color{blue} \mathrm{gui.currsim.collision.collide\_bounds()}\\ &\color{blue} \mathrm{gui.currsim.collision.collide\_all()}\\ &\color{green} \mathrm{gui.currsim.collision.response.PGS()}\\ &\color{green} \mathrm{gui.currsim.collision.response.apply\_impulses()}\\ &\color{magenta}\mathrm{gui.currsim.scene.update\_posn()}\\ &\cdots \end{align*} \]

The magenta part shows the numerical integration functions you will implement later. At the start of the loop, we will advance the linear velocity and angular velocity using the scene’s gravity value.

Then, we enter the blue block of code that is related to collision detection. The Collision.clearCollision function clears all the contact points detected in the previous iteration. Collision.collide_bounds detects all the collision points between boxes and fixed boundaries and Collision.collide_all detects all the collision points between boxes and boxes. In both Collision.collide_bounds and Collision.collide_all , the function CollisionResponse.addcontact will be called so that you can store the information related to these contact points, which defines the linear complementarity problem (LCP) you will solve later.

Thereafter, the green code CollisionResponse.PGS and CollisionResponse.apply_impulses are collision response codes that you will implement later. Assuming that you already constructed the linear complementarity system in the process of calling your CollisionResponse.addContact for each collision point, Your CollisionResponse.PGS implementation will be responsible for using the Projected Gauss-Seidel Algorithm (PGS) to solve the LCP. Then we call your CollisionResponse.apply_impulses implementation to update the velocities of the rigid boxes, using the impulses that were determined by PGS.

Finally, we use the updated velocities to update the position and orientation of the boxes, which you implement in Scene.update_posn.

The rest of the code in the main loop is for GUI update and we skip over those here.

2.3 Rigid Body Data

As for the details regarding how the rigid-body related data is being stored: Scene is the class where we initialize different tests and store the states of the rigid boxes. There are Scene.N boxes being stored in a Taichi field Scene.boxes, which holds objects of a Taichi dataclass called BoxState. For detailed information about the usage of Taichi dataclasses, reference the Taichi documentation. Briefly, to access the information related to the ith box, use Scene.boxes[i].xx, where xx is the attribute you would like to access. For details about the information stored in each box, see the BoxState data class in scene.py. In general, if you need to update information related to the rigid boxes, always go back to Scene.boxes. This is a field that stores most of the information you need.

2.4 Collision Data

Collision detection is run by calling Collision.collide_bounds and Collision.collide_all. These methods will determine which boxes are colliding, compute contact points and normals, and call the function ContactState.addContact in response.py for each contact. For detailed information on the arguments passed to addContact (which it is your job to write), please reference the function’s doc-string in response.py.

Also, detection of each contact point causes some information to be stored for visualization. Specifically, a boolean value for each rigid body that indicates whether it is involved in any collisions is stored in Collision.coll, and the positions of contact points are stored in the Taichi field Collision.collPs. Both of these two fields will be used for visualization later. The colliding objects will have white edges and collision points will be highlighted as magenta circles.

In your implementation process, feel free to borrow linear algebra code from util.py. Both 2D cross product and 2D transformation code all reside in this file.

3 Numerical Integration

If you run the starter code and try out the different tests, you should be able see to some boxes frozen in the air. They are not moving since the numerical integration code is left empty for now.

Thus, let’s start with implementing a numerical integrator that should be familiar to you now: Symplectic Euler.

Specifically, the Symplectic Euler can be divided into two parts: update_vel and update_posn in scene.py. update_vel is called before update_posn in the main loop, so you should update the velocities stored in Scene.boxes first using the gravity stored at Scene.g, and then later update the position in update_posn using the new velocities stored in Scene.boxes (which by then will have been modified to account for contacts). Note that you should update both the orientation and the position, based on the angular velocity and the linear velocity respectively.

To test your numerical integration, you can try out the tests below with \(C_r = 1\) and \(\beta = 0\).

You can check the reference videos of these three tests at this phase by clicking on the links attached to the test names and use them to debug your implementation. The box-wall collision test could be useful for debugging your angular velocity update, since the box in that test has nonzero initial angular velocity.

4 Collision Response

Now let’s move on to implement the collision response so that your boxes would bounce back when they hit the fixed boundaries or other boxes. To achieve this, you need to implement the CollisionResponse class stored in response.py .

In this assignment, we would like to iteratively solve the contact impulse by using the Projected Gauss-Seidel Solver (PGS). The PGS will be used to solve a linear complementary problem as follows:

\[ \begin{align*} \mathbf{v}_{\textrm{res}} &= JM^{-1}J^T \boldsymbol{\gamma} + ( 1 + C_r) JV + \beta \mathbf{d}\\ \boldsymbol{\gamma} &\geq \mathbf{0}\\ \mathbf{v}_{\textrm{res}} &\geq \mathbf{0}\\ \mathbf{v}_{\textrm{res}} \odot \boldsymbol{\gamma} &= \mathbf{0}\\ \end{align*} \]

where \(C_r\) is the restitution ratio and \(\beta\) is the coefficient for overlap-repair impulse.

To compute and use the impulses, there are three required steps: construction of the linear complementary problem, solving the problem and apply the impulses so that the velocities of the rigid bodies would be updated respectively.

The first step is covered in the functions CollisionResponse.init_contact and CollisionResponse.addContact. init_contact would be used at the start of a simulation to initialize all the fields that are required for the linear complementary problem.addContact will be called in the main update loop whenever a collision point is detected. You can use it to procedurally construct the J and M matrices in the problem.

Then, to solve the complementary problem, you can implement the Projected Gauss-Seidel algorithm in PGS to solve the problem. Your PGS iterative process could terminate when the impulse vector \(\boldsymbol{\gamma}\) converges. To determine if the impulse vector converges, you can check if the \(\|\Delta\boldsymbol{\gamma}\| < 10^{-4}\).

Eventually, to apply the impulses to the rigid bodies, update the velocities based using the solved impulses in the function apply_impulses .

If you have implemented these three functions correctly, you could run your implementation against the first three tests: box-wall collision, box-box collision, and box-corner collision with \(C_r = 1\) and \(\beta = 0\). You should be able to see that your boxes bouncing off walls and each other now. Also, we help you check if the kinetic energy of all boxes is conserved before and after applying impulses. If your implementation is not conserving kinetic energy, an error message saying that

the energy before applying impulse xxx and the energy after applying impulse xxx are significantly different.

would show up. Here are the three reference videos for the box-wall collision, box-box collision, and box-corner collision. Use those reference videos and the kinetic energy conservation error message to debug your implementation.

Eventually, you can test your nonoverlap-repair impulse implementation by switching the coefficients to \(C_r = 0.7\) and \(\beta = 0.8 /dt\) (These two coefficients are commented in pa2_main.py). Then compare your simulation results against our reference results for the box stacking test and the random box falling test.

5 Extra Credit

You can add friction to your CollisionResponse class and expand your linear complementary problem respectively. You can set \(\mu = 0.8\) and use the box stacking test to test your friction code. You should observe that the boxes are no longer sliding off each other over time.

6 Submission

Use the parameters \(C_r = 1\) , \(\beta = 0\) and \(\mu = 0\) to run the box-wall collision test, box-box collision test and box-corner collision test. Record your simulation by changing your A2GUI instantiation code in pa2_main.py to gui = A2GUI(width, rigidSims, xx , record=True), where xx = Init.BW_C, Init.BC_C, Init.BB_C .

Use the parameters \(C_r = 0.7\) , \(\beta = 0.8 / dt\) and \(\mu = 0\) to run the box stacking test, random box test. Record your simulation by changing your A2GUI instantiation code in pa2_main.py to gui = A2GUI(width, rigidSims, xx , record=True), where xx = Init.STACK, Init.RANDOM.

Combine all of your recordings into one demo video and submit it. Also, as what you did for previous assignments, submit a pdf file including the link to the commit of your submission