[SOLVED] Subset mapping & collisions

[timp]

Hi! I am building a scene in which I connect two deformable objects (A & B) via a bilateral interaction constraint (BIC). Since the surface on which those objects are to be connected to one another is quite complex, I use a separate mesh of this surface to easily find all the relevant points for my BIC. The surface mesh is a truncated piece of the mesh of object A. I now “connect” the surface mesh to object A and B via a subset mapping and a barycentric mapping, respectively.
So far, everything works as I want it to and the two objects deform and move together.
If, however, I add a collision model with a collision pipeline at the root node, my simulation gets stuck in the very first simulation step (-> the objects disappear and SOFA becomes unresponsive). If I use a BIC without the subset mapping (e.g., by manually picking a few points of objects A & B to be connected to each other) the simulation works as expected.
Does any of you know what could lead to the objects disappearing?
Any help is appreciated, thank you! Let me know if you need more details.

[Hugo]

Hi @timp

NaN could explain it or wrong pointer, but it is weird.
Could you please share with me the full scene (in PM if you prefer) so that I can investigate this?

An usual way to investigate would be to compile in debug and see the trace.

Best,

Hugo

[timp]

Hi @hugo

Thanks for your reply!
Here’s the full scene:

import Sofa
import os

def createScene(rootNode):

        ##########################################
        # Scene setup                            #
        ##########################################

        rootNode.createObject('RequiredPlugin', name='SofaPython')
        rootNode.createObject('RequiredPlugin', name='SofaSparseSolver')

        rootNode.createObject('VisualStyle', displayFlags='showVisualModels showBehaviorModels hideCollisionModels hideBoundingCollisionModels')
        rootNode.createObject('BackgroundSetting', color='1 1 1')
        rootNode.createObject('OglSceneFrame', style='Arrows', alignment='TopRight')

        rootNode.findData('gravity').value='0 0 0' 
        
        rootNode.createObject('FreeMotionAnimationLoop')
        rootNode.createObject('GenericConstraintSolver', maxIterations='1000', tolerance = '0.001')

        rootNode.createObject('DefaultPipeline')
        rootNode.createObject('BruteForceDetection', name='N2')
        rootNode.createObject('DefaultContactManager', response='FrictionContact', responseParams='mu=0.6')
        rootNode.createObject('LocalMinDistance', name='Proximity', alarmDistance='2.0', contactDistance='0.5')

        ##########################################
        # Cube                                   #
        ##########################################

        cube = rootNode.createChild('cube')
        cube.createObject('EulerImplicitSolver', name='odesolver_cube')
        cube.createObject('SparseLDLSolver', name='LDLsolver')

        cube.createObject('MeshVTKLoader', name='loader', filename='mesh/cube.vtk', rotation='180 0 0', translation='0 30 0')
        cube.createObject('TetrahedronSetTopologyContainer', src='@loader', name='container')
        cube.createObject('TetrahedronSetTopologyModifier')
        cube.createObject('TetrahedronSetTopologyAlgorithms', template='Vec3d')
        cube.createObject('TetrahedronSetGeometryAlgorithms', template='Vec3d')

        cube.createObject('MechanicalObject', name='MO_cube', template='Vec3d')
        cube.createObject('TetrahedronFEMForceField', template='Vec3d', name='FEM_cube', method='large', poissonRatio='0.3',  youngModulus='200', drawAsEdges='false')
        cube.createObject('UniformMass', totalMass='0.0008')
        cube.createObject('LinearSolverConstraintCorrection', solverName='LDLsolver')
        
        cube.createObject('BoxROI', name='cubeROI_fix', box='20 -2 -20 -20 2 20', drawBoxes='true', position='@MO_cube.rest_position', tetrahedra='@container.tetrahedra')
        cube.createObject('FixedConstraint', indices='@cubeROI_fix.indices')

        #cube/cubeVisu
        cubeVisu = cube.createChild('visualModel')
        cubeVisu.createObject('MeshSTLLoader', filename='mesh/cube.stl', name='loader', rotation='180 0 0', translation='0 30 0')
        cubeVisu.createObject('OglModel', name='visual', src='@loader', color='0.5 0.8 0.2 0.5')
        cubeVisu.createObject('BarycentricMapping')  

        #cube/collision
        cubeCollision = cube.createChild('collisionModel')
        cubeCollision.createObject('MeshSTLLoader', name='loader', filename='mesh/cube.stl', triangulate='1', rotation='180 0 0', translation='0 30 0')
        cubeCollision.createObject('MeshTopology', src='@loader')
        cubeCollision.createObject('MechanicalObject', src='@loader')
        cubeCollision.createObject('TTriangleModel', selfCollision='0') 
        cubeCollision.createObject('TLineModel', selfCollision='0')
        cubeCollision.createObject('TPointModel', selfCollision='0')
        cubeCollision.createObject('BarycentricMapping')  

        #cube/cubeInterfacePoints
        cubeInterPts = cube.createChild('cubeInterfacePoints')
        cubeInterPts.createObject('MeshSTLLoader', filename='mesh/cube_surface.stl', name='loader', rotation='180 0 0', translation='0 30 0')
        cubeInterPts.createObject('MeshTopology', position='@loader.position', triangles='@loader.triangles', name='mesh_cubeInterface')
        cubeInterPts.createObject('MechanicalObject', name='pointCloud_cube', template='Vec3d', showObject='1', showObjectScale='15')
        cubeInterPts.createObject('TTriangleModel')
        cubeInterPts.createObject('TLineModel')
        cubeInterPts.createObject('TPointModel')
        cubeInterPts.createObject('SubsetMapping')

        ##########################################
        # Cuboid                                 #
        ##########################################                        
        
        cuboid = rootNode.createChild('cuboid')
        cuboid.createObject('EulerImplicitSolver', name='odesolver_cuboid')
        cuboid.createObject('SparseLDLSolver', name='LDLsolver')

        cuboid.createObject('MeshVTKLoader', name='loader', filename='mesh/cuboid.vtk', rotation='0 0 90', translation='15 45 0')
        cuboid.createObject('TetrahedronSetTopologyContainer', src='@loader', name='container')
        cuboid.createObject('TetrahedronSetTopologyModifier')
        cuboid.createObject('TetrahedronSetTopologyAlgorithms', template='Vec3d')
        cuboid.createObject('TetrahedronSetGeometryAlgorithms', template='Vec3d')

        cuboid.createObject('MechanicalObject', name='MO_cuboid', template='Vec3d')
        cuboid.createObject('TetrahedronFEMForceField', template='Vec3d', name='FEM_cuboid', method='large', poissonRatio='0.3',  youngModulus='500', drawAsEdges='false')
        cuboid.createObject('UniformMass', totalMass='0.0008')
        #cuboid.createObject('UncoupledConstraintCorrection')
        cuboid.createObject('LinearSolverConstraintCorrection', solverName='LDLsolver')

        #cuboid/cuboidVisu
        cuboidVisu = cuboid.createChild('visualModel')
        cuboidVisu.createObject('MeshSTLLoader', filename='mesh/cuboid.stl', name='loader', rotation='0 0 90', translation='15 45 0')
        cuboidVisu.createObject('OglModel', name='visual', src='@loader', color='0.1 0.6 0.5 0.5')
        cuboidVisu.createObject('BarycentricMapping')  

        #cuboid/collision
        cuboidCollision = cuboid.createChild('collisionModel')
        cuboidCollision.createObject('MeshSTLLoader', name='loader', filename='mesh/cuboid.stl', triangulate='1', rotation='0 0 -90', translation='-45 45 0')
        cuboidCollision.createObject('MeshTopology', src='@loader')
        cuboidCollision.createObject('MechanicalObject', src='@loader')
        cuboidCollision.createObject('TTriangleModel', selfCollision='0') 
        cuboidCollision.createObject('TLineModel', selfCollision='0')
        cuboidCollision.createObject('TPointModel', selfCollision='0')
        cuboidCollision.createObject('BarycentricMapping')  

        #cuboid/cuboidConnection
        cuboidConn = cuboid.createChild('cuboidConnection')
        cuboidConn.createObject('MeshSTLLoader', filename='mesh/cube_surface.stl', name='loader', rotation='180 0 0', translation='0 30 0')
        cuboidConn.createObject('MeshTopology', position='@loader.position', triangles='@loader.triangles')
        cuboidConn.createObject('MechanicalObject', name='pointCloud_cuboid', template='Vec3d', showObject='1', showObjectScale='5')
        cuboidConn.createObject('TTriangleModel')
        cuboidConn.createObject('TLineModel')
        cuboidConn.createObject('TPointModel')
        cuboidConn.createObject('BarycentricMapping')

        ##########################################
        # Bilateral Interaction Constraint       #
        ##########################################
        
        id_pts=""
        for i in range(len(cubeInterPts.getObject('mesh_cubeInterface').findData('position').value)):
            id_pts= id_pts+str(i)+" "
        rootNode.createObject(
            'BilateralInteractionConstraint', 
            name='BIC', 
            template='Vec3d', 
            object1='@cube/cubeInterfacePoints/pointCloud_cube', 
            object2='@cuboid/cuboidConnection/pointCloud_cuboid', 
            first_point=id_pts, 
            second_point=id_pts,
        )
        
        
        return rootNode

Thank you!

[timp]

…oh, and here are the mesh files in case you want to run the scene.

[Hugo]

Hi @timp

I am really the sorry the transfer expired. Could you restart it please? I will try to give a look within the next days! My apologies!
Best,

Hugo

[timp]

Hi @hugo

The files are here

Thanks for your help!

[Hugo]

Hi @timp

I am sorry it took me some time to find out the issue (without mentioning the fact that I am just coming back from some days off).

The problem comes from the fact that you can not have collision computed on the attached area. Otherwise the collision and the bilateral constraint would be contradictory: one is attaching the two objects, while the other is repulsing the same two objects.

Here is a version of the scene working for my version of SOFA (latest master branch):

import Sofa
import os

def createScene(rootNode):

        ##########################################
        # Scene setup                            #
        ##########################################

        rootNode.createObject('RequiredPlugin', name='SofaPython')
        rootNode.createObject('RequiredPlugin', name='SofaSparseSolver')

        rootNode.createObject('VisualStyle', displayFlags='showVisualModels showBehaviorModels hideCollisionModels hideBoundingCollisionModels')
        rootNode.createObject('OglSceneFrame', style='Arrows', alignment='TopRight')

        rootNode.findData('gravity').value='0 0 0' 
        
        rootNode.createObject('FreeMotionAnimationLoop')
        rootNode.createObject('GenericConstraintSolver', maxIterations='1000', tolerance = '0.001')

        rootNode.createObject('DefaultPipeline')
        rootNode.createObject('BruteForceDetection', name='N2')
        rootNode.createObject('DefaultContactManager', response='FrictionContact', responseParams='mu=0.6')
        rootNode.createObject('LocalMinDistance', name='Proximity', alarmDistance='2.0', contactDistance='0.5')

        ##########################################
        # Cube                                   #
        ##########################################

        cube = rootNode.createChild('cube')
        cube.createObject('EulerImplicitSolver', name='odesolver_cube')
        cube.createObject('SparseLDLSolver', name='LDLsolver')

        cube.createObject('MeshVTKLoader', name='loader', filename='mesh/cube.vtk', rotation='180 0 0', translation='0 30 0')
        cube.createObject('TetrahedronSetTopologyContainer', src='@loader', name='container')
        cube.createObject('TetrahedronSetTopologyModifier')
        cube.createObject('TetrahedronSetGeometryAlgorithms', template='Vec3d')

        cube.createObject('MechanicalObject', name='MO_cube', template='Vec3d')
        cube.createObject('TetrahedronFEMForceField', template='Vec3d', name='FEM_cube', method='large', poissonRatio='0.3',  youngModulus='200')
        cube.createObject('UniformMass', totalMass='0.0008')
        cube.createObject('LinearSolverConstraintCorrection', solverName='LDLsolver')
        
        cube.createObject('BoxROI', name='cubeROI_fix', box='20 -2 -20 -20 2 20', drawBoxes='true', position='@MO_cube.rest_position', tetrahedra='@container.tetrahedra')
        cube.createObject('FixedConstraint', indices='@cubeROI_fix.indices')

        #cube/cubeVisu
        cubeVisu = cube.createChild('visualModel')
        cubeVisu.createObject('MeshSTLLoader', filename='mesh/cube.stl', name='loader', rotation='180 0 0', translation='0 30 0')
        cubeVisu.createObject('OglModel', name='visual', src='@loader', color='0.5 0.8 0.2 0.5')
        cubeVisu.createObject('BarycentricMapping')  

        #cube/cubeInterfacePoints
        cubeInterPts = cube.createChild('cubeInterfacePoints')
        cubeInterPts.createObject('MeshSTLLoader', filename='mesh/cube_surface.stl', name='loader', rotation='180 0 0', translation='0 30 0')
        cubeInterPts.createObject('MeshTopology', position='@loader.position', triangles='@loader.triangles', name='mesh_cubeInterface')
        cubeInterPts.createObject('MechanicalObject', name='pointCloud_cube', template='Vec3d', showObject='1', showObjectScale='15')
        cubeInterPts.createObject('SubsetMapping')

        ##########################################
        # Cuboid                                 #
        ##########################################                        
        
        cuboid = rootNode.createChild('cuboid')
        cuboid.createObject('EulerImplicitSolver', name='odesolver_cuboid')
        cuboid.createObject('SparseLDLSolver', name='LDLsolver')

        cuboid.createObject('MeshVTKLoader', name='loader', filename='mesh/cuboid.vtk', rotation='0 0 90', translation='15 45 0')
        cuboid.createObject('TetrahedronSetTopologyContainer', src='@loader', name='container')
        cuboid.createObject('TetrahedronSetTopologyModifier')
        cuboid.createObject('TetrahedronSetGeometryAlgorithms', template='Vec3d')

        cuboid.createObject('MechanicalObject', name='MO_cuboid', template='Vec3d')
        cuboid.createObject('TetrahedronFEMForceField', template='Vec3d', name='FEM_cuboid', method='large', poissonRatio='0.3',  youngModulus='500')
        cuboid.createObject('UniformMass', totalMass='0.0008')
        #cuboid.createObject('UncoupledConstraintCorrection')
        cuboid.createObject('LinearSolverConstraintCorrection', solverName='LDLsolver')

        #cuboid/cuboidVisu
        cuboidVisu = cuboid.createChild('visualModel')
        cuboidVisu.createObject('MeshSTLLoader', filename='mesh/cuboid.stl', name='loader', rotation='0 0 90', translation='15 45 0')
        cuboidVisu.createObject('OglModel', name='visual', src='@loader', color='0.1 0.6 0.5 0.5')
        cuboidVisu.createObject('BarycentricMapping')  

        #cuboid/cuboidConnection
        cuboidConn = cuboid.createChild('cuboidConnection')
        cuboidConn.createObject('MeshSTLLoader', filename='mesh/cube_surface.stl', name='loader', rotation='180 0 0', translation='0 30 0')
        cuboidConn.createObject('MeshTopology', position='@loader.position', triangles='@loader.triangles')
        cuboidConn.createObject('MechanicalObject', name='pointCloud_cuboid', template='Vec3d', showObject='1', showObjectScale='5')
        cuboidConn.createObject('BarycentricMapping')

        ##########################################
        # Bilateral Interaction Constraint       #
        ##########################################
        
        id_pts=""
        for i in range(len(cubeInterPts.getObject('mesh_cubeInterface').findData('position').value)):
            id_pts= id_pts+str(i)+" "
        rootNode.createObject(
            'BilateralInteractionConstraint', 
            name='BIC', 
            template='Vec3d', 
            object1='@cube/cubeInterfacePoints/pointCloud_cube', 
            object2='@cuboid/cuboidConnection/pointCloud_cuboid', 
            first_point=id_pts, 
            second_point=id_pts,
        )
        
        
        return rootNode

I hope this helps.
Best,

Hugo

[timp]

Hi @hugo,

Thanks for having a look at my scene! Your explanation for why the simulation crashes makes sense, I changed the scene accordingly and now it works.
Thank you!

[Hugo]

Great to read !
Thanks for having closed the topic. Feel free to open a new one anytime!

Best

Hugo

[Author]

Posts