[Mapping] Implement tetrahedral-topological changes in SubsetTopologicalMapping

[ScheiklP]

Bonjour,
I am currently adding TETRA topology changes in SubsetTopologicalMapping.

This first commit is just a WIP commit for the REMOVE case to get your feedback.
If the style and content is ok, I will also implement the ADD case.
In particular, I am unsure about the line
toTetrahedronMod->removeTetrahedra(tetra_array_buffer_destination, false);
The removeTriangle distinguishes between different topology types of isolated items, while removeTetrahedra does not.
Should this thus just set to be false?

I also think that the naive implementation of following the TRIANGLESREMOVED is invalid.

[ERROR]   [SubsetTopologicalMapping(SubsetTopologicalMapping)] Invalid d_triangleS2D size : 13943 != 13944
[ERROR]   [SubsetTopologicalMapping(SubsetTopologicalMapping)] Invalid d_triangleS2D size : 13942 != 13944
[ERROR]   [SubsetTopologicalMapping(SubsetTopologicalMapping)] Invalid d_tetrahedronS2D size : 6315 != 6314
[ERROR]   [SubsetTopologicalMapping(SubsetTopologicalMapping)] Invalid d_triangleS2D size : 13942 != 13943
[ERROR]   [SubsetTopologicalMapping(SubsetTopologicalMapping)] Invalid d_tetrahedronS2D size : 6315 != 6314
[ERROR]   [SubsetTopologicalMapping(SubsetTopologicalMapping)] Invalid d_triangleS2D size : 13942 != 13943
[ERROR]   [SubsetTopologicalMapping(SubsetTopologicalMapping)] Invalid d_tetrahedronS2D size : 6315 != 6314

Cheers,
Paul

[ScheiklP]

Test 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]
]


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=5, xmin=-10.0, xmax=10.0, ymin=-10.0, ymax=10.0, zmin=-10.0, zmax=10.0)
    topology_node.addObject("TetrahedronSetTopologyContainer")
    topology_node.addObject("TetrahedronSetTopologyModifier")
    topology_node.addObject("Hexa2TetraTopologicalMapping", input="@RegularGridTopology", output="@TetrahedronSetTopologyContainer")

    topology_node.init()
    positions = topology_node.RegularGridTopology.position.array()
    tetrahedra = topology_node.TetrahedronSetTopologyContainer.tetrahedra.array()

    first_positions = [point.tolist() for point in positions if point[1] <= 0.1]
    first_tetrahedra = []
    for tetra in tetrahedra:
        tetra_positions = [positions[index] for index in tetra]
        if all([point[1] <= 0.1 for point in tetra_positions]):
            # find the corresponding indices in the first_positions
            new_tetra = [first_positions.index(point.tolist()) for point in tetra_positions]
            first_tetrahedra.append(new_tetra)

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

    second_positions = [point.tolist() for point in positions if point[1] >= -0.1]
    second_tetrahedra = []
    for tetra in tetrahedra:
        tetra_positions = [positions[index] for index in tetra]
        if all([point[1] >= -0.1 for point in tetra_positions]):
            # find the corresponding indices in the second_positions
            new_tetra = [second_positions.index(point.tolist()) for point in tetra_positions]
            second_tetrahedra.append(new_tetra)

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

    hybrid_object_node = scene_node.addChild("hybrid_object")
    hybrid_object_node.addObject("TetrahedronSetTopologyContainer", position=positions, tetrahedra=tetrahedra)
    hybrid_object_node.addObject("TetrahedronSetTopologyModifier")
    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, -10.0, -10.0, 10.0, -9.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(
        "TetrahedronSetTopologyContainer",
        position=first_positions,
        tetrahedra=first_tetrahedra,
    )
    first_part_node.addObject("TetrahedronSetTopologyModifier")
    first_part_node.addObject("TetrahedronSetGeometryAlgorithms", template="Vec3")
    first_part_node.addObject(
        "SubsetTopologicalMapping",
        input="@../TetrahedronSetTopologyContainer",
        output="@TetrahedronSetTopologyContainer",
        samePoints=False,
        handleTriangles=True,
        handleTetrahedra=True,
    )
    first_part_node.addObject("MechanicalObject")
    first_part_node.addObject("TetrahedralCorotationalFEMForceField", youngModulus=1000, poissonRatio=0.3)
    first_part_node.addObject("UniformMass", totalMass=1.0)
    first_part_node.addObject("SubsetMapping", indices=first_subset_indices, handleTopologyChange=True)

    second_part_node = hybrid_object_node.addChild("second_part")
    second_part_node.addObject(
        "TetrahedronSetTopologyContainer",
        position=second_positions,
        tetrahedra=second_tetrahedra,
    )
    second_part_node.addObject("TetrahedronSetTopologyModifier")
    second_part_node.addObject("TetrahedronSetGeometryAlgorithms", template="Vec3")
    second_part_node.addObject(
        "SubsetTopologicalMapping",
        input="@../TetrahedronSetTopologyContainer",
        output="@TetrahedronSetTopologyContainer",
        samePoints=False,
        handleTriangles=True,
        handleTetrahedra=True,
    )
    second_part_node.addObject("MechanicalObject")
    second_part_node.addObject("TetrahedralCorotationalFEMForceField", youngModulus=1000, poissonRatio=0.3)
    second_part_node.addObject("UniformMass", totalMass=1.0)
    second_part_node.addObject("SubsetMapping", indices=second_subset_indices, handleTopologyChange=True)

[bakpaul]

It feels like you are iterating twice over this array although you could have done it only once latter by appending the tetra_array_buffer_destination in a scoped if line 997

[ScheiklP]

Sure, I can change that. However, right now the changes are not applied correctly.
Do you see anything obvious that is missing? See error messages above.

[bakpaul]

I don't know if that is the reason, but you are working with the datas tS2D ans tD2S but those are for triangles. You should be modifying the datas teS2D and teD2S

[ScheiklP]

Ohhhhh, I see. That makes much more sense. :D
Thanks! :)

[ScheiklP]

That solved some of the issues, but uncovered a much larger issue underneath. :D
See #5127

[bakpaul]

Ok, thank you about that, we will take a look, this seems quite serious indeed.