<!-- GENERATED by tools/generators/examplesDocsEmitter.py from python/ - do not edit -->
(testmodels-flexiblependulumbeamcomparison)=
# flexiblePendulumBeamComparison.py

You can view and download this file on Github: [flexiblePendulumBeamComparison.py](https://github.com/jgerstmayr/EXUDYN/tree/master/python/TestModels/flexiblePendulumBeamComparison.py)

```{code-block} python
:linenos:

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  A highly flexible pendulum, released from horizontal under gravity, modelled three times
#           with 10 elements: ObjectBeamGeometricallyExact2D, ObjectBeamGeometricallyExact (3D, moving in
#           the x-y plane) and ObjectANCFCable2D (#2730). Gravity acts as nodal forces, half the weight
#           of an element at each of its nodes, the same for all three. Measured: with 10 elements the
#           three differ from a 40-element GE2D solution by 9 mm (GE2D), 11 mm (GE3D) and 5 mm (ANCF2D)
#           in the tip position; with 40 elements GE3D agrees with GE2D to 0.1 mm and ANCF2D to 0.7 mm.
#           The result is the sum of the tip positions at the end, so that a change of any of the three
#           moves it.
#
# Author:   Johannes Gerstmayr
# Date:     2026-09-29
#
# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

import exudyn as exu
from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities

import numpy as np

testIsActive = exu.sys.get('testIsActive', False)

L = 0.5                 #length of the pendulum
nElements = 10
E = 1e8                 #a very soft elastomer
rho = 1000
h = 0.002               #height of the cross section (in y)
b = 0.01                #width (in z)
A = b*h
Izz = b*h**3/12         #bending in the x-y plane
Iyy = h*b**3/12
nu = 0.3
ks = 10*(1+nu)/(12+11*nu)   #shear correction
G = E/(2*(1+nu))
g = 9.81
lElement = L/nElements
tEnd = 1
stepSize = 0.0025


def Pendulum(kind):
    """the tip position over time"""
    SC = exu.SystemContainer()
    mbs = SC.AddSystem()
    nodes = []
    for i in range(nElements+1):
        x = i*lElement
        if kind == 'GE2D':
            nodes.append(mbs.AddNode(NodeRigidBody2D(referenceCoordinates=[x, 0, 0])))
        elif kind == 'GE3D':
            nodes.append(mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[x, 0, 0] + eulerParameters0)))
        else:
            nodes.append(mbs.AddNode(NodePoint2DSlope1(referenceCoordinates=[x, 0, 1, 0])))
    for i in range(nElements):
        if kind == 'GE2D':
            mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers=[nodes[i], nodes[i+1]], length=lElement,
                          massPerLength=rho*A, crossSectionInertia=rho*Izz, bendingStiffness=E*Izz,
                          axialStiffness=E*A, shearStiffness=ks*G*A))
        elif kind == 'GE3D':
            section = exu.BeamSection()
            section.stiffnessMatrix = np.diag([E*A, ks*G*A, ks*G*A, G*(Iyy+Izz), E*Iyy, E*Izz])
            section.inertia = np.diag([rho*(Iyy+Izz), rho*Iyy, rho*Izz])
            section.massPerLength = rho*A
            mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[nodes[i], nodes[i+1]], length=lElement,
                                                       sectionData=section))
        else:
            mbs.AddObject(ObjectANCFCable2D(nodeNumbers=[nodes[i], nodes[i+1]], length=lElement,
                          massPerLength=rho*A, bendingStiffness=E*Izz, axialStiffness=E*A))
    #gravity as nodal forces: half the weight of an element at each of its nodes
    for i in range(nElements+1):
        weight = rho*A*lElement*(0.5 if i in [0, nElements] else 1)*g
        mbs.AddLoad(LoadForceVector(markerNumber=mbs.AddMarker(MarkerNodePosition(nodeNumber=nodes[i])),
                                    loadVector=[0, -weight, 0]))
    #the hinge at the first node
    mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=mbs.AddNode(NodePointGround()), coordinate=0))
    for i in range(3 if kind == 'GE3D' else 2):
        mbs.AddObject(ObjectConnectorCoordinate(markerNumbers=[mGround,
                      mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nodes[0], coordinate=i))]))
    sTip = mbs.AddSensor(SensorNode(nodeNumber=nodes[-1], outputVariableType=exu.OutputVariableType.Position,
                                    storeInternal=True, writeToFile=False))
    mbs.Assemble()
    simulationSettings = exu.SimulationSettings()
    simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
    simulationSettings.timeIntegration.endTime = tEnd
    simulationSettings.timeIntegration.verboseMode = 0
    simulationSettings.solution.file.write = False
    simulationSettings.solution.sensors.writePeriod = stepSize
    simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
    mbs.SolveDynamic(simulationSettings)
    return mbs.GetSensorStoredData(sTip)


tips = dict((kind, Pendulum(kind)) for kind in ['GE2D', 'GE3D', 'ANCF2D'])
for kind in ['GE3D', 'ANCF2D']:
    difference = np.abs(tips[kind][:, 1:3] - tips['GE2D'][:, 1:3]).max(axis=1)
    exu.Print(kind + ' against GE2D: max tip difference', difference.max(), 'until t=0.5:', difference[tips['GE2D'][:, 0] <= 0.5].max())
exu.Print('tip at t=1:', dict((kind, tips[kind][-1, 1:3].round(6).tolist()) for kind in tips))

testResult = sum(tips[kind][-1, 1] + tips[kind][-1, 2] for kind in tips)
exu.Print('solution of flexiblePendulumBeamComparison=', testResult)
exu.sys['testResult'] = testResult
```