Attaching 2 rods to a platform with different objects

[Benjamin]

Hello everyone,

I am working with SOFA 19.06. I would like to create 2 different beams using BeamInterpolation and adaptiveBeamForceFieldAndMass connect to a platform. I already know how to do this with the dof0TransformNode0 and dof1TransformNode1 parameters but with this, I cannot get the force and moment at the end of the clamping points. That’s why, I would like to connect differently the rod to the platform. I tried to use AttachConstraint but I did not suceed to create it. I think the problem come from the different rotations and positions of the frames between the platform and the ends of the rods.

I have several questions relatively to this.
1- Are solvers necessary for each object or only for the first/main child ?
2- Is AttachConstraint is the best way to somehow rewrite the dof0/1TransformNode0/1 effect/parameters, or do you think RigidMapping or another function could be more suitable ?
3- If I create two other points on the platform object representating the clamping points. How can link them together to respect rigid conditions even if those points do not have the same orientation ?

Here is the code I am working on:

import Sofa, numpy
from splib.constants import Key
from math import sin, cos, sqrt
from splib.animation import AnimationManager
stepTime=1.0e-3
gravity=0.0#-9.81e3
NumberElement=10

edgesList=['0 1','1 2','2 3','3 4','4 5','5 6','6 7','7 8','8 9','9 10']
positionB1=['110.0 0.0 0.0    0.0 0.0 1.0 0.0', '100.0 0.0 0.0    0.0 0.0 1.0 0.0', '90.0 0.0 0.0    0.0 0.0 1.0 0.0','80.0 0.0 0.0      0.0 0.0 1.0 0.0', '70.0 0.0 0.0    0.0 0.0 1.0 0.0', '60.0 0.0 0.0      0.0 0.0 1.0 0.0',
    '50.0 0.0 0.0      0.0 0.0 1.0 0.0','40.0 0.0 0.0    0.0 0.0 1.0 0.0', '30.0 0.0 0.0      0.0 0.0 1.0 0.0','20.0 0.0 0.0    0.0 0.0 1.0 0.0', '10.0 0.0 0.0      0.0 0.0 1.0 0.0']
positionB2=['-110.0 0.0 0.0    0.0 0.0 0.0 1.0','-100.0 0.0 0.0    0.0 0.0 0.0 1.0', '-90.0 0.0 0.0    0.0 0.0 0.0 1.0', '-80.0 0.0 0.0      0.0 0.0 0.0 1.0', '-70.0 0.0 0.0    0.0 0.0 0.0 1.0', '-60.0 0.0 0.0      0.0 0.0 0.0 1.0',
    '-50.0 0.0 0.0      0.0 0.0 0.0 1.0', '-40.0 0.0 0.0    0.0 0.0 0.0 1.0', '-30.0 0.0 0.0      0.0 0.0 0.0 1.0', '-20.0 0.0 0.0    0.0 0.0 0.0 1.0', '-10.0 0.0 0.0      0.0 0.0 0.0 1.0']
DOF0TransformNode0T=['0.0 0.0 0.0    0.0 0.0 0.0 1.0', '0.0 0.0 0.0    0.0 0.0 0.0 1.0', '0.0 0.0 0.0    0.0 0.0 0.0 1.0', '0.0 0.0 0.0    0.0 0.0 0.0 1.0', '0.0 0.0 0.0    0.0 0.0 0.0 1.0',
    '0.0 0.0 0.0    0.0 0.0 0.0 1.0', '0.0 0.0 0.0    0.0 0.0 0.0 1.0','0.0 0.0 0.0    0.0 0.0 0.0 1.0', '0.0 0.0 0.0    0.0 0.0 0.0 1.0', '0.0 0.0 0.0    0.0 0.0 0.0 1.0']
DOF1TransformNode1T=['0.0 0.0 0.0    0.0 0.0 0.0 1.0', '0.0 0.0 0.0    0.0 0.0 0.0 1.0', '0.0 0.0 0.0    0.0 0.0 0.0 1.0', '0.0 0.0 0.0    0.0 0.0 0.0 1.0', '0.0 0.0 0.0    0.0 0.0 0.0 1.0',
    '0.0 0.0 0.0    0.0 0.0 0.0 1.0', '0.0 0.0 0.0    0.0 0.0 0.0 1.0','0.0 0.0 0.0    0.0 0.0 0.0 1.0', '0.0 0.0 0.0    0.0 0.0 0.0 1.0', '0.0 0.0 0.0    0.0 0.0 0.0 1.0']

def createScene(rootNode):
    rootNode.createObject('RequiredPlugin',pluginName='SoftRobots BeamAdapter SofaPython SofaSparseSolver')
    rootNode.createObject('VisualStyle',displayFlags='showVisualModels showBehaviorModels showCollisionModels hideBoundingCollisionModels showForceFields showInteractionForceFields hideWireframe')
    rootNode.createObject('BackgroundSetting', color=[0.0, 0.0, 0.0, 1.0])
    AnimationManager(rootNode)
    rootNode.dt = stepTime
    rootNode.gravity = [0.0, -gravity, 0.0]
    
    RPCModelNode=rootNode.createChild('RPCModel')
    RPCModelNode.createObject('EulerImplicitSolver', rayleighStiffness=0.0, printLog=False, rayleighMass=0.0)
    RPCModelNode.createObject('SparseLDLSolver', name='ldl')
    RPCModelNode.createObject('GenericConstraintCorrection', name='GCC', solverName='ldl')

    # Creation of the platform
    PlatformNode=RPCModelNode.createChild('Platform')
    # PlatformNode.createObject('EulerImplicitSolver', rayleighStiffness=0.0, printLog=False, rayleighMass=0.0)
    # PlatformNode.createObject('SparseLDLSolver', name='ldl')
    # PlatformNode.createObject('GenericConstraintCorrection', name='GCC', solverName='ldl')
    PlatformNode.createObject('MechanicalObject', name='PlatformStructure', template='Rigid3d', showObject='true', position='0 0 0  0 0 0 1', showObjectScale='0.01')
    PlatformNode.createObject('MeshTopology', edges=0)
    PlatformNode.createObject('PartialFixedConstraint', name='ConstraintPlat', indices='0', fixedDirections='0 0 1 1 1 0')
    PlatformNode.createObject('UniformMass', name='MassPlat', indices='0', totalMass='0.004', showAxisSizeFactor='0.001')

    Beam1Node=RPCModelNode.createChild('Beam1')
    Beam1Node.createObject('MechanicalObject', name='BeamStructure1', template='Rigid3d', showObject='true', position=positionB1, showObjectScale='0.01')
    Beam1Node.createObject('MeshTopology', edges=edgesList)
    Beam1Node.createObject('UniformMass', name='MassAct1', indices='0', totalMass='0.004', showAxisSizeFactor='0.001')
    Beam1Node.createObject('BeamInterpolation', name='interpolation', crossSectionShape='circular', radius='0.125', innerRadius='0.0', 
        dofsAndBeamsAligned='false', defaultYoungModulus='72e9', DOF0TransformNode0=DOF0TransformNode0T, DOF1TransformNode1=DOF1TransformNode1T) 
    Beam1Node.createObject('AdaptiveBeamForceFieldAndMass', massDensity=2201e-9, name='LinkForceField', interpolation='@interpolation', reinforceLength=1)
    Beam1Node.createObject('PartialFixedConstraint', name='ConstraintAct1', indices=0, fixedDirections='1 1 1 1 1 1') 

    Beam2Node=RPCModelNode.createChild('Beam2')
    # Beam2Node.createObject('EulerImplicitSolver', rayleighStiffness=0.0, printLog=False, rayleighMass=0.0)
    # Beam2Node.createObject('SparseLDLSolver', name='ldl')
    # Beam2Node.createObject('GenericConstraintCorrection', name='GCC', solverName='ldl')
    Beam2Node.createObject('MechanicalObject', name='BeamStructure', template='Rigid3d', showObject='true', position=positionB2, showObjectScale='0.01')
    Beam2Node.createObject('MeshTopology', edges=edgesList)
    Beam2Node.createObject('UniformMass', name='MassAct2', indices='0', totalMass='0.004', showAxisSizeFactor='0.001')
    Beam2Node.createObject('BeamInterpolation', name='interpolation', crossSectionShape='circular', radius='0.125', innerRadius='0.0', 
        dofsAndBeamsAligned='false', defaultYoungModulus='72e9', DOF0TransformNode0=DOF0TransformNode0T, DOF1TransformNode1=DOF1TransformNode1T) 
    Beam2Node.createObject('AdaptiveBeamForceFieldAndMass', massDensity=2201e-9, name='LinkForceField', interpolation='@interpolation', reinforceLength=1)
    Beam2Node.createObject('PartialFixedConstraint', name='ConstraintAct2', indices=0, fixedDirections='1 1 1 1 1 1') 

    RPCModelNode.createObject('AttachConstraint',name='ClampingPoint1', object1="@Platform", object2="@Beam1", indices1="0", indices2="10", constraintFactor="1")
    RPCModelNode.createObject('AttachConstraint',name='ClampingPoint2', object1="@Platform", object2="@Beam2", indices1="0",indices2="10", constraintFactor="1")
    return rootNode

Cheers,
Benjamin

[Benjamin]

Hello everyone,

Maybe my post was not as clear as I thought so I created something easier. I created 2 beams. What I want is to connect the first indices of both beams.

I tried two things:

– I tried to use AttachConstraint but it did not work. The two beams were linked but when I moved one the other one did not follow even if there was the link between them.
I tried to take some inspiration from

https://www.sofa-framework.org/community/forum/topic/how-to-do-attachconstraint/

I tried to update the position and rest position of the beam but it did not work.

– I also tried to use BilateralInteractionConstraint as propose in the following post:

https://www.sofa-framework.org/community/forum/topic/best-way-to-implement-ligaments-attaching-finger-bones-together/

But error saying that the determinant of a matrix is too small but I do not see how I can solve that

You can find the BeamMultiCompMoving.py file where both options are implemented at the following link.

Can you help me on this connection between two points of different objects, please ?

Cheers,
Benjamin

[Hugo]

Hi @benjamin

Firstly, sorry for the late reply.
Regarding your SOFA version, v19.06 starts being really old. You might want to update sometimes soon (using SofaPython3!).

Second, thanks a lot for specifying your first post it helps.

Regarding the AttachConstraint approach, could you share with us all the options (data fields) you gave to the AttachConstraint please?

Regarding the BilateralIC, is the error on the determinant occurring also without the constraint ?

NB : your link on BeamMultiCompMoving.py is not properly set

Best,

Hugo

[Benjamin]

Hi @Hugo,

Thank you for your answer. Understandably, you take time to answer due to the number of people asking for your help.

That’s true that my SOFA version starts being old. You gave me that version (containing the soft robot plugin and the inverse plugin) on a GDR day. If I want to update it, I don’t know if it will be easy to get the same plugins.

– For the AttachConstraint and the BilaterallC, I’m using the same data fields (two objects (beams) with only one index).
– If I commented on the BilaterallC, there is no error, just a warning about the denominator threshold of GC that is reached at the first iteration.

Sorry for the link to GitHub. So here is the code that I’m using.

import Sofa
from splib.constants import Key
# mm and kg units
# BEAM_LENGTH = 100  # in mm
# BEAM_RADIUS = 0.5  # in mm
# BEAM_DENSITY = 6.45e-6  # in kg/mm^3
# YOUNG_MODULUS = 160e6  # 160 GPa = 160e9 Pa = 160e9 N/(m.s^2) = 160e6 N/(mm.s^2)
# ELEMENT_NUMBER = 10
# GRAVITATIONAL_ACCELERATION = 9806.65  # in mm/s^2
# STEP_TIME = 0.001  # in s
# LOAD_MASS = 0.5  # in kg

# SI Units
BEAM_LENGTH = 0.1  # in m
BEAM_RADIUS = 0.5e-3  # in m
BEAM_DENSITY = 6450  # in kg/m^3
YOUNG_MODULUS = 160e9  # in Pa
ELEMENT_COUNT = 5
GRAVITATIONAL_ACCELERATION = 9.80665  # in m/s^2
STEP_TIME = 0.01  # in s
LOAD_MASS = 0.1  # in kg

# OpenGL
SLICE_COUNT = 100
SECTOR_COUNT = 5

def createScene(rootNode):
    print('---------- Entering createScene ----------')

    rootNode.createObject('RequiredPlugin',pluginName='SoftRobots BeamAdapter SofaPython SofaSparseSolver')

    rootNode.createObject('VisualStyle',
                          displayFlags='showVisualModels showBehaviorModels showCollisionModels hideBoundingCollisionModels showForceFields showInteractionForceFields hideWireframe')

    rootNode.createObject('BackgroundSetting', color=[0.0, 0.0, 0.0, 1.0])

    rootNode.dt = STEP_TIME
    rootNode.gravity = [0.0, -GRAVITATIONAL_ACCELERATION, 0.0]
    rootNode.createObject('GenericConstraintSolver', name="genericConstraintSolver1", maxIterations="10000", tolerance="1e-12")
    rootNode.createObject('FreeMotionAnimationLoop')
    # rootNode.createObject('DefaultPipeline', verbose="0")
    # rootNode.createObject('BruteForceDetection', name="N2")
    # rootNode.createObject('DefaultContactManager', response="FrictionContact", responseParams="mu=0.6")
    # rootNode.createObject('LocalMinDistance', name="Proximity", alarmDistance="0.05", contactDistance="0.001", angleCone="0.00")
    rootNode.createObject('OglSceneFrame', style="Arrows", alignment="TopRight")
    # rootNode.createObject('DefaultPipeline', verbose="0")
    # rootNode.createObject('BruteForceDetection', name="N2")
    # rootNode.createObject('DefaultPipeline', verbose="0")
    rootNode.createObject('PythonScriptController', classname="Controller")
    
    ##### Compute the beam topology using python
    edgesList = []
    for i in range(0, ELEMENT_COUNT):
        edgesList.append([i, i + 1])

    print("edgesList = ", edgesList)

    positionsList = []
    for i in range(0, ELEMENT_COUNT + 1):
        dx = BEAM_LENGTH / ELEMENT_COUNT
        positionsList.append([dx * i, 0.0, 0.0, 0.0, 0.0, 0.0, 1])

    print("positionsList = ", positionsList)

    wall = rootNode.createChild('wall')
    wall.createObject('EulerImplicitSolver', name="odesolverwall", rayleighStiffness="0.1", rayleighMass="0.1")
    # wall.createObject('SparseLDLSolver')
    wall.createObject('CGLinearSolver', iterations="100", tolerance="1.0e-9", threshold="1.0e-9" )
    # wall.createObject('BTDLinearSolver', printLog="false", verbose="false")
    wall.createObject('MechanicalObject', name='wallFrame', template="Rigid3", position=['0.0 -0.05 0.0  0.0 0.0 0.707 0.707','0.0 0.0 0.0  0.0 0.0 0.707 0.707','0.0 0.05 0.0  0.0 0.0 0.707 0.707'])
    wall.createObject('MeshTopology', edges=['0 1','1 2'])
    # beam.createObject('BeamInterpolation', name='interpolation', crossSectionShape='rectangular', lengthY=2*BEAM_RADIUS, lengthZ=2*BEAM_RADIUS, defaultYoungModulus=YOUNG_MODULUS)
    wall.createObject('BeamInterpolation', name='wallinter', crossSectionShape='circular', radius=BEAM_RADIUS,
                      innerRadius=0.0, defaultYoungModulus=YOUNG_MODULUS, dofsAndBeamsAligned=1)
    wall.createObject('AdaptiveBeamForceFieldAndMass', name='wallForceField', computeMass=1, massDensity=BEAM_DENSITY)
    wall.createObject('PartialFixedConstraint', name='clampingProxEnd', indices=['0','1','2'], fixedDirections='0 1 1 1 1 1')
    wall.createObject('UniformMass', totalMass="0.1", showAxisSizeFactor="0.01")
    # wall.createObject('LinearSolverConstraintCorrection')
    wall.createObject('UncoupledConstraintCorrection')
    
    ##### Beam model
    beam = rootNode.createChild('beam')
    beam.createObject('EulerImplicitSolver', name="odesolverbeam", rayleighStiffness="0.1", rayleighMass="0.1")
    beam.createObject('CGLinearSolver', iterations="100", tolerance="1.0e-9", threshold="1.0e-9" )
    # beam.createObject('SparseLDLSolver')
    # beam.createObject('BTDLinearSolver', printLog="false", verbose="false")
    beam.createObject('MechanicalObject', name='frame', template="Rigid3", position=positionsList)
    beam.createObject('MeshTopology', edges=edgesList)
    # beam.createObject('BeamInterpolation', name='interpolation', crossSectionShape='rectangular', lengthY=2*BEAM_RADIUS, lengthZ=2*BEAM_RADIUS, defaultYoungModulus=YOUNG_MODULUS)
    beam.createObject('BeamInterpolation', name='interpolation', crossSectionShape='circular', radius=BEAM_RADIUS,
                      innerRadius=0.0, defaultYoungModulus=YOUNG_MODULUS, DOF0TransformNode0=['0.0 0.0 0.0  0.0 0.0 0.0 0.0','0.0 0.0 0.0  0.0 0.0 0.0 0.0','0.0 0.0 0.0  0.0 0.0 0.0 0.0'],
                      DOF1TransformNode1=['0.0 0.0 0.0  0.0 0.0 0.0 0.0','0.0 0.0 0.0  0.0 0.0 0.0 0.0','0.0 0.0 0.0  0.0 0.0 0.0 0.0'], dofsAndBeamsAligned=1.0)
    beam.createObject('AdaptiveBeamForceFieldAndMass', name='BeamForceField', computeMass=1, massDensity=BEAM_DENSITY)
    beam.createObject('UniformMass', totalMass="0.1", showAxisSizeFactor="0.01" )
    # beam.createObject('LinearSolverConstraintCorrection')
    beam.createObject('UncoupledConstraintCorrection')
    ##### Clamping
    # beam.createObject('FixedConstraint', name='clamping', indices='0')
    beam.createObject('PartialFixedConstraint', name='clampingProxEnd', indices='0', fixedDirections='0 1 1 1 1 1')
    beam.createObject('FixedConstraint', name='clampingDisEnd', indices=ELEMENT_COUNT)

    # rootNode.createObject('AttachConstraint', name='constforPal', object1="@wall", object2="@beam", indices1="0", indices2="0", constraintFactor="1",twoWays='1')
    # rootNode.createObject('BilateralInteractionConstraint', template="Rigid3", object1="@wall", object2="@beam", first_point="1", second_point="0")
    
    print('---------- Exiting createScene ----------')

    return rootNode

class Controller(Sofa.PythonScriptController):

    def createGraph(self, node):
        print('---------- Entering createGraph ----------')
        self.node = node
        print('---------- Exiting createGraph ----------')

    def reset(self):
        print('---------- Entering reset ----------')
        self.time = 0.0
        print('---------- Exiting reset ----------')
        return 0

    def onBeginAnimationStep(self, dt):
        return 0

    def onEndAnimationStep(self, dt):
        return 0

    def onKeyPressed(self, key):
        print("Key Pressed")
        if key == Key.U: # mouvements suivant l'axe X
            self.node.wall.wallFrame.rest_position=self.node.wall.wallFrame.position
        return 0

[Author]

Posts