Propagate Changes in SubsetTopologicalMapping

[ScheiklP]

Hi, while implementing TETRAHEDRAREMOVED for SubsetTopologicalMapping (https://github.com/sofa-framework/sofa/pull/5106) I saw a lot of errors in the ENDING_EVENT.

[ERROR]   [SubsetTopologicalMapping(SubsetTopologicalMapping)] Invalid d_pointS2D size : 22 != 23
[ERROR]   [SubsetTopologicalMapping(SubsetTopologicalMapping)] Invalid d_pointS2D size : 22 != 23
[ERROR]   [SubsetTopologicalMapping(SubsetTopologicalMapping)] Invalid d_triangleS2D size : 19 != 20
[ERROR]   [SubsetTopologicalMapping(SubsetTopologicalMapping)] Invalid d_triangleS2D size : 19 != 20
[ERROR]   [SubsetTopologicalMapping(SubsetTopologicalMapping)] Invalid d_pointS2D size : 21 != 22
[ERROR]   [SubsetTopologicalMapping(SubsetTopologicalMapping)] Invalid d_pointS2D size : 21 != 22

This was a bit suspicious, because I just followed the implementation of TRIANGLESREMOVED. So I created a dummy scene with only triangles, and used v24.06 of SOFA (so without any of my changes). Sadly, the error messages also exist there. And it is not only the message, there is actually something wrong with the topology containers.

Video:

https://github.com/user-attachments/assets/41c09909-a885-4c14-800b-4fd3679081dc

Example scene

plugin_list = [
    "MultiThreading",  # Needed to use components [ParallelBVHNarrowPhase,ParallelBruteForceBroadPhase]
    "Sofa.Component.AnimationLoop",  # Needed to use components [FreeMotionAnimationLoop]
    "Sofa.Component.Collision.Detection.Algorithm",  # Needed to use components [CollisionPipeline]
    "Sofa.Component.Collision.Detection.Intersection",  # Needed to use components [MinProximityIntersection]
    "Sofa.Component.Collision.Geometry",  # Needed to use components [TriangleCollisionModel]
    "Sofa.Component.Collision.Response.Contact",  # Needed to use components [CollisionResponse]
    "Sofa.Component.Constraint.Lagrangian.Correction",  # Needed to use components [LinearSolverConstraintCorrection]
    "Sofa.Component.Constraint.Lagrangian.Solver",  # Needed to use components [GenericConstraintSolver]
    "Sofa.Component.LinearSolver.Direct",  # Needed to use components [SparseLDLSolver]
    "Sofa.Component.Mass",  # Needed to use components [UniformMass]
    "Sofa.Component.ODESolver.Backward",  # Needed to use components [EulerImplicitSolver]
    "Sofa.Component.SolidMechanics.FEM.Elastic",  # Needed to use components [TetrahedralCorotationalFEMForceField]
    "Sofa.Component.StateContainer",  # Needed to use components [MechanicalObject]
    "Sofa.Component.Topology.Container.Dynamic",  # Needed to use components [TetrahedronSetTopology]
    "Sofa.Component.Visual",  # Needed to use components [VisualStyle]
    "Sofa.Component.Constraint.Projective",  # Needed to use components [FixedProjectiveConstraint]
    "Sofa.Component.Engine.Select",  # Needed to use components [BoxROI]
    "Sofa.Component.Topology.Container.Grid",  # Needed to use components [RegularGridTopology]
    "Sofa.Component.Mapping.Linear",  # Needed to use components [SubsetMapping]
    "Sofa.Component.LinearSolver.Iterative",  # Needed to use components [CGLinearSolver]
    "Sofa.Component.Topology.Mapping",  # Needed to use components [Quad2TriangleTopologicalMapping,SubsetT
]

def createScene(root):

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

    root.gravity = [0.0, -9.81, 0.0]
    root.dt = 0.01

    root.addObject("FreeMotionAnimationLoop")
    root.addObject("VisualStyle", displayFlags=["showForceFields", "showBehaviorModels", "showCollisionModels"])

    root.addObject("CollisionPipeline")
    root.addObject("ParallelBruteForceBroadPhase")
    root.addObject("ParallelBVHNarrowPhase")
    root.addObject("MinProximityIntersection", alarmDistance=5.0, contactDistance=1.0)

    root.addObject("CollisionResponse", response="FrictionContactConstraint")
    root.addObject("GenericConstraintSolver")

    scene_node = root.addChild("scene")

    topology_node = scene_node.addChild("topology")
    topology_node.addObject("RegularGridTopology", nx=5, ny=5, nz=1, xmin=-10.0, xmax=10.0, ymin=-10.0, ymax=10.0, zmin=0.0, zmax=0.0)
    topology_node.addObject("TriangleSetTopologyContainer")
    topology_node.addObject("TriangleSetTopologyModifier")
    topology_node.addObject("Quad2TriangleTopologicalMapping", input="@RegularGridTopology", output="@TriangleSetTopologyContainer")

    topology_node.init()
    positions = topology_node.RegularGridTopology.position.array()
    triangles = topology_node.TriangleSetTopologyContainer.triangles.array()

    first_positions = [point.tolist() for point in positions if point[0] <= 0.1]
    first_triangles = []
    for triangle in triangles:
        triangle_positions = [positions[index] for index in triangle]
        if all([point[0] <= 0.1 for point in triangle_positions]):
            # find the corresponding indices in the first_positions
            new_triangle = [first_positions.index(point.tolist()) for point in triangle_positions]
            first_triangles.append(new_triangle)

    first_subset_indices = []
    for index, point in enumerate(positions):
        if point[0] <= 0.1:
            first_subset_indices.append(index)

    second_positions = [point.tolist() for point in positions if point[0] >= -0.1]
    second_triangles = []
    for triangle in triangles:
        triangle_positions = [positions[index] for index in triangle]
        if all([point[0] >= -0.1 for point in triangle_positions]):
            # find the corresponding indices in the second_positions
            new_triangle = [second_positions.index(point.tolist()) for point in triangle_positions]
            second_triangles.append(new_triangle)

    second_subset_indices = []
    for index, point in enumerate(positions):
        if point[0] >= -0.1:
            second_subset_indices

    hybrid_object_node = scene_node.addChild("hybrid_object")
    hybrid_object_node.addObject("TriangleSetTopologyContainer", position=positions, triangles=triangles)
    hybrid_object_node.addObject("TriangleSetTopologyModifier")
    hybrid_object_node.addObject("EulerImplicitSolver")
    hybrid_object_node.addObject("SparseLDLSolver", template="CompressedRowSparseMatrixMat3x3d")
    hybrid_object_node.addObject("MechanicalObject", showObject=True)
    hybrid_object_node.addObject("BoxROI", box=[-10.0, 9.0, -10.0, 10.0, 10.0, 10.0])
    hybrid_object_node.addObject("FixedProjectiveConstraint", indices="@BoxROI.indices")
    hybrid_object_node.addObject("TriangleCollisionModel")
    hybrid_object_node.addObject("LinearSolverConstraintCorrection")

    first_part_node = hybrid_object_node.addChild("first_part")
    first_part_node.addObject(
        "TriangleSetTopologyContainer",
        position=first_positions,
        triangles=first_triangles,
    )
    first_part_node.addObject("TriangleSetTopologyModifier")
    first_part_node.addObject("TriangleSetGeometryAlgorithms", template="Vec3")
    first_part_node.addObject(
        "SubsetTopologicalMapping",
        input="@../TriangleSetTopologyContainer",
        output="@TriangleSetTopologyContainer",
        samePoints=False,
        handleTriangles=True,
    )
    first_part_node.addObject("MechanicalObject")
    first_part_node.addObject("TriangleFEMForceField", youngModulus=1000, poissonRatio=0.49)
    first_part_node.addObject("UniformMass", totalMass=1.0)
    # first_part_node.addObject("TriangleCollisionModel")
    first_part_node.addObject("SubsetMapping", indices=first_subset_indices, handleTopologyChange=True)

    second_part_node = hybrid_object_node.addChild("second_part")
    second_part_node.addObject(
        "TriangleSetTopologyContainer",
        position=second_positions,
        triangles=second_triangles,
    )
    second_part_node.addObject("TriangleSetTopologyModifier")
    second_part_node.addObject("TriangleSetGeometryAlgorithms", template="Vec3")
    second_part_node.addObject(
        "SubsetTopologicalMapping",
        input="@../TriangleSetTopologyContainer",
        output="@TriangleSetTopologyContainer",
        samePoints=False,
        handleTriangles=True,
    )
    second_part_node.addObject("MechanicalObject")
    second_part_node.addObject("TriangleFEMForceField", youngModulus=1000, poissonRatio=0.49)
    # second_part_node.addObject("TriangleCollisionModel")
    second_part_node.addObject("UniformMass", totalMass=1.0)
    second_part_node.addObject("SubsetMapping", indices=second_subset_indices, handleTopologyChange=True)

[bakpaul]

Let's isolate the problem. if you remove all the subet topological mapping and only have a simple triangleSetTopology (only keep the hybrid object), does this still happens ?

[ScheiklP]

Oh now it is even worse. :D It just segfaults.

Adapted scene

plugin_list = [
    "MultiThreading",  # Needed to use components [ParallelBVHNarrowPhase,ParallelBruteForceBroadPhase]
    "Sofa.Component.AnimationLoop",  # Needed to use components [FreeMotionAnimationLoop]
    "Sofa.Component.Collision.Detection.Algorithm",  # Needed to use components [CollisionPipeline]
    "Sofa.Component.Collision.Detection.Intersection",  # Needed to use components [MinProximityIntersection]
    "Sofa.Component.Collision.Geometry",  # Needed to use components [TriangleCollisionModel]
    "Sofa.Component.Collision.Response.Contact",  # Needed to use components [CollisionResponse]
    "Sofa.Component.Constraint.Lagrangian.Correction",  # Needed to use components [LinearSolverConstraintCorrection]
    "Sofa.Component.Constraint.Lagrangian.Solver",  # Needed to use components [GenericConstraintSolver]
    "Sofa.Component.LinearSolver.Direct",  # Needed to use components [SparseLDLSolver]
    "Sofa.Component.Mass",  # Needed to use components [UniformMass]
    "Sofa.Component.ODESolver.Backward",  # Needed to use components [EulerImplicitSolver]
    "Sofa.Component.SolidMechanics.FEM.Elastic",  # Needed to use components [TetrahedralCorotationalFEMForceField]
    "Sofa.Component.StateContainer",  # Needed to use components [MechanicalObject]
    "Sofa.Component.Topology.Container.Dynamic",  # Needed to use components [TetrahedronSetTopology]
    "Sofa.Component.Visual",  # Needed to use components [VisualStyle]
    "Sofa.Component.Constraint.Projective",  # Needed to use components [FixedProjectiveConstraint]
    "Sofa.Component.Engine.Select",  # Needed to use components [BoxROI]
    "Sofa.Component.Topology.Container.Grid",  # Needed to use components [RegularGridTopology]
    "Sofa.Component.Mapping.Linear",  # Needed to use components [SubsetMapping]
    "Sofa.Component.LinearSolver.Iterative",  # Needed to use components [CGLinearSolver]
    "Sofa.Component.Topology.Mapping",  # Needed to use components [Quad2TriangleTopologicalMapping,SubsetT
]

def createScene(root):

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

    root.gravity = [0.0, -9.81, 0.0]
    root.dt = 0.01

    root.addObject("FreeMotionAnimationLoop")
    root.addObject("VisualStyle", displayFlags=["showForceFields", "showBehaviorModels"])

    root.addObject("CollisionPipeline")
    root.addObject("ParallelBruteForceBroadPhase")
    root.addObject("ParallelBVHNarrowPhase")
    root.addObject("MinProximityIntersection", alarmDistance=5.0, contactDistance=1.0)

    root.addObject("CollisionResponse", response="FrictionContactConstraint")
    root.addObject("GenericConstraintSolver")

    scene_node = root.addChild("scene")

    topology_node = scene_node.addChild("topology")
    topology_node.addObject("RegularGridTopology", nx=5, ny=5, nz=1, xmin=-10.0, xmax=10.0, ymin=-10.0, ymax=10.0, zmin=0.0, zmax=0.0)
    topology_node.addObject("TriangleSetTopologyContainer")
    topology_node.addObject("TriangleSetTopologyModifier")
    topology_node.addObject("Quad2TriangleTopologicalMapping", input="@RegularGridTopology", output="@TriangleSetTopologyContainer")

    topology_node.init()
    positions = topology_node.RegularGridTopology.position.array()
    triangles = topology_node.TriangleSetTopologyContainer.triangles.array()

    hybrid_object_node = scene_node.addChild("hybrid_object")
    hybrid_object_node.addObject("TriangleSetTopologyContainer", position=positions, triangles=triangles)
    hybrid_object_node.addObject("TriangleSetTopologyModifier")
    hybrid_object_node.addObject("TriangleSetGeometryAlgorithms")
    hybrid_object_node.addObject("EulerImplicitSolver")
    hybrid_object_node.addObject("SparseLDLSolver", template="CompressedRowSparseMatrixMat3x3d")
    hybrid_object_node.addObject("MechanicalObject", showObject=True)
    hybrid_object_node.addObject("TriangleFEMForceField", youngModulus=1000, poissonRatio=0.49)
    hybrid_object_node.addObject("BoxROI", box=[-10.0, 9.0, -10.0, 10.0, 10.0, 10.0])
    hybrid_object_node.addObject("FixedProjectiveConstraint", indices="@BoxROI.indices")
    hybrid_object_node.addObject("UniformMass", totalMass=1.0)
    hybrid_object_node.addObject("TriangleCollisionModel")
    hybrid_object_node.addObject("LinearSolverConstraintCorrection")

https://github.com/user-attachments/assets/8f742837-94ea-424c-8ced-e2f3a0adf2b4

[ScheiklP]

Ah, but it is actually the TriangleForcField's fault. Does not segfault with MeshSpringForceField. So maybe that is another bug to put in the backlog. :D

Without the topological subset mapping (just hybrid) there are no issues. However, I think not every edge that should be removed is actually removed.

As you can see, on the right, there is still an edge, that should have been removed, when the triangle was removed.

https://github.com/user-attachments/assets/072aefa8-c1ca-4142-ac48-458d506b1cbe