[SOLVED] RigidMapping from Rigid3d to only some particles of FEMbody

[Simon]

RigidMapping/Connection between Rigid3d and ONLY SOME particles of “FEM body”
Hello alltogether,

I try to model a rotating structure (3D FEM “beam”). At one end its surface should be connected rigidly to a “Reference Particle” which is given a rotational boundary condition.

Using stuff from this post.
https://www.sofa-framework.org/community/forum/topic/rotating-an-object-with-a-angular-velocity-and-conserving-collisions/

I build my scene
https://github.com/swunsch/sofa-help/tree/master/05_rotating_beam/quest_1
see file: beam_rot_06_rigi_map_part_contrl_q_1_1.scn

<?xml version="1.0" ?>
<Node name="lroot" gravity="0 0 0" dt="1e-3">
    <VisualStyle displayFlags="showForceFields showBehavior" />
    <RequiredPlugin name="SofaPython" pluginName="SofaPython" />
    
    <Node name="ParticleControl" >
        <TransformEngine template="Rigid3d" name="engine" input_position="0 0 0  0 0 0 1" />
        <EulerImplicitSolver name="EulerImplicit"  rayleighStiffness="0.1" rayleighMass="0.1" />
        <CGLinearSolver name="CG Solver" iterations="25" tolerance="1e-5" threshold="1e-5"/>
        <MechanicalObject template="Rigid3d" name="myRigidParticleForRotation" position="@engine.output_position" showObject="true" showObjectScale="5"/>
        <Node name="beam">
            <EulerImplicitSolver name="cg_odesolver"   rayleighStiffness="0.5e-5" rayleighMass="0.0" />
            <CGLinearSolver name="linear solver" iterations="25" tolerance="1e-09" threshold="1e-09" />
            <MeshGmshLoader name="meshLoader" filename="beam_hexa_coarse_01_mod.msh" />
            <HexahedronSetTopologyContainer name="topo" src="@meshLoader" />
            <MechanicalObject name="dofs" src="@meshLoader" showObject="true" />
            <HexahedronSetGeometryAlgorithms template="Vec3d" name="GeomAlgo" />
            <DiagonalMass  name="computed using mass density" massDensity="7.85e-09" />
            <HexahedronFEMForceField template="Vec3d" name="FEM" method="large" poissonRatio="0.3" youngModulus="2.1e5" />
            <BoxROI name="slaveParticles" box="0 -16 -11  0 16 11" drawBoxes="true"/>
<!--             <RigidMapping input="@.." output="@." /> -->
            <RigidMapping input="@.." output="@." />
        </Node>
    </Node>
    <PythonScriptController name="angularVel" filename="angularVelBeam.py" classname="rotate"/>
</Node>

Problem 1: The Rigid mapping maps the to the whole beam, not just the particles from BoxROI but using

<RigidMapping input="@.." output="@." />

I think @. uses of the whole MechanicalObject beam so all the mechanicalDOFs of the beam. Okay.

What I need is the Mapping to only some of the DOFs of the MechanicalObject beam. I can select the particles needed with

<BoxROI name="slaveParticles" box="0 -16 -11 0 16 11" drawBoxes="true"/>

but in the
file: beam_rot_06_rigi_map_part_contrl_q_1_2.scn
I cannot map to those points/particles/DOFs from BoxROI selection

both
<RigidMapping input="@.." output="@slaveParticles" />
or
<RigidMapping input="@.." output="@slaveParticles.indices" />

are not working.

How is it possible to connect (rigid connection) only some points/particles/DOFs to the Rigid3d?

Maybe something with a SubSet… (which I haven’t really understood yet).
Thanks in advance.

[Hugo]

Hi @simon !

Looking at your scene, it seems you are more looking for boundary conditions (fixing your 3D beam) or a constraint than for a mapping.
Right?
Do you want a part of your deformable to be rigidly attached to a rigid frame?

Did you take a look at the following examples?

• examples/Components/constraint/BilateralInteractionConstraint.scn
• examples/Components/constraint/FixedConstraint.scn
• examples/Components/constraint/AttachConstraint.scn

Best

Hugo

[Simon]

oops..
and from https://www.sofa-framework.org/community/forum/topic/rigid-multibody-modeling-in-sofa/

maybe I had an overall solver missing.

https://github.com/swunsch/sofa-help/blob/master/05_rotating_beam/quest_1_1/beam_rot_12_ctrl_rigid_attach.scn

Looks much better now. Now I guess I have to adjust my rotation speed.

[Hugo]

Hey @simon

In your current scene, it is strange to have so many solvers and integration schemes. Normally you have one per object, or one for several objects. But here you have 3 solvers and 2 objects!

Moreover, I think the BilateralInteractionConstraint does correspond better to your needs. Here is an example:

<?xml version="1.0" ?>
<Node name="lroot" gravity="0 -10 0" dt="1e-5">
    <VisualStyle displayFlags="showForceFields showBehavior" />
    <RequiredPlugin name="SofaPython" pluginName="SofaPython" />

    <FreeMotionAnimationLoop />
    <GenericConstraintSolver maxIterations="1000" tolerance="0.001"/>

    <Node name="ParticleControl" >
        <TransformEngine template="Rigid3d" name="engine" input_position="0 0 0  0 0 0 1" />
        <EulerImplicitSolver name="EulerImplicit"  rayleighStiffness="0.1" rayleighMass="0.1" />
        <CGLinearSolver name="CG Solver" iterations="25" tolerance="1e-5" threshold="1e-5"/>
        <MechanicalObject template="Rigid3d" name="myRigidParticleForRotation" position="@engine.output_position" showObject="true" showObjectScale="5"/>
        <UncoupledConstraintCorrection />

        <Node name="Constraints">
            <MechanicalObject name="points" template="Vec3d" position="0 15 10 0 15 -10 0 -15 -15 0 -15 10" showObject="true" showObjectScale="5"/>
            <RigidMapping />
        </Node>
    </Node>
    <Node name="beam">
        <EulerImplicitSolver name="cg_odesolver"   rayleighStiffness="0.5e-5" rayleighMass="0.0" />
        <CGLinearSolver name="linear solver" iterations="25" tolerance="1e-09" threshold="1e-09" />
        <MeshGmshLoader name="meshLoader" filename="beam_hexa_coarse_01_mod.msh" />
        <HexahedronSetTopologyContainer name="topo" src="@meshLoader" />
        <MechanicalObject name="dofs" src="@meshLoader" showObject="true" />
        <HexahedronSetGeometryAlgorithms template="Vec3d" name="GeomAlgo" />
        <DiagonalMass  name="computed using mass density" massDensity="7.85e-09" />
        <HexahedronFEMForceField template="Vec3d" name="FEM" method="large" poissonRatio="0.3" youngModulus="2.1e5" />
        <UncoupledConstraintCorrection />

    </Node>

     <BilateralInteractionConstraint object1="@ParticleControl/Constraints/points" object2="@beam/dofs" first_point="0" second_point="2" />
     <BilateralInteractionConstraint object1="@ParticleControl/Constraints/points" object2="@beam/dofs" first_point="1" second_point="3" />
     <BilateralInteractionConstraint object1="@ParticleControl/Constraints/points" object2="@beam/dofs" first_point="2" second_point="1" />
     <BilateralInteractionConstraint object1="@ParticleControl/Constraints/points" object2="@beam/dofs" first_point="3" second_point="0" />
</Node>

[Simon]

Thanks @Hugo,

it’s going in the right direction. Sadly my earlier have not been posted (somehow the FORUM and me having some digitial issues) but this is not a problem because we have skipped alot of bad scenes from my side I think.

So I am saying it again now. Good to have you back in the forum, without you it’s not the same.

And on goes the journey.
I modified the scene:
https://github.com/swunsch/sofa-help/blob/master/05_rotating_beam/quest_1_1/beam_rot_16_rigi_bilat_face.scn

And I am seeing too much movement/flexibility in the constraint:

I want to make the constraint stiffer or somehow reduce this violation. Which parameters can be adjusted?

[Damien Marchal]

Hi Simon,

At Defrost we are using Rigidification to do so. This will associate a rigid frame to a group of deformable dofs. As this is rather complex to implment we are using some helper function to do so. The functions are implemented in STLIB.
Here: https://github.com/SofaDefrost/STLIB/blob/master/python/stlib/physics/mixedmaterial/rigidification.py

The example scene look like:

def createScene(rootNode):
        """
        """
        from stlib.scene import MainHeader
        from stlib.physics.deformable import ElasticMaterialObject
        from stlib.physics.mixedmaterial import Rigidify
        from splib.objectmodel import setData

        MainHeader(rootNode, plugins=["SofaSparseSolver"])
        rootNode.VisualStyle.displayFlags = "showBehavior"

        modelNode = rootNode.createChild("Modeling")
        elasticobject = ElasticMaterialObject(modelNode, "mesh/liver.msh", "ElasticMaterialObject")

        # Rigidification of the elasticobject for given indices with given frameOrientations.
        o = Rigidify(modelNode,
                     elasticobject,
                     name="RigidifiedStructure",
                     frames=[[0., 0., 0], [0., 0., 0]],
                     groupIndices=[[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], [48, 49, 50, 51]])

        # Activate some rendering on the rigidified object.
        setData(o.RigidParts.dofs, showObject=True, showObjectScale=1, drawMode=2)
        setData(o.RigidParts.RigidifiedParticules.dofs, showObject=True, showObjectScale=0.1,
                drawMode=1, showColor=[1., 1., 0., 1.])
        setData(o.DeformableParts.dofs, showObject=True, showObjectScale=0.1, drawMode=2)
        o.RigidParts.createObject("FixedConstraint", indices=0)

        simulationNode = rootNode.createChild("Simulation")
        simulationNode.createObject("EulerImplicitSolver")
        simulationNode.createObject("CGLinearSolver")
        simulationNode.addChild(o)
return rootNode

[Hugo]

Hey @simon

I think what you actually want to achieve is called : rigidification.
You have one (not very clean) example: examples/Components/animationloop/MechanicalMatrixMapper.pyscn

Here is another example in STLIB.

Let me know if this helps.
Best

Hugo

[Hugo]

Oups sorry I did not refresh the page Damien;)
Thanks for your answer!

[Simon]

Hi @Damien, Hi @Hugo,

thanks for your help, I will have a try at your suggestions.

[Simon]

sorry for the delay…

@damien-marchaluniv-lille1-fr

with the precompiled “DeFrost” Version of SOFA I get following erros, when running your example code. The “plugins/SoftRobots/docs/tutorials/FirstSteps” Tutorial however runs with no problem. What could be the problem?

[ERROR]   [PythonScript] ImportError: cannot import name _distributor_init
  File "stlib_example_damien.pyscn", line 8, in createScene
    from stlib.physics.mixedmaterial import Rigidify
  File "C:/Users/swtouch/SOFA/SOFA_v19.06.99_custom_Win64_v8.1/plugins/STLIB/python\stlib\physics\mixedmaterial\__init__.py", line 17, in <module>
    from rigidification import Rigidify
  File "C:/Users/swtouch/SOFA/SOFA_v19.06.99_custom_Win64_v8.1/plugins/STLIB/python\stlib\physics\mixedmaterial\rigidification.py", line 25, in <module>
    from splib.numerics import Vec3, Quat, sdiv
  File "C:/Users/swtouch/SOFA/SOFA_v19.06.99_custom_Win64_v8.1/plugins/STLIB/python\splib\numerics\__init__.py", line 38, in <module>
    import numpy
  File "C:\Users\swtouch\Anaconda2\lib\site-packages\numpy\__init__.py", line 140, in <module>
    from . import _distributor_init

[Simon]

@damien-marchaluniv-lille1-fr

could it be the / and \ in the path strings?

[Simon]

Hi @Hugo,

again I’ve Problems with the forum, I’ve written a reply but it’s not listed in this post. But I get an E-mail, that Simon has posted a new reply. And in this mail I see my post. So I try again with this reply.

I’ve adopted your suggested MechanicalMatrixMapper.pyscn example to my needs.

See:
https://github.com/swunsch/sofa-help/tree/master/05_rotating_beam/quest_1_2

But I don’t get the fixedConstraint to work now. I get the error:

[WARNING] [FixedConstraint(Fix)] Can not find the topology, wont be able to handle topological changes

[Hugo]

Hi @simon

Thank you for noticing us the issue with the forum.
When it happens (to anyone of you all), let us know as soon as possible using the contact form.

It is no error but just a warning. It should not be of a concern in your case.

I was able to run the scene, even if the result is strange:

Best

Hugo

[Simon]

Hi @Hugo,

aha ok, I understand it is not an error it’s a warning. So what I’ve build so far is the right way for my needs? And the fixedConstraint is not working properly in this case? Or have I done something wrong?

Maybe the TransFormEngine and the Constraint don’t work together well. I’ll try something in this direction, but in general I will need a way to prescribe the rotation of the Rigid3D in a time dependent manner, while its position is fixed.

I can’t see the picture, but I know the behaviour is strange.

Thanks
Simon

[Simon]

maybe I’ve messed up my TransformEngine I’ll have to check that first…

[Hugo]

Hi @simon

It works for me:

Best,

Hugo

[Simon]

Hi @Hugo,

okay cool, can I see your corresponding .pyscn file + .py(for the rotation, if you do so)?

I prepare another question because I’ve come up with something where my problem originates. maybe…

Thank you

[Simon]

I assume what I need to change in my use of TransformEngine is something from that post:

https://www.sofa-framework.org/community/forum/topic/rotation-orientation-of-a-rigid-model/

In my case Rigid at (2,0,0). I have to adopt my .py to get a local rotation, as this is what I want.

Question: The TransformEngine is always a transformation referring to the global reference frame at (0,0,0)?

This is what I assume from the behaviour of:

rigid_frame_rotation.pyscn in:

https://github.com/swunsch/sofa-help/tree/master/05_rotating_beam/quest_transform_engine

[Simon]

Hi @damien-marchaluniv-lille1-fr,

now it is working. I have forgotten to get a python27 + numpy for the DeFrost SOFA version to work.

Now I can progress.
Thanks

[Author]

Posts