# Assignment 3

### Table of Contents

- 1 Installation
- 2 Code Structure
- 3 Numerical Integration
- 4 Collision Response
- 5 Extra Credit
- 6 Submission

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 A3 starter code from: https://github.coecis.cornell.edu/cs5643/assignments. You should be able to see 7 files in total:

`scene.py`

`collision.py`

`response.py`

`rigid_sim.py`

`GUI.py`

`util.py`

`pa3_main.py`

`pa3_main.py`

is the main file you need to run for this
assignment. To debug and visualize your code, run ```
python
pa3_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 `pa3_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 = A3GUI(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 `SimGUI`

s is fed to `A3GUI`

, 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 `A3GUI`

about the test it should
show first.

## 2.2 Main Loop

We move on to the rest of the code in `pa3_main.py`

to
understand that how the three classes `Scene`

,
`Collision`

, and `CollisionResponse`

interact to
update the rigid body states.

\[ \begin{align*}\texttt{while }& \texttt{ gui.is_window_running:}\\ &\color{magenta} \texttt{gui.currsim.scene.update_vel()}\\ &\color{blue} \texttt{gui.currsim.collision.clearCollision()}\\ &\color{blue} \texttt{gui.currsim.collision.collide_bounds()}\\ &\color{blue} \texttt{gui.currsim.collision.collide_all()}\\ &\color{green} \texttt{gui.currsim.collision.response.PGS()}\\ &\color{green} \texttt{gui.currsim.collision.response.apply_impulses()}\\ &\color{magenta}\texttt{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
`i`

th 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 `contact.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 `contact.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\).

- box-wall collision: a rotating box translating towards a wall
- box-box collision: a non-rotating box collides translating towards another box staying still
- box-corner collision: a box translating towards a corner formed by boundaries.

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 will 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 complementarity 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: constructing the linear complementarity problem, solving the problem, and applying the impulses to update the velocities of the rigid bodies.

The first step is covered in the functions
`CollisionResponse.init_contact`

and
`CollisionResponse.addContact`

. `init_contact`

is
called at the start of a simulation to initialize all the fields that
are required for the linear complementarity 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 should terminate when the impulse
vector \(\boldsymbol{\gamma}\)
converges. To determine when the solver has converged, you can check if
\(\|\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 can 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 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

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

will 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 `pa3_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 A3GUI instantiation code in
`pa3_main.py`

to ```
gui = A3GUI(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 A3GUI
instantiation code in `pa3_main.py`

to ```
gui =
A3GUI(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