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