How to use BeamAdapter with a self intersecting collision model?

[fwhotto]

I want to develop a scene where I have a beam being deployed as in examples/SingleBeamDeployment.scn. But in addition to collision with the carotids/vascalature I want simulate the beam colliding with itself. As I understand one can add a linecollision model with the selfCollision flag. But the behaviour is not as I would expect, the self intersection is not active. If you have any directions it would be much appreciated.
Current scene

<?xml version="1.0"?>
<Node name="root" gravity="0 0 0" dt="0.01" bbox="0 0 0 10 10 10">
    <RequiredPlugin name="BeamAdapter"/> <!-- Needed to use components [WireBeamInterpolation WireRestShape] -->
    <RequiredPlugin name="Sofa.Component.Constraint.Projective"/> <!-- Needed to use components [FixedProjectiveConstraint] -->
    <RequiredPlugin name="Sofa.Component.LinearSolver.Direct"/> <!-- Needed to use components [BTDLinearSolver] -->
    <RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> <!-- Needed to use components [EulerImplicitSolver] -->
    <RequiredPlugin name="Sofa.Component.SolidMechanics.Spring"/> <!-- Needed to use components [RestShapeSpringsForceField] -->
    <RequiredPlugin name="Sofa.Component.StateContainer"/> <!-- Needed to use components [MechanicalObject] -->
    <RequiredPlugin name="Sofa.Component.Topology.Container.Dynamic"/> <!-- Needed to use components [EdgeSetGeometryAlgorithms EdgeSetTopologyContainer EdgeSetTopologyModifier] -->
    <RequiredPlugin name="Sofa.Component.Topology.Container.Grid"/> <!-- Needed to use components [RegularGridTopology] -->
    <RequiredPlugin name="Sofa.Component.Visual"/> <!-- Needed to use components [VisualStyle] -->
    <RequiredPlugin name="Sofa.Component.IO.Mesh"/> <!-- Needed to use components [MeshVTKLoader] -->
    <RequiredPlugin name="Sofa.Component.Collision.Geometry"/>
    <RequiredPlugin name="Sofa.Component.Collision.Detection.Algorithm"/>
    <RequiredPlugin name="Sofa.Component.Collision.Detection.Intersection"/>
    <RequiredPlugin name="Sofa.Component.Collision.Response.Contact"/>

    <VisualStyle displayFlags="showVisualModels showBehaviorModels showCollisionModels hideMappings showForceFields" />

    <DefaultAnimationLoop />
    <DefaultVisualManagerLoop />

    <!-- Collision pipeline -->
    <DefaultPipeline depth="15" verbose="0" draw="0" />
    <BruteForceBroadPhase name="N2" />
    <BVHNarrowPhase/>
    <MinProximityIntersection name="Proximity" alarmDistance="3." contactDistance="0.5" />
    <DefaultContactManager name="Response" response="FrictionContactConstraint" />

    <Node name="EdgeTopology">
        <MeshOBJLoader name="loader" filename="polyline.obj" />

        <RodMeshSection name="ComplexPreshape" youngModulus="20000" radius="0.9" massDensity="0.00000155" />
        
        <WireRestShape template="Rigid3d" name="BeamRestShape" printLog="1" wireMaterials="@ComplexPreshape"/>
        <EdgeSetTopologyContainer name="meshLinesBeam"/>
        <EdgeSetTopologyModifier name="Modifier"/> 
        <EdgeSetGeometryAlgorithms name="GeomAlgo" template="Rigid3d" /> 
        <MechanicalObject name="dofTopo2" template="Rigid3d" />
    </Node>  
    
    <Node name="BeamModel">
        <EulerImplicitSolver rayleighStiffness="0.2" rayleighMass="0.1" printLog="false" />
        <BTDLinearSolver verbose="0"/>
        <RegularGridTopology name="MeshLines" nx="60" ny="1" nz="1"
                                    xmax="0.0" xmin="0.0" ymin="0" ymax="0" zmax="0" zmin="0"
                                    p0="0 0 0" drawEdges="1"/>
        <MechanicalObject template="Rigid3d" name="DOFs Container" ry="-90" /> 

        <WireBeamInterpolation name="BeamInterpolation" WireRestShape="@../EdgeTopology/BeamRestShape" />
        <AdaptiveBeamForceFieldAndMass name="BeamForceField" interpolation="@BeamInterpolation"/>
        <InterventionalRadiologyController name="DeployController" template="Rigid3d" instruments="BeamInterpolation"
                                    startingPos="0 0 0 0 0 0 1" xtip="0 0 0" printLog="1" 
                                    rotationInstrument="0 0 0" step="0.5" speed="2.0" 
                                    listening="1" controlledInstrument="0"/>
        <FixedProjectiveConstraint name="FixedConstraint" indices="0" />
        <RestShapeSpringsForceField name="RestSPForceField" points="@DeployController.indexFirstNode" angularStiffness="1e8" stiffness="1e8"/> 

        <!-- SIMPLE self-collision for the rod -->
        <Node name="CollisionModel">
            <EdgeSetTopologyContainer name="topo"/>
            <EdgeSetTopologyModifier />
            <EdgeSetGeometryAlgorithms />
            <MechanicalObject template="Vec3d" name="x"/>
            <!-- Enable self-collision on the line model -->
            <LineCollisionModel selfCollision="1" proximity="0.0"/>
            <!-- Map the beam frames to the Vec3 collision points/lines -->
            <BeamLinearMapping template="Rigid3d,Vec3d" input="@../mstate" output="@x"/>
        </Node>


    </Node>  
    
</Node>

[bakpaul]

Hello !

You are right, adding the parameter 'selfCollision="true"' should do the trick. What seems to be the issue here is that you are setting a proximity to 0.0 on the collision model. Along with some final solving error, because the constraint solver might not be converging, it is easy to cross the beam as if nothing was there... Given that your rod seems to have a radius of 0.9, why not trying this with a proximity of the same size ?

You could also try to change the MinProximityIntersection by a NewProximityIntersection. Not all the intersector support all pairs of elements, the one I've mentioned is the most complete one to date.

[fwhotto]

I also wrote a scene without dependencies if its easier for you to understand what is going on this way?

<?xml version="1.0"?>
<Node name="root" gravity="0 0 0" dt="0.01" bbox="0 0 0 10 10 10">
    <RequiredPlugin name="BeamAdapter"/> <!-- Needed to use components [WireBeamInterpolation WireRestShape] -->
    <RequiredPlugin name="Sofa.Component.Constraint.Projective"/> <!-- Needed to use components [FixedProjectiveConstraint] -->
    <RequiredPlugin name="Sofa.Component.LinearSolver.Direct"/> <!-- Needed to use components [BTDLinearSolver] -->
    <RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> <!-- Needed to use components [EulerImplicitSolver] -->
    <RequiredPlugin name="Sofa.Component.SolidMechanics.Spring"/> <!-- Needed to use components [RestShapeSpringsForceField] -->
    <RequiredPlugin name="Sofa.Component.StateContainer"/> <!-- Needed to use components [MechanicalObject] -->
    <RequiredPlugin name="Sofa.Component.Topology.Container.Dynamic"/> <!-- Needed to use components [EdgeSetGeometryAlgorithms EdgeSetTopologyContainer EdgeSetTopologyModifier] -->
    <RequiredPlugin name="Sofa.Component.Topology.Container.Grid"/> <!-- Needed to use components [RegularGridTopology] -->
    <RequiredPlugin name="Sofa.Component.Visual"/> <!-- Needed to use components [VisualStyle] -->
    <RequiredPlugin name="Sofa.Component.IO.Mesh"/> <!-- Needed to use components [MeshVTKLoader] -->
    <RequiredPlugin name="Sofa.Component.Collision.Geometry"/>
    <RequiredPlugin name="Sofa.Component.Collision.Detection.Algorithm"/>
    <RequiredPlugin name="Sofa.Component.Collision.Detection.Intersection"/>
    <RequiredPlugin name="Sofa.Component.Collision.Response.Contact"/>
    <RequiredPlugin name="Sofa.Component.SceneUtility"/>
    <RequiredPlugin name="Sofa.Component.Mapping.Linear"/> <!-- Needed to use components [BeamLinearMapping] --> 
    <RequiredPlugin name="Sofa.Component.AnimationLoop"/> <!-- Needed to use components [FreeMotionAnimationLoop] -->  
    <RequiredPlugin name="Sofa.Component.Constraint.Lagrangian.Solver"/> <!-- Needed to use components [GenericConstraintSolver] -->  


    <VisualStyle displayFlags="showAll showCollisionModels" />

    <FreeMotionAnimationLoop />
    <DefaultVisualManagerLoop />
    <GenericConstraintSolver maxIterations="1000" tolerance="1e-3" printLog="1"/>

    <!-- Collision pipeline -->
    <CollisionPipeline name="pipeline" printLog="1"/>
    <CollisionResponse name="response" response="FrictionContactConstraint" responseParams="mu=0.0" printLog="1" />
    <BruteForceBroadPhase name="N2" />
    <BVHNarrowPhase/>
    <NewProximityIntersection alarmDistance="0.1" contactDistance="0.01" printLog="1"/>


    <Node name="EdgeTopology">
        <RodStraightSection name="StraightSection" youngModulus="20000" radius="0.9" massDensity="0.00000155" nbEdgesCollis="30" nbEdgesVisu="196" length="980.0"/>
        <RodSpireSection name="SpireSection" youngModulus="20000" radius="0.9" massDensity="0.00000155" nbEdgesCollis="5" nbEdgesVisu="4" length="20.0" spireDiameter="25" spireHeight="0.0"/>
        
        <WireRestShape template="Rigid3d" name="BeamRestShape" printLog="1" wireMaterials="@StraightSection @SpireSection"/>
        <EdgeSetTopologyContainer name="meshLinesBeam"/>
        <EdgeSetTopologyModifier name="Modifier"/> 
        <EdgeSetGeometryAlgorithms name="GeomAlgo" template="Rigid3d" /> 
        <MechanicalObject name="dofTopo2" template="Rigid3d" />
    </Node>  
    
    <Node name="BeamModel">
        <EulerImplicitSolver rayleighStiffness="0.2" rayleighMass="0.1" printLog="false" />
        <BTDLinearSolver verbose="0"/>
        <RegularGridTopology name="MeshLines" nx="60" ny="1" nz="1"
                                    xmax="0.0" xmin="0.0" ymin="0" ymax="0" zmax="0" zmin="0"
                                    p0="0 0 0" drawEdges="1"/>
        <MechanicalObject template="Rigid3d" name="DOFs Container" ry="-90" /> 

        <WireBeamInterpolation name="BeamInterpolation" WireRestShape="@../EdgeTopology/BeamRestShape" />
        <AdaptiveBeamForceFieldAndMass name="BeamForceField" interpolation="@BeamInterpolation"/>
        <InterventionalRadiologyController name="DeployController" template="Rigid3d" instruments="BeamInterpolation"
                                    startingPos="0 0 0 0 0 0 1" xtip="0 0 0" printLog="1" 
                                    rotationInstrument="0 0 0" step="0.5" speed="0.5" 
                                    listening="1" controlledInstrument="0"/>
        <FixedProjectiveConstraint name="FixedConstraint" indices="0" />
        <RestShapeSpringsForceField name="RestSPForceField" points="@DeployController.indexFirstNode" angularStiffness="1e8" stiffness="1e8"/> 

        <!-- SIMPLE self-collision for the rod -->
        <Node name="CollisionModel">

            <!-- edge topology for lines (fed from the same loader as the rest shape) -->
            <EdgeSetTopologyContainer name="topo"
                                    src="@../../EdgeTopology/meshLinesBeam"
                                    checkTopology="1"/> 
            <EdgeSetTopologyModifier/>
            <EdgeSetGeometryAlgorithms/>

            <MechanicalObject template="Vec3d" name="x"
                            src="@../../EdgeTopology/meshLinesBeam"/> 

            <!-- Enable self-collision on the line model -->
            <LineCollisionModel selfCollision="1" proximity="0.9" contactFriction="0.0" topology="@topo" printLog="1"/>

            <!-- Map the beam frames to the Vec3 collision points/lines -->
            <BeamLinearMapping template="Rigid3d,Vec3d" input="@../DOFs_Container" output="@x" printLog="1"/>
        </Node>


    </Node>  
    
</Node>

[bakpaul]

I will summon the two main contributors/users of BeamAdapter : @fredroy and @epernod do you have an idea on that ??

[fredroy]

I dont think the selfCollision would bode well with the "InterventionalRadiologyController" component.
The principle behind this controller is to deploy the DOFs from a starting position, so you deploy lets say 10% of your wire, 90% of the DOFs of the wire will be positioned on startPos (which is 0 0 0 in your caase).
So if the collision detection is allowed to detect collision on the same object, it will find that those 90% are in collision, as the distance between them is < of your alarmDistance. So we can presume it will be messy after this 🙃

I guess the trick would be to disable collision detection on self for the DOFs which are not deployed. It seems there is some code about it in IBeamAdapter, as WireBeamInterpolation has a beamCollision field, where the "deployed" edges are listed (apparently). So ideally, one would want to supply this list of edges to the list of edges in the topology in the collision node (in your case `CollisionModel/topo) but I did not try it 😅

[fwhotto]

Hi Frederick,

Thank you for the prompt reply. That makes sense. As I understand this field is not exposed in the .scn files? Would it be possible to try this using SofaPy or it would require forking and modifying the beamadapter plugin and building sofa from source?

[fredroy]

Hi Frederick,

Thank you for the prompt reply. That makes sense. As I understand this field is not exposed in the .scn files? Would it be possible to try this using SofaPy or it would require forking and modifying the beamadapter plugin and building sofa from source?

Actually this field should be exposed, as it is a "Data".

There was (and still is) a project to convert all the code in the InterventionalRadiologyController component into a python3 script, as its code is an aggregate of C++ for the last 15 years 👀 (student code over student code over.....).
If you want a "more" modern code/approach, you may use the Cosserat plugin which was started much more recently and is using heavily python3. But there is not the same level of "features" so it may not be possible in your case.

[github-actions[bot]]

⚠️ ⚠️ ⚠️
@fwhotto
Feedback has been given to you by the project reviewers, however we have not received a response from you. Without further news in the coming weeks, this discussion will be automatically closed in order to keep this forum clean and fresh 🌱 Thank you for your understanding

[github-actions[bot]]

⚠️ ⚠️ ⚠️
@fwhotto
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.