Unexpected reduction of DOFs when using SI units (WireRestShape / WireBeamInterpolation)

[MarkLeipi]

Hello everyone,

I am currently working with the BeamAdapter plugin and noticed an unexpected behavior when switching my scene from mm-based units to SI units (meters).

To isolate the issue, I created an example scene, where you can choose between two instruments "sheath" (in mm) and "sheath_SI" (in m)

Here is the scene:

import sys
import Sofa
import Sofa.Core

# [INSTRUMENTEN-DICTIONARY]
DOF_straight = 20
DOF_wire = 10

instrument_specs = {
    "sheath_SI": {
        "length_straight": 0.403,
        "length_wire": 0.02334,
        "innerRadius" : 0.005,
        "radius" : 0.01,
        "YM" : 820e6,
        "density" : 1350,
        "DOFs": DOF_straight,
        "DOFs_2" : DOF_wire,
        "speed" : 0.05,
        "ROI": "0 0 0 1 1 1",
        "spireDia": 0.04050,
        "xtip": 0.001

    },
    "sheath": {
        "length_straight": 403-23.34,
        "length_wire": 23.34,
        "innerRadius" : 5,
        "radius" : 10,
        "YM" : 820,
        "density" : 1350e-9,
        "DOFs": DOF_straight,
        "DOFs_2" : DOF_wire,
        "speed": 20,
        "ROI": "0 0 0 500 500 500",
        "spireDia": 40.50,
        "xtip": 1 
        }
}

# [INSTRUMENTENAUSWAHL]

instrument_choice = "sheath_SI"   # "schleuse", "guide", "needle"
spec = instrument_specs[instrument_choice]


# [SZENE]

def createScene(rootNode):

    rootNode.dt = 0.1
    rootNode.gravity = [0, 0, 0]

    rootNode.addObject(
        'RequiredPlugin',
        pluginName='Sofa.Component.IO.Mesh '
                   'BeamAdapter '
                   'Sofa.Component.Constraint.Projective '
                   'Sofa.Component.LinearSolver.Direct '
                   'Sofa.Component.ODESolver.Backward '
                   'Sofa.Component.StateContainer '
                   'Sofa.Component.Topology.Container.Constant '
                   'Sofa.Component.Topology.Container.Grid '
                   'Sofa.Component.Visual '
                   'Sofa.Component.SolidMechanics.Spring ' 
                   'Sofa.Component.Topology.Container.Dynamic'
    )

    rootNode.addObject(
        'VisualStyle',
        displayFlags='showVisualModels showBehaviorModels hideMappings showForceFields'
    )

    rootNode.addObject('DefaultAnimationLoop')
    rootNode.addObject('DefaultVisualManagerLoop')

    # --------------------------------------------------------
    # TOPOLOGIE
    # --------------------------------------------------------

    topoLines = rootNode.addChild('EdgeTopology')

    topoLines.addObject(
        'RodStraightSection',
        name='StraightSection',
        length=spec["length_straight"],
        radius=spec["radius"],
        innerRadius = spec["innerRadius"],
        nbBeams= spec["DOFs"],
        youngModulus=spec["YM"],
        massDensity= spec["density"])
    

    topoLines.addObject("RodSpireSection",
        name='SpireSection',
        length=spec["length_wire"],
        radius=spec["radius"],
        innerRadius = spec["innerRadius"],
        nbBeams= spec["DOFs_2"],
        youngModulus= spec["YM"],
        massDensity= spec["density"],
        spireDiameter = spec["spireDia"],
        spireHeight = 0)
    
    topoLines.addObject(
        'WireRestShape',
        name='BeamRestShape',
        template='Rigid3d',
        wireMaterials='@StraightSection @SpireSection'
    )

    topoLines.addObject('EdgeSetTopologyContainer', name='meshLinesBeam')
    topoLines.addObject('EdgeSetTopologyModifier')
    topoLines.addObject('EdgeSetGeometryAlgorithms', template='Rigid3d')
    topoLines.addObject('MechanicalObject', name='dofTopo2', template='Rigid3d')

    # --------------------------------------------------------
    # BEAMMECHANICS
    # --------------------------------------------------------

    BeamMechanics = rootNode.addChild('BeamModel', bbox = spec["ROI"])

    BeamMechanics.addObject(
        'EulerImplicitSolver',
        rayleighStiffness=0.2,
        rayleighMass=0.1
    )

    BeamMechanics.addObject('BTDLinearSolver')
 
    BeamMechanics.addObject(
        'RegularGridTopology',
        nx=DOF_straight + DOF_wire + 1, ny=1, nz=1,
        xmin=0, xmax=0,
        ymin=0, ymax=0,
        zmin=0, zmax=0,
    )

    BeamMechanics.addObject(
        'MechanicalObject',
        name='DOFs Container',
        template='Rigid3d',
        ry=-90,
        showIndices=True,
        showIndicesScale=0.01
    )

    BeamMechanics.addObject(
        'WireBeamInterpolation',
        name='BeamInterpolation',
        WireRestShape='@../EdgeTopology/BeamRestShape'
    )

    BeamMechanics.addObject(
        'AdaptiveBeamForceFieldAndMass',
        interpolation='@BeamInterpolation'
    )

    BeamMechanics.addObject(
        'InterventionalRadiologyController',
        name='DeployController',
        template='Rigid3d',
        instruments='BeamInterpolation',
        startingPos=[0, 0, 0, 0, 0, 0, 1],
        xtip=spec["xtip"],
        rotationInstrument=0,        
        speed=spec["speed"],      
        listening=True,
        controlledInstrument=0
    )

    BeamMechanics.addObject(
        'FixedProjectiveConstraint',
        indices=0
    )

    BeamMechanics.addObject(
        'RestShapeSpringsForceField',
        points='@DeployController.indexFirstNode',
        angularStiffness=1e8,
        stiffness=1e8
    )

# [MAIN]

def main():
    import SofaRuntime
    import Sofa.Gui

    root = Sofa.Core.Node('root')
    createScene(root)
    Sofa.Simulation.init(root)

    Sofa.Gui.GUIManager.Init('deploy_scene', 'qglviewer')
    Sofa.Gui.GUIManager.createGUI(root, __file__)
    Sofa.Gui.GUIManager.SetDimension(1080, 1080)
    Sofa.Gui.GUIManager.MainLoop(root)
    Sofa.Gui.GUIManager.closeGUI()

if __name__ == '__main__':
    main()

What i observed for the case of SI-units ("sheath_SI"):

Increasing nbBeams for the RodSpireSection does not lead to a corresponding increase in DOFs. Regardless of how large I choose nbBeams, only up to 3 DOFs are created.

Inspecting the scene graph shows that in the length_list of WireBeamInterpolation, the beam lengths of each beam seem to never get below ~0.01, even though it should do so, because of the small total length.

In contrast, when using mm ("sheath") units, the discretization behaves as expected.

My question:

Is there some sort of threshhold I can adapt to avoid this behavior, or is there another way to correctly model a beam with small lengths like in this example?

Thanks for any kind of support :)

[MarkLeipi]

Update/Solution:

In the InterventionalRadiologyController there is a parameter threshold which is homogeneous to the length unit used in the simulation.
The default value is: threshold = 0.01
When using SI units (meters), this corresponds to a minimum allowed beam segment length of 1 cm.

Setting the threshold to a smaller SI-consistent value fixed my issue:
threshold = 0.01e-3 # 10 µm

[hugtalbot]

Hey @MarkLeipi
Thanks for sharing 👏