Original Post
Thanks in advance. I'm studying about implicit integration method. I can understand how it works, however, I failed to implement it. I've coded for only one spring for two pointmasses. But it failed. Just only an identity jacobian matrix is calculated. Anyone has an example or code segments for implicit spring? Here's my code, what's wrong with it?
ComputeJacobian()
{
Vector3f oD = m_poPM1->GetPoint() - m_poPM2->GetPoint();
Matrix3f oDDT(
oD.x*oD.x, oD.x*oD.y, oD.x*oD.z,
oD.y*oD.x, oD.y*oD.y, oD.y*oD.z,
oD.z*oD.x, oD.z*oD.y, oD.z*oD.z );
float fLengthInv = oD.Length();
if( !IsZero( fLengthInv ) )
fLengthInv = 1.0f/fLengthInv;
oDDT = oDDT * fLengthInv * fLengthInv;
m_oJX = m_fSpringConst * ( oDDT + ( Matrix3f::IDENTITY - oDDT )*( 1.0f - m_fOrgLength*fLengthInv ) );
m_oJV = Matrix3f::IDENTITY * m_fDampingConst;
}
Update( float fDT )
{
// -------------------------------------------------------------------------
// implicit integration( backward euler method )
// -------------------------------------------------------------------------
ComputeJacobian();
Matrix3f M = Matrix3f::IDENTITY * m_poPM1->GetMass();
Matrix3f A = M - fDT*m_oJV - fDT*fDT*m_oJX;
m_poPM1->SetForce( m_poPM1->GetForce() + fDT*m_oJX*m_poPM1->GetVelocity() );
Vector3f b = fDT * m_poPM1->GetForce();
if( A.MakeInverse() == false ) // not invertible
return;
Vector3f oDV = A * b;
Vector3f oDX = fDT * ( m_poPM1->GetVelocity() + oDV );
m_poPM1->SetVelocity( m_poPM1->GetVelocity() + oDV );
m_poPM1->SetPoint( m_poPM1->GetPoint() + oDX );
m_poPM2->SetVelocity( m_poPM2->GetVelocity() - oDV );
m_poPM2->SetPoint( m_poPM2->GetPoint() - oDX );
}