HexahedronCompositeFEMForceFieldAndMass Solver

[ScheiklP]

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.

FEM-2024-03-18_18.24.42.mp4

SOFA: v23.12

Cheers,
Paul

[alxbilger]

Hi,

The message is located here:

sofa/Sofa/Component/SolidMechanics/FEM/NonUniform/src/sofa/component/solidmechanics/fem/nonuniform/HexahedronCompositeFEMMapping.h

Line 105 in b5aefb4

msg_warning() << "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.";

I don't think that the message is true. @epernod Do you remember why you added this message in 0170a07#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

[ScheiklP]

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.

FEM-2024-03-19_10.41.08.mp4

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

FEM-2024-03-19_10.54.23.mp4

This is without any interaction with the scene.

Cheers,
Paul

[ScheiklP]

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