Understanding Linear Velocities in a Rotational System Simulation Using Simbody

[Amit Misra]

First of all, sorry to post this question during the holiday period, but I'm hoping someone might be able to help.

I'm working on a simulation using Simbody where I have a chassis connected to the ground via a free joint, with two motors and wheels attached. The motors apply torque to the wheels, causing them to rotate. I expected that, since there are no external forces and the system should conserve linear momentum, the chassis would not exhibit any net linear motion - only rotational motion due to the internal torques.

However, when I run the simulation, I'm observing non-zero linear velocities for the chassis. Here is a snippet of the output data I'm getting:

Time Angular Velocities (deg/sec) Linear Velocities (m/sec)
0 0 0 0 ~[0,0,0]
0.1 3.38813e-21 0.108405 -0.308656 ~[-0.00477689,5.21783e-23,1.05879e-22]
0.2 1.89735e-19 0.216809 -0.617312 ~[-0.00955379,-5.25913e-21,-4.23516e-21]
0.3 1.0842e-18 0.325214 -0.925967 ~[-0.0143307,1.51674e-20,-1.86347e-20]


I understand that for any point on a rotating body, the linear velocity is given by v = omega x r, where omega is the angular velocity and r is the position vector from the axis of rotation. This means that even if the center of mass isn't translating, points on the chassis should have linear velocities due to rotation.

My questions are:

Should the chassis exhibit linear velocity in this scenario? Is the observed linear velocity a result of rotational motion (v = omega x r), or is it indicating an unintended translation of the chassis's center of mass?

How can I determine whether the linear velocities are due to rotation or translation? What steps can I take to verify and resolve this issue in my simulation?

Are there any best practices in Simbody to prevent unintended linear motion when simulating rotational systems? For example, should I constrain the chassis's translation or adjust integrator settings to improve numerical stability?

Any insights or suggestions would be greatly appreciated!

Additional Information:

The simulation does not include external forces like gravity or friction.
The chassis is connected to the ground using a MobilizedBody::Free joint.
The motors are applying constant torque to the wheels via a custom force implementation.

Code: Select all

#include <Simbody.h>

using namespace SimTK;


class ShowEnergy : public DecorationGenerator {
public:
    explicit ShowEnergy(const MultibodySystem& mbs,
                        const MobilizedBody::Free& shaft)
        : m_mbs(mbs), m_shaft(shaft) {}
    void generateDecorations(const State&                state,
                             Array_<DecorativeGeometry>& geometry) override;
private:
    const MultibodySystem& m_mbs;
    const MobilizedBody::Free& m_shaft;
};

// Custom motor force class
class MyMotorForce : public Force::Custom::Implementation {
private:
    MobilizedBody::Pin m_wheel;
    Real m_forceMagnitude;

public:
    MyMotorForce(const MobilizedBody::Pin& wheel, Real forceMagnitude)
        : m_wheel(wheel), m_forceMagnitude(forceMagnitude) {
    }

    void calcForce(const State& state,
                   Vector_<SpatialVec>& bodyForces,
                   Vector_<Vec3>& particleForces,
                   Vector& mobilityForces) const override {
        m_wheel.applyOneMobilityForce(state, 0, m_forceMagnitude,
                                      mobilityForces);
    }


    Real calcPotentialEnergy(const State&) const override { return 0; }
};

int main() {
    try {
        // Create the system and add gravity
        MultibodySystem system;
        SimbodyMatterSubsystem matter(system);
        GeneralForceSubsystem forces(system);
        // Force::Gravity gravity(forces, matter, Vec3(0, 0, 0));

        // matter.setShowDefaultGeometry(false);

        // Define the bodies
        Body::Rigid chassisBody(MassProperties(1.0, Vec3(0), UnitInertia(1)));
        Body::Rigid motorBody(MassProperties(0.1, Vec3(0), UnitInertia(0.01)));
        Body::Rigid wheelBody(MassProperties(1, Vec3(0), UnitInertia(0.1)));
        Body::Rigid motorBody2(MassProperties(0.1, Vec3(0), UnitInertia(0.01)));
        Body::Rigid wheelBody2(MassProperties(1, Vec3(0), UnitInertia(0.1)));


        // Decorate bodies
        const Vec3 halfLengths(.02, .04, .3);
        chassisBody.addDecoration(Transform(),
                                  DecorativeBrick(Vec3(1, 1, 1)).setColor(Red));

        Rotation YtoZ(Pi / 2, XAxis);
        motorBody.addDecoration(YtoZ,
                                DecorativeCylinder(0.01, 0.01).setColor(Blue));
        // Rotation YtoX(Pi, ZAxis);

        wheelBody.addDecoration(YtoZ,
                                DecorativeCylinder(1, 0.05).setColor(Green));

        motorBody2.addDecoration(Transform(),
                                DecorativeCylinder(0.01, 0.01).setColor(Orange));

        Rotation YtoZ2(Pi / 2, XAxis);
        wheelBody2.addDecoration(Transform(),
                                DecorativeCylinder(1, 0.05).setColor(Orange));

        MobilizedBody::Free chassis(matter.Ground(), Transform(Vec3(0)),
                                    chassisBody, Transform(Vec3(0)));
        MobilizedBody::Weld motor(chassis, Transform(Vec3(0, 0, 1)), motorBody,
                                  Transform(Vec3(0)));

        MobilizedBody::Pin wheel(motor, Transform(Vec3(0, 0, 1)), wheelBody,
                                 Transform(Vec3(0)));

        MobilizedBody::Weld motor2(chassis, Transform(Vec3(0, 1, 0)), motorBody2,
                                  Transform(Vec3(0)));

        MobilizedBody::Pin wheel2(motor2, Transform(Vec3(0, 1, 0)), wheelBody2,
                                 Transform(Vec3(0)));
        // wheelBody.addDecoration(Transform(), DecorativeCylinder(0.5, 0.1).setColor(Blue));

        Visualizer viz(system);
        viz.setDesiredFrameRate(60);
        viz.setShowFrameRate(true);
        viz.setBackgroundType(Visualizer::BackgroundType::SolidColor);
        viz.addDecorationGenerator(new ShowEnergy(system, chassis));
        system.addEventReporter(new Visualizer::Reporter(viz, 1.0 / 60.0));

        // Add the motor force
        Force::Custom(forces, new MyMotorForce(wheel, 0.174533));

        // Initialize the system and state
        State state = system.realizeTopology();
        // State state = system.getDefaultState();

        // Integrate the system over time
        RungeKuttaMersonIntegrator integrator(system);
        integrator.setAccuracy(1.0e-6);
        integrator.setMaximumStepSize(1e-3);
        integrator.initialize(state);

        TimeStepper ts(system, integrator);
        ts.initialize(state);
        // ts.stepTo(100.0);

        while (integrator.getTime() < 2000.0) {
            integrator.stepTo(integrator.getTime() + 0.1);
            const State& state = integrator.getState();
            viz.report(integrator.getState());
            // TODO: Save or analyze state
        }
    } catch (const std::exception& e) {
        std::cout << "EXCEPTION: " << e.what() << std::endl;
        return 1;
    }
    return 0;


}


void ShowEnergy::generateDecorations(const State&                state,
                                     Array_<DecorativeGeometry>& geometry)
{
    m_mbs.realize(state, Stage::Dynamics);
    const Real E=m_mbs.calcEnergy(state);
    DecorativeText energy;
    energy.setTransform(Vec3(-.2,0,.5))
            .setText("Energy: " + String(E, "%.6f"))
            .setScale(.09)
            .setColor(Black);
    geometry.push_back(energy);

    // Shaft S angular velocity w_GS comes back expressed in G; we want to see it in S.
    const Rotation& R_GS = m_shaft.getBodyRotation(state);
    const auto liner_vel = m_shaft.getBodyVelocity(state);
    const auto& R_SG = ~R_GS; // Transpose is inverse.
    const auto w_GS_G = 180/Pi * m_shaft.getBodyAngularVelocity(state);
    const auto w_GS_S = R_SG * w_GS_G;
    const auto liner_vel_body = R_SG * liner_vel;
    std::cout << state.getTime() << '\t' << w_GS_S[0] << ' ' << w_GS_S[1] << ' ' << w_GS_S[2] << ' '<< liner_vel_body[1]<<'\n';
}

[Michael Sherman]

Sorry for the delay in responding to this.

I ran your program to see what's happening. My guess is you didn't realize that when you apply a mobility force it is applied equal and opposite to the two connected bodies. So a torque is being applied both to the actuated wheel and the free-floating chassis. You can see this clearly if you reduce the applied torque by a factor of 1000:

Code: Select all

// Add the motor force
Force::Custom(forces, new MyMotorForce(wheel, .001 * 0.174533));

In that case you can see the green wheel accelerate in one direction and the chassis react in the opposite direction.

If you want to apply force only to the wheel you can do that with a body force, but a mobility force always acts on both bodies.

Regards,
Sherm

[Amit Misra]

Thank you for your detailed explanation, Sherm. I appreciate the insights about how mobility forces behave and the equal and opposite torques applied to connected bodies. That cleared up a lot of confusion.
Based on your explanation, I have updated the code and applied the torque directly to the wheel body. Here is the added code snippet for reference

Code: Select all

Vec3 torqueInG(0, 0, m_forceMagnitude);

// A SpatialVec is {angular, linear}. To apply a pure torque:
SpatialVec torqueOnly(torqueInG, Vec3(0));

// Apply this torque directly to the wheel body.
m_wheel.applyBodyForce(state, torqueOnly, bodyForces);

Could you kindly review this approach and let me know if this correctly isolates the force application to the wheel body? Your feedback would be invaluable.


Complete Code for Reference

Code: Select all

#include <Simbody.h>

using namespace SimTK;


class ShowEnergy : public DecorationGenerator {
public:
    explicit ShowEnergy(const MultibodySystem& mbs,
                        const MobilizedBody::Free& shaft)
        : m_mbs(mbs), m_shaft(shaft) {}
    void generateDecorations(const State&                state,
                             Array_<DecorativeGeometry>& geometry) override;
private:
    const MultibodySystem& m_mbs;
    const MobilizedBody::Free& m_shaft;
};

// Custom motor force class
class MyMotorForce : public Force::Custom::Implementation {
private:
    MobilizedBody::Pin m_wheel;
    Real m_forceMagnitude;

public:
    MyMotorForce(const MobilizedBody::Pin& wheel, Real forceMagnitude)
        : m_wheel(wheel), m_forceMagnitude(forceMagnitude) {
    }

    void calcForce(const State& state,
                   Vector_<SpatialVec>& bodyForces,
                   Vector_<Vec3>& particleForces,
                   Vector& mobilityForces) const override {
        // m_wheel.applyOneMobilityForce(state, 0, m_forceMagnitude,
        //                               mobilityForces);

        Vec3 torqueInG(0, 0, m_forceMagnitude);

        // A SpatialVec is {angular, linear}. To apply a pure torque:
        SpatialVec torqueOnly(torqueInG, Vec3(0));

        // Apply this torque directly to the wheel body.
        m_wheel.applyBodyForce(state, torqueOnly, bodyForces);

        // No mobilityForces are applied here, so no equal and opposite reaction
        // is generated via a joint coordinate. Just the wheel gets the torque.
    }


    Real calcPotentialEnergy(const State&) const override { return 0; }
};

int main() {
    try {
        // Create the system and add gravity
        MultibodySystem system;
        SimbodyMatterSubsystem matter(system);
        GeneralForceSubsystem forces(system);
        // Force::Gravity gravity(forces, matter, Vec3(0, 0, 0));

        // matter.setShowDefaultGeometry(false);

        // Define the bodies
        Body::Rigid chassisBody(MassProperties(1.0, Vec3(0), UnitInertia(1)));
        Body::Rigid motorBody(MassProperties(0.1, Vec3(0), UnitInertia(0.01)));
        Body::Rigid wheelBody(MassProperties(1, Vec3(0), UnitInertia(0.1)));
        Body::Rigid motorBody2(MassProperties(0.1, Vec3(0), UnitInertia(0.01)));
        Body::Rigid wheelBody2(MassProperties(1, Vec3(0), UnitInertia(0.1)));


        // Decorate bodies
        const Vec3 halfLengths(.02, .04, .3);
        chassisBody.addDecoration(Transform(),
                                  DecorativeBrick(Vec3(1, 1, 1)).setColor(Red));

        Rotation YtoZ(Pi / 2, XAxis);
        motorBody.addDecoration(YtoZ,
                                DecorativeCylinder(0.01, 0.01).setColor(Blue));
        // Rotation YtoX(Pi, ZAxis);

        wheelBody.addDecoration(YtoZ,
                                DecorativeCylinder(1, 0.05).setColor(Green));

        motorBody2.addDecoration(Transform(),
                                DecorativeCylinder(0.01, 0.01).setColor(Orange));

        Rotation YtoZ2(Pi / 2, XAxis);
        wheelBody2.addDecoration(Transform(),
                                DecorativeCylinder(1, 0.05).setColor(Orange));

        MobilizedBody::Free chassis(matter.Ground(), Transform(Vec3(0)),
                                    chassisBody, Transform(Vec3(0)));
        MobilizedBody::Weld motor(chassis, Transform(Vec3(0, 0, 1)), motorBody,
                                  Transform(Vec3(0)));

        MobilizedBody::Pin wheel(motor, Transform(Vec3(0, 0, 1)), wheelBody,
                                 Transform(Vec3(0)));

        MobilizedBody::Weld motor2(chassis, Transform(Vec3(0, 1, 0)), motorBody2,
                                  Transform(Vec3(0)));

        MobilizedBody::Pin wheel2(motor2, Transform(Vec3(0, 1, 0)), wheelBody2,
                                 Transform(Vec3(0)));
        // wheelBody.addDecoration(Transform(), DecorativeCylinder(0.5, 0.1).setColor(Blue));

        Visualizer viz(system);
        viz.setDesiredFrameRate(60);
        viz.setShowFrameRate(true);
        viz.setBackgroundType(Visualizer::BackgroundType::SolidColor);
        viz.addDecorationGenerator(new ShowEnergy(system, chassis));
        system.addEventReporter(new Visualizer::Reporter(viz, 1.0 / 60.0));

        // Add the motor force
        Force::Custom(forces, new MyMotorForce(wheel, 1.174533));

        // Initialize the system and state
        State state = system.realizeTopology();
        // State state = system.getDefaultState();

        // Integrate the system over time
        RungeKuttaMersonIntegrator integrator(system);
        integrator.setAccuracy(1.0e-6);
        integrator.setMaximumStepSize(1e-3);
        integrator.initialize(state);

        TimeStepper ts(system, integrator);
        ts.initialize(state);
        // ts.stepTo(100.0);

        while (integrator.getTime() < 2000.0) {
            integrator.stepTo(integrator.getTime() + 0.1);
            const State& state = integrator.getState();
            viz.report(integrator.getState());
            // TODO: Save or analyze state
        }
    } catch (const std::exception& e) {
        std::cout << "EXCEPTION: " << e.what() << std::endl;
        return 1;
    }
    return 0;


}


void ShowEnergy::generateDecorations(const State&                state,
                                     Array_<DecorativeGeometry>& geometry)
{
    m_mbs.realize(state, Stage::Dynamics);
    const Real E=m_mbs.calcEnergy(state);
    DecorativeText energy;
    energy.setTransform(Vec3(-.2,0,.5))
            .setText("Energy: " + String(E, "%.6f"))
            .setScale(.09)
            .setColor(Black);
    geometry.push_back(energy);

    // Shaft S angular velocity w_GS comes back expressed in G; we want to see it in S.
    const Rotation& R_GS = m_shaft.getBodyRotation(state);
    const auto liner_vel = m_shaft.getBodyVelocity(state);
    const auto& R_SG = ~R_GS; // Transpose is inverse.
    const auto w_GS_G = 180/Pi * m_shaft.getBodyAngularVelocity(state);
    const auto w_GS_S = R_SG * w_GS_G;
    const auto liner_vel_body = R_SG * liner_vel;
    std::cout << state.getTime() << '\t' << w_GS_S[0] << ' ' << w_GS_S[1] << ' ' << w_GS_S[2] << ' '<< liner_vel_body[1]<<'\n';
}

[Amit Misra]

Dear Sherm,

Sorry to bother you again, but I am still struggling to fully understand and visualize the issue you explained. I deeply appreciate your earlier response, but I think I may still be misunderstanding some key points.

My main goal is to verify the conservation of momentum in the system. If I rotate the wheel in one direction, the chassis should rotate in the opposite direction to conserve angular momentum. This works perfectly when the chassis and motor rotation are symmetrical.

However, when I introduce asymmetry into the system by adding another wheel or offsetting components to make the chassis unsymmetrical, I expected there to be some wobbling or vibrations due to the imbalance but not linear motion of the center of mass. Instead, I observed a linear velocity in the chassis, and I am confused about where my reasoning might be flawed.

You mentioned earlier that "a mobility force is applied equal and opposite to the two connected bodies," and I thought this might mean I was inadvertently introducing a linear force along with the torque, causing the linear motion. To address this, I modified the force calculation to directly apply the torque to the wheel body instead of using mobility forces. However, while this resolved the linear motion, I ended up with no motion at all, which is also incorrect.

Interestingly, when I replace the free joint with a ball joint, the simulation produces the expected behavior. I understand that a ball joint only allows rotational motion, which may explain why linear velocities are not observed. Here is the output from the simulation using a ball joint:

0 0 0 0 ~[0,0,0]
0.1 0 0 -2.24601 ~[0,0,0]
0.2 0 0 -4.49203 ~[0,0,0]
0.3 0 0 -6.73804 ~[0,0,0]
0.4 0 0 -8.98405 ~[0,0,0]
0.5 0 0 -11.2301 ~[0,0,0]
0.6 0 0 -13.4761 ~[0,0,0]
0.7 0 0 -15.7221 ~[0,0,0]


This behavior aligns with my expectations, but it has left me more confused about why the free joint case produces linear motion. Could you help clarify why linear motion arises in the asymmetrical case and whether my approach to applying the torque directly to the wheel is correct? Perhaps I am overlooking something fundamental about how forces and torques interact in Simbody or how they propagate through the system.

Thank you again for your time and patience.
I truly appreciate your help.

[Michael Sherman]

Hi, Amit.

I'm not sure I understand what you're trying to do, physically. The mobility force you had before very accurately models a motor at the joint -- motors apply equal and opposite torques to the armature and case. The current set up with a body force is a kind of external "hand of god" where a torque is applied to the wheel with no obvious physical mechanism. I suppose you could achieve this physically by attaching a couple of rockets to the wheel (like a pinwheel firework).

I applied a torque also to "wheel2" and noticed that you have the joint axis pointed incorrectly -- a Pin joint always rotates around the common z axis of the two joint frames (F and M in the documentation). You'll need to rotate those frames to align z with the wheel axis.

Regarding conservation of angular momentum, that applies to isolated systems but you are applying external forces so I'm not sure whether your expectations makes sense. My suggestion would be to start by simulating some simple textbook systems to gain confidence that you are seeing correct physics (I can guarantee that Simbody is doing that correctly for a simple system like this one). Also, it's important to verify that the model you built is the one you wanted -- the odd axis direction for the orange wheel suggests there might still be a modeling bug.

BTW it's hard to get a sense of what's happening visually with the huge torque that you're using. I had to reduce it by a factor of 1000 so the wheel wouldn't spin so fast that I'd see a strobe effect with the 60fps visualization. I'd suggest using small forces until you are confident everything is set up correctly.

Regards,
Sherm

[Amit Misra]

Hi Sherm,

Thank you for your detailed feedback and observations. I truly appreciate you taking the time to review my setup.

What I am trying to model is a simple reaction wheel system. In this setup, one body (the base) is free-floating in space, and the other (the wheel) is controlled by a motor. The motor applies torque to the wheel, which in turn causes the base to rotate in the opposite direction, enabling orientation control of the free-floating body.

The "hand of god" torque you mentioned is intended to represent the torque applied by the motor to the wheel in my simulation.

Regarding the higher torque values, I used them to exaggerate the effect for faster visualization, but I understand how this might complicate interpreting the results. I will adjust to more realistic values for improved accuracy.

While testing this setup, I noticed some linear motion in the system, which seems unusual to me. I am not sure why it is occurring and would appreciate your thoughts on what might be causing it.

I hope this explanation clarifies my setup and observations. Thank you again for your guidance!

[Michael Sherman]

Hi, Amit. If a reaction wheel is activated by a motor, you should switch back to using a mobility force on a Pin joint to model that, so that you correctly get torques on both bodies. As I mentioned, the current implementation is similar to a fireworks pinwheel where you have a pair of rockets attached to the wheel. I don't think that's a good model for a reaction wheel. Also, those are external forces but I think you want only internal ones.

I didn't notice the linear motion in the visualization, maybe I just needed some object fixed to Ground to see the chassis motion. Numerical integration does introduce errors depending on the accuracy setting -- you could check that by running at a very tight tolerance like 1e-8 to see whether the motion is a numerical artifact or is part of the physics as modeled.

Regards, Sherm

[Amit Misra]

Hi Sherm,

Thank you for your feedback and for pointing out the distinction between external and internal forces. I am a bit confused about the comment: "Also, those are external forces but I think you want only internal ones."

Currently, I apply the torque using the following code:

Code: Select all

// Custom Motor Force Implementation
class MotorForce : public Force::Custom::Implementation {
public:
    MotorForce(SimbodyMatterSubsystem& matter, const MobilizedBody::Pin& wheel,
               Real maxTorque, Real frequency)
        : matter(matter), wheel(wheel), maxTorque(maxTorque), frequency(frequency) {}

    // Required virtual methods implementation
    void calcForce(const State& state, Vector_<SpatialVec>& bodyForces,
                  Vector_<Vec3>& particleForces, Vector& mobilityForces) const override {
        Real t = state.getTime();
        Real torque = maxTorque * std::sin(2 * Pi * frequency * t);


         wheel.applyOneMobilityForce(state,0,torque,mobilityForces);
    }

    Real calcPotentialEnergy(const State& state) const override {
        return 0;
    }

    // Additional methods for parameter adjustment
    void setMaxTorque(Real torque) {
        maxTorque = torque;
    }

    void setFrequency(Real freq) {
        frequency = freq;
    }

private:
    SimbodyMatterSubsystem& matter;
    MobilizedBody::Pin wheel;
    Real maxTorque;
    Real frequency;
};

Code: Select all

Force::Custom motorForce(forces, new MotorForce(matter, wheel, maxTorque, frequency));

Would this be considered an external force? If so, could you please provide guidance on how to apply this as an internal force? For example, would using a mobility force through the pin joint, like pinJoint.applyMobilityForce, correctly represent an internal force?

I will also test with a tighter tolerance (1e-8) as suggested to check if the linear motion is a numerical artifact or part of the physics as modeled.

Thank you again for your insights. I greatly appreciate your help!

complete code

Code: Select all

#include <Simbody.h>
#include <iostream>

using namespace SimTK;

// Custom Motor Force Implementation
class MotorForce : public Force::Custom::Implementation {
public:
    MotorForce(SimbodyMatterSubsystem& matter, const MobilizedBody::Pin& wheel,
               Real maxTorque, Real frequency)
        : matter(matter), wheel(wheel), maxTorque(maxTorque), frequency(frequency) {}

    void calcForce(const State& state, Vector_<SpatialVec>& bodyForces,
                  Vector_<Vec3>& particleForces, Vector& mobilityForces) const override {
        Real t = state.getTime();
        Real torque = maxTorque * std::sin(2 * Pi * frequency * t);
        wheel.applyOneMobilityForce(state, 0, torque, mobilityForces);
    }

    Real calcPotentialEnergy(const State& state) const override {
        return 0;
    }

    void setMaxTorque(Real torque) {
        maxTorque = torque;
    }

    void setFrequency(Real freq) {
        frequency = freq;
    }

private:
    SimbodyMatterSubsystem& matter;
    MobilizedBody::Pin wheel;
    Real maxTorque;
    Real frequency;
};

int main() {
    MultibodySystem system;
    SimbodyMatterSubsystem matter(system);
    GeneralForceSubsystem forces(system);

    Body::Rigid groundBody(MassProperties(1, Vec3(0), Inertia(1)));
    MobilizedBody::Ground ground(matter.Ground());

    const Real chassisMass = 100.0;
    Body::Rigid chassisBody(MassProperties(chassisMass, Vec3(0),
        UnitInertia(1.0, 1.0, 1.0)));

    chassisBody.addDecoration(Transform(),
       DecorativeBrick(Vec3(0.5, 0.2, 0.1)).setColor(Blue));

    MobilizedBody::Free chassis(ground, Transform(), chassisBody, Transform());

    const Real wheelMass = 10.0;
    Body::Rigid wheelBody(MassProperties(wheelMass, Vec3(0), UnitInertia(1.0, 1.0, 1.0)));

    const Real motorBlockMass = 2.0;
    Body::Rigid motorBlockBody(MassProperties(motorBlockMass, Vec3(0), UnitInertia(0.1, 0.1, 0.1)));

    Rotation wheelRotation(Pi / 2, YAxis);

    motorBlockBody.addDecoration(wheelRotation, DecorativeBrick(Vec3(0.1, 0.1, 0.1)).setColor(Green));

    MobilizedBody::Weld motorBlock(chassis, Transform(wheelRotation, Vec3(0, 0.2, 0)),
        motorBlockBody, Transform());

    wheelBody.addDecoration(wheelRotation,
        DecorativeCylinder(0.2, 0.2).setColor(Red));

    MobilizedBody::Pin wheel(chassis, Transform(wheelRotation, Vec3(1, 0, 0)),
                           wheelBody, Transform(wheelRotation));

    Real maxTorque = 1.0;
    Real frequency = 1.0;
    Force::Custom motorForce(forces, new MotorForce(matter, wheel, maxTorque, frequency));

    Visualizer viz(system);
    viz.setShowFrameRate(true);
    viz.setBackgroundType(Visualizer::SolidColor);
    system.addEventReporter(new Visualizer::Reporter(viz, 1.0 / 30));

    system.realizeTopology();
    State state = system.getDefaultState();

    RungeKuttaMersonIntegrator integ(system);
    integ.setAccuracy(1e-12);
    TimeStepper ts(system, integ);
    ts.initialize(state);

    for (int i = 0; i < 100; ++i) {
        ts.stepTo(i);
        auto newstate = ts.getState();
        system.realize(newstate, Stage::Dynamics);
        std::cout << chassis.getBodyVelocity(newstate) << std::endl;
        std::cout << wheel.getBodyVelocity(newstate) << std::endl;
    }

    return 0;
}

[Michael Sherman]

Hi, Amit. Yes, applying a mobility forces is what I meant by "internal". It models a motor on the wheel axle. It was the earlier body force that was external (the rocket-powered pinwheel).

Regards,
Sherm