Collision rigid-deformable objects

[etagliabue]

Hi everybody,

I have been struggling with some issues related to collisions for a while, I guess since I still do not have a complete understanding on how things work in SOFA.
What I am trying to do is to simulate the interaction between a probe (rigid object) and a breast (deformable object).

1) I decided to start from a simplified scenario, and I can simulate the collision of a rigid sphere (simplified “probe”) with a deformable hemisphere with breast properties (“breast”). Here is the xml I am using:

<Node name="root" dt="0.005" gravity="0 0 0">
    <VisualStyle displayFlags="showCollisionModels hideVisualModels" />
    <DefaultAnimationLoop />
    <CollisionPipeline verbose="0" draw="0"/>
    <BruteForceDetection name="N2"/>
    <NewProximityIntersection name="Proximity" alarmDistance="1" contactDistance="0.4"/>
    <CollisionResponse name="Response" response="default"/>

    <MeshGmshLoader name="meshLoader" filename="Data/emisphere.msh"/>
    <Node name="Hemisphere">
        <EulerImplicitSolver rayleighStiffness="0"/>
        <CGLinearSolver iterations="200" tolerance="1e-06" threshold="1e-06"/>
        <TetrahedronSetTopologyContainer name="topo" src="@../meshLoader" />
        <TetrahedronSetGeometryAlgorithms template="Vec3d" name="GeomAlgo" />
        <MechanicalObject template="Vec3d" name="myLiver" showObject="1"/>
        <TetrahedronFEMForceField name="FEM" youngModulus="1681" poissonRatio="0.49"  method="large" />
        <DiagonalMass name="myDiagonalMass" massDensity="1" />
        <UniformVelocityDampingForceField template="Vec3d" name="uniformVelocityDampingFF0" />
        <BoxROI name="ROI" box="-10 -0.5 -10 10 0.5 10" drawBoxes="0" drawSize="0"/>
        <FixedConstraint name="Boundaries" indices="@ROI.indices"/>   
        <Node name="Collision">
            <MeshGmshLoader name="meshLoader" filename="Data/emisphere.msh" />
            <Mesh src="@meshLoader"/>
            <MechanicalObject scale="1.0"/>
            <Triangle name="CollisionModel" contactStiffness="100"/>
            <BarycentricMapping name="CollisionMapping" />
        </Node>
    </Node>

    <Node name="Sphere">
        <EulerImplicitSolver rayleighStiffness="0"/>
        <CGLinearSolver iterations="200" tolerance="1e-06" threshold="1e-06"/>
        <MechanicalObject template="Rigid" name="myParticle" position="-2 10 0    0 0 0 1" showObject="1" />
        <UniformMass totalmass="1" />
        <Sphere name="Floor" listRadius="1" simulated="1" moving="1" contactStiffness="100" />
    </Node>

</Node> 

First basic question: is there a way of imposing that when the simulation starts, the rigid sphere displaces from its current position to another one and freezes there? I was thinking about setting an initial velocity and run the simulation for a fixed number of steps, but maybe there’s another way.

2) However, when I try to make the scenario more complex, things stop working.
In particular: if I keep the deformable hemisphere and add the probe model with these lines:

<Node name="Probe">
        <EulerImplicitSolver rayleighStiffness="0"/>
        <CGLinearSolver iterations="200" tolerance="1e-06" threshold="1e-06"/>
        <MeshGmshLoader name="meshLoader" filename="Data/probe.msh" />
        <Mesh src="@meshLoader"/>
        <MechanicalObject name="myParticle" showObject="1" scale3d="0.1 0.1 0.1" translation="0 8 0" rotation="0 0 0" />
        <UniformMass totalmass="10" />
        <Node name="CollisionNode">
            <MeshSTLLoader name="meshLoader2" filename="Data/probe.stl" />
            <Mesh src="@meshLoader2"/>
            <MechanicalObject scale3d="0.1 0.1 0.1" translation="@../myParticle.translation" rotation="@../myParticle.rotation" position="@../myParticle.position" />
<!-- I tried to add all these 3 contact models since the probe has vertices just on its sides, but I need collision detection also along faces and edges -->
            <Triangle name="ProbeCollision" contactStiffness="1000"/>
            <LineModel name="Lines" contactStiffness="1000" />
            <PointModel name="Points" contactStiffness="1000" />           
        </Node>
    </Node>

The contact is detected but the probe disappears from the scene immediately after collision. What may be the cause? Am I not defining the collision model properly?

3) I also tried to do the opposite: keep the rigid sphere and use the MSH of the breast (it has 700 vertices) in place of the hemisphere. Collision is detected but the deformation induced to the breast cannot be visually seen (even though the elastic properties are the same as those of the hemisphere). Which could be the cause of this? Here is the code for loading the breast:

<MeshGmshLoader name="meshLoader" filename="Data/Breast.msh"/>
    <Node name="Breast">
        <EulerImplicitSolver rayleighStiffness="0"/>
        <CGLinearSolver iterations="200" tolerance="1e-06" threshold="1e-06"/>
        <TetrahedronSetTopologyContainer name="topo" src="@../meshLoader" />
        <TetrahedronSetGeometryAlgorithms template="Vec3d" name="GeomAlgo" />
        <MechanicalObject template="Vec3d" name="myLiver" showObject="1"/>
        <TetrahedronFEMForceField name="FEM" youngModulus="1681" poissonRatio="0.49"  method="large" />
        <DiagonalMass name="myDiagonalMass" massDensity="1" />
        <UniformVelocityDampingForceField template="Vec3d" name="uniformVelocityDampingFF0" />
        <FixedConstraint name="Boundaries" indices="0 10 11 12 24 25 26 27 43 55 63 65 66 67 76 77 78 86 87 88 89 90 91 96 98 100 "/>         
        <Node name="Collision">
            <MeshGmshLoader name="meshLoader" filename="Data/Breast.msh" />
            <Mesh src="@meshLoader"/>
            <MechanicalObject scale="1.0"/>
            <!-- I tried to add Line and Point models as well but nothing changed -->
            <Triangle name="CollisionModel" contactStiffness="200"/>
            <BarycentricMapping name="CollisionMapping" />
        </Node>
    </Node>

Do you have any suggestions to have the whole simulation working properly?
Sorry for my multiple questions, and thank you in advance for any help!

Thank you!
Eleonora

[Hugo]

Hi Eleonora,

Welcome to the SOFA forum and into the community!
I will shortly reply point by point but I am still investigating:

1) to constraint a motion of your object, I would advise to use the LinearMovementConstraint : this class makes possible to impose displacements depending on the simulation steps/time. You will find an example in examples/Components/constraint/

2/3) I am currently investigating on this and get back to you as soon as I have info.

Best,

Hugo

[Hugo]

Hi Eleonora !

Thank you for your patience. I really think your project is promising.
I dug into your issue with your scenes and data. There is actually general remarks:

• first, note that your mesh of the probe is not properly meshed (stretched, non-regular triangles), you need a good and regular mesh for collision.
• second, your breast mesh (Breast.msh) is in strange unit. You need to be careful with units, since if you use a breast in meters then you can use the usual SI unit (Pascal Pa) for the Young modulus. If not in meter, you need to adapt the Young modulus.

Note that in your case I think that :

• the contact stiffness was really high and to run your simulation without instabilities you needed to decrease the time step.
• moreover, you can simplify your simulation by considering the probe as rigid (with regards to the breast)

Here is a suggestion:

<Node name="root" dt="0.002" gravity="0 -9.8 0">
    <VisualStyle displayFlags="hideCollisionModels showVisualModels" />
    <DefaultAnimationLoop />
    <CollisionPipeline verbose="0" draw="0"/>
    <BruteForceDetection name="N2"/>
    <NewProximityIntersection name="Proximity" alarmDistance="1" contactDistance="0.4"/>
    <CollisionResponse name="Response" response="default"/>

    <MeshGmshLoader name="meshLoader" filename="Data/emisphere.msh"/>
    <Node name="Hemisphere">
        <EulerImplicitSolver rayleighStiffness="0"/>
        <CGLinearSolver iterations="200" tolerance="1e-09" threshold="1e-09"/>
        <TetrahedronSetTopologyContainer name="topo" src="@../meshLoader" />
        <TetrahedronSetGeometryAlgorithms template="Vec3d" name="GeomAlgo" />
        <MechanicalObject template="Vec3d" name="myLiver" showObject="1"/>
        <TetrahedronFEMForceField name="FEM" youngModulus="1681" poissonRatio="0.49"  method="large" />
        <DiagonalMass name="myDiagonalMass" massDensity="1" />
        <UniformVelocityDampingForceField template="Vec3d" name="uniformVelocityDampingFF0" />
        <BoxROI name="ROI" box="-10 -0.5 -10 10 0.5 10" drawBoxes="0" drawSize="0"/>
        <FixedConstraint name="Boundaries" indices="@ROI.indices"/>  
        <Node name="Collision">
            <MeshGmshLoader name="meshLoader" filename="Data/emisphere.msh" />
            <Mesh src="@meshLoader"/>
            <MechanicalObject scale="1.0"/>
            <Triangle name="CollisionModel" contactStiffness="100"/>
            <BarycentricMapping name="CollisionMapping" />
            <Node name="VisualNode">
                <OglModel />
                <IdentityMapping name="VisualMapHemiSphere" input="@../" output="@." />
            </Node>
        </Node>
    </Node>

    <!-- Probe -->
    <Node name="Probe">
        <EulerImplicitSolver rayleighStiffness="0"/>
        <CGLinearSolver iterations="200" tolerance="1e-09" threshold="1e-09"/>
        <MechanicalObject name="myParticle" template="Rigid" showObject="1" scale3d="0.1 0.1 0.1" position="0 0 0 0 0 0 0" translation="0 8 0" rotation="0 0 0" />
        <UniformMass totalmass="2" />
        <LinearMovementConstraint template="Rigid" keyTimes="0 0.05 0.2 10" movements="  0 0 0   0 0 0
                                                                                      0 0 0   0 0 0
                                                                                      0 -1.3 0  0 0 0
                                                                                      0 -1.3 0  0 0 0" />
        <Node name="CollisionNode">
            <MeshObjLoader name="meshLoader2" filename="Data/probe.obj" />
            <Mesh src="@meshLoader2"/>
            <MechanicalObject scale3d="0.1 0.1 0.1" name="CollModel" />
            <Triangle name="ProbeCollision" contactStiffness="10"/>
            <LineModel name="Lines" contactStiffness="10" />
            <PointModel name="Points" contactStiffness="10" />
            <RigidMapping name="RigidMap" input="@../" output="@CollModel" />

            <Node name="VisualNode">
                <OglModel />
                <IdentityMapping name="VisualMapProbe" input="@../" output="@." />
            </Node>

        </Node>
    </Node>

</Node>

I hope this helps

Hugo

[etagliabue]

Hi Hugo,

Thank you for your answer. I tried the solution you posted and it works correctly. Anyway, if I replace the hemisphere with the breast mesh, SOFA crashes as soon as collision with the probe is detected. I tried to consider as collision model just the subset of the mesh lying close to the probe, and I tried some different models (at different resolutions) but the behaviour is always the same. Do you think this has to deal with the breast mesh model itself? It doesn’t seem that irregular to me.

Thank you!
Eleonora

[Hugo]

Hi Eleonora,

There is clearly a question of units.
What is the size of your breast mesh ? Is it in millimeter ?
What is the size of the probe ?

Everything should use the same units. Moreover, if your meshes are in millimeter, your mass is in kg, then : the young modulus (~stiffness) should be give in MPa = N/mm2.

I will try to make your scene (breast+probe) work.
Best,

Hugo

[etagliabue]

Hi there,

Both the breast and the probe are in mm.
However, if I use a Young modulus in MPa (~ 0.01), the breast does not behave as I expect (when deformed only by gravity). I will try to think about this measurement units issue and let you know if I make any progress.

Thank you again.
Eleonora

[Hugo]

Hi Eleonora,

What units did you use for the gravity (SI units are m/s²), since you use mm it should also be coherently set.

Best,

Hugo

[etagliabue]

Hi Hugo,

Thanks for your suggestion. I’m still working on it, I’ll let you know if I solve the issue.
I have another question: can you confirm that the TetrahedralCorotationalFEMForceField assumes that my model has a linear elastic behaviour? Is there a way to use other constitutive models in SOFA?

Many thanks!
Eleonora

[Hugo]

Hi Eleonora,

That’s a good question : TetrahedralCorotationalFEMForceField implements a linear elastic constitutive law indeed. However, it is decomposing the actual element configuration into rotational and deformational components. Larger displacement (with small deformations) can therefore be handled still using a linear elastic material. More info here.

Other materials are available in SOFA. You can take a look at the TetrahedronHyperelasticityFEMForceField. This class implements hyperelastic materials.

Best regards,

Hugo

[Author]

Posts