Assignment 2: Rigid Bodies

**Posted:** Monday, 27 August 2018

**Updated:** Tuesday, 28 August 2018

**Due:** Monday, 10 September 2018

In this assignment you will implement a basic rigid body simulator, with support for simple collisions and/or constraints.

**Note:** Here is a new update of the OpenGL base code that fixes some bugs and adds some helper functions for drawing rigid bodies. Pay special attention to the comments at the top section of `draw.hpp`

.

Implement a simulation of a single rigid body, using the principles discussed in class. Typically, the rigid body representation is separated into two components:

The dynamics model, which stores the body’s physical properties and time-varying kinematic state. This is what we have talked about most of the time in class.

The shape representation, which describes the fixed geometrical shape of the body in body space. This is used for computing the inertia matrix, for rendering, and for collision detection. (In practice the rendered shape can be much more detailed than the collision shape, but we will use the same shape for both in this assignment.)

I have provided a partially implemented `Shape`

class (also included in the base code above) that takes care of the latter functionality. You should fill in the missing pieces, and then implement a rigid body object that tracks the dynamical state, and contains a `Shape`

as a member variable.

The rigid body should store its mass \(m\) and body-space inertia matrix \(\mathbf I_0\), its current position and rotation \(\mathbf q = (\mathbf x,\boldsymbol q)\), and its current linear and angular velocity \(\mathbf u = (\mathbf v,\boldsymbol\omega)\) (or angular momentum \(\mathbf L\) instead of \(\boldsymbol\omega\)). You may want to also precompute the inverse of the body-space inertia matrix to simplify the computation at run time.

As discussed in class, rotations are best represented using quaternions. If you are using C++, the Eigen linear algebra library also provides a `Quaternion`

class. You can also implement quaternions yourself using the description in the Witkin and Baraff notes.

Collision handling is hard to do in a generic way independent of the time integrator, so we will fix symplectic Euler as our time integration scheme. Implement the velocity update and position update of the rigid body as two separate functions, so that a symplectic Euler step consists of calling the velocity update followed by the position update (and later, we’ll be able to add collision handling by inserting it in between the two). Remember to normalize the rotation at the end of the position update.

**Notes:** (i) You may find that a spinning object gradually increases in energy due to the explicit integration of the gyroscopic force \(-\boldsymbol\omega\times\mathbf I\boldsymbol\omega\) in the velocity update. Try using an explicit midpoint scheme for this term instead. (ii) If you use Eigen’s quaternion class, it won’t let you do \(\boldsymbol q \mathrel{{+}{=}} \frac12\boldsymbol\omega\boldsymbol q\,\Delta t\) in the position update. A workaround is to update its `coeffs()`

instead.

To simplify the implementation of later extensions, you may also want to implement helper functions to transform a given point from body space to world space and vice versa, and to compute the world-space velocity of a given point.

Test your implementation on the following examples:

Create a unit-mass box of size \(0.4\times0.8\times0.1\), and give it an impulse of strength \(\mathbf j = (0,0,-1)\) at its right edge, i.e. \((0.2,0,0)\) relative to the center of mass. Verify that the resulting linear and angular velocity are consistent with the analytical solution. Initialize it rotated \(90^\circ\) about the \(z\)-axis and do the same test; verify that the result is identical to an unrotated box of size \(0.8\times0.4\times0.1\).

Demonstrate the tennis racket theorem on the \(0.4\times0.8\times0.1\) box: Give it an initial angular velocity slightly perturbed from a coordinate axis, for example \(\boldsymbol\omega = (0.1, 5.0, 0.1)\). Try this for all three coordinate axes. The rotation should be stable about when the box spins about the principal axes with maximum and minimum inertia, but it should tumble in an interesting way when spun about the intermediate axis. Draw arrows proportional to the angular velocity and angular momentum, to show that angular momentum remains nearly constant (up to time integration error).

Use zero gravity for both these test cases. In your submission, include videos for the three axes in the tennis racket theorem example.

Completing the base requirement is worth 60% of the points for this assignment. Each of the following four extensions is worth an additional 20%, so you need to do at least two of them to have the possibility of getting 100% points. Note that the total points for the assignment are capped at 100%. The extensions are mostly independent and do not have to be done in sequential order, i.e. you can choose any two or more to implement in your submission, not necessarily the first two. Your implemented extensions also do not have to be mutually compatible.

This component does not have any dependencies other than the base.

Support multiple rigid bodies in your program. Give each one a restitution coefficient \(\epsilon\) and friction coefficient \(\mu\). At each time step, perform collision detection between all pairs of bodies. If any pair is colliding, apply a collision impulse between them to resolve collisions.

To do this, it is useful to create a `World`

object that contains all of the rigid bodies, and takes care of time integration and collision handling. This way, you should not need to call the rigid body update functions directly from your main program.

To simplify collision detection, you only need to implement collision detection between sphere/sphere and sphere/ground pairs in this component. Discrete collision detection is fine, i.e. just check if the spheres are currently intersecting. Then you should set the collision point to be somewhere in their intersection, and the collision normal to be along the line joining the sphere centers. Do something similar for collisions between a sphere and a ground plane.

As suggested by Guendelman et al., perform collision response before the velocity and position updates. Do this via the following loop: Pick a collision pair whose relative normal velocity is non-separating. For the collision parameters, combine the coefficients of the two bodies as \(\epsilon=\min(\epsilon_1,\epsilon_2)\) and \(\mu=\max(\mu_1,\mu_2)\). Apply a collision impulse as described in the lecture slides, or in the Witkin and Baraff notes. Repeat until all collision pairs are separating, or until a maximum iteration count is reached.

Make sure your collision impulses are correct by checking that they result in the expected post-impulse relative normal velocity.

Demonstrate your implementation by making an animation of a billiards-type scenario, where a collection of balls on a table gets struck by another fast-moving ball. Friction with the table should cause the balls should roll realistically as well as translating.

This component depends on the first one.

Generalize your collision handling method to boxes. This will primarily require changing the collision detection method to test for box/box, box/sphere, and box/floor intersections.

For box/floor collisions, it is sufficient to only test the vertices of the box.

A simple, though inexact, way to test for intersections between general shapes is to test only a selected collection of sample points on one shape for interpenetration into the other (and vice versa if necessary). The provided `Shape`

code creates a set of sample points for an arbitrary box.

The sample points are defined in body space. To test a sample point for interpenetration into another body, you will have to (i) transform it into world space, then into the body space of the other body, (ii) test it for collision, and if so (iii) transform the collision normal back into world space using the rotation quaternion.

For box/sphere collisions, it is sufficient to take the *center* of the sphere as the sample point, test it against the box, and count it as a collision if its distance is less than the *radius* of the sphere. You will have to modify the collision point to be on the sphere surface instead of the center, though, to get an appropriate frictional response.

For box/box collisions, use the sample points of one box to test against the other, and vice versa as well.

Implementing reliable collision detection is a notoriously error-prone task, so you will probably need to debug your code a lot. It can be useful to visualize the collision pairs at each time step, by drawing a little sphere at the collision point and a little arrow coming out of it pointing along the collision normal. You will want to draw your spheres and boxes in wireframe so that you can see the actual collision points, and pause the animation whenever you see anything odd.

Demonstrate your more general collision detection implementation by creating an animation of multiple balls and boxes colliding in interesting ways, such as a bowling ball hitting a collection of tall boxes representing bowling pins.

This component does not depend on any other extensions.

Allow two or more rigid bodies to be attached together to create an articulated body. To do this, for each joint, specify a point in the body space of body 1 and a point in the body space of body 2, and constrain them to have the same world-space position. By adding constraints between multiple rigid bodies, you can connect them into chains, trees, and networks.

To enforce the constraints, you may use any of the techniques discussed in class, such as penalty forces or projection methods. In the latter case, see the lecture slides for the constraint Jacobian, and verify that your implementation is correct by using finite differences. A small amount of drift in the constraints is acceptable, but try not to let it get too large. If you’re using a penalty method, you may have to increase the strength of the penalty forces and use multiple small time steps per frame.

To demonstrate your implementation, show a multiple pendulum made of at least three rigid bodies, with the topmost one also connected to a fixed point.

Optionally, you could also create an articulated human-like character out of boxes, like a “ragdoll” used for unconscious characters in games. If you enable collision handling, disable it for bodies connected to each other.

This component depends on the first one, and can combine well with the second.

You may find that the pairwise collision resolution scheme is not sufficient to prevent bodies in stacks from sinking into each other. To fix this problem, implement shock propagation as described by Guendelman et al. This requires (i) building the contact graph and perfoming topological sorting, and (ii) iteratively applying zero-restitution impulses on layers of the contact graph from the bottom up. Run shock propagation as an extra step in the `World`

object’s symplectic time integrator, in between the velocity and position updates.

If constructing the contact graph is too much work, you could try simply sorting by the \(y\)-coordinate of the center of mass for a minor penalty in points.

When computing contact impulses for a body in the shock propagation pass, remember to only consider its collisions with other bodies below it (which would have already been processed by shock propagation). Those bodies should be treated as immovable, i.e. the contact impulse should be computed ignoring their mass and inertia matrix, and their linear and angular velocities should not be affected by it.

Demonstrate a stack of at least 10 objects that remains stable for several seconds. If you did not implement collision detection with boxes, simply create a stack of spheres, and reset their \(x\)- and \(z\)-coordinates to zero at the end of the time step to keep them lined up, as if they are in a tube of tennis balls. After a few seconds, hit it halfway up with another sphere (releasing the \(x\)- and \(z\)-coordinate constraints if any) to prove that the objects are not simply glued together.

Another useful test is to create a stack of blocks with overhangs given by the harmonic series (see block-stacking problem) and show that it is stable, while a stack with constant overhangs for each block is not. However, this is not a required part of the assignment.

A submission form will be opened on Moodle for this assignment. There, you should submit a zip file containing the following items:

a report in PDF format describing which extensions you implemented in this assignment, and a brief description of your submitted video(s),

a video or set of videos showing the results of your implementation, and

a copy of your implementation code.

Any additional information included in the report about how your implementation works will also be useful for my grading, and will help me give you more detailed feedback on your assignment.