| CS 5643
Rigid Body Contact
Due date: Sun,
Apr 5, 2009 (before
In this assignment you
will implement a 2D rigid-body system with broad
and narrow phase
collision detection and velocity-level contact resolution based on a
Gauss-Seidel complementarity constraint solver. The starter code
provides you with a simple frictionless penalty-force simulator, and
allows you to draw simulation scenarios and input them into your
simulator as images.
You will test the scalability of your simulator on various image-based
simulation scenarios, such as towers of stacked rigid blocks.
on your own, or in a group of two people.
Starter Code (cs5643.rigidbody):
This project has significant starter code, primarily to support
rigid-body dynamics, all-pairs penalty-based contact, image-based scene
construction, and OpenGL
rendering. It is available via CMS.
In this assignment,
you will modify
this package as needed. Feel free to borrow code
from your previous
assignments, but all code used should be your own.
- RigidImageSimulation is
access point which accepts one optional argument specifying the image
simulate. Images must be
uncompressed 24-bit TGA format, and specify
non-white objects on a white background. Pin constraints immobilize rigid
bodies, and occur for image objects that are pure shades of
blue, (0,0,b), or lighter, (a,a,b). Test
simulation images are provided and described below. The default
simplest simulation is
"lcp.tga" to start. You'll also find all keyboard and
mouse controls in this
object; a summary of default key assignments is as follows:
- SPACE : toggles simulation.
- 'r' : Resets simulations, and pauses it.
- 'w' : Toggles anti-aliased wireframe rendering.
- 'b' : Toggles rendering of simulated contact
Block is bounded by a disk.
- 'l' : Toggles between single-step and multi-step
Euler integration per frame.
- 'j' : Jiggle--applies random accelerations to objects.
- 'e' : Toggles png-frame exporting to the "frames"
- Test Images: I've
included a number of helpful test images (see below). These are
uncompressed TGA images, with empty space specified by white
(255,255,255) colors. The intensity of the color specifies the mass
density. Blue objects are pinned. Feel free to draw your own elaborate
simulations, and submit them to us for others to try.
- Penalty contact between
two Block objects is currently implemented for you using a simple
damped spring. Familiarize yourself with how contact is handled, and
how the narrow and broad phase detection decides which Block-Block
pairs to apply forces to. Notice that for difficult scenarios, the
contact stiffness must be high, and the time-step size much be small,
yet your objects might still interpenetrate especially with high
rotational velocities. You will also use the simple Block bounding
disks for a robust disk-disk contact model in your velocity-level
As before, the starter code will compile and run using JDK 1.5 or
recommend you get Sun's latest JDK here. In
addition to Java the starter code also uses JOGL/OpenGL and Vecmath
familiar from previous assignments. Although you're welcome to
modify the OpenGL portion of the assignment, it is not necessary to
Assignments Steps: Here
are the steps you need to address, most of which involve modifying the
object to support narrow and broad phase tests:
Profiling: In addition to timing your broad and narrow phase
collision times, you can use the "java
-Xprof ..." option to help identify and remove program
bottlenecks. You should also use the JVM's -server
option to encourage aggressive runtime optimizations.
1. Rigidbody dynamics are
mostly implemented, but some minor things are missing such as
rotational dynamics. As a first step in familiarizing yourself
with the code, add torque in RigidBody.applyContactForceW(Point2d
contactPointW, Vector2d contactForceW) to make the
objects spin; you can do this by taking the z-component of the cross
product of the xy-plane force and relative position, i.e., the triple
scalar product, zHat-dot-(x-p)-cross-force. Also, finish implementing
Symplectic Euler integrator so that it works. (You can find both parts
by searching for "TODO")
2. Narrow phase collision detection: Even with only a couple letters, e.g.,
collision processing is slow due to both the small time-steps used and
the naive all-pairs
collision processing, i.e., each boundary Block on
object "i" is tested against every boundary Block on
object "j" for overlap of their bounding disks. Your first step
is to acclerate
object-object collision detection of overlapping Block/Disk primitives,
and thus speed up the
penalty-based simulator. You
could do this in various ways. First, you could build a bounding
volume hierarchy on each RigidBody's Block
primitives, with the most natural hierarchy for Block primitives being
one with axis-aligned bounding boxes (an AABB tree). Be careful
that the bounds contain the support of any disk-like penalty forces,
etc. Build the narrow phase data structure in
each object's RigidBody constructor.
- Improve mouse picking and
dragging: Use your
narrow phase test to improve the mouse picking to select a (nearby)
on the object rather than just its center of mass. You can do this by
implementing a better RigidBody.intersectsW(Point2d) method to see if a
point intersects a rigid body's blocks, then modifying the mouse spring
force interaction appropriately.
- Support object-object overlap
queries to identify the overlapping Block disks for
object-object pairs. This involves recursively testing one object's AABB-Tree against
another object's tree. Your penalty-based
should now be much faster, and support much more complex scenes.
3. Velocity-level contact
solver: The main part of this assignment is the
implementation of the velocity-level complementarity constraint
solver. You will implement the project Gauss-Seidel solver for
the nonlinear complementarity problem discussed in class [Erleben
2007]. The steps you need to pursue are as follows:
the simple disk-disk contact geometry from the penalty-based solver to
generate your velocity-level contacts. Given two overlapping disks and
the line segment connecting their center points, use the line's
mispoint as the collision point and its orientation as the contact
normal. You will use discrete-time collision detection to find
overlapping pairs (as opposed to continuous/interval-based collision
detection like in assignment #2), which will work reasonably well for
the disk-disk contacts and sufficiently small time-steps.
- Contact Jacobian: Given
the contacts, you should implement basic functionality to evaluate the
contact Jacobian to map between rigidbody velocity and
normal/tangential contact velocity. Note that you need not explicitly
compute a sparse matrix data structure and store the Jacobian. However,
you will want to evaluate columns/rows, and know which contacts each
rigidbody has, etc.
- Setting up w = A lambda + b:
In preparation for projected Gauss-Seidel, you will need to compute the
b vector of unconstrained
contact velocities, and support evaluation of (block) rows of the A matrix. Exploit sparsity by using
the fact that contact k is directly row-influenced by
other contacts involving body i and j.
- Projected Gauss-Seidel (PGS)
complementarity solver: Iteratively solve the complementarity
problem by performing Gauss-Seidel iterations with appropriate clamping
of the Lagrange multipliers. Use a fixed number of PGS iterations, such
as 10 or 100 or 1000. At each iteration you will require a single row
of the A matrix, however you
will find it more convenient to evaluate/access both normal &
tangential rows for contact k at the same time. Large numbers of
iterations will benefit from efficient access and multiplication of the
rows of A.
- Pin constraints: Use
inverse-mass-matrix filtering to support pinned rigid bodies (RigidBody.isPinned()),
by treating pinned objects as having infinite linear and rotational
mass. Observe that these filters will essentially zero-out certain
contributions to the A
matrix---notably its diagonal entries.
"Room for resting beans only, please.":
To make it extra challenging, your friction forces need to keep as many
of your beans from moving as possible. That means no implausibly
sliding, jittering, or chattering beans---these aren't coffee
beans! The simulator will track how many beans are resting at any
and the historical maximum (nRestingMax
It will also estimate a not-so-robust height of your resting jelly bean
) and its historical maximum (Best Height
modify the Jelly Bean Factory code, or
your contest submission will be void. Also, you are not allowed
to simply set the velocity of beans to zero
(which will void your submission). You must rely on the same general-purpose
projected Gauss-Seidel (or better) complementarity-based velocity-level
contact solver used in all examples.
thing... if you can achieve a plausibly impressive >256 resting jelly beans
your frictional contact solver, then you will attain the impressive OCD
Distinction. This will be hard. Very hard. No one has even done it
before. Heck, they won't even fit in the jar. Come to think of it, it's
probably not even possible. You might as well not even try...
|Jelly Bean Factory Contest
(Email png image of simulation to email@example.com for inclusion)
| #Beans | BestHeight
| 346 | 0.51
| 290 | 0.54
Penal T. Bounce
6. Other Things To Try:
In order to simulate a very large pile of jelly beans, or make your
castle destruction animation, you may find that you need a little more
sophistication or other functionality. Here are some things to
Broad phase collision
Once you have an efficient narrow phase processing and the
velocity-level constraint solver working, you can collide
detailed objects together efficiently and your bottleneck will shift to
all-pairs broad phase test for certain scenarios. You can implement any
broad phase collision
detection scheme provided it yields decent performance on the large
examples. Schemes you might consider are uniform spatial subdivision
and related hashing schemes, hierarchical grids (better performance on
variable object sizes), octrees, kd-trees, as well as sweep and prune.
Feel free to use your own code from the second (Spaghetti Factory)
Adaptive Time Stepping:
Given that contacts are detected using discrete collision tests, you
may want to monitor the maximum
Block speed to reduce the time-step size to avoid missing collisions.
Exploiting temporal coherence:
In addition to adaptive time-stepping, you can exploit temporal
coherence by reusing previous values from your PGS solve to "warm
start" solutions at the next time step. This is particularly
useful for quasi-static examples, such as stacks, but it can also be
useful for less contrived examples. One challenge is establishing
contact-contact correspondence between time steps.
Restitution coefficient: The
contact conditions derived in class impose a nonnegative normal
velocity at contacts, however for new impacting contacts you may wish
to impose an impact condition on the normal velocity based on a
Stacking examples are particularly challenging for PGS, and shock
propagation techniques are commonly used to accelerate convergence or
provide more plausible approximate solutions [Erleben 2007]. Try
building a contact graph and using shock propagation to improve the
stability of your stacks, and achieve that extra special OCD
Erleben mentions various enhancements to improve the performance of the
PGS solver. Feel free to incorporate these into your submission.
Other forces: Feel free to add
springs, joint constraints or other rigid body system elements to allow
you to model more interesting mechanisms.
submit a brief
(in txt or PDF format) describing your approach and any findings, in
addition to your Java implementation. Provide videos to
document any results you want us to see, any creative artifacts, and
your best Jelly Bean Factory run, and your best stacking example
runs. If you are working with a partner, be sure to form and
submit your zip file as a group.
in a portable format such as QuickTime, mpg, or divx, but not native
formats, e.g., not the FRAPS codec for your machine.
Please submit videos with a fixed
resolution of 720-by-720 (the default
On collaboration and academic
integrity: You are allowed to collaborate on the assignments to
the extent of formulating ideas as a group, and derivation of physical
equations. However, you must conduct your programming and write up
completely on your own (or with your partner), and understand what you
are writing. You may not use code from the web. Please also list the names of everyone that
you discussed the assignment with. You are expected to
maintain the utmost level of academic integrity in the course. Any
violation of the code of academic integrity will be penalized
- K. Erleben, Velocity-based
shock propagation for multibody dynamics animation, ACM Trans.
Graph. 26, 2, Jun. 2007. (most similar to Mar09 lecture)
A., and Baraff, D., Eds. 2001. Physically Based Modeling: Principles
and Practice. Course Notes. ACM SIGGRAPH '01.
Guendelman, Robert Bridson, Ronald P. Fedkiw, Nonconvex Rigid Bodies With Stacking,
ACM Transactions on Graphics, 22(3), July 2003, pp. 871-878.
- Collision Detection
(discusses broad and narrow phases):
- Bounding volume hierarchies (BVHs):
rigid bodies, you can precompute each rigid body's bounding volume
hierarchy (BVH) in
the frame of reference of the rigid body. Some common choices are:
- Sphere Trees are described in
- AABB Trees are described (although
for deformable models) in:
den Bergen, Efficient Collision Detection of Complex
Deformable Models using AABB Trees, Journal of Graphics
Tools: JGT, 2(4), 1--14, 1997.
rigid bodies there is no need to update
the bounding volume hierarchy every frame, unlike with the deformable
model case in [van den Bergen 1997].
the bounds are axis-aligned in each object's material frame, when
testing two objects' AABB-rees against each other, bounds are rotated
due to the relative rigid transformation between the two objects.
- k-DOPs (generalized box bounds):
- OBB-Trees (oriented box bounds):
Copyright Doug James, March 2009.