Simple stretch test

[ateixeira0163]

Hi there! I was recently interested by the behavior of a simple stretch simulation using different solvers. Interestingly enough, I found some weird behaviors while using some solvers and I could not understand why. To exemplify, I used a scene described below:

root.gravity = [0, 0, 0]
root.dt = .1
root.addObject('DefaultAnimationLoop') [or] root.addObject('FreeMotionAnimationLoop')

beam = root.addChild('beam')
beam.addObject('EulerImplicitSolver', solveConstraint=True)
# --- LinearSolvers --- # 
beam.addObject('CG' [or] 'SparseLDLSolver' [or] 'SparseCholeskySolver' [or] 'ShewchukPCGLinearSolver' + 'SparseLDL')
# --- Rest --- # 
beam.addObject('RegularGridTopology', name="beam_grid_topology", n=[5, 5, round(L / (D / 4))],
               min=[-D / 2, -D / 2, 0], max=[D / 2, D / 2, L])
beam.addObject('HexahedronSetTopologyContainer', name='topology', src='@beam_grid_topology')
beam.addObject('HexahedronSetGeometryAlgorithms')
beam.addObject('MechanicalObject', template='Vec3d', name='mech_obj', src='@beam_grid_topology')
beam.addObject('MeshMatrixMass', massDensity=1160, topology='@topology')
beam.addObject('HexahedronFEMForceField', name='ff', src='@topology', poissonRatio=.45,
                            youngModulus=1.25e6)
beam.addObject('BoxROI', name='box_fix', box=[-D, -D, -e, D, D, e], template='Vec3d')
beam.addObject('FixedConstraint', indices='@box_fix.indices')
beam.addObject('BoxROI', name='box_end', box=[-D, -D, L - e, D, D, L + e], template='Vec3d')
beam.addObject('LinearMovementConstraint', indices='@box_end.indices',
                                            keyTimes=[0, 100],
                                            movements=[0, 0, 0,
                                                       0, 0, .1],
                                            relativeMovements=True, showMovement=True)
beam.addObject('LinearSolverConstraintCorrection', solverName='linear_solver')

So the idea being just a beam being stretched by a LinearMovementConstraint while fixed on the other end using Hexahedron elements.
But, the thing is, for example, using:
A. DefaultAnimationLoop + EulerImplicitSolver + CG + HexahedronFF + MeshMatrixMass = I get a plausible result as expected.

But if I change the linear solver for a direct solver like SparseLDLSolver for example, the beam starts to oscillate quite quickly and it only gets worse (img down below). I've tried others direct solvers SparseLU, SparseCholesky but with no success (always the same result) even when trying to reduce the tolerance. Even the ShewchukPCG that is also an iterative solver (but uses a direct solver to precondition) didn't manage to keep the beam from oscillating.

I also did some tests with Caribou, but I sent the questions to Jean Nicolas.

After that, I tried to do the same tests but using FreeMotionAnimationLoop (this is kinda important for my simulation goal at the end). The first result using the combo EulerImplicitSolver + CG + HexaFF + MeshMatrixMass result in the image below:

Where it seems that the deformation is not "communicated" with the rest of the body. As I read, the FreeMotionAnimationLoop needs a ConstraintCorrection object. For that I need to use LinearSolverConstraintCorrection (the other ones are too approximated for my end case). The problem with that being that the direct solvers doesn't solve correctly the problem but they work with the ConstraintCorrection and CG does not work with LinearSolverConstraintCorrection. I was trying to evaluate the condition number of the stiffness matrix but I only found the matrix by elements, is there a way to assemble and get all the stiffness/mass matrix of the system before calling the linear solver ?

So I would like to know if there's any tests made for the stretching and any thoughts on this.
Thanks! (I hope it was not very confusing, because I, myself, am confused with these results)
Alexandre T

[alxbilger]

Hi @ateixeira0163,

in your stretch test, can you test if there is a difference between
node.addObject('CGLinearSolver', template='CompressedRowSparseMatrixd')
and
node.addObject('CGLinearSolver', template='CompressedRowSparseMatrixMat3x3d')
?

Regarding the assembly of the stiffness/mass matrices, it is still a work in progress. Some works are available in pull request (sofa-framework/SofaPython3#262 and sofa-framework/SofaPython3#260).

FreeMotionAnimationLoop requires a constraint solver AND constraint corrections.

[alxbilger]

With SparseLDLSolver could you test disabling the permutations by setting applyPermutation to false?

[ateixeira0163]

I did not found an applyPermutation in the SparseLDLSolver class.

[alxbilger]

You are probably not up-to-date (

sofa/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLDLSolverImpl.h

Line 153 in feb6926

, d_applyPermutation(initData(&d_applyPermutation, true ,"applyPermutation", "If true the solver will apply a fill-reducing permutation to the matrix of the system."))

)

[ateixeira0163]

I have the 21.12 version. But indeed this line does not exist in the source code that I have.
I'll do a pull and recompile to test it out.

[ateixeira0163]

Indeed I was not on the master branch. Some compiling and debugging time after, I've got Sofa working and I did the test with applyPermutation=False but the result is still the same.

[hugtalbot]

Hi @ateixeira0163

(edit)
I tested your scene with :

• L = 10
• D = 1
• e = 0.2

and using the following cases:

• DefaultAnimationLoop
  • CGLinearSolver (without LinearSolverConstraintCorrection since no Lagrangian constraint are computed with DefaultAL)
  • SparseLUSolver (without LinearSolverConstraintCorrection since no Lagrangian constraint are computed with DefaultAL)
• FreeMotionAnimationLoop (with GenericConstraintCorrection)
  • CGLinearSolver + LinearSolverConstraintCorrection: it works only because there is actually no Lagrangian constraint in the scene. If there is any constraint, the CGLinearSolver must have a template different than GraphScattered which corresponds to unbuilt computation (use e.g. template="CompressedRowSparseMatrixMat3x3d" instead of GraphScattered). In other words, if there is any constraint, the CGLinearSolver must build the A matrix.
  • SparseLDLSolver + LinearSolverConstraintCorrection

everything works here.
Note that:

• with the FreeMotionAnimationLoop, only one direct solver seems to be permitted
• I changed slight the LinearMovementConstraint in my scene, here is what I wrote:

    beam.addObject('LinearMovementConstraint', indices='@box_end.indices',
                                                keyTimes=[0, 100, 999],
                                                movements=[0, 0, 0,
                                                           0, 0, .1,
                                                           0, 0, .1],
                                                relativeMovements=True, showMovement=True)

Let me know if this helps

[alxbilger]

This scene works without any oscillation and without any crash.
Let me know if this scene is not a good example of your problem. Please provide a full scene if you want me to try a new scene.
I am on SOFA master.

[ateixeira0163]

Hi!
I've managed to make this work before, but the main difference being the scaling.
If I change the physical parameters on the scene you sent to:
L = 3.5 and D = .083, in the 18,6s time mark, the beam will start to oscillate.
I understand that the system is quite unstable, but since there's no perturbation, it shouldn't be a problem.
Like if use the TetrahedronElasticForce from Caribou, the system does not oscillate.

[alxbilger]

Ok, with the parameters L and D that you provided I also observe the oscillations

[alxbilger]

I am not totally convinced the issue comes from the solver. It may come from the model, but the CG attenuates the oscillations.
To verify that, you can try to change the model. For example, you could try to have a linear geometry (method small in HexahedronFEMForceField). That would exclude the computation of the rotations. You could also switch to a tetra topology to try other components: FastTetrahedralCorotationalForceField, TetrahedralCorotationalFEMForceField, or TetrahedronHyperelasticityFEMForceField.

In Caribou, how many Newton iterations did you perform?

[ateixeira0163]

I don't think either.
When I said I tried with Caribou, it was only the ForceField that I changed (HexahedronElasticForce for elastic model). Also, using the HexahedronFEMForceFieldAndMass yields a better result, with oscillations that only appears if the dt is too big. It seems that the problem is in the building of the Mass Matrix in some sort.

[hugtalbot]

We would need additional eyes on this topic. Instabilities (vibrations/oscillations) are noticed when a bar is stretched. Here are the current conclusions of the investigation:

• iterative CG solver results in a stable stretching until it's released
• while a direct solver leads to instabilities

@ateixeira0163 My only interpretation is that the time step is way to large regarding the mesh size (about dx=0.001 while `dt=0.1") which IMO contradicts the statibility principle according which a limit time step cannot be exceeded depending on the mesh size.
Apparently you succeeded to get it stable in another software (with same time step and mesh), do you confirm?

A 🍺 for the problem solver ! @alxbilger @jnbrunet @ChristianDuriez @courtecuisse @StephaneCotin !

A basic scene to test (click to expand)

(E = 1.25e6 ν = .45 ρ = 1160 dt = .1 gravity = 0 displacement = 0.001/step) :

# Required import for python
import Sofa
import numpy as np

L = 0.35
D = 0.0083
e = 1e-5

def main():
    import SofaRuntime
    import Sofa.Gui

    root = Sofa.Core.Node("root")
    createScene(root)
    Sofa.Simulation.init(root)

    Sofa.Gui.GUIManager.Init("myscene", "qglviewer")
    Sofa.Gui.GUIManager.createGUI(root, __file__)
    Sofa.Gui.GUIManager.SetDimension(1080, 600)
    Sofa.Gui.GUIManager.MainLoop(root)
    Sofa.Gui.GUIManager.closeGUI()


def createScene(root):
    root.gravity = [0, 0, 0]
    root.dt = .1
    # root.addObject('DefaultAnimationLoop')
    root.addObject('FreeMotionAnimationLoop')
    root.addObject('GenericConstraintSolver', maxIterations=1000, tolerance=1e-6)
    root.addObject('DefaultVisualManagerLoop')
    root.addObject('VisualStyle', displayFlags="showForceFields showBehaviorModels")

    beam = root.addChild('beam')
    beam.addObject('EulerImplicitSolver', solveConstraint=True)
    # --- LinearSolvers --- # 
    # beam.addObject('CGLinearSolver', template="CompressedRowSparseMatrixMat3x3d", name='linear_solver_correction', iterations=1000, tolerance=1e-5, threshold=1e-5)
    beam.addObject('SparseLDLSolver', name='linear_solver_correction', template="CompressedRowSparseMatrixMat3x3d")
    # --- Rest --- # 
    beam.addObject('RegularGridTopology', name="beam_grid_topology", n=[5, 5, round(L / (D / 4))],
                   min=[-D / 2, -D / 2, 0], max=[D / 2, D / 2, L])
    beam.addObject('HexahedronSetTopologyContainer', name='topology', src='@beam_grid_topology')
    beam.addObject('HexahedronSetGeometryAlgorithms')
    beam.addObject('MechanicalObject', template='Vec3d', name='mech_obj', src='@beam_grid_topology', showObject=True, showObjectScale=2)
    beam.addObject('MeshMatrixMass', massDensity=1160, topology='@topology', geometryState='@mech_obj')
    beam.addObject('HexahedronFEMForceField', name='ff', src='@topology', poissonRatio=.45,
                                youngModulus=1.25e6)
    beam.addObject('BoxROI', name='box_fix', box=[-D, -D, -e, D, D, e], template='Vec3d', drawBoxes=True)
    beam.addObject('FixedConstraint', indices='@box_fix.indices')
    beam.addObject('BoxROI', name='box_end', box=[-D, -D, L - e, D, D, L + e], template='Vec3d')
    beam.addObject('LinearMovementConstraint', indices='@box_end.indices',
                                                keyTimes=[0, 100, 999],
                                                movements=[0, 0, 0,
                                                           0, 0, .1,
                                                           0, 0, .1],
                                                relativeMovements=True, showMovement=True)
    beam.addObject('LinearSolverConstraintCorrection', solverName='linear_solver_correction')

    return root



# Function used only if this script is called from a python environment
if __name__ == '__main__':
    main()

[hugtalbot]

and I add @Sidaty1 who has a recent experience to the topic!

[github-actions[bot]]

⚠️ ⚠️ ⚠️
@ateixeira0163
We are setting a new forum management policy: topics like this one, in which the last reply is more than 4 months old, will automatically be closed. Without further news in the coming weeks, our bot will proceed to automatic closure, thus keeping this forum clean and fresh 🌱 Thank you for your understanding

[github-actions[bot]]

⚠️ ⚠️ ⚠️
@ateixeira0163
In accordance with our forum management policy, the last reply is more than 4 months old and is therefore closed. Our objective is to keep the forum up to date and offer the best support experience.

Please feel free to reopen it if the topic is still active by providing us with an update. Please feel free to open a new thread at any time - we'll be happy to help where and when we can.

[github-actions[bot]]

⚠️ ⚠️ ⚠️
@ateixeira0163
In accordance with our forum management policy, the last reply is more than 4 months old and is therefore closed. Our objective is to keep the forum up to date and offer the best support experience.

Please feel free to reopen it if the topic is still active by providing us with an update. Please feel free to open a new thread at any time - we'll be happy to help where and when we can.