Particle constraint Jacobian derivative

[hop]

I have been following through Witkin's "An Introduction to Physically Based Modelling: Constrained Dynamics" http://www.cs.cmu.edu/~baraff/sigcourse/notesf.pdf so that I fully understand constraint formulation and solution for a particle system before attempting to apply it to a rigid body system.

I have followed through to the end of section 3 and I am happy with the derivation of the matrix equation for the unknown Lagrange multipliers (page F6 eqn 10) :

J * M^-1 * J^T *lamda = -Jdot * qdot - J * M^1 * Q

where if q is the state vector for all the particles, then qdot is the time derivative of the state vector of the particles, Q is the vector of external forces, J = del_C / del_q, and Jdot = del_Cdot / del_q

However when I try to apply this theory to the a simple test case of a particle constrained to a circle on paper I quickly realise that it is straight forward to find J (by inspection) but I cannot see an easy way to determine Jdot for this system, or indeed for any other (e.g. 2-particle distance constraint).

Looking at the equivalent rigid body matrix equation it seems that Jdot is nowhere to be seen! Have I perhaps mis-interpreted its meaning? Thanks for any help.

Following through Cattos derivation form http://www.continuousphysics.com/ftp/pu ... namics.pdf Section 6 I get the expression:

J * M^-1 * J^T *lamda = -J * qdot / dt - J * M^1 * Q

If this is equivalent then this implies that Jdot = del_Cdot / del_t rather than Jdot = del_Cdot / del_q (as in the first derivation. Can anyone offer any insight into this?

[Dirk Gregorius]

Erin solves on the so called "velocity" level while Witkin solves on the "acceleration" level.

Erin uses Newton in impulse form and combines the position and velocity constraint in a Baumgarte fashion (see the paper you mentioned). I assume that external forces are already integrated here to define the sequential impulse scheme.

M * ( v2 - v1 ) = JT * lambda
J * v2 + beta * C / dt = 0 (where dC/dt = J * v)

-> J * W * JT * lambda = -beta * C / dt - J * v1


Witking builds the acceleration constraint by building the time derivative of the *velocity* constraint:

dC/dt = J * v -> d2C/dt^2 = dJ/dt * v + J * dv/dt = 0

You can plug this into the Newton-Euler equation and get the Witking forumulas:

dx/dt = v
dv/dt = M^-1 * ( f_ext + JT * lambda )
dJ/dt * v + J * dv/dt = 0

Usually you don't need to build the time-derivative of the Jacobian which becomes a rank 3 tensor (basically an array of Hessian matrices). So you identify the dJ/dt * v by inspection. The Baraff papers explain this stuff quite in detail, but I would call this deprecated and even modern engineering software solves on the velocity level.
The acceleration formulation is pretty uncommon nowadays. See e.g. the PhD of Anitescu.


HTH,
-Dirk

[hop]

Thanks.

When you write M * ( v2 - v1 ) = JT * lambda
did you mean M * ( v2 - v1 ) / dt = JT * lambda ?

which gives the correct units of acceleration for Newton's 2nd.

-> J * W * JT * lambda = -beta * C / dt - J * v1 / dt

It looks like Witkin's derivation with Catto's expression for acceleration gives the matrix equation I expect.

Regarding the Baumgarte term: I have found that using this term instead of a separate Box2d style position correction post velocity constraint solve gives less pleasing results. i.e. a rope of these distance constraints swings less naturally.

Code: Select all

	m_invMass = m_pParticle[ 0 ]->GetInverseMass() + m_pParticle[ 1 ]->GetInverseMass();
	m_mass = 1.0f / m_invMass;
	const Vector d = m_pParticle[ 1 ]->GetPos() - m_pParticle[ 0 ]->GetPos();
	m_jacobian = Normalise( d );

       // loop N times....
	const Vector dv = m_pParticle[ 1 ]->GetVel() - m_pParticle[ 0 ]->GetVel();
	// definition of velocity constraint cDot = Jv
	const float cDot = Dot( m_jacobian, dv );

	// calculate the Lagrange Multiplier lambda i.e. the constraint space impulse magnitude
	const float beta = 0.0f;
	const float lambdaDt = ( -cDot - beta * c ) / m_invMass * Dot( m_jacobian, m_jacobian );

	// F = J^T * lambda   therefore  J^T * lambda * dt = F * dt = impulse
	const Vector impulse = m_jacobian * lambdaDt;

	m_pParticle[ 0 ]->ApplyImpulse( -impulse );
	m_pParticle[ 1 ]->ApplyImpulse( impulse );

       // end loop


       // position correction loop N times - don't do this if beta > 0

	// position constraint is error in length
	const float c = Length( d ) - m_length;

	// calculate the Lagrange Multiplier lambda
	const float lambda = -c / m_invMass;

	// calculate the corrective impulse
	Vector impulse = m_jacobian * lambda;

	// under-relax the impulse
	impulse *= s_positionCorrectionRelaxationFactor;

	// apply instantaneous (impulse-like) change in positions to the particles
	const Vector dx0 = impulse * m_pParticle[ 0 ]->GetInverseMass();
	const Vector dx1 = impulse * m_pParticle[ 1 ]->GetInverseMass();

	ASSERT( LengthSqr( dx0 ) < 100.0f * 100.0f, "Very large position correction attempted on particle 0. Solution is unstable." );
	ASSERT( LengthSqr( dx1 ) < 100.0f * 100.0f, "Very large position correction attempted on particle 1. Solution is unstable." );

	m_pParticle[ 0 ]->SetPos( m_pParticle[ 0 ]->GetPos() - dx0 );
	m_pParticle[ 1 ]->SetPos( m_pParticle[ 1 ]->GetPos() + dx1 );

        // end position loop

I'm not sure I see this, or indeed what the justification for a position correction of this form is.

[Dirk Gregorius]

Nope I didn't meant this. In my equation lambda is an impulse. Anyway, the difference is only subtle.

Yeah, the unnatural swinging motion is indeed the reason for the post-stabilization. Baumgarte stabilization is energetic. Mathematically Baumgarte stabilizes a differential-algebraic system (DAE) while post-stabilization stabilizes the ODE. The later is more pleasing. For more details on this topic I suggest reading the papers of Ury Asher on constraint stabilization, but they are mathematically quite heavy.

[fishboy82]

Hi:
could you give us some link about Uri Asher's paper on Post Stabilize,Dirk. I googled for them but couldn't find any

[hop]

http://www.cs.rutgers.edu/~dpai/papers/ClinePai03.pdf compares Baumgarte stabilization to post-stabilization and references 2 Ascher papers. I presume these are the relevant documents.

[fishboy82]

Thanks a lot!