Equation of Motion

[Youne]

Hi everyone!

I'm new in bullet physics engine. I hope you will be able to help me because I am a mechanical engineer and I'm not a computer programmer.

In my model I have to calculate the force applied on each object for each step simulation. In a simple model which I made I had to manually calculate the equation of motion and then implement it in the code. I was able to do this, however when the model is made up of several objects linked one to the other with constraintes it's more complicated to extract the value of the force applied to a given object.

How can I extract the value of the force (or impulse) and with which function?
Where in the code can I find the equation of motion?

Thanks for your help!

Youn

[Erwin Coumans]

For unconstrained motion, Bullet uses the Euler integrator:
Given position P, linear velocity V and delta time, new position P' is:
P' = P + V * dt;
Orientation is updated using either quaternion derivative or exponential map)

See implementation in integrateTransform method in Bullet/src/LinearMath/btTransformUtil/btTransformUtil.h. This is called from btDiscreteDynamicsWorld in void btDiscreteDynamicsWorld::integrateTransforms(btScalar timeStep).

For linked multibody constraint motion, Bullet uses pairwise constraints (full coordinates method). All constraints are derived from btTypedConstraint (see Bullet/src/BulletDynamics/ConstraintSolver/btTypedConstraint.h). You can use constraint->getAppliedImpulse(); to get the applied constraint impulse after calling dynamicsWorld->stepSimulation (...).

Hope this helps,
Erwin

[Youne]

Thanks a lot for your replay
I have one more question !

That is BULLET makes a calculation of forces? it means that each time we have a possibility of asking the value of the applied force at a given point of a given object? (like Mechanical calculation code "ANSYS.Inc, SOLIDWORKS")

Thanks

Y

[Dirk Gregorius]

Bullet solves the following constrained systems which builds a system of differential and algebraic equations (DAE).

dv/dt = M^-1 * (f_ext + f_c)
C(x) = 0

The geometric constraints result in constraint forces f_c that act onto the system and can be computed using Lagrange multipliers f_c = JT * lambda where JT is the transposed Jacobian matrix J = dC/dx. The classical mechanical way to solve such a system is two build the second time derivative of the constraint equation and plug those into the equations of motion. This allows you to solve for lambda and the constraint forces. Once you have the constraint forces you can use any integrator you want. This approach suffers from infinite friction forces though. If you want a sorrow explanation of this approach you might look into the work of Baraff.

Bullet (and all other physic libraries for games) go a slightly different approach. First the system is integrated forward in time without considering the constraint forces. Then in a second phase we correct velocities to satisfy each individual constraint in an iterative loop using sequential impulses. Since correcting one constraint may violate another one we need to do this for several iterations. This is called relaxation and is similar to the Gauss-Seidel algorithm for solving systems of linear equations. Once the velocities are within some threshold or if the maximum iteration count is exceeded we update the position using the corrected velocities. Mathematically Bullet solves the following system similar to the work of Anitescu:

We approximate the equations of motion using the differential quotient:

(1) dv/dt ~= ( v(t+dt) - v(t) ) / dt = M^-1 * ( f_ext + JT * lambda )


The velocity constraints can be found by building the time derivative of position constraints C(x)

(2) dC/dt = J * v = 0


If you combine both equations you end up with ( solve (1) for v(t+dt) and plug it into (2) ):

J * W * JT * lamdba * dt = -J * ( v(t) + M^-1 * f_ext * dt )

If you look at the right hand side you notice that v(t) + M^-1 * f_ext * dt is the new velocity without considering the constraints. The described algorithm above origins in solving this system using an iterative projected Gauss-Seidel algorith for a Mixed Linear Complementaty Problem (MLCP). For equality constraints only the above equations would build a system of linear equations. In the context of inequatlity constraints like for contacts, friction, limits or motors you end up with a with the MLCP.

W.r.t. to you final question:
Yes, you can query for the constraint forces, but notice that due to the iterative nature of the solver you get only an approximation. Also note that the contact and constraint forces might only be correct in the sense that the system will behave dynamically correct. E.g. if you have a table standing on the floor the solver should give you a reaction force f = mg/4 for each leg. Instead it might return f = mg/2 for two diagonal legs and f = 0 for the other two. F1 = 1/8 mg and F2 = 3/8 mg are also possible for diagonal pairs. So in this sense Bullet is not comparable to programms from mechanical engineering like e.g. Ansys or Solidworks. Personally I am constructional engineer and I wouldn't use any of the game physic libraries (e.g. PhysX, ODE or Bullet) for mechanical engineering problems. Note that Bullet was designed for realtime simulation like found in games or maybe also movies and not so much with mechanical engineering as focus.


HTH,
-Dirk

[Youne]

Hi Dirk !

I would like to sincerely thank you as you answered my question straight away whilst it took me two weeks of reflection and trying to understand the code and i was still not able to find the answer. It's exacly what I wanted to know. Now my objective is to extract the value of the applied force which will act on each rigidbody. For the moment, I understood very well the equation of motion which they use and the method they use to solve them. But I can't find in which part of the code the applied force is calculated?

Thank you

Youne

[Erwin Coumans]

Now my objective is to extract the value of the applied force which will act on each rigidbody. For the moment, I understood very well the equation of motion which they use and the method they use to solve them. But I can't find in which part of the code the applied force is calculated?

Bullet operates on the velocity level, so it is an impulse that is calculated and applied.

See 'resolveSingleCollisionCombinedCacheFriendly' called from called from btSequentialImpulseConstraintSolver::solveGroupCacheFriendlyIterations in bullet\src\BulletDynamics\ConstraintSolver\btSequentialImpulseConstraintSolver.cpp.

This applied impulse is temporarily data, so it is not available to the user right now, unless you disable 'SOLVER_CACHE_FRIENDLY' in the solver mode. In that case, the applied impulse is stored in 'btConstraintPersistentData::m_appliedImpulse' , as part of the contact point information. You can check in btCollisionInterfaceDemo how to iterate over contact points, to get to this data.

The applied impulse will be easier available in Bullet 2.67, within a week or so.

Hope this helps,
Erwin

[Dirk Gregorius]

In order to satisfy a velocity constraint you apply incremental impulses. Since you repeat this in several iterations let's call these impulses a "delta impulse". You can accumulate these delta impulse and their sum gives you the total impulse. Since in Bullet forces and impulses are linear related F = P / dt this total impulse is you constraint force. A generic example for a bilateral constraint goes something like this:

// Compute "delta impulse"
float dv = J * v
float dP = -dv / (J*W*JT)

// Apply impulse to bodies...

// Accumulate
mAccumulatedImpulse += dP;


So all you need to do is to store the accumulated impulse with your joints and contacts and give them some accessor. The accumulated impulse is basically your constraint force. Note that this accumulated impulse can also be used to warmstart the iterative solver, but in presense of Baumgarte stabilization you get stability issues with unbalanced torques then. This information in only for completeness. To my knowledge Bullet doesn't warmstart at the moment.

[Erwin Coumans]

To my knowledge Bullet doesn't warmstart at the moment.

Actually Bullet provides the option for warmstarting. Just use

Code: Select all

constraintSolver->setSolverMode(btSequentialImpulseConstraintSolver::SOLVER_USE_WARMSTARTING);

The funny thing is that it is disabled by accident in latest Bullet 2.66, and noone noticed It will be enabled again by default in Bullet 2.67, this week.

Thanks,
Erwin

[Youne]

Hi Erwin,

I just download the new version of bullet 2.67, and i try to see the modification on SequentialImpulseConstraintolver but I have not noticed any change to allow access to impulses applied.

Can you explain to me, How i have to do to extract the impulses applied?

thank you very much

[Erwin Coumans]

There are two cases, contact constraints and other constraints derived from btTypedConstraint (btHingeConstraint, btGeneric6DofConstraint, btPoint2PointConstraint, btConeTwistConstraint).

For contacts, you can read the applied impulse in the new member of btManifoldPoint called m_appliedImpulse.
For your information: the following code fragment inside the btSequentialImpulseConstraintSolver (around line 840) copies the applied impulse after all iterations:

Code: Select all

{
		int numPoolConstraints = m_tmpSolverConstraintPool.size();
		int j;
		for (j=0;j<numPoolConstraints;j++)
		{
			const btSolverConstraint& solveManifold = m_tmpSolverConstraintPool[j];
			btManifoldPoint* pt = (btManifoldPoint*) solveManifold.m_originalContactPoint;
			btAssert(pt);
			pt->m_appliedImpulse = solveManifold.m_appliedImpulse;
			//do a callback here?
		}
	}

For all btTypedConstraint-derived constraints, you can get the applied impulse using getAppliedImpulse().

Hope this helps,
Erwin