plotSensorExamples.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  This example serves as demonstration for PlotSensor
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2022-02-19
  8#
  9# 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.
 10#
 11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 12
 13import exudyn as exu
 14from exudyn.itemInterface import *
 15from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 16import exudyn.graphics as graphics #only import if it does not conflict
 17
 18import numpy as np #for postprocessing
 19from math import pi
 20
 21L=0.5
 22mass = 1.6          #mass in kg
 23spring = 4000       #stiffness of spring-damper in N/m
 24damper = 8          #damping constant in N/(m/s)
 25
 26u0=-0.08            #initial displacement
 27v0=1                #initial velocity
 28f =80               #force on mass
 29x0=f/spring         #static displacement
 30
 31SC = exu.SystemContainer()
 32mbs = SC.AddSystem()
 33
 34#node for 3D mass point:
 35nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
 36
 37#add rigid body for sensor tests:
 38iCube0 = InertiaCuboid(density=5000, sideLengths=[0.2,0.1,0.5])
 39iCube0 = iCube0.Translated([0.1,0.2,0.3])
 40dictCube0 = mbs.CreateRigidBody(
 41              inertia=iCube0,  # includes COM
 42              nodeType=exu.NodeType.RotationRxyz,
 43              initialAngularVelocity=[4, 0.1, 0.1],
 44              returnDict=True)
 45[n0, b0] = [dictCube0['nodeNumber'], dictCube0['bodyNumber']]
 46
 47#add spring damper system
 48n1=mbs.AddNode(NodePoint(referenceCoordinates = [L,0,0],
 49             initialCoordinates = [u0,0,0],
 50             initialVelocities= [v0,0,0]))
 51
 52
 53#add mass point (this is a 3D object with 3 coordinates):
 54massPoint = mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = n1))
 55
 56#marker for ground (=fixed):
 57groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
 58#marker for springDamper for first (x-)coordinate:
 59nodeMarker  =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0))
 60
 61#spring-damper between two marker coordinates
 62nC = mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker],
 63                                          stiffness = spring, damping = damper))
 64
 65
 66#add load:
 67mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker,
 68                                          load = f))
 69
 70#add sensor:
 71sForce = mbs.AddSensor(SensorObject(objectNumber=nC,
 72                           storeInternal=True,
 73                           outputVariableType=exu.OutputVariableType.Force))
 74
 75sDisp = mbs.AddSensor(SensorNode(nodeNumber=n1, storeInternal=True, fileName='solution/sDisp.txt',
 76                           outputVariableType=exu.OutputVariableType.Displacement))
 77sVel = mbs.AddSensor(SensorNode(nodeNumber=n1, storeInternal=True,
 78                           outputVariableType=exu.OutputVariableType.Velocity))
 79
 80sOmega = mbs.AddSensor(SensorNode(nodeNumber=n0, storeInternal=True,
 81                           outputVariableType=exu.OutputVariableType.AngularVelocity))
 82sPos = mbs.AddSensor(SensorNode(nodeNumber=n0, storeInternal=True,
 83                           outputVariableType=exu.OutputVariableType.Position))
 84sRot = mbs.AddSensor(SensorNode(nodeNumber=n0, storeInternal=True,
 85                           outputVariableType=exu.OutputVariableType.Rotation))
 86#dummy sensor, writes only zeros
 87sDummy= mbs.AddSensor(SensorNode(nodeNumber=nGround, storeInternal=True,
 88                           outputVariableType=exu.OutputVariableType.Displacement))
 89
 90#%%++++++++++++++++++++
 91mbs.Assemble()
 92
 93tEnd = 4     #end time of simulation
 94h = 0.002    #step size; leads to 1000 steps
 95
 96simulationSettings = exu.SimulationSettings()
 97simulationSettings.solutionSettings.solutionWritePeriod = 0.005  #output interval general
 98simulationSettings.solutionSettings.writeSolutionToFile = False
 99simulationSettings.solutionSettings.sensorsWritePeriod = 1*h  #output interval of sensors
100
101simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h) #must be integer
102simulationSettings.timeIntegration.endTime = tEnd
103simulationSettings.timeIntegration.verboseMode = 1
104simulationSettings.displayComputationTime = True
105
106simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
107
108# exu.StartRenderer()              #start graphics visualization
109#mbs.WaitForUserToContinue()    #wait for pressing SPACE bar to continue
110
111#start solver:
112mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
113dispExplicit=mbs.GetSensorStoredData(sDisp)
114velExplicit=mbs.GetSensorStoredData(sVel)
115omegaExplicit=mbs.GetSensorStoredData(sOmega)
116
117mbs.SolveDynamic(simulationSettings)#, solverType=exu.DynamicSolverType.ExplicitEuler)
118
119#SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit
120# exu.StopRenderer()               #safely close rendering window!
121
122#evaluate final (=current) output values
123u = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
124print('displacement=',u)
125
126# data=mbs.GetSensorStoredData(0)
127# print('sensor data=',data)
128
129
130
131
132# import matplotlib.pyplot as plt
133mbs.PlotSensor(sensorNumbers=sDisp, components=0, closeAll=True)
134
135mbs.PlotSensor(sVel, 0) #SIMPLEST command to plot x-coordinate of velocity sensor
136
137#compare difference of sensors:
138mbs.PlotSensor(sensorNumbers=sVel, components=0, newFigure=False, colorCodeOffset=1,
139            offsets=[-velExplicit], labels='difference of velocity \nof expl./impl. integrator')
140
141mbs.PlotSensor(sensorNumbers=sForce, components=0, newFigure=False, factors=[1e-3], colorCodeOffset=2)
142
143#internal data and file names; compute difference to external data:
144extData = np.loadtxt('solution/sDisp.txt', comments='#', delimiter=',')
145mbs.PlotSensor(sensorNumbers=['solution/sDisp.txt',sDisp,sDisp], components=0, xLabel='time in seconds',
146            offsets=[0,0,-extData],
147            markerStyles=['','x',''], lineStyles=['-','','-'], markerDensity=0.05,
148            labels=['Displacement from file','Displacement internal','diff between file and \ninternal data (precision)'])
149
150mbs.PlotSensor(sensorNumbers=sOmega, components=[0,1,2],
151          yLabel='angular velocities with offset 0\nand scaled with $\\frac{180}{\pi}$',
152          factors=180/pi, offsets=0,fontSize=12,title='angular velocities',
153          lineWidths=[3,5,1], lineStyles=['-',':','-.'], colors=['r','g','b'])
154
155mbs.PlotSensor(sensorNumbers=[sRot]*3+[sOmega]*3, components=[0,1,2]*2,
156          colorCodeOffset=3, newFigure=True, fontSize=14,
157          yLabel='Tait-Bryan rotations $\\alpha, \\beta, \\gamma$ and\n angular velocities around $x,y,z$',
158          title='compare rotations and angular velocities')
159
160mbs.PlotSensor(sensorNumbers=sRot, components=[0,1,2], markerStyles=['* ','x','^ '], #add space after marker symbol to draw empty
161            lineWidths=2, markerSizes=12, markerDensity=15)
162
163
164#create subplots:
165subs=[3,2]
166mbs.PlotSensor(sensorNumbers=sOmega, components=0, newFigure=True,  subPlot=[*subs,1])
167mbs.PlotSensor(sensorNumbers=sOmega, components=1, newFigure=False, subPlot=[*subs,2])
168mbs.PlotSensor(sensorNumbers=sOmega, components=2, newFigure=False, subPlot=[*subs,3])
169mbs.PlotSensor(sensorNumbers=sPos,   components=0, newFigure=False, subPlot=[*subs,4])
170mbs.PlotSensor(sensorNumbers=sPos,   components=1, newFigure=False, subPlot=[*subs,5])
171mbs.PlotSensor(sensorNumbers=sPos,   components=2, newFigure=False, subPlot=[*subs,6])
172
173#compare different simulation results (could also be done with stored files ...):
174omegaImplicit=mbs.GetSensorStoredData(sOmega)
175mbs.PlotSensor(sensorNumbers=[sOmega,sOmega], components=[0,0], newFigure=True,  subPlot=[1,3,1],
176           offsets=[0.,omegaExplicit-omegaImplicit], sizeInches=[12,4], labels=['omegaX impl.','omegaX expl.'])
177mbs.PlotSensor(sensorNumbers=[sOmega,sOmega], components=[1,1], newFigure=False, subPlot=[1,3,2],
178           offsets=[0.,omegaExplicit-omegaImplicit], sizeInches=[12,4], labels=['omegaX impl.','omegaX expl.'])
179mbs.PlotSensor(sensorNumbers=[sOmega,sOmega], components=[2,2], newFigure=False, subPlot=[1,3,3],
180           offsets=[0.,omegaExplicit-omegaImplicit], sizeInches=[12,4], labels=['omegaY impl.','omegaY expl.'],
181           fileName='solution/fig_omega.pdf')
182
183
184#PHASE Plot, more complicated ...; using dummy sensor with zero values
185data = 0.*mbs.GetSensorStoredData(sDisp) #create data set
186data[:,1] = mbs.GetSensorStoredData(sDisp)[:,1] #x
187data[:,2] = mbs.GetSensorStoredData(sVel)[:,1]  #y
188mbs.PlotSensor(sensorNumbers=[sDummy], componentsX=[0], components=[1], xLabel='Position', yLabel='Velocity',
189           offsets=[data], labels='velocity over displacement', title='Phase plot',
190           rangeX=[-0.01,0.04],rangeY=[-1,1], majorTicksX=6, majorTicksY=6)
191
192##plot y over x:
193#mbs.PlotSensor(sensorNumbers=s0, componentsX=[0], components=[1], xLabel='x-Position', yLabel='y-Position')