geometricallyExactBeamMarkerTest.py

You can view and download this file on Github: geometricallyExactBeamMarkerTest.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Body markers on ObjectBeamGeometricallyExact (#2730, RG4.8.7 of revision2026b): loads through
  5#           MarkerBodyMass, MarkerBodyPosition and MarkerBodyRigid on the elements give the same result as
  6#           the equivalent loads on the nodes. (1) a flexible pendulum under gravity through
  7#           LoadMassProportional, against gravity as nodal forces; (2) a cantilever with a tip force
  8#           and a tip torque on the last element's end, against the same on the last node; (3) a force
  9#           at the middle of an element, against half of it on each of its nodes.
 10#
 11# Author:   Johannes Gerstmayr
 12# Date:     2026-09-29
 13#
 14# 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.
 15#
 16#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 17
 18import exudyn as exu
 19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 20
 21import numpy as np
 22
 23testIsActive = exu.sys.get('testIsActive', False)
 24
 25L = 0.5; nElements = 8; E = 1e8; rho = 1000; h = 0.002; b = 0.01; nu = 0.3
 26A = b*h; Izz = b*h**3/12; Iyy = h*b**3/12; G = E/(2*(1+nu)); g = 9.81
 27lElement = L/nElements
 28
 29
 30def Beam(mbs):
 31    """nodes and elements of a straight beam along x, and a marker for the first node's coordinates"""
 32    nodes = [mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[i*lElement, 0, 0] + eulerParameters0))
 33             for i in range(nElements+1)]
 34    section = exu.BeamSection()
 35    section.stiffnessMatrix = np.diag([E*A, G*A, G*A, G*(Iyy+Izz), E*Iyy, E*Izz])
 36    section.inertia = np.diag([rho*(Iyy+Izz), rho*Iyy, rho*Izz])
 37    section.massPerLength = rho*A
 38    elements = [mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[nodes[i], nodes[i+1]], length=lElement,
 39                                                           sectionData=section)) for i in range(nElements)]
 40    mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=mbs.AddNode(NodePointGround()), coordinate=0))
 41    return (nodes, elements, mGround)
 42
 43
 44def Fix(mbs, mGround, node, coordinates):
 45    for i in coordinates:
 46        mbs.AddObject(ObjectConnectorCoordinate(markerNumbers=[mGround, mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=node, coordinate=i))]))
 47
 48
 49#(1) the pendulum: gravity through the elements or on the nodes
 50def Pendulum(throughElements):
 51    SC = exu.SystemContainer(); mbs = SC.AddSystem()
 52    (nodes, elements, mGround) = Beam(mbs)
 53    if throughElements:
 54        for e in elements:
 55            mbs.AddLoad(LoadMassProportional(markerNumber=mbs.AddMarker(MarkerBodyMass(bodyNumber=e)), loadVector=[0, -g, 0]))
 56    else:
 57        for i in range(nElements+1):
 58            weight = rho*A*lElement*(0.5 if i in [0, nElements] else 1)*g
 59            mbs.AddLoad(LoadForceVector(markerNumber=mbs.AddMarker(MarkerNodePosition(nodeNumber=nodes[i])), loadVector=[0, -weight, 0]))
 60    Fix(mbs, mGround, nodes[0], [0, 1, 2])
 61    mbs.Assemble()
 62    s = exu.SimulationSettings()
 63    s.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
 64    s.timeIntegration.numberOfSteps = 200; s.timeIntegration.endTime = 0.5
 65    s.solution.file.write = False
 66    s.linearSolver.solverType = exu.LinearSolverType.EigenSparse
 67    mbs.SolveDynamic(s)
 68    return mbs.GetNodeOutput(nodes[-1], exu.OutputVariableType.Position)
 69
 70
 71#(2), (3) the cantilever, clamped (position and e1..e3; e0 follows from the norm constraint)
 72def Cantilever(case, throughElements):
 73    SC = exu.SystemContainer(); mbs = SC.AddSystem()
 74    (nodes, elements, mGround) = Beam(mbs)
 75    F = [0, -0.5*E*Izz/L**2, 0]; M = [0, 0, 0.3*E*Izz/L]
 76    if case == 'tip':
 77        if throughElements:
 78            mTip = mbs.AddMarker(MarkerBodyRigid(bodyNumber=elements[-1], localPosition=[0.5*lElement, 0, 0]))
 79        else:
 80            mTip = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodes[-1]))
 81        mbs.AddLoad(LoadForceVector(markerNumber=mTip, loadVector=F))
 82        mbs.AddLoad(LoadTorqueVector(markerNumber=mTip, loadVector=M))
 83    else:   #the middle of the last element, or half on each of its nodes
 84        if throughElements:
 85            mbs.AddLoad(LoadForceVector(markerNumber=mbs.AddMarker(MarkerBodyPosition(bodyNumber=elements[-1])), loadVector=F))
 86        else:
 87            for n in nodes[-2:]:
 88                mbs.AddLoad(LoadForceVector(markerNumber=mbs.AddMarker(MarkerNodePosition(nodeNumber=n)), loadVector=0.5*np.array(F)))
 89    Fix(mbs, mGround, nodes[0], [0, 1, 2, 4, 5, 6])
 90    mbs.Assemble()
 91    s = exu.SimulationSettings()
 92    s.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
 93    s.staticSolver.numberOfLoadSteps = 10
 94    s.staticSolver.newton.relativeTolerance = 1e-10
 95    mbs.SolveStatic(s)
 96    return mbs.GetNodeOutput(nodes[-1], exu.OutputVariableType.Position)
 97
 98
 99testResult = 0
100for (name, solve) in [('pendulum', lambda through: Pendulum(through)),
101                      ('cantilever, tip', lambda through: Cantilever('tip', through)),
102                      ('cantilever, element middle', lambda through: Cantilever('middle', through))]:
103    pElements = solve(True)
104    pNodes = solve(False)
105    exu.Print(name + ': tip', pElements.round(8), ', difference element loads - node loads:', np.linalg.norm(pElements - pNodes))
106    testResult += sum(pElements)
107
108exu.Print('solution of geometricallyExactBeamMarkerTest=', testResult)
109exu.sys['testResult'] = testResult