pendulumIftommBenchmark.py
You can view and download this file on Github: pendulumIftommBenchmark.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Pendulum according to IFToMM multibody benchmarks
5# https://www.iftomm-multibody.org/benchmark/problem/Planar_simple_pendulum/
6#
7# Author: Johannes Gerstmayr
8# Date: 2022-06-15
9#
10# 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.
11#
12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
13
14import exudyn as exu
15from exudyn.utilities import *
16
17SC = exu.SystemContainer()
18mbs = SC.AddSystem()
19
20L = 1 #distance
21mass = 1
22g = 9.81
23
24r = 0.05 #just for graphics
25graphicsBackground = GraphicsDataRectangle(-1.2*L,-1.2*L, 1.2*L, 0.2*L, [1,1,1,1]) #for appropriate zoom
26graphicsSphere = GraphicsDataSphere(point=[0,0,0], radius=r, color=[1.,0.2,0.2,1], nTiles = 16)
27
28#add ground object and mass point:
29oGround = mbs.AddObject(ObjectGround(referencePosition = [0,0,0],
30 visualization = VObjectGround(graphicsData = [graphicsBackground])))
31nMass = mbs.AddNode(NodePoint2D(referenceCoordinates=[-L,0],
32 initialCoordinates=[0,0],
33 initialVelocities=[0,0]))
34oMass = mbs.AddObject(MassPoint2D(physicsMass = mass, nodeNumber = nMass,
35 visualization = VObjectMassPoint2D(graphicsData = [graphicsSphere])))
36
37mMass = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
38mGround = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition = [0,0,0]))
39oDistance = mbs.AddObject(DistanceConstraint(markerNumbers = [mGround, mMass], distance = L))
40
41#add loads:
42mbs.AddLoad(Force(markerNumber = mMass, loadVector = [0, -mass*g, 0]))
43
44addSensors = True
45if addSensors:
46 sDist = mbs.AddSensor(SensorObject(objectNumber=oDistance, storeInternal=True,
47 outputVariableType=exu.OutputVariableType.Distance))
48 sPos = mbs.AddSensor(SensorNode(nodeNumber=nMass, storeInternal=True,
49 outputVariableType=exu.OutputVariableType.Position))
50
51 sVel = mbs.AddSensor(SensorNode(nodeNumber=nMass, storeInternal=True,
52 outputVariableType=exu.OutputVariableType.Velocity))
53
54
55 def UFenergy(mbs, t, sensorNumbers, factors, configuration):
56 pos = mbs.GetSensorValues(sensorNumbers[0])
57 vel = mbs.GetSensorValues(sensorNumbers[1])
58 Ekin = 0.5*mass*(vel[0]**2 + vel[1]**2)
59 Epot = mass * g * pos[1]
60 return [Ekin+Epot] #return total energy
61
62 sEnergy = mbs.AddSensor(SensorUserFunction(sensorNumbers=[sPos,sVel], storeInternal=True,
63 sensorUserFunction=UFenergy))
64
65#print(mbs)
66
67mbs.Assemble()
68
69simulationSettings = exu.SimulationSettings()
70
71#for performance (energy error < 5e-5J):
72#without sensors, takes 0.037 seconds on i7 surface book laptop
73endTime = 10
74stepSize = 0.8e-3
75
76#accuracy:
77# endTime = 9.99 #in benchmark results, values are only given until 9.99 seconds
78# stepSize = 1e-5
79
80
81simulationSettings.timeIntegration.numberOfSteps = int(endTime/stepSize)
82simulationSettings.timeIntegration.endTime = endTime
83simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/5
84simulationSettings.solutionSettings.sensorsWritePeriod = simulationSettings.timeIntegration.endTime/100
85#simulationSettings.displayComputationTime = True
86simulationSettings.timeIntegration.verboseMode = 1
87#simulationSettings.timeIntegration.verboseModeFile = 0
88
89#these Newton settings are slightly faster than full Newton:
90simulationSettings.timeIntegration.newton.useModifiedNewton = True
91simulationSettings.timeIntegration.newton.modifiedNewtonJacUpdatePerStep = True
92
93#simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.75
94#simulationSettings.timeIntegration.adaptiveStep = False
95
96#simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = True
97#simulationSettings.solutionSettings.coordinatesSolutionFileName= "coordinatesSolution.txt"
98
99simulationSettings.displayStatistics = True
100#simulationSettings.solutionSettings.recordImagesInterval = 0.04
101
102SC.visualizationSettings.nodes.defaultSize = 0.05
103useGraphics = False
104
105if useGraphics:
106 exu.StartRenderer()
107
108#mbs.WaitForUserToContinue()
109#exu.InfoStat()
110mbs.SolveDynamic(simulationSettings,
111 # solverType=exu.DynamicSolverType.TrapezoidalIndex2
112 )
113#exu.InfoStat()
114
115if useGraphics:
116 SC.WaitForRenderEngineStopFlag()
117 exu.StopRenderer() #safely close rendering window!
118
119
120#plot constraint error:
121if addSensors:
122
123 mbs.PlotSensor(sensorNumbers=sDist, offsets=[-L], closeAll=True)
124 mbs.PlotSensor(sensorNumbers=sPos, components=[0,1], newFigure=True)
125
126 mbs.PlotSensor(sensorNumbers=sEnergy, yLabel='total energy', newFigure=True)