Surface Pressure Constraint causes unrealistic sliding

[MikeP]

I have been evolving the geometries and layouts of simple soft robots which are made up of soft parts, rigid parts and actuating parts. I have successfully managed to do this using Sofa and SofaPython3 where I use and control the pressure value in a Surface Pressure Forcefield component to actuate the actuating parts. Unfortunately, I was not getting the desired result. As the pressures used would be sufficient when the actuating part was large but would barely do anything when the actuating part was small. I then attempted to use the Surface pressure constraint from the Soft Robots plugin instead which seemed to solve the issue of smaller volumes not actuating as much as larger volumes, and also allowed me to specifically actuate the volume growth iteself. Unfortunately, it created a new issue – when the soft robot’s actuator inflates the soft robot begins to slide along the surface, despite friction contact being implemented. It appears that when the actuator inflates, a reaction force is created causing the robot to slide.

1. Please could someone explain how both the Surface Pressure forcefield and the surface pressure constraint works and why they don’t seem to produce similar results (for example a pressure value of 150 for the forcefield component seems to equate to a pressure value of 1.5 for the contraint component).

2. How can I prevent the unrealistic sliding caused by using the surface constraint component?

The code illustrating an example is provided below:

import sys, os
import math

sys.path.append(os.path.abspath("./bindings/Sofa/package"))
sys.path.append(os.path.abspath("./bindings/SofaRuntime/package"))
sys.path.append(os.path.abspath("./bindings/SofaTypes/package"))

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

class FingerController(Sofa.Core.Controller):
    """ This is a custom controller to perform actions when events are triggered """

    def __init__(self, *args, **kwargs):
        # These are needed (and the normal way to override from a python class)
        Sofa.Core.Controller.__init__(self, *args, **kwargs)
        self.node = kwargs["node"]
        self.rootNode = kwargs["rootNode"]
        self.act1_size = kwargs["act1_size"]
        return

    def onAnimateBeginEvent(self, event):
        if self.rootNode.time.value > 1.5:
            self.node.Surface.surfaceConstraint.value.value = [0.6*math.sin(self.rootNode.time.value-1.5-math.pi/2)+0.6]

class FingerController2(Sofa.Core.Controller):
    """ This is a custom controller to perform actions when events are triggered """

    def __init__(self, *args, **kwargs):
        # These are needed (and the normal way to override from a python class)
        Sofa.Core.Controller.__init__(self, *args, **kwargs)
        self.node = kwargs["node"]
        self.rootNode = kwargs["rootNode"]
        self.act2_size = kwargs["act2_size"]
        return

    def onAnimateBeginEvent(self, event):
        if self.rootNode.time.value > 1.5+math.pi:
            self.node.Surface.surfaceConstraint.value.value = [-0.6*math.sin(self.rootNode.time.value-1.5-math.pi/2)+0.6]

def createScene(root):
    # Set time step and gravity of the Scene
    root.gravity=[0, -9.81, 0]
    root.dt=0.01
    root.addObject('RequiredPlugin', pluginName="SoftRobots SofaGeneralObjectInteraction SofaConstraint SofaGeneralSimpleFem SofaOpenglVisual SofaGeneralLoader SofaImplicitOdeSolver SofaLoader SofaMiscCollision SofaSimpleFem")
    # Scene must now include a VisualLoop
    root.addObject('DefaultVisualManagerLoop')
    # Scene must now include a AnimationLoop
    root.addObject('FreeMotionAnimationLoop')
    root.addObject('GenericConstraintSolver', tolerance=1e-3, maxIt=1000)

    root.addObject('DefaultPipeline', verbose=0)
    root.addObject('BruteForceDetection', name="N2")
    root.addObject('DefaultContactManager', name="Response", response="FrictionContact", responseParams="mu=1")
    root.addObject('NewProximityIntersection', name="Proximity", alarmDistance=0.8, contactDistance=0.5)
    root.addObject('DefaultCollisionGroupManager', name="Group")

    floor = root.addChild('Floor')
    floor.addObject('MeshGmshLoader', name="floorLoader", filename="/home/mike/SOFA_v20.12.01_Linux/share/sofa/mesh/test_meshes/floor.msh")
    floor.addObject('QuadSetTopologyContainer', name="quadContainer", src="@floorLoader")
    floor.addObject('TriangleSetTopologyContainer', name='triangleContainer')
    floor.addObject('TriangleSetTopologyModifier', name='triangleModifier')
    floor.addObject('Quad2TriangleTopologicalMapping', input="@quadContainer", output="@triangleContainer")
    floor.addObject('MechanicalObject', name="FloorDOF", src="@triangleContainer", scale=50, translation=[-50.0, 0.0, -50.0], rotation=[0.0, 0.0, 0.0])
    floor.addObject('TriangleCollisionModel', src="@triangleContainer", simulated=0, moving=0, bothSide=1, group=1)
    floor.addObject('LineCollisionModel', src="@triangleContainer", simulated=0, moving=0, group=1)
    floor.addObject('PointCollisionModel', src="@triangleContainer", simulated=0, moving=0, group=1)
    floor.addObject('OglModel', name="FloorVis", color="white")

    soft = root.addChild('soft')
    soft.addObject('EulerImplicitSolver', rayleighStiffness=0)
    soft.addObject('CGLinearSolver', iterations=200, tolerance=1e-06, threshold=1e-06)
    soft.addObject('MeshGmshLoader', name="meshLoader", filename="/home/mike/meshes/FourPart/3*3*3/5softpart.msh")
    soft.addObject('HexahedronSetTopologyContainer', name="topo", src="@meshLoader")
    soft.addObject('MechanicalObject', name="softDOF", template="Vec3d", translation=[0.0, 1.0, 0.0])
    soft.addObject('UniformMass', vertexMass=0.3)
    soft.addObject('PrecomputedConstraintCorrection')
    soft.addObject('HexahedralFEMForceField', name="FEM", youngModulus=500, poissonRatio=0.3, method="large")
    surface = soft.addChild('Surface')
    surface.addObject('QuadSetTopologyContainer', name="QuadContainer")
    surface.addObject('QuadSetTopologyModifier', name="QuadModifier")
    surface.addObject('Hexa2QuadTopologicalMapping', input="@../topo", output="@QuadContainer")
    collision = surface.addChild('Collsion')
    collision.addObject('TriangleSetTopologyContainer', name="TriContainer")
    collision.addObject('TriangleSetTopologyModifier', name="TriModifier")
    collision.addObject('Quad2TriangleTopologicalMapping', input="@../QuadContainer", output="@TriContainer")
    collision.addObject('TriangleCollisionModel', name="CollisionMO", template="Vec3d", group=2)
    collision.addObject('LineCollisionModel', group=2)
    collision.addObject('PointCollisionModel', group=2)
    collision.addObject('OglModel', name="softVis", color="red")
    collision.addObject('BarycentricMapping', input="@../../softDOF", output="@softVis")

    actuator1 = root.addChild('actuator1')
    actuator1.addObject('EulerImplicitSolver', rayleighStiffness=0)
    actuator1.addObject('CGLinearSolver', iterations=200, tolerance=1e-06, threshold=1e-06)
    actuator1.addObject('MeshGmshLoader', name="meshLoader", filename="/home/mike/meshes/FourPart/3*3*3/5actuator1part.msh")
    actuator1.addObject('HexahedronSetTopologyContainer', name="topo", src="@meshLoader")
    actuator1.addObject('MechanicalObject', name="actuator1DOF", translation=[0.0, 1.0, 0.0])
    actuator1.addObject('UniformMass', vertexMass=0.3)
    actuator1.addObject('PrecomputedConstraintCorrection')
    surface = actuator1.addChild('Surface')
    surface.addObject('QuadSetTopologyContainer', name="QuadContainer")
    surface.addObject('QuadSetTopologyModifier', name="QuadModifier")
    surface.addObject('Hexa2QuadTopologicalMapping', input="@../topo", output="@QuadContainer")
    surface.addObject('MeshSpringForceField', name="Springs", stiffness=200, damping=1)
    # surface.addObject('SurfacePressureForceField', name="Pressure1", pressure=0, pulseMode=0)
    surface.addObject('SurfacePressureConstraint', name='surfaceConstraint', template='Vec3d', value=0, valueType='pressure')
    collision = surface.addChild('Collision')
    collision.addObject('TriangleSetTopologyContainer', name="TriContainer")
    collision.addObject('TriangleSetTopologyModifier', name="TriModifier")
    collision.addObject('Quad2TriangleTopologicalMapping', input="@../QuadContainer", output="@TriContainer")
    collision.addObject('TriangleCollisionModel', name="CollisionMO", template="Vec3d", group=2)
    collision.addObject('LineCollisionModel', group=2)
    collision.addObject('PointCollisionModel', group=2)
    collision.addObject('OglModel', name="actuator1Vis", color="green")
    collision.addObject('BarycentricMapping', input="@../../actuator1DOF", output="@actuator1Vis")

    actuator2 = root.addChild('actuator2')
    actuator2.addObject('EulerImplicitSolver', rayleighStiffness=0)
    actuator2.addObject('CGLinearSolver', iterations=200, tolerance=1e-06, threshold=1e-06)
    actuator2.addObject('MeshGmshLoader', name="meshLoader", filename="/home/mike/meshes/FourPart/3*3*3/5actuator2part.msh")
    actuator2.addObject('HexahedronSetTopologyContainer', name="topo", src="@meshLoader")
    actuator2.addObject('MechanicalObject', name="actuator2DOF", translation=[0.0, 1.0, 0.0])
    actuator2.addObject('UniformMass', vertexMass=0.3)
    actuator2.addObject('PrecomputedConstraintCorrection')
    surface = actuator2.addChild('Surface')
    surface.addObject('QuadSetTopologyContainer', name="QuadContainer")
    surface.addObject('QuadSetTopologyModifier', name="QuadModifier")
    surface.addObject('Hexa2QuadTopologicalMapping', input="@../topo", output="@QuadContainer")
    surface.addObject('MeshSpringForceField', name="Springs", stiffness=200, damping=1)
    # surface.addObject('SurfacePressureForceField', name="Pressure2", pressure=2, pulseMode=0)
    surface.addObject('SurfacePressureConstraint', name='surfaceConstraint', template='Vec3d', value=0, valueType='pressure')
    collision = surface.addChild('Collision')
    collision.addObject('TriangleSetTopologyContainer', name="TriContainer")
    collision.addObject('TriangleSetTopologyModifier', name="TriModifier")
    collision.addObject('Quad2TriangleTopologicalMapping', input="@../QuadContainer", output="@TriContainer")
    collision.addObject('TriangleCollisionModel', name="CollisionMO", template="Vec3d", group=2)
    collision.addObject('LineCollisionModel', group=2)
    collision.addObject('PointCollisionModel', group=2)
    collision.addObject('OglModel', name="actuator2Vis", color="yellow")
    collision.addObject('BarycentricMapping', input="@../../actuator2DOF", output="@actuator2Vis")

    root.addObject('AttachConstraint', object1='@soft', object2='@actuator1', twoWay=True, indices1=[0, 1, 2, 3, 4, 5], indices2=[6, 7, 8, 9, 10, 11])
    root.addObject('AttachConstraint', object1='@soft', object2='@actuator2', twoWay=True, indices1=[2, 5, 8, 11], indices2=[4, 6, 8, 10])
    root.addObject('AttachConstraint', object1='@actuator1', object2='@actuator2', twoWay=True, indices1=[2, 5, 8, 11], indices2=[0, 2, 4, 6])

    root.addObject(FingerController(name="ActuationController", rootNode=root, node=actuator1, act1_size=2))
    root.addObject(FingerController2(name="ActuationController2", rootNode=root, node=actuator2, act2_size=2))

def main():
    import SofaRuntime
    import Sofa.Gui
    # Make sure to load all SOFA libraries
    SofaRuntime.importPlugin("SofaComponentAll")
    #Create the root node
    root = Sofa.Core.Node("root")
    # Call the below 'createScene' function to create the scene graph
    createScene(root)
    Sofa.Simulation.init(root)

    if not USE_GUI:
        for iteration in range(1000):
            Sofa.Simulation.animate(root, root.dt.value)
    else:
        # Find out the supported GUIs
        print ("Supported GUIs are: " + Sofa.Gui.GUIManager.ListSupportedGUI(","))
        # 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("GUI was closed")

    print("Simulation is done.")
    softTotal = [0.0, 0.0, 0.0]
    softPos = root.soft.softDOF.position.value
    for vertex in softPos:
        softTotal += vertex
    actuator1Pos = root.actuator1.actuator1DOF.position.value
    for vertex in actuator1Pos:
        softTotal += vertex
    softAverage = softTotal/(len(softPos)+len(actuator1Pos))
    print(softAverage)
    fitness = ((softAverage[0])**2 + (softAverage[1])**2 + (softAverage[2])**2)**0.5
    print(fitness)

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

The mesh files can be found here

[Hugo]

Hey @quantanovo

I am so sorry for the delay of my reply! I am having hard time to reply all questions, sorry for keeping you waiting.

1. The SurfacePressureForceField is an external force applied to your object, applying a specific force on the surface, which is spread over the surface points depending on the surface area.
On the other hand, the SurfacePressureConstraint is a BaseConstraint which have to be solved using the Lagrange Multiplier method. This component will therefore constraint the object (find the according forces) to deform so that the pressure inside the volume is ensured.

The difference of forces remind me of this topic. The forces resulting from the constraint resolution (forces ensuring that the pressure inside the volume is conserved) with SurfacePressureConstraint are affected be the time steps. It looks like there is a problem here. Could you try to change the time step and see if it affects the ratio of force (1.5/150, if you change the time step with dt=0.001 should become 0.15/150)

could you keep me posted ?

2. Regarding the sliding, have you tried increasing again the mu friction coefficient to see it this creates the desired effect?

Best wishes,

Hugo

[Hugo]

Hi @quantanovo

Even if it took me some time, I carried on the investigation. It appears that the SurfacePressureConstraint is timestep dependent and the pressure set as input has to actually be related to pressure/dt as explained in the associated GitHub issue.

You may want to follow closely the ongoing discussion there.
Best wishes,

Hugo

[Author]

Posts