Benchmark test case: cantilever subjected to end shear force (volumetric FEM)

[julien-siguenza]

NB: A similar topic has already been opened here considering the Shell FEM Elastic plugin. This is a brand new topic where the same test case is studied, but considering the original volumetric FEM modeling of SOFA Framework.

---

Hi there,

I create this topic to implement a benchmark test case of cantilever subjected to end shear force, as presented in the section 3.1 of the publication of Sze et al. (2004):

I thus implemented the following python scene, which attempts to reproduce this test case within the SOFA environment:

    import Sofa
    import Sofa.Gui

    import os
    import sys
    import numpy as np
    import pandas as pd
    import pyvista as pv
    from pathlib import Path

    def main():

        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, 1080)
        Sofa.Gui.GUIManager.MainLoop(root)
        Sofa.Gui.GUIManager.closeGUI()

    def createScene(root):

        path_mesh = "plate.msh"
        path_model = "plate.vtk"
        young = 1.2e6
        poisson = 0.0

        body = pv.read(path_model)

        mask = body.point_data["Constraint"] == 1
        ids_constraint = np.where(mask == True)[0].tolist()

        mask = body.point_data["Load"] == 1
        ids_load = np.where(mask == True)[0].tolist()

        root.gravity=[0.0, 0.0, 0.0]
        root.dt=0.05
        root.bbox="-20 -20 -20 20 20 20"

        root.addObject('RequiredPlugin', name="Sofa.Component.AnimationLoop")
        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.Lagrangian.Correction")
        root.addObject('RequiredPlugin', name="Sofa.Component.Constraint.Lagrangian.Solver")
        root.addObject('RequiredPlugin', name="Sofa.Component.Constraint.Projective")
        root.addObject('RequiredPlugin', name="Sofa.Component.IO.Mesh")
        root.addObject('RequiredPlugin', name="Sofa.Component.LinearSolver.Direct")
        root.addObject('RequiredPlugin', name="Sofa.Component.ODESolver.Backward")
        root.addObject('RequiredPlugin', name="Sofa.Component.SolidMechanics.FEM.Elastic")
        root.addObject('RequiredPlugin', name="Sofa.Component.MechanicalLoad")
        root.addObject('RequiredPlugin', name="Sofa.Component.StateContainer")
        root.addObject('RequiredPlugin', name="Sofa.Component.Topology.Container.Constant")
        root.addObject('RequiredPlugin', name="Sofa.Component.Topology.Container.Dynamic")
        root.addObject('RequiredPlugin', name="Sofa.Component.Visual")
        root.addObject('RequiredPlugin', name="Sofa.Component.ODESolver.Forward")
        root.addObject('RequiredPlugin', name="Sofa.GL.Component.Rendering3D")
        root.addObject('RequiredPlugin', name="Sofa.Component.LinearSolver.Preconditioner")
        
        root.addObject('VisualStyle',
            displayFlags=[
                "hideVisualModels",
                "hideBehaviorModels",
                "showMappings",
                "showForceFields"
            ]
        )
        
        root.addObject('DefaultAnimationLoop')

        root.addObject('EulerImplicitSolver',
            name="ODE solver",
            rayleighStiffness=0.0,
            rayleighMass=0.0,
        )
        # root.addObject('StaticSolver',
        #     name="ODE solver",
        # )
        
        root.addObject('SparseLDLSolver',
            name="Linear solver",
            template="CompressedRowSparseMatrixMat3x3d",
        )
        # root.addObject('CGLinearSolver',
        #     name="Linear solver",
        #     iterations=5000,
        #     tolerance=1e-20,
        #     threshold=1e-20,
        # )
        
        plate = root.addChild('plate')
        plate.addObject('MeshTopology',
            name="mesh",
            filename=str(path_mesh),
        )
        plate.addObject('MeshGmshLoader',
            name="loader",
            filename=str(path_mesh),
        )
        plate.addObject('MechanicalObject',
            template="Vec3",
            name="dofs",
            src="@loader",
        )
        plate.addObject('HexahedronSetTopologyContainer',
            src="@loader",
            name="topo",
        )
        plate.addObject('HexahedronFEMForceField',
            template="Vec3",
            name="FEM",
            youngModulus=young,
            poissonRatio=poisson,
            method="large",
        )
        # plate.addObject('DiagonalMass',
        #    name="mass",
        #    massDensity=1.0e-20,
        #    topology="@topo",
        # )
        plate.addObject('FixedConstraint',
            name="FixedConstraint",
            indices=ids_constraint,
        )
        plate.addObject('ConstantForceField',
            name="Load",
            indices=ids_load,
            totalForce=[0.0, 0.0, 0.0],
            showArrowSize=10.0,
        )
        plate.addObject(TimeDependentLoad(node=plate))

        return root


    class TimeDependentLoad(Sofa.Core.Controller):

        def __init__(self, *args, **kwargs):
            
            Sofa.Core.Controller.__init__(self, *args, **kwargs)
            self.node = kwargs.get('node')
            self.load = self.node.getObject('Load')

        def onAnimateBeginEvent(self, event):
            
            current_time = self.node.time.value
            self.updateLoad(current_time)

        def updateLoad(self, current_time):

            ramp = 1.0
            if current_time <= ramp:
                self.load.totalForce = [0.0, 0.0, 4.0*current_time/ramp]
            else:
                self.load.totalForce = [0.0, 0.0, 4.0]


    if __name__ == '__main__':
        main()

Two different meshes (both composed of hexahedrons) are here considered to define the rectangular plate with the dimensions (L, b, h):

• Mesh1, having a spacial discretization dx = 0.1 (at the left on the image below).
• Mesh2, having a spacial discretization dx = 0.05 (at the right on the image below).

For each mesh, the group of nodes located at the plate extremities are flagged with specific data fields, which are then recovered in the python scene above to later apply the boundary conditions in the simulation:

The maximum force Pmax is applied with a temporal ramp of 1.0, using the TimeDependentLoad class. Given that the timestep is set to 0.05, this means that the force Pmax is decomposed into 20 incremental loads.

The ODE solver is here initially set to EulerImplicitSolver (backward), and the linear solver is set to SparseLDLSolver (direct).

The HexahedronFEMForceField is finally used to define the mechanical behavior of the plate, setting the Young's modulus to 1.2e6, and the Poisson ratio to 0.0.

Results

In the first place, I am only trying to reach the equilibrium state of the plate when the maximum force Pmax is applied:

• I am able to achieve this objective for the Mesh1:

Enregistrement.2024-08-19.170525.mp4

• However, the Mesh2 is giving some diverging oscillations and the equilibrium state is never reached:

Enregistrement.2024-08-19.171054_compressed.mp4

---

FYI: @hugtalbot @bakpaul

[bakpaul]

Hi ! thank you for this work, this really help to get some very important insights on the current state of the Framework and will help for regression testing in the future !

I am actually wondering if this oscillation comes from the projective constraint you are using at the base ?

As you might know, projective constraints are weak, in the way that they are not explicitly expressed in the equation but more corrected a posteriori, leading to the fact that you'll get the effect of such projection on the internal forces in the next step. This might result in an oscillating behavior as you are experiencing, and is magnified here by the fact that your beam has way more element along the beam than the first one. This leads to getting less effect on the high deformation of the basis on the last node, thus increasing the delay between its position correction and the base correction.

I can think of three ways of checking if this is effectively the case :

1. Try to keep the same number of element along the beam than the more coarse one while keeping the same number along the other dimensions, to see if we go back to the other behavior
2. Try to add Rayleigh damping, if what I expect is happening, the oscillations should be modified, but I cannot predict the changes.
3. Replace the projective constraints by either springs or langrangian based constraints.

[julien-siguenza]

Hi Paul!

Thank you for your feedback 🙏

I followed your third advice and tried to replace the FixedConstraint by a RestShapeSpringsForceField:

plate.addObject('RestShapeSpringsForceField',
    stiffness=1e6,
    angularStiffness=1e6,
    points=ids_constraint,
)

Unfortunately, I end up with exactly the same behavior as described in my previous message. So I guess the issue does not come from the projective constraint.

Any other idea?

Thanks!