🔒 Locked

Joint issues (constraints)

Started by Black-Panther Jul 29, 2009 at 5:13 AM 10 replies 1.9k views

Original Post

Jul 29, 2009 at 5:13 AM

Hi all! Again I have a question! I now implemented a PGS (after Erin Catto) as constraint-solver and modelled all collision through normal and friction constraints. So far, collision response works well --> therefore the solver seems to be implemented correctly. Now I wanted to add some joints. Started with simply once, like only-plane-movement. worked fine. Then I went on to ball-socket-joints. In the special-case the anchor-point is placed at each center of both bodies, it seems to work. But if not: eg for the pendulum-case it behaves very strange. Here some code how I implemented the ball-joint: Code snippet Annotations: - The c'tor of ogpVelConstraint takes the two involved bodies and the 4 vector-entries of the jacobian as argument. - I tried various values for fBias. The basic-behaviour of the sim. did not change - For those who wonder why the forth, and not the second jacobian-entry is negated, this is because of my special definition of the rotationmatrix and quaternions. Point-velocity in my case is computed after: vp = v + r x w To show you how my simulation acts, I uploaded two vids. (usually I simulate at 60Hz --> same problem). The big box is static and is used as anchor-point for the second-little box --> pendulum Without error correction With error correction (fBias == 0.1f here) I'm thankful for any suggestions what might be wrong! Thanks Black-Panther

Dirk Gregorius

Jul 29, 2009 at 7:24 AM

const og3DVector vError = -fBias / dt * ((v1 - v2) - m_vPosDiff);

Should be: const og3DVector vError = -( v1 - v2 ) - fBias * m_vPosDiff / dt

Don't multiplyscale the velocity error with dt and the bias! Signs could be different. I didn't check how you setup your Jacobians. Besides this the code looks good :)

Check whether you makes the same mistake in you conact constraints!

Cheers,
-Dirk

Black-Panther

Jul 29, 2009 at 8:30 AM

Are u sure about that? In my eyes this can't be. Because m_vPosDiff is the position difference between the anchor-points at the beginning. so to say the CORRECT distance between them. And if I want to compute the error between the should-value and the is-value, I've to compute:

const og3DVector vDeltaPos = (v1 - v2) - m_vPosDiff;

Afterwards I take this delta and compute a Xi proportional to it:
const og3DVector vError = fFactor * vDeltaPos;

where factor is nothing else then:
const float fFactor = -fBias / dt;

All together leads to my formula:
const og3DVector vError = -fBias / dt * ((v1 - v2) - m_vPosDiff);

Otherwise I only would scale the initial-distance between the anchor-points.

NOTE: Perhaps you didn't see, that I don't force the anchor-points to be at the same position, but I try to force them to keep their initial-distance to each other!

Dirk Gregorius

Jul 29, 2009 at 9:16 AM

I am 100% sure. Whether the points are coincident or not doesn't matter. You only remove the relative velocity at the anchor points.

The formula goes like this:
J * M^-1 JT * lambda = -0.1 * C / dt - J * v
P = JT * lambda

Check Erin's code and presentations if you don't trust me. Without position correction you iteratively seek impulses that remove the relative velocities at the anchor points. Since you suffer from drift you put back some stabilization term (called Baumgarte stabilization). Using the Baumgarte term you can bring the anchors to what ever distance you like.

Here is some demo I programmed lately if you want to check that this works in praxis (iterate through the demos a bit until you come to the joints):
www.dondickied.de/download/demo.zip

Did you actually check my suggestion before going into a theoretical discussion here :)

Cheers,
-Dirk

Numsgil

Jul 29, 2009 at 10:45 AM

Just to back up Dirk:

Here's the term plus baumgarte from my engine for collisions:

( 1 + getRestitution() ) * ( deltaVel.dot3(getNormal()) ) + depth / deltaTime * baumgarteUnitFraction

where depth is the error term

Dirk Gregorius

Jul 29, 2009 at 10:50 AM

A small tip Numsgil. This way you will gain energy for coefficient of restitution e = 1 since you basically additionally add the bias. You can look at the demo with the 5 balls in my demo above. I use somethink like this:

plVector3 v_rel;
plComputeRelativeVelocityAt( v_rel, pBody1, r1, pBody2, r2 );
plReal vn_rel = plVec3Dot( v_rel, n );

const plReal kMinVelocity = 1.0f;
if ( vn_rel < -kMinVelocity )
{
plReal ve = restitution * vn_rel;
if ( ve < out.mBias )
{
out.mBias = ve;
}
}

Cheers,
-Dirk

Black-Panther

Jul 29, 2009 at 11:18 AM

Sure I tried it... It doesn't work either:
const float fBias = 0.2f;
const og3DVector vError = -(v1 - v2) - (fBias / dt) * m_vPosDiff;
DDDtest

------------------------------------------------------
Just to be sure:
My velocity-constraint looks like this: (' means the first time derivative)
C' = v1 + r1 x w1 - v2 - r2 x w2 = 0

coming from the position-constraint:
C = p1 - p2 - diff = 0
where diff is p1 - p2 at the constraint-initialisation. Therefore the two anchor-points are to keep Length(diff) far apart from each other.

As you wrote the baumgarte stabilization-term is proportional to C:
Xi = -(ERP / dt) * C

So the equation becomes:
J*V = Xi

_________________________________

If this all is true, then my code was correct. Maybe the name fBias is confusing. Call it ERP. And perhaps the name m_vPosDiff too. It's not the current position error. This is the position difference at the beginning of the simulation.

If this is wrong, then I don't understand where I'm wrong. Hopefully some of you could point that out, please.

Thanks a lot

Dirk Gregorius

Jul 29, 2009 at 11:37 AM

Look correct. Did you try it without stabilization. Does it work in this case? Also, how do you compute the effective mass?

Black-Panther

Jul 29, 2009 at 11:58 AM

Yes, I did. See first video.

I don't know what you mean by effective mass. I compute the jacobian and the constraint-solver (PGS) does the rest. The effective mass should be computed there indirectly... Or am I wrong? The joint only uses masses and inertia tensor to compute the B-matrix (like in Cattos paper):

See method computeMatrixB()
ogpVelConstraint

Now the ball-joint looks like this: (seems more than it is)
ogpJoint

Numsgil

Jul 29, 2009 at 2:44 PM

Quote:
Original post by DonDickieD
A small tip Numsgil. This way you will gain energy for coefficient of restitution e = 1 since you basically additionally add the bias. You can look at the demo with the 5 balls in my demo above. I use somethink like this:

plVector3 v_rel;
plComputeRelativeVelocityAt( v_rel, pBody1, r1, pBody2, r2 );
plReal vn_rel = plVec3Dot( v_rel, n );

const plReal kMinVelocity = 1.0f;
if ( vn_rel < -kMinVelocity )
{
plReal ve = restitution * vn_rel;
if ( ve < out.mBias )
{
out.mBias = ve;
}
}

Cheers,
-Dirk

Thanks.

Actually, thinking about it, if I were to work on a "real" "production" level physics engine I'd probably dump Baumgarte (it's basically a hack anyway) and do a three pass constraint solver: first perform a force level constraint solver for resting contacts. Then do a velocity level constraint solver for colliding objects. Then do a position level constraint solver to make things not interpenetrate visually.

To the best of my knowledge I don't think anyone's ever tried that before. it's usually just pick a level (velocity being the most popular from what I can see) and formulate everything in terms of velocity constraints.

Dirk Gregorius

Jul 29, 2009 at 3:44 PM

Solving on the acceleration level can yield infinite forces in presence of friction. E.g. a simple falling rod. The velocity level is not as bad as some people always state here. In the demo above you have a two pass solver if you choose "Projection" as stabilization method. First impulses and then positions to resolve penetrations and joint separation. This improves articulated structures quite a bit in my opinion.

An advantage of solving on the acceleration level is that you know the constraint forces and can use every integrator you like. This solves some other problems.

Numsgil

Jul 29, 2009 at 4:20 PM

Quote:
Original post by DonDickieD
Solving on the acceleration level can yield infinite forces in presence of friction. E.g. a simple falling rod.

I was thinking more in terms of static resting contact forces, or something like a bridge. Let's say you're making a game like bridge builder. If you modeled the bridge as a constrained group of rigid bodies, you could reuse the exact same Jacobian across multiple frames, assuming that you don't allow any bridge deformation (either it holds or it breaks). So you could directly decompose the effective mass matrix when the level first loaded (that's the part that takes n^3 time), and use it to easily determine stress and strain values for some given forces the bridge is supporting in n^2 time (if it's sparse it might be more like linear time), which is fast enough to do on a per frame basis as the train drives over.

In a more general physics engine, it means your normal constraint solver doesn't have to convolute resting forces (which change little per frame and are physically forces, not impulses) with collision responses (which change probably every single frame (if not, they're resting contact forces) and are physically impulses) you should get better/faster results. Especially if you can precompute the effective mass matrix during level load. Or even inside your tool pipeline.

Quote:

The velocity level is not as bad as some people always state here. In the demo above you have a two pass solver if you choose "Projection" as stabilization method. First impulses and then positions to resolve penetrations and joint separation. This improves articulated structures quite a bit in my opinion.

Yes, that was what I was imagining.

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.