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')