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

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

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

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  The 3-node (quadratic) ObjectBeamGeometricallyExact2D (#2208): nodes at the start, the middle
#           and the end of the element, in this order. (1) A cantilever with a small tip force: one
#           quadratic element gives the Timoshenko deflection F*L^3/(3*EI) + F*L/GA = -3.43333e-4 to 1e-7,
#           the size of the geometrically nonlinear part (the linear element: -3.4203e-4 with 8). (2) The same cantilever bent into a half
#           circle by a tip moment pi*EI/L, tip at y = 2L/pi = 0.636620: 0.616191, 0.635647, 0.636563 and
#           0.636616 with 1, 2, 4 and 8 quadratic elements - fourth order - against 0.707107, 0.653281 and
#           0.640729 with 2, 4 and 8 linear elements. (3) A pendulum of 4 quadratic elements falling under
#           gravity through MarkerBodyMass, hinged at one end: tip y -0.001862 at t = 1, 8 and 16 elements
#           give -0.001827 and -0.001830.
#
# Author:   Johannes Gerstmayr
# Date:     2026-09-30
#
# 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 = 1; EI = 100; GA = 1e4; EA = 1e5


def Beam(nElements, force=0, moment=0, dynamic=False):
    """the tip position of a beam of quadratic elements: a clamped cantilever, or a pendulum under gravity"""
    SC = exu.SystemContainer(); mbs = SC.AddSystem()
    nodes = [mbs.AddNode(NodeRigidBody2D(referenceCoordinates=[L*i/(2*nElements), 0, 0])) for i in range(2*nElements+1)]
    elements = [mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers=nodes[2*k:2*k+3], length=L/nElements,
                              massPerLength=1, crossSectionInertia=0.01, bendingStiffness=EI,
                              axialStiffness=EA, shearStiffness=GA)) for k in range(nElements)]
    mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=mbs.AddNode(NodePointGround()), coordinate=0))
    for i in ([0, 1] if dynamic else [0, 1, 2]):
        mbs.AddObject(ObjectConnectorCoordinate(markerNumbers=[mGround, mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nodes[0], coordinate=i))]))
    simulationSettings = exu.SimulationSettings()
    simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
    simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
    if dynamic:
        for e in elements:
            mbs.AddLoad(LoadMassProportional(markerNumber=mbs.AddMarker(MarkerBodyMass(bodyNumber=e)), loadVector=[0, -9.81, 0]))
        mbs.Assemble()
        simulationSettings.timeIntegration.numberOfSteps = 1000
        simulationSettings.timeIntegration.endTime = 1
        simulationSettings.solution.file.write = False
        mbs.SolveDynamic(simulationSettings)
    else:
        mTip = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodes[-1]))
        mbs.AddLoad(LoadForceVector(markerNumber=mTip, loadVector=[0, force, 0]))
        mbs.AddLoad(LoadTorqueVector(markerNumber=mTip, loadVector=[0, 0, moment]))
        mbs.Assemble()
        simulationSettings.staticSolver.numberOfLoadSteps = 10 if moment != 0 else 1
        mbs.SolveStatic(simulationSettings)
    return mbs.GetNodeOutput(nodes[-1], exu.OutputVariableType.Position)


force = -0.1
tipForce = Beam(1, force=force)[1]
exu.Print('small tip force, one element: deflection', tipForce, ', Timoshenko', force*L**3/(3*EI)+force*L/GA)
testResult = 1000*tipForce

for nElements in [2, 4]:
    tipMoment = Beam(nElements, moment=np.pi*EI/L)[1]
    exu.Print('half circle,', nElements, 'elements: tip y', tipMoment, ', exact', 2*L/np.pi)
    testResult += tipMoment

tipPendulum = Beam(4, dynamic=True)[1]
exu.Print('pendulum, 4 elements: tip y at t = 1:', tipPendulum)
testResult += 100*tipPendulum

exu.Print('solution of geometricallyExactBeam2DquadraticTest=', testResult)
exu.sys['testResult'] = testResult
```