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)