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
 15#a local run must not open a plot window when windows are suppressed (#2477)
 16if exu.special.userInterface.suppressPlots: import matplotlib; matplotlib.use('Agg')
 17from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 18import numpy as np
 19import exudyn.graphics as graphics
 20
 21SC = exu.SystemContainer()
 22mbs = SC.AddSystem()
 23
 24L = 0.8 #distance
 25mass = 2.5
 26g = 9.81
 27
 28r = 0.05 #just for graphics
 29graphicsBackground = graphics.Lines([[-1.2*L,-1.2*L,0], [1.2*L,-1.2*L,0], [1.2*L,0.2*L,0], [-1.2*L,0.2*L,0], [-1.2*L,-1.2*L,0]], color=[1,1,1,1]) #for appropriate zoom
 30graphicsSphere = graphics.Sphere(point=[0,0,0], radius=r, color=[1.,0.2,0.2,1], nTiles = 16)
 31
 32#add ground object and mass point:
 33oGround = mbs.AddObject(ObjectGround(referencePosition = [0,0,0], 
 34                           visualization = VObjectGround(graphicsData = [graphicsBackground])))
 35nMass = mbs.AddNode(NodePoint2D(referenceCoordinates=[L,0], 
 36                                initialCoordinates=[0,0],
 37                                initialVelocities=[0,0]))
 38oMass = mbs.AddObject(MassPoint2D(mass = mass, nodeNumber = nMass, 
 39                                  visualization = VObjectMassPoint2D(graphicsData = [graphicsSphere])))
 40
 41mMass = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
 42mGround = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition = [0,0,0]))
 43oDistance = mbs.AddObject(DistanceConstraint(markerNumbers = [mGround, mMass], distance = L))
 44
 45#add loads:
 46mbs.AddLoad(Force(markerNumber = mMass, loadVector = [0, -mass*g, 0])) 
 47
 48sDist = mbs.AddSensor(SensorObject(objectNumber=oDistance, storeInternal=True, 
 49                                   outputVariableType=exu.OutputVariableType.Distance))
 50
 51#print(mbs)
 52
 53mbs.Assemble()
 54
 55simulationSettings = exu.SimulationSettings()
 56
 57f = 1000000
 58simulationSettings.timeIntegration.numberOfSteps = int(1*f)
 59simulationSettings.timeIntegration.endTime = 0.001*f
 60simulationSettings.solution.file.writePeriod = simulationSettings.timeIntegration.endTime/5000
 61simulationSettings.solution.sensors.writePeriod = simulationSettings.timeIntegration.endTime/50000
 62#simulationSettings.show.computationTime = True
 63simulationSettings.timeIntegration.verboseMode = 1
 64simulationSettings.timeIntegration.verboseModeFile = 0
 65
 66#these Newton settings are slightly faster than full Newton:
 67simulationSettings.timeIntegration.newton.useModifiedNewton = True
 68simulationSettings.timeIntegration.newton.modifiedNewtonJacUpdatePerStep = True
 69
 70simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.60 #0.62 is approx. the limit
 71simulationSettings.timeIntegration.adaptiveStep = False
 72
 73simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = True
 74simulationSettings.solution.file.name= "solution/coordinatesSolution.txt"
 75
 76simulationSettings.show.statistics = True
 77#simulationSettings.solution.recordImagesInterval = 0.04
 78
 79SC.visualizationSettings.nodes.defaultSize = 0.05
 80SC.renderer.Start()
 81
 82#SC.renderer.DoIdleTasks()
 83mbs.SolveDynamic(simulationSettings, 
 84                 # solverType=exu.DynamicSolverType.TrapezoidalIndex2
 85                 )
 86
 87SC.renderer.DoIdleTasks()
 88SC.renderer.Stop() #safely close rendering window!
 89
 90nODE2 = len(mbs.systemData.GetODE2Coordinates())
 91print("ODE2=",nODE2)
 92
 93#plot constraint error:
 94
 95mbs.PlotSensor(sensorNumbers=sDist, offsets=[-L], closeAll=True)
 96
 97#old way, better use PlotSensor:
 98import matplotlib.pyplot as plt
 99import matplotlib.ticker as ticker
100
101#plot y-acceleration:
102data = np.loadtxt(OutputFilePath('solution/coordinatesSolution.txt'), comments='#', delimiter=',')
103plt.figure()
104plt.plot(data[:,0], data[:,1+2*nODE2+1], 'b-')
105
106ax=plt.gca() # get current axes
107ax.grid(True, 'major', 'both')
108ax.xaxis.set_major_locator(ticker.MaxNLocator(10)) 
109ax.yaxis.set_major_locator(ticker.MaxNLocator(10)) 
110plt.tight_layout()
111plt.show()