Assistance with mapping an internal stiff spring forcefield to a soft actuator

[tolstoys19]

Hello.

I have been trying to model a fiber reinforced soft actuator for my phD project. I have succesfully added the internal air pressure caivty using the surfacepressure constraint of the SoftRobots plugin. To model the soft actuators themselves I am basing myself on the work done on the following paper:

https://ieeexplore.ieee.org/abstract/document/9981637
10.1109/IROS47612.2022.9981637

Where the fibers are modeled as parallel stiff springs and then a barycentric mapping is used to map the stiffness of the fibers to the actuator. For now in my case my fiber is just a single line that goes through the middle of the bottom layer of my actuator.

I am able to add the fiber as a stiff spring and map it to the soft actuator body and the simulation can run using the CGlinear solver. However when I add the soft actuator cavity I need to change to the GenericConstraintSolver and add a FreeAnimationLoop. When i do this my simulation crashes and the only I wan get it to run is to stop mapping the forces of the fiber node which then makes the reinforcement useless.

Do you have anyway to solve or workaround this problem?

I have attached my files here.

Thank you in advance for your help
Question.zip

[EulalieCoevoet]

Hello,

Currently the forces under a mapped node are not correctly mapped when using a SparseLDLSolver, which you need, to use the combination of GenericConstraintSolver/FreeAnimationLoop, if I'm correct.
There's a WIP developpement for this problem (sorry I didn't find any references) and a temporary solution ; you can use the component MechanicalMatrixMapper.

Here's a working scene:

def createScene(rootNode):
    rootNode.addObject('RequiredPlugin', pluginName=['SoftRobots'])
    rootNode.addObject('VisualStyle', displayFlags='showBehavior')
    rootNode.gravity.value = [9810, 0, 0]
    rootNode.addObject('FreeMotionAnimationLoop')
    rootNode.addObject('GenericConstraintSolver', tolerance=1e-12, maxIterations=10000)
    rootNode.addObject('EulerImplicitSolver', name='odesolver', rayleighStiffness=0.1, rayleighMass=0.1)

    finger = rootNode.addChild('finger')
    finger.addObject('SparseLDLSolver')
    finger.addObject('MeshVTKLoader', name='loader', filename='data/mesh/actwithconnector.vtk')
    finger.addObject('MeshTopology', src='@loader', name='container')
    finger.addObject('MechanicalObject', name='tetras', template='Vec3', showObject=True, showObjectScale=1)
    finger.addObject('TetrahedronFEMForceField', template='Vec3', name='FEM', method='large', poissonRatio=0.36,
                     youngModulus=500)
    finger.addObject('UniformMass', totalMass=3)
    finger.addObject('BoxROI', name='boxROISubTopo', box=[28, -10, -8, 31, 10, 8], strict=False, drawBoxes=True)
    finger.addObject('BoxROI', name='boxROISubTopo2', box=[0, -10, -4, 28, 10, -2], strict=False, drawBoxes=True)
    finger.addObject('BoxROI', name='boxROI', box=[-1.5, -10, -20, 1.5, 10, 20], drawBoxes=True)
    finger.addObject('RestShapeSpringsForceField', points='@boxROI.indices', stiffness=1e12, angularStiffness=1e12)
    finger.addObject('LinearSolverConstraintCorrection')

    fiber = finger.addChild("fiber")
    nbDOFs = 10
    fiber.addObject("MechanicalObject", template="Vec3", name="DOF",
                    position=[[i*3, 0, -4] for i in range(nbDOFs)],
                    showObject=True, showObjectScale=1)
    fiber.addObject('MeshTopology', name='lines', lines=[[i, i + 1] for i in range(nbDOFs-1)])
    fiber.addObject('UniformMass', totalMass=0.8)
    fiber.addObject("FixedConstraint", name="FixedConstraint", indices=[0])

    Ks = 1e6
    Kd = 5
    L = 3
    fiber.addObject("StiffSpringForceField", template="Vec3d", name="springs", showArrowSize=0.5, drawMode=1,
                    spring=[[i, i + 1, Ks, Kd, L] for i in range(nbDOFs-1)])
    fiber.addObject('BarycentricMapping', name='mapping')

    # Temporary solution
    finger.addObject('MechanicalMatrixMapper', template="Vec3,Vec3",
                         nodeToParse=fiber.linkpath,  # where to find the forces to map
                         object1=finger.tetras.linkpath,  # where to map the forces
                         object2=finger.tetras.linkpath)  # in case of multi-mapping, here you can give the second parent
                        

    cavity = finger.addChild('cavity')
    cavity.addObject('MeshVTKLoader', name='cavityLoader', filename='data/mesh/CHAMBER.vtk')
    cavity.addObject('MeshTopology', src='@cavityLoader', name='cavityMesh')
    cavity.addObject('MechanicalObject', name='cavity')
    cavity.addObject('SurfacePressureConstraint', name='SurfacePressureConstraint', template='Vec3', value=0.5,
                     triangles='@cavityMesh.triangles', valueType='pressure')
    cavity.addObject('BarycentricMapping', name='mapping')

[tolstoys19]

Thank You Very much for your answer.　i was able to run the simulation but I needed to change some lines from your code.
nodeToParse=fiber.linkpath -> fiber.getLinkPath()
object1=finger.tetras.linkpath -> finger.tetras.getLinkPath()
object2=finger.tetras.linkpath-> finger.tetras.getLinkPath()

I also uploaded my new files here in case someone needs them as reference later on.

again, thank you very much for your help
Question (1).zip

[EulalieCoevoet]

I'm glad you managed to make it work. Thanks for your feedback and for sharing your code.

Just so you know, linkpath is a new (and prefered) feature of SofaPython3. "It's better than getLinkPath() in terms of performance and it allows to make links between objects that are not yet in the same graph." I'm quoting a SofaPython3 dev here ^^
I think it should be in the last release of SOFA (v22.12), I've been using it for a while now.

[hugtalbot]

Hi @tolstoys19

EDIT: I see that @EulalieCoevoet already answered before me! Anyway I post my short reply

Thank you for your question!
First of all, I will poke @sescaida since he's one of the authors of the paper which you refer to.

Second, it seems to me that there is some confusion here "when I add the soft actuator cavity I need to change to the GenericConstraintSolver and add a FreeAnimationLoop". AnimationLoops and ConstraintSolvers are two different components:

• the AnimationLoop is in charge of defining the step to be computed within one time step (and their order)
• the ConstraintSolver is a solver in charge of solving the constraint problem using the Lagrange multiplier method. To use such constraint resolution (as Eulalie wrote), a FreeMotionAnimationLoop is however the only AnimationLoop supporting constraint solving using the Lagrange multiplier.

I hope this helps