Rigid multibody modeling in SOFA

[Christiane]

Hello,
I’m trying to build a rigid multibody model of a spine in SOFA. Unfortunately, it’s not going well…

I created three rigid mechanical objects (those are my vertebrae) which I try to join using JointSpringForceField. I then fix the bottom vertebra and apply a force on the top one. Problem: nothing moves. The top vertebra behaves as if it were fixed.

I’m obviously doing something wrong, but quite frankly don’t know what, even after looking at a lot of the examples. I’m pasting my code below:

<?xml version=”1.0″?>
<Node name=”root” dt=”0.1″ gravity=”0 0 -9.81″>
<VisualStyle displayFlags=”showForceFields showBehaviorModels showVisual showInteractionForceFields” />
<DefaultAnimationLoop/>

<!– Vertèbre T1 –>
<Node name=”T1″ >
<EulerSolver />
<CGLinearSolver />
<MechanicalObject name=”vertT1″ template=”Rigid” position=”-46.730 13.813 594.049 0.1264 -0.1761 0.0215 0.9760″ />

<!– Application d’une force –>
<ConstantForceField name=”forceT1″ indices=”0″ forces=”5 0 0 0 0 0″ showColor=”1 0 1 1″ arrowSizeCoef=”1.0″ />
<!LinearForceField points=”0″ forces=”5 0 0 0 0 0 1 ” />

<!– Export des positions et des forces –>
<WriteState name=”StateT1″ filename=”positionsT1.data” writeX=”1″ writeF=”1″ DOFsX=”0″ />

</Node>

<!– Vertèbre T2 –>
<Node name=”T2″ >
<EulerSolver />
<CGLinearSolver />
<MechanicalObject name=”vertT2″ template=”Rigid” position=”-54.695 10.630 578.688 -0.1270 0.1552 0.0380 0.9789″ />

<!– Export des positions et des forces –>
<WriteState name=”StateT2″ filename=”positionsT2.data” writeX=”1″ writeF=”1″ DOFsX=”0″ />

</Node>

<!– Vertèbre T3 –>
<Node name=”T3″ >
<EulerSolver />
<CGLinearSolver />
<MechanicalObject name=”vertT3″ template=”Rigid” position=”-60.588 3.019 560.031 -0.2165 0.1218 0.0763 0.9656″ />

<!– Fixation de la vertèbre –>
<!FixedConstraint template=”Rigid” name=”fixedVertebra” />

<!– Export des positions et des forces –>
<WriteState name=”StateT3″ filename=”positionsT3.data” writeX=”1″ writeF=”1″ DOFsX=”0″ />
</Node>

<!– Disques IVs –>
<JointSpringForceField template=”Rigid3d” name=”T1-T2″ object1=”@T1/vertT1″ object2=”@T2/vertT2″ spring=”BEGIN_SPRING 0 0 FREE_AXIS 1 1 1 1 1 1 KS_T 500 5000 KS_R 500 5000 REST_T -7.9650 -3.1830 -15.3610 REST_R -0.25771 0.316322 0.0132941 0.912828 END_SPRING ”/>
<JointSpringForceField template=”Rigid3d” name=”T2-T3″ object1=”@T2/vertT2″ object2=”@T3/vertT3″ spring=”BEGIN_SPRING 0 0 FREE_AXIS 1 1 1 1 1 1 KS_T 500 5000 KS_R 500 5000 REST_T 5.8930 7.6110 18.6570 REST_R -0.0820912 -0.0291694 0.0561321 0.994571 END_SPRING ”/>

</Node>

If anyone can point me in the right direction, I would appreciate it.

Regards,
Christiane

[Hugo]

Hi @chcaouette

Regarding your scene, the issue comes from the fact that each rigid position is solved independently with its own solver. Without any constraint or mapping linking the rigids, they won’t affect each other and therefore won’t move.

What you can do is to you one single solver for all rigid dof:

<?xml version="1.0"?>
<Node name="root" dt="0.01" gravity="0 0 -9.81">
    <VisualStyle displayFlags="showForceFields showBehaviorModels showVisual showInteractionForceFields" />
    <DefaultAnimationLoop/>
    <EulerImplicitSolver />
    <CGLinearSolver />

    <!-- Vertèbre T1 -->
    <Node name="T1" >
        <MechanicalObject name="vertT1" template="Rigid" position="-46.730 13.813 594.049 0.1264 -0.1761 0.0215 0.9760" />

        <!-- Application d’une force -->
        <ConstantForceField name="forceT1" indices="0" forces="5 0 0 0 0 0" showColor="1 0 1 1" arrowSizeCoef="1.0" />
        <!-- <LinearForceField points="0" forces="5 0 0 0 0 0 1 " /> -->

        <!-- Export des positions et des forces -->
        <WriteState name="StateT1" filename="positionsT1.data" writeX="1" writeF="1" DOFsX="0" />
    </Node>

    <!-- Vertèbre T2 -->
    <Node name="T2" >
        <MechanicalObject name="vertT2" template="Rigid" position="-54.695 10.630 578.688 -0.1270 0.1552 0.0380 0.9789" />

        <!-- Export des positions et des forces -->
        <WriteState name="StateT2" filename="positionsT2.data" writeX="1" writeF="1" DOFsX="0" />
    </Node>

    <!-- Vertèbre T3 -->
    <Node name="T3" >
        <MechanicalObject name="vertT3" template="Rigid" position="-60.588 3.019 560.031 -0.2165 0.1218 0.0763 0.9656" />

        <!-- Fixation de la vertèbre -->
        <FixedConstraint template="Rigid" name="fixedVertebra" />

        <!-- Export des positions et des forces -->
        <WriteState name="StateT3" filename="positionsT3.data" writeX="1" writeF="1" DOFsX="0" />
    </Node>

    <!-- Disques IVs -->
    <JointSpringForceField template="Rigid3d" name="T1-T2" object1="@T1/vertT1" object2="@T2/vertT2" spring="BEGIN_SPRING 0 0 FREE_AXIS 1 1 1 1 1 1 KS_T 500 5000 KS_R 500 5000 REST_T -7.9650 -3.1830 -15.3610 REST_R -0.25771 0.316322 0.0132941 0.912828 END_SPRING "/>
    <JointSpringForceField template="Rigid3d" name="T2-T3" object1="@T2/vertT2" object2="@T3/vertT3" spring="BEGIN_SPRING 0 0 FREE_AXIS 1 1 1 1 1 1 KS_T 500 5000 KS_R 500 5000 REST_T 5.8930 7.6110 18.6570 REST_R -0.0820912 -0.0291694 0.0561321 0.994571 END_SPRING "/>

</Node>

Moreover, there is maybe other alternative for the modeling, like the ArticulatedSystemMapping.scn available in examples/Components/mapping folder.

Cheers,

Hugo

[Christiane]

Hello Hugo,

Thanks, solving with the same solver for all rigids seems to solve the problem.

The next problem I’m having is with the JointSpringForceField. I’m able to specify the initial length using REST_T in its definition, which gets rid of any initial force in the springs; however, I’m unable to specify its initial orientation. I use REST_R to specify it, but the interface doesn’t seem to read it properly and it’s not applied. It just initializes its rotation by itself, and then when I solve I have an initial rotational load that changes my vertebrae orientations.

Any thoughts on that? Is it a bug, or do I have a typo someplace?

Regards,
Christiane

[Hugo]

Hi @chcaouette

I never used the JointSpringFF. What I can tell is that if you want to define rest_position and rotation of your DOFs, the best is to specify them in the rest_position=”x y z qw qx qy qz”. Moreover, in the code it is noted that the REST_R corresponds to the offset orientation of the spring.

Does this fit what you want to do?

Hugo

[Author]

Posts