Open Dynamics Engine v0.035 Pre-release User Guide
Russell Smith, Tuesday 29 October, 2002

This document is Copyright © 2001,2002 Russell Smith.
This is an UNFINISHED document, there is still a lot of stuff to write.

Contents


1. Introduction

The Open Dynamics Engine (ODE) is a free, industrial quality library for simulating articulated rigid body dynamics. For example, it is good for simulating ground vehicles, legged creatures, and moving objects in VR environments. It is fast, flexible and robust, and it has built-in collision detection. ODE is being developed by Russell Smith.

If ``rigid body simulation'' does not make much sense to you, check out What is a Physics SDK?.

This is the user guide for the fourth prototype release of ODE (version 0.03). This release is to help continue testing the core algorithms and solidifying the API. There will probably only be a few (hopefully minor) changes to the API in the future.

1.1. Current features

ODE is best for simulating articulated rigid body structures. An articulated structure is created when rigid bodies of various shapes are connected together with joints of various kinds. Examples are ground vehicles (where the wheels are connected to the chassis) or legged creatures (where the legs are connected to the body).

ODE is designed to be used in interactive or real-time simulation. It is particularly good for simulating moving objects in changeable virtual reality environments. This is because it is fast, robust and stable, and the user has complete freedom to change the structure of the system even while the simulation is running.

ODE uses a highly stable integrator, so that the simulation errors should not grow out of control. The physical meaning of this is that the simulated system should not "explode" for no reason (believe me, this happens a lot with other simulators if you are not careful). ODE emphasizes speed and stability over physical accuracy.

ODE has hard contacts. This means that a special non-penetration constraint is used whenever two bodies collide. The alternative, used in many other simulators, is to use virtual springs to represent contacts. This is difficult to do right and extremely error-prone.

ODE has a built-in collision detection system. However you can ignore it and do your own collision detection if you want to. The current collision primitives are sphere, box, capped cylinder and plane - more collision objects will come later. ODE's collision system provides fast identification of potentially intersecting objects, through the concept of ``spaces''.

Here are the features:

1.2. New in the next release

1.3. New in the 0.03 release

1.4. New in the 0.025 release

1.5. New in the 0.02 release

1.6. Future features

1.6.1. Plan for 0.035 release

1.6.2. Plan for 0.04 release

1.6.3. Other stuff for later releases

1.7. ODE's license

ODE is Copyright © 2001,2002 Russell L. Smith. All rights reserved.

This library is free software; you can redistribute it and/or modify it under the terms of EITHER:

  1. The GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. The text of the GNU Lesser General Public License is included with this library in the file LICENSE.TXT.
  2. The BSD-style license that is included with this library in the file LICENSE-BSD.TXT.
This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the files LICENSE.TXT and LICENSE-BSD.TXT for more details.

1.8. The ODE community

Do you have questions or comments about ODE? Think you can help? Please write to the ODE mailing list.


2. How to install and use ODE

2.1. Installing ODE

Step 1: Unpack the ODE archive.

Step 2: Get the GNU make tool. Many Unix platforms come with this, although sometimes it is called gmake. I have provided a version of GNU make for windows here.

Step 3: Edit the settings in the file config/user-settings. The list of supported platforms is given in that file.

Step 4: Run make to configure and build ODE and the graphical test programs. Be sure you are using the correct version of GNU make. The make targets for building the parts of ODE are:

All of these targets will do an implicit make configure. If the configurator screws up then you can edit the settings directly in include/ode/config.h.

Step 5: To install the ODE library onto your system you should copy the lib/ and include/ directories to a suitable place, e.g. on Unix:

There is currently no built-in support to build Windows DLLs or Unix shared libraries --- although it is not hard to add this yourself.

ODE has been verified to build on at least the following platforms:

2.2. Using ODE

The best way to understand how to use ODE is to look at the test/example programs that come with it. Note the following things:


3. Concepts

3.1. Background

[Here is where I will write some background information about rigid body dynamics and simulation. But in the meantime, please refer to Baraff's excellent SIGGRAPH tutorial].

3.2. Rigid bodies

A rigid body has various properties from the point of view of the simulation. Some properties change over time:

Other body properties are usually constant over time: Conceptually each body has an x-y-z coordinate frame embedded in it, that moves and rotates with the body:

The origin of this coordinate frame is the body's point of reference. Some values in ODE (vectors, matrices etc) are relative to the body coordinate frame, and others are relative to the global coordinate frame.

Note that the shape of a rigid body is not a dynamical property (except insofar as it influences the various mass properties). It is only collision detection that cares about the detailed shape of the body.

3.2.1. Islands and disabled bodies

Bodies are connected to each other with joints. An ``island'' of bodies is a group that can not be pulled apart - in other words each body is connected somehow to every other body in the island.

Each island in the world is treated separately when the simulation step is taken. This is useful to know: if there are N similar islands in the simulation then the step computation time will be O(N).

Each body can be enabled or disabled. Disabled bodies are effectively ``turned off'' and are not updated during a simulation step. Disabling bodies is an effective way to save computation time when it is known that the bodies are motionless or otherwise irrelevent to the simulation.

If there are any enabled bodies in an island then every body in the island will be enabled at the next simulation step. Thus to effectively disable an island of bodies, every body in the island must be disabled. If a disabled island is touched by another enabled body then the entire island will be enabled, as a contact joint will join the enabled body to the island.

3.3. Integration

The process of simulating the rigid body system through time is called integration. Each integration step advances the current time by a given step size, adjusting the state of all the rigid bodies for the new time value. There are two main issues to consider when working with any integrator:

ODE's current integrator is very stable, but not particularly accurate unless the step size is small. For most uses of ODE this is not a problem -- ODE's behavior still looks perfectly physical in almost all cases. However ODE should not be used for quantitative engineering until this accuracy issue has been addressed in a future release.

3.4. Force accumulators

Between each integrator step the user can call functions to apply forces to the rigid body. These forces are added to "force accumulators" in the rigid body object. When the next integrator step happens, the sum of all the applied forces will be used to push the body around. The forces accumulators are set to zero after each integrator step.

3.5. Joints and constraints

In real life a joint is something like a hinge, that is used to connect two objects. In ODE a joint is very similar: It is a relationship that is enforced between two bodies so that they can only have certain positions and orientations relative to each other. This relationship is called a constraint -- the words joint and constraint are often used interchangeably. The following picture shows three different constraint types:

The first is a ball and socket joint that constraints the ``ball'' of one body to be in the same location as the ``socket'' of another body. The second is a hinge joint that constraints the two parts of the hinge to be in the same location and to line up along the hinge axle. The third is a slider joint that constraints the ``piston'' and ``socket'' to line up, and additionally constraints the two bodies to have the same orientation.

Each time the integrator takes a step all the joints are allowed to apply constraint forces to the bodies they affect. These forces are calculated such that the bodies move in such a way to preserve all the joint relationships.

Each joint has a number of parameters controlling its geometry. An example is the position of the ball-and-socket point for a ball-and-socket joint. The functions to set joint parameters all take global coordinates, not body-relative coordinates. A consequence of this is that the rigid bodies that a joint connects must be positioned correctly before the joint is attached.

3.6. Joint groups

A joint group is a special container that holds joints in a world. Joints can be added to a group, and then when those joints are no longer needed the entire group of joints can be very quickly destroyed with one function call. However, individual joints in a group can not be destroyed before the entire group is emptied.

This is most useful with contact joints, which are added and remove from the world in groups every time step.

3.7. Joint error and the error reduction parameter (ERP)

When a joint attaches two bodies, those bodies are required to have certain positions and orientations relative to each other. However, it is possible for the bodies to be in positions where the joint constraints are not met. This ``joint error'' can happen in two ways:

  1. If the user sets the position/orientation of one body without correctly setting the position/orientation of the other body.
  2. During the simulation, errors can creep in that result in the bodies drifting away from their required positions.
Here is an example of error in a ball and socket joint (where the ball and socket do not line up):

There is a mechanism to reduce joint error: during each simulation step each joint applies a special force to bring its bodies back into correct alignment. This force is controlled by the error reduction parameter (ERP), which has a value between 0 and 1.

The ERP specifies what proportion of the joint error will be fixed during the next simulation step. If ERP=0 then no correcting force is applied and the bodies will eventually drift apart as the simulation proceeds. If ERP=1 then the simulation will attempt to fix all joint error during the next time step. However, setting ERP=1 is not recommended, as the joint error will not be completely fixed due to various internal approximations. A value of ERP=0.1 to 0.8 is recommended (0.2 is the default).

A global ERP value can be set that affects most joints in the simulation. However some joints have local ERP values that control various aspects of the joint.

3.8. Soft constraint and constraint force mixing (CFM)

Most constraints are by nature ``hard''. This means that the constraints represent conditions that are never violated. For example, the ball must always be in the socket, and the two parts of the hinge must always be lined up. In practice constraints can be voilated by unintentional introduction of errors into the system, but the error reduction parameter can be set to correct these errors.

Not all constraints are hard. Some ``soft'' constraints are designed to be violated. For example, the contact constraint that prevents colliding objects from penetrating is hard by default, so it acts as though the colliding surfaces are made of steel. But it can be made into a soft constraint to simulate softer materials, thereby allowing some natural penetration of the two objects when they are forced together.

There are two parameters that control the distinction between hard and soft constraints. The first is the error reduction parameter (ERP) that has already been introduced. The second is the constraint force mixing (CFM) value, that is described below.

3.8.1. Constraint force mixing (CFM)

What follows is a somewhat technical description of the meaning of CFM. If you just want to know how it is used in practice then skip to the next section.

Traditionally the constraint equation for every joint has the form

J * v = c

where v is a velocity vector for the bodies involved, J is a ``Jacobian'' matrix with one row for every degree of freedom the joint removes from the system, and c is a right hand side vector. At the next time step, a vector lambda is calculated (of the same size as c) such that the forces applied to the bodies to preserve the joint constraint are

force = JT * lambda

ODE adds a new twist. ODE's constraint equation has the form

J * v = c + CFM * lambda

where CFM is a square diagonal matrix. CFM mixes the resulting constraint force in with the constraint that produces it. A nonzero (positive) value of CFM allows the original constraint equation to be violated by an amount proportional to CFM times the restoring force lamda that is needed to enforce the constraint. Solving for lambda gives

(J M-1 JT + CFM/h) lambda = c/h

Thus CFM simply adds to the diagonal of the original system matrix. Using a positive value of CFM has the additional benefit of taking the system away from any singularity and thus improving the factorizer accuracy.

3.8.2. How to use ERP and CFM

ERP and CFM can be independently set in many joints. They can be set in contact joints, in joint limits and various other places, to control the spongyness and springyness of the joint (or joint limit).

If CFM is set to zero, the constraint will be hard. If CFM is set to a positive value, it will be possible to violate the constraint by ``pushing on it'' (for example, for contact constraints by forcing the two contacting objects together). In other words the constraint will be soft, and the softness will increase as CFM increases. What is actually happeneng here is that the constraint is allowed to be violated by an amount proportional to CFM times the restoring force that is needed to enforce the constraint. Note that setting CFM to a negative value can have undesirable bad effects, such as instability. Don't do it.

By adjusting the values of ERP and CFM, you can achieve various effects. For example you can simulate springy constraints, where the two bodies oscillate as though connected by springs. Or you can simulate more spongy constraints, without the oscillation. In fact, ERP and CFM can be selected to have the same effect as any desired spring and damper constants. If you have a spring constant kp and damping constant kd, then the corresponding ODE constants are:

ERP = h k_p / (h k_p + k_d)
CFM = 1 / (h k_p + k_d)

where h is the stepsize. These values will give the same effect as a spring-and-damper system simulated with implicit first order integration.

Increasing CFM, especially the global CFM, can reduce the numerical errors in the simulation. If the system is near-singular, then this can markedly increase stability. In fact, if the system is mis-behaving, one of the first things to try is to increase the global CFM.

3.9. Collision handling

[There is a lot that needs to be written about collision handling.]

Collisions between bodies or between bodies and the static environment are handled as follows:

  1. Before each simulation step, the user calls collision detection functions to determine what is touching what. These functions return a list of contact points. Each contact point specifies a position in space, a surface normal vector, and a penetration depth.

  2. A special contact joint is created for each contact point. The contact joint is given extra information about the contact, for example the friction present at the contact surface, how bouncy or soft it is, and various other properties.

  3. The contact joints are put in a joint "group", which allows them to be added to and removed from the system very quickly. The simulation will speed goes down as the number of contacts goes up, so various strategies can be used to limit the number of contact points.

  4. A simulation step is taken.

  5. All contact joints are removed from the system.
Note that the built-in collision functions do not have to be used - other collision detection libraries can be used as long as they provide the right kinds of contact point information.

3.10. Typical simulation code

A typical simulation will proceed like this:

  1. Create a dynamics world.
  2. Create bodies in the dynamics world.
  3. Set the state (position etc) of all bodies.
  4. Create joints in the dynamics world.
  5. Attach the joints to the bodies.
  6. Set the parameters of all joints.
  7. Create a collision world and collision geometry objects, as necessary.
  8. Create a joint group to hold the contact joints.
  9. Loop:
    1. Apply forces to the bodies as necessary.
    2. Adjust the joint parameters as necessary.
    3. Call collision detection.
    4. Create a contact joint for every collision point, and put it in the contact joint group.
    5. Take a simulation step.
    6. Remove all joints in the contact joint group.
  10. Destroy the dynamics and collision worlds.

3.11. Physics model

The various methods and approximations that are used in ODE are discussed here.

3.11.1. Friction approximation

[We really need more pictures here.]

The Coulomb friction model is a simple, but effective way to model friction at contact points. It is a simple relationship between the normal and tangential forces present at a contact point (see the contact joint section for a description of these forces). The rule is:

| fT | <= mu * | fN |

where fN and fT are the normal and tangential force vectors respectively, and mu is the friction coefficient (typically a number around 1.0). This equation defines a "friction cone" - imagine a cone with fN as the axis and the contact point as the vertex. If the total friction force vector is within the cone then the contact is in "sticking mode", and the friction force is enough to prevent the contacting surfaces from moving with respect to each other. If the force vector is on the surface of the cone then the contact is in "sliding mode", and the friction force is typically not large enough to prevent the contacting surfaces from sliding. The parameter mu thus specifies the maximum ratio of tangential to normal force.

ODE's friction models are approximations to the friction cone, for reasons of efficiency. There are currently two approximations to chose from:

  1. The meaning of mu is changed so that it specifies the maximum friction (tangential) force that can be present at a contact, in either of the tangential friction directions. This is rather non physical because it is independent of the normal force, but it can be useful and it is the computationally cheapest option. Note that in this case mu is a force limit an must be chosen appropriate to the simulation.
  2. The friction cone is approximated by a friction pyramid aligned with the first and second friction directions [I really need a picture here]. A further approximation is made: first ODE computes the normal forces assuming that all the contacts are frictionless. Then it computes the maximum limits fm for the friction (tangential) forces from

    fm = mu * | fN |

    and then proceeds to solve for the entire system with these fixed limits (in a manner similar to approximation 1 above). This differs from a true friction pyramid in that the "effective" mu is not quite fixed. This approximation is easier to use as mu is a unitless ratio the same as the normal Coloumb friction coefficient, and thus can be set to a constant value around 1.0 without regard for the specific simulation.


4. Data types and conventions

4.1. The basic data types

The ODE library can be built to use either single or double precision floating point numbers. Single precision is faster and uses less memory, but the simulation will have more numerical error that can result in visible problems. You will get less accuracy and stability with single precision.

[must describe what factors influence accuracy and stability].

The floating point data type is dReal. Other commonly used types are dVector3, dVector4, dMatrix3, dMatrix4, dQuaternion.

4.2. Objects and IDs

There are various kinds of object that can be created:

Functions that deal with these objects take and return object IDs. The object ID types are dWorldID, dBodyID, etc.

4.3. Argument conventions

All 3-vectors (x,y,z) supplied to ``set'' functions are given as individual x,y,z arguments.

All 3-vector result arguments to get() function are pointers to arrays of dReal.

Larger vectors are always supplied and returned as pointers to arrays of dReal.

All coordinates are in the global frame except where otherwise specified.

4.4. C versus C++

The ODE library is written in C++, but its public interface is made of simple C functions, not classes. Why is this?

4.5. Debugging

The ODE library can be compiled in "debugging" or "release" mode. Debugging mode is slower, but function arguments are checked and many run-time tests are done to ensure internal consistency. Release mode is faster, but no checking is done.


5. World

The world object is a container for rigid bodies and joints. Objects in different worlds can not interact, for example rigid bodies from two different worlds can not collide.

All the objects in a world exist at the same point in time, thus one reason to use separate worlds is to simulate systems at different rates.

Most applications will only need one world.


dWorldID dWorldCreate();
Create a new, empty world and return its ID number.


void dWorldDestroy (dWorldID);
Destroy a world and everything in it. This includes all bodies, and all joints that are not part of a joint group. Joints that are part of a joint group will be deactivated, and can be destroyed by calling, for example, dJointGroupEmpty().


void dWorldSetGravity (dWorldID, dReal x, dReal y, dReal z);
void dWorldGetGravity (dWorldID, dVector3 gravity);
Set and get the world's global gravity vector. The units are m/s/s, so Earth's gravity vector would be (0,0,-9.81), assuming that +z is up. The default is no gravity, i.e. (0,0,0).


void dWorldSetERP (dWorldID, dReal erp);
dReal dWorldGetERP (dWorldID);
Set and get the global ERP value, that controls how much error correction is performed in each time step. Typical values are in the range 0.1--0.8. The default is 0.2.


void dWorldSetCFM (dWorldID, dReal cfm);
dReal dWorldGetCFM (dWorldID);
Set and get the global CFM (constraint force mixing) value. Typical values are in the range 1e-9 -- 1. The default is 1e-5 if single precision is being used, or 1e-10 if double precision is being used.


void dWorldStep (dWorldID, dReal stepsize);
Step the world.


void dWorldImpulseToForce (dWorldID, dReal stepsize,
			   dReal ix, dReal iy, dReal iz, dVector3 force);
If you want to apply a linear or angular impulse to a rigid body, instead of a force or a torque, then you can use this function to convert the desired impulse into a force/torque vector before calling the dBodyAdd... function.

This function is given the desired impulse as (ix,iy,iz) and puts the force vector in force. The current algorithm simply scales the impulse by 1/stepsize, where stepsize is the step size for the next step that will be taken.

This function is given a dWorldID because, in the future, the force computation may depend on integrator parameters that are set as properties of the world.


void dCloseODE();
This deallocates some extra memory used by ODE that can not be deallocated using the normal destroy functions, e.g. dWorldDestroy(). You can use this function at the end of your application to prevent memory leak checkers from complaining about ODE.


6. Rigid body functions

6.1. Creating and destroying


dBodyID dBodyCreate (dWorldID);
Create a body in the given world with default mass parameters at position (0,0,0). Return its ID.


void dBodyDestroy (dBodyID);
Destroy a body. All joints that are attached to this body will be put into limbo (i.e. unattached and not affecting the simulation, but they will NOT be deleted).

6.2. Position and orientation


void dBodySetPosition   (dBodyID, dReal x, dReal y, dReal z);
void dBodySetRotation   (dBodyID, const dMatrix3 R);
void dBodySetQuaternion (dBodyID, const dQuaternion q);
void dBodySetLinearVel  (dBodyID, dReal x, dReal y, dReal z);
void dBodySetAngularVel (dBodyID, dReal x, dReal y, dReal z);
const dReal * dBodyGetPosition   (dBodyID);
const dReal * dBodyGetRotation   (dBodyID);
const dReal * dBodyGetQuaternion (dBodyID);
const dReal * dBodyGetLinearVel  (dBodyID);
const dReal * dBodyGetAngularVel (dBodyID);
These functions set and get the position, rotation, linear and angular velocity of the body. After setting a group of bodies, the outcome of the simulation is undefined if the new configuration is inconsistent with the joints/constraints that are present. When getting, the returned values are pointers to internal data structures, so the vectors are valid until any changes are made to the rigid body system structure.

Hmmm. dBodyGetRotation returns a 4x3 rotation matrix.

6.3. Mass and force


void dBodySetMass (dBodyID, const dMass *mass);
void dBodyGetMass (dBodyID, dMass *mass);
Set/get the mass of the body (see the mass functions).


void dBodyAddForce            (dBodyID, dReal fx, dReal fy, dReal fz);
void dBodyAddTorque           (dBodyID, dReal fx, dReal fy, dReal fz);
void dBodyAddRelForce         (dBodyID, dReal fx, dReal fy, dReal fz);
void dBodyAddRelTorque        (dBodyID, dReal fx, dReal fy, dReal fz);
void dBodyAddForceAtPos       (dBodyID, dReal fx, dReal fy, dReal fz,
                                        dReal px, dReal py, dReal pz);
void dBodyAddForceAtRelPos    (dBodyID, dReal fx, dReal fy, dReal fz,
                                        dReal px, dReal py, dReal pz);
void dBodyAddRelForceAtPos    (dBodyID, dReal fx, dReal fy, dReal fz,
                                        dReal px, dReal py, dReal pz);
void dBodyAddRelForceAtRelPos (dBodyID, dReal fx, dReal fy, dReal fz,
                                        dReal px, dReal py, dReal pz);
Add forces to bodies (absolute or relative coordinates). The forces are accumulated on to each body, and the accumulators are zeroed after each time step.

The ...RelForce and ...RelTorque functions take force vectors that are relative to the body's own frame of reference.

The ...ForceAtPos and ...ForceAtRelPos functions take an extra position vector (in global or body-relative coordinates respectively) that specifies the point at which the force is applied. All other functions apply the force at the center of mass.


const dReal * dBodyGetForce  (dBodyID);
const dReal * dBodyGetTorque (dBodyID);
Return the current accumulated force and torque vector. The returned pointers point to an array of 3 dReals. The returned values are pointers to internal data structures, so the vectors are only valid until any changes are made to the rigid body system.


void dBodySetForce  (dBodyID b, dReal x, dReal y, dReal z);
void dBodySetTorque (dBodyID b, dReal x, dReal y, dReal z);
Set the body force and torque accumulation vectors. This is mostly useful to zero the force and torque for deactivated bodies before they are reactivated, in the case where the force-adding functions were called on them while they were deactivated.

6.4. Utility


void dBodyGetRelPointPos (dBodyID, dReal px, dReal py, dReal pz,
                          dVector3 result);
void dBodyGetRelPointVel (dBodyID, dReal px, dReal py, dReal pz,
                          dVector3 result);
void dBodyGetPointVel    (dBodyID, dReal px, dReal py, dReal pz,
                          dVector3 result);
Utility functions that take a point on a body (px,py,pz) and return that point's position or velocity in global coordinates (in result). The dBodyGetRelPointXXX functions are given the point in body relative coordinates, and the dBodyGetPointVel function is given the point in global coordinates.


void dBodyGetPosRelPoint (dBodyID, dReal px, dReal py, dReal pz,
                       	  dVector3 result);
This is the inverse of dBodyGetRelPointPos(). It takes a point in global coordinates (x,y,z) and returns the point's position in body-relative coordinates (result).


void dBodyVectorToWorld   (dBodyID, dReal px, dReal py, dReal pz,
                           dVector3 result);
void dBodyVectorFromWorld (dBodyID, dReal px, dReal py, dReal pz,
                           dVector3 result);
Given a vector expressed in the body (or world) coordinate system (x,y,z), rotate it to the world (or body) coordinate system (result).

6.5. Miscellaneous


void  dBodySetData (dBodyID, void *data);
void *dBodyGetData (dBodyID);
Get and set the body's user-data pointer.


void dBodyEnable (dBodyID);
void dBodyDisable (dBodyID);
int dBodyIsEnabled (dBodyID);
Enable and disable a body. Disabled bodies are effectively ``turned off'' and are not updated during a simulation step. However, if a disabled body is connected to island containing one or more enabled bodies then it will be re-enabled at the next simulation step.

dBodyIsEnabled return 1 if a body is enabled or 0 if it is disabled. New bodies are created in the enabled state.


void dBodySetFiniteRotationMode (dBodyID, int mode);
This function control the way a body's orientation is updated at each time step. The mode argument can be:


int dBodyGetFiniteRotationMode (dBodyID);
Return the current finite rotation mode of a body (0 or 1).


void dBodySetFiniteRotationAxis (dBodyID, dReal x, dReal y, dReal z);
This sets the finite rotation axis for a body. This is axis only has meaning when the finite rotation mode is set (see dBodySetFiniteRotationMode()).

If this axis is zero (0,0,0), full finite rotations are performed on the body.

If this axis is nonzero, the body is rotated by performing a partial finite rotation along the axis direction followed by an infitesimal rotation along an orthogonal direction.

This can be useful to alleviate certain sources of error caused by quickly spinning bodies. For example, if a car wheel is rotating at high speed you can call this function with the wheel's hinge axis as the argument to try and improve its behavior.


void dBodyGetFiniteRotationAxis (dBodyID, dVector3 result);
Return the current finite rotation axis of a body.


int dBodyGetNumJoints (dBodyID b);
Return the number of joints that are attached to this body.


dJointID dBodyGetJoint (dBodyID, int index);
Return a joint attached to this body, given by index. Valid indexes are 0 to n-1 where n is the value returned by dBodyGetNumJoints().


void dBodySetGravityMode (dBodyID b, int mode);
int dBodyGetGravityMode (dBodyID b);
Set/get whether the body is influenced by the world's gravity or not. If mode is nonzero it is, if mode is zero, it isn't. Newly created bodies are always influenced by the world's gravity.


7. Joint types and joint functions

7.1. Creating and destroying


dJointID dJointCreateBall (dWorldID, dJointGroupID);
dJointID dJointCreateHinge (dWorldID, dJointGroupID);
dJointID dJointCreateSlider (dWorldID, dJointGroupID);
dJointID dJointCreateContact (dWorldID, dJointGroupID,
                              const dContact *);
dJointID dJointCreateUniversal (dWorldID, dJointGroupID);
dJointID dJointCreateHinge2 (dWorldID, dJointGroupID);
dJointID dJointCreateFixed (dWorldID, dJointGroupID);
dJointID dJointCreateAMotor (dWorldID, dJointGroupID);
Create a new joint of a given type. The joint is initially in "limbo" (i.e. it has no effect on the simulation) because it does not connect to any bodies. The joint group ID is 0 to allocate the joint normally. If it is nonzero the joint is allocated in the given joint group. The contact joint will be initialized with the given dContact structure.


void dJointDestroy (dJointID);
Destroy a joint, disconnecting it from its attached bodies and removing it from the world. However, if the joint is a member of a group then this function has no effect - to destroy that joint the group must be emptied or destroyed.


dJointGroupID dJointGroupCreate (int max_size);
Create a joint group. The max_size argument is now unused and should be set to 0. It is kept for backwards compatibility.


void dJointGroupDestroy (dJointGroupID);
Destroy a joint group. All joints in the joint group will be destroyed.


void dJointGroupEmpty (dJointGroupID);
Empty a joint group. All joints in the joint group will be destroyed, but the joint group itself will not be destroyed.

7.2. Miscellaneous


void dJointAttach (dJointID, dBodyID body1, dBodyID body2);
Attach the joint to some new bodies. If the joint is already attached, it will be detached from the old bodies first. To attach this joint to only one body, set body1 or body2 to zero - a zero body refers to the static environment. Setting both bodies to zero puts the joint into "limbo", i.e. it will have no effect on the simulation.

Some joints, like hinge-2 need to be attached to two bodies to work.


void dJointSetData (dJointID, void *data);
void *dJointGetData (dJointID);
Get and set the joint's user-data pointer.


int dJointGetType (dJointID);
Get the joint's type. One of the following constants will be returned:

dJointTypeBallA ball-and-socket joint.
dJointTypeHingeA hinge joint.
dJointTypeSliderA slider joint.
dJointTypeContactA contact joint.
dJointTypeUniversalA universal joint.
dJointTypeHinge2A hinge-2 joint.
dJointTypeFixedA fixed joint.
dJointTypeAMotorAn angular motor joint.


dBodyID dJointGetBody (dJointID, int index);
Return the bodies that this joint connects. The index must be either 0 or 1, referring to either the ``first'' or ``second'' body. If one of these returned body IDs is zero, the joint connects the other body to the static environment. If both body IDs are zero, the joint is in ``limbo'' and has no effect on the simulation.


void dJointSetFeedback (dJointID, dJointFeedback *);
dJointFeedback *dJointGetFeedback (dJointID);
During the world time step, the forces that are applied by each joint are computed. These forces are added directly to the joined bodies, and the user normally has no way of telling which joint contributed how much force.

If this information is desired then the user can allocate a dJointFeedback structure and pass its pointer to the dJointSetFeedback() function. The feedback information structure is defined as follows:

typedef struct dJointFeedback {
  dVector3 f1;       // force that joint applies to body 1
  dVector3 t1;       // torque that joint applies to body 1
  dVector3 f2;       // force that joint applies to body 2
  dVector3 t2;       // torque that joint applies to body 2
} dJointFeedback;

During the time step any feedback structures that are attached to joints will be filled in with the joint's force and torque information. The dJointGetFeedback() function returns the current feedback structure pointer, or 0 if none is used (this is the default). dJointSetFeedback() can be passed 0 to disable feedback for that joint.

Now for some API design notes. It might seem strange to require that users perform the allocation of these structures. Why not just store the data statically in each joint? The reason is that not all users will use the feedback information, and even when it is used not all joints will need it. It will waste memory to store it statically, especially as this structure could grow to store a lot of extra information in the future.

Why not have ODE allocate the structure itself, at the user's request? The reason is that contact joints (which are created and destroyed every time step) would require a lot of time to be spent in memory allocation if feedback is required. Letting the user do the allocation means that a better allocation strategy can be provided, e.g simply allocating them out of a fixed array.

The alternative to this API is to have a joint-force callback. This would work of course, but it has a few problems. First, callbacks tend to pollute APIs and sometimes require the user to go through unnatural contortions to get the data to the right place. Second, this would expose ODE to being changed in the middle of a step (which would have bad consequences), and there would have to be some kind of guard against this or a debugging check for it - which would complicate things.


int dAreConnected (dBodyID, dBodyID);
Utility function: return 1 if the two bodies are connected together by a joint, otherwise return 0.

7.3. Joint parameter setting functions

7.3.1. Ball and socket

A ball and socket joint looks like this:


void dJointSetBallAnchor (dJointID, dReal x, dReal y, dReal z);
Set the joint anchor point.


void dJointGetBallAnchor (dJointID, dVector3 result);
Get the joint anchor point.

7.3.2. Hinge

A hinge joint looks like this:


void dJointSetHingeAnchor (dJointID, dReal x, dReal y, dReal z);
void dJointSetHingeAxis (dJointID, dReal x, dReal y, dReal z);
Set hinge anchor and axis parameters.


void dJointGetHingeAnchor (dJointID, dVector3 result);
void dJointGetHingeAxis (dJointID, dVector3 result);
Get hinge anchor and axis parameters.


dReal dJointGetHingeAngle (dJointID);
dReal dJointGetHingeAngleRate (dJointID);
Get the hinge angle and the time derivative of this value. The angle is measured between the two bodies, or between the body and the static environment. The angle will be between -pi..pi.

When the hinge anchor or axis is set, the current position of the attached bodies is examined and that position will be the zero angle.

7.3.3. Slider

A slider joint looks like this:


void dJointSetSliderAxis (dJointID, dReal x, dReal y, dReal z);
Set the slider axis parameter.


void dJointGetSliderAxis (dJointID, dVector3 result);
Get the slider axis parameter.


dReal dJointGetSliderPosition (dJointID);
dReal dJointGetSliderPositionRate (dJointID);
Get the slider linear position (i.e. the slider's ``extension'') and the time derivative of this value.

When the axis is set, the current position of the attached bodies is examined and that position will be the zero position.

7.3.4. Universal

A universal joint looks like this:

A universal joint is like a ball and socket joint that constrains an extra degree of rotational freedom. Given axis 1 on body 1, and axis 2 on body 2 that is perpendicular to axis 1, it keeps them perpendicular. In other words, rotation of the two bodies about the direction perpendicular to the two axes will be equal.

In the picture, the two bodies are joined together by a cross. Axis 1 is attached to body 1, and axis 2 is attached to body 2. The cross keeps these axes at 90 degrees, so if you grab body 1 and twist it, body 2 will twist as well.

Universal joints show up in cars, where the engine causes a shaft, the drive shaft, to rotate along its own axis. At some point you'd like to change the direction of the shaft. The problem is, if you just bend the shaft, then the part after the bend won't rotate about its own axis. So if you cut it at the bend location and insert a universal joint, you can use the constraint to force the second shaft to rotate about the same angle as the first shaft.

Another use of this joint is to attach the arms of a simple virtual creature to its body. Imagine a person holding their arms straight out. You may want the arm to be able to move up and down, and forward and back, but not to rotate about its own axis.

Here are the universal joint functions:


void dJointSetUniversalAnchor (dJointID, dReal x, dReal y, dReal z);
void dJointSetUniversalAxis1 (dJointID, dReal x, dReal y, dReal z);
void dJointSetUniversalAxis2 (dJointID, dReal x, dReal y, dReal z);
void dJointGetUniversalAnchor (dJointID, dVector3 result);
void dJointGetUniversalAxis1 (dJointID, dVector3 result);
void dJointGetUniversalAxis2 (dJointID, dVector3 result);
FIXME: I should eventually write a description of these functions here.

7.3.5. Hinge-2

A hinge-2 joint looks like this:

The hinge-2 joint is the same as two hinges connected in series, with different hinge axes. An example, shown in the above picture is the steering wheel of a car, where one axis allows the wheel to be steered and the other axis allows the wheel to rotate.

The hinge-2 joint has an anchor point and two hinge axes. Axis 1 is specified relative to body 1 (this would be the steering axis if body 1 is the chassis). Axis 2 is specified relative to body 2 (this would be the wheel axis if body 2 is the wheel).

Axis 1 can have joint limits and a motor, axis 2 can only have a motor.

Axis 1 can function as a suspension axis, i.e. the constraint can be compressible along that axis.


void dJointSetHinge2Anchor (dJointID, dReal x, dReal y, dReal z);
void dJointSetHinge2Axis1 (dJointID, dReal x, dReal y, dReal z);
void dJointSetHinge2Axis2 (dJointID, dReal x, dReal y, dReal z);
Set hinge-2 anchor and axis parameters. Axis 1 and axis 2 must not lie along the same line.


void dJointGetHinge2Anchor (dJointID, dVector3 result);
void dJointGetHinge2Axis1 (dJointID, dVector3 result);
void dJointGetHinge2Axis2 (dJointID, dVector3 result);
Get hinge-2 anchor and axis parameters.


dReal dJointGetHinge2Angle1 (dJointID);
dReal dJointGetHinge2Angle1Rate (dJointID);
dReal dJointGetHinge2Angle2Rate (dJointID);
Get the hinge-2 angles (around axis 1 and axis 2) and the time derivatives of these values.

When the anchor or axis is set, the current position of the attached bodies is examined and that position will be the zero angle.

7.3.6. Fixed

The fixed joint maintains a fixed relative position and orientation between two bodies, or between a body and the static environment. Using this joint is almost never a good idea in practice, except when debugging. If you need two bodies to be glued together it is better to represent that as a single body.

Currently the fixed joint does not support a non-identity relative rotation between two bodies, it only supports a relative offset.


void dJointSetFixed (dJointID);
Call this on the fixed joint after it has been attached to remember the current desired relative offset between the bodies.

7.3.7. Contact

A contact joint looks like this:

The contact joint prevents body 1 and body 2 from inter-penetrating at the contact point. It does this by only allowing the bodies to have an ``outgoing'' velocity in the direction of the contact normal. Contact joints typically have a lifetime of one time step. They are created and deleted in response to collision detection.

Contact joints can simulate friction at the contact by applying special forces in the two friction directions that are perpendicular to the normal.

When a contact joint is created, a dContact structure must be supplied. This has the following definition:
struct dContact {
  dSurfaceParameters surface;
  dContactGeom geom;
  dVector3 fdir1;
};
geom is a substructure that is set by the collision functions. It is described in the collision section.

fdir1 is a "first friction direction" vector that defines a direction along which frictional force is applied. It must be of unit length and perpendicular to the contact normal (so it is typically tangential to the contact surface). It should only be defined if the dContactFDir1 flag is set in surface.mode. The "second friction direction" is a vector computed to be perpendicular to both the contact normal and fdir1.

surface is a substructure that is set by the user. Its members define the properties of the colliding surfaces. It has the following members:

7.3.8. Angular Motor

An angular motor (AMotor) allows the relative angular velocities of two bodies to be controlled. The angular velocity can be controlled on up to three axes, allowing torque motors and stops to be set for rotation about those axes (see the ``Stops and motor parameters'' section below). This is mainly useful in conjunction will ball joints (which do not constrain the angular degrees of freedom at all), but it can be used in any situation where angular control is needed. To use an AMotor with a ball joint, simply attach it to the same two bodies that the ball joint is attached to.

The AMotor can be used in different modes. In dAMotorUser mode, the user directly sets the axes that the AMotor controls. In dAMotorEuler mode, AMotor computes the euler angles corresponding to the relative rotation, allowing euler angle torque motors and stops to be set. An AMotor joint with euler angles looks like this:

In this diagram, a0, a1 and a2 are the three axes along which angular motion is controlled. The green axes (including a0) are anchored to body 1. The blue axes (including a2) are anchored to body 2. To get the body 2 axes from the body 1 axes the following sequence of rotations is performed:

There is an important restriction when using euler angles: the theta1 angle must not be allowed to get outside the range -pi/2 ... pi/2. If this happens then the AMotor joint will become unstable (there is a singularity at +/- pi/2). Thus you must set the appropriate stops on axis number 1.


void dJointSetAMotorMode (dJointID, int mode);
int dJointGetAMotorMode (dJointID);
Set (and get) the angular motor mode. The mode parameter must be one of the following constants:

dAMotorUserThe AMotor axes and joint angle settings are entirely controlled by the user. This is the default mode.
dAMotorEulerEuler angles are automatically computed. The axis a1 is also automatically computed. The AMotor axes must be set correctly when in this mode, as described below. When this mode is initially set the current relative orientations of the bodies will correspond to all euler angles at zero.


void dJointSetAMotorNumAxes (dJointID, int num);
int dJointGetAMotorNumAxes (dJointID);
Set (and get) the number of angular axes that will be controlled by the AMotor. The argument num can range from 0 (which effectively deactivates the joint) to 3. This is automatically set to 3 in dAMotorEuler mode.


void dJointSetAMotorAxis (dJointID, int anum, int rel,
			  dReal x, dReal y, dReal z);
void dJointGetAMotorAxis (dJointID, int anum, dVector3 result);
int dJointGetAMotorAxisRel (dJointID, int anum);
Set (and get) the AMotor axes. The anum argument selects the axis to change (0,1 or 2). Each axis can have one of three ``relative orientation'' modes, selected by rel: The axis vector (x,y,z) is always specified in global coordinates regardless of the setting of rel. There are two GetAMotorAxis functions, one to return the axis and one to return the relative mode.

For dAMotorEuler mode:


void dJointSetAMotorAngle (dJointID, int anum, dReal angle);
Tell the AMotor what the current angle is along axis anum. This function should only be called in dAMotorUser mode, because in this mode the AMotor has no other way of knowing the joint angles. The angle information is needed if stops have been set along the axis, but it is not needed for axis motors.


dReal dJointGetAMotorAngle (dJointID, int anum);
Return the current angle for axis anum. In dAMotorUser mode this is simply the value that was set with dJointSetAMotorAngle. In dAMotorEuler mode this is the corresponding euler angle.


dReal dJointGetAMotorAngleRate (dJointID, int anum);
Return the current angle rate for axis anum. In dAMotorUser mode this is always zero, as not enough information is available. In dAMotorEuler mode this is the corresponding euler angle rate.

7.4. General

The joint geometry parameter setting functions should only be called after the joint has been attached to bodies, and those bodies have been correctly positioned, otherwise the joint may not be initialized correctly. If the joint is not already attached, these functions will do nothing.

For the parameter getting functions, if the system is out of alignment (i.e. there is some joint error) then the anchor/axis values will be correct with respect to body 1 only (or body 2 if you specified body 1 as zero in the dJointAttach() function).

The default anchor for all joints is (0,0,0). The default axis for all joints is (1,0,0).

When an axis is set it will be normalized to unit length. The adjusted axis is what the axis getting functions will return.

When measuring a joint angle or position, a value of zero corresponds to the initial position of the bodies relative to each other.

Note that there are no functions to set joint angles or positions (or their rates) directly, instead you must set the corresponding body positions and velocities.

7.5. Stop and motor parameters

When a joint is first created there is nothing to prevent it from moving through its entire range of motion. For example a hinge will be able to move through its entire angle, and a slider will slide to any length.

This range of motion can be limited by setting stops on the joint. The joint angle (or position) will be prevented from going below the low stop value, or from going above the high stop value. Note that a joint angle (or position) of zero corresponds to the initial body positions.

As well as stops, many joint types can have motors. A motor applies a torque (or force) to a joint's degree(s) of freedom to get it to pivot (or slide) at a desired speed. Motors have force limits, which means they can apply no more than a given maximum force/torque to the joint.

Motors have two parameters: a desired speed, and the maximum force that is available to reach that speed. This is a very simple model of real life motors, engines or servos. However, is it quite useful when modeling a motor (or engine or servo) that is geared down with a gearbox before being connected to the joint. Such devices are often controlled by setting a desired speed, and can only generate a maximum amount of power to achieve that speed (which corresponds to a certain amount of force available at the joint).

Motors can also be used to accurately model dry (or Coulomb) friction in joints. Simply set the desired velocity to zero and set the maximum force to some constant value - then all joint motion will be impeded by that force.

The alternative to using joint stops and motors is to simply apply forces to the affected bodies yourself. Applying motor forces is easy, and joint stops can be emulated with restraining spring forces. However applying forces directly is often not a good approach and can lead to severe stability problems if it is not done carefully.

Consider the case of applying a force to a body to achieve a desired velocity. To calculate this force you use information about the current velocity, something like this:

force = k * (desired speed - current speed)

This has several problems. First, the parameter k must be tuned by hand. If it is too low the body will take a long time to come up to speed. If it is too high the simulation will become unstable. Second, even if k is chosen well the body will still take a few time steps to come up to speed. Third, if any other ``external'' forces are being applied to the body, the desired velocity may never even be reached (a more complicated force equation would be needed, which would have extra parameters and its own problems).

Joint motors solve all these problems: they bring the body up to speed in one time step, provided that does not take more force than is allowed. Joint motors need no extra parameters because they are actually implemented as constraints. They can effectively see one time step into the future to work out the correct force. This makes joint motors more computationally expensive than computing the forces yourself, but they are much more robust and stable, and far less time consuming to design with. This is especially true with larger rigid body systems.

Similar arguments apply to joint stops.

7.5.1. Parameter functions

Here are the functions that set stop and motor parameters (as well as other kinds of parameters) on a joint:


void dJointSetHingeParam (dJointID, int parameter, dReal value);
void dJointSetSliderParam (dJointID, int parameter, dReal value);
void dJointSetHinge2Param (dJointID, int parameter, dReal value);
void dJointSetAMotorParam (dJointID, int parameter, dReal value);
dReal dJointGetHingeParam (dJointID, int parameter);
dReal dJointGetSliderParam (dJointID, int parameter);
dReal dJointGetHinge2Param (dJointID, int parameter);
dReal dJointGetAMotorParam (dJointID, int parameter);
Set/get limit/motor parameters for each joint type. The parameter numbers are:

dParamLoStopLow stop angle or position. Setting this to -dInfinity (the default value) turns off the low stop. For rotational joints, this stop must be greater than -Pi to be effective.
dParamHiStopHigh stop angle or position. Setting this to dInfinity (the default value) turns off the high stop. For rotational joints, this stop must be less than Pi to be effective. If the high stop is less than the low stop then both stops will be ineffective.
dParamVelDesired motor velocity (this will be an angular or linear velocity).
dParamFMaxThe maximum force or torque that the motor will use to achieve the desired velocity. This must always be greater than or equal to zero. Setting this to zero (the default value) turns off the motor.
dParamFudgeFactorThe current joint stop/motor implementation has a small problem: when the joint is at one stop and the motor is set to move it away from the stop, too much force may be applied for one time step, causing a ``jumping'' motion. This fudge factor is used to scale this excess force. It should have a value between zero and one (the default value). If the jumping motion is too visible in a joint, the value can be reduced. Making this value too small can prevent the motor from being able to move the joint away from a stop.
dParamBounceThe bouncyness of the stops. This is a restitution parameter in the range 0..1. 0 means the stops are not bouncy at all, 1 means maximum bouncyness.
dParamCFMThe constraint force mixing (CFM) value used when not at a stop.
dParamStopERPThe error reduction parameter (ERP) used by the stops.
dParamStopCFMThe constraint force mixing (CFM) value used by the stops. Together with the ERP value this can be used to get spongy or soft stops. Note that this is intended for unpowered joints, it does not really work as expected when a powered joint reaches its limit.
dParamSuspensionERPSuspension error reduction parameter (ERP). Currently this is only implemented on the hinge-2 joint.
dParamSuspensionCFMSuspension constraint force mixing (CFM) value. Currently this is only implemented on the hinge-2 joint.

If a particular parameter is not implemented by a given joint, setting it will have no effect.

These parameter names can be optionally followed by a digit (2 or 3) to indicate the second or third set of parameters, e.g. for the second axis in a hinge-2 joint, or the third axis in an AMotor joint. A constant dParamGroup is also defined such that:

dParamXi = dParamX + dParamGroup * (i-1)


8. Support functions

8.1. Rotation functions

Rigid body orientations are represented with quaternions. A quaternion is four numbers [cos(theta/2) sin(theta/2)*u] where theta is a rotation angle and `u' is a unit length rotation axis.

Every rigid body also has a 3x3 rotation matrix that is derived from the quaternion. The rotation matrix and the quaternion always match.

Some information about quaternions:

The following are utility functions for dealing with rotation matrices and quaternions.


void dRSetIdentity (dMatrix3 R);
Set R to the identity matrix (i.e. no rotation).


void dRFromAxisAndAngle (dMatrix3 R,
                         dReal ax, dReal ay, dReal az, dReal angle);
Compute the rotation matrix R as a rotation of angle radians along the axis (ax,ay,az).


void dRFromEulerAngles (dMatrix3 R,
                        dReal phi, dReal theta, dReal psi);
Compute the rotation matrix R from the three Euler rotation angles.


void dRFrom2Axes (dMatrix3 R, dReal ax, dReal ay, dReal az,
                  dReal bx, dReal by, dReal bz);
Compute the rotation matrix R from the two vectors `a' (ax,ay,az) and `b' (bx,by,bz). `a' and `b' are the desired x and y axes of the rotated coordinate system. If necessary, `a' and `b' will be made unit length, and `b' will be projected so that it is perpendicular to `a'. The desired z axis is the cross product of `a' and `b'.


void dQSetIdentity (dQuaternion q);
Set q to the identity rotation (i.e. no rotation).


void dQFromAxisAndAngle (dQuaternion q, dReal ax, dReal ay, dReal az,
                         dReal angle);
Compute q as a rotation of angle radians along the axis (ax,ay,az).


void dQMultiply0 (dQuaternion qa,
                  const dQuaternion qb, const dQuaternion qc);
void dQMultiply1 (dQuaternion qa,
                  const dQuaternion qb, const dQuaternion qc);
void dQMultiply2 (dQuaternion qa,
                  const dQuaternion qb, const dQuaternion qc);
void dQMultiply3 (dQuaternion qa,
                  const dQuaternion qb, const dQuaternion qc);
Set qa = qb*qc. This is that same as qa = rotation qc followed by rotation qa. The 0/1/2 versions are analogous to the multiply functions, i.e. 1 uses the inverse of qb, and 2 uses the inverse of qc. Option 3 uses the inverse of both.


void dQtoR (const dQuaternion q, dMatrix3 R);
Convert quaternion q to rotation matrix R.


void dRtoQ (const dMatrix3 R, dQuaternion q);
Convert rotation matrix R to quaternion q.


void dWtoDQ (const dVector3 w, const dQuaternion q, dVector4 dq);
Given an existing orientation q and an angular velocity vector w, return in dq the resulting dq/dt.

8.2. Mass functions

The mass parameters of a rigid body are described by a dMass structure:
typedef struct dMass {
  dReal mass;   // total mass of the rigid body
  dVector4 c;   // center of gravity position in body frame (x,y,z)
  dMatrix3 I;   // 3x3 inertia tensor in body frame, about POR
} dMass;

The following functions operate on this structure:


void dMassSetZero (dMass *);
Set all the mass parameters to zero.


void dMassSetParameters (dMass *, dReal themass,
                         dReal cgx, dReal cgy, dReal cgz,
                         dReal I11, dReal I22, dReal I33,
                         dReal I12, dReal I13, dReal I23);
Set the mass parameters to the given values. themass is the mass of the body. (cx,cy,cz) is the center of gravity position in the body frame. The Ixx values are the elements of the inertia matrix:
    [ I11 I12 I13 ]
    [ I12 I22 I23 ]
    [ I13 I23 I33 ]


void dMassSetSphere (dMass *, dReal density, dReal radius);
Set the mass parameters to represent a sphere of the given radius and density, with the center of mass at (0,0,0) relative to the body.


void dMassSetCappedCylinder (dMass *, dReal density, int direction,
                             dReal a, dReal b);
Set the mass parameters to represent a capped cylinder of the given parameters and density, with the center of mass at (0,0,0) relative to the body. The radius of the cylinder (and the spherical cap) is a. The length of the cylinder (not counting the spherical cap) is b. The cylinder's long axis is oriented along the body's x, y or z axis according to the value of direction (1=x, 2=y, 3=z).


void dMassSetBox (dMass *, dReal density,
                  dReal lx, dReal ly, dReal lz);
Set the mass parameters to represent a box of the given dimensions and density, with the center of mass at (0,0,0) relative to the body. The side lengths of the box along the x, y and z axes are lx, ly and lz.


void dMassAdjust (dMass *, dReal newmass);
Given mass parameters for some object, adjust them so the total mass is now newmass. This is useful when using the above functions to set the mass parameters for certain objects - they take the object density, not the total mass.


void dMassTranslate (dMass *, dReal x, dReal y, dReal z);
Given mass parameters for some object, adjust them to represent the object displaced by (x,y,z) relative to the body frame.


void dMassRotate (dMass *, const dMatrix3 R);
Given mass parameters for some object, adjust them to represent the object rotated by R relative to the body frame.


void dMassAdd (dMass *a, const dMass *b);
Add the mass b to the mass a.

8.3. Math functions

[There are quite a lot of these, but they're not standardized enough to document yet].

8.4. Error and memory functions

[Document these later].


9. Collision functions

This chapter describes the built-in collision detection system of ODE. Using ODE's collision detection is optional - an alternative collision detection system can be used as long as it can supply the right kinds of information to ODE.

9.1. Contact points

If two bodies touch, or if a body touches a static feature in its environment, the contact is represented by one or more "contact points". Each contact point has a corresponding dContactGeom structure:
struct dContactGeom {
  dVector3 pos;       // contact position
  dVector3 normal;    // normal vector
  dReal depth;        // penetration depth
  dGeomID g1,g2;      // the colliding geoms
};
pos records the contact position, in global coordinates.

depth is the depth to which the two bodies inter-penetrate each other. If the depth is zero then the two bodies have a grazing contact, i.e. they "only just" touch. However, this is rare - the simulation is not perfectly accurate and will often step the bodies too far so that the depth is nonzero.

normal is a unit length vector that is, generally speaking, perpendicular to the contact surface.

g1 and g2 are the geometry objects that collided.

The convention is that if body 1 is moved along the normal vector by a distance depth (or equivalently if body 2 is moved the same distance in the opposite direction) then the contact depth will be reduced to zero. This means that the normal vector points "in" to body 1.

In real life, contact between two bodies is a sophisticated thing. Representing contacts by contact points is only an approximation. Contact "patches" or "surfaces" might be more physically accurate, but representing these things in high speed simulation software is a challenge.

Each extra contact point added to the simulation will slow it down some more, so sometimes we are forced to ignore contact points in the interests of speed. For example, when two boxes collide many contact points may be needed to properly represent the geometry of the situation, but we may choose to keep only the best three. Thus we are piling approximation on top of approximation.

9.2. Geometry objects

To use ODE's collision detection, geometry objects must be associated with the rigid bodies. A geometry object represents a rigid shape in space. Geometry objects are distinct from rigid bodies in that a geometry object has geometrical properties (size, shape, position and orientation) but no dynamical properties (such as velocity or mass).

Every geometry object is an instance of a class, such as sphere, plane, or box. You can define your own classes as well.

Every geometry object has a position vector and a 3*3 rotation matrix. If a geometry object is associated with a rigid body then its position and rotation is actually the position and rotation of that body. The point of reference for the standard classes usually corresponds to their centers of mass. This makes them particularly easy to connect to dynamics objects. If other points of reference are required, transformation objects should be used to encapsulate the primitives.

9.2.1. Sphere class


dGeomID dCreateSphere (dSpaceID space, dReal radius);
Create a sphere geometry object of the given radius, and return its ID. If space is nonzero, insert it into that space. The point of reference for a sphere is its center.


void dGeomSphereSetRadius (dGeomID sphere, dReal radius);
Set the radius of the given sphere.


dReal dGeomSphereGetRadius (dGeomID sphere);
Return the radius of the given sphere.

9.2.2. Box class


dGeomID dCreateBox (dSpaceID space, dReal lx, dReal ly, dReal lz);
Create a box geometry object of the given x/y/z side lengths (lx,ly,lz), and return its ID. If space is nonzero, insert it into that space. The point of reference for a box is its center.


void dGeomBoxSetLengths (dGeomID box, dReal lx, dReal ly, dReal lz);
Set the side lengths of the given box.


void dGeomBoxGetLengths (dGeomID box, dVector3 result);
Return in result the side lengths of the given box.

9.2.3. Plane class


dGeomID dCreatePlane (dSpaceID space,
                      dReal a, dReal b, dReal c, dReal d);
Create a plane geometry object of the given parameters, and return its ID. If space is nonzero, insert it into that space. The plane equation is
a*x+b*y+c*z = d
The plane's normal vector is (a,b,c), and it must have length 1. Unlike other geometry objects, planes disregard their assigned position and rotation, i.e. the parameters are always in global coordinates. In other words it is assumed that the plane is always part of the static environment and not tied to any movable object.


void dGeomPlaneSetParams (dGeomID plane, dReal a, dReal b, dReal c, dReal d);
Set the parameters of the given plane.


void dGeomPlaneGetParams (dGeomID plane, dVector4 result);
Return in result the parameters of the given plane.

9.2.4. Capped cylinder class


dGeomID dCreateCCylinder (dSpaceID space, dReal radius, dReal length);
Create a capped cylinder geometry object of the given parameters, and return its ID. If space is nonzero, insert it into that space.

A capped cylinder is like a normal cylinder except it has half-sphere caps at its ends. This feature makes collision detection particularly easy. The cylinder's length, not counting the caps, is given by length. The cylinder is aligned along the geometry's local Z axis. The radius of the caps, and of the cylinder itself, is given by radius.


void dGeomCCylinderSetParams (dGeomID ccylinder,
                              dReal radius, dReal length);
Set the parameters of the given capped cylinder.


void dGeomCCylinderGetParams (dGeomID ccylinder,
                              dReal *radius, dReal *length);
Return in radius and length the parameters of the given capped cylinder.

9.2.5. Geometry Group class

Geometry object groups (GeomGroups) are a way to speed up collision detection in some situations. They represent a collection of geometry objects that are located together in space.

GeomGroups solve the following problem: Suppose you have two cars driving over some terrain. Each car is made up of many bodies. Without GeomGroups, the collision computation time between the two cars is proportional to the total number of parts (or even to the square of this number, depending on which space container is used). To speed up the process, a GeomGroup object is created to represent each car. The GeomGroup is inserted into the space, but the individual geometry objects of the cars are not. The bounding box of each GeomGroup is computed to encompass the car geometry objects. At each time step, the dSpaceCollide() function does a single intersection test between the GeomGroups. If there is an intersection, the contents of the GeomGroup will be individually tested against each other. In this way the best case performance (when the cars are far away from each other) is improved.

GeomGroups have the same interface as any other geometry object.

The geometry objects inside a GeomGroup will not self-intersect if they are not a member of a space. If you want them to self-intersect, insert them into a separate space and call dSpaceCollide() on that. Do not insert them into the same space the GeomGroup is in.


dGeomID dCreateGeomGroup (dSpaceID space);
Create an empty geometry group object, and return its ID. If space is nonzero, insert it into that space.


void dGeomGroupAdd (dGeomID group, dGeomID x);
Add the geometry object x to the given GeomGroup. Note that when the GeomGroup is destroyed, the objects that it contains will not be.


void dGeomGroupRemove (dGeomID group, dGeomID x);
Remove the geometry object x from the given GeomGroup. If it is not in the GeomGroup then nothing is done.


int dGeomGroupGetNumGeoms (dGeomID group);
Return the number of geometry object in the given GeomGroup.


dGeomID dGeomGroupGetGeom (dGeomID group, int i);
Return a specific geometry object from the given GeomGroup. i is the index of the object, it must be between 0 and n-1 where n is the number of objects in the group.


int dGeomGroupQuery (dGeomID group, dGeomID x);
Return 1 if x is contained in the given GeomGroup, or 0 if it is not.

9.2.6. Geometry Transform class

A geometry transform object (T) encapsulates another geometry object (E), allowing E to be positioned and rotated arbitrarily with respect to its point of reference.

Most geometry objects (like the sphere and box) have their point of reference corresponding to their center of mass, allowing them to be easily connected to dynamics objects. Transform objects give you more flexibiliy - for example, you can offset the center of a sphere, or rotate a cylinder so that its axis is something else.

T mimics the object E that it encapsulates: T is inserted into a space and attached to a body as though it was E. E itself must not be inserted into a space or attached to a body. E's position and rotation are set to constant values that say how it is transformed relative to T. If E's position and rotation are left at their default values, T will behave exactly like E would have if you had used it directly.


dGeomID dCreateGeomTransform (dSpaceID space);
Create a new geometry transform object, and return its ID. If space is nonzero, insert it into that space. On creation the encapsulated geometry is set to 0.


void dGeomTransformSetGeom (dGeomID g, dGeomID obj);
Set the geometry object that the geometry transform g encapsulates. The object obj must not be inserted into any space, and must not be associated with any body.

If g has its clean-up mode turned on, and it already encapsulates an object, the old object will be destroyed before it is replaced with the new one.


dGeomID dGeomTransformGetGeom (dGeomID g);
Get the geometry object that the geometry transform g encapsulates.


void dGeomTransformSetCleanup (dGeomID g, int mode);
int dGeomTransformGetCleanup (dGeomID g);
Set and get the clean-up mode of geometry transform g. If the clean-up mode is 1, then the encapsulated object will be destroyed when the geometry transform is destroyed. If the clean-up mode is 0 this does not happen. The default clean-up mode is 0.


void dGeomTransformSetInfo (dGeomID g, int mode);
int dGeomTransformGetInfo (dGeomID g);
Set and get the "information" mode of geometry transform g. The mode can be 0 or 1. The default mode is 0.

With mode 0, when a transform object is collided with another object (using dCollide (tx_geom,other_geom,...)), the g1 field of the dContactGeom structure is set to the geometry object that is encapsulated by the transform object. This value of g1 allows the caller to interrogate the type of the geometry that is transformed, but it does not allow the caller to determine the position in global coordinates or the associated body, as both of these properties are used differently for encapsulated geometry objects.

With mode 1, the g1 field of the dContactGeom structure is set to the transform object itself. This makes the object appear just like any other kind of geometry object, as dGeomGetBody() will return the attached body, and dGeomGetPosition() will return the global position. To get the actual type of the encapsulated geometry in this case, dGeomTransformGetGeom() must be used.

9.2.7. General geometry object functions


void dGeomDestroy (dGeomID);
Destroy a geometry object, removing it from any space it is in first.


void dGeomSetData (dGeomID, void *);
void *dGeomGetData (dGeomID);
These functions set and get the user-defined data pointer stored in the geometry object.


void dGeomSetBody (dGeomID, dBodyID);
dBodyID dGeomGetBody (dGeomID);
These functions set and get the body associated with the geometry object. Setting the body automatically attaches the position vector and rotation matrix of the body to the geometry object. Setting a body ID of zero gives the geometry object its own position and rotation, independent from any body.


void dGeomSetPosition (dGeomID, dReal x, dReal y, dReal z);
void dGeomSetRotation (dGeomID, const dMatrix3 R);
Set the position vector and rotation matrix of the geometry object. These functions are analogous to dBodySetPosition() and dBodySetRotation(). If the geometry object is attached to a body, the body's position / rotation will also be changed.


const dReal * dGeomGetPosition (dGeomID);
const dReal * dGeomGetRotation (dGeomID);
Return pointers to the geometry object's position vector and rotation matrix. The returned values are pointers to internal data structures, so the vectors are valid until any changes are made to the geometry object. If the geometry object is attached to a body, the body's position / rotation pointers will be returned, i.e. the result will be identical to calling dBodyGetPosition() or dBodyGetRotation().


void dGeomGetAABB (dGeomID, dReal aabb[6]);
Return in aabb the axis aligned bounding box for the given geometry object.


dReal *dGeomGetSpaceAABB (dGeomID);
Return a pointer to the axis aligned bounding box that was precomputed for the given geometry object by the dSpaceCollide() function. This returns 0 if the dSpaceCollide() function is not currently being called. This is much faster than calling dGeomGetAABB().

9.2.8. Composite objects

Consider the following collision objects:

If these objects are meant to be rigid then it is necessary to use a single rigid body to represent each of them. But there are no single geometry classes that can represent things like a table or a molecule. The solution is to use a composite collision object that is a combination of several geometry objects.

No extra functions are needed to manage composite objects: simply create each component geometry object and attach it to the same body. To move and rotate the component objects with respect to each other, use geometry transforms to encapsulate them. That's all there is to it!

However there is one caveat: You should never create a composite object that will result in collision points being generated very close together. For example, consider a table that is made up of a box for the top and four boxes for the legs. If the legs are flush with the top, and the table is lying on the ground on its side, then the contact points generated for the boxes may coincide where the legs join to the top. ODE does not currently optimize away coincident contact points, so this situation can lead to numerical errors and strange behavior.

In this example the table geometry should be adjusted so that the legs are not flush with the sides, making it much more unlikely that coincident contact points will be generated. In general, avoid having different contact surfaces that overlap, or that line up along their edges.

9.3. Collision detection


int dCollide (dGeomID o1, dGeomID o2, int flags,
              dContactGeom *contact, int skip);
Given two geometry objects that potentially touch (o1 and o2), generate contact information for them. Internally, this just calls the correct class-specific collision functions for o1 and o2.

flags specifies how contacts should be generated if the objects touch. Currently the lower 16 bits of flags specifies the maximum number of contact points to generate. If this number is zero, this function just pretends that it is one - in other words you can not ask for zero contacts. All other bits in flags must be zero. In the future the other bits may be used to select other contact generation strategies.

contact must point to an array of contact geometry information that can hold at least the maximum number of contacts. Note: the elements of the contact array do not necessarily have to be contiguous. skip is the number of bytes between each dContactGeom structure in the contact array. If skip is sizeof(dContactGeom) then contact points to a "normal" (C-style) contact array. If skip is larger than this, then the dContactGeom structures are embedded in some other larger structures. It is an error for skip to be smaller than sizeof(dContactGeom).

If the objects touch, this returns the number of contact points generated (and updates the contact array), otherwise it returns 0 (and the contact array is not touched).

9.4. User defined classes

You can define your own geometry classes using the functions in this section. The standard geometry classes do not have any special access to the internals of ODE, they use the public functions exactly as you would.

Every geometry class has a unique number (0,1,2, etc...). A new geometry class (call it `X') must provide the following to ODE:

  1. Functions that will handle collision detection and contact generation between X and one or more other classes. These functions must be of type dColliderFn, which is defined as
    typedef int dColliderFn (dGeomID o1, dGeomID o2, int flags,
                             dContactGeom *contact, int skip);
    	
    This has exactly the same interface as dCollide(). Each function will handle a specific collision case, where o1 has type X and o2 has some other known type.

  2. A "selector" function, of type dGetColliderFnFn, which is defined as
    typedef dColliderFn * dGetColliderFnFn (int num);
    	
    This function takes a class number (num), and returns the collider function that can handle colliding X with class num. It should return 0 if X does not know how to collide with class num. Note that if classes X and Y are to collide, only one needs to provide a function to collide with the other.

    This function is called infrequently - the return values are cached and reused.

  3. A function that will compute the axis aligned bounding box (AABB) of instances of this class. This function must be of type dGetAABBFn, which is defined as
    typedef void dGetAABBFn (dGeomID g, dReal aabb[6]);
    	
    This function is given g, which has type X, and returns the axis-aligned bounding box for g. The aabb array has elements (minx, maxx, miny, maxy, minz, maxz). If you don't want to compute tight bounds for the AABB, you can just supply a pointer to dInfiniteAABB(), which returns +/- infinity in each direction.

  4. The number of bytes of "class data" that instances of this class need. For example a sphere stores its radius in the class data area, and a box stores its side lengths there.
The following things are optional for a geometry class:
  1. A function that will destroy the class data. Most classes will not need this function, but some will want to deallocate heap memory or release other resources. This function must be of type dGeomDtorFn, which is defined as
    typedef void dGeomDtorFn (dGeomID o);
    	
    The argument o has type X.

  2. A function that will test whether a given AABB intersects with an instance of X. This is used as an early-exit test in the space collision functions. This function must be of type dAABBTestFn, which is defined as
    typedef int dAABBTestFn (dGeomID o1, dGeomID o2, dReal aabb2[6]);
    	
    The argument o1 has type X. If this function is provided it is called by dSpaceCollide() when o1 intersects geometry object o2, which has an AABB given by aabb2. It returns 1 if aabb2 intersects o1, or 0 if it does not.

    This is useful, for example, for large terrains. Terrains typically have very large AABBs, which are not very useful to test intersections with other objects. This function can test another object's AABB against the terrain without going to the computational trouble of calling the specific collision function. This has an especially big savings when testing against GeomGroup objects.

Here are the functions used to manage custom classes:


int dCreateGeomClass (const dGeomClass *classptr);
Register a new geometry class, defined by classptr. The number of the new class is returned. The convention used in ODE is to assign the class number to a global variable with the name dXxxClass where Xxx is the class name (e.g. dSphereClass).

Here is the definition of the dGeomClass structure:
struct dGeomClass {
  int bytes;                  // bytes of custom data needed
  dGetColliderFnFn *collider; // collider function
  dGetAABBFn *aabb;           // bounding box function
  dAABBTestFn *aabb_test;     // aabb tester, can be 0 for none
  dGeomDtorFn *dtor;          // destructor, can be 0 for none
};


void * dGeomGetClassData (dGeomID);
Given a geometry object, return a pointer to the class's custom data (this will be a block of the required number of bytes). Internal ODE classes may use a slightly faster macro instead.


dGeomID dCreateGeom (int classnum);
Create a geometry object of the given class number. The custom data block will initially be set to 0. This object can be added to a space using dSpaceAdd().


int dGeomGetClass (dGeomID);
Given a geometry object, this returns its class number.

When you implement a new class you will usually write a function that does the following:

  1. If the class has not yet been created, create it. You should be careful to only ever create the class once.
  2. Call dCreateGeom() to make an instance of the class.
  3. Set up the custom data area.
This is what dCreateSphere() and the other geometry creation functions do.

9.5. Utility functions


void dClosestLineSegmentPoints (dVector3 const a1, dVector3 const a2,
				dVector3 const b1, dVector3 const b2,
				dVector3 cp1, dVector3 cp2)
Given two line segments A and B with endpoints a1-a2 and b1-b2, return the points on A and B that are closest to each other (in cp1 and cp2). In the case of parallel lines where there are multiple solutions, a solution involving the endpoint of at least one line will be returned. This will work correctly for zero length lines, e.g. if a1==a2 and/or b1==b2.


int dBoxTouchesBox (const dVector3 _p1, const dMatrix3 R1,
                    const dVector3 side1, const dVector3 _p2,
                    const dMatrix3 R2, const dVector3 side2);
Given boxes (p1,R1,side1) and (p2,R2,side2), return 1 if they intersect or 0 if not. p is the center of the box, R is the rotation matrix for the box, and side is a vector of x/y/z side lengths.


void dInfiniteAABB (dGeomID geom, dReal aabb[6]);
This function can be used as the AABB-getting function in a geometry class, if you don't want to compute tight bounds for the AABB. It returns +/- infinity in each direction.

9.6. Space

A space is a container for geometry objects. It is similar to the rigid body "world", except for collision instead of dynamics.

The space does high level collision culling, which means that it can identify which pairs of geometry objects are potentially touching. You can safely call dCollide() for only those pairs, instead of having to call dCollide() for every object-object pair. This can save a huge amount of time.

Two collision culling algorithms are currently available:

All spaces guarantee that if a geometry object's collider function is called, the object's AABB function will have been called prior to that without the geometry position having being changed in between the calls. This protocol can be useful if the AABB function wants to precompute some data that is used by the collider function. The collider function can check if it is being called inside a space by checking if the return value of dGeomGetSpaceAABB() is nonzero.

Here are the functions used for spaces:


dSpaceID dSimpleSpaceCreate();
dSpaceID dHashSpaceCreate();
Create a space, either of the simple or multi-resolution hash table kind.


void dHashSpaceSetLevels (dSpaceID space, int minlevel, int maxlevel);
This sets some parameters for a multi-resolution hash table space. The smallest and largest cell sizes used in the hash table will be 2^minlevel and 2^maxlevel respectively.


void dSpaceDestroy (dSpaceID);
Destroy a space. When a space is destroyed, all the geometry objects in that space are automatically destroyed as well.


void dSpaceAdd (dSpaceID, dGeomID);
void dSpaceRemove (dSpaceID, dGeomID);
int dSpaceQuery (dSpaceID, dGeomID);
The first two functions add and remove geometry objects to/from the space. They may be called by the geometry object creation/deletion functions. dSpaceQuery() will return 1 if the given geometry object is in the given space, 0 if not. It is provided for convenience.


void dSpaceCollide (dSpaceID space, void *data,
                    dNearCallback *callback);
Call a callback function one or more times, for all potentially intersecting objects in the space. The callback function is of type dNearCallback, which is defined as:
typedef void dNearCallback (void *data, dGeomID o1, dGeomID o2);
The data variable is passed from
dSpaceCollide() to the callback function. Its meaning is user defined. The o1 and o2 arguments are the geometry objects that may be near each other. The callback function can call dCollide() on o1 and o2, perhaps first determining whether to collide them at all based on other information.


10. How to make good simulations

[just notes for now]

10.1. Integrator accuracy and stability

10.2. Behavior may depend on step size

10.3. Making things go faster

What factors does execution speed depend on? Each joint removes a number of degrees of freedom (DOFs) from the system. For example the ball and socket removes three, and the hinge removes five. For each separate group of bodies connected by joints, where:

then the computing time per step for the group is proportional to:
k1 O(m1) + k2 O(m23) + k2 O(n)

ODE currently relies on factorization of a ``system'' matrix that has one row/column for each DOF removed (this is where the O(m23) comes from). In a 10 body chain that uses ball and socket joints, roughly 30-40% of the time is spent filling in this matrix, and 30-40% of the time is spent factorizing it.

Thus, to speed up your simulation you might consider:

In the future ODE will implement techniques that scale better with the number of joints.

10.4. Making things stable

Increasing the global CFM will make the system more numerically robust and less susceptable to stability problems. It will also make the system look more ``spongy'', so a tradeoff has to be found.

Redundant constraints (two or more constraints that ``try and do the same job'') will fight each other and cause stability problems. The numerical cause of this problem is singularity in the system matrix. One example of this is if two contacts joints connect the same pair of bodies at the same point. Another example is if a virtual hinge joint is created between two bodies by connecting them with two ball joints, spaced apart along the hinge axis (this is bad because the two ball joints try to remove six degrees of freedom from the system, but a real hinge joint would only remove five).

Redundant constraints fight each other and generate strange forces in the system that can swamp the normal forces. For example, an affected body might fly around as though it has a life of its own, with complete disregard for gravity.

10.5. Using constraint force mixing (CFM)

10.6. Avoiding singularities

10.7. Other stuff


11. FAQ

11.1. How do I connect a body to the static environment with a joint?

Use dJointAttach with arguments (body,0) or (0,body).

11.2. Does ODE need or use graphics library X ?

No. ODE is a computational engine, and is completely independent of any graphics library. However the examples that come with ODE use OpenGL, and most interesting uses of ODE will need some graphics library to make the simulation visible to the user. But that's your problem.

11.3. Why do my rigid bodies bounce or penetrate on collision? My restitution is zero!

Sometimes when rigid bodies collide without restitution, they appear to inter-penetrate slightly and then get pushed apart so that they only just touch. The problem gets worse as the time step gets larger. What is going on?

The contact joint constraint is only applied after the collision is detected. If a fixed time step is being used, it is likely that the bodies have already penetrated when this happens. The error reduction mechanism will push the bodies apart, but this can take a few time steps (depending on the value of the ERP parameter).

This penetration and pushing apart sometimes makes the bodies look like they are bouncing, although it is completely independent of whether restitution is on or not.

Some other simulators have individual rigid bodies take variable sized timesteps to make sure bodies never penetrate much. However ODE takes fixed size steps, as automatically choosing a non-penetrating step size is problematic for an articulated rigid body simulator (the entire ARB structure must be stepped to account for the first penetration, which may result in very small steps).

There are three fixes for this problem:

11.4. How can an immovable body be created?

In other words, how can you create a body that doesn't move, but that interacts with other bodies? The answer is to create a geometry object only, without the corresponding rigid body object. The geometry object is associated with a rigid body ID of zero. Then in the contact callback when you detect a collision between two geometry objects with a nonzero body ID and a zero body ID, you can simply pass those two IDs to the dJointAttach() function as normal. This will create a contact between the rigid body and the static environment.

Don't try to get the same effect by setting a very high mass/inertia on the ``motionless'' body and then resetting it's position/orientation on each time step. This can cause unexpected simulation errors.

11.5. Why would you ever want to set ERP less than one?

From the definition of the ERP value, it seems than setting it to one is the best approach, because then all joint errors will be fully corrected at each time step. However, ODE uses various approximations in its integrator, so ERP=1 will not usually fix 100% of the joint error. ERP=1 can work in some cases, but it can also result in instability in some systems. In these cases you have the option of reducing ERP to get a better behaving system.

11.6. Is it advisable to set body velocities directly, instead of applying a force or torque?

You should only set body velocities directly if you are setting the system to some initial configuration. If you are setting body velocities every time step (for example from motion capture data) then you are probably abusing your physical model, i.e. forcing the system to do what you want rather than letting it happen naturally.

The preferred method of setting body velocities during the simulation is to use joint motors. They can set body velocities to a desired value in one time step, provided that the force/torque limit is high enough.

11.7. Why, when I set a body's velocity directly, does it come up to speed slower when joined to other bodies?

What is likely happening is that you are setting the velocity of one body without also setting the velocity of the bodies that it is joined to. When you do this, you cause error in the system in subsequent time steps as the bodies come apart at their joints. The error reduction mechanism will eventually correct for this and pull the other bodies along, but it may take a few time steps and it will cause a noticeable "drag" on the original body.

Setting the velocity of a body will affect that body alone. If it is joined to other bodies, you must set the velocity of each one separately (and correctly) to prevent this behavior.

11.8. Should I scale my units to be around 1.0 ?

Say you need to simulate some behavior on the scale of a few millimeters and a few grams. These small lengths and masses will usually work in ODE with no problem. However occasionally you may experience stability problems that are caused by lack of precision in the factorizer. If this is the case, you can try scaling the lengths and masses in your system to be around 0.1..10. The time step should also be be scaled accordingly. The same guideline applies when large lengths and masses are being used.

In general, length and mass values around 0.1..1.0 are better as the factorizer may not lose so much precision. This guideline is especially helpful when single precision is being used.

11.9. I've made a car, but the wheels don't stay on properly!

If you are building a car simulation, typically you create a chassis body and attach four wheel bodies. However, you may discover that when you drive it around the wheels rotate in incorrect directions, as though the joint was somehow becoming ineffective. The problem is observed when the car is moving fast (so the wheels are rotating fast), and the car tries to turn a corner. The wheels appear to rotate off their proper constraints as though the ``axles'' had become bent. If the wheels are rotating slowly, or the turn is made slowly, the problem is less apparent.

The problem is that numerical errors are being caused by the high rotation speed of the wheels. Two functions are provided to fix this problem: dBodySetFiniteRotationMode() and dBodySetFiniteRotationAxis(). The wheel bodies should have their finite rotation mode set, and the wheel's finite rotation axes should be set every time step to match their hinge axes. This will hopefully fix most of the problem.

11.10. How do I make ``one way'' collision interaction

Suppose you need to have two bodies (A and B) collide. The motion of A should affect the motion of B as usual, but B should not influence A at all. This might be necessary, for example, if B is a physically simulated camera in a VR environment. The camera needs collison response so that it doesn't enter into any scene objects by mistake, but the motion of the camera should not affect the simulation. How can this be achieved?

Here is a good solution: when the collision is detected, don't create a contact joint between A and B as you normally would. Instead, attach the contact joint between B and 0 (the static environment). That way the body A will appear to B as though it is static and unmovable. This approach may result in some penetration between A and B, but this will not be a problem in many applications.

11.11. The Windows version of ODE crashes with large systems

ODE requires stack space roughly on the order of O(n)+O(m2), where n is the number of bodies and m is the sum of all the joint constraint dimensions. If m is large, this can be a lot of space!

Unix-like operating systems typically allocate stack space as it is needed, with an upper limit that might be in the hundreds of Mb. Windows compilers normally allocate a much smaller stack. If you experience crashes when running large systems, try increasing the stack size. For example, the MS VC++ command line compiler accepts the /Stack:num flag to set the upper limit.

11.12. My simple rotating bodies are unstable!

If you have a box whose sides have different lengths, and you start it rotating in free space, you should observe that it just tumbles at the same speed forever. But sometimes in ODE the box will gain speed by itself, spinning faster and faster until it ``explodes'' (disappears off to infinity). Here is the explanation:

ODE uses a first order semi-implicit integrator. The ``semi implicit'' means that some forces are calculated as though an implicit integrator is being used, and other forces are calculated as though the integrator is explicit. The constraint forces (applied to bodies to keep the constraints together) are implicit, and the "external" forces (applied by the user, and due to rotational effects) are explicit. Now, inaccuracy in implicit integrators is manifested as a reduction in energy - in other words the integrator damps the system for you. Inaccuracy in explicit integrators has the opposite effect - it increases the system energy. This is why systems simulated with explicit first order integrators can explode.

So, a single body tumbling in space is effectively explicitly integrated. If the body's moments of inertia were equal (e.g. if it is a sphere) then the rotation axis will remain constant, and the integrator error will be small. If the body's moments of inertia are unequal then the rotation axis wobbles as momentum is transferred between different rotation directions. This is the correct physical behavior, but it results in higher integrator error. The integrator in this case is explicit so the error increases the energy, which causes faster and faster rotation, causing more and more error - leading to the explosion. The problem is particularly evident with long thin objects, where the 3 moments of inertia are highly unequal.

To prevent this, do one or more of the following:

In the future I may add a feature to ODE to modify the rotational dynamics of selected bodies so that they exhibit no rotational error with ODEs integrator.

11.13. My rolling bodies (e.g. wheels) sometimes get stuck between geoms

Consider a system where rolling bodies roll over an environment made up of multiple geometry objects. For example, this might be a car driving over a terrain (the rolling bodies are the wheels). If you find that the rolling bodies mysteriously come to a stop when they roll from one geometry object to another, or when they receive multiple contact points, then you may need to use a different contact friction model. This section explains the problem and the solution.

11.13.1. The problem

An example of such a system is shown in the following picture of a ball that has just rolled down a ramp and touched the ground:

Normally, the ball should continue rolling along the ground, towards the right. However, if ODE's default contact friction mode is being used then the ball will come to a complete stop when it hits the ground. Why?

ODE has two ways to approximate friction: the default way (called the constant-force-limit approximation, or ``box friction'') and an improved way (called ``friction pyramid approximation 1'') which is obtained by setting the dContactApprox1 flag in the contact joint's surface mode field.

Consider the above picture. There are two contact points, one between the ball and the ramp, the other between the ball and the ground. If the box friction mode is used in both contacts and the mu parameter is set to dInfinity then the ball can not slip against the ramp or ground at either contact.

If no slip is possible at a ball contact point, then the center of the ball must move along a path that is an arc around the contact point. Thus the center of the ball is required to simultaneously move along the path ``Arc 1'' and ``Arc 2''. The only way to satisfy both paths at once is for the ball to stop moving altogether.

This is not a bug in ODE - so what is going on here? Objects in real life do not get stuck like this. The problem is that, in the simple ``box'' approximation of friction the tangential force available at a contact constraint to stop it slipping is independent of the normal force that prevents penetration. This is not real-life physics, so we should not be surprised that non-real-life motion results.

Note that this problem does not occur if mu is set to zero, but this is not a helpful solution because we need some amount of friction to model the real world.

11.13.2. The solution

The solution is to use the dContactApprox1 flag in the contact's surface mode field, and set mu to some appropriate value between 0 and infinity. This mode ensures that there will only be a tangential anti-slipping force at the contact point if the contact normal force is nonzero. In the above example it turns out that contact-1 will have a zero normal force, so there will be no force applied at contact-1 at all, and the problem is solved! (the ball will roll along the ground properly.)

The dContactApprox1 mode may not be appropriate in all situations, which is why it is optional. It is important to remember that, although it is a better friction approximation, it is not true Coulomb friction. Thus it is still possible that you may encounter some examples of non-physical behavior.


12. Known BUGS


13. ODE internals

[only notes for now]

Why don't I implement a proper friction pyramid or friction cone (e.g. Baraff's version) ?} Because I have to factor non-symmetric (and possibly indefinite) matrices, for either static or dynamic friction. Speed was considered more important - the current friction approximation only needs a symmetric factorization, which is twice as fast.

13.1. Matrix storage conventions

Matrix operations like factorization are expensive, so we must store the data in a way that is most useful to the matrix code. I want to do 4-way SIMD optimizations later, so the format is this: store the matrix by rows, and each row is rounded up to a multiple of 4 elements. The extra "padding" elements at the end of each row/column must be set to 0. This is called the "standard format". Hopefully this decision will remain good in the future, as more and more processors have 4-way SIMD (especially for fast 3D graphics).

The exception: matrices that have only one column or row (vectors), are always stored as consecutive elements in standard row format, i.e. there is no interior padding, only padding at the end.

Thus: all 3x1 floating point vectors are stored as 4x1 vectors: (x,x,x,0).

13.2. Internals FAQ

13.2.1. Why do some structures have a dx prefix and some have a d prefix?

The dx prefix is used for internal structures that should never be visible externally. The d prefix is used for structures that are part of the public interface.

13.2.2. Returned vectors

There seem to be 2 ways of returning vectors in ODE, e.g.:
    const dReal* dBodyGetPosition (dxBodyID);
    void dWorldGetGravity (dxWorldID, dVector3);
Why? The second way is the 'official' way. The first way returns pointers to volatile internal data structures and is less clean API-wise. For a stable API I feel that filling in vectors is cleaner than returning pointers to vectors, for two reasons:

  1. The returned vector values may have to be calculated somehow, so there is no internal ``cache'' to return a pointer to.
  2. The internal data structures may be moved, which is a problem if the user keeps the returned pointer and uses it later.
As it happens these two cases don't currently happen in ODE - most returned vector data is cached and always at the same address. But having the freedom to change things in the future is useful. The current API shouldn't slow you down because the cases where you need to be fast (i.e. getting body transforms) return pointers anyway - breaking my own rule.
Russell Smith, Tuesday 29 October, 2002