[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.

[epernod]

wow wow wow.... ok I have to admit I've never used this component and wish I had never looked at it!
So for me your code is correct in respect to the rest of the component... but for me the component should be totally rewritten.

I don't understand why in the init method we retrieve the same tetrahedra to create the map. For me this subsetTopologyMapping should feed the toModel topology using input tetrahedra and then handle change on this list.

The component is inheriting from public sofa::core::topology::TopologicalMapping but is not using its maps

    // Map which gives for each index (global index) of an element in the INPUT topology
    // the corresponding index (local index) of the same element in the OUTPUT topology :
    std::map<Index, Index> Glob2LocMap;

    std::map<Index, sofa::type::vector<Index> > In2OutMap;

which is the main concept of topological mappings.

You code is working when you remove a Tetrahedron in the input topology and also present in the output topology but in fact you should also handle changes when a tetrahedron is remove from the input topology and NOT present in the output topology. Indeed in this case the buffer of tetrahedra is renumbered in the input topology and you need to propagate those reumbering in the maps of your mapping.

[epernod]

So regarding your PR, does it works for you case?
if yes, I would suggest to merge it and open an issue to keep in mind that this component should be totally changed.

[epernod]

I would suggest to use tetra_source_id to show that we are working on a set of indices and not on tetrahedra arrays

[damienmarchal]

remove commented code ?

[ScheiklP]

So regarding your PR, does it works for you case? if yes, I would suggest to merge it and open an issue to keep in mind that this component should be totally changed.

Hi @epernod sadly it does not work for me. :(
I still get the initial tetra remove errors from #5127