pendulum2Dconstraint.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Mathematical pendulum with constraint;
  5#           Remark: update from pendulum.py example
  6#
  7# Author:   Johannes Gerstmayr
  8# Date:     2019-12-26
  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 = 0.8 #distance
 22mass = 2.5
 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
 45sDist = mbs.AddSensor(SensorObject(objectNumber=oDistance, storeInternal=True,
 46                                   outputVariableType=exu.OutputVariableType.Distance))
 47
 48#print(mbs)
 49
 50mbs.Assemble()
 51
 52simulationSettings = exu.SimulationSettings()
 53
 54f = 1000000
 55simulationSettings.timeIntegration.numberOfSteps = int(1*f)
 56simulationSettings.timeIntegration.endTime = 0.001*f
 57simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/5000
 58simulationSettings.solutionSettings.sensorsWritePeriod = simulationSettings.timeIntegration.endTime/50000
 59#simulationSettings.displayComputationTime = True
 60simulationSettings.timeIntegration.verboseMode = 1
 61simulationSettings.timeIntegration.verboseModeFile = 0
 62
 63#these Newton settings are slightly faster than full Newton:
 64simulationSettings.timeIntegration.newton.useModifiedNewton = True
 65simulationSettings.timeIntegration.newton.modifiedNewtonJacUpdatePerStep = True
 66
 67simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.60 #0.62 is approx. the limit
 68simulationSettings.timeIntegration.adaptiveStep = False
 69
 70simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = True
 71simulationSettings.solutionSettings.coordinatesSolutionFileName= "coordinatesSolution.txt"
 72
 73simulationSettings.displayStatistics = True
 74#simulationSettings.solutionSettings.recordImagesInterval = 0.04
 75
 76SC.visualizationSettings.nodes.defaultSize = 0.05
 77exu.StartRenderer()
 78
 79#mbs.WaitForUserToContinue()
 80#exu.InfoStat()
 81mbs.SolveDynamic(simulationSettings,
 82                 # solverType=exu.DynamicSolverType.TrapezoidalIndex2
 83                 )
 84#exu.InfoStat()
 85
 86SC.WaitForRenderEngineStopFlag()
 87exu.StopRenderer() #safely close rendering window!
 88
 89nODE2 = len(mbs.systemData.GetODE2Coordinates())
 90print("ODE2=",nODE2)
 91
 92#plot constraint error:
 93
 94mbs.PlotSensor(sensorNumbers=sDist, offsets=[-L], closeAll=True)
 95
 96#old way, better use PlotSensor:
 97import matplotlib.pyplot as plt
 98import matplotlib.ticker as ticker
 99
100#plot y-acceleration:
101data = np.loadtxt('coordinatesSolution.txt', comments='#', delimiter=',')
102plt.figure()
103plt.plot(data[:,0], data[:,1+2*nODE2+1], 'b-')
104
105ax=plt.gca() # get current axes
106ax.grid(True, 'major', 'both')
107ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
108ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
109plt.tight_layout()
110plt.show()