Problem with soft and rigid body collision

[Thanit09]

Hello everyone,

I am new to SOFA, and my goal is to simulate the interaction between a rigid and soft body, and to extract the surface topology of the soft body during the interaction (think of an object lying on a mattress).

To start off I tried to modify the liver example of SOFA python.
I replaced the sphere with an obj model of a sphere, which is part of SOFA, and defined it as rigid. Additionally the sphere collision model of the liver was replaced by an Triangle/Point/LineCollisionModel, same for the sphere.
Since I need an intersection detection between two meshes I added the MinProximityIntersection object.

Running code can be seen here:

import Sofa

USE_GUI = True


def main():
    import SofaRuntime
    import Sofa.Gui

    root = Sofa.Core.Node("root")
    createScene(root)
    Sofa.Simulation.init(root)
    SofaRuntime.importPlugin("Sofa.Component.Collision.Detection.Algorithm")

    if not USE_GUI:
        for iteration in range(10):
            Sofa.Simulation.animate(root, root.dt.value)
    else:
        Sofa.Gui.GUIManager.Init("myscene", "qglviewer")
        Sofa.Gui.GUIManager.createGUI(root, __file__)
        Sofa.Gui.GUIManager.SetDimension(1080, 1080)
        Sofa.Gui.GUIManager.MainLoop(root)
        Sofa.Gui.GUIManager.closeGUI()


def createScene(root):
    root.gravity=[0, -9.81, 0]
    root.dt=0.02

    root.addObject('DefaultVisualManagerLoop')
    root.addObject('DefaultAnimationLoop')

    root.addObject('VisualStyle', displayFlags="showCollisionModels showVisualModels showForceFields")
    root.addObject('RequiredPlugin', pluginName="SofaImplicitOdeSolver SofaLoader SofaOpenglVisual SofaBoundaryCondition SofaGeneralLoader SofaGeneralSimpleFem Sofa.Component.Collision.Detection.Algorithm Sofa.Component.Collision.Detection.Intersection Sofa.Component.LinearSolver.Iterative Sofa.Component.Mass") 
    root.addObject('DefaultPipeline', name="CollisionPipeline")
    root.addObject('BruteForceDetection', name="N2")
    root.addObject('MinProximityIntersection', useLineLine=True, usePointPoint=True, useLinePoint=True)
    #root.addObject('LocalMinDistance', name='proximity', alarmDistance=0.3, contactDistance=0.15, angleCone=0.0)
    root.createObject('DefaultContactManager', name="CollisionResponse", response="PenalityContactForceField")
    root.addObject('DiscreteIntersection')

    root.addObject('MeshObjLoader', name="LiverSurface", filename="mesh/liver-smooth.obj")
    root.addObject('MeshObjLoader', name="SphereSurface", filename="mesh/ball.obj", triangulate="true", scale=1.0, translation='0 10 0')
    liver = root.addChild('Liver')
    liver.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness=0.1, rayleighMass=0.1)
    liver.addObject('CGLinearSolver', name="linear_solver", iterations=25, tolerance=1e-09, threshold=1e-09)
    liver.addObject('MeshGmshLoader', name="meshLoader", filename="mesh/liver.msh")
    liver.addObject('TetrahedronSetTopologyContainer', name="topo", src="@meshLoader")
    liver.addObject('MechanicalObject', name="dofs", src="@meshLoader")
    liver.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="GeomAlgo")
    liver.addObject('DiagonalMass', name="Mass", massDensity=1.0)
    liver.addObject('TetrahedralCorotationalFEMForceField', template="Vec3d", name="FEM", method="large", poissonRatio=0.3, youngModulus=3000, computeGlobalMatrix=False)
    liver.addObject('FixedConstraint', name="FixedConstraint", indices="3 39 64")

    visu = liver.addChild('Visu')
    visu.addObject('OglModel', name="VisualModel", src="@../../LiverSurface")
    visu.addObject('BarycentricMapping', name="VisualMapping", input="@../dofs", output="@VisualModel")

    surf = liver.addChild('Surf')
    #surf.addObject('MeshSTLLoader', name="LiverSurface", filename="/resources/Axle_long_vertical.STL")COPY meshes/Axle_long_vertical_small.stl /resources/Axle_long_vertical_small.STL
    surf.addObject('MeshTopology', name = 'Axle_long', src="@../../LiverSurface")
    #surf.addObject('MechanicalObject')
    surf.addObject('MechanicalObject', name="spheres", position="@../../LiverSurface.position", src='@../../LiverSurface')
    surf.addObject('TriangleCollisionModel', name='triangleLoader',contactStiffness='100', selfCollision='0', group='1', moving=True, simulated=True)
    surf.addObject('LineCollisionModel', moving=True, simulated=True)
    surf.addObject('PointCollisionModel', moving=True, simulated=True)
    surf.addObject('BarycentricMapping', name="CollisionMapping", input="@../dofs", output="@spheres")

    root.addObject(KeyPressedController(name = "SphereCreator"))

    return root


class KeyPressedController(Sofa.Core.Controller):
    """ This controller monitors new sphere objects.
    Press ctrl and the L key to make spheres falling!
    """
    def __init__(self, *args, **kwargs):
        Sofa.Core.Controller.__init__(self, *args, **kwargs)
        self.iteration = 0

    def onKeypressedEvent(self, event):
        # Press L key triggers the creation of new objects in the scene
        if event['key']=='L':
            self.createNewSphere()
            
    def createNewSphere(self):
        root = self.getContext()
        newSphere = root.addChild('FallingSphere-'+str(self.iteration))
        newSphere.addObject('EulerImplicitSolver')
        newSphere.addObject('CGLinearSolver', threshold='1e-09', tolerance='1e-09', iterations='200')
        newSphere.addObject('MeshTopology', src="@../SphereSurface")
        MO = newSphere.addObject('MechanicalObject', showObject=True, position=[-2, 600+self.iteration, 0, 0, 0, 0, 1], name=f'Particle-{self.iteration}', template='Rigid3d')
        Mass = newSphere.addObject('UniformMass', totalMass=1)
        Force = newSphere.addObject('ConstantForceField', name="CFF", totalForce=[0, -1, 0, 0, 0, 0] )
        #Sphere = newSphere.addObject('SphereCollisionModel', name="SCM", radius=1.0, group='1')
        
        surf = newSphere.addChild('Surf')
        surf.addObject('MeshTopology', name = 'sphere', src="@../../SphereSurface")
        surf.addObject('MechanicalObject', name=f"spheres-{self.iteration}", position=f"@../Particle-{self.iteration}.position", src='@../../SphereSurface')
        surf.addObject('TriangleCollisionModel', name='triangleLoader',contactStiffness='100', selfCollision='0', group='1', moving=True, simulated=True)
        surf.addObject('LineCollisionModel', moving=True, simulated=True)
        surf.addObject('PointCollisionModel', moving=True, simulated=True)
        surf.addObject('RigidMapping', name=f"CollisionMapping-{self.iteration}", input=f"@../Particle-{self.iteration}", output=f"@spheres-{self.iteration}")


        newSphere.init()
        self.iteration = self.iteration+1


# Function used only if this script is called from a python environment
if __name__ == '__main__':
    main()

As soon as the sphere approaches the liver a long orange line appears extending upwards from the sphere. With continuation of the simulation the liver starts deforming and at some point the meshes disappear first partially (you can see the tetrahedral mesh inside the liver) and then completely and you have to zoom out to see them again, even though the objects don't appear to be that much larger.

I assume I have done something wrong with the collision model.
Could someone please give me pointers on how to solve this problem?

Thank you very much for your help!

[hugtalbot]

Hi @Thanit09

So sorry for the delay, most of the team was in summer vacations!

Regarding your rigid sphere, since it is a rigid body, no need to use a full mesh for the physics : all you need is one rigid frame in charge of the rigid motion. You could still map on this frame any mesh as a collision mesh.

The behavior that you noticed with the collision models is due to the fact that you are using penalities and this needs fine tuning.
To avoid this, you could go for a constraint approach to solve the collision detection. But the scene would be different:

• a FreeMotionAnimationLoop should be used
• a ConstaintSolver should be used
• the collision response should be response="FrictionContactConstraint"
• for each of your object a ConstraintCorrection will be needed for the corrective motion to apply

An example of this is available in the FrictionContact example.

In your case, a working example would be:

# Required import for python
import Sofa


# Choose in your script to activate or not the GUI
USE_GUI = True


def main():
    import SofaRuntime
    import Sofa.Gui

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

    if not USE_GUI:
        for iteration in range(10):
            Sofa.Simulation.animate(root, root.dt.value)
    else:
        Sofa.Gui.GUIManager.Init("myscene", "qglviewer")
        Sofa.Gui.GUIManager.createGUI(root, __file__)
        Sofa.Gui.GUIManager.SetDimension(1080, 1080)
        Sofa.Gui.GUIManager.MainLoop(root)
        Sofa.Gui.GUIManager.closeGUI()


def createScene(root):
    root.gravity=[0, -9.81, 0]
    root.dt=0.02

    root.addObject("RequiredPlugin", pluginName=[    'Sofa.Component.Collision.Detection.Algorithm',
    'Sofa.Component.Collision.Detection.Intersection',
    'Sofa.Component.Collision.Geometry',
    'Sofa.Component.Collision.Response.Contact',
    'Sofa.Component.Constraint.Projective',
    'Sofa.Component.IO.Mesh',
    'Sofa.Component.LinearSolver.Iterative',
    'Sofa.Component.Mapping.Linear',
    'Sofa.Component.Mass',
    'Sofa.Component.ODESolver.Backward',
    'Sofa.Component.SolidMechanics.FEM.Elastic',
    'Sofa.Component.StateContainer',
    'Sofa.Component.Topology.Container.Dynamic',
    'Sofa.Component.Visual',
    'Sofa.GL.Component.Rendering3D'
    ])

    root.addObject('VisualStyle', displayFlags="showCollisionModels hideVisualModels showForceFields")

    root.addObject('FreeMotionAnimationLoop')
    root.addObject('GenericConstraintSolver', maxIterations=1000, tolerance=1e-6)

    root.addObject('CollisionPipeline', name="CollisionPipeline")
    root.addObject('BruteForceBroadPhase', name="BroadPhase")
    root.addObject('BVHNarrowPhase', name="NarrowPhase")
    root.addObject('DefaultContactManager', name="CollisionResponse", response="FrictionContactConstraint")
    root.addObject('MinProximityIntersection', useLineLine=True, usePointPoint=True, alarmDistance=0.3, contactDistance=0.15, useLinePoint=True)

    root.addObject('MeshOBJLoader', name="LiverSurface", filename="mesh/liver-smooth.obj")

    liver = root.addChild('Liver')
    liver.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness=0.1, rayleighMass=0.1)
    liver.addObject('SparseLDLSolver', name="linear_solver", template="CompressedRowSparseMatrixMat3x3d")
    liver.addObject('MeshGmshLoader', name="meshLoader", filename="mesh/liver.msh")
    liver.addObject('TetrahedronSetTopologyContainer', name="topo", src="@meshLoader")
    liver.addObject('MechanicalObject', name="dofs", src="@meshLoader")
    liver.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="GeomAlgo")
    liver.addObject('DiagonalMass', name="Mass", massDensity=1.0)
    liver.addObject('TetrahedralCorotationalFEMForceField', template="Vec3d", name="FEM", method="large", poissonRatio=0.3, youngModulus=3000, computeGlobalMatrix=False)
    liver.addObject('FixedConstraint', name="FixedConstraint", indices="3 39 64")
    liver.addObject('LinearSolverConstraintCorrection')

    LiverSurf = liver.addChild('ExtractSurface')
    LiverSurf.addObject('TriangleSetTopologyContainer', name="Container", position="@../topo.position")
    LiverSurf.addObject('TriangleSetTopologyModifier', name="Modifier")
    LiverSurf.addObject('Tetra2TriangleTopologicalMapping', name="SurfaceExtractMapping", input="@../topo", output="@Container")

    LiverCollision = LiverSurf.addChild('Surf')
    LiverCollision.addObject('TriangleSetTopologyContainer', name="Container", src="@../Container")
    LiverCollision.addObject('MechanicalObject', name="surfaceDOFs")
    LiverCollision.addObject('PointCollisionModel', name="CollisionModel")
    LiverCollision.addObject('IdentityMapping', name="CollisionMapping", input="@../../dofs", output="@surfaceDOFs")

    newSphere = root.addChild('FallingSphere-0')
    newSphere.addObject('EulerImplicitSolver')
    newSphere.addObject('CGLinearSolver', threshold='1e-09', tolerance='1e-09', iterations='200')
    MO = newSphere.addObject('MechanicalObject', showObject=True, position=[-2, 10, 0, 0, 0, 0, 1], name=f'Particle-0', template='Rigid3d')
    Mass = newSphere.addObject('UniformMass', totalMass=1)
    Force = newSphere.addObject('ConstantForceField', name="CFF", totalForce=[0, -1, 0, 0, 0, 0] )
    Sphere = newSphere.addObject('SphereCollisionModel', name="SCM", radius=1.0 )
    newSphere.addObject('UncoupledConstraintCorrection' )

    return root



# Function used only if this script is called from a python environment
if __name__ == '__main__':
    main()

Hope this helps

PS: remember that we do provide training sessions to get started and progress on SOFA, more in this discussion here

[github-actions[bot]]

⚠️ ⚠️ ⚠️
@Thanit09
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]]

⚠️ ⚠️ ⚠️
@Thanit09
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]]

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