COMS W4167: Physically Based Computer Animation
Theme III, Milestone II – Rigid Body Simulation
Due 22:00:00 PM Tuesday 16th Nov 2021 (EST)
Academic Honesty Policy
You are permitted and encouraged to discuss your work with other students. Except where explicitly stated
otherwise, you may work out equations in writing on paper or a whiteboard. You are encouraged to use
Piazza to converse with other students, the TAs, and the instructor.
HOWEVER, you may NOT share source code or hard copies of source code. Refrain from activities or
the sharing materials that could cause your source code to APPEAR TO BE similar to another student’s
enrolled in this or previous years. We will be monitoring source code for individuality. Cheating will be
dealt with severely. Source code should be yours and yours only. Do not cheat. For more details, please
refer to the full academic honesty policy on the departmental website and on Piazza.
1 Introduction
In Theme III Milestone I we derived the equations of motion for a rigid body and explored the simulation
of rigid body dynamics in the absence of collisions. In Theme II we explored various methods for resolving
collisions, and encountered many shortcomings with these methods (recall the problems associated with the
rigidification failsafe). In this milestone we will bring together ideas from these previous milestones and
derive a method for rigid body contact. In particular, the response methods in this milestone can gracefully
handle breaking contact in the face of multiple co-dependent collisions.
The methods in this milestone depend on more advanced numerical techniques than previous assignments;
we will make use of both a linear complementarity problem (LCP) solver, as well as a convex quadratic-
program solver.
Note: This assignment has proven very tricky in the past. However, everything is in the writeup. Many
issues from previous years have been resolved by reading the writeup more closely.
2 New XML Features
This milestone introduces the following new simulation features:
1. The rigidbodycollisionhandling feature specifies the detection and response methods to use:
The detection attribute specifies which collision detection method to use; the only valid value for this
assignment is all-pairs. The response method specifies which response method to use; the only valid
values for this assignment are lcp and velocity-projection.
2. We have added a fixed attribute to the rigidbody feature:
If fixed is set to 1, forces and contact impulses will not effect the body.
1
Figure 1: Breaking contacts in a rigid body simulation. Active contacts are denoted by small stars, breaking
contacts by small circles. Naively resolving collisions in an iterative fashion can produce an incorrect solution.
3 The Linear Complementarity Problem
3.1 Problem Formulation
Consider a system of rigid bodies in contact at k points. In order to prevent interpenetration, we would like
to compute a set of k impulses applied at these contact points that yield separating relative velocities. That
is, for each contact i, we seek impulses that yield ġi = n̂i · (ṙ+ip− ṙ
+
iq) ≥ 0. We denote the relative velocity at
the ith contact point by ġi, the contact normal at the i
th contact by n̂i, and the post-impulse world-space
velocities of rigid bodies p and q at the ith contact by ṙ+ip and ṙ
+
iq (a − superscript defnotes pre-impulse).
Denote the change in velocity due to the action of all impulses at the ith contact point on each body by δṙip
and δṙiq. With this notation, the i
th constraint becomes ġi = n̂i · (ṙ−ip − ṙ
−
iq) + n̂i · (δṙip − δṙiq) ≥ 0 where
ṙ−ip and ṙ
−
iq are the known pre-impulse velocities. Let us examine the change in velocities in more detail.
Denote the ith impulse by λin̂i. Expanding δṙip in terms of the effect of all impulses on the center
of mass velocity and the effect of all impulse-induced torques on the angular velocity, we find that δṙip =
δvp + δωp× r′ip (recall that prime denotes world-space positions relative to the center of mass). To compute
the change in velocity, we simply sum the effect of all impulses on body p (let sjp indicate whether impulse
j acts on body p, and if so capture the correct sign): δvp = M
−1
p
∑k−1
j=0 sjpλjn̂j . To compute the impulse-
induced torque, we sum impulses crossed with their point of action relative to the center of mass of body
p: δωp = I
−1
p
∑k−1
j=0 r
′
jp × sjpλjn̂j . Making these substitutions, we have: δṙip = M
−1
p
∑k−1
j=0 sjpλjn̂j +
I−1p
∑k−1
j=0 r
′
jp × sjpλjn̂j × r
′
ip.
Armed with the knowledge of how all impulses change the velocity at the points of contact, we can expand
the constraint as:
ġi = n̂i · (ṙ−ip − ṙ
−
iq) + n̂i · (δṙip − δṙiq)
=
k−1∑
j=0
λj
(
sjpn̂
T
i M
−1
p n̂j
)
+
k−1∑
j=0
λj
(
sjpn̂
T
i I
−1
p r
′
jp × n̂j × r
′
ip
)
−
k−1∑
j=0
λj
(
sjqn̂
T
i M
−1
q n̂j
)
−
k−1∑
j=0
λj
(
sjqn̂
T
i I
−1
q r
′
jq × n̂j × r
′
iq
)
+ n̂i · (ṙ−ip − ṙ
−
iq)
Observe that this constraint is linear in all impulse magnitudes λj . We can therefore express the contact
constraints as a linear function of the impulse magnitudes: Aλ + b. The ith entry of this vector gives the
post-impulse relative velocity at the ith contact. The non-penetration constraints are thus:
Aλ+ b ≥ 0 (1)
2
n̂0 n̂1 g
Figure 2: Block resting on table with center of mass over the edge.
Our work is not yet complete, however. As currently formulated, nothing prevents constraint impulses
from pulling. Unless we are modeling an adhesive surface, this behavior is undesirable, so we impose the
additional constraint that contact forces can only push:
λ ≥ 0 (2)
We need to introduce one final constraint before our work is complete, however. We want to enforce the
fact that conditions that aren’t currently violated should cause no impulse (i.e. λ = 0). If the motion of
two objects is breaking away, or they are not in contact, our formulation should not apply an impulse. This
often happens when resolving one constraint eliminates the need for another constraint to act (Figure 1).
Otherwise, if we ignore the fact that no impulses should be applied along breaking contacts, we could
potentially apply a response that is too strong and add energy to the system. Thus, we introduce the
condition that:
(Aλ+ b)iλi = 0 ∀i (3)
What does this condition say? If λi > 0, then (Aλ + b)i = 0. That is, if an impulse is acting at some
point, the bodies must remain in contact at that point. If (Aλ + b)i > 0, then λi = 0. That is, if bodies
are separating, no impulse can act. This celebrated result is known as the Signorini-Fichera Condition.
Together, these three conditions on contact constraints can be solved as a linear complementarity problem.
3.2 Example System
Before proceeding, we will first examine an example problem and its solution.
3.3 Block On Table Edge With Two Contacts
Consider a block sitting on a fixed table with 2 contact points, and with the block’s center of mass positioned
beyond the edge of the table. See Figure 2. Labeling the rigid body A, we can write the constraints on the
relative velocity as:
ġ0 = n̂0 · ṙ+0A = n̂0 · ṙ
−
0A + n̂0 · δṙ0A ≥ 0
ġ1 = n̂1 · ṙ+1A = n̂1 · ṙ
−
1A + n̂1 · δṙ1A ≥ 0
Expanding δṙ0A and δṙ1A we find
δṙ0A = δvA + δωA × r′0A
δṙ1A = δvA + δωA × r′1A
where primes denote the positions of the contacts relative to the body’s center of mass. δvA is the change
in the center of mass’ velocity due to both impulses, and δωA is the change in angular velocity about the
3
Figure 3: A table under gravity in two dimensions with three contacts. Observe that multiple sets of impulses
give the same solution – we could apply only the large impulse to the center leg, or the two smaller impulses
to the outer legs.
center of mass due to both impulses. Expanding δvA and δωA we find:
δvA = λ0n̂0/MA + λ1n̂1/MA
δωA = r
′
0A × λ0n̂0/IA + r
′
1A × λ1n̂1/IA
Substituting these expressions into the constraints, we find
ġ0 = n̂0 · ṙ−0A + λ0 (n̂0 · n̂0/MA + n̂0 · r
′
0A × n̂0 × r
′
0A/IA) + λ1 (n̂0 · n̂1/MA + n̂0 · r
′
1A × n̂1 × r
′
0A/IA) ≥ 0
ġ1 = n̂1 · ṙ−1A + λ0 (n̂1 · n̂0/MA + n̂1 · r
′
0A × n̂0 × r
′
1A/IA) + λ1 (n̂1 · n̂1/MA + n̂1 · r
′
1A × n̂1 × r
′
1A/IA) ≥ 0
or in matrix form(
n̂0 · n̂0/MA + n̂0 · r′0A × n̂0 × r
′
0A/IA n̂0 · n̂1/MA + n̂0 · r
′
1A × n̂1 × r
′
0A/IA
n̂1 · n̂0/MA + n̂1 · r′0A × n̂0 × r
′
1A/IA n̂1 · n̂1/MA + n̂1 · r
′
1A × n̂1 × r
′
1A/IA
)(
λ0
λ1
)
+
(
n̂0 · ṙ−0A
n̂1 · ṙ−1A
)
≥
(
0
0
)
Let us substitute in specific values and examine the solution. Let n̂0 = (0, 1, 0)
T , n̂1 = (0, 1, 0)
T , MA = 4,
IA = 20, r
′
0A = (−2,−1, 0), r
′
1A = (−1,−1, 0), ṙ
−
0A = (0,−2, 0), and ṙ
−
1A = (0,−2, 0). This gives the system:(
9/20 7/20
7/20 3/20
)(
λ0
λ1
)
+
(
−2
−2
)
≥
(
0
0
)
which is solved by λ0 = 0 and λ1 ≈ 13.333 (as a sanity check, verify that these values satisfy the three
constraints in the LCP). What does this mean? λ0 = 0 implies that the leftmost impulse is not active and
the contact is breaking. λ1 ≈ 13.333 implies that the second point remains in contact with the table, and
that an impulse is acting to prevent interpenetration. The net effect is that the block will pivot about the
end of the table, as expected.
4 Redundant Constraints
While our formulation is enough to obtain a correct solution, it does not pin down a unique solution. For
example, Figure 3 illustrates a table with three legs resting on the ground. There are a (uncountably
infinite) number of solutions to this contact problem as we have phrased it. We could apply a single impulse
of magnitude I to the center leg and no impulses to the side legs, we could apply two impulses of magnitude
I/2 to the side legs, or we could apply any combination of (pushing) impulses with total magnitude I that
induce no torque.
What does this redundancy mean in linear algebra terminology? For the three legged table, the matrix A
is singular. That is, Aλ+b = 0 has multiple solutions. Many of the standard LCP algorithms, unfortunately,
will fail when presented with a singular matrix. While one can construct a number of interesting problems
that do not have this singular matrix, if we desire a truly robust method we will have to work a little harder.
4
It is important to stress that this is a numerical issue. Conceptually (numerics aside), even though λ is not
uniquely defined, the final velocities are uniquely defined. Luckily, this problem can be solved by a number
of numerical methods, some of which are more robust to this redundancy than the LCP. We will repose this
problem as a minimization, but first let us formalize some of the concepts that we explored above, using a
notation that applies equally to the LCP and minimization viewpoints.
5 More General Notation
The notation employed thus far has been traditionally used in the graphics literature. While the derivations
so far have been straightforward in that they follow immediately from writing down the physics of rigid
bodies and ‘turning the crank’, they can be a tad unwieldy.
Before rephrasing the contact problem to address the singular matrix A, let us take a moment to construct
a more general notation than the one employed above. Rewriting the notation will have a few advantages.
First, the more general notation will make our results easier to carry over to other settings (e.g. simulation
of deformable bodies). Second, the more general notation will make higher-level algebraic manipulations
easier, and makes it easier to see higher-level properties of the system, such as the symmetry of A.
5.1 Notation
Take each rigid body’s reduced or generalized coordinate representation (center of mass position, orientation),
and concatenate them into a single column vector x. Similarly, concatenate the rate of change of the reduced
coordinates into a single column vector ẋ. As an example, for our 2D representation or rigid bodies, x ∈ R3N
and ẋ ∈ R3N , where N is the number of rigid bodies. For 3 rigid bodies these vectors look like:
x =
(
X0x X
0
y θ
0 X1x X
1
y θ
1 X2x X
2
y θ
2
)
ẋ =
(
V 0x V
0
y ω
0 V 1x V
1
y ω
1 V 2x V
2
y ω
2
)
When applying forces and impulses to rigid bodies, we need some way to compute world-space positions
and velocities from reduced positions and velocities. Let ri(x) be a function that computes the world space
position of the ith point from the reduced representation. For example, if the seventh particle is attached to
the first rigid body, in our representation this function takes the form (tilde denotes body space):
r7(x) = R(θ
1)r̃7 + X
1
To compute the world-space velocity of the ith particle, we simply employ the chain rule. That is, ṙi(x) =
∇riẋ. Continuing with the above example, ∇r7 ∈ R2×9 is given by:
∇r7 =
(
0 0 0 1 0 − sin θ1r̃7x − cos θ1r̃7y 0 0 0
0 0 0 0 1 cos θ1r̃7x − sin θ1r̃7y 0 0 0
)
Convince yourself that the multiplication ∇riẋ gives the same velocity we derived in Theme III Milestone I.
Consider now a set of contacts C. For some contact k ∈ C, the contact occurs between two points ri ∈ R2
and rj ∈ R2. Let the normal for this contact, n̂k ∈ R2, point from i to j. We can thus express the relative
velocity at this contact as n̂Tk (ṙi − ṙj) = n̂
T
k (∇ri −∇rj)ẋ ≡ n̂
T
k Γkẋ. For our simulations, Γk ∈ R
2×3N , and
contains 2 non-zero R2×3 blocks. Continuing our example of three rigid bodies, if the 3rd contact is between
the 7th point (on body 1) and the 9th point (on body 2), Γ3 is given by:
Γ3 =
(
0 0 0 1 0 − sin θ1r̃7x − cos θ1r̃7y −1 0 sin θ2r̃9x + cos θ2r̃9y
0 0 0 0 1 cos θ1r̃7x − sin θ1r̃7y 0 −1 − cos θ2r̃9x + sin θ2r̃9y
)
If we want to apply an equal and opposite impulse y to points i and j at contact k, the resulting impulse
on the reduced coordinates is given by ΓTk y. Similarly, we can map the contact normals to their reduced
coordinate counterparts by ηk = Γ
T
k n̂k ∈ R
3N . Finally, concatenating these normals into one matrix, we
have N =
(
η0 η1 . . . η|C|−1
)
∈ R3N×|C|.
5
qx
qy q̇y
q̇x
Figure 4: On the left, unit-mass particles moving in configuration space constrained to lie above the horizon-
tal. On the right, the same particles in velocity space and the constraint impulses. Post-impulse velocities
are denoted by stars. Notice that the post-impulse velocities are given by the closest point in the admissible
region.
5.2 A Revisited
Let us revisit the construction of the system Aλ+ b for the LCP. The post-impulse relative velocity at all
contact points in world-space can be expressed in terms of the reduced coordinates as NT ẋ+. Mapping the
impulse magnitudes in world-space to the reduced coordinate space gives Nλ, and the corresponding changes
in velocities are expressed as M−1Nλ. Thus, we find that NT ẋ+ = NT (ẋ− + M−1Nλ) = NTM−1Nλ +
NT ẋ− ≥ 0. That is, A = NTM−1N and b = NT ẋ−. The symmetry of A now follows immediately from
the symmetry of M (take the transpose of NTM−1N and see what you get).
6 Velocity-Projection for Inelastic Contact
When solving the LCP, we encountered problems with a singular matrix that arose due to redundant contact
directions. However there is a way to rephrase our problem as a minimization problem. More specifically, a
convex quadratic-program, an optimization problem for which algorithms exist that are robust in the face
of redundant constraints (see [1]).
How do we pose our problem as a minimization? We want our contact to be resolved in a way that
changes velocity the least. This will preserve the physical properties of the system as much as possible while
ensuring no penetration occurs. We do this by minimizing the kinetic energy of the velocity difference. Think
about it like this: we can’t conserve energy, but we can try to change it as little as necessary.
If our initial velocity is v− and our resulting velocity is v+, let us define the difference as δv (ref. Figure
5). Minimizing the kinetic energy of the difference then is: argmin
δv
1
2
δvTMδv.
The constraint that we have is non-penetration. This means 0 ≤ v+Tn = v−Tn+ δvTn = c(δv).
The Lagrangian function for this problem is L(δv, λ) = 1
2
δvTMδv − λ(v−Tn+ δvTn).
Optimization theory then tells us that the necessary conditions for a minimum are the so called KKT
(Karush-Kuhn-Tucker) conditions (ref. [2, p. 321]). They are:
6
V+
V-
V
n
Figure 5: A sketch of our minimization problem.
∇δvL = 0
c(δv) ≥ 0
λ ≥ 0
λc(δv) = 0
(4)
Substituting our definitions in, this means Mδv − λn = 0, so:
δv = M−1λn (5)
v−
T
n+ δvTn ≥ 0 (6)
λ ≥ 0 (7)
λ(v−
T
n+ δvTn) = 0 (8)
6.1 Connection to LCP
In fact we see that this gives the same result as the LCP formulation. (5) is equivalent to the LCP definition of
the impulses, (6) is equivalent to (1) (the non-penetration condition), (7) is equivalent to (2) (non-adhesion)
and (8) is equivalent to (3) (the Signorini-Fichera Condition).
As we would expect, the KKT conditions guarantee that impulses are only active for velocities that
are bringing particles into the inadmissible region. Further, notice that impulses only act in the positive
constraint gradient direction, that is they are only pushing. Finally, observe that all post-impulse relative
velocities are separating or zero. Thus, all of the content of the LCP formulation is reflected in the velocity-
projection solution.
6.2 General formulation of the problem
As a general formulation, the post-collision velocity is given by the minimizer of δT = 1
2
(v− ẋ−)TM(v− ẋ−)
subject to NTv ≥ 0. Expanding this multiplication, we find that δT = 1
2
vTMv− 1
2
vTMẋ−− 1
2
(ẋ−)TMv−
1
2
(ẋ−)TMẋ− = 1
2
vTMv− vTMẋ− − 1
2
(ẋ−)TMẋ− where we have used the symmetry of M to combine the
7
cross terms. We can drop the constant term − 1
2
(ẋ−)TMẋ−, as constant offsets will not change the minimum
of the function. Thus, the post-collision velocities are given by:
ẋ+ = argmin
v
(
1
2
vTMv − vTMẋ− : NTv ≥ 0
)
7 Required Features
7.1 Rigid-Body Rigid-Body All-Pairs Detection
To detect collisions, please complete the detectCollisions method of RigidBodyAllPairsCollisionDetector.cpp.
As input this method takes a vector of all rigid bodies in the simulation, and as output constructs a set of
template type RigidBodyCollision containing all collisions in the system. For each pair of rigid bodies, detect
any vertex-edge collisions. A vertex and an edge are considered colliding if the distance between the two is
less than (keep this inequality strict!) the sum of the radii of each body. You should be able to adapt this
code from Theme II. When inserting a RigidBodyCollision into the collision set, you will have to compute:
1. The index of the first rigid body of the collision (i in the constructor)
2. The index of the second rigid body of the collision (j in the constructor)
3. The vector from the center of mass of the first body to the point of contact (r0 in the constructor)
4. The vector from the center of mass of the second body to the point of contact (r1 in the constructor)
5. The unit-length contact normal (nhat in the constructor)
7.2 LCP
Please complete the resolveCollisions method of RigidBodyLCPCollisionResolver.cpp. As input this method
takes both a vector of all rigid bodies and a set of all collisions between rigid bodies, and as output modifies
the velocities and angular velocities of the rigid bodies so as to prevent penetration.
You will have to compute the matrix A and the vector b described in Sections 3 and 5.2. We have
provided a function lcputils::solveLCPwithODE that takes as inputs A and b and returns a vector λ of
impulse magnitudes (please see the code for an example). You will then have to apply the resulting impulses
to the rigid bodies.
We have provided the expected values for A and b in the scene files collisionslcp/test01.xml and col-
lisionslcp/test03.xml. This code uses the implementation of Dantzig’s method from the open source rigid
body engine ODE.
Note on fixed rigid bodies: Impulses have no effect on fixed rigid bodies. Therefore, to handle fixed rigid
bodies in the response, when you are constructing A ignore any terms corresponding to an impulses’ effect
on a fixed body.
7.3 Velocity-Projection
Please complete the resolveCollisions method of RigidBodyVelocityProjectionCollisionResolver. As input
this method takes both a vector of all rigid bodies and a set of all collisions between rigid bodies, and as
output modifies the velocities and angular velocities of the rigid bodies so as to prevent penetration.
You will have to compute the mass matrix M, the product −Mẋ−, and the generalized constraint matrix
N described in Sections 5 and 6. We have provided a function solve quadprog that takes these quantities as
inputs and returns the solution to the QP. Please see the code for an example.
We have provided the expected values for M, −Mẋ−, and N in the scene files collisionsvelocityprojec-
tion/test01.xml and collisionsvelocityprojection/test05.xml. This code uses the open-source QuadProg++
implementation of the Goldfarb-Idnani algorithm.
8
Note on fixed rigid bodies: The easiest way to handle fixed bodies in the solve is to simply not expose
those bodies’ degrees of freedom to the QP solver. That is, shrink M, Mẋ−, and N to not include fixed
degrees of freedom, and do not add fixed bodies’ contributions to N.
8 Extra Credits: Elastic Response
The LCP method as we formulated above attempts to find a valid (non-pulling, normal) impulse such
that the post-response velocity does not violate the collision constraint. The Velocity Projection method
above attempts to project the velocity back to the admissible domain too, removing any constraint-violating
component but nothing more. As a result, both responses are fully inelastic by construction, which is
apparent from the NT q̇+ ≥ 0 constraint common to both methods. An elastic response can be obtained
by changing this constraint from “the post-response velocity cannot be approaching” to requiring that “the
post-response velocity is at least as much departing as the pre-response velocity was approaching”, i.e.
NT q̇+ ≥ −NT q̇−. We refer the readers to recent work by [3] for more details.
However, as soon as we move into the elastic regime, the multiple simultaneous contact problem be-
comes much more complicated, as pointed out in [3]. It is straightforward to list a small set of fundamental
properties that a desired collision response method should possess, such as symmetry preservation, energy
conservation and not creating artificial sticking between objects, but none of the existing methods prior to [3]
was able to achieve all of these desiderata. Fortunately, with a critical observation on the property of LCP,
it is possible to design a method satisfying all the desiderata, through relatively minor changes to what we
have built in the previous sections. Simply speaking, LCP/Velocity Projection should be applied prudently
and repeated when necessary.
For this section, please carefully read [3], and
1. implement fully elastic response in both LCP and Velocity Projection methods. The formulation in [3]
is general for any coefficient of restitution between 0 and 1, but for this milestone we are only interested
in the fully elastic case.
2. implement the Generalized Reflection algorithm (algorithm 2 in [3]) on top of both LCP and Velocity
Projection methods. The resulting collision response algorithm should have all the desired properties.
All of the implementation in this section should be done in source files RigidBodyGRLCPCollisionRe-
solver.cpp/.h and RigidBodyGRVelocityProjectionCollisionResolver.cpp/.h, separate from the inelastic LCP
and Velocity Projection algorithms in the previous sections. Two new collision response methods, “gr-lcp”
and “gr-velocity-projection”, have been added to the “response” property of the rigidbodycollisionhandling
XML tag.
You can find the tests in the t3m2 extra credit folder to test the various properties of your collision
response implementation. You will be awarded 15% extra credit for a correct Generalized Reflection imple-
mentation that passes all these tests. Most of these tests are taken from [3], therefore, in addition to the
oracle’s visual output, you also have the paper’s video as a qualitative reference to compare your animation
to: http://www.cs.columbia.edu/cg/rosi/rosi.mov. You may find it useful to watch a lecture given by
the author here.
Note on line 7 of Algorithm 2 in [3]: Since our LCP solver and Quadratic Programming solver, which
GR uses as building blocks, have limited numerical precision, a small epsilon should be used for this test.
Specifically, please only proceed into the true branch if the pre-response normal velocity (left hand side of
the inequality on line 7) is smaller than −10−14.
9 Creative Scene
As part of your final submission for this milestone, please include a scene of your design that best shows off
your program. Based on the quality of your scene, you will have the opportunity to earn up to 10% extra
9
http://www.cs.columbia.edu/cg/rosi/rosi.mov
credit. Your scene will be judged by a secret committee of top scientists using the highly refined criteria of:
1. How well the scene shows off this milestone’s ‘magic ingredients’ (a la Iron Chef ).
2. Aesthetic considerations. The more beautiful, the better.
3. Originality.
Top examples will be posted to Piazza.
To submit this scene, place the XML file in the assets/Creative/ directory of your submission, and
submit a video of the scene to Courseworks.
To export PNGs of your scene, run
path/to/FOSSSim -s path/to/creativescene.xml -m path/to/output/dir/
This will output a PNG file for each frame of the simulation. The program will output black frames if you
pass it -d 0 so leave that flag out.
Then, create a video from your PNGs using ffmpeg
ffmpeg -r 24 -f image2 -i pngs/frame%05d.png -vcodec libx264 -pix_fmt yuv420p out.mp4
If you are having issue with saving using svgs, please refer to Piazza post
9.1 Creative scene submission
Please download and submit the mp4 file to Courseworks. Remember to generate videos of your creative
scene before submitting, as the Codio box will become read-only after you submit. The due date for creative
scene will be on 12:00pm (noon) Oct. 21. Please refer to announcement on Courseworks
10 FAQ
Theme 3 Milestone 2 FAQ:
Q: What is the Matrix M referred to in the writeup?
A: This is the diagonal ’Mass’ matrix. In our case the matrix M is a diagonal matrix with repeating sets of
three entries per body: M, M, I (so the third entry per body is the moment of inertia).
Q: What are the vectors r and ṙ?
A: r = R(θ)r̃ + CoM (like milestone 1, so no third dimension)
ṙ is the corresponding (∇r)q̇ (which has a third dimension)
Q: Should we include fixed bodies in our mass matrix M?
A: No. The mass matrix should be square with length 3(N −K), with N being the number of bodies, and
K being the number of fixed bodies.
Q: Whats the best way to debug penetration for Velocity Projection?
A: Velocity projection method is a constrained optimization. There are naturally two kinds of bugs: ones
that violate the constraints, and ones that fails the optimization. The first kind of bugs can be spotted
by observing that your solution doesn’t satisfy the constraint inequalities. The second kind of bugs can be
spotted by checking the direction of the gradient of the objective function (in our case the norm of velocity
change).
If you are seeing penetration it is a violation of the constraint “normal velocity ≥ 0”. You can try to find
a minimal example that reproduces this violation and verify your implementation manually.
Notes:
10
https://ffmpeg.org/ffmpeg.html
https://piazza.com/class/kta7z4zpa6l1ob?cid=189
• Your Mass matrix should not have negative entries in it.
• Make sure your impulses only get applied once. People have had issues in the past were they applied
it multiple times on the same collision point.
• Be very very careful about what calculations are done in body space coordinates and which are world
space.
• A key to this milestone is being consistent with a choice of sign in collision detection normal direction.
All will work out in the math if you stay consistent.
• A zero rotation is an identity matrix, not all zeros.
• The outgoing velocity of two colliding objects is dependent on their individual masses. Conserve
Momentum.
• It is extremely useful for understanding and debugging purposes to work out the dimensions of your
matrix algebra by hand and confirm the code matches your results.
• Be sure to note whether the variables passed into a function you are using are using body-coordinates
or world-space as a reference frame. This is often a cause of error.
• Having a velocity of zero is not the same as being a fixed body.
• Lambda only tells you the strength of the impulse that should be applied at that contact point, so you
should ensure that it’s being applied in the correct direction for each body.
References
[1] D Goldfarb and A Idnani. A numerically stable dual method for solving strictly convex quadratic
programs. Mathematical Programming, (27), 1983.
[2] Jorge Nocedal and Stephen J Wright. Nonlinear Equations. Springer, 2006.
[3] Breannan Smith, Danny M Kaufman, Etienne Vouga, Rasmus Tamstorf, and Eitan Grinspun. Reflections
on simultaneous impact. ACM Transactions on Graphics (TOG), 31(4):106, 2012.
11
Introduction
New XML Features
The Linear Complementarity Problem
Problem Formulation
Example System
Block On Table Edge With Two Contacts
Redundant Constraints
More General Notation
Notation
A Revisited
Velocity-Projection for Inelastic Contact
Connection to LCP
General formulation of the problem
Required Features
Rigid-Body Rigid-Body All-Pairs Detection
LCP
Velocity-Projection
Extra Credits: Elastic Response
Creative Scene
Creative scene submission
FAQ