COMS W4167: Physically Based Computer Animation
Theme IV, Milestone I – Elastic Simulation
Due 22:00:00 PM Tuesday 23th 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 IV, we’ll tackle the elastic body problem. Simulation of the elastic material is based on formulation
of an elastic energy that is quadratic in a strain that describes how the object is deformed. Using different
aspects of the deformation as our strain, we can capture different modes of deformation in the energies, and
thus obtain the forces that counter-act these modes. In this milestone we’ll examine three elastic forces:
spring force, bending force, and the constant-strain-triangle (CST) force.
2 New XML Features
This milestone introduces the following new simulation features:
1. The elasticbodyspringforce feature specifies a spring force:
The i1 and i2 attributes specify between which particles this force acts; note that different than the
spring force we implemented in the first theme, this force does not require an edge between the two
endpoints. The alpha attribute specifies the material’s elastic modulus (equivalent to EA where E is
the young’s modulus and A is the area of the cross-section of the spring). The l0 attribute specifies
the rest length which means the same as the spring in theme I.
2. The elasticbodybendingforce feature specifies a bending force:
The i1, i2 and i3 attributes specify between which particle this force acts. To represent a bending
deformation we need three particles, where the second one is the hinge. The alpha attribute specifies
the material’s elastic modulus in resisting the bending mode. The theta0 attribute specifies the rest
angle in radian.
3. The elasticbodycstforce feature specifies a CST force:
1
Figure 1: The deformation map.
The i1, i2 and i3 attributes specify the three vertices of the triangle where this force acts. The
youngsmodulus and poisonratio attributes are the material parameters. ebxx, ebxy and ebyy encodes
the “resting strain” �̄ that is discussed in the next section. xb1x and xb1y specify the undeformed
position of vertex 1, and xb2x, xb2y, xb3x and xb3y are similar.
3 Kinematics
3.1 The Deformation Map and Deformation Gradient
The kinematics of a deforming material are represented by a deformation mapping Φ : R2 → R2 taking every
‘undeformed’ material point x on the plane to its corresponding deformed position Φ(x). In general, Φ is
an arbitrary nonlinear map taking 2 coordinates denoting a position to a new set of 2 coordinates denoting
some other position.
Now we’d like to know how much this deformation is stretching, compressing or shearing the material.
Since a generic deformation is an arbitrary nonlinear map, we only look at each point locally (by examining
its infinitesimal neighborhood). Consider a point x̄ in the undeformed space, and a nearby point x̄ + δx̄.
The deformed position are x and x + δx, and by definition of Φ we have:
x = Φ(x̄)
x + δx = Φ(x̄ + δx̄)
= Φ(x̄) +∇Φ(x̄)δx̄ +O(‖δx̄‖2)
δx = ∇Φ(x̄)δx̄ +O(‖δx̄‖2)
∇Φ(x̄) is the Jacobian of the transformation Φ evaluated at x̄, often called the Deformation Gradient.
The high order terms in the end can be omitted because we assume δx̄ is infinitesimal.
3.2 The Strain
Now we are ready to evaluate how much lengths change after such a transformation. Under the assumption
that δx̄ is infinitesimal, the line from x̄ to x̄+δx̄ in the undeformed space remains a line after the deformation,
2
going from x to x+ δx. Therefore the lengths of this piece of material before and after deformation are ‖δx̄‖
and ‖δx‖ respectively. For the ease of computation we look at the change in squared length (so as to avoid
taking square roots):
δxT δx− δx̄T δx̄ = (∇Φ(x̄)δx̄)T (∇Φ(x̄)δx̄)− δx̄T δx̄
= δx̄T∇Φ(x̄)T∇Φ(x̄)δx̄− δx̄T δx̄
= δx̄T [∇Φ(x̄)T∇Φ(x̄)− I]δx̄
The reason we can use change in squared length instead of change in length itself as a measure of stretching
is that, to the first order, they are equivalent as long as the deformation is small. In order to make this
change in squared length a scale-invariant representation of relative stretching/compression, we extract out
the term in the middle, [∇Φ(x̄)T∇Φ(x̄)− I], and call it the strain, or �.
This strain is a second-order tensor, unlike the 1D case where a single scalar fully describes how much a
spring is stretched as shown in class two weeks ago. In 2D materials can be stretched differently in different
directions, so a second-order tensor is needed. The physical meaning of this tensor is that when you hit
it on both sides with a direction δx̄, it returns a change in squared length in that direction. In its matrix
representation, the eigenvectors of this matrix are the two directions that are stretched or compressed the
most, and the corresponding eigenvalues tell you how much the stretching/compression is in that direction.
Another interesting thing to note is that � is rotation-invariant, even though ∇Φ is not. As a measure
of deformation, the strain should not depend on any rigid motion such as translation and rotation. It is
obvious that ∇Φ, which is the result of a spatial differential operator, stays the same under translation, but
what about rotation? To test this, we apply a constant rotation Q to any deformation Φ to obtain a new
deformation Φ̃, and see how � responds:
Φ̃ = QΦ
∇Φ̃ = Q∇Φ
�̃ = ∇Φ̃T∇Φ̃− I
= (Q∇Φ)T (Q∇Φ)− I
= ∇ΦTQTQ∇Φ− I
= ∇ΦT∇Φ− I
= �
Thus we’ve convinced ourselves that our strain � is invariant under any rigid motions, therefore it qualifies
as a measure of deformation.
3.3 Example Deformation
Let’s look at an example with real numbers. Consideration an affine deformation of a rectangular sheet of
rubber as in Figure 2. This rectangle is stretched along the green axis by 80%, and compressed along the
red axis by 40%. This means our deformation map Φ is multiplication by a matrix with two eigenvalues: 1.8
corresponding to the eigenvector (3, 1) (green axis), and 0.6 corresponding to the eigenvector (-1, 3) (red
axis). Therefore Φ looks like this:
Φ(x) =
1
√
10
(
3 −1
1 3
)(
1.8 0
0 0.6
)
1
√
10
(
3 1
−1 3
)
x
=
(
1.68 0.36
0.36 0.72
)
x
3
Figure 2: An example deformation of square, with stretching in one direction and compression in the other.
The two arrows denote the two eigenvectors of the strain tensor with corresponding eigenvalues (λ1 and λ2)
marked on the arrows.
We can compute the strain tensor � from Φ:
� = ∇Φ(x̄)T∇Φ(x̄)− I
=
(
1.68 0.36
0.36 0.72
)T (
1.68 0.36
0.36 0.72
)
− I
=
(
1.952 0.864
0.864 0.352
)
It is easy to verify that this 2×2 matrix has eigenvalues 2.24 and -0.64, corresponding to a 80% stretching
(1.82 − 1 = 2.24) and a 40% compression (0.62 − 1 = −0.64).
4 Physics
The strain tensor gives us a complete description of the deformation of the material, but it doesn’t tell us
how the material will respond. That is out of the scope of kinematics, and we have to put in some ingredients
from physics. Here in this section we’ll attempt to write down the energy that arises from the deformation,
which governs how the simulation will take place in time. For this purpose, we can proceed in either a
mechanics-oriented way, or a mathematics-oriented way.
4.1 The Mathematical Way
Imagine we have somehow obtained the energy function w of strain �, we can take its Taylor expansion:
w(�) = w(0) +
∂w(0)
∂�
�+
1
2
∂2w(0)
∂�2
�2 +O(‖�‖3)
w(0) is a constant independent of � and can be conveniently assume to be 0. The linear term is also
0 because the energy should be minimized at zero strain (the definition of “resting shape”). If we adopt
the small strain assumption, the higher order terms can be omitted since the energy is dominated by the
quadratic term. This is basically saying, no matter how complicated the material laws are, under the small
strain assumption we can always write down the energy as a quadratic function of the strain.
4.2 The Mechanical Way
Physically, the energy is equal to the product between the stress and the strain:
4
w(�) =
1
2
2∑
i=1
2∑
j=1
�ijσij
and the stress σ is defined as the elastic moduli times the strain:
σkl =
2∑
i=1
2∑
j=1
cijkl�ij
therefore
w(�) =
1
2
2∑
i=1
2∑
j=1
2∑
k=1
2∑
l=1
cijkl�ij�kl
Here we obtain an energy that is quadratic in strain again. c is a fourth-order tensor that encodes all
the material properties such as stiffness upon stretching/compression in each direction and shearing etc.
Although a generic fourth-order tensor requires a 16 scalars to represent, symmetry in both the strain tensor
and the stress tensor reduces this number to 9. If we further assumes that the material is isotropic, then we
only need two parameters to fully describe this material: Young’s modulus E and Poisson ratio ν:
σxxσyy
σxy
= E
1− ν2
1 ν 0
ν 1 0
0 0
1− ν
2
�xx�yy
2�xy
4.3 Total Energy
The w(�) defined above gives the energy density at a point in the material. Since � is a function of undeformed
space position x̄, so is w. To obtain the total elastic energy that is stored in an elastic body, we integrate
this energy density over the entire body:
W =
∫
Ω
w(x̄)dΩ
For this milestone, we simply assume that the strain is constant for any element (discussed in the next
section), therefore the energy is simplified to W = w(x̄)A where A is the 2D element’s area.
4.4 Discretization
All the analysis up to now has been performed on a continuum, in the form of differential equations. Continua
have an infinite number of degrees of freedom and cannot be manipulated directly in computers – we need
to discretize.
The discretization here is in the spatial dimensions, but it’s in the same spirit with the temporal dis-
cretization we’ve implemented in the first Theme with integrators. By spatial discretization we represent the
body of interest with a set of nodes, connected by edges, polygon faces or polyhedron volumes, depending
on the dimension of the body. These building blocks are called elements. There are a fixed number of nodes
in each element, and due to the limited amount of information these nodes can provide, we have to assume
5
certain properties on the element itself. For example, a triangle element with only one node in each corner
(one of the simplest and most commonly used element types) simply does not have any information about
how the deformation mapping changes over the triangle to the second or higher order, so we have to assume
that the deformation Φ is linear, and when we need information about a point inside the triangle, we linearly
interpolate between the three nodes. Linear deformation map implies constant deformation gradient, and as
a result this element type is called constant strain triangle (discussed later in more details).
Of course the interpolation won’t be able to faithfully reproduce the true property values at an arbitrary
point in the element, so a discretization error is introduced; but since our hardware technology isn’t advanced
enough to handle infinite degrees of freedom, this error is inevitable. As long as the discretization error
approaches zero as we refine our mesh further and further, we say the discretization converges and thus is
valid. For our milestone, we do not consider more complicated (but potentially more accurate) elements.
5 Required Features
5.1 Elastic Body Spring Force
The spring force is similar to the spring force we implemented in Theme 1, with the only difference being
that we use a shape-independent material property as the stiffness. Following the concept from above, we
define the strain as:
� =
‖xi − xj‖ − l0
l0
The energy density is:
w =
1
2
α�2
=
α
2l20
(‖xi − xj‖ − l0)2
Assuming that the strain is constant along the spring, the total energy is
W = wl0
=
α
2l0
(‖xi − xj‖ − l0)2
Taking the gradient of this energy gives the force. Please implement the computation for energy, force, and
hessian in addEnergyToTotal(), addGradEToTotal() and addHessEToTotal() in ElasticBodySpringForce.cpp.
Note that implementation of the energy is optional; it won’t be graded. Implementing the energy will help
implementing the force by providing more intuitions and more debugging information.
5.2 Elastic Body Bending Force
When a thin rod is bent, why does it try to restore its initial straight configuration? It is because the
bending introduces differential compression/stretching along the cross section of the rod, and thus increases
the elastic energy. If we explicitly model the cross section of the rod, we will have an array of axial springs
side-by-side, connected by short and highly stiff radial-directed springs at endpoints. When bent springs on
the inner layer are compressed and springs on the outer layer are stretched, therefore creating the restoring
force. However this modeling scheme has two serious disadvantages: it is too stiff and has bad aspect ratios,
both of which lead to numerical difficulties. In practice we model thin rods as a single polyline, augmented
with a bending force to resist bending.
Since the bending energy arises from axial compression and stretching of the material, it is naturally
quadratic in the curvature κ of the material centerline:
6
Figure 3: The bending at the common node between two edges
w =
α
2
κ2
where α captures all other constants, including the material’s elasticity modulus and cross section area. As
our discretization refines, curvature κ can be approximated by
θ
‖ēij‖+ ‖ējk‖
, where θ is the angle by which
the next edge (eij) “deviates” from the direction of the previous edge (ejk):
θ = atan2(eij × ejk, eij · ejk)
where function atan2(y, x), a C library function, returns the direction angle of vector 〈x, y〉.
Therefore the total bending energy is:
W = w(‖ēij‖+ ‖ējk‖)
=
α
2
κ2(‖ēij‖+ ‖ējk‖)
≈
α
2
(
θ
‖ēij‖+ ‖ējk‖
)2(‖ēij‖+ ‖ējk‖)
=
α
2(‖ēij‖+ ‖ējk‖)
θ2
Allowing non-straight (non-zero curvature or non-zero bending angle) resting shape, we define the energy
as
W =
α
2(‖ēij‖+ ‖ējk‖)
(θ − θ0)2
Please calculate the gradient and hessian, and implement addEnergyToTotal(),anddGradEToTotal() and
addHessEToTotal() in ElasticBodyBendingForce.cpp. Note that implementation of the energy is optional; it
won’t be graded. Implementing the energy will help implementing the force by providing more intuitions
and more debugging information. The following formula may be useful:
7
Figure 4: A triangle undergoing a deformation Φ
∂arctan(x)
∂x
=
1
1 + x2
(
∂a× b
∂a
)T =
(
0 1
−1 0
)
b
5.3 Elastic Body Constant Strain Triangle Force
Both the spring force and the bending force essentially assumes an infinitely thin rod-like structure, where
the strain is fully captured by a single scalar (either the edge length change or the bending angle). For 2D
elastic bodies we need to use the full-fledged treatment in the previous sections. Now consider the triangle
involving particle i, j and k:
We assume that the strain is constant over the interior of the triangle (hence the name “constant strain
triangle”), and so the deformation map Φ is affine (without loss of generality, from here on we assume i =
1, j = 2 and k = 3):
x = ∇Φ(x̄− x̄1) + x1
∇Φ =
(
x2 − x1 x3 − x1
y2 − y1 y3 − y1
)(
x̄2 − x̄1 x̄3 − x̄1
ȳ2 − ȳ1 ȳ3 − ȳ1
)−1
where xi denotes the x component of node xi, yi denotes its y component, and the overhead bar denotes
the undeformed configuration as before. Denoting the components in ∇Φ as ∇Φ =
(
a c
b d
)
, we define the
strain:
8
� =
(
�xx �xy
�xy �yy
)
= ∇ΦT∇Φ− I
=
(
a2 + b2 − 1 ac+ bd
ac+ bd c2 + d2 − 1
)
and the energy:
W =
Ā
2
σxxσyy
σxy
T
�xx�yy
2�xy
=
Ā
2
(
E
1− ν2
1 ν 0
ν 1 0
0 0
1− ν
2
�xx�yy
2�xy
)T
�xx�yy
2�xy
=
EĀ
2(1− ν2)
[�2xx + 2ν�xx�yy + �
2
yy + 2(1− ν)�
2
xy]
where Ā is the undeformed triangle area:
Ā =
1
2
‖(x̄3 − x̄1)× (x̄2 − x̄1)‖
In order to find the force, we take the gradient of the energy above:
−F = ∇xW
=
EĀ
1− ν2
[(�xx + ν�yy)∇x�xx + (�yy + ν�xx)∇x�yy + 2(1− ν)�xy∇x�xy]
∇x�xx = 2a∇xa+ 2b∇xb
∇x�yy = 2c∇xc+ 2d∇xd
∇x�xy = a∇xc+ c∇xa+ b∇xd+ d∇xb
Now letting
(
x̄2x − x̄1x x̄3x − x̄1x
x̄2y − x̄1y x̄3y − x̄1y
)−1
=
(
e g
f h
)
we have
∇Φ =
(
(−e− f)x1 + ex2 + fx3 (−g − h)x1 + gx2 + hx3
(−e− f)y1 + ey2 + fy3 (−g − h)y1 + gy2 + hy3
)
∇a =
−e− f
0
e
0
f
0
,∇b =
0
−e− f
0
e
0
f
,∇c =
−g − h
0
g
0
h
0
,∇d =
0
−g − h
0
g
0
h
,
9
With these derivation, you should be able to implement the elastic energy, force, and hessian of a
constant strain triangle. Please do so in addEnergyToTotal(), addGradEToTotal(), and addHessEToTotal()
in ElasticBodyCSTForce.cpp. Note that implementation of the energy is optional; it won’t be graded.
Implementing the energy will help implementing the force by providing more intuitions and more debugging
information.
Debugging help
If you have problems implementing the elasticity, try debugging a few simple cases where you know the exact
solution. Then test the output of your algorithm for the deformation of each point to the exact solution that
you know to be true.
An easy example to use for this is the triangle with vertices (0, 0), (0, 1) and (1, 0) (ref. the first image
below). Another simple test case is the square with vertices (0, 0), (0, 1), (0, 1) and (1, 1) (ref. the second
image below)
A triangle with vertices (0, 0), (0, 1) and (1, 0)
A square with vertices (0, 0), (0, 1), (0, 1) and (1, 1)
10
6 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
6.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) Nov. 24 Please refer to announcement on Courseworks
11
https://ffmpeg.org/ffmpeg.html
https://piazza.com/class/kta7z4zpa6l1ob?cid=189
Introduction
New XML Features
Kinematics
The Deformation Map and Deformation Gradient
The Strain
Example Deformation
Physics
The Mathematical Way
The Mechanical Way
Total Energy
Discretization
Required Features
Elastic Body Spring Force
Elastic Body Bending Force
Elastic Body Constant Strain Triangle Force
Creative Scene
Creative scene submission