Softbody to softbody collision model breaking into each other when force applied

[thundermms]

Problem

Description
I am new to SOFA and finite element modeling, and I would love to get some input and opinion from you guys!
Basically I am trying to simulate a human muscle,fat and skin (layered on top of each other as attached in image) All of these have a Tetrahedron .msh, visual and collision model setup.
I am trying to get soft body to soft body collision set up but I am facing an issue where the layers tear into another one if i drag it more than an extent or if i drag it fast. ( as shown in the video).
Is this behaviour because I haven't set up the collision model properly yet or is this an expected outcome due to the Poisson ratio and Youngs modulus of each layer.
I have attached my code as well for reference

---

Environment

Context

• System: Ubuntu 20.04
• Version of SOFA: SOFA V24.06.00
• State: Installed from the website

Env vars

python -c "exec( \"import os, sys\nprint('#################')\nprint('--- sys.version ---')\nprint(sys.version)\nprint('--- PATH ---')\ntry:\n  print(os.environ['PATH'])\nexcept Exception:\n  pass\nprint('--- SOFA_ROOT ---')\ntry:\n  print(os.environ['SOFA_ROOT'])\nexcept Exception:\n  pass\nprint('--- PYTHONPATH ---')\ntry:\n  print(os.environ['PYTHONPATH'])\nexcept Exception:\n  pass\nprint('--- sys.path ---')\ntry:\n   print(str(sys.path))\nexcept Exception:\n   pass\nprint('#################')\" )"

--- sys.version ---
3.10.12 (main, Nov  6 2024, 20:22:13) [GCC 11.4.0]
--- PATH ---
/home/sun/sofa/build/bin:/home/sun/.local/bin:/usr/bin/python3:/usr/local/bin:/home/sun/sofa/build/bin:/home/sun/.local/bin:/usr/bin/python3:/usr/local/bin:/opt/ros/humble/bin:/home/sun/.local/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games:/snap/bin:/snap/bin
--- SOFA_ROOT ---
/home/sun/SOFA/v24.06.00:/home/sun/SOFA/v24.06.00:
--- PYTHONPATH ---
/home/sun/SOFA/v24.06.00/plugins/SofaPython3/lib/python3/site-packages:/home/sun/SOFA/v24.06.00/plugins/SofaPython3/lib/python3/site-packages:/opt/ros/humble/lib/python3.10/site-packages:/opt/ros/humble/local/lib/python3.10/dist-packages
--- sys.path ---
['', '/home/sun/SOFA/v24.06.00/plugins/SofaPython3/lib/python3/site-packages', '/opt/ros/humble/lib/python3.10/site-packages', '/opt/ros/humble/local/lib/python3.10/dist-packages', '/usr/lib/python310.zip', '/usr/lib/python3.10', '/usr/lib/python3.10/lib-dynload', '/home/sun/.local/lib/python3.10/site-packages', '/home/sun/.local/lib/python3.10/site-packages/rerun_sdk', '/usr/local/lib/python3.10/dist-packages', '/usr/lib/python3/dist-packages']

---

WhatsApp.Video.2024-11-18.at.11.42.36.AM.mp4

    def createScene(root):
    # root.gravity=[0, -9.81, 0]
    root.dt=0.01
    root.addObject("VisualGrid", nbSubdiv=10, size=10)
    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',
    'Sofa.Component.Constraint.Lagrangian.Correction',
    'Sofa.Component.Constraint.Lagrangian.Solver',
    'Sofa.Component.AnimationLoop'
    ])
    root.addObject('OglSceneFrame', style="Arrows", alignment="TopRight")

    root.addObject('DefaultAnimationLoop')

    root.addObject('VisualStyle', displayFlags="showCollisionModels")
    root.addObject('CollisionPipeline', name="CollisionPipeline")
    root.addObject('BruteForceBroadPhase', name="BroadPhase")
    root.addObject('DefaultContactManager', name="CollisionResponse", response="PenalityContactForceField")
    # root.addObject('DefaultContactManager', name="CollisionResponse", response="FrictionContactConstraint")
    root.addObject('LocalMinDistance', name="LMD-proximity" , alarmDistance=0.1, contactDistance=0.05, useLMDFilters=0)
    root.addObject('DiscreteIntersection')
    root.addObject('CollisionPipeline', verbose="0", draw="0", depth="6")
    root.addObject('BVHNarrowPhase', name="NarrowPhase")
    root.addObject('DefaultAnimationLoop', name="AnimationLoop")
    root.addObject('GenericConstraintSolver', resolutionMethod="false", tolerance="1e-4", maxIterations="250", multithreading="true")


     ######################################################################################################################################################################
    ##
    ##                                          This part of the Scene defines Skin Layer of Tissue Phantom
    ##
    ######################################################################################################################################################################


    ##---------------------------Skin Mechanical Model---------------------------------

    skin = root.addChild('skin')
    
    skin.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness="0.1", rayleighMass="0.1")
    skin.addObject('CGLinearSolver', name="linear_solver", iterations="25", tolerance="1e-09", threshold="1e-09")
    skin.addObject('MeshGmshLoader', name="meshLoader", filename="D:\SOFA Learning\My Learning\mesh\Phantom Surface Meshes\ReducedMesh\Skin1000.msh", triangulate="true", scale=0.1)
    #skin.addObject('MeshGmshLoader', name="meshLoader", filename="D:\SOFA Learning\My Learning\mesh\Phantom Surface Meshes\ReducedMesh\skin600.msh", triangulate="true")
    skin.addObject('TetrahedronSetTopologyContainer', name="topo", src="@meshLoader")
    skin.addObject('MechanicalObject', name="dofs", src="@meshLoader", translation2=[0., 1., 0.], rotation2=[90., 180., 180.], showObjectScale=0.1)
    #skin.addObject('MechanicalObject', name="dofs", src="@meshLoader")
    skin.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="GeomAlgo")
    #skin.addObject('DiagonalMass', name="Mass", massDensity="1.0")
    ####skin.addObject('MeshMatrixMass', name="mass", totalMass= "0.00751", showAxisSizeFactor="0.005")
    skin.addObject('MeshMatrixMass', name="mass", totalMass= "0.00676", showAxisSizeFactor="0.005")
    ####skin.addObject('TetrahedralCorotationalFEMForceField', template="Vec3d", name="FEM", method="large", poissonRatio="0.35", youngModulus="130000", computeGlobalMatrix="0")
    skin.addObject('TetrahedralCorotationalFEMForceField', template="Vec3d", name="FEM", method="large", poissonRatio="0.45", youngModulus="12000", computeGlobalMatrix="0")
    #skin.addObject('FixedConstraint', name="FixedConstraint", indices="5 4 20 21 50")
    ####skin.addObject('FixedConstraint', indices=" 128 129 203 34 310 177 457 570 168 182")
    #skin.addObject('GenericConstraintCorrection', name="skin_ConstraintCorrection")

    ##---------------------Skin Collision Model---------------------------------------------

     # Skin COLLISION MODEL
    skinCollision = skin.addChild('skinCollision')
    skinCollision.addObject('MeshOBJLoader', name="loader", filename="D:\SOFA Learning\My Learning\mesh\Phantom Surface Meshes\ReducedMesh\Skin1000.obj")
    skinCollision.addObject('MeshTopology', src="@loader")
    skinCollision.addObject('MechanicalObject', name="skin_model", translation2=[0., 1., 0.], rotation2=[90., 180., 180.])
    skinCollision.addObject('TriangleCollisionModel')
    skinCollision.addObject('LineCollisionModel')
    skinCollision.addObject('PointCollisionModel')
    skinCollision.addObject('TetrahedronCollisionModel')
    skinCollision.addObject('BarycentricMapping', name="CollisionMapping", input="@../dofs", output="@skin_model")


    ##--------------------Skin Visual Model--------------------------------------------------

    # Skin VISUAL MODEL
    skinVisual = skin.addChild('skinVisual')
    
    skinVisual.addObject('MeshOBJLoader', name="loader", filename="D:\SOFA Learning\My Learning\mesh\Phantom Surface Meshes\Curved Phantom skin.obj")
    skinVisual.addObject('OglModel', name="VisualModel", translation=[0., 1., 0.], material="Default Diffuse 1 1 0 0 1 Ambient 1 0.2 0 0 1 Specular 0 1 0 0 1 Emissive 0 1 0 0 1 Shininess 0 45",
                          src="@loader",color=[0.921, 0.65, 0.101], rotation = [90, 180, 180])
    #skinVisual.addObject('OglModel', name="VisualModel", src="@loader", texturename="/home/darkn117/Desktop/texture-blood-vessel_36888-212.jpg")
    skinVisual.addObject('BarycentricMapping', name="VisualMapping", input="@../dofs", output="@VisualModel")
    

    ######################################################################################################################################################################
    ##
    ##                                          This part of the Scene defines Fat Layer of Tissue Phantom
    ##
    ######################################################################################################################################################################


    ##---------------------------Fat Mechanical Model---------------------------------

    fat = root.addChild('Fat')
    
    fat.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness="0.1", rayleighMass="0.1")
    fat.addObject('CGLinearSolver', name="linear_solver", iterations="25", tolerance="1e-09", threshold="1e-09")
    fat.addObject('MeshGmshLoader', name="meshLoader", filename="D:\SOFA Learning\My Learning\mesh\Phantom Surface Meshes\ReducedMesh\Fat1000.msh", triangulate="true", scale=0.1)
    #fat.addObject('MeshGmshLoader', name="meshLoader", filename="D:\SOFA Learning\My Learning\mesh\Phantom Surface Meshes\ReducedMesh\fat600.msh", triangulate="true")
    fat.addObject('TetrahedronSetTopologyContainer', name="topo", src="@meshLoader")
    fat.addObject('MechanicalObject', name="dofs", src="@meshLoader", translation2=[0., 0.5, 0.], rotation2=[90., 180., 180.], showObjectScale=0.1)
    #fat.addObject('MechanicalObject', name="dofs", src="@meshLoader")
    fat.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="GeomAlgo")
    #fat.addObject('DiagonalMass', name="Mass", massDensity="1.0")
    fat.addObject('MeshMatrixMass', name="mass", totalMass= "0.00691", showAxisSizeFactor="0.005")
    fat.addObject('TetrahedralCorotationalFEMForceField', template="Vec3d", name="FEM", method="large", poissonRatio="0.45", youngModulus="12000", computeGlobalMatrix="0")
    #fat.addObject('FixedConstraint', name="FixedConstraint", indices="5 4 20 21 50")
    ####fat.addObject('FixedConstraint', indices="0 3 4 67 68 71 72 75 77 139 140 142 144 147 148 211 212 215 216 219 221 283 284 286 288 289")
    #fat.addObject('GenericConstraintCorrection', name="fat_ConstraintCorrection")

    ##---------------------Fat Collision Model---------------------------------------------

     # FAT COLLISION MODEL
    fatCollision = fat.addChild('fatCollision')
    fatCollision.addObject('MeshOBJLoader', name="loader", filename="D:\SOFA Learning\My Learning\mesh\Phantom Surface Meshes\ReducedMesh\Fat1000.obj")
    fatCollision.addObject('MeshTopology', src="@loader")
    fatCollision.addObject('MechanicalObject', name="Fat_model", translation2=[0., 0.5, 0.], rotation2=[90., 180., 180.])
    fatCollision.addObject('TriangleCollisionModel')
    fatCollision.addObject('LineCollisionModel')
    fatCollision.addObject('PointCollisionModel')
    fatCollision.addObject('TetrahedronCollisionModel')
    fatCollision.addObject('BarycentricMapping', name="CollisionMapping", input="@../dofs", output="@Fat_model")


    ##--------------------Fat Visual Model--------------------------------------------------

    # FAT VISUAL MODEL
    fatVisual = fat.addChild('fatVisual')
    
    fatVisual.addObject('MeshOBJLoader', name="loader", filename="D:\SOFA Learning\My Learning\mesh\Phantom Surface Meshes\Curved Phantom Fat.obj")
    fatVisual.addObject('OglModel', name="VisualModel", translation=[0., 0.5, 0.], material="Default Diffuse 1 1 0 0 1 Ambient 1 0.2 0 0 1 Specular 0 1 0 0 1 Emissive 0 1 0 0 1 Shininess 0 45",
                          src="@loader",color=[0.941, 0.858, 0.313], rotation = [90, 180, 180])
    #fatVisual.addObject('OglModel', name="VisualModel", src="@loader", texturename="/home/darkn117/Desktop/texture-blood-vessel_36888-212.jpg")
    fatVisual.addObject('BarycentricMapping', name="VisualMapping", input="@../dofs", output="@VisualModel")


    ######################################################################################################################################################################
    ##
    ##                                          This part of the Scene defines Muscle Layer of Tissue Phantom
    ##
    ######################################################################################################################################################################


    ##---------------------------Muscle Mechanical Model---------------------------------

    muscle = root.addChild('Muscle')
    
    muscle.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness="0.1", rayleighMass="0.1")
    muscle.addObject('CGLinearSolver', name="linear_solver", iterations="50", tolerance="1e-09", threshold="1e-09")
    muscle.addObject('MeshGmshLoader', name="meshLoader", filename="mesh/Muscle600.msh", triangulate="true", scale=0.1)
    #muscle.addObject('MeshGmshLoader', name="meshLoader", filename="D:\SOFA Learning\My Learning\mesh\Phantom Surface Meshes\ReducedMesh\Muscle600.msh", triangulate="true")
    muscle.addObject('TetrahedronSetTopologyContainer', name="topo", src="@meshLoader")
    muscle.addObject('MechanicalObject', name="dofs", src="@meshLoader", showObjectScale=0.1,translation2=[0, 0, 0], rotation2=[90, 180, 180.])
    #muscle.addObject('MechanicalObject', name="dofs", src="@meshLoader")
    muscle.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="GeomAlgo")
    #muscle.addObject('DiagonalMass', name="Mass", massDensity="1.0")
    muscle.addObject('MeshMatrixMass', name="mass", totalMass= "0.04659", showAxisSizeFactor="0.005")
    muscle.addObject('TetrahedralCorotationalFEMForceField', template="Vec3d", name="FEM", method="large", poissonRatio="0.49", youngModulus="5000", computeGlobalMatrix="0")
    #muscle.addObject('FixedConstraint', name="FixedConstraint", indices="5 4 20 21 50")     246 399 574 
    muscle.addObject('FixedConstraint', indices="179 284 137 392 228 231 57 19 366 306 27 398  220 27 112 113 268 15 138 207 16 57 65 138 57 314 63 305 355 65 386 120 123 200 363 543")
    # muscle.addObject('FixedConstraint', indices="0 3 4 67 68 71 72 75 77 139 140 142 144 147 148 211 212 215 216 219 221 283 284 286 288 289")
    #muscle.addObject('GenericConstraintCorrection', name="Muscle_ConstraintCorrection")

    ##---------------------Muscle Collision Model---------------------------------------------

     # Pumonary Vein COLLISION MODEL
    muscleCollision = muscle.addChild('muscleCollision')
    muscleCollision.addObject('MeshOBJLoader', name="loader", filename="mesh/Muscle600.obj")
    muscleCollision.addObject('MeshTopology', src="@loader")
    muscleCollision.addObject('MechanicalObject', name="Muscle_model",translation2=[0, 0, 0], rotation2=[90, 180, 180.])
    muscleCollision.addObject('TriangleCollisionModel')
    muscleCollision.addObject('LineCollisionModel')
    muscleCollision.addObject('PointCollisionModel')
    muscleCollision.addObject('BarycentricMapping', name="CollisionMapping", input="@../dofs", output="@Muscle_model")


    ##--------------------Muscle Visual Model--------------------------------------------------

    # Pumonary Vein VISUAL MODEL
    muscleVisual = muscle.addChild('muscleVisual')
    
    muscleVisual.addObject('MeshOBJLoader', name="loader", filename="mesh/Muscle Body.obj",translation=[0, 0, 0], rotation=[90, 180, 180.])
    muscleVisual.addObject('OglModel', name="VisualModel", material="Default Diffuse 1 1 0 0 1 Ambient 1 0.2 0 0 1 Specular 0 1 0 0 1 Emissive 0 1 0 0 1 Shininess 0 45",
                          src="@loader",color=[1, 0.05, 0.3])
    #muscleVisual.addObject('OglModel', name="VisualModel", src="@loader", texturename="/home/darkn117/Desktop/texture-blood-vessel_36888-212.jpg")
    muscleVisual.addObject('BarycentricMapping', name="VisualMapping", input="@../dofs", output="@VisualModel")

[bakpaul]

Hello !

Many parameters influence this behavior. We can start with the ones with direct influence :

• The LocalMinDistance has two parameters : the alarmDistance and the contactDistance.
  • The contactDistance is a kind of skinning width, the contact doesn't arise at the exact point of the model, but a bit in front. This helps to avoid interpenetration due to numerical noise.
  • The alarmDistance is the distance between objects at which we start to add constraints, even if there is no contact. This parameter is very important because it avoid to miss contact when object are subject to too much forces (might be the case here). Here I find it a bit small with respect to the contactDistance.
  • The collision model also have a distance parameter which is proximity that will add more distance between the contact point and the actual point used for the solving. this plays a role of "hysteresis": by solving the contact while enforcing overshooting of the response (more space than necessary) you make sure that the modle will still be on the right side of the other object's surface at the end of the solving

You are using Penality-based constraints which are very weak and doesn't grantee to have no interpenetration at the end of the time step. So for contact with a lot of stress like this this is really not recommended. You should use Lagrangian-based constraints instead.

[bakpaul]

And by the way you have two animationLoop and A constraint solver : only one animationLoop is needed and because you are not using Lagrangian-based constraint, you don't need the GenericConstraintSolver