How to extract the global stiffness matrix of a system of multiple beam?

[Matyas-D]

Hi everyone !

I am currently working on a system composed of multiple 1D beams, and I need to extract its global stiffness matrix. I tried to follow the example given in access_stiffness_matrix.py but if I use the "ConstantForceField" the matrix is empty and if I use the "AdaptiveBeamForceFieldAndMass" the K matrix is not built. Please see an example code below.

I also tried to use access_matrix.py but this does give me the global system matrix and not the stiffness matrix.

Can anyone clarify how to proceed?

Thanks,
Matyas

def createScene(rootNode):

    parallelRobot = rootNode.addChild('ParallelRobot')

    parallelRobot.addObject('EulerImplicitSolver')
    linear_solver = parallelRobot.addObject('SparseLDLSolver')

    positions = [[-50,0,0,0,0,0,1],[-25,0,0,0,0,0,1],[0,0,0,0,0,0,1],[0,25,0,0,0,0.7071,0.7071],[0,50,0,0,0,0.7071,0.7071]]
    parallelRobot.addObject('MechanicalObject', name = 'frame',  template = "Rigid3", position = positions)
    parallelRobot.addObject('FixedConstraint', name = 'clamping1',  indices = 0)
    parallelRobot.addObject('FixedConstraint', name = 'clamping2',  indices = 4)

    externalForce = parallelRobot.addObject('ConstantForceField',  name="PointForce",  indices=2,
                                             force=[-30.0, 30.0, 0.0, 0.0, 0.0, 0.0, 1.0], showArrowSize=0.005)

    beam1 = parallelRobot.addChild('Beam1')
    beam1.addObject('MeshTopology', edges = [[0,1],[1,2]])
    beam1.addObject('BeamInterpolation', name = 'interpolation', crossSectionShape = 'circular',  radius = 0.5,
                                       dofsAndBeamsAligned = 'false',
                                       DOF0TransformNode0 = [[0,0,0,0,0,0,1],[0,0,0,0,0,0,1]],
                                       DOF1TransformNode1 = [[0,0,0,0,0,0,1],[0,0,0,0,0,0,1]])
    beamForceField1 = beam1.addObject('AdaptiveBeamForceFieldAndMass', name = 'BeamForceField', computeMass = 1,
                                       massDensity = 9.58e-7, interpolation = '@interpolation', reinforceLength = 1)

    visuNode1 = parallelRobot.addChild('VisuNode1')
    visuNode1.addObject('CylinderGridTopology', name='topology', radius=0.5, n=[5, 5, 100], length=50,
                        axis=[1.0, 0.0, 0.0])
    visuNode1.addObject('MechanicalObject', name='cylinder', src='@topology')
    visuNode1.addObject('AdaptiveBeamMapping', interpolation="@../Beam1/interpolation", input="@../frame",
                        output="@cylinder")
    oglNode = visuNode1.addChild('ogl')
    oglNode.addObject('OglModel')
    oglNode.addObject('IdentityMapping')

    beam2 = parallelRobot.addChild('Beam2')
    beam2.addObject('MeshTopology', edges = [[2,3],[3,4]])
    beam2.addObject('BeamInterpolation', name = 'interpolation', crossSectionShape = 'circular', radius = 0.5,
                                       dofsAndBeamsAligned = 'false',
                                       DOF0TransformNode0 = [[0,0,0,0,0,0.7071,0.7071],[0,0,0,0,0,0,1]],
                                       DOF1TransformNode1 = [[0,0,0,0,0,0,1],[0,0,0,0,0,0,1]])
    beamForceField2 = beam2.addObject('AdaptiveBeamForceFieldAndMass', name = 'BeamForceField', computeMass = 1,
                                       massDensity = 9.58e-7,  interpolation = '@interpolation', reinforceLength = 1)

    rootNode.addObject(MyController(name='Controller',  externalForce=externalForce, linear_solver=linear_solver,
                                    beamForceField = beamForceField1))

    visuNode2 = parallelRobot.addChild('VisuNode2')
    visuNode2.addObject('CylinderGridTopology', name='topology', radius=0.5, n=[5, 5, 100], length=50,
                        axis=[1.0, 0.0, 0.0])
    visuNode2.addObject('MechanicalObject', name='cylinder', src='@topology')
    visuNode2.addObject('AdaptiveBeamMapping', interpolation="@../Beam2/interpolation", input="@../frame",
                        output="@cylinder")
    oglNode = visuNode2.addChild('ogl')
    oglNode.addObject('OglModel')
    oglNode.addObject('IdentityMapping')

    return rootNode


class MyController(Sofa.Core.Controller):

    def __init__(self, *args, **kwargs):

        Sofa.Core.Controller.__init__(self, *args, **kwargs)
        self.externalForce = kwargs["externalForce"]
        self.linear_solver = kwargs.get("linear_solver")
        self.beamForceField = kwargs["beamForceField"]

    def onAnimateEndEvent(self, event):
        print(self.externalForce.assembleKMatrix().toarray())
        print(self.beamForceField)
        print(self.linear_solver.A().toarray())

[alxbilger]

Hi,

in the case of the constant force field: f=constant so K=df/dx=0

For the beam force field, what happens exactly?

[Matyas-D]

Hi @alxbilger,

Thank you for your answer.

The purpose of my question is to characterize the way any node would move when a force is applied somewhere on the system. For example if I change the constant force from f=0 to f=30N all the nodes will have a slithly different position and df/dx =/= 0. I suppose there is a non zero stiffness matrix even at static state.

For the beam force field there is no method called "assembleKMatrix", if I try it anyway SOFA shuts down.

[alxbilger]

I confirm there is a bug in AdaptiveBeamForceFieldAndMass preventing to call assembleKMatrix properly. I am trying to debug it. I'll let you know.

[alxbilger]

sofa-framework/SofaPython3#325 should be a fix allowing you to use assembleKMatrix on a AdaptiveBeamForceFieldAndMass.

Let us know if it works.

[Matyas-D]

Thank you for fixing this. I currently have a ready-to-use version of SOFA and I don't have the material to compile it but I will try later, please don't expect a quick feedback.

[alxbilger]

Sure! If the pull request is merged, a binary will be generated and you won't have to compile. Just replace the SofaPython3 binary by the new one.

[github-actions[bot]]

⚠️ ⚠️ ⚠️
@Matyas-D
We are setting a new forum management policy: topics like this one, in which the last reply is more than 4 months old, will automatically be closed. Without further news in the coming weeks, our bot will proceed to automatic closure, thus keeping this forum clean and fresh 🌱 Thank you for your understanding

[github-actions[bot]]

⚠️ ⚠️ ⚠️
@Matyas-D
In accordance with our forum management policy, the last reply is more than 4 months old and is therefore closed. Our objective is to keep the forum up to date and offer the best support experience.

Please feel free to reopen it if the topic is still active by providing us with an update. Please feel free to open a new thread at any time - we'll be happy to help where and when we can.

[github-actions[bot]]

⚠️ ⚠️ ⚠️
@Matyas-D
In accordance with our forum management policy, the last reply is more than 4 months old and is therefore closed. Our objective is to keep the forum up to date and offer the best support experience.

Please feel free to reopen it if the topic is still active by providing us with an update. Please feel free to open a new thread at any time - we'll be happy to help where and when we can.