sofa-framework
[Mapping] Implement tetrahedral-topological changes in SubsetTopologicalMapping
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
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)
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.
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.
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 https://github.com/sofa-framework/sofa/issues/5127