Interaction between two soft objects

[CelSag]

Hello team,
I'm trying to simulate the behavior of a soft object when it's being pressed with a soft finger actuator.

I haven't found any related topic/ resources with such example, so I just tried to do collisions like it's done between soft/rigid objects, with

    rootNode.addObject('FreeMotionAnimationLoop')
    rootNode.addObject('DefaultPipeline')
    rootNode.addObject('BruteForceBroadPhase')
    rootNode.addObject('BVHNarrowPhase')
    rootNode.addObject('LocalMinDistance', name='Proximity', alarmDistance=0.2, contactDistance=0.1)
    rootNode.addObject('DefaultCollisionGroupManager',  name="Group")

    rootNode.addObject('DefaultContactManager', response="FrictionContactConstraint")
    rootNode.addObject('GenericConstraintSolver', tolerance=1e-12, maxIterations=10000)

and both of my objects defined with a TetrahedronFEMForceField (and everything that goes with it), EulerImplicitSolver, SparseLDLSolver, FixedConstraint, LinearSolverConstraintCorrection and a collision model.

The finger is also controlled with a SurfacePressureConstraint.
The simulation works normally at the beginning, but after the collision, the model of the finger explodes for some reason.

I would really appreciate it if you could guide me toward a related example or advise on how to solve this problem. Here's the full code in case you need it:

import Sofa
import os
import Sofa.Core
import Sofa.Simulation
import SofaRuntime
import numpy as np

radius = 50
translateFinger1 = [2,  -6.5 +radius, -7.5]

def createScene_collision(rootNode):
    PATH = os.path.abspath(os.getcwd())

    rootNode.addObject('RequiredPlugin',
                   pluginName='SofaComponentAll SofaOpenglVisual SoftRobots SofaPython3 SofaLoader SofaSimpleFem SofaEngine SofaDeformable SofaImplicitOdeSolver SofaConstraint SofaSparseSolver')

    rootNode.addObject('VisualStyle', displayFlags=' showBehaviorModels hideMappings showForceFields')
    rootNode.findData('gravity').value = [0,0, 0] #[9806, 0, 0]
    rootNode.addObject('OglSceneFrame', style='Arrows', alignment='TopRight')

    rootNode.addObject('FreeMotionAnimationLoop')
    rootNode.addObject('DefaultPipeline')
    rootNode.addObject('BruteForceBroadPhase')
    rootNode.addObject('BVHNarrowPhase')
    rootNode.addObject('LocalMinDistance', name='Proximity', alarmDistance=0.2, contactDistance=0.1)
    rootNode.addObject('DefaultCollisionGroupManager',  name="Group")

    rootNode.addObject('DefaultContactManager', response="FrictionContactConstraint") 
    rootNode.addObject('GenericConstraintSolver', tolerance=1e-12, maxIterations=10000)
   
    
    # # ### FINGER
    # # ##########################################
    ###getting indices to block
    finger = rootNode.addChild('finger_init')
    finger.addObject('MeshGmshLoader', name="loader", filename= PATH + '/models/simple_pneunet/for_sofa/basic_finger_tetra.msh', translation = translateFinger1) 
    finger.addObject('MeshTopology', src='@loader', name='container')
    finger.addObject('MechanicalObject', name='tetras', template='Vec3', showIndices=False, showIndicesScale=4e-5)
    
    finger.addObject('BoxROI', name='boxROI', box=[-1, radius -10, -7.5, 1 ,radius + 10, 7.5], drawBoxes = True)
    finger.init()

    idx_blocked = rootNode.finger_init.boxROI.indices
    finger_init = rootNode.finger_init
    rootNode.removeChild(finger_init)

    ### -------------------------------------
    finger = rootNode.addChild('finger')
    finger.addObject('EulerImplicitSolver', name='odesolver', rayleighStiffness=0.1, rayleighMass=0.1)
    finger.addObject('SparseLDLSolver', name='preconditioner', template="CompressedRowSparseMatrixd")

    loader = finger.addObject('MeshGmshLoader', name="loader", filename= PATH + '/models/simple_pneunet/for_sofa/basic_finger_tetra.msh') 
    finger.addObject('MeshTopology', src='@loader', name='container')
    MO = finger.addObject('MechanicalObject', name='tetras', template='Vec3', showIndices=False, showIndicesScale=4e-5)
    finger.addObject('DiagonalMass', name="Mass", massDensity="8.8e-10") 
    finger.addObject('TetrahedronFEMForceField', template='Vec3d', name='FEM', method='large', poissonRatio=0.4,  youngModulus=0.6e3)
    finger.init()
    finger.addObject('FixedConstraint', indices=idx_blocked)
    finger.addObject('LinearSolverConstraintCorrection', solverName='preconditioner')

    # Constraint                             
    cavity = finger.addChild('cavity')
    cavity.addObject('MeshGmshLoader', name='loader_cavity', filename='models/simple_pneunet/for_sofa/Cavity_triangle.msh')
    cavity.addObject('MeshTopology', src='@loader_cavity', name='topo')
    cavity.addObject('MechanicalObject', name='cavity')
    cavity.addObject('SurfacePressureConstraint', name='SurfacePressureConstraint', template='Vec3', value=0.0001, triangles='@topo.triangles', valueType='pressure')
    cavity.addObject('BarycentricMapping', name='mapping',  mapForces=False, mapMasses=False)

    # Collision                           
    collisionFinger = finger.addChild('collisionFinger')
    collisionFinger.addObject('MeshGmshLoader', name='loader_colli', filename='models/simple_pneunet/for_sofa/basic_finger_triangle.msh')
    collisionFinger.addObject('MeshTopology', src='@loader_colli', name='topo')
    collisionFinger.addObject('MechanicalObject', name='collisMech')
    collisionFinger.addObject('TriangleCollisionModel', selfCollision=False)
    collisionFinger.addObject('LineCollisionModel', selfCollision=False)
    collisionFinger.addObject('PointCollisionModel', selfCollision=False)
    collisionFinger.addObject('BarycentricMapping')

    # Controller                             
    finger.addObject(PressureController(name="PressureController", node = rootNode))

    # # ### OBJECT
    # # ##########################################
    fruit = rootNode.addChild('fruit')
    fruit.addObject('EulerImplicitSolver', name='odesolver')
    fruit.addObject('SparseLDLSolver', name='linearSolver', template="CompressedRowSparseMatrixd")

    fruit.addObject('MeshGmshLoader', name="loader", filename= PATH + '/fruit_tetra.msh', translation=[70, -80, 0])
    fruit.addObject('MeshTopology', src='@loader', name='container')
    fruit.addObject('MechanicalObject', name='tetras', template='Vec3', showIndices=False, showIndicesScale=4e-5)
    fruit.addObject('TetrahedronFEMForceField', template='Vec3d', name='FEM', method='large', poissonRatio=0.4,  youngModulus=0.8e3, computeVonMisesStress="2", showVonMisesStressPerElement="true")

    fruit.addObject('DiagonalMass', name="Mass", massDensity="1.15e-6")
    fruit.addObject('BoxROI', name='boxROI', box=[50, -110, -30, 100 ,-50, 30], drawBoxes = True)
    fruit.init()
    fruit.addObject('FixedConstraint', indices ='@boxROI.indices' )
    fruit.addObject('LinearSolverConstraintCorrection', solverName='linearSolver')

    #collision
    fruitCollis = fruit.addChild('fruitCollis')
    fruitCollis.addObject('MeshGmshLoader', name='loader', filename=PATH + '/fruit_triangle.msh', translation=[70, -80, 0])
    fruitCollis.addObject('MeshTopology', src='@loader')
    fruitCollis.addObject('MechanicalObject')
    fruitCollis.addObject('TriangleCollisionModel')
    fruitCollis.addObject('LineCollisionModel')
    fruitCollis.addObject('PointCollisionModel')
    fruitCollis.addObject('BarycentricMapping')

    return rootNode

class PressureController(Sofa.Core.Controller):
    def __init__(self, *a, **kw):
        Sofa.Core.Controller.__init__(self, *a, **kw)
        self.node = kw["node"]
        self.constraints = []
        self.dofs = []
        self.dofs.append(self.node.getChild('finger').tetras)
        self.constraints.append(self.node.getChild('finger').cavity.SurfacePressureConstraint)

    def onKeypressedEvent(self,e):
        increment = 0.01
        key = e['key']
        if ord(key) == 4: ##enter key 
            i = 0
            pressureValue = self.constraints[i].value.value[0] + increment
            if pressureValue > 1.5:
                pressureValue = 1.5
            print(pressureValue)
            self.constraints[i].value = [pressureValue]

        if ord(key) == 3: ##delete key 
            i = 0
            pressureValue = self.constraints[i].value.value[0] - increment
            if pressureValue < 0:
                pressureValue = 0
            print(pressureValue)
            self.constraints[i].value = [pressureValue]

def main():
    import Sofa.Gui

    # Make sure to load all SOFA libraries and plugins
    SofaRuntime.importPlugin("SofaBaseMechanics")
    SofaRuntime.importPlugin("SofaComponentAll")
    SofaRuntime.importPlugin("Sofa.Component.IO.Mesh")

    #Create the root node
    root = Sofa.Core.Node("root")
    root = createScene_collision(root)
    Sofa.Simulation.init(root)

    # Launch the GUI (qt or qglviewer)
    Sofa.Gui.GUIManager.Init("myscene", "qglviewer")
    Sofa.Gui.GUIManager.createGUI(root, __file__)
    Sofa.Gui.GUIManager.SetDimension(1080, 1080)
    # Initialization of the scene will be done here
    Sofa.Gui.GUIManager.MainLoop(root)
    Sofa.Gui.GUIManager.closeGUI()
    print("Simulation is done.")

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

I can provide the meshes is needed.
Any help is welcome!!

[hugtalbot]

Hi @CelSag

First of all welcome and thank you creating a full and clear post!

Is the simulation crashing / diverging exactly at the moment of the contact ?
How do you define your time step dt exactly (defined in the root node)? You might have a too-large time step.

Another thing to investigate would be whether the GenericConstraintSolver does converge. You could activate the printLog data to read more about this.

Best wishes.

[CelSag]

Hi Hugo,
Thank you for your reply!

To answer your questions;

• the simulation is crashing a few iterations after the contact
• i'm working with dt = 0.01, i tried with a smaller one but the same issue occurs
• following your advice, i activated the GenericConstraintSolver printlog; it seems that it does converge, but here are the results:

At contact, which seems to work fine:

Then after a few iterations, when it starts acting up:

It seems that the number of constraints alternates between 13, 10, 7 and 4 from that point on. Is there a reason why the number of constraints suddenly goes up and down? Could it be why it crashes?

Thanks again for your help!

[hugtalbot]

Hi @CelSag

Thanks for the feedback.
The number of constraints may indeed change since the number of contact points may itself change depending on the relative motions of both objects. Nothing strange here.

i'm working with dt = 0.01, i tried with a smaller one but the same issue occurs

When you say smaller, how much smaller?

Did you check the coherency of your units?

Best wishes

[github-actions[bot]]

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

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

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