accessFunctionsTest.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  A body marker places a force, a torque or a connector at a local position of a body, and the
 5#           body computes the Jacobians there (its access functions, #2744). Some bodies define them at
 6#           restricted local positions only: ObjectANCFCable at the centerline, ObjectANCFThinPlate at the
 7#           midsurface, ObjectBeamGeometricallyExact2D at the axis. Assemble() checks the markers through which
 8#           a connector or a load acts:
 9#           (1) a rotating table, ObjectRotationalMass1D, driven by a torque on its axis and by a force at its
10#               rim, fixed to the table and tangential, whose moment about the axis it takes (#2775);
11#           (2) a cantilever ANCF cable under a force at its tip, on the centerline;
12#           (3) what Assemble() refuses: a force off the centerline of the cable.
13#
14# Author:   Johannes Gerstmayr
15# Date:     2026-10-02
16#
17# 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.
18#
19#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
20
21import exudyn as exu
22from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
23import numpy as np
24
25testIsActive = exu.sys.get('testIsActive', False)
26
27def TableModel():
28    """a rotating table with inertia 0.5, driven by a torque 1 about its axis and a tangential force 1, fixed to the
29    table, at its rim at radius 0.2"""
30    SC = exu.SystemContainer()
31    mbs = SC.AddSystem()
32    nTable = mbs.AddNode(Node1D(referenceCoordinates=[0]))
33    oTable = mbs.AddObject(ObjectRotationalMass1D(inertia=0.5, nodeNumber=nTable, referencePosition=[0, 0, 0]))
34    mAxis = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oTable, localPosition=[0, 0, 0]))
35    mRim = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oTable, localPosition=[0.2, 0, 0]))
36    mbs.AddLoad(LoadTorqueVector(markerNumber=mAxis, loadVector=[0, 0, 1]))
37    mbs.AddLoad(LoadForceVector(markerNumber=mRim, loadVector=[0, 1, 0], bodyFixed=True))
38    sRim = mbs.AddSensor(SensorMarker(markerNumber=mRim, outputVariableType=exu.OutputVariableType.Position, storeInternal=True))
39    mbs.Assemble()
40    return (mbs, sRim)
41
42def CableModel(offset):
43    """a cantilever ANCF cable, one element of length 1 with EI = 1, a force -0.01 in y at the tip, at offset from
44    the centerline; clamped: the position and the lateral slopes of node 0 fixed"""
45    SC = exu.SystemContainer()
46    mbs = SC.AddSystem()
47    nodes = [mbs.AddNode(NodePointSlope1(referenceCoordinates=[x, 0, 0, 1, 0, 0])) for x in [0, 1]]
48    oCable = mbs.AddObject(ObjectANCFCable(nodeNumbers=nodes, length=1, massPerLength=1,
49                                           bendingStiffness=1, axialStiffness=1e4))
50    oGround = mbs.AddObject(ObjectGround())
51    nGround = mbs.AddNode(NodePointGround())
52    mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
53    for coordinate in [0, 1, 2, 4, 5]:
54        mCoordinate = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nodes[0], coordinate=coordinate))
55        mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround, mCoordinate]))
56    mTip = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oCable, localPosition=[1, offset, 0]))
57    mbs.AddLoad(LoadForceVector(markerNumber=mTip, loadVector=[0, -0.01, 0]))
58    mbs.Assemble()
59    return (mbs, nodes[-1])
60
61#(1) the torque 1 + 0.2*1 turns the table by phi = t^2/2 * 1.2/0.5 = 1.2 rad after 1 s; the rim follows
62(mbs, sRim) = TableModel()
63simulationSettings = exu.SimulationSettings()
64simulationSettings.timeIntegration.numberOfSteps = 100
65simulationSettings.timeIntegration.endTime = 1
66mbs.SolveDynamic(simulationSettings)
67rim = mbs.GetSensorValues(sRim)
68exu.Print('position of the rim after 1 s:', rim, ', expected', 0.2*np.array([np.cos(1.2), np.sin(1.2), 0]))
69
70#(2) the tip deflection, F L^3/(3 EI) = 0.01/3 for the Euler-Bernoulli beam
71(mbs, nTip) = CableModel(offset=0)
72mbs.SolveStatic()
73tip = mbs.GetNodeOutput(nTip, exu.OutputVariableType.Displacement)
74exu.Print('tip deflection of the cable:', tip[1], ', expected', -0.01/3)
75
76#(3) refused at Assemble(), with the reason
77refused = 0
78for (name, Build) in [('force off the centerline of the cable', lambda: CableModel(offset=0.1))]:
79    try:
80        Build()
81        exu.Print(name + ': accepted')
82    except Exception as error:
83        exu.Print(name + ': refused,', str(error).split(' [Python file')[0].strip())
84        refused += 1
85
86testResult = rim[0] + rim[1] + tip[1] + refused
87exu.Print('solution of accessFunctionsTest=', testResult)
88exu.sys['testResult'] = testResult