[SOLVED] BilateralInteractionConstraint reduces gravity by 100x

[Joshua]

Hello,
I’m attempting to model a soft body attached to a rigid body. As suggested in several posts on the forum, I rigidly map the points of a surface mesh to the 6dof rigid body, and use a BilateralInteractionConstraint to attach the soft body to the rigid surface. This works as a method of attachment, but the constraint seems to dilute gravitational (and possibly other) forces on the soft body by a factor of 100x as they are experienced by the rigid body. In other words, an upward force of only 1% of the soft body’s weight applied to the rigid body is enough to compensate for the additional weight of the entire soft body. When the upward force is applied directly to the soft body, a force equal to the total weight of that body is needed (as expected). Can anyone explain this discrepancy?

import Sofa
from stlib.scene import Scene
from stlib.visuals import VisualModel

import numpy as np

def createScene(rootNode):
    scene = Scene(rootNode, plugins=["SoftRobots",
                                     "SofaSparseSolver",
                                     "SofaMiscCollision",
                                     "SofaPython"])
    # scene.VisualStyle.displayFlags = 'showForceFields showBehaviorModels showInteractionForceFields'
    rootNode.createObject("BackgroundSetting", color=[1., 1., 1.])

    rootNode.dt = 0.01
    rootNode.gravity = [0., 0., -981.]

    rootNode.createObject("FreeMotionAnimationLoop")
    rootNode.createObject("GenericConstraintSolver", maxIterations=1e3, tolerance=1e-5)

    # Rigid object

    rigid_node = rootNode.createChild("rigid")

    totalMass = .01
    inertiaMatrix = [1., 0., 0., 0., 1., 0., 0., 0., 1.]
    volume = 1.0
    uniformScale = 20

    rigid_node.createObject('EulerImplicitSolver', name='odesolver')
    rigid_node.createObject('CGLinearSolver', name='Solver')

    rigid_node.createObject("MechanicalObject", template="Rigid3", name="dofs", showObjectScale=uniformScale)
    mass = rigid_node.createObject('UniformMass', name="mass", vertexMass=[totalMass, volume, inertiaMatrix[:]])
    mass.totalMass = totalMass

    rigid_node.createObject('UncoupledConstraintCorrection', compliance=1e-8)

    points = rigid_node.createChild('points')
    loader = points.createObject('MeshObjLoader', name='loader', filename='mesh/cube.obj', scale=uniformScale)
    points.createObject('Mesh', src=loader, name="MeshPoints")
    points.createObject('MechanicalObject', name='dofs', template='Vec3d', showObject=True, showObjectScale=1,
                        src=loader)
    # geomAlgo = points.createObject('TriangleSetGeometryAlgorithms', template="Vec3d", name='GeomAlgo')
    # mass = points.createObject('DiagonalMass', massDensity=.00001, name='mass')

    points.createObject('RigidMapping', name="mapping", input=rigid_node, output=points)
    points.createObject('UncoupledConstraintCorrection', compliance=0.002)

    visual = VisualModel(rigid_node, 'mesh/cube.obj', scale=[uniformScale]*3,
                         color=[0.75, 0.75, 0.25, 1.])
    visual.writeZTransparent = True
    visual.createObject('RigidMapping', name="mapping")

    # elastic object

    eobject = rootNode.createChild("eobject")

    loader = eobject.createObject('MeshVTKLoader', name='loader', filename="finger.vtk")

    integration = eobject.createObject('EulerImplicitSolver', name='integration')
    solver = eobject.createObject('SparseLDLSolver', name="solver")

    container = eobject.createObject('TetrahedronSetTopologyContainer', src='@loader', name='container')
    dofs = eobject.createObject('MechanicalObject', template='Vec3d', name='dofs', src='@loader')

    geomAlgo = eobject.createObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name='GeomAlgo',
                                    drawTetrahedra=True)

    forcefield = eobject.createObject('TetrahedronFEMForceField', template='Vec3d',
                                      method='large', name='forcefield',
                                      poissonRatio=0.25, youngModulus=20000)

    constraintSolver = eobject.createObject('LinearSolverConstraintCorrection', solverName=solver.name)

    mass = eobject.createObject('DiagonalMass', massDensity=.00001, name='mass')

    # constraint
    rootNode.createObject('BilateralInteractionConstraint', template='Vec3d',
                          name="BilateralInteractionConstraint",
                          object1=rigid_node.points,
                          object2=eobject,
                          first_point=5,
                          second_point=83,
                          )

    # # force
    eobject.init()

    ### THe EFFECTIVE MASS OF THE ATTACHED BODY SEEMS TO BE ONLY 1% OF ITS ACTUAL MASS ### 
    net_mass = rigid_node.getObject("mass").totalMass + eobject.getObject("mass").totalMass/100

    rigid_node.createObject("ConstantForceField", indices=[0], force=[0, 0, 981.*net_mass, 0, 0, 0],
                            arrowSizeCoef=.1/net_mass)

    return rootNode

Thanks for your help!
Josh

[Hugo]

Hi @joshuaf

I am really sorry for this delayed reply.
I wanted to find some time to look at your simulation.

Your simulation includes constraints, solved following the Lagrange multiplier method. You are right, the force that must be applied to compensate the gravity must be scaled by 0.01. It actually corresponds to the time step.

I think this comes from the Lagrange constraint resolution which involves the time step. I need to get through the equations to provide you a complete explanation. I will get back to you next week.

Sorry again for our latency!
Best wishes,

Hugo

[Joshua]

Hi Hugo,

I really appreciate your reply, and look forward to hearing what you uncover when you get a chance.

Best,
Josh

[Joshua]

Following up, when I replace the attached soft body with a rigid one, the factor of 100 (or 1/dt based on your earlier post) disappears. If anything this leaves me more confused, but I hope it provides someone more knowledgeable with some clue what’s happening!

import Sofa
from stlib.scene import Scene
from stlib.visuals import VisualModel

import numpy as np

def createScene(rootNode):
    scene = Scene(rootNode, plugins=["SoftRobots",
                                     "SofaSparseSolver",
                                     "SofaMiscCollision",
                                     "SofaPython"])
    # scene.VisualStyle.displayFlags = 'showForceFields showBehaviorModels showInteractionForceFields'
    rootNode.createObject("BackgroundSetting", color=[1., 1., 1.])

    rootNode.dt = 0.01
    rootNode.gravity = [0., 0., -981.]

    rootNode.createObject("FreeMotionAnimationLoop")
    rootNode.createObject("GenericConstraintSolver", maxIterations=1e3, tolerance=1e-5)

    # Rigid object 1

    rigid_node1 = rootNode.createChild("rigid1")

    totalMass = .01
    inertiaMatrix = [1., 0., 0., 0., 1., 0., 0., 0., 1.]
    volume = 1.0
    uniformScale = 1
    rayleighMass = .5
    rayleighStiffness = 1

    rigid_node1.createObject('EulerImplicitSolver', name='odesolver',
                             rayleighMass=rayleighMass, rayleighStiffness=rayleighStiffness)
    rigid_node1.createObject('CGLinearSolver', name='Solver')

    rigid_node1.createObject("MechanicalObject", template="Rigid3", name="dofs", showObjectScale=uniformScale)
    rigid_mass1 = rigid_node1.createObject('UniformMass', name="mass", vertexMass=[totalMass, volume, inertiaMatrix[:]])
    rigid_mass1.totalMass = totalMass

    rigid_node1.createObject('UncoupledConstraintCorrection', compliance=1e-8)

    points1 = rigid_node1.createChild('points')
    loader1 = points1.createObject('MeshObjLoader', name='loader', filename='mesh/cube.obj', scale=uniformScale)
    points1.createObject('Mesh', src=loader1, name="MeshPoints")
    points1.createObject('MechanicalObject', name='dofs', template='Vec3d', showObject=True, showObjectScale=1,
                         src=loader1)

    points1.createObject('RigidMapping', name="mapping", input=rigid_node1, output=points1)
    points1.createObject('UncoupledConstraintCorrection', compliance=0.002)

    visual1 = VisualModel(rigid_node1, 'mesh/cube.obj', scale=[uniformScale]*3,
                          color=[0.75, 0.75, 0.25, 1.])
    visual1.writeZTransparent = True
    visual1.createObject('RigidMapping', name="mapping")

    # Rigid object 2

    rigid_node2 = rootNode.createChild("rigid2")

    rigid_node2.createObject('EulerImplicitSolver', name='odesolver',
                             rayleighMass=rayleighMass, rayleighStiffness=rayleighStiffness)
    rigid_node2.createObject('CGLinearSolver', name='Solver')

    rigid_node2.createObject("MechanicalObject", template="Rigid3", name="dofs", showObjectScale=uniformScale,
                             translation=[0, 0, -2])
    rigid_mass2 = rigid_node2.createObject('UniformMass', name="mass", vertexMass=[totalMass, volume, inertiaMatrix[:]])
    rigid_mass2.totalMass = totalMass

    rigid_node2.createObject('UncoupledConstraintCorrection', compliance=1e-8)

    points2 = rigid_node2.createChild('points')
    loader2 = points2.createObject('MeshObjLoader', name='loader', filename='mesh/cube.obj', scale=uniformScale,
                                   translation=[0, 0, 0])
    points2.createObject('Mesh', src=loader2, name="MeshPoints")
    points2.createObject('MechanicalObject', name='dofs', template='Vec3d', showObject=True, showObjectScale=1,
                         src=loader2)

    points2.createObject('RigidMapping', name="mapping", input=rigid_node2, output=points2)
    points2.createObject('UncoupledConstraintCorrection', compliance=0.002)

    visual2 = VisualModel(rigid_node2, 'mesh/cube.obj', scale=[uniformScale]*3,
                          color=[0.75, 0.75, 0.25, 1.])
    visual2.writeZTransparent = True
    visual2.createObject('RigidMapping', name="mapping")

    # constraint
    rootNode.createObject('BilateralInteractionConstraint', template='Vec3d',
                          name="BilateralInteractionConstraint",
                          object1=rigid_node1.points,
                          object2=rigid_node2.points,
                          first_point=[0, 3, 4, 7],
                          second_point=[1, 2, 5, 6],
                          )

    net_mass = rigid_node1.getObject("mass").totalMass + rigid_node2.getObject("mass").totalMass
    rigid_node2.createObject("ConstantForceField", indices=[0], force=[0, 0, 981.*net_mass, 0, 0, 0],
                             arrowSizeCoef=.01/net_mass)

    return rootNode

[Hugo]

Hi @joshuaf

Sorry for my latency, I am in the middle of different administrative things!

I really must thank you for investigating further with a full-rigid approach. It definitely helps: since it shows that the problem comes from the computation of the implicit part of the internal forces (mechanical model) with these constraints.

It is therefore a very good clue.
Let me further dig into this and I get back to you asap.

Best

Hugo

[Hugo]

Dear @joshuaf

Sorry for the delay of this reply. Your question got me stuck for some time and I needed to get to some constraint resolution experts. They actually did not get why at a first glance the forces were not balancing the gravity acceleration.

I think I just found the point.

In the UncoupledConstraintCorrection, it exist (in a galaxy far far away) a data named “useOdeSolverIntegrationFactors“. This data allows to use as factor to apply the corrective motion (coming from the lamda multipliers, solution of the constraint problem) either:

* (default) the correctionVelocityFactor and correctionPositionFactor (two other data in the UncoupledConstraintCorrection) –> useOdeSolverIntegrationFactors=0
* the integration factors of the ODE solver (EulerImplicitSolver in your case) –> useOdeSolverIntegrationFactors=1

On the other hand, the LinearSolverConstraintCorrection offers only one solution for applying the corrective motion: it is using the integration factors of the ODE solver. Therefore, if you do not specify the data useOdeSolverIntegrationFactors=1 in your UncoupledConstraintCorrection, the constraints won’t be applied in the same way on the rigid () and the deformable object.

Your fix: just setting the data to true in the UncoupledConstraintCorrection.

Damned! all this just for a flag. I will investigate more what is the interest of using correctionVelocityFactor and correctionPositionFactor, provided by the user.
Best,

Hugo

[Joshua]

Looks like that worked! Thanks so much, Hugo!

[jnbrunet]

Hey @hugo,

Nice find!

Do you think we could add some documentation (maybe here) on these integration factors? In fact, all of those methods that need to be overridden by ODE solvers:

double OdeSolver::getPositionIntegrationFactor (int /*inputDerivative*/, int /*outputDerivative*/);
double OdeSolver::getSolutionIntegrationFactor (int /*outputDerivative*/);
double OdeSolver::getVelocityIntegrationFactor();
double OdeSolver::getPositionIntegrationFactor();

are very difficult to grasp…and the doxygen isn’t super clear.

Jean-Nico

[Author]

Posts