Runge Kutta Merson t1>t0 assert

[Jean-Philippe Jodoin]

Hi Sherman, first I want to thank you for the great help you have provided me with my other problems. It is really appreciated. I have a really recurrent assert in RungeKuttaMersonIntegratorRep::attemptAStep. It is this line: assert(t1 > t0); (t1 seems to be equal to t0). In this case, I have a a couple of pendulum and after a couple of seconds the system start to become really slow and this assert will popup. I also wanted to know if you prefer me to post these problem I encounter on the forum, or if you prefer something else. Either case, let me know.

So here is a simple code:

MultibodySystem system;
SimbodyMatterSubsystem matter(system); //Subsystem that contains all the bodies of the system
GeneralForceSubsystem forces(system); //Subsystem that contains all the forces of the system

if(usegravity)
Force::UniformGravity gravity(forces, matter, Vec3(0, -9.81, 0)); //We add a gravity force acting on all the body of the system


Body::Rigid BetaBox(MassProperties(35000.0, Vec3(0), Inertia(80000)));
BetaBox.addDecoration(Transform(), DecorativeBrick().setOpacity(0.5));

SimTK::Rotation Rot(SimTK::Quaternion(1,0,0,0));
SimTK::Rotation Rot45(SimTK::Quaternion(0.923,0,0.382,0));
MobilizedBody::Ball box1(matter.Ground(), Transform(Vec3(0,5,0)), BetaBox, Transform(Rot,Vec3(3, 8, 0)));
MobilizedBody::Ball box2(matter.Ground(), Transform(Vec3(0,5,0)), BetaBox, Transform(Rot45, Vec3(3, -8, 0)));
MobilizedBody::Ball box3(box1, Transform(Vec3(0,0,0)), BetaBox, Transform(Vec3(3, 3, 0))); //Transform(Rot2, Vec3(-7, -7, -7)).invert())
MobilizedBody::Ball box4(box3, Transform(Vec3(0,0,0)), BetaBox, Transform(Vec3(4, 4, 0)));
MobilizedBody::Ball box5(box4, Transform(Vec3(0,0,0)), BetaBox, Transform(Vec3(3, 3, 0)));
SimTK::Constraint::ConstantAngle(box1, UnitVec3(1,0,0), box2, UnitVec3(0,1,0), -Pi);

system.updDefaultSubsystem().addEventReporter( new VTKEventReporter(system, 0.01));

// Initialize the system and state.
system.realizeTopology();
State state = system.getDefaultState();
// Simulate it.
RungeKuttaMersonIntegrator integ(system);
TimeStepper ts(system, integ);
ts.initialize(state);
ts.stepTo(25);

[Michael Sherman]

Hi, Jean-Philippe. Sorry for the delay -- it's grant-writing season here.

The problem is that you are using the ConstantAngle constraint outside its usable range -- it cannot be used efficiently at angles very close to 0 or +/-180, and you have it set to -Pi. The reason is that this is only a single constraint equation, but aligning two axes requires *two* constraint equations (because there is only one rotation left).

Your system works fine at -Pi/2, and even works at -.99*Pi although it is noticeable slower.

I wish I could say "you should have read the documentation". In fact I found this comment:

// We would like to enforce the constraint that cos(theta) is a constant. This can be done
// with a single constraint equation as long as theta is sufficiently far away from 0 and
// 180, with the numerically best performance at theta=90 degrees where cos(theta)==0.
//
// If you want to enforce that two axes are aligned with one another (that is, the angle
// between them is 0 or 180), that takes *two* constraint equations since the only remaining
// rotation is about the common axis.

But it was in the hidden implementation code rather than in the API -- my bad! I've fixed the API comment in Simbody 1.7 (current svn trunk).

Best regards,
Sherm

[Michael Sherman]

Here is the newly checked-in documentation for the ConstantAngle API:
/**
* This constraint consists of a single constraint equation that enforces that
* a unit vector v1 fixed to one body (the "base body") must maintain a fixed
* angle theta with respect to a unit vector v2 fixed on the other body (the
* "follower body"). This can be done with a single constraint equation as long
* as theta is sufficiently far away from 0 and +/-Pi (180 degrees), with the
* numerically best performance at theta=Pi/2 (90 degrees).
*
* @warning
* Do not use this constraint to \e align the vectors, that is for angles near
* 0 or +/- Pi; performance will noticeably degrade within a few degrees of
* these limits and numerical integration will eventually fail at the limits.
*
* If you want to enforce that two axes are aligned with one another (that
* is, the angle between them is 0 or +/-Pi), that takes \e two constraint
* equations since the only remaining rotation is about the common axis. (That
* is, two rotational degrees of freedom are removed; that can't be done with
* one constraint equation -- the situation is analogous to the inability of
* a Rod (distance) constraint to keep to points at 0 distance.) Instead,
* you can use two ConstantAngle constraints on pairs of vectors perpendicular
* to the aligned ones, so that each ConstantAngle is set to the optimal 90 degrees.
*
* This constraint is enforced by an internal scalar torque applied equal and
* opposite on each body, about the mutual perpendicular to the two vectors.
*
* The assembly condition is the same as the run-time constraint: the
* bodies must be rotated until the vectors have the right angle between them.
*/

[Michael Sherman]

er, that's "two" points, not "to" points

[Jean-Philippe Jodoin]

Thanks a lot, that clarify it. I'll update my Simbody SVN trunk to get the extra comments.


Jean-Philippe Jodoin