Benchmark test case: cantilever subjected to end shear force

[julien-siguenza]

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:

    import Sofa
    import Sofa.Gui

    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(rootNode):

        rootNode.addObject('RequiredPlugin', name="Shell")
        rootNode.addObject('RequiredPlugin', name="Sofa.Component.Constraint.Projective")
        rootNode.addObject('RequiredPlugin', name="Sofa.Component.MechanicalLoad")
        rootNode.addObject('RequiredPlugin', name="Sofa.Component.IO.Mesh")
        rootNode.addObject('RequiredPlugin', name="Sofa.Component.LinearSolver.Direct")
        rootNode.addObject('RequiredPlugin', name="Sofa.Component.Mapping.Linear")
        rootNode.addObject('RequiredPlugin', name="Sofa.Component.Mass")
        rootNode.addObject('RequiredPlugin', name="Sofa.Component.ODESolver.Backward")
        rootNode.addObject('RequiredPlugin', name="Sofa.Component.StateContainer")
        rootNode.addObject('RequiredPlugin', name="Sofa.Component.Topology.Container.Constant")
        rootNode.addObject('RequiredPlugin', name="Sofa.Component.Visual")
        rootNode.addObject('RequiredPlugin', name="Sofa.GL.Component.Rendering3D")
        rootNode.addObject('RequiredPlugin', name="Sofa.GUI.Component")

        rootNode.dt = 0.01
        rootNode.gravity = [0, 0, 0]
        rootNode.addObject('DefaultAnimationLoop')

        rootNode.addObject('VisualStyle',
            displayFlags=[
                "hideVisualModels",
                "hideBehaviorModels",
                "showMappings",
                "showForceFields"
            ]
        )

        plate = rootNode.addChild('plate')
        plate.addObject('EulerImplicitSolver',
            rayleighStiffness=0.0,
            rayleighMass=5.0,
        )
        plate.addObject('SparseLDLSolver',
            template="CompressedRowSparseMatrixMat3x3d",
        )
        plate.addObject('MeshGmshLoader',
            name="loader",
            filename="mesh/plate.msh",
        )
        plate.addObject('MeshTopology',
            name="topology",
            src="@loader",
        )
        plate.addObject('MechanicalObject',
            template="Rigid3",
        )
        plate.addObject('UniformMass',
            totalMass=0.005,
        )
        plate.addObject('FixedConstraint',
            indices=[0, 3, 85, 86, 87],
            # indices=[0, 3, 85, 86, 87, 4, 84, 125, 202, 236, 238],
        )
        plate.addObject('ConstantForceField',
            indices=[1, 2, 43, 44, 45],
            totalForce=[0, 0, 0.04],
            showArrowSize=100,
            printLog=True,
        )
        plate.addObject('TriangularBendingFEMForceField',
            youngModulus=1.2e6,
            poissonRatio=0.0,
            thickness=0.1,
        )

        visu = plate.addChild('Visu')
        visu.addObject('OglModel',
            src="@../topology",
        )
        visu.addObject('IdentityMapping')

        return

    if __name__ == '__main__':
        
        main()

I give the following mesh as input to the MeshGmshLoader, that defines a rectangular plate with the dimensions (L, b):

I fix the left extremity of the plate using a FixedConstraint and load the right extremity using the ConstantForceField, which I first set with a tiny totalForce of 0.04 (1% of the maximum force Pmax) in the direction tangential to the plate.

I define a big amount of rayleighMass (5.0) when setting the EulerImplicitSolver, in order to reach more quickly the steady deformation of the plate.

The TriangularBendingFEMForceField is finally used to define the shell, setting the Young's modulus to 1.2e6, the Poisson ratio to 0.0 and the thickness (dimension h) to 0.1.

For the considered force, I am thus expecting a moderate deformation of the plate.
Instead, I am ending up with a result that I did not expect at all:

I suspect that this unexpected behavior is due to the FixedConstraint which probably constrains only the translational DOFs of the shell, but not the rotational DOFs of the shell on the selected nodes.
Could anyone please confirm this and tell me how I can constrain the rotational DOFs?

A dirty workaround I here found to constrain the rotational DOFs on the desired extremal nodes was to constrain also the neighboring nodes:

I finally end up with the following result:

The behavior of the shell is much softer than what I am expecting looking at the reference result of this test case.
Another worrying thing is that increasing either the Young's modulus or the thickness of the shell both lead to a very similar result...

Can anyone please help me see more clearly in all these inconsistencies?
I guess I defined something wrong in my scene, but cannot find out what.

If it can be of any help, here are the messages that are displayed in my terminal when I am running my scene:

Thanks for helping 🙏

FYI: @hugtalbot @bakpaul