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

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

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

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  Lateral buckling of the right-angle frame (Argyris et al. 1979; Simo and Vu-Quoc 1986;
#           Nachbagauer and Gerstmayr 2013), with ObjectANCFBeam or ObjectBeamGeometricallyExact:
#           two arms of length 0.24 m with a thin rectangular cross section 30 x 0.6 mm, clamped at the
#           end of the first arm, joined rigidly at the corner. The tip of the second arm is DRIVEN by
#           displacement: a coordinate constraint prescribes its displacement along the first arm, and
#           its reaction force is the load P. So the static solver passes the buckling point, where a
#           load-driven computation stops at the critical load 1.088 N, and follows the post-buckling
#           path, on which P stays close to 1.088 N while the tip moves out of the plane. A small
#           out-of-plane force is the imperfection.
#           With ObjectANCFBeam (useGeometricallyExact = False) the static solver stops at about 3 % of
#           the drive (P = 0.31 N): there the Newton iteration converges only from very close to the
#           solution, although its Jacobian is consistent and the constrained tangent stiffness is
#           positive definite; neither more elements, another penalty factor, another linear solver
#           nor more load steps change it, and without the imperfection the drive is completed (#2763).
#           K. Nachbagauer, J. Gerstmayr. Structural and Continuum Mechanics Approaches for a 3D Shear
#           Deformable ANCF Beam Finite Element: Application to Buckling and Nonlinear Dynamic Examples,
#           Journal of Computational and Nonlinear Dynamics, Vol. 9(1), pp. 011013-1 - 011013-8, 2013.
#
# Author:   Johannes Gerstmayr
# Date:     2023-05-02, displacement-driven 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
import exudyn.graphics as graphics

import numpy as np

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

useGeometricallyExact = True    #False: ObjectANCFBeam, whose static solve stops at about 3 % of the drive (#2763, see Details)
nElements = 16                  #per arm
uxMax = 2e-3                    #prescribed displacement of the tip along the first arm (x); 200 steps of 1e-5 follow the post-buckling path, much larger steps may jump to another branch
imperfection = 1e-4             #out-of-plane force at the tip [N]
numberOfLoadSteps = 200

L = 0.24                #length of each arm
w = 0.0006              #thickness, out of the plane of the frame
h = 0.03                #height, in the plane of the frame
Em = 7.124e10; nu = 0.31; Gm = Em/(2*(1+nu))
A = h*w; Iyy = h*w**3/12; Izz = w*h**3/12
J = h*w**3/3            #torsion of a thin rectangle
rho = 1e2               #not used by the static solver

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

sectionData = exu.BeamSection()
sectionData.stiffnessMatrix = np.diag([Em*A, Gm*A, Gm*A, Gm*J, Em*Iyy, Em*Izz])
sectionData.inertia = rho*J*np.eye(3)
sectionData.massPerLength = rho*A

sectionGeometry = exu.BeamSectionGeometry()
sectionGeometry.polygonalPoints = exu.Vector2DList()
for p in [[h, -w], [h, w], [-h, w], [-h, -w]]:
    sectionGeometry.polygonalPoints.Append(p)

lElem = L/nElements
rotZ = RotationMatrixZ(-np.pi/2)

def Arm(position0, direction, rotation, slopes):
    """the nodes and elements of one arm; returns the first and the last node"""
    nodes = []
    for k in range(nElements+1):
        p = list(np.array(position0) + k*lElem*np.array(direction))
        if useGeometricallyExact:
            nodes += [mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=p+list(RotationMatrix2EulerParameters(rotation))))]
        else:
            nodes += [mbs.AddNode(NodePointSlope23(referenceCoordinates=p+slopes))]
    for k in range(nElements):
        if useGeometricallyExact:
            mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[nodes[k], nodes[k+1]], length=lElem,
                                                       sectionData=sectionData,
                                                       visualization=VBeam3D(sectionGeometry=sectionGeometry)))
        else:
            mbs.AddObject(ObjectANCFBeam(nodeNumbers=[nodes[k], nodes[k+1]], length=lElem,
                                         sectionData=sectionData, crossSectionPenaltyFactor=[1, 1, 1],
                                         visualization=VANCFBeam(sectionGeometry=sectionGeometry)))
    return nodes[0], nodes[-1]

nInit, n1 = Arm([0, 0, 0], [1, 0, 0], np.eye(3), [0, 1, 0, 0, 0, 1])      #first arm, along x
nInitB, n1B = Arm([L, 0, 0], [0, -1, 0], rotZ, [1, 0, 0, 0, 0, 1])       #second arm, along -y

nGround = mbs.AddNode(NodePointGround(visualization=VNodePointGround(show=False)))
mCoordinateGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))

#clamping at the start of the first arm, and the rigid corner
if useGeometricallyExact:
    mbs.AddObject(GenericJoint(markerNumbers=[mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround)),
                                              mbs.AddMarker(MarkerNodeRigid(nodeNumber=nInit))],
                               visualization=VGenericJoint(axesRadius=2*w, axesLength=2*w)))
else:
    for i in range(9):
        mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround,
                                                          mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nInit, coordinate=i))]))
mbs.AddObject(GenericJoint(markerNumbers=[mbs.AddMarker(MarkerNodeRigid(nodeNumber=n1)),
                                          mbs.AddMarker(MarkerNodeRigid(nodeNumber=nInitB, localHT=exu.HT(rotation=rotZ)))],
                           visualization=VGenericJoint(axesRadius=2*w, axesLength=2*w)))

#the tip is driven along x; the reaction force of the constraint is the load P
def UFoffset(mbs, t, itemNumber, lOffset):
    return uxMax*t

oDrive = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround,
                                                           mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n1B, coordinate=0))],
                                            offsetUserFunction=UFoffset))
mbs.AddLoad(Force(markerNumber=mbs.AddMarker(MarkerNodePosition(nodeNumber=n1B)), loadVector=[0, 0, imperfection]))

sForce = mbs.AddSensor(SensorObject(objectNumber=oDrive, outputVariableType=exu.OutputVariableType.Force,
                                    storeInternal=True, writeToFile=False))
sDisp = mbs.AddSensor(SensorNode(nodeNumber=n1B, outputVariableType=exu.OutputVariableType.Displacement,
                                 storeInternal=True, writeToFile=False))
mbs.Assemble()

simulationSettings = exu.SimulationSettings()
simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
simulationSettings.staticSolver.useLoadFactor = False   #the drive follows the user function of t in [0, 1]
simulationSettings.staticSolver.numberOfLoadSteps = numberOfLoadSteps
simulationSettings.staticSolver.newton.relativeTolerance = 1e-8
simulationSettings.staticSolver.newton.absoluteTolerance = 1e-7
simulationSettings.solution.file.write = False
simulationSettings.solution.sensors.writePeriod = 1/numberOfLoadSteps

SC.visualizationSettings.bodies.beams.axialTiling = 50
SC.visualizationSettings.view0.scene.drawWorldBasis = True
SC.visualizationSettings.view0.scene.worldBasisSize = 0.1

if not testIsActive:
    SC.renderer.Start()

mbs.SolveStatic(simulationSettings)

if not testIsActive:
    SC.renderer.DoIdleTasks()
    SC.renderer.Stop()

force = mbs.GetSensorStoredData(sForce)[:, 1]       #the load P on the tip: the force of the drive
uz = mbs.GetSensorStoredData(sDisp)[:, 3]           #out-of-plane displacement of the tip
testResult = 0
for displacement in [1e-3, 5e-3, 1e-2]:
    forceAt = np.interp(displacement, np.abs(uz), force)
    exu.Print('out-of-plane tip displacement', displacement, 'at P =', round(forceAt, 4))
    testResult += forceAt
exu.Print('at the end: P =', round(force[-1], 4), ', out-of-plane tip displacement', round(uz[-1], 5))
testResult += force[-1] + abs(uz[-1])

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

if not testIsActive:
    mbs.PlotSensor(sensorNumbers=[np.column_stack((np.abs(uz), force)), np.array([[0, 1.088], [np.abs(uz).max(), 1.088]])],
                   components=[0, 0], componentsX=[-1, -1], labels=['computed', 'critical load 1.088 N'],
                   xLabel='out-of-plane tip displacement', yLabel='P')
```