sofa-framework
HexahedronCompositeFEMForceFieldAndMass Solver
Problem
Hi, the example for the HexahedronCompositeFEMForceFieldAndMass complains, that
[WARNING] [HexahedronCompositeFEMMapping(HexahedronCompositeFEMMapping)] This object only support Direct Solving but an Indirect Solver in the scene is calling method applyJT(constraint) which is
not implemented. This will produce un-expected behavior.
Even when using the SparseLDLSolver instead of the CGLinearSolver.
Further, when I try to use the component with the liver model, with an embedded sphere (same components, but with the Free motion animation loop), the model just diverges.
SOFA: v23.12
Cheers, Paul
Hi,
The message is located here: https://github.com/sofa-framework/sofa/blob/b5aefb456472e05579913a806ee88bd171457f7a/Sofa/Component/SolidMechanics/FEM/NonUniform/src/sofa/component/solidmechanics/fem/nonuniform/HexahedronCompositeFEMMapping.h#L105
I don't think that the message is true. @epernod Do you remember why you added this message in https://github.com/sofa-framework/sofa/commit/0170a07db4edefb9d82c3c90ca9aa48e09da4e20#diff-a123e50be73be01f8474ca9ef41fb5da89fecf57c07cde10d88cfb39ec8f0f2a?
This function is called even with a direct solver. Here for example.
According to the history of this function in the component HexahedronCompositeFEMMapping, I would say that nobody took the time to implement it. Unfortunately, without this function, the mapping will not be able to map stiffness matrices. In the case of your example, the only contribution of a mapped stiffness matrix would come from the collision. It means that the collision forces don't contribute to the stiffness matrix, i.e. they are considered explicit. The scene may behave correctly even in this case, especially if you don't have collision.
If the message bothers you, and if you build SOFA yourself, you can comment it as a temporary solution.
Unfortunately, I don't know the components of your scene, so I don't know why it diverges with the liver. You can always try the usual tricks: decrease the time step size, increase the nb of CG iterations, increase the number of elements etc. Does it work with DefaultAnimationLoop? Perhaps the function applyJT is important after all?..
Alex
Hi @alxbilger , I ported the scene to python and the FreeMotionAnimation loop. Some weird observations:
- the TriangleCollisionModels did not work at all -> had to switch to SphereCollisionModel
- the collision of the super heavy ball with the composite object does nothing to the deformation
import Sofa
import Sofa.Core
def createScene(root):
root.gravity = [0, -700, 0]
root.dt = 0.02
root.addObject("RequiredPlugin", name="Sofa.Component.Collision.Detection.Algorithm")
root.addObject("RequiredPlugin", name="Sofa.Component.Collision.Detection.Intersection")
root.addObject("RequiredPlugin", name="Sofa.Component.Collision.Geometry")
root.addObject("RequiredPlugin", name="Sofa.Component.Collision.Response.Contact")
root.addObject("RequiredPlugin", name="Sofa.Component.Constraint.Projective")
root.addObject("RequiredPlugin", name="Sofa.Component.Engine.Select")
root.addObject("RequiredPlugin", name="Sofa.Component.IO.Mesh")
root.addObject("RequiredPlugin", name="Sofa.Component.LinearSolver.Iterative")
root.addObject("RequiredPlugin", name="Sofa.Component.Mapping.Linear")
root.addObject("RequiredPlugin", name="Sofa.Component.ODESolver.Backward")
root.addObject("RequiredPlugin", name="Sofa.Component.SolidMechanics.FEM.NonUniform")
root.addObject("RequiredPlugin", name="Sofa.Component.StateContainer")
root.addObject("RequiredPlugin", name="Sofa.Component.Topology.Container.Constant")
root.addObject("RequiredPlugin", name="Sofa.Component.Topology.Container.Grid")
root.addObject("RequiredPlugin", name="Sofa.Component.Visual")
root.addObject("RequiredPlugin", name="Sofa.GL.Component.Rendering3D")
root.addObject("RequiredPlugin", name="Sofa.Component.LinearSolver.Direct")
root.addObject("RequiredPlugin", name="Sofa.Component.AnimationLoop") # Needed to use components [FreeMotionAnimationLoop]
root.addObject("RequiredPlugin", name="Sofa.Component.Constraint.Lagrangian.Correction") # Needed to use components [GenericConstraintCorrection,LinearSolverConstraintCorrection]
root.addObject("RequiredPlugin", name="Sofa.Component.Constraint.Lagrangian.Solver") # Needed to use components [GenericConstraintSolver]
root.addObject("RequiredPlugin", name="Sofa.Component.Mass") # Needed to use components [UniformMass]
root.addObject("RequiredPlugin", name="Sofa.Component.Mapping.NonLinear") # Needed to use components [RigidMapping]
root.addObject("FreeMotionAnimationLoop")
root.addObject("VisualStyle", displayFlags=["showVisual", "showForceFields", "showCollisionModels"])
plane_node = root.addChild("plane")
plane_node.addObject("EulerImplicitSolver")
plane_node.addObject("CGLinearSolver")
plane_node.addObject(
"MeshOBJLoader",
name="meshLoader_0",
filename="mesh/plane_loop_2.obj",
scale=0.2,
handleSeams=True,
rotation=[90, 0, 90],
translation=[0, -2.01, 0],
)
plane_node.addObject("MechanicalObject", template="Rigid3d")
plane_node.addObject("UniformMass", totalMass=10.0)
plane_node.addObject("FixedConstraint")
plane_node.addObject("UncoupledConstraintCorrection")
plane_collision_node = plane_node.addChild("collision")
plane_collision_node.addObject("MeshTopology", src="@../meshLoader_0")
plane_collision_node.addObject("MechanicalObject")
# plane_collision_node.addObject("TriangleCollisionModel")
plane_collision_node.addObject("SphereCollisionModel", radius=0.2, group=0)
plane_collision_node.addObject("UncoupledConstraintCorrection")
plane_collision_node.addObject("RigidMapping")
plane_visual_node = plane_node.addChild("visual")
plane_visual_node.addObject(
"OglModel",
name="plane",
src="@../meshLoader_0",
material="Default Diffuse 1 1 0.4 0.4 1 Ambient 1 0.8 0.8 0.8 1 Specular 0 1 1 1 1 Emissive 0 1 1 1 1 Shininess 0 45",
)
plane_visual_node.addObject("RigidMapping")
root.addObject("CollisionPipeline", depth=6, verbose=False, draw=False)
root.addObject("BruteForceBroadPhase")
root.addObject("BVHNarrowPhase")
root.addObject("MinProximityIntersection", name="Proximity", alarmDistance=0.5, contactDistance=0.3)
root.addObject("CollisionResponse", response="FrictionContactConstraint")
root.addObject("GenericConstraintSolver")
composite_node = root.addChild("composite")
composite_node.addObject(
"SparseGridMultipleTopology",
n=[6, 3, 3],
fileTopology="mesh/bubille_out.obj",
fileTopologies=["mesh/bubille_out.obj", "mesh/bubille_in1.obj", "mesh/bubille_in2.obj"],
nbVirtualFinerLevels=3,
finestConnectivity=False,
stiffnessCoefs=[1, 0.0001, 50],
massCoefs=[1, 1, 1],
)
composite_node.addObject("EulerImplicitSolver", vdamping=0, rayleighMass=0, rayleighStiffness=0)
composite_node.addObject("SparseLDLSolver")
composite_node.addObject("MechanicalObject")
composite_node.addObject(
"HexahedronCompositeFEMForceFieldAndMass",
drawType="0",
lumpedMass=False,
nbVirtualFinerLevels=2,
youngModulus=600,
poissonRatio=0.3,
method="polar",
density=0.1,
updateStiffnessMatrix=False,
printLog=False,
)
composite_node.addObject("BoxROI", box="-5 -2.1 -10 10 -1.9 10")
composite_node.addObject("FixedConstraint", indices="@BoxROI.indices")
composite_node.addObject("LinearSolverConstraintCorrection")
collision_node = composite_node.addChild("collision")
collision_node.addObject("MeshOBJLoader", name="loader", filename="mesh/bubille_out.obj")
collision_node.addObject("MeshTopology", src="@loader")
collision_node.addObject("MechanicalObject", src="@loader")
collision_node.addObject("HexahedronCompositeFEMMapping")
# collision_node.addObject("TriangleCollisionModel", group=0)
collision_node.addObject("SphereCollisionModel", group=0, radius=0.3)
visual_node = collision_node.addChild("visual")
visual_node.addObject("MeshOBJLoader", name="meshLoader_2", filename="mesh/bubille_out.obj", handleSeams=True)
visual_node.addObject("OglModel", name="VisualBody", src="@meshLoader_2", normals="0", color=[0.1, 0.8, 0.3, 0.6])
visual_node.addObject("IdentityMapping", input="@..", output="@VisualBody")
soft_bead_node = composite_node.addChild("soft bead")
soft_bead_node.addObject("MeshOBJLoader", name="meshLoader_1", filename="mesh/bubille_in1.obj", handleSeams=True)
soft_bead_node.addObject("OglModel", name="VisualBody1", src="@meshLoader_1", normals="0", color=[1, 0, 0, 1])
soft_bead_node.addObject("HexahedronCompositeFEMMapping", input="@..", output="@VisualBody1")
stiff_bead_node = composite_node.addChild("stiff bead")
stiff_bead_node.addObject("MeshOBJLoader", name="meshLoader_3", filename="mesh/bubille_in2.obj", handleSeams=True)
stiff_bead_node.addObject("OglModel", name="VisualBody2", src="@meshLoader_3", normals="0", color=[0, 0, 1, 1])
stiff_bead_node.addObject("HexahedronCompositeFEMMapping", input="@..", output="@VisualBody2")
ball_node = root.addChild("ball")
ball_node.addObject("EulerImplicitSolver", vdamping=0, rayleighMass=0, rayleighStiffness=0)
ball_node.addObject("CGLinearSolver", iterations=100, tolerance=1e-5, threshold=1e-5)
ball_node.addObject("MechanicalObject", template="Rigid3d", position=[0, 5, 0, 0, 0, 0, 1], showObject=True, showObjectScale=2.0)
ball_node.addObject("UniformMass", totalMass=10000.0)
ball_node.addObject("SphereCollisionModel", radius=0.5, group=1)
ball_node.addObject("UncoupledConstraintCorrection")
As a side question: Am I even using the right components? How would you model this scene of a liver with an embedded tumor? I also tested the Heterogeneous-TetrahedronFEMForceField.scn example, but that is even more unstable. When you interact with the object through the mouse, it applies a huge force in the opposite direction.
I also simplified the liver scene to just the SOFA liver. Same problem with the instability.
import Sofa
import Sofa.Core
PLUGINS = [
"Sofa.Component.AnimationLoop",
"Sofa.Component.Collision.Detection.Algorithm",
"Sofa.Component.Collision.Detection.Intersection",
"Sofa.Component.Collision.Response.Contact",
"Sofa.Component.Constraint.Lagrangian.Solver",
"Sofa.Component.Visual",
"Sofa.Component.Collision.Geometry",
"Sofa.Component.Constraint.Projective",
"Sofa.Component.LinearSolver.Iterative",
"Sofa.Component.Mapping.NonLinear",
"Sofa.Component.Mass",
"Sofa.Component.ODESolver.Backward",
"Sofa.Component.StateContainer",
"Sofa.GL.Component.Rendering3D",
"Sofa.Component.Constraint.Lagrangian.Correction",
"Sofa.Component.Topology.Container.Dynamic",
"MultiThreading",
"Sofa.Component.SolidMechanics.FEM.NonUniform",
"Sofa.Component.Topology.Container.Grid",
"Sofa.Component.IO.Mesh",
"Sofa.Component.LinearSolver.Direct",
"Sofa.Component.Mapping.Linear",
"Sofa.Component.Topology.Container.Constant",
]
def createScene(root: Sofa.Core.Node):
plugin_set = set(PLUGINS)
for plugin in plugin_set:
root.addObject("RequiredPlugin", name=plugin)
root.gravity = [0.0, 0.0, -9.81]
root.dt = 0.02
root.addObject("FreeMotionAnimationLoop")
root.addObject(
"VisualStyle",
displayFlags=[
"showVisual",
"showForceFields",
"showBehaviorModels",
],
)
root.addObject("CollisionPipeline", depth=6, verbose=False, draw=False)
root.addObject("BruteForceBroadPhase")
root.addObject("BVHNarrowPhase")
root.addObject("MinProximityIntersection", name="Proximity", alarmDistance=0.5, contactDistance=0.3)
root.addObject("CollisionResponse", response="FrictionContactConstraint")
root.addObject("GenericConstraintSolver")
scene_node = root.addChild("scene")
composite_node = scene_node.addChild("composite")
mesh_files = [
"mesh/liver.obj",
]
composite_node.addObject(
"SparseGridMultipleTopology",
n=[6, 6, 6],
fileTopology=mesh_files[0],
fileTopologies=mesh_files,
nbVirtualFinerLevels=2,
finestConnectivity=False,
stiffnessCoefs=[1] * len(mesh_files),
massCoefs=[1] * len(mesh_files),
)
composite_node.addObject("EulerImplicitSolver", vdamping=0, rayleighMass=0, rayleighStiffness=0)
composite_node.addObject("SparseLDLSolver")
composite_node.addObject("MechanicalObject")
composite_node.addObject(
"HexahedronCompositeFEMForceFieldAndMass",
drawType=0,
lumpedMass=False,
nbVirtualFinerLevels=2,
youngModulus=600,
poissonRatio=0.3,
method="polar",
density=0.1,
updateStiffnessMatrix=False,
)
composite_node.addObject("BoxROI", box=[-1, -5, -1, 5, 5, 5])
composite_node.addObject("FixedConstraint", indices="@BoxROI.indices")
composite_node.addObject("LinearSolverConstraintCorrection")
liver_visual = composite_node.addChild("visual")
liver_visual.addObject("MeshOBJLoader", filename=mesh_files[0])
liver_visual.addObject(
"OglModel",
src="@MeshOBJLoader",
material="Transparent Diffuse 1 1 0 1 0.45 Ambient 0 1 1 1 1 Specular 1 0 0 1 1 Emissive 0 1 0 0 1 Shininess 1 100",
)
liver_visual.addObject("HexahedronCompositeFEMMapping")
return root
This is without any interaction with the scene.
Cheers, Paul
Hi @alxbilger ,
quick update: The main reason behind the instability is the nbVirtualFinerLevels value.
If we set that to 1 -> no problem. Also for strictly convex objects, it is also no problem.
But I would still love to get deformation through collision working... Any idea on what would have to be done there?
Cheers, Paul
Hi @ScheiklP,
I'm encountering the same issue you described and was hoping to get your insights.
Here's a quick overview of what I'm trying to achieve:
Simulate a realistic human finger composed of multiple materials and meshes (I have all the necessary meshes: dermis, subcutaneous tissue, and bones).
Enable the finger to interact with an external soft material and observe the resulting deformation.
If you've managed to resolve this or have any tips on how to approach such a setup, I'd greatly appreciate your guidance.
Thanks in advance!
Hi @Paolo-Susini , Sadly I never got it to work properly. So pinging @alxbilger again. :D
Cheers, Paul