CS计算机代考程序代写 arm COMS W4167: Physically Based Computer Animation

COMS W4167: Physically Based Computer Animation
Theme III, Milestone I – Rigid Body Simulation

Due 22:00:00 PM Tue 9nd 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 we will explore methods for simulating rigid bodies – objects in which the distance between
any two points is fixed for all time. Rigid bodies have a number of interesting applications ranging from
their use in robot motion planning to studying the time evolution of tops. Rigid bodies are also one of the
most common primitives simulated in video game physics engines.

Rigid bodies see heavy use as they posses many favorable numeric properties; simulating similar systems
with masses and springs or finite elements will lead to very stiff systems and hence strict time step limits.
As a reduced coordinate representation, rigid bodies also have the potential to yield large savings due to the
reduction in the number of degrees of freedom.

The use of reduced or generalized coordinates can be viewed as a constraint enforcement method. For
example, one can encode the position of a pendulum with one variable that changes with time – the deflection
of the pendulum arm from the vertical. The position of the pendulum bob can be recovered using simple
trigonometry, the angle of deflection, and the length of the arm. By construction, the inextensibility of
the arm cannot be violated. In contrast, if one were to use cartesian coordinates, the inextensibility of the
constraint would have to be enforced in some manner. Usually reduced coordinates are difficult to write
down for complex constrained systems, but for rigid bodies they prove both manageable and advantageous.

In this first milestone we will explore the dynamics of rigid bodies in the absence of collisions.

2 New XML Features

This milestone introduces the following new simulation commands:

1. The simtype node specifies the type of this simulation:

The valid values for type are “particle-system” and “rigid-body”. Setting the type to “particle-system”
indicates that this scene file describes a particle simulation as implemented in the previous two themes.
In this theme, we will use a setting of “rigid-body”.

2. The rigidbodyintegrator node specifies both the time integration technique and the time-step to use for
this rigid body simulation:

1

0

1

2

34

Figure 1: The hull of a rigid body. Note the clockwise specification of vertices.

The type attribute specifies which time integration technique to use: valid values are “explicit-euler”
and “symplectic-euler”. The dt attribute specifies the time-step to use. dt expects a positive scalar.

3. The rigidbodyvertex node creates a new vertex that the user can assign to a rigid body:

The x and y attributes specify the position of the vertex. Both x and y expect scalar values. The m
attribute specifies the mass of the vertex and must be a positive scalar.

4. The rigidbody node creates a new rigid body:

Each p attribute references a rigidbodyvertex and adds a new vertex to this rigidbody ’s hull. The
vertices that define the (generally non-convex) hull of this rigidbody are specified in clockwise order.
The first and final vertex are also connected. See Figure 1. The center of mass, the total mass, and
the moment of inertia of the rigid body are computed from the positions and masses of these vertices.
The vx and vy attributes specify the initial velocity of the center of mass of this rigidbody. The omega
attribute specifies the initial angular velocity about the center of mass of this rigidbody. The r attribute
specifies the radius of this rigidbody for rendering and collision-detection purposes.

5. The rigidbodyspringforce node creates a spring force that acts on rigid bodies:

The i and j nodes specify which rigid bodies this force acts between. If either value is set to −1, this
indicates that the corresponding endpoint is fixed in space and not attached to a rigid body. If an
endpoint is attached to a rigid body, the pix, piy, pjx, and pjy nodes specify where on each rigid body
the spring is attached relative to the center of mass (assuming the original orientation). If instead
the corresponding endpoint is specified as fixed, these values specify the position in space where the
spring endpoint is fixed to. The k node specifies the stiffness of the spring and the l0 node specifies
the rest-length of the spring.

6. The rigidbodygravityforce node creates a (constant near-earth) gravity force that acts on all rigid bodies:

2

r 0i

r i

r i

Figure 2: The body space and world space frames. r′i = R(t)r
0
i and ri = R(t)r

0
i + X(t).

The fx and fy attributes specify the x and y components of the force, respectively.

7. The rigidbodywindforce node creates a wind-like force that acts on all rigid bodies:

The positive scalar beta attribute scales the overall strength of the force. The fx and fy attributes
specify the x and y components of the force, respectively. The integer pointsperedge attributes specifies
the number of quadrature points per edge to use when computing this force (see description below).

This milestone introduces the following new rendering commands:

1. The rigidbodycolor node sets the color of an entire rigid body:

The body attribute specifies which rigid body this color applies to. The r, g, and b nodes specify the
color of the body. Their values must be between 0 and 1.

2. The rigidbodyspringcolor node sets the color of a spring:

The spring attribute specifies which spring this color applies to. The r, g, and b nodes specify the
color of the spring. Their values must be between 0 and 1.

3 Rigid Body Dynamics

A rigid body is a collection of masses (or a continuum, but here we consider only a finite number of points)
in which the distance between any two masses remains fixed for all time. Intuitively, the two transformations
that preserve this invariant are translations and rotations. Therefore, not surprisingly, we will find that we
can decompose the dynamics of a rigid body into two components: a component associated with the center
of mass that behaves like a point mass, and a component associated with rotations about the center of mass.

3.1 Body Space, World Space, and Center of Mass

Consider a 2D rigid body composed of N points masses lying in a 2D plane. Recall that we define the center
of mass of a collection of masses as

X =


imiri∑
imi

=


imiri
M

3

ω

r i

pperp

ppara

ṙ i

φ

Figure 3: A rigid body rotating about its center of mass and axis ω.

where mi is the mass of the i
th particle, ri is the position in world space of the i

th particle, and M =

imi
is the total mass.

Define body space as a 2D orthonormal coordinate system with the origin placed at the body’s center of
mass and in which the body is considered to have no rotation (see Figure 2). If r0i is the position of any
point on the rigid body in body space, then the point’s position in world space is given by

ri = R(t)r
0
i + X(t) (1)

where R(t) is a 2× 2 rotation matrix describing the orientation of the rigid body at time t, and X(t) is the
position of the body’s center of mass at time t. As R is purely a rotation, we can simply store a single angle
θ to reconstruct R.

What does this construction gain us? At any time t, if we know the orientation of the body θ and
the center of mass’ position X, we can compute the position of any particle on the rigid body. Instead of
simulating 2N degrees of freedom, we need only simulate 3. Furthermore, as our representation only admits
rigid motions, that is rotations and translations, we do not need to enforce the rigidity constraint with soft
constraints (e.g. stiff forces), hard constraints (e.g. Lagrange multipliers), or some other technique.

Before we can simulate a rigid body with this reduced coordinate representation, however, we need to
ascertain how the velocity of a particle on the body relates to this reduced representation, and how forces
act on this reduced representation.

3.2 Linear and Angular Velocity

We can compute the velocity of a particle attached to a rigid body by taking the time derivative of (1):

ṙi = Ṙ(t)r
0
i + Ẋ(t)

Note that r0i is fixed throughout the simulation and thus its time derivative is 0. Ẋ is simply the velocity of
the rigid body’s center off mass, but how do we compute Ṙ? Before proceeding, let us examine the columns
of R.

Component-wise, we can write R as: (
Rxx Ryx
Rxy Ryy

)
Consider the action of R on the cartesian x-axis:

Rx̂ =

(
Rxx Ryx
Rxy Ryy

)(
1
0

)
=

(
Rxx
Rxy

)

4

That is, the first column of R is simply the body space x-axis rotated into world space! Similarly, we find
that the second column of R is the body space y-axis rotated into world space. Therefore, if we can compute
how these vectors evolve instantaneously, we can compute Ṙ. In order to compute these quantities, we will
first derive a more general 3D result.

Consider a rigid body rotating about its center of mass and some axis ω with angular velocity ω = |ω|
(see Figure 3). Consider a vector pointing from the center of mass to some point on the body, say r′i =
Rr0i = (ri − X) = p

perp + ppara, where pperp is perpendicular to ω and ppara is parallel to ω. ppara is
unaffected by the rotation, while pperp will trace out a circle. Simple trigonometry gives that radius of the
circle is |ri −X| sinφ, where φ is the angle between (ri −X) and ω. Thus, the instantaneous speed of any
point on this circle is given by |ω||ri −X| sinφ. Furthermore, the direction of the velocity of this point is
perpendicular to both ω and (ri −X). What operation gives this magnitude and this direction? The cross
product! Therefore, the velocity of ri −X is given by ω × (ri −X).

As a notational convenience, note that we can represent a cross product as a matrix-vector multiplication:

a× b =


aybz − azbyazbx − axbz
axby − aybx


 =


 0 −az ayaz 0 −ax
−ay ax 0




bxby
bz


 = a∗b

Returning to the problem of computing Ṙ, first note that translations of the rigid body do not change
R, and thus we need only consider changes to R due to the angular velocity. Invoking the above result:

Ṙ =

(
ω ×

(
Rxx
Rxy

)
ω ×

(
Ryx
Ryy

))
= ω∗R

Thus, the world space velocity of ith particle in the rigid body is ṙi = ω(t)
∗R(t)r0i + Ẋ(t). That is, if

in addition to X(t) and θ(t) we know Ẋ(t) and ω(t), we can recover the velocity of any point on the rigid
body! Therefore, we can represent the entire state of the rigid body with these four quantities.

If we wanted to simulate rigid bodies without any forces, the above would be enough; Ẋ and ω would
remain constant throughout the simulation, and we would just update X(t) and θ(t) using these constant
values. This situation is not terribly interesting, however, so we will now see how Newton’s second law
(F = ma) looks when phrased in terms of our reduced coordinates.

3.3 Forces and Torques

3.3.1 Center of Mass Acceleration

To determine how a force effects the center of mass’ velocity, Ẋ, we will start by computing the total
momentum of the rigid body. The total momentum P of a rigid body is given by:

P =

i

miṙi = (

i

mi)


imiṙi∑
imi

= MẊ

That is, the momentum of the system can be computed by viewing all of the body’s mass as moving with
the center of mass. Taking the derivative of this quantity gives us the center of mass’ acceleration:

MẌ =

i

mir̈i =

i

fi

fi is the total force acting on mass i. We have to be a little careful going forward, as this force includes both
external forces (e.g. gravity) as well as the internal forces acting within the rigid body. Let us assume that
each point j exerts a force on each other point i to maintain the rigidity constraint. That is, we can write
the force on particle i as fi = f

ext
i +


j 6=i fij where f

ext
i is the sum of the external forces acting on i, and fij

is the internal force on i from point j. Substituting this sum into the rate of change of momentum:

MẌ =

i

fi =

i

(fexti +

j 6=i

fij) =

i

fexti +

i


j 6=i

fij

5

Let us examine

i


j 6=i fij in more detail. We can rewrite this sum and invoke Newton’s third law (‘for

every action there is an opposite and equal reaction’), giving:∑
i


j 6=i

fij =

i


j>i

(fij + fji) =

i


j>i

(fij − fij) = 0

That is, provided the internal constraint forces obey Newton’s third law, to compute the acceleration of the
center of mass, we simply sum the external forces acting on each particle of the body:

MẌ =

i

fexti = F

3.3.2 Angular Acceleration

To determine the angular acceleration of a rigid body, we will first compute the time derivative of the angular
momentum measured with respect to the origin. The total angular momentum of rigid body is given by

L =

i

ri × pi

and thus we compute the time derivative as:

L̇ =

i

ṙi × pi +

i

ri × ṗi

In the first term, the velocity of a particle is parallel to its momentum, so this term is 0. In the second term,
ṗi = fi = f

ext
i +


j 6=i fij is the total force acting on the point, which as in the previous section includes both

the external forces on the point and the internal constraint-maintaining forces. Making this substitution, we
find:

L̇ =

i

ri × ṗi =

i

ri × (fexti +

j 6=i

fij) =

i

ri × fexti +

i

ri ×

j 6=i

fij

Let us examine the second term in more detail. We can rearrange the sum and invoke Newton’s third law
to obtain: ∑

i

ri ×

j 6=i

fij =

i


j 6=i

ri × fij =

i


j>i

(ri × fij + rj × fji) =

i


j>i

((ri − rj)× fij)

When is this sum always 0? When ri − rj is parallel to fij . That is, when all constraint forces are central
and obey Newton’s third law, the time rate of change in angular momentum is simply the sum of all external
torques:

L̇ =

i

ri × fexti

We can further decompose the angular momentum into a component associated with the center of mass’
motion and a component associated with motion (spin) about the center of mass. Let r′i = R(t)r

0
i denote

the body space coordinate of particle i rotated into world space by the body’s orientation. Substituting a
particle’s position expressed relative to the center of mass into the formula for angular momentum:

L =

i

ri × pi =

i

(X + r′i)×mi(Ẋ + ṙ

i)

=

i

X×miẊ +

i

X×miṙ′i +

i

r′i ×miẊ +

i

r′i ×miṙ

i

= X× (

i

mi)Ẋ + X× (

i

miṙ

i) + (


i

mir

i)× Ẋ +


i

r′i ×miṙ

i

= X×P + X× (

i

miṙ

i) + (


i

mir

i)× Ẋ +


i

r′i × p

i

6

Let us look more closely at the two terms in parenthesis. Expanding

imir

i gives:∑

i

mir

i =


i

mi(ri −X) =

i

miri −

i

miX = (

i

mi)


imiri∑
imi

− (

i

mi)X = MX−MX = 0


imiṙ


i is simply the time derivative of this function, and hence is also 0. Thus, the total angular momentum

of a rigid body with respect to the origin is given by:

L = X×P +

i

r′i × p

i = L

point + Lspin

Conveniently, the angular momentum decomposes into a component associated with the center of mass and
the total momentum of the system, and a component associated with velocity relative to the center of mass.
We will now compute L̇point, and use this to obtain a simple formula for L̇spin.

L̇point = Ẋ×P + X× Ṗ = X× F

Here we have used that Ẋ is parallel to P and that Ṗ = F as proved in the previous section. We now
compute L̇spin as

L̇spin = L̇− L̇point =

i

ri × fexti −X× F

=

i

(X + r′i)× f
ext
i −X× F = X×


i

fexti +

i

r′i × f
ext
i −X× F

= X× F +

i

r′i × f
ext
i −X× F =


i

r′i × f
ext
i

This is a remarkable result! If we compute the angular momentum in the non-inertial (accelerating) reference
frame of the center of mass, the rate of change of the angular momentum takes the simple form L̇spin =∑

i r

i × f

ext
i .

We now restrict ourselves to 2D, and without loss of generality consider a body undergoing only a
rotational motion with instantaneous angular velocity ω (as we know how to compute torques in the non-
inertial frame of the center of mass). We can compute the velocity of the ith point as:

ṙ′i =

∣∣∣∣∣∣
x̂ ŷ ẑ
0 0 ω
xi yi 0

∣∣∣∣∣∣ = ω

−yixi

0


The angular momentum of this particular point is given by (in 2D a scalar) mir

i × ṙ


i. Expanding this cross

product we find that:

mir

i × ṙ


i = mi

∣∣∣∣∣∣
x̂ ŷ ẑ
xi yi 0
−yi xi 0

∣∣∣∣∣∣ω = mi

 00
x2i + y

2
i


ω = mir′i · r′iω

We can thus express the total angular momentum in the center of mass frame as

Lspin =

i

mir

i × ṙ


i =

(∑
i

mir

i · r

i

)
ω = Iω

where I =

imir

i · r

i is called the moment of inertia. In 2D, this value is a constant throughout the simu-

lation. To see that I is constant, rotate the entire body about the center of mass: I =

imi(Rr

i)

TRr′i =∑
imi(r


i)

TRTRr′i =

imi(r

i)

T r′i, as a rotation is an orthogonal matrix (R
T = R−1, and thus RTR = Id).

In conclusion, we find that L̇spin = Iω̇ =

i r

i × f

ext
i = Γ, which gives us a simple expression for ω̇.

Note that here we rely on the fact that I is constant in time (a fact not true in 3D).

7

3.4 Equations of Motion

To summarize, the equations of motion for our reduced coordinate representation of a rigid body are given
by:

d

dt



X
θ


ω


 =





ω

F/M
Γ/I




In two dimensions, X, Ẋ, and F are vectors of length 2. θ, ω, Γ, M , and I are scalars. From here, one
can utilize his or her favorite time integration technique.

4 Required Features for Theme III

4.1 Computation of Total Mass, Center of Mass, Moment of Inertia

Please edit RigidBodies/RigidBody.cpp and complete the computeTotalMass, computeCenterOfMass, and
computeMomentOfInertia methods. The interface for these methods is documented in the header file. The
oracle will check and grade these values for all scenes. If your values are incorrect, the oracle will print the
correct values.

4.2 Computation of Momentum

Please edit RigidBodies/RigidBody.cpp and complete the computeTotalMomentum method. The interface
for this method is documented in the header file. The oracle will check and grade this value for all scenes.
If your value is incorrect, the oracle will print the correct value.

4.3 Computation of Angular Momentum

In Section 3.3.2 we decomposed the total angular momentum of a rigid body into L = Lpoint +Lspin. Please
edit RigidBodies/RigidBody.cpp and complete the computeCenterOfMassAngularMomentum and compute-
SpinAngularMomentum methods. The interface for these methods is documented in the header file. The
oracle will check and grade these values for all scenes. If your values are incorrect, the oracle will print the
correct values.

4.4 Computation of Kinetic Energy Components

In Section 3.3.2, we separated the angular momentum into the components L = Lpoint + Lspring. Please
derive a similar decomposition for the kinetic energy of a rigid body; that is, separate the kinetic energy
into components T = T point + T spin. Please implement these functions in the computeCenterOfMassKinet-
icEnergy and computeSpinKineticEnergy functions of RigidBody.cpp. These values are serialized and graded.
If your values are incorrect, the oracle will print the correct values.

4.5 Explicit and Symplectic Euler

Please edit the stepScene methods of RigidBodies/RigidBodyExplicitEuler.cpp and
RigidBodies/RigidBodySymplecticEuler.cpp to implement explicit Euler and symplectic Euler respectively.
The documentation for these methods is included in the source files. After implementing these methods,
you should be able to pass all tests under theme3assets/inertiatests/. As in the first theme, your symplectic
Euler implementation should be implicit in velocity during the position update.

8

v 0

v 1

v 2

v 3n̂ i

Figure 4: Subdivision of an edge for computing the wind force. Here the number of sample points is four.

4.6 Near-Earth Gravity Force

Please implement a force corresponding to the (pointwise) potential U(ri) = −mig · ri. As this force acts on
all points of the rigid body, to simplify computations, work out the total potential energy, force, and torque
on the rigid body due to this force:

U body =

i

U(ri)

Fbody =

i

−∇U(ri)

Γbody =

i

r′i ×−∇U(ri)

The resulting formulas are very simple and intuitive.
Please fill in the implementations of computePotentialEnergy and computeForceAndTorque in Rigid-

Bodies/RigidBodyGravityForce.cpp (note that the potential energy is serialized and graded). After imple-
menting these methods, you should be able to pass all tests under theme3assets/gravitytests/.

4.7 Spring Force

Please implement a force corresponding to the potential U(ri, rj) =
1
2
k(|rj − ri| − l0)2, where ri and rj are

either points on a rigid body or fixed points in space. If the endpoint is a rigid body, you will have to compute
both the force and the corresponding torque. Please fill in the implementations of computePotentialEnergy
and computeForceAndTorque in RigidBodies/RigidBodySpringForce.cpp (note that the potential energy is
serialized and graded). The interfaces for these functions are documented in the corresponding header. After
implementing these methods, you should be able to pass all tests under theme3assets/springtests/.

Note: Please be careful with the special case where the spring length is 0 and the rest length l0 is 0 (can
you see where in the force computation this will cause trouble?). You will have to explicitly check for this
case and exert no force when it is detected.

4.8 Wind Force

We will implement a force that approximates the action of ‘wind’ on a rigid body. The magnitude of this
force should scale with the difference between the wind’s velocity and the velocity of some point on the edge
in the normal direction, that is with (vwind − vj) · n̂i. The velocity along an edge of a rigid body is not
constant, so we will have to subdivide the edge into at least two regions when computing the action of this
force. See Figure 4. To summarize, given a rigid body with N edges and P sample points per edge:

Fwind =

N−1∑
i=0

P−1∑
j=0

β
li
P

(vwind − vj) · n̂in̂i

9

β scales the magnitude of the force, and li/P is the fraction of length assigned to each sample point of the i
th

edge. As the vertices of the rigid body are specified in clockwise order, one can easily compute the outward
normal n̂i for a given edge.

When implementing this force, you will have to compute the torque at each sample point on each edge.
Please implement this force in the source file RigidBodies/RigidBodyWindForce.cpp.

5 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
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

5.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

6 Food For Thought

Answers to these questions are not required, we have included them only for your personal edification.

1. Did we have to choose the center of mass to encode the body’s translation? What would happen if we
selected a different point in body space?

2. Can you identify where in this derivation we used the assumption of rigidity?

7 Complications in moving to 3D

The complexity of a rigid body simulation jumps significantly from 2D to 3D, be warned! Complications
include:

1. The axis of rotation is now a three vector.

10

https://ffmpeg.org/ffmpeg.html
https://piazza.com/class/kta7z4zpa6l1ob?cid=189

2. The moment of inertia is dependent on orientation (it is a tensor).

3. The derivative of angular momentum has extra terms because the inertia tensor changes over time.

4. The representation of orientation is no longer a single angle. There are many options, all with tradeoffs:

(a) Euler angles: Simple (only three scalars), but suffers from ‘gimbal lock’ (who has watched Apollo
13?). There are also 24 possible conventions for representing orientation with Euler angles, and
this ambiguity can be confusing.

(b) Rotation matrix : Simple, but suffers from numeric drift, and requires the storage of nine scalars.

(c) Quaternion: Requires only four scalars, easy to ‘renormalize’ to correct for drift. Renormalization
is not free, however. Conceptually more difficult to work with.

(d) Exponential map: No renormalization required. Updating orientations fit ‘naturally’ into time
stepping routines such as explicit Euler.

11

Introduction
New XML Features
Rigid Body Dynamics
Body Space, World Space, and Center of Mass
Linear and Angular Velocity
Forces and Torques
Center of Mass Acceleration
Angular Acceleration

Equations of Motion

Required Features for Theme III
Computation of Total Mass, Center of Mass, Moment of Inertia
Computation of Momentum
Computation of Angular Momentum
Computation of Kinetic Energy Components
Explicit and Symplectic Euler
Near-Earth Gravity Force
Spring Force
Wind Force

Creative Scene
Creative scene submission

Food For Thought
Complications in moving to 3D