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

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

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

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  Kinetic and potential (elastic) energy of flexible bodies as output variables (#2202): free,
#           undamped oscillations of two-element models of each beam and plate element, and a kinematic
#           tree under gravity with joint springs and constant joint forces. Each model gets sensors of
#           KineticEnergy and PotentialEnergy for every element, and the total is kept by the
#           generalized-alpha method without numerical damping (spectralRadius = 1) - up to the time
#           discretization of the nonlinear oscillation, which the printed relative change shows.
#           The elastic energy of an element is integrated with the strains and the integration rules
#           of its elastic forces, so that its derivative is exactly the elastic force.
#           (1) ANCFCable2D, (2) ANCFCable, (3) BeamGeometricallyExact2D, (4) BeamGeometricallyExact,
#           (5) ANCFBeam, (6) ANCFThinPlate, (7) KinematicTree: a revolute and a prismatic joint with
#           P control springs, constant joint forces and gravity.
#
# Author:   Johannes Gerstmayr
# Date:     2026-10-01
#
# 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
from exudyn.shells import ShellMesh

import numpy as np

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

SC = exu.SystemContainer()
mbs = SC.AddSystem()

KE = exu.OutputVariableType.KineticEnergy
PE = exu.OutputVariableType.PotentialEnergy
models = {}     #name: list of energy sensors; their sum is the total energy of the model

def EnergySensors(name, bodies):
    models[name] = [mbs.AddSensor(SensorBody(bodyNumber=b, outputVariableType=ovt, storeInternal=True))
                    for b in bodies for ovt in [KE, PE]]

L = 0.5         #element length
section = exu.BeamSection()
section.stiffnessMatrix = np.diag([100, 60, 80, 1, 1.5, 2])
section.inertia = np.diag([2e-3, 1e-3, 1e-3])
section.massPerLength = 2

#(1) ANCFCable2D: three nodes, the middle one moving sideways
nodes = [mbs.AddNode(NodePoint2DSlope1(referenceCoordinates=[i*L, 0, 1, 0], initialVelocities=[0, 0.2 if i == 1 else 0, 0, 0]))
         for i in range(3)]
EnergySensors('ANCFCable2D', [mbs.AddObject(ObjectANCFCable2D(nodeNumbers=[nodes[i], nodes[i+1]], length=L,
              massPerLength=2, bendingStiffness=1, axialStiffness=100)) for i in range(2)])

#(2) ANCFCable (3D)
nodes = [mbs.AddNode(NodePointSlope1(referenceCoordinates=[i*L, 1, 0, 1, 0, 0],
                                     initialVelocities=[0, 0.2 if i == 1 else 0, 0.1 if i == 2 else 0, 0, 0, 0]))
         for i in range(3)]
EnergySensors('ANCFCable', [mbs.AddObject(ObjectANCFCable(nodeNumbers=[nodes[i], nodes[i+1]], length=L,
              massPerLength=2, bendingStiffness=1, axialStiffness=100)) for i in range(2)])

#(3) BeamGeometricallyExact2D
nodes = [mbs.AddNode(NodeRigidBody2D(referenceCoordinates=[i*L, 2, 0], initialVelocities=[0, 0.2 if i == 1 else 0, 0.5 if i == 0 else 0]))
         for i in range(3)]
EnergySensors('BeamGeometricallyExact2D', [mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers=[nodes[i], nodes[i+1]],
              length=L, massPerLength=2, crossSectionInertia=2e-3, bendingStiffness=1,
              axialStiffness=100, shearStiffness=60)) for i in range(2)])

#(4) BeamGeometricallyExact (3D), Euler parameter nodes
nodes = [mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[i*L, 3, 0] + list(eulerParameters0),
                                     initialVelocities=[0, 0.2 if i == 1 else 0, 0.1 if i == 2 else 0] + list(
                                         AngularVelocity2EulerParameters_t([0.3, 0, 0] if i == 0 else [0, 0, 0], eulerParameters0))))
         for i in range(3)]
EnergySensors('BeamGeometricallyExact', [mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[nodes[i], nodes[i+1]],
              length=L, sectionData=section)) for i in range(2)])

#(5) ANCFBeam (3D), with the cross-section deformation as a penalty
nodes = [mbs.AddNode(NodePointSlope23(referenceCoordinates=[i*L, 4, 0, 0, 1, 0, 0, 0, 1],
                                      initialVelocities=[0, 0.2 if i == 1 else 0, 0.1 if i == 2 else 0, 0, 0, 0, 0, 0, 0]))
         for i in range(3)]
EnergySensors('ANCFBeam', [mbs.AddObject(ObjectANCFBeam(nodeNumbers=[nodes[i], nodes[i+1]], length=L, sectionData=section,
              crossSectionPenaltyFactor=[1, 1, 1])) for i in range(2)])

#(6) ANCFThinPlate: 2x1 elements, one corner moving out of the plane
plate = ShellMesh(vertices=[[0, 5, 0], [1, 5, 0], [1, 5.5, 0], [0, 5.5, 0]], numberOfElementsX=2, numberOfElementsY=1,
                  youngsModulus=2e5, poissonsRatio=0.3, density=500, thickness=0.02)
plate.CreateANCFThinPlateElements(mbs)
mbs.SetNodeParameter(plate.vertexNodeNumbers[2], 'initialVelocities', [0, 0, 0.1] + [0]*6)
EnergySensors('ANCFThinPlate', plate.elementNumbers)

#(7) KinematicTree: a revolute joint about z with a link along x, then a prismatic joint along x; P control springs
#towards offsets, constant joint forces and gravity in -y; no D control, so the total of T and V is kept
nKT = mbs.AddNode(NodeGenericODE2(referenceCoordinates=[0, 0.2], initialCoordinates=[0.3, 0], initialCoordinates_t=[0, 0.1],
                                  numberOfODE2Coordinates=2))
oKT = mbs.AddObject(ObjectKinematicTree(nodeNumber=nKT, baseOffset=[0, 7, 0], linkParents=[-1, 0],
                                        jointTypes=[exu.JointType.RevoluteZ, exu.JointType.PrismaticX],
                                        jointTransformations=exu.Matrix3DList([np.eye(3)]*2),
                                        jointOffsets=exu.Vector3DList([[0, 0, 0], [0.5, 0, 0]]),
                                        linkInertiasCOM=exu.Matrix3DList([np.eye(3)*0.01]*2),
                                        linkCOMs=exu.Vector3DList([[0.25, 0, 0], [0.1, 0, 0]]),
                                        linkMasses=[1, 2], gravity=[0, -9.81, 0],
                                        jointPositionOffsetVector=[0.1, 0.2], jointPControlVector=[50, 200],
                                        jointVelocityOffsetVector=[0, 0], jointDControlVector=[0, 0], jointForceVector=[2, -1]))
EnergySensors('KinematicTree', [oKT])

mbs.Assemble()

simulationSettings = exu.SimulationSettings()
simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
simulationSettings.timeIntegration.numberOfSteps = 1000
simulationSettings.timeIntegration.endTime = 1
simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
#the cross-section deformation of ObjectANCFBeam carries no mass, so its mass matrix is singular for initial accelerations
simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = False
simulationSettings.timeIntegration.verboseMode = 0
simulationSettings.solution.file.write = False
simulationSettings.solution.sensors.writePeriod = 0.01
mbs.SolveDynamic(simulationSettings)

testResult = 0
for (name, sensors) in models.items():
    total = sum(mbs.GetSensorStoredData(s)[:, 1] for s in sensors)
    elastic = sum(mbs.GetSensorStoredData(s)[:, 1] for s in sensors[1::2])
    exu.Print('%-28s total energy %.8f -> %.8f, relative change %.2e (largest elastic energy %.4f)' % (
        name, total[0], total[-1], (total[-1]-total[0])/total[0], elastic.max()))
    testResult += total[-1]

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