* if |x| exceeds limit, set x to limit with the sign of x
*/
if(fabs(x)>limit){// limit delVdotDelR to a fraction of speed of sound
x=limit*copysign(1,x);
}
}
/*
* deviator of a tensor
*/
staticinlineMatrix3dDeviator(constMatrix3dM){
Matrix3deye;
eye.setIdentity();
eye*=M.trace()/3.0;
returnM-eye;
}
/*
* Polar Decomposition M = R * T
* where R is a rotation and T a pure translation/stretch matrix.
*
* The decomposition is achieved using SVD, i.e. M = U S V^T,
* where U = R V and S is diagonal.
*
*
* For any physically admissible deformation gradient, the determinant of R must equal +1.
* However, scenerios can arise, where the particles interpenetrate and cause inversion, leading to a determinant of R equal to -1.
* In this case, the inversion direction is heuristically identified with the eigenvector of the smallest entry of S, which should work for most cases.
* The sign of this corresponding eigenvalue is flipped, the original matrix M is recomputed using the flipped S, and the rotation and translation matrices are
* obtained again from an SVD. The rotation should proper now, i.e., det(R) = +1.