12,244,322 members (46,060 online)
Article
alternative version

6.6K views
8 bookmarked
Posted

# Linear Rigid Body Dynamics with OpenGL

, 25 Jul 2013 CPOL
 Rate this:
Simulating displacement, velocity, and acceleration in code.

## Introduction

To animate 3-D models realistically, it's important to understand the principles of dynamics, which describe the motion of objects in the real world. Computational physics has become commonplace in games and simulation tools, and this article explains how to combine linear dynamics with OpenGL rendering.

Scientists and engineers have understood the laws of dynamics since the age of Newton, and programmers have become increasingly interested in using these principles to animate objects in 3-D environments. If you search on the web, you'll find plenty of material on the topic. In particular, I recommend Chris Hecker's articles and the lecture notes written by Andrew Witkin, David Baraff, and Michael Kass.

Many sources discuss the theory behind dynamics and some show how the equations can be implemented in code. But this article goes further and shows how rigid body dynamics can be implemented in a real-world application based on OpenGL. This application, whose code can be obtained here, assigns a velocity and acceleration to a sphere, and then renders the path of the sphere as it moves through space. Figure 1 shows what this looks like on my Linux system.

##### Figure 1: A Sphere's Motion Subject to Linear Dynamics

This may not look particularly impressive, but the methodology behind the code can be extended to more complex circumstances. But before I can discuss the code, I need to present the theory underlying linear dynamics. The first section discusses this in detail.

NOTE: This article focuses solely on linear dynamics, and does not discuss rotational dynamics. In addition, this presentation deals with difference equations only, not differential or integral equations.

## 1.  The Basics of Linear Dynamics

In many video games, objects always seem to move at the same speed. But in the real world, the speed of an object changes over time, particularly under the influence of gravity. To explain how this works, this section discusses position and velocity first, and then discusses the important concepts of acceleration and force.

### 1.1  Position and Velocity

If you're familiar with graphical modeling, you know that the position of a point is given as a set of coordinates that identify a location relative to another point called the origin. In three-dimensional Euclidean space, coordinates are given in (x, y, z) triples, where x identifies the distance from the origin along the x-axis, y identifies the distance from the origin along the y-axis, and z identifies the distance from the origin along the z-axis.

It's important to distinguish distance from position. If a point's position, denoted r, is given as (x, y, z), it's distance from the origin, denoted |r|, equals sqrt(x2 + y2 + z2). Distance identifies how far away a point is, and because it has a single value, it's called a scalar. In contrast, position identifies both how far away a point is and its direction. A quantity that identifies magnitude and direction is called a vector.

When you create animated designs, it's not enough to have a point's position. You have to know how its position changes from one frame to the next. We'll refer to the time interval between rendering frames as Δt, which is measured in seconds. For example, if an application's frame rate is 120 fps (frames per second), the time interval Δt = 1/120 = 0.008333 seconds.

The change in position over time is referred to as velocity. If a point's position changes by (Δx, Δy, Δz) over the time interval Δt, the average velocity over the interval is expressed as (Δx/Δt, Δy/Δt, Δz/Δt). To simplify the notation, we'll refer to the average velocity vector as v and its components will be given as vx = Δx/Δt, vy = Δy/Δt, and vz = Δz/Δt.

Note: This discussion deals exclusively with the average velocity over a time interval instead of the instantaneous velocity at a specific time. This is a big deal for physicists but not for programmers. Therefore, I'm going to stop using the phrase "average velocity over an interval."

In physics, velocity is expressed in meters per second. That is, if a point's velocity is (5, 5, 0), its position changes five meters in the +x direction and five meters in the +y direction every second. Meters aren't commonly used by programmers, so we'll refer to velocity in terms of generic units per second.

Suppose you know a point's velocity but you don't know how it's position changes from one frame to the next. In this situation, the following equation becomes helpful:

This states that, to find the new position, add the old position to the product of the velocity and time interval. This should make sense. If you drive at a speed of 64 mph (103 kph) for half an hour, you will travel a distance of 64 × 0.5 = 32 miles or 103 × 0.5 = 51.5 kilometers.

Figure 2 presents another example. The point's initial location is y0 = 0.5 and it travels upward at velocity vy = 1. Each time interval is 0.75s, so the first updated location, y1 = 0.5 + 1(0.75) = 1.25. The second location, y2 = 1.25 + 1(0.75) = 2. Each location exceeds the previous location by 0.75.

##### Figure 2: A Point's Trajectory with Constant Velocity

The path taken by a point over time is called its trajectory, and if the velocity is constant, the trajectory will be a straight line. The next discussion will explore what happens when an object's velocity changes over time.

### 1.2.  Acceleration and Force

Just as velocity measures the change in position over time, acceleration measures the change in velocity over time. Like position and velocity, acceleration is a vector. The following equation shows how it relates to velocity:

An important question arises. We know how to update a point's position given its velocity and we know how to update the velocity given the acceleration. But how do we update a point's position given its initial velocity and acceleration? That is, if we know ri, vi, and a, how do we obtain ri+1?

To obtain the answer, remember that the velocity changes from vi to vi+1 over the course of the time interval Δt. The average velocity over the interval can be computed as follows:

With this expression, we can determine the updated position as follows:

Now let's look at the trajectory of a point with constant acceleration. Suppose a point has r0 = 0.5, v0 = 2 and a = -1. After the first time interval, the updated position, r1, equals 0.5 + 2(0.75) + (0.5)(-1)(0.75)2 = 1.72. The updated velocity, v1, equals 2 + (-1)(0.75) = 1.25. After the second time interval, r2 = 1.72 + (1.25)(0.75) + (0.5)(-1)(0.75)2 = 2.375 and v2 = 1.25 + (-1)(0.75) = 0.5.

By keeping track of the cumulative time interval, we can obtain the position at each step without computing the velocity. For example, we can proceed from r0 to r2 by setting the time interval in the equation to 2Δt. That is, r2 = r0 + v0(2Δt) + (0.5)(a)(2Δt)2 = 2.375. Figure 3 shows the trajectory of this point from Δt to 6Δt.

##### Figure 3: A Point's Trajectory with Constant Acceleration

As shown in the figure, the trajectory of a point with initial velocity and constant acceleration is a parabola. This should make sense, as the equation for position is quadratic with respect to time and the graph of a quadratic equation forms a parabola.

Programmers and engineers generally limit their interest to position, velocity, and acceleration, and don't take into account the change in acceleration over time. This is because force plays a very large role in physics —from friction to gravity to electromagnetism, the phenomena that affect physical motion are typically quantified using forces. When a force acts on an object and causes it to move, the object's acceleration equals the force divided by the object's mass. In equation form, F = ma, where F is the force in Newtons, m is the mass in kilograms, and a is the acceleration in meters/second2.

At this point, you should have a clear understanding of the relationships between position, velocity, and acceleration. These are all given as vectors, and the next section will show how they can be used to constrain the motion of an object in an OpenGL rendering.

## 2.  Linear Dynamics and OpenGL

Modern OpenGL applications store vertex positions in structures called vertex buffer objects, or VBOs. As the rendering process proceeds, a vertex shader processes the elements of the VBO to determine the final position of the vertex. For this reason, we'll rely on the vertex shader to update the VBO data with the result of physical computation.

In particular, the goal of this discussion is to present an OpenGL application that animates a sphere so that it follows the trajectory depicted in Figure 1. Like the point in Figure 2, this sphere has an upward initial velocity but downward acceleration.

To update the positions of the sphere's vertices, three pieces of information must be processed: the elapsed time, the velocity, and the acceleration. In this section, I'll discuss how the application obtains the time and then proceed to show how the velocity, acceleration, and updated position are computed.

### 2.1  Obtaining Elapsed Time

To determine how vertex positions should be updated, the first step is to obtain the elapsed time. Different windowing systems and operating systems provide different timing functions, but the code in this example application is based on the OpenGL Utility Toolkit, or GLUT. The example application relies on two functions to obtain timing information:

1. `glutIdleFunc(void (*func)(void))` - identifies the function to be called when the application is idle (no events are being received or processed)
2. `glutGet(GLUT_ELAPSED_TIME)` - returns the number of milliseconds that have passed since `glutInit` was called or the first call to `glutGet(GLUT_ELAPSED_TIME)`

To configure timing in code, the main function in the animate_sphere.cpp file invokes the following lines of code:

```glutIdleFunc(update_vertices);
start_time = glutGet(GLUT_ELAPSED_TIME);```

The first line states that the `update_vertices` function should be called when the application is idle. The second line computes `start_time`, which sets the starting time in milliseconds.

The `update_vertices` function obtains the current time and uses it to compute the elapsed time. The following code shows how this works:

```void update_vertices() {
current_time = glutGet(GLUT_ELAPSED_TIME);
delta_t = (current_time - start_time)/1000.0f;
...
} ```

The `delta_t` variable contains the elapsed time that will be used in our physics equations. Note that this isn't the time between frames, but the total time that has passed since the application started.

### 2.2  Processing Velocity and Acceleration

After determining the timing data, the application computes the change in position brought about by the figure's velocity and acceleration. In animate_sphere.cpp, the dynamic parameters are initialized with the following code:

```#define INIT_POSITION 0.5f
#define INIT_VELOCITY 0.8f
#define ACCELERATION -0.4f
...
glm::vec3 init_position = glm::vec3(0.0f, INIT_POSITION, 0.0f);
glm::vec3 init_velocity = glm::vec3(INIT_VELOCITY, INIT_VELOCITY, 0.0f);
glm::vec3 acceleration = glm::vec3(0.0f, ACCELERATION, 0.0f);```

As shown, r0 is initialized to (0.0, 0.5, 0.0), v0 is initialized to (0.8, 0.8, 0.0), and a is set to (0.0, -0.4, 0.0). Each is represented by a `glm::vec3` data type, which is provided by the OpenGL Math Library.

The update_vertices function computes the updated position vector with the following code:

`delta_r = init_position + delta_t * init_velocity + 0.5f * delta_t * delta_t * acceleration;`

This is the same equation that we derived earlier, and should look familiar. The result identifies how each vertex position should be updated, and it can be used in one of at least two ways. First, `delta_r` can be added to the elements of the VBO containing vertex locations. Second, `delta_r` can be sent as a uniform to the vertex shader, which will add `delta_r` to each vertex's final location. In OpenGL, a uniform is a quantity that doesn't change from vertex to vertex.

The animate_sphere application adopts the second method, and the following code sets the uniform's data equal to `delta_r`:

`glUniform3fv(delta_location, 1, &(delta_r[0]));`

When the vertex shader executes, it adds the updated position to each vertex in the figure. You can see this by examining the vertex shader (animate_sphere.vert), which contains the following line of code:

`gl_Position = mvp * vec4(in_coords + delta, 1.0);`

This shows that the vertex's position is set equal to the original VBO position (`in_coords`) plus the updated position (`delta`) transformed by the modelview-projection matrix (`mvp`). As the shader executes, it updates the vertex locations and the result is presented way back in Figure 1.

## 3.  Conclusion

If you want objects in a graphical model to move realistically, you need to have a basic understanding of dynamics. This article has discussed the theory of linear dynamics for rigid bodies, and has focused on the relationships between position, velocity, and acceleration. The second part of this article discussed how this theory can be implemented in code. In particular, the animate_sphere application sets a sphere's initial velocity and acceleration, computes the dynamic equations in C, and uses OpenGL to render its motion through space.

## Share

 United States
No Biography provided