Simulating deformation of a simple set of beams via BeamFEMForceField

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