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