Page 1 of 2

Prescribed motion

Posted: Tue Oct 27, 2009 9:17 am
by mpoursina
Dear,
Within my research, I want to apply the prescribed motion on the inward and outward boundaries of RNA.
I have started with a single pendulum problem and tried to use
Vector_<Real> coefficients(3);
coefficients[0] = 0.0;
coefficients[1] = 1.0;
coefficients[2] = 0.0;
const Function::Polynomial func3(coefficients);
Constraint::PrescribedMotion(matter, &func3, MobilizedBodyIndex(0), MobilizerQIndex(0));

It compiles fine without any error. But at
ts.initialize(state);
it gives the following message:
SimTK::Integrator::InitializationFailed at memory location 0x0018f640
and stops the simulation.
It would be grateful to have your comments on this problem.
My second question is that, if I want to apply motion or force at the boundaries which one do you think can be implemented easier?
I thank again in advance for your time.
Best regards
Mohammad Poursina



RE: Prescribed motion

Posted: Tue Oct 27, 2009 9:33 am
by sherm
Hi, Mohammad.

In your example above it appears that you applied the PrescribedMotion constraint to Ground (MobilizedBodyIndex(0)) -- but mobilities go with the moving (child) body, not the parent. You need to specify the MobilizedBodyIndex of a body that has some mobility; that index would have been returned at the time you added the body to the system.

Note that prescribed motion is a constraint that applies *only* to mobilities (that is, joint degrees of freedom). If you want to use prescribed motion you must first build your model so that it has mobilities corresponding to the motion you want to control. That means if you want to use prescribed motion at both ends of a chain, both the first and last bodies would have to have joints back to ground; somewhere in the middle you would have to split a body or a joint and add a Weld or Ball constraint.

It is always easier to add forces because they can be applied to any model anywhere. However, if what you are actually trying to control is motion, it is not always easy to find the right forces. Using prescribed motion for that purpose can be very convenient but you do have to plan your model to enable that.

Regards,
Sherm

RE: Prescribed motion

Posted: Tue Oct 27, 2009 9:58 am
by mpoursina
Dear Sherm,

Thanks for your reply. I implemented your comment in the code using:
MobilizedBodyIndex(1). The same message still appears. When I follow it, it seems that the problem is raised from the following line:
getSystem().realizeModel(updAdvancedState());
in
void IntegratorRep::initialize(const State& initState)

I appreciate to have your comments on it.
I thank you in advance for your time.

Best regards
Mohammad


RE: Prescribed motion

Posted: Tue Oct 27, 2009 10:42 am
by sherm
If you run exactly the same model, but without the prescribed motion constraint, does that work?

Sherm

RE: Prescribed motion

Posted: Tue Oct 27, 2009 10:51 am
by mpoursina
Dear Sherm,

It works when I comment out the prescribed motion from my code.

Mohammad

RE: Prescribed motion

Posted: Tue Oct 27, 2009 11:18 am
by sherm
Hi, Mohammad.

I'm not sure whether this is the problem, but this code is definitely wrong:
const Function::Polynomial func3(coefficients);
Constraint::PrescribedMotion(matter,&func3,...);

You have allocated the Function on the stack and then passed its address (&func3) to a parameter that expects to be able to "take ownership" of the object, meaning it must be a heap-allocated object. Try replacing those two statements with this one and let me know what happens:
Constraint::PrescribedMotion
(matter,
new Function::Polynomial(coefficients),
...);

That will create a heap-allocated object for the PrescribedMotion constraint to own.

Sherm

RE: Prescribed motion

Posted: Tue Oct 27, 2009 2:15 pm
by mpoursina
Dear Sherm,

I still have the same problem.

Regards

Mohammad

RE: Prescribed motion

Posted: Tue Oct 27, 2009 3:26 pm
by sherm
OK, let me see if I can reproduce it here.

My understanding:
This is just a 1dof pendulum on a pin joint, with a function-based prescribed motion on its one coordinate q, with the function being the polynomial q= 0*t^2 + 1*t + 0; i.e., q=t.

Is that right?

Sherm

RE: Prescribed motion

Posted: Tue Oct 27, 2009 3:48 pm
by mpoursina
Dear Sherm,

What you said is exactly what I am implementing.
The following is the related code:

#include "SimTKsimbody.h"
#include "SimTKsimbody_aux.h"
#include <cstdio>
#include <iomanip>
#include <vector>
using namespace SimTK;
using namespace std;

int main() {

// Create the system.
MultibodySystem system;
SimbodyMatterSubsystem matter(system);
GeneralForceSubsystem forces(system);
Force::UniformGravity gravity(forces, matter, Vec3(0, -9.8, 0));
Body::Rigid pendulumBody(MassProperties(1.0, Vec3(0), Inertia(1)));
pendulumBody.addDecoration(Transform(), DecorativeSphere(0.1));
MobilizedBody::Pin pendulum1(matter.Ground(), Transform(Vec3(0)),pendulumBody, Transform(Vec3(0, 1, 0)));
system.realizeTopology();
State state = system.getDefaultState();
Vector_<Real> coefficients(3);
coefficients[0] = 0.0;
coefficients[1] = 1.0;
coefficients[2] = 0.0;
Constraint::PrescribedMotion(matter,new Function::Polynomial(coefficients),MobilizedBodyIndex(1), MobilizerQIndex(0));
VerletIntegrator integ(system);
TimeStepper ts(system, integ);
ts.initialize(state);
ts.stepTo(10.0);
}

Thanks for your time.

Mohammad

RE: Prescribed motion

Posted: Tue Oct 27, 2009 4:00 pm
by sherm
Mohammad, the problem is that you didn't use a try/catch block so that you can see error messages! The code was trying to tell you what was wrong. When I ran it I got this message:
... must call realizeTopology() first

Try/catch is the only way Simbody has to report errors. You should never use it without that.
In this case you added the PrescribedMotion constraint *after* you had already realized the topology (meaning you "finalized" the system).

But before you fix the problem, surround the entire main program with try/catch and make sure you are able to see the error messages. Otherwise you'll waste a lot of time chasing minor typos and the like. This is all you have to do:

main() {
try
{
// your existing code goes here

} catch (const std::exception& e) {
std::printf("EXCEPTION THROWN: %s\n",
e.what());
exit(1);

} catch (...) {
std::printf("UNKNOWN EXCEPTION THROWN\n");
exit(1);
}

return 0; // normal return
}


Please let me know (1) whether you are able to see the error message, and (2) whether moving the realizeTopology() call fixes the problem.

Regards,
Sherm