Programming Assignment 4: Physical Simulation

Description

Physical simulation is used in movies and video games to animate a variety of phenomena: explosions, car crashes, water, cloth, and so on. Such animations are very di cult to keyframe, but relatively easy to simulate given the physical laws that govern their motion. In this assignment, you will build simulations of increasing complexity, ending up with a cloth model built from a large number of particles connected to each other with springs.

Requirements (maximum 10p) on top of which you can do extra credit

1. Euler integrator (1p)

2. Spring system (2p)

3. Trapezoid integrator (2p)

4. Pendulum system (2p)

5. Cloth system (3p)

• Getting Started

In this assignment, you’ll be both solving ordinary di erential equations (ODEs) numer-ically by writing di erent kinds of integrators, and formulate dynamical systems of your own (springs, cloth, etc.) that are driven by forces such as gravity, springs, and viscosity.

1.1 Simple Newtonian Mechanics and Di erential Equations

The motion of many physically-based systems is governed by Newton’s laws of motion, which are second order di erential equations that, as you know, relate forces to acceler-ations. For a single point-like particle at position x, this means

m	d2x			= F;	(1)
		2
		dt

i.e., the familiar \force equals mass times acceleration”. The force F is a vector sum over all the di erent forces acting upon the particle. The forces can vary with time, and also other particles may have an e ect. We’ll talk more about forces below.

As we saw in class, it pays o to reduce higher-order di erential equations to rst order, so that we can write a solver once and use it for many di erent kinds of problems. For our simple point dynamics, this means adding an auxiliary velocity variable v to each particle, and writing

dt		v	=	F=m	:	(2)
d		x		v

This equation is of rst order | i.e., it only involves the rst time derivatives of the system state that consists of the position and velocity of the particle | and it says the same thing as the one above. Solving this di erential equations means that starting from a given initial state (x₀; v₀), we trace out the curve determined by the in nitesimal displacements

• v

d _v = dt _F=m

that say \in the in nitely small time step dt, my current position changes by dt times my current velocity, and my velocity changes by dt times the current forces divided by my mass”.

When we have a larger number n > 1 of point-like particles in our system, we write x_i for the position and v_i for the velocity of the ith particle. Because the above equations of motion apply to all particles separately (the acceleration for a single particle is determined by the forces acting on that particle), we can stack all the positions and velocities in a single big vector X(t)

• 1

x₁(t)

Bv₁(t)C

• C

Bx₂(t)C

X(t) = ^B ._. ^C ; (3)

B C

B . C

• C @x_n(t)A

v_n(t)

and following the above, write this bigger particle system abstractly as

• 1

v₁

• F₁(X; t)=m₁ C

dtX = B		v...2	C = f(X; t):		(4)
d			C	def
	B
	B	Fn(X; t)=mn	C
	B		C
@			A

The vector X has length 6n for n particles (all positions are 3D, as are all velocities), and it lives in what is called phase space. We denote the force acting on the ith particle by | you guessed it | F_i(X; t). The derivative d=dtX depends potentially on time | for example, a changing wind force | and the current state | e.g. spring forces are determined by the current positions of the particles | so we explicitly write it out as a function f(X; t) that returns the derivative vector given the state vector X and time t.

1.2 Forces

The core component of particle system simulations are forces. Suppose we are given a particle’s position x_i, velocity v_i, and mass m_i. We can then express forces such as

gravity

F(x_i; v_i; m_i) = m_ig; where g = 0 g 0 ^T

or perhaps viscous drag that opposes motion (given a drag constant k):

F(x_i; v_i; m_i) = kv_i:

We can also express forces that involve other particles as well. For instance, if we con-nected particles i and j with an with an undamped spring of rest length r and spring constant k, it would yield a force of:

F(xi; vi; mi) = k(jjdjj	r)	d		, where d = xi		xj :
		jjdjj

Think back to your high school physics: what is the expression of the gravitational force between two massive bodies that are so far away from each other that they can be considered point-like?

Summing over all forces yields the net force, and dividing the net force by the mass gives the acceleration a_i = _d^d_t v_i.

• Numerical integrators

It is important for you to understand the abstraction between numerical integrators and particle systems. The numerical integrator does not know anything about the physics of the system. It can request the particle system to compute the derivative f(X; t) using

the system’s evalF method. This function is the critical communication channel between a system and an integrator. It takes as input a state vector and returns a force vector for this particular input state, which are both represented as 1D arrays regardless of the precise type of particle system (only the size of the array varies). This allows integrators to be general and reusable.

A particle system stores its current state X, but the numerical integrator might request the forces for a di erent state, in particular for the trapezoidal method. It is critical that the particle system uses the correct state, the one requested by the call to evalF, to compute the forces. Make sure you understand the di erence between this internal state of the system and that requested by the integrator.

2.1 Refresher on Euler and Trapezoidal Rule

The simplest integrator is the explicit Euler method. For an Euler step, given state X, we examine f(X; t) at X, then step to the new state value. This requires to pick a step size h, and we take the following step based on f(X; t), which depends on our system.

X(t + h) =X + hf(X; t)

Look back to Equation (4) and see if you can spot a similarity; we’ve sort of replaced dt with a small but nite number and multiplied it over to the other side. This technique, while easy to implement, can be unstable for all but the simplest particle systems. As a result, one must use small step sizes (h) to achieve reasonable results.

There are numerous other methods that provide greater accuracy and stability. For this problem set, you will implement the Trapezoidal approach, which works by using the av-erage of the force f₀ at the current state and the force f₁ after an Euler step of step size h:

f₀ =f(X; t)

f₁ =f(X + hf₀; t + h)

h

X(t + h) =X + ₂ (f₀ + f₁)

This method makes it critical for the particle system to be able to evaluate the forces at a state X + hf₀ other than its current state X. The state of the particle system should not be updated until reaching the last equation.

2.2 Simple example system

We have supplied you with a rst-order ODE similar to one we saw in class (SimpleSystem).

This system has a single particle and its state is de ned by its x-y-z coordinates:

• 1

x

X = @ y A

z

and the right-hand side is

0 1

y

f(X; t) = @ x A :

0

This is only a rst-order system and the right-hand side of the ODE does not describe physical forces. (Using pen and paper, take a moment to plot some points and their associated derivatives on paper to understand what is going on. It is not hard!)

The arrays passed between the system and the integrator consist of a single Vec3f. The z coordinate does not do anything interesting but we kept it to make the various systems in this problem set more consistent.

The exact solution of the ODE is a circle with the equation

rcos(t + k)

X(t) =

rsin(t + k)

but that is not quite what we get from our implementation. Test this with example.exe. Pick the simple particle system, Euler integrator and a large step size. You’ll see the particle spiral outwardly in a 2D space, similar to the image below.

While large step size allows you to see the divergence faster, it will eventually occur regardless of how small you make the step size. The other integrators – midpoint, trape-zoidal and Runge-Kutta – are also unstable, but diverge at a much slower rate.

• Detailed instructions

Since we have supplied you with one working particle system (the simple system) and one working integrator (midpoint), you can test your work regardless of in which order you complete the requirements. Still, you’ll probably want to do R2, R4 and R5 in that order since they are progressively more di cult applications of the same techniques.

R1 Euler integrator (1p)

Implement an Euler integrator into the eulerStep function within integrators.cpp. Once it is correct, you’ll see the particle systems moving when you select the Euler integrator in the application.

R2 Spring system (2p)

For this requirement, you’ll implement a simple physically modeled particle system where one particle is xed in place and another hangs from it attached with a spring. Unlike the simple system, this equation involves a real spring force and is hence a 2nd-order system that needs to be reduced to 1st order by introducing the velocity variable like we saw above.

Head over to particle systems.cpp to ll in the missing parts of SpringSystem class implementation. First, tackle the reset function which sets the initial state of the particle system. You essentially just have to set the positions of the particles and store the data about the spring in the system. You need to pay attention to the indices where the components of position and velocity go in the state vector.

The real work happens in the evalF function which computes the forces a ecting the state given as a parameter. You can start by applying only gravity; you should see one particle staying put at origin and the other one falling. Then add the drag and spring forces. There are helper functions for the forces at the beginning of particle systems.cpp; fGravity is already implemented, while the other two, fDrag and fSpring are up to you.

Fixing a particle in place is just a matter of setting the net forces a ecting it to zero.

R3 Trapezoid integrator (2p)

Just like R1, only this time you implement a trapezoid integrator within the trapezoidStep function. The trapezoid integrator is considerably better than the Euler integrator, and should give results roughly equivalent to the midpoint integrator we supplied you with.

R4 Pendulum system (2p)

Here you’ll repeat what you did in R2 but in a more general form, building a particle system where an arbitrary amount of particles are connected with springs in a chain and suspended from the topmost one. Appropriate rendering could make such a system look like a length of rope.

You again need to ll in the reset and evalF functions. It is substantially no di erent from R2, but there are more particles (note that PendulumSystem’s constructor takes the particle count as a parameter) and it’s easier to make mistakes while indexing into the state array. You might nd it useful to de ne small helper functions, e.g. pos and vel, which take the index of a particle and return the index at which that particle’s position and velocity are stored in the state array.

At this point it’s also useful to consider forming the proper springs in the PendulumSystem::reset function. After this, in PendulumSystem::evalF you can loop rst over the particles to compute independent forces (gravity) and then loop over the springs, adding their re-sultant force to both a ected particles. This makes the implementation simpler; you

only have to think about the particle indices once instead of guring out the neighbors of each particle separately { just be careful to only add each spring once into the list. Furthermore, with this approach your evalF can be essentially the same in R5.

R5 Cloth system (3p)

What you have learned in R2 and R4 is su cient to implement a simulation of cloth in ClothSystem. We merely need to build a more elaborate con guration of particles and springs than we did before.

Figure 1: Left to right: structural springs, shear springs, and ex springs

We begin with a uniform grid of particles and connect them to their vertical and horizonal neighbors with springs. These springs keep the particle mesh together and are known as structural springs. Then, we add additional shear springs to prevent the cloth from collapsing diagonally. Finally, we add exion springs to prevent the cloth from folding over onto itself. Notice that the ex springs are only drawn as curves to make it clear that they skip a particle|they are still \straight” springs.

Write your code very carefully here; it is easy to make mistake and connect springs to particles that don’t exist. Again we recommend that you create helper methods, like the previously mentioned pos and vel, but this time taking the x and y grid coordinates of the particle and returning indices in the state vector belonging to that particle. As already suggested for R4, here it’s extremely useful to generate the springs in reset and simply evaluate their result in evalF instead of guring which neighboring particles exist for which particles.

First, implement structural springs. Make sure it looks as you expect before moving on. Run the simulation and you should obtain something that looks like a net. As usual, viscous drag helps prevent explosions. For faster debugging, use small meshes of e.g. 3 3 particles.

Once you’ve made sure your structural springs are correct, add in the shear springs, and then the ex springs. When inserting a new set of springs you can temporarily remove the previous ones to see whether the initial shape of the new ones looks correct.

Don’t be too discouraged if your rst tests blow up because of instability. Use the most stable integrator you have, stick to low step sizes, and use the example binary as a yardstick for how low the step sizes need to be. Doing the easy Runge-Kutta extra credit

or one of the harder extras will help a lot to keep the cloth together.

• Extra Credit

As always, you are free to do other extra credit work than what’s listed here | just make sure to describe what you did in your README. We’ll be fair with points, but if you want to attempt something grand, better talk to us beforehand. On some items, you have the choice to implement more or less features; we’ll give you points based on how much you did.

4.1 Recommended

Easy: Implement RK4 Runge-Kutta integrator (2p)

A RK4 integrator is nearly as easy to implement as the trapezoid integrator; it’s just more code to write. RK4 allows using substantially larger step sizes without your system blowing up. See http://en.wikipedia.org/wiki/Runge-Kutta methods

Easy: Spray/waterfall/sprinkler system (1-3p)

Create your own particle system that keeps spraying out particles which disappear after a while (think \water hose” or \rising smoke”). For full points, also have the particles bounce o a plane or other object (you don’t have to render the object). If you’d like to also render your particles using semi-transparent textures or other cool e ects, there’s an extra for that in the Medium section.

4.2 Easy

Wind (2p)

Add a random wind force with an on/o toggle to your cloth simulation.

Mouse drag/poke UI (2-3p)

Provide a mouse-based interface for users to interact with the cloth. You may, for instance, allow the user to click on certain parts of the cloth and drag parts around.

Frictionless collisions (2p)

Implement frictionless collisions of cloth with a simple primitive such as a sphere. This is simpler than it may sound at rst: just check whether a particle is \inside” the sphere; if so, just project the point back to the surface. (You don’t need to draw the primitive.)

Cloth tearing (1+p) Simulate the cloth breaking up by removing springs that have stretched beyond some threshold. Simply doing this with an uniform treshold is worth one point, and additional features like non-uniform tresholds (anisotropic or spatially varying) and generating a reasonable mesh with the boundaries handled nicely are worth more.

4.3 Medium

Particle editor (2-3p)

Give the user a mouse UI to modify some particle parameters using spline curves | for instance, the emission rate and dispersal of a geyser, or wind strength and direction over time. Provide a way to save and load these parameters from les.

Particle rendering (1-4p)

Create a cool rendering mode for some particle system. You could make your extra credit sprinkler/waterfall system emit particles that look like uid/ re/smoke, do a visualization of kinetic and potential energy or forces in one of your spring-based systems, or render the cloth mesh. A crude cloth mesh as in the example binary is worth 1p, but you can apply better shading or texture mapping for more points. We award points based on your best system, so invest in quality rather than quantity.

Adaptive solver scheme

Implement an adaptive solver scheme (look up adaptive Runge-Kutta-Felhberg tech-niques or check out the MATLAB ode45 function).

GPU simulation (4+p)

Transfer the computational load of the simulation to run on the graphics processor using OpenGL compute shaders or other suitable API (OpenCL, CUDA) { we have provided starter code for OpenGL so it might be the easiest especially if you have no former experience in GPU computing. The classroom computers have GL compute support in case your own computer doesn’t. You might want to see the slides for the Programming parallel computers course, lectures 4 and 5 are relevant. Although the slides focus on CUDA, many concepts apply to GL compute as well, just with di erent names. This OpenGL Wiki page gives a good overall description of GL compute shaders as well.

To get started, rst you’ll have to uncomment the line including ComputeCloth.hpp from the top of App.hpp. You’ll have to ll in function bodies to BxcPlusD.glsl, evalF.glsl and reset.glsl, as well as the body of the CPU-side method Advance in ComputeCloth.cpp. Each le contains some helpful comments and starter code that should help you get started.

Points are awarded based on which systems, integrators and additional features (like wind and tearing) you implement. You’ll get four points for the basic cloth system with no extra features.

4.4 Hard

Implicit integration (8-10p)

Implement an implicit integration scheme, such as implicit Euler. Such techniques allow much greater stability for sti systems of di erential equations, such as the ones that arise from cloth simulation. (Even though implicit techniques su er from numerical error just like explicit ones, the inaccuracy tends to bias the solution towards stable solutions, thus allowing far greater step sizes.) This will require

the ability to evaluate the Jacobian of the force vector with respect to its own components, and a Newton solver. We recommend you use the Eigen C++ linear algebra template library for solving the large linear systems that arise. Be sure to look at the Bara & Witkin notes. You get more points for implementing more than one kind of implicit solver (e.g., implicit midpoint rule).

Implicit optimization (4p)

The Jacobian matrix of the force function is required for implicit integration. While the numerical values of the matrix change for each iteration in the Newton method, the sparsity pattern stays the same; only certain entries of the matrix are ever non-zero. See if you can use Eigen’s \sparse direct solvers” to perform a symbolic factorization of the matrix beforehand. Then, given the numerical values in the current step, this allows you to quickly compute a matrix factorization that lets you solve for the Newton step quickly.

Constraints (8p)

Extend your particle system to support constraints, as described here. This means that you remove the springs, and using so-called Lagrange multipliers, project the remaining forces (gravity, etc.) onto a subspace that does not violate the con-straints.

Robust model of cloth (15p) Implement the technique described in http://graphics.stanford.edu/papers/cloth-sig02/

Simulate rigid-body dynamics http://www.cs.cmu.edu/ bara /sigcourse/notesd1.pdf or deformable models: http://graphics.cs.cmu.edu/projects/stvk/ or

uids: http://physbam.stanford.edu/~fedkiw/

In theory, particle systems can be used to achieve similar e ects. However, greater accuracy and e ciency can be achieved through more complex physical and math-ematical models. You will receive credit according to how large and hard your implementation is.

• Submission

Make sure your code compiles and runs both in Release and Debug modes on Visual Studio 2019, preferably in the VDI environment. Comment out any func-tionality that is so buggy it would prevent us seeing the good parts. Check that your README.txt (which you hopefully have been updating throughout your work) accurately describes the nal state of your code.

Fill in whatever else is missing, including any feedback you want to share. We were not kidding when we said we prefer brutally honest feedback.

Package all the code, project and solution les required to build your submission, the README.txt and any screenshots, logs or other les you want to share into a ZIP archive. Note: the solution le itself does not contain any code, make sure the contents of src/ are present in the nal submission.

Sanity check: look inside your ZIP archive. Are the les there? (Better yet, unpack the archive into another folder, and see if you can still open the solution, compile the code and run.)

Submit your archive in MyCourses folder:

“Assignment 4: Physical Simulation”.

12