Drift in beam deformation example with tetrahedral mesh and SurfacePressureForceField

[th-skam]

While working on validation examples for a material model, I encountered a behaviour I was not expecting.

The setup involves a beam aligned with the x-axis that gets deflected in the y-direction due to a load applied to its bottom surface. The load is applied using the SurfacePressureForceField.

The figures show a comparison of two deformed beams. They are simulated using the same parameters, only the meshes are different: tetrahedral (blue) and hexahedral (green/red)

I see 2 problems here but this issue is mainly about the second one.

1. Apparent "extra stiffness" in the tetrahedral beam. We have discussed in the past about locking with tetrahedral meshes and, as I am using a near-incompressible model here ($\nu = 0.495$), this was not a big surprise.

2. The tetrahedral beam gets deflected laterally. At first, I thought this was an isolated problem caused solely by the interpolation of the SurfacePressureForceField from the triangles to the vertices. However, by reducing Poisson's ratio, the deflection decreases, so these two issues could be related.

What do you think about this? Could this be caused by a malfunctioning SurfacePressureForceField or do you think it's a problem inherent in the linear tetrahedral FEM?

The scene file for reference

import Sofa

def createScene(root):
    root.gravity=[0, 0, 0]
    root.dt = 1e-3
    root.addObject("RequiredPlugin",pluginName=['Sofa.Component.Constraint.Projective',
                                                'Sofa.Component.Engine.Select',
                                                'Sofa.Component.LinearSolver.Direct',
                                                'Sofa.Component.Mapping.Linear',
                                                'Sofa.Component.Mass',
                                                'Sofa.Component.MechanicalLoad',
                                                'Sofa.Component.ODESolver.Backward',
                                                'Sofa.Component.SolidMechanics.FEM.Elastic',
                                                'Sofa.Component.StateContainer',
                                                'Sofa.Component.Topology.Container.Dynamic',
                                                'Sofa.Component.Topology.Container.Grid',
                                                'Sofa.Component.Topology.Mapping',
                                                'Sofa.Component.Visual',
                                                'Sofa.GL.Component.Rendering3D'
                                                ])
  
    root.addObject("DefaultAnimationLoop")
    createBeam(root, femType="Hex", zOffset=0)
    createBeam(root, femType="Tet", zOffset=1.5)
    
def createBeam(root, femType, zOffset):

    # Beam Geometry (in mm)
    nPtsLat = 4
    nPtsLength = 10*nPtsLat - 9
    grid = root.addChild("Grid"+femType)
    grid.addObject("RegularGridTopology", name="gridTopo", 
        min= [0 , 0, 0 + zOffset], 
        max= [10, 1, 1 + zOffset], 
        n  = [nPtsLength, nPtsLat, nPtsLat]
    )

    # Beam Model
    beam = root.addChild("Beam"+femType)
    beam.addObject("VisualStyle", displayFlags=["showForceFields"])

    # Solver Setup
    beam.addObject("NewtonRaphsonSolver", name="newtonSolver", warnWhenLineSearchFails=False,
                   maxNbIterationsNewton=200, maxNbIterationsLineSearch=5,
                   absoluteResidualStoppingThreshold  = 1e-12,
                   relativeInitialStoppingThreshold   = 1e-12,
                   absoluteEstimateDifferenceThreshold= 1e-12,
                   relativeEstimateDifferenceThreshold= 1e-12
    )
    beam.addObject("StaticSolver", name="staticSolver", newtonSolver="@newtonSolver")
    beam.addObject("SparseLDLSolver", name="linearSolver", template="CompressedRowSparseMatrixMat3x3d") 

    # Mesh & FEM Setup
    match femType:
        case "Hex":
            beam.addObject("HexahedronSetTopologyContainer", name="hexContainer", 
                           position="@Grid"+femType+"/gridTopo.position", 
                           hexahedra="@Grid"+femType+"/gridTopo.hexahedra")
            beam.addObject("HexahedronSetTopologyModifier", name="hexModifier")
            beam.addObject("HexahedronSetGeometryAlgorithms", name="hexGeomAlgo")
            beam.addObject("MechanicalObject", name="mechObj", template="Vec3d")
            beam.addObject("DiagonalMass", massDensity=1e-6)
            beam.addObject("HexahedronFEMForceField", name="FEM", 
                           youngModulus=0.008, 
                           poissonRatio=0.495, 
                           method="large"
            )
        case "Tet":
            beam.addObject("TetrahedronSetTopologyContainer", name="tetContainer", 
                           position="@Grid"+femType+"/gridTopo.position")
            beam.addObject("TetrahedronSetTopologyModifier", name="tetModifier")
            beam.addObject("TetrahedronSetGeometryAlgorithms", name="tetGeomAlgo")
            beam.addObject("Hexa2TetraTopologicalMapping" , name="hexa2tetMap", 
                           input="@Grid"+femType+"/gridTopo", output="@tetContainer")
            beam.addObject("MechanicalObject", name="mechObj", template="Vec3d")
            beam.addObject("DiagonalMass", massDensity=1e-6)
            beam.addObject("TetrahedronFEMForceField", name="FEM", 
                           youngModulus=0.008, 
                           poissonRatio=0.495, 
                           method="large"
            )
        case _:
            return "Invalid femType"

    # Boundary Conditions - Fixed End
    beam.addObject("BoxROI", name="fixedROI", drawBoxes=False, doUpdate=False,
                   box=[-0.01, -0.01, -0.01+zOffset, 0.01, 1.01, 1.01+zOffset])
    beam.addObject("FixedProjectiveConstraint", name="fixedConstraint", indices="@fixedROI.indices")

    # Surface Load - Pressure on Bottom Faces
    surf = beam.addChild("Surface")
    surf.addObject("VisualStyle", displayFlags=["showForceFields"])

    match femType:
        case "Hex":
            surf.addObject("QuadSetTopologyContainer", name="quadContainer", 
                           position="@../hexContainer.position")
            surf.addObject("Hexa2QuadTopologicalMapping", name="hex2quadMap", 
                           input="@../hexContainer", output="@quadContainer")
            surf.addObject("MechanicalObject", name="surfMechObj", template="Vec3d")
            surf.addObject("BoxROI", name="loadROI", drawBoxes=False, doUpdate=False,
                           box=[0, 0, 0+zOffset,  10, 0.2, 1+zOffset])
            surf.addObject("SurfacePressureForceField", name="load", 
                           quadIndices="@loadROI.quadIndices", 
                           pressure = -4e-6,
                           drawForceScale=1e7
            )
        case "Tet":
            surf.addObject("TriangleSetTopologyContainer", name="triContainer", 
                           position="@../tetContainer.position")
            surf.addObject("Tetra2TriangleTopologicalMapping", name="tet2triMap", 
                           input="@../tetContainer", output="@triContainer")
            surf.addObject("MechanicalObject", name="surfMechObj", template="Vec3d")
            surf.addObject("BoxROI", name="loadROI", drawBoxes=False, doUpdate=False,
                           box=[0, 0, 0+zOffset,  10, 0.2, 1+zOffset])
            surf.addObject("SurfacePressureForceField", name="load", 
                           triangleIndices="@loadROI.triangleIndices", 
                           pressure = -4e-6,
                           drawForceScale=1e7
            )
        case _:
            return "Invalid femType"

    surf.addObject("IdentityMapping", name="idMap", input="@../mechObj", output="@surfMechObj")

[hugtalbot]

Hi @th-skam

Thanks for the issue. The drift that you notice for your tetrahedral mesh comes from the fact that you mesh is "oriented". As you can notice all tetrahedra have an edge in the following direction: ↘ diffusing errors/forces in an anisotropic way.

I would therefore recommend to remesh your beam.
For instance, you could use edge swapping.

EDIT:
you can use in the Hexa2TetraTopologicalMapping, the data swapping=True

[th-skam]

Hello @hugtalbot

Yes, swapping the edges solved the drift:

Thanks for the assist!