Open ScheiklP opened 1 month ago
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:
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
Problem
Hi, the example for the HexahedronCompositeFEMForceFieldAndMass complains, that
Even when using the
SparseLDLSolver
instead of theCGLinearSolver
.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