Assignment 6 - Intro. to Simulation and Animation

Due Wednesday Nov. 27, 2019 at 3:00 PM

Overview

This is a two part assignment that focuses on concepts from (1) simulation and (2) animation.

The first part on simulation has you derive variational time integrators for three physical systems: a single spring pendulum, a double spring pendulum, and an elastic mesh. You are given programs to simulate all three systems, but each program is missing the lines of code needed to perform the integration step. Your task in this first part is to complete these programs by implementing the missing update rules. Of course, to ensure a stable simulation, you need to make sure that these update rules are symplectic. Hence, you need to derive the update rules by using the discrete Lagrangian and applying the discrete Euler-Lagrange (DEL) equations.

The second part on animation has you write two programs that each examines a set of keyframes and interpolates across the set to generate appropriate frames between keyframes. The interpolation process involves the use of splines – Catmull-Rom splines to be exact.

Before you begin, you may want to review the assignment material with these lecture notes. The files for this assignment can be downloaded here. More in-depth details for the assignment will be given in the sections below. There is also an extra credit section at the bottom.

By the way, if you end up liking the material from this assignment a lot, then we highly recommend taking CS176 Introduction to Computer Graphics Research, taught by Prof. Mathieu Desbrun and Prof. Peter Schröder. The course does vary a bit each year, but it usually covers more advanced topics related to time integration and splines. The course is usually offered during the winter term.

Your Assignment: Part 1 (50 Points Total)

Within the Part1 folder are three more folders. Each folder contains a different simulation program that you need to complete:

1.

The Single Spring Pendulum (20 Points): The Single_Spring_Pendulum folder contains a nearly-complete program for simulating a 2D single spring pendulum. The program is missing the lines of code for the integration step – i.e. the update rules needed to transition the system from one time step to the next. Your task is to derive these update rules and implement them in the TODO section of the single_spring_pendulum_demo.cpp file. To ensure a stable simulation, you need to make sure that your update rules are symplectic; that is, you should derive them using the discrete Lagrangian and the discrete Euler-Lagrange (DEL) equations.

Regarding the single spring pendulum, imagine a simple pendulum of mass $m$ , except replace the string that attaches to the mass with a spring of stiffness $κ$ and rest length $l$ . As gravity pulls on the mass, the spring stretches and contracts into different lengths from its rest length. This generates a spring force that causes the pendulum system to behave more chaotically than a simple pendulum.

To help you start off, let us do the first few steps of the derivation together. First, we are going to use Cartesian coordinates for this derivation as opposed to polar coordinates. While both coordinate systems can work, Cartesian coordinates may prove to be more convenient when dealing with springs. A little introductory physics work gives us the kinetic energy of the system as:

$T = \frac{1}{2} m (\overset{\cdot}{x}^{2} + \overset{\cdot}{y}^{2})$

and the potential energy as:

$U = \frac{1}{2} κ {(\sqrt{x^{2} + y^{2}} - l)}^{2} - m g y$

We can then write our continuous Lagrangian as:

$L = T - U = \frac{1}{2} m (\overset{\cdot}{x}^{2} + \overset{\cdot}{y}^{2}) - \frac{1}{2} κ {(\sqrt{x^{2} + y^{2}} - l)}^{2} - (- m g y)$

Since we are working in Cartesian coordinates, we characterize the spring pendulum system with two position-state variables: $x$ and $y$ . Note that this is a bit different from the lecture notes’ simple pendulum example, where we only have the angle $θ$ as our only state variable. Formally, we represent our position-state variables with the vector $q = [x, y]$ . Let us see how this works out.

We need to write the discrete analog to our continuous Lagrangian. This process essentially boils down introducing time steps and replacing the velocity terms, $\overset{\cdot}{x}$ and $\overset{\cdot}{y}$ , with their finite differences. Letting our time step be $Δ t$ , our state at the current time step be $q_{k}$ , and our state at the next time step be $q_{k + 1}$ , we can express our discrete Lagrangian as:

$L_{d} (q_{k}, q_{k + 1}, Δ t) = Δ t [\frac{1}{2} m ({(\frac{x_{k + 1} - x_{k}}{Δ t})}^{2} + {(\frac{y_{k + 1} - y_{k}}{Δ t})}^{2}) - \frac{1}{2} κ {(\sqrt{x_{k}^{2} + y_{k}^{2}} - l)}^{2} + m g y_{k}]$

Take a second to verify this yourself. Refer back to the lecture notes for a refresher of the discrete Lagrangian if necessary.

We now need to compute the DEL equations with respect to $q_{k}$ and to $q_{k + 1}$ . But since $q_{k} = [x_{k}, y_{k}]$ is a vector of two position-state variables, we need to compute equations for all of: $p_{x_{k}}$ , $p_{x_{k + 1}}$ , $p_{y_{k}}$ , and $p_{y_{k + 1}}$ :

$p_{x_{k}} = - D_{1} L_{d} (x_{k}, x_{k + 1})$

$p_{y_{k}} = - D_{1} L_{d} (y_{k}, y_{k + 1})$

$p_{x_{k + 1}} = D_{2} L_{d} (x_{k}, x_{k + 1})$

$p_{y_{k + 1}} = D_{2} L_{d} (y_{k}, y_{k + 1})$

At this point, all that is left is to take the correct derivatives and rearrange terms with some algebra. We let you take it from here.

We show further below a movie of how the single spring pendulum program should perform once the correct symplectic update rules are coded. The command line arguments for the program are:

./single_pendulum [xres] [yres] [x_start] [y_start]

where the x_start and y_start variables are the starting positions of the pendulum mass. The arguments we use for the movie are:

./single_pendulum 500 500 4 -2

The top-left corner of the screen displays the kinetic and potential energies at each time step. The top-right corner displays the total energy. The bottom-left corner displays the smallest and largest total energies the program has seen throughout the course of the simulation. Recall from the lecture notes that variational time integrators are supposed to follow the law of conservation of energy. Due to numerical imprecision, we cannot perfectly conserve energy in any discrete simulation; but we can get close by keeping the energy bounded. If you were to let the program run for an extended period of time, then you will see that the “Min Total” and “Max Total” numbers barely change throughout the entire simulation.

2.

The Double Spring Pendulum (20 Points): The Double_Spring_Pendulum folder contains a nearly-complete program for simulating a 2D double spring pendulum. Like the single spring pendulum program, the double spring pendulum program is also missing the lines of code for the integration step. Your task is to derive and code the symplectic update rules for this system in the TODO section of the double_spring_pendulums_demo.cpp file.

The double spring pendulum system is similar to the single spring pendulum system, except a second pendulum mass is attached to the first mass via another spring. This results in a system that acts quite chaotically, as shown in the left video at the top of this page. Just like with the single spring pendulum system, start this derivation by writing down the continuous Lagrangian using Cartesian coordinates. Then write down the discrete analog and apply the DEL equations. The derivatives and algebra get quite nasty here, so we recommend using symbolic mathematics software like Wolfram Mathematica to aid your derivation.

In the end, you should have 8 update rules – 4 for each pendulum mass. You should have an update rule for each of these variables:

m1.x
m1.y
m2.x
m2.y
m1.px
m1.py
m2.px
m2.py

The left movie at the top of this page shows the double spring pendulum program running when the symplectic update rules are coded correctly. The program takes the following command line arguments:

./double_pendulum [xres] [yres] [x_start_1] [y_start_1] [x_start_2] [y_start_2]

and the arguments that we use for the movie are:

./double_pendulum 500 500 0 -1 4 2

3.

The Elasticity Demo (10 Points): The Elasticity folder contains a nearly-complete program for simulating a 2D elastic mesh. Like the pendulum programs, the elasticity demo is also missing the lines of code for the integration step. Your task is to code the symplectic update rules for this sytem in the TODO section of the elastic_demo.cpp file.

Unlike the pendulum exercises, this part does not require you to do a full derivation. This is more of a test of how comfortable you are with the general concept of the Lagrangian and its discrete counterpart. The actual mathematics needed to simulate an elastic surface or volume is unfortunately outside the scope of this course (check out the extra credit section though if you’re curious).

The elastic system is governed by general Lagrangian mechanics, meaning we can write a continuous Lagrangian of the form:

$L = T - U = \frac{1}{2} m (\overset{\cdot}{x}^{2} + \overset{\cdot}{y}^{2}) - U (x, y)$

for some potential function $U$ that depends on $x$ and $y$ . A pair, $(x, y)$ , represents the coordinates of a vertex on our elastic mesh. You do not need to worry about the exact details of $U$ . For this exercise, assume that you know how to compute:

$\frac{\partial U}{\partial x} = - f_{x}$

$\frac{\partial U}{\partial y} = - f_{y}$

where $f_{x}$ and $f_{y}$ are the forces acting on vertex $(x, y)$ . The given elasticity demo program already computes $f_{x}$ and $f_{y}$ for every vertex for you (see the comments in the TODO section). Knowing that you already have values for $f_{x}$ and $f_{y}$ for every vertex, can you do a quick derivation for the update rules of this system?

As always, start by writing down the discrete analog to the continuous Lagrangian and then apply the DEL equations. If you begin to wonder how you’re supposed to handle $U$ , remember that you already “know” $f_{x}$ and $f_{y}$ .

The elasticity demo program just takes in one argument: an .obj file of a 2D mesh. We have provided one called man.obj from [1]. You can run the program with:

./simulate man.obj

The program allows you to click and drag the vertices of the mesh. If you move the vertices a little, then you will see the mesh start to “jiggle” due to the elastic interactions among the vertices. Those of you who wish to learn about the actual process for simulating elasticity should refer to the extra credit section. The following movie gives a visual demonstration of the program at work:

As a note, the vertices are not meant to be dragged a large distance. If you drag the vertices outside the bounds of a face, then you will cause triangles to overlap with each other. All the force calculations afterward will go crazy and cause the mesh to explode. The program unfortunately has no protection against this, but you may (for extra credit) want to try adding something to prevent the vertices from being dragged too much.

Your Assignment: Part 2 (50 Points Total)

Within the Part2 folder are two more folders. These folders contain test files for the two programs that you will write:

1.

The “I-Bar” Animator (30 Points): In the I_Bar folder, there is a file named I-bar_code.cpp and three other files named test.script, test2.script, and test3.script. The .script files specify the keyframes for a single point. The format of the file is as so (using test.script as an example):

75                         # total number of frames in the animation
Frame 0                    # first keyframe is frame 0 in the animation
translation -5 0 0         # translation of the point in the scene for this keyframe
scale       1 1 1          # scaling of the point
rotation    0 0 1 0        # rotation of the point, where the angle is the last numerical value in degrees
Frame 15                   # second keyframe is frame 15 in the animation
translation -1 24 0
scale       2 2 2
rotation    0 1 1 120
...                        # and so on...

Your assignment is to write a program that processes a given script input file and uses OpenGL to create an animation loop based off the keyframes given in the script. Here are the specifications for the program:

Call the program keyframe.
Use the following specifications for the glFrustum camera parameters (you may hardcode these):
- near: 1.0
- far: 60.0
- left: -1.0
- right: 1.0
- bottom: -1.0
- top: 1.0
Use [0, 0, 40] for the camera position. Let the camera angle be 0 (i.e. no need for a call to glRotatef for the camera this time).
The program should initially just display the “I-bar” object drawn by the code snippet in I-bar_code.cpp. The center of the I-bar should be positioned at the point described by frame 0. Without any transformation specifications, the given code snippet in I-bar_code.cpp draws the I-bar centered at the origin.
When the user presses a key on the keyboard (pick any key), the program should step forward one frame. The point described by the frames should always be the center of the I-bar.
Use Catmull-Rom splines to interpolate the translations, scalings, and rotations across keyframes. Remember that you must represent rotations as quaternions, and interpolate the quaternions component-wise. Don’t forget to normalize the quaternions after interpolating them!
The last keyframe loops back to frame 0.

You may process the script files however you like, though we recommend simply using the ifstream functions. You might also find it convenient to use Eigen for the matrix and vector operations.

Your code from Assignment 3 and the code from Part1 should be good references for you to write this program. It does not have to be too long.

The following movies show the animation loops that you should get for each of the three test scripts. The movies are for test.script, test2.script, and test3.script in that order. You do not need to display the frame number in your program. It exists in these videos to show you how each of your frames should look like.

2.

Interpolating a Smoothing Bunny (20 Points): In the Bunny_Frames folder, there are is a folder called keyframes with 5 .obj files, each one representing a keyframe of the Stanford bunny as it undergoes implicit smoothing. The keyframes are of frames 0, 5, 10, 15, and 20 of the smoothing process. Your task is to write a program that interpolates across these keyframes to generate the missing 16 frames.

The order in which the vertices are listed is consistent across the 5 given .obj files. You can treat each vertex position as a translation and simply use Catmull-Rom splines to interpolate a vertex across the keyframes. For the earlier frames (before keyframe 5), you should use keyframe 0 for both the first and second control points of your splines. For the final few frames (after keyframe 15), you should use keyframe 20 for both the third and fourth control points of your splines.

You may take your liberties on how you write this program. For instance, you may choose not to write an OpenGL program and instead write a regular program that takes in the 5 .obj keyframe files and outputs 16 .obj files for the missing frames. We have provided the correct, interpolated frames, in a folder called interpolated_frames. You can use these to check the output of your program (e.g. compare each file in this folder with each file that your program outputs and see if they match to within floating point precision). You can also load each .obj file output of your program into your Assignment 5 program to visually check each frame. The movie on the right, at the top of this page, shows the correct bunny render for each of the 20 frames. And of course, if you really want to, you can always write an OpenGL program that interpolates all the frames and steps through them as shown in the movie.

Extra Credit

Here are a few extra credit ideas related to simulation and animation:

Create a program that simulates the double spring pendulum in 3D. Design an appropriate way for a user to see the behavior of the system as time goes on.
Modify the elasticity demo to incorporate gravity and a surface (e.g. floor) below the elastic mesh. Correctly model the behavior of the elastic mesh as it falls due to gravity onto the floor and bounces around. Adding friction is another idea.
For those of you who are curious about simulating elastic surfaces and volumes, we direct you to the following papers: [2] and [3]. The first resource is a paper that introduces and tests more advanced variational time integrators using elastic meshes. The second resource is a set of notes that detail the discretization of elastic forces. Look through these and try to implement an idea from them. For instance, try modifying the elasticity demo to load in a different mesh and have it still display correct energy conservation. This can be a tricky, since all the meshes that we’ve been dealing with are 3D and don’t easily convert over to 2D.
Modify the elasticity demo to prevent the user from dragging a vertex in the mesh beyond the bounds of a surrounding face (so that we can’t make the mesh explode by dragging the vertices too much). You will probably want to use a halfedge implementation for this to find the distances from the dragged vertex to the opposite edges of incident faces. You might then pick a distance to be the upper bound of how much the vertex can travel based on which face the user moves the vertex towards.

Of course, you can always pursue your own idea as well. Please consult with the TA’s ahead of time to make sure your extra credit topic is appropriate. Please also include details regarding your extra credit in the README. Extra credit points will be awarded based on the complexity of the work.

TA’s are not liable to help with extra credit.

What to Submit

Before submitting, please comment your code clearly and appropriately, making sure to give details regarding any non-trivial parts of your submission.

Submit a .zip or .tar.gz file containing the files that we need for all parts of the assignment. In addition, please submit in the .zip or .tar.gz a README with instructions on how to compile and execute your programs as well as any comments and information that you would like us to know. The README is especially important if your submission has significant issues, since we want to know what you have tried to debug; this will help us figure out how close you were at getting the programs to work and give you a grade faster.

We ask that you name your .zip or .tar.gz file lastname_firstname_hw6.zip or lastname_firstname_hw6.tar.gz respectively, where you replace the “firstname” and “lastname” fields with your actual first and last name.

Please submit your zip or tar.gz to moodle in the area for Assignment 6.

References

[1] http://www.dgp.toronto.edu/~rms/software/Deform2D/

[2] http://www.geometry.caltech.edu/pubs/KWTKMSD06.pdf

[3] http://www.femdefo.org/

Written by Kevin (Kevli) Li (Class of 2016).
Links: Home Assignments Contacts Policies Resources