Runge Kutta Merson t1>t0 assert

Simbody is useful for internal coordinate and coarse grained molecule modeling, large scale mechanical models like skeletons, and anything else that can be modeled as bodies interconnected by joints, acted upon by forces, and restricted by constraints.
POST REPLY
User avatar
Jean-Philippe Jodoin
Posts: 14
Joined: Thu Sep 10, 2009 6:32 am

Runge Kutta Merson t1>t0 assert

Post by Jean-Philippe Jodoin » Thu Oct 01, 2009 6:59 am

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);

User avatar
Michael Sherman
Posts: 807
Joined: Fri Apr 01, 2005 6:05 pm

RE: Runge Kutta Merson t1>t0 assert

Post by Michael Sherman » Thu Oct 01, 2009 9:02 pm

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

User avatar
Michael Sherman
Posts: 807
Joined: Fri Apr 01, 2005 6:05 pm

RE: Runge Kutta Merson t1>t0 assert

Post by Michael Sherman » Thu Oct 01, 2009 9:35 pm

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.
*/

User avatar
Michael Sherman
Posts: 807
Joined: Fri Apr 01, 2005 6:05 pm

RE: Runge Kutta Merson t1>t0 assert

Post by Michael Sherman » Thu Oct 01, 2009 9:39 pm

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

User avatar
Jean-Philippe Jodoin
Posts: 14
Joined: Thu Sep 10, 2009 6:32 am

RE: Runge Kutta Merson t1>t0 assert

Post by Jean-Philippe Jodoin » Fri Oct 02, 2009 4:43 am

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


Jean-Philippe Jodoin

POST REPLY