RigidPendulumTest.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  A `RigidBody2D` hanging in a `RevoluteJoint2D` and swinging under gravity, integrated
 5#           over 0.5 s: the vertical position of the tip at the end.
 6#           The model compares against a reference value written into it, so its
 7#           result is that difference and its reference solution is 0 (#2632).
 8#
 9# Author:   Johannes Gerstmayr
10# Date:     2019-11-01, reworked 2026-09-24
11#
12# 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.
13#
14#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15
16import exudyn as exu
17from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
18import exudyn.graphics as graphics
19
20testIsActive = exu.sys.get('testIsActive', False)
21exu.sys['testTolerance'] = 4e-13 #the tolerance RunAllModelUnitTests used for these ten
22
23SC = exu.SystemContainer()
24mbs = SC.AddSystem()
25
26rect = [-2,-2,2,2] #xmin,ymin,xmax,ymax
27background = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[rect[0],rect[1],0, rect[2],rect[1],0, rect[2],rect[3],0, rect[0],rect[3],0, rect[0],rect[1],0]} #background
28oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
29a = 0.5     #half x-dim of pendulum
30b = 0.05    #half y-dim of pendulum
31massRigid = 12
32inertiaRigid = massRigid/12*(2*a)**2
33g = 9.81    # gravity
34
35graphics2 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[-a,-b,0, a,-b,0, a,b,0, -a,b,0, -a,-b,0]} #background
36nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[a,0,0], initialVelocities=[0,0,0*2]));
37oRigid = mbs.AddObject(RigidBody2D(mass=massRigid, inertia=inertiaRigid,nodeNumber=nRigid,visualization=VObjectRigidBody2D(graphicsData= [graphics2])))
38
39mRigidSupport = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[-a,0.,0.])) #support point
40mRigidMidPoint = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ 0., 0.,0.])) #mid point
41mRigidMass = mbs.AddMarker(MarkerBodyMass(bodyNumber=oRigid)) 
42
43mGround = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=[0,0,0.]))
44mbs.AddObject(RevoluteJoint2D(markerNumbers=[mGround,mRigidSupport]))
45
46mbs.AddLoad(Force(markerNumber = mRigidMidPoint, loadVector = [0, -0.5*massRigid*g, 0])) #split force into two parts: gravity and force ...
47mbs.AddLoad(Gravity(markerNumber = mRigidMass, loadVector = [0, -0.5*g, 0])) 
48
49mbs.Assemble()
50#mbs.systemData.Info()
51
52simulationSettings = exu.SimulationSettings() #takes currently set values or default values
53simulationSettings.solution.file.write=False
54
55simulationSettings.timeIntegration.numberOfSteps = 1000
56simulationSettings.timeIntegration.endTime = 0.5
57simulationSettings.timeIntegration.newton.relativeTolerance = 1e-10 
58simulationSettings.timeIntegration.verboseMode = 1 
59simulationSettings.show.statistics = False
60
61simulationSettings.solution.file.writePeriod = 1e-4
62
63simulationSettings.timeIntegration.newton.useModifiedNewton = True
64
65#simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
66#simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
67simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
68
69if not testIsActive: 
70    SC.renderer.Start()
71
72mbs.SolveDynamic(simulationSettings)
73
74if not testIsActive: 
75    SC.renderer.DoIdleTasks()
76    SC.renderer.Stop() #safely close rendering window!
77
78
79u = mbs.GetNodeOutput(nRigid, exu.OutputVariableType.Position) #tip node
80# if True:
81#     errorRigidPendulum = u[1] - (-0.49796067298096375 ) #2021-02-06: --0.49796067298096375 
82if True:
83    errorRigidPendulum = u[1] - (-0.4979606729809297) #2021-02-04: -0.4979606729809297
84else:
85    errorRigidPendulum = u[1] - (-0.4979662392961769) #2019-12-26(new initial acc): -0.4979662392961769; 15.12.2019: (-0.4980200584148534); before 15.12.2019: (-0.4980200584133354) # old test (load at tip)   0*(- 0.4905431986572512) #0*(-0.037344780490849015) #yield-displacement
86exu.Print('solution rigid pendulum=',u[1])
87
88testResult = abs(errorRigidPendulum)
89exu.Print('solution of RigidPendulumTest=', testResult)
90exu.sys['testResult'] = testResult