Simulating deformation of a simple set of beams via BeamFEMForceField

[Andrija]

Hi everyone,

I am trying to simulate the deformation of simple set of beams, where beams are oriented differently in space, by using BeamFEMForceField. The end nodes of the beams are modeled as MechanicalObject with template Rigid3d. I am using "MeshTopology in form of lines. Although it seems like the simulation is working as it should (deformation of the beams looks fine), I am having problem to orient the lines to follow the direction how my nodes (the Rigid3d degrees of freedom are connected) are connected. I pasted the code below. How it is now the MeshToplogy lines are orineted in x direction, but they should be oriented so to follow the direction how my Rigid3d degrees of freedom are connected. Do you know how I can change the orientation of lines or some other suggestions how to do it? In all of the documentation examples the BeamFEMForceField is demonstrated on a simple horizontal beam, but how to simulate a beam that is under some angle in space (beam node coordinates are so that the beam is under some angle in space)?

This is how it looks now:

This is how it should look (connected with light color):

Here is the code (disregard the Visual model for now):

# Importing modules
import Sofa
import os

def createScene(rootNode):

    rootNode.createObject('RequiredPlugin', name='OpenVis', pluginName='SofaOpenglVisual')

    # Solvers
    rootNode.createObject('EulerImplicitSolver', name="cg odesolver", printLog="false", rayleighStiffness="0.1", rayleighMass="0.1")
    rootNode.createObject('CGLinearSolver', iterations="100", name="linear solver", threshold="1e-20", tolerance="1e-20")

    mechanics = rootNode.createChild("Mechanics")
    mechanics.createObject("MechanicalObject", template="Rigid3d", name="DOFs", position=[                                                                                      
        [2, 2, 22, 0, 0, 0, 1],
        [12.5, 2, 10, 0, 0, 0, 1],                                                                                       
        [23, 2, -2, 0, 0, 0, 1],
        [35, 9.5, 8, 0, 0, 0, 1],                                                                                        
        [48, 17, 18, 0, 0, 0, 1]])                                                                                      
    mechanics.createObject("BeamFEMForceField", name="FEM", radius=0.1, youngModulus=80000000, poissonRatio=0.4)
    mechanics.createObject("MeshTopology", name="lines", lines="0 1 1 2 2 3 3 4")
    mechanics.createObject("UniformMass", template="Rigid3d", name="mass", vertexMass="0.1 0.1 [1 0 0,0 1 0,0 0 1]")
    mechanics.createObject("FixedConstraint", template="Rigid3d", name="fixOrigin", indices="0")

    # Defining viusal model
    Visu = mechanics.createChild('Visu')
    Visu.createObject('OglModel', name="Visual", position="-1 -0.5 -0.5 -1 0.5 -0.5 -1 0.5 0.5 -1 -0.5 0.5 1 -0.5 -0.5 1 0.5 -0.5 1 0.5 0.5 1 -0.5 0.5 -1 -0.5 -0.5 -1 0.5 -0.5 -1 0.5 0.5 -1 -0.5 0.5 1 -0.5 -0.5 1 0.5 -0.5 1 0.5 0.5 1 -0.5 0.5 -1 -0.5 -0.5 -1 0.5 -0.5 -1 0.5 0.5 -1 -0.5 0.5 1 -0.5 -0.5 1 0.5 -0.5 1 0.5 0.5 1 -0.5 0.5 -1 -0.5 -0.5 -1 0.5 -0.5 -1 0.5 0.5 -1 -0.5 0.5 1 -0.5 -0.5 1 0.5 -0.5 1 0.5 0.5 1 -0.5 0.5", triangles="3 1 0 3 2 1 3 6 2 3 7 6 7 5 6 7 4 5 4 1 5 4 0 1 5 1 2 2 6 5 4 7 3 4 3 0 11 9 8 11 10 9 11 14 10 11 15 14 15 13 14 15 12 13 12 9 13 12 8 9 13 9 10 10 14 13 12 15 11 12 11 8 19 17 16 19 18 17 19 22 18 19 23 22 23 21 22 23 20 21 20 17 21 20 16 17 21 17 18 18 22 21 20 23 19 20 19 16 27 25 24 27 26 25 27 30 26 27 31 30 31 29 30 31 28 29 28 25 29 28 24 25 29 25 26 26 30 29 28 31 27 28 27 24")
    Visu.createObject('RigidMapping', template="Rigid3d,ExtVec3f", rigidIndexPerPoint="1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4", input="@..", output="@Visual")

Many thanks!

Best regards,
Andrija

[Hugo]

Hi @andrija987

Could you tell me what you mean by “problem to orient the lines to follow the direction how my nodes”?
I am not sure to understand the problem.

To visualize the lines I just modified the MeshTopology line as follows:

mechanics.createObject("MeshTopology", name="lines", lines="0 1 1 2 2 3 3 4", position="@DOFs.position", drawEdges="1")

Note that the orientation of each frame connecting the beam is defined in the quaternion part in the degrees of freedom : x y z + quaternion (=Rigid3d)

Best

Hugo

[Andrija]

Hi @hugo,

Many thanks for the reply and great help!!!

Actually that did the trick and solved what I was searching for. This is what I asked. Sorry for some reason the links to figures I provided didn’t work out, and I provided the figures to better illustrate what I was asking. But the line of code you wrote worked quite nicely for me.

Sorry I have another question in relation to this. How it is possible to fix all degrees of freedom in only certain desired indices? I mean how to constrain motion, both translation and rotation degrees of freedom? In the example and code above, how I can fix both translation and rotation of indices 0?

Using the FixedConstraint, seems to constraint only translation degrees of freed but not rotation. Then when I try to use FixedRotationConstraint, there is no option to constrain the rotation of only indices I want (e.g. indices 0), but rather when I use it it constraints rotation of all the indices. Or maybe I don’t know how to apply this constrain in the right way. Can you maybe help you with some suggestions what I can use?

If you are free/interested I can share more about my project in SOFA and what we are trying to realize, maybe we can promote SOFA through our project, and vice verse. Feel free to reach me out for possible organizing meeting in the future.

Thanks!

Best,
Andrija

[Hugo]

Hi @andrija987

I am not myself working extensively with rigids, and I never noticed that the FixedRotationConstraint was applied on all points. From what I understand, you would like to couple both the FixedConstraint and FixedRotationConstraint but only on some indices. Is this correct?
If so, you (if you feel to do it, it’s the best way to learn) or I (if I have some time) could implement this.

Promotion of SOFA, citations are always most welcome within the open-source community. The more we are, the better the project will become!

We could chat with the next weeks, if you’d like. I would be interested to learn more about your project.
Best wishes,

Hugo

[Andrija]

Hi @hugo,

Many thanks again for the reply.

Yes that is what I would like. To couple or use bothFixedConstraint and FixedRotationConstraint but apply them only to desired indices. In my case I would like to fix both translation and rotation of the indices 0. Is seems like there is no option that I use FixedRotationConstraint but apply only to the desired indices. When I write this in my code/scene:
mechanics.createObject("FixedRotationConstraint", template="Rigid3d", name="fixRotation", indices="0", FixedXRotation="1", FixedYRotation="1", FixedZRotation="1")
I get the warning message meaning that this attribute for selecting incidences was not used.
Unused Attribute: "indices" with value: "0" (/Mechanics/fixRotation)
I really don’t know how to apply this FixedRotationConstraint to only for example indices 0? Or if there is some other ways to constrain all degrees of freedom of indices 0 (desired indices in general)?

It would be really great to have a chat with you then I can share in more depth what we are trying to do at the end, and what is our project about. Think it would be interested for you and SOFA. Of course we are more then happy to promote SOFA with citations in our future papers, in general as a really great framework to use.

I am free in the next weeks, feel free to write me or reach me out at: Andrija.Milojevic@lut.fi, or use my email address I use and registered here at SOFA forum. Then we can correspond about the meeting. It would be great to share what we are doing in SOFA, or trying to do 🙂 .

I am looking forward hearing from you!

Best regards,
Andrija

[Author]

Posts