[SOLVED] BilateralInteractionConstraint reduces gravity by 100x

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