Rotation parameters do not account for large rotations
It seems that the rotation parameters are computed from the rotation matrix by assuming only small rotations, because of the incremental rotation handling in MBDyn. However, if we are using automatic differentiation and line search based nonlinear solvers for example, we could use larger time steps with larger changes in the rotation matrix during the nonlinear solution phase. This could violate the assumption of small rotations.
class Param_Manip : public Vec3_Manip {
public:
Param_Manip() {};
inline void Manipulate(Vec3& v, const Mat3x3& m) const {
// singularity test
doublereal d = 1. + m.Trace();
if (fabs(d) < std::numeric_limits<doublereal>::epsilon()) {
silent_cerr("Param_Manip(): divide by zero, "
"probably due to singularity in rotation parameters" << std::endl);
throw ErrDivideByZero(MBDYN_EXCEPT_ARGS);
}
v = m.Ax()*(4./d);
};
};
So, I would like to propose a new algorithm for the calculation of the rotation parameters in order to eliminate this limitation of small rotations. See the attached file: rotation_matrix_to_rotation_param.m