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