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

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

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

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  The section forces, moments and strains of ObjectBeamGeometricallyExact (#2753): a cantilever
#           bent by large tip loads. (1) A bending moment, then a torsional moment at the tip: every element
#           carries it as TorqueLocal, and CurvatureLocal is the moment divided by the stiffness, exactly.
#           (2) A tip force: in every element, the section force ForceLocal rotated by the rotation matrix
#           at the middle of the element gives the tip force of 4.7 - to 5.6e-3 with 8 elements, 3.9e-4 with 32:
#           the section forces of an element are constant, the mean of a field that varies along it.
#           Also the other output variables at the ends of an element, against its nodes.
#
# 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
EA = 1e5; GA = 4e4; GJ = 2; EIy = 5; EIz = 3


def Cantilever(force, torque, nElements=8):
    """solve statically; return mbs, the elements and the nodes"""
    lElement = L/nElements
    SC = exu.SystemContainer(); mbs = SC.AddSystem()
    nodes = [mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[i*lElement, 0, 0]+list(eulerParameters0)))
             for i in range(nElements+1)]
    section = exu.BeamSection()
    section.stiffnessMatrix = np.diag([EA, GA, GA, GJ, EIy, EIz])
    section.inertia = np.diag([2e-3, 1e-3, 1e-3])
    section.massPerLength = 1
    elements = [mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[nodes[i], nodes[i+1]], length=lElement,
                                                           sectionData=section)) for i in range(nElements)]
    mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=mbs.AddNode(NodePointGround()), coordinate=0))
    for i in [0, 1, 2, 4, 5, 6]:
        mbs.AddObject(ObjectConnectorCoordinate(markerNumbers=[mGround, mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nodes[0], coordinate=i))]))
    mTip = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodes[-1]))
    mbs.AddLoad(LoadForceVector(markerNumber=mTip, loadVector=force))
    mbs.AddLoad(LoadTorqueVector(markerNumber=mTip, loadVector=torque, bodyFixed=True))
    mbs.Assemble()
    simulationSettings = exu.SimulationSettings()
    simulationSettings.staticSolver.numberOfLoadSteps = 10
    simulationSettings.staticSolver.newton.relativeTolerance = 1e-10
    simulationSettings.staticSolver.newton.absoluteTolerance = 1e-12
    mbs.SolveStatic(simulationSettings)
    return (mbs, elements, nodes)


testResult = 0
Out = lambda mbs, e, t, x=0: mbs.GetObjectOutputBody(e, t, localPosition=[x, 0, 0])

#(1) tip moments, one at a time, given in the frame of the tip: constant along the beam
for torque in [[0, 0, 0.5*EIz/L], [0.4*GJ/L, 0, 0]]:
    (mbs, elements, nodes) = Cantilever([0, 0, 0], torque)
    errorMoment = max(np.linalg.norm(Out(mbs, e, exu.OutputVariableType.TorqueLocal) - torque) for e in elements)
    errorCurvature = max(np.linalg.norm(Out(mbs, e, exu.OutputVariableType.CurvatureLocal) - np.array(torque)/[GJ, EIy, EIz]) for e in elements)
    exu.Print('tip moment', torque, ': largest error of TorqueLocal', errorMoment, ', of CurvatureLocal', errorCurvature)
    testResult += sum(Out(mbs, elements[-1], exu.OutputVariableType.CurvatureLocal))

#(2) tip force: every section carries it
force = [0, -1.5*EIz/L**2, 0.5*EIy/L**2]
for nElements in [32, 8]:
    (mbs, elements, nodes) = Cantilever(force, [0, 0, 0], nElements)
    errorForce = 0
    for e in elements:
        rotation = Out(mbs, e, exu.OutputVariableType.RotationMatrix).reshape(3, 3)
        errorForce = max(errorForce, np.linalg.norm(rotation @ Out(mbs, e, exu.OutputVariableType.ForceLocal) - force))
    exu.Print(nElements, 'elements, tip force: largest error of the section force, rotated to global', errorForce)
testResult += sum(Out(mbs, elements[0], exu.OutputVariableType.TorqueLocal))
testResult += sum(Out(mbs, elements[0], exu.OutputVariableType.StrainLocal))

#the output at the ends of an element is the one of its nodes
errorEnds = 0
for (x, node) in [(-0.5*L/nElements, nodes[-2]), (0.5*L/nElements, nodes[-1])]:
    for outputVariable in [exu.OutputVariableType.Position, exu.OutputVariableType.RotationMatrix,
                           exu.OutputVariableType.AngularVelocity]:
        errorEnds = max(errorEnds, np.linalg.norm(Out(mbs, elements[-1], outputVariable, x) - mbs.GetNodeOutput(node, outputVariable)))
exu.Print('element ends against the nodes: largest difference', errorEnds)
rotation = Out(mbs, elements[-1], exu.OutputVariableType.Rotation, 0.5*L/nElements)
testResult += sum(rotation)

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