[Topology] Topological changes on HexahedronSetTopology and EdgeSetTopology Container

[ScheiklP]

Problem

Hi, I would like to create scene where a SphereCollisionModel cuts

• HexahedronSetTopology
• EdgeSetTopology

For the HexahedronSetTopology, there seems to be a bug in mapping the correct indices from a triangle in collision to the corresponding hexahedron through Quad2TriangleTopologicalMapping

For the EdgeSetTopology, there is no cutting happening on LineCollisionModel, probably because there is no corresponding implementation in TopologicalChangeManager::removeItemsFromCollisionModel?

The minimal scene for reproducing the behavior (also works through the mouse manager) is:

Click to expand!

import Sofa
import Sofa.Core


def createScene(
    root_node: Sofa.Core.Node,
) -> Sofa.Core.Node:

    # test_cutting = "hexa"
    test_cutting = "edges"

    plugin_list = [
        "Sofa.Component.Collision.Detection.Algorithm",  # [BVHNarrowPhase, BruteForceBroadPhase, DefaultPipeline]
        "Sofa.Component.Collision.Detection.Intersection",  # [NewProximityIntersection]
        "Sofa.Component.Collision.Geometry",  # [LineCollisionModel, PointCollisionModel]
        "Sofa.Component.Collision.Response.Contact",  # [DefaultContactManager]
        "Sofa.Component.Setting",  # [BackgroundSetting]
        "Sofa.Component.Topology.Container.Dynamic",  # [HexahedronSetTopologyContainer, HexahedronSetTopologyModifier]
        "Sofa.Component.Topology.Container.Grid",  # [RegularGridTopology]
        "Sofa.Component.Visual",  # [VisualStyle]
        "Sofa.Component.Constraint.Projective",  # [FixedConstraint]
        "Sofa.Component.LinearSolver.Iterative",  # [CGLinearSolver]
        "Sofa.Component.Mass",  # [UniformMass]
        "Sofa.Component.ODESolver.Backward",  # [EulerImplicitSolver]
        "Sofa.Component.SolidMechanics.FEM.Elastic",  # [HexahedronFEMForceField]
        "Sofa.Component.Topology.Mapping",  # [Hexa2TetraTopologicalMapping]
        "SofaCarving",  # [CarvingManager]
        "Sofa.Component.SolidMechanics.Spring",  # [StiffSpringForceField]
    ]

    plugin_node = root_node.addChild("Plugins")

    for plugin in plugin_list:
        plugin_node.addObject("RequiredPlugin", pluginName=plugin, name=plugin)

    root_node.addObject("DefaultAnimationLoop")
    root_node.addObject("DefaultVisualManagerLoop")
    root_node.addObject(
        "VisualStyle",
        displayFlags=["showVisual", "showForceFields", "showCollisionModels", "showBehaviorModels", "showInteractionForceFields"],
    )

    root_node.addObject("DefaultPipeline")
    root_node.addObject("BruteForceBroadPhase")
    root_node.addObject("BVHNarrowPhase")
    root_node.addObject("DefaultContactManager", response="PenalityContactForceField")

    root_node.addObject(
        "NewProximityIntersection",
        alarmDistance=3.0,
        contactDistance=1.1,
    )

    root_node.gravity = [0.0, -98.1, 0.0]

    scene_node = root_node.addChild("scene")

    ######################
    # Cuttable HEXA Object
    ######################
    if test_cutting == "hexa":
        hexa_node = scene_node.addChild("hexas")
        topology_creation_node = scene_node.addChild("topology_creation")

        grid = topology_creation_node.addObject("RegularGridTopology", nx=10, ny=2, nz=2, xmin=0.0, xmax=10.0, ymin=0.0, ymax=2.0, zmin=0.0, zmax=2.0)

        topology_creation_node.init()

        grid_positions = grid.position.array()
        grid_hexahedra = grid.hexahedra.array()

        hexa_node.addObject("CGLinearSolver")
        hexa_node.addObject("EulerImplicitSolver")

        hexa_node.addObject("HexahedronSetTopologyContainer", hexahedra=grid_hexahedra, position=grid_positions)
        hexa_node.addObject("HexahedronSetTopologyModifier")
        hexa_node.addObject("MechanicalObject", showObject=True, showObjectScale=2.0)
        hexa_node.addObject("UniformMass", totalMass=0.001)
        hexa_node.addObject("HexahedronFEMForceField", youngModulus=100, poissonRatio=0.0)

        # There are a total of 40 points
        hexa_node.addObject("FixedConstraint", indices=[0, 10, 20, 30, 9, 19, 29, 39])

        triangle_node = hexa_node.addChild("triangles")
        triangle_node.addObject("TriangleSetTopologyContainer")
        triangle_node.addObject("TriangleSetTopologyModifier")
        triangle_node.addObject("Quad2TriangleTopologicalMapping")
        triangle_node.addObject("TriangleCollisionModel", group=0, tags="CarvingSurface")

        cutting_sphere_node = scene_node.addChild("cutting_sphere")
        cutting_sphere_node.addObject("CGLinearSolver")
        cutting_sphere_node.addObject("EulerImplicitSolver")
        cutting_sphere_node.addObject("MechanicalObject", template="Rigid3d", position=[5.0, 5.0, 1.0, 0.0, 0.0, 0.0, 1.0])
        cutting_sphere_node.addObject("UniformMass", totalMass=0.001)
        cutting_sphere_node.addObject("SphereCollisionModel", radius=1.0, group=1, tags="CarvingTool")

    elif test_cutting == "edges":
        ######################
        # Cuttable EdgeSet Object
        ######################
        edge_node = scene_node.addChild("edges")

        edge_node.addObject("CGLinearSolver")
        edge_node.addObject("EulerImplicitSolver")

        positions = [[x, 1.0, 1.0] for x in range(15, 25)]
        edges = [[x, x + 1] for x in range(len(positions) - 1)]
        edge_node.addObject("EdgeSetTopologyContainer", position=positions, edges=edges)
        edge_node.addObject("EdgeSetTopologyModifier")

        edge_node.addObject("MechanicalObject", showObject=True, showObjectScale=2.0)
        edge_node.addObject("UniformMass", totalMass=0.001)
        springs = [[first, second, 100, 0, 0] for first, second in edges]
        edge_node.addObject("StiffSpringForceField", spring=springs)
        edge_node.addObject("FixedConstraint", indices=[0, len(positions) - 1])
        edge_node.addObject("LineCollisionModel", group=0, tags="CarvingSurface")

        cutting_sphere_node_2 = scene_node.addChild("cutting_sphere_2")
        cutting_sphere_node_2.addObject("CGLinearSolver")
        cutting_sphere_node_2.addObject("EulerImplicitSolver")
        cutting_sphere_node_2.addObject("MechanicalObject", template="Rigid3d", position=[20.0, 5.0, 1.0, 0.0, 0.0, 0.0, 1.0])
        cutting_sphere_node_2.addObject("UniformMass", totalMass=0.001)
        cutting_sphere_node_2.addObject("SphereCollisionModel", radius=1.0, group=1, tags="CarvingTool")

    root_node.addObject("CarvingManager", active=True, carvingDistance=1.0)

I would be happy to implement the bug fixes and missing functions myself, but I would need some help on where to do what. :D
Maybe @epernod and @hugtalbot could help me there? :)

---

Environment

Context

• System: Ubuntu 20.04
• Version of SOFA: master branch at commit 95e9357da1a457a2a158dfc26012b214ac0aa889
• State: "Install directory"

[epernod]

Quickly looking at your example, there is a Hexa2QuadTopologyMapping missing. You should have:

- HexahedronToplogy (volume of your mesh with FEM)
    |
    - -> QuadTopology  (surface of your mesh, only quad on the surface will be present)
           Hexa2QuadTopologyMapping
            |
            - -> TriangleTopology  (dividing each quad into 2 triangles)
                    Quad2TriangleTopologyMapping

Note that I think TriangleCollisionModel directly work with QuadTopology. operating the division of each quad in 2 internally.

Let me know if it changes the behavior

[ScheiklP]

Version A (Quad with TriangleCollision):

        hexa_node.addObject("HexahedronSetTopologyContainer", hexahedra=grid_hexahedra, position=grid_positions)
        hexa_node.addObject("HexahedronSetTopologyModifier")
        hexa_node.addObject("MechanicalObject", showObject=True, showObjectScale=2.0)
        hexa_node.addObject("UniformMass", totalMass=0.001)
        hexa_node.addObject("HexahedronFEMForceField", youngModulus=100, poissonRatio=0.0)

        hexa_node.addObject("FixedConstraint", indices=[0, 60, 120, 180, 59, 119, 179, 239])

        quad_node = hexa_node.addChild("quads")
        quad_node.addObject("QuadSetTopologyContainer")
        quad_node.addObject("QuadSetTopologyModifier")
        quad_node.addObject("Hexa2QuadTopologicalMapping")
        quad_node.addObject("TriangleCollisionModel", group=0, tags="CarvingSurface")

-> If I remove elements through the mouse, it removes the incorrect indices (but it does not crash / warn) https://user-images.githubusercontent.com/29635054/184161506-3e617100-2bc0-4923-a838-2396039f89ea.mp4

-> If I let the SphereCollisionModel fall, it does not correctly detect the collision. Maybe the index order is incorrect -> flipped triangle? https://user-images.githubusercontent.com/29635054/184161484-242044af-1fe0-4a04-982a-946be839e861.mp4

Version B (Hexa -> Quad -> Triangle with TriangleCollision): Same as before.

[hugtalbot]

To close this issue: what is expected is a message in XX2XXTopologicalMapping which checks whether the input topology includes the appropriate container