martes, 20 de abril de 2010

An Overview on Constraint Enforced Formulations of Variational Dynamics

A couple of weeks ago a very interesting paper entered my hard drive: Review of classical approaches for constraint enforcement in multibody systems.
It is a fairly clarifying overview on solid grounds of certain methodologies to solve multibody systems dynamics.
These methodologies differ from the ones I have encountered in the spreaded dynamics engines in their apparently more robust and simple formulation.
Here is a synthesis of this paper:

Multibody systems present two distinguishable features:
  1. Bodies undergo finite relative rotations, which introduce nonlinearities
  2. Bodies are connected by mechanical joints that impose restrictions, which mean a set of governing equations that combine differential and algebraic equations (ODEs and DAEs, respectively).
Lagrange's equation of the first kind has the following aspect:
  • M=M(q,t) is the mass matrix
  • q is the generalized coordinates vector
  • B is the constraint matrix
  • λ is the array of Lagrange multipliers
  • F is the dynamic externally applied forces
If all the constraints are holonomic (velocity independent), B is called the Jacobian matrix, and the generalized coordinates, q, are linked by m algebraic constraints.
Lagrange's equations of the first kind form a set of (m+n) Differential Algebraic Equations.
The approach of the following methods is to use algebraic procedures in order to eliminate Lagrange's multipliers and then obtain a set of ODEs.
  • Maggi's formulation: implies the creation of a vector containing the so called kinematic parameters, generalized speeds or independent quasi-velocities by the analyst in order to obtain a Γ matrix. This matrix spans the null space of the constraint matrix and allows for the elimination of the Lagrange multipliers.
  • Index-1 formulation: requires that initial condition of the problem be subjected to the constraint conditions. Then, it is possible to obtain a system of 2nd order ODEs that is solvable by rearranging the previous equation, extracting the Lagrange multipliers and hence obtaining:
  • Null space formulation: this method solves the system of second order ODEs by premultiplying the first part of the equation by the transposed null space matrix thus eliminating the Lagrange multipliers.
  • Udwadia-Kalaba formulation: this method represents a more compact and general form of solving the DAEs by means of the Moore-Penrose generalized inverse. It is based on Gauss' Principle of Minimum Constraint, which establishes that the explicit equations of motion be expressed as the solution of a quadratic minimization problem subjected to constraints, but at the acceleration level.

All these formulations transform the (2n+m) first order DAEs into ODEs by eliminating Lagrange multipliers.
Maggi's formulation yields (2n-m) first order ODEs.
Index-1, null space and Udwadia-Kalaba form sets of (n) second order ODEs which could be alternatively recast into (2n) first order ODEs for the n generalized coordinates and the n generalized velocities.
The main advantage of these methodologies is not so much the reduction in the number of equations but rather in the change from DAEs to ODEs.
There is a warning on the constraint drift phenomenon, for which these method will be more affected (not so much in Maggi's formulation), and that would require constraint stabilization techniques.
My conclusion is that it seems the way to go, not only for the claims of more stable and quick numerical techniques available to get them working, but also for an apparently more clear approach in the theoretical field.
Particularly, I have done some research into the Udwadia-Kalaba formulation, and definitively is a very promising one.

jueves, 15 de abril de 2010


After seeing my own simulator running there is a whole lot of things that I wish it had.
All this stuff is more or less complicated to achieve, but from my previous research I have found some Open Source engines from which I can learn a lot.

The best documented one so far is ODE (, and also the one I intend to merge with Blender.
So, I present a very quick and thorough overview on how this works, in order to have some scope.

The first thing to disentangle is how it gets working. For that, just looking at the user's manual one gets the following algorithm:

1.-Create a Dynamics world
  • This basically means to make an instance of the dxWorld struct.
2.-Create the bodies
  • Attach the bodies to the world (this means adding data to the dxBody* array of the dynamics world).
  • Set their properties (position and orientation of point of reference, its linear and angular velocities, the mass of the bodies and some other stuff of the like).
3.-Create the joints
  • Attach the joints to the world (by adding data to the dxJoints* array of the dynamics world)
  • Set their properties (depending of the selected type of joint, one has to provide different details).
4.-Manage collisions
  • Create a new collision world (just by making an instance to the dxSpace struct).
  • Create a joint group where collisions will be stored temporarily for every frame step.
  • Apply forces to bodies.
  • Adjust joint parameters.
  • Call collision detection.
  • Create a contact joint for every detected collision point and add it to the collision joint group.
  • Take a simulation step.
  • Clear the collision joint group.
6.-Destroy the dynamics and the collision worlds (wow, that sounds evil...hehehe).

All of this seems fairly easy to do but then one has to get to know where and how things have been implemented. This is when things become a bit harder.
I am sure the file structure and the class definitions make perfect sense for the programmers of this engine.
Fortunately, they have been careful enough so as to comment everything in an understandable way, for which I feel deeply grateful.
After some digging, I have managed to restructure the ode/src folder into the following topics:
  • Accessories
    • Memory management
    • Math
    • Matrix handling
    • External applications
  • Collision
  • Core
  • Joints
  • Primitives
  • Solver
Which means one can tackle the engine in an ordered manner and find things when needed.

From a theoretical point of view, this engine presents the following features:
  • A Lagrange multiplier velocity based model from Trinkle and Stewart
  • A friction model from Baraff
  • A Danzig LCP solver

martes, 13 de abril de 2010

First Ogre+Verlet+Gauss-Seidel simulation

Once the Ogre3D engine is ready to draw what we numerically compute, it was about time to start having some fun.

I found this excellent article from Thomas Jakobsen: where a very neat and simple engine is implemented (thanks Mr. Jakobsen).

It is very well explained so adapting it to my little lab was not very hard.
Here I proudly present my very first simulation!

I am perfectly conscious that it lacks of a lot of things...namely collision detection, a proper stable, precise integrator, an optimized implementation...and a long etcetera...but it's just a baby!

Here is a zip file with the code (an extension from what I presented in the last post):

domingo, 11 de abril de 2010

Setting Ogre3D as Graphic Environment

The last days I have been working on getting a Graphics Environment for my Toy Physics Engine.
Last December and January I had already fiddled around with excellent engines Ogre ( and Bullet (, trying to put them to work together.I have a post from february where some simulations on youtube can be seen and some preliminary introductions are made.

However, now I am starting to understand things, and whenever I have the drive to test them always encounter the problem on how visualize them.

So, I came across with Blender (, another fantastic tool, and got surprised to find their Google Summer of Code mentorship for 2010 (
I have applied for it. Deadline was 9th April, which meant preparing everything in a hurry and some embarrasing trouble with their wiki page (I have already apologised in their bf-commiters mailing list...ooops).
Here is my proposal for such an important event: (
Hopefully they find it interesting and contact me to complete it!

In the meantime, I have prepared a short and basic Ogre framework to begin doing some tests.
It is composed of a few files and classes integrating Ogre.

I generated it basically copying, cutting and pasting from the very good tutorials available at the Ogre web page (

It goes as follows:

The first step is to intall and cofigure Ogre in your computer. This is the best tutorial I found:

Then, using whatever IDE of your taste (I am on EasyEclipse), just get the following files to compile.

The core of the application is the BaseApp object, from which I will later extend further as needed. It is basically in charge of doing everything Ogre requires to be done to put things up in a rendering window.


The BaseFrameListener is another very important piece of code that allows one to control the rendering flow.
In order to manage whatever happens on every rendering step, Ogre provides the FrameListener object.
With the installed package comes an ExampleFrameListener.h file, which overrides the basic Ogre FrameListener and allows for a fairly good control over the inputting devices such as keyboard and mouse.I have slightly customised it for my purposes.
It has been a bit tricky however, and has made my computer crash a few times.
Anyway, here it is, domesticated at last!


Once the beast is under control, the rest is a piece of cake...hehehe... 
You can make as many instances of a FrameListener as you need in an Ogre application.
I have chosen to create a separate framelistener for managing whatever goes under the physics, to keep things up neat and tidy.
It is the PhysicsFrameListener class:


The elementary brick on which everything else will be constructed is likely to be called PhysicsEntity.
So I have defined my PhysicsEntity Class.
In this code is still at its very bones, just presenting a couple of functions to make things visible.
Hopefully in the near future some more flesh will get added on it!



Last, but not least, is the main.cpp file, where everything gets blended:


And this is it. Now...Let's get busy!