mainSystemUserFunctionsTest.py
You can view and download this file on Github: mainSystemUserFunctionsTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test of 2 mass-spring-dampers with system-level user functions;
5# these are only tests and should not necessarily replicated, as there are more efficient solutions for some cases - see the comments.
6#
7# Author: Johannes Gerstmayr
8# Date: 2024-10-17
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
15from exudyn.utilities import *
16import exudyn.graphics as graphics
17
18import numpy as np
19
20useGraphics = True #without test
21#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
22#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
23try: #only if called from test suite
24 from modelUnitTests import exudynTestGlobals #for globally storing test results
25 useGraphics = exudynTestGlobals.useGraphics
26except:
27 class ExudynTestGlobals:
28 pass
29 exudynTestGlobals = ExudynTestGlobals()
30#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31useGraphics = False
32
33SC = exu.SystemContainer()
34mbs = SC.AddSystem()
35
36
37caseList = ['SpringDamper','BreakSpring','StopMass','ReplicateSpringDamper']
38result = 0
39
40# case = caseList[3]
41# if True:
42for case in caseList:
43
44 L=1.2 #m
45 mass = 8 #mass in kg
46 spring = 12000 #stiffness of spring-damper in N/m
47 damper = 20 #damping constant in N/(m/s)
48 load = 80 #N
49
50
51 #ground for attachments
52 ground = mbs.CreateGround()
53
54 #mass point
55 mass0 = mbs.CreateMassPoint(referencePosition=[L,0,0],
56 initialDisplacement=[0.,0.5,0],
57 initialVelocity=[50,0.,0.],
58 physicsMass=mass,
59 graphicsDataList=[graphics.Sphere(radius=0.2,
60 color=graphics.color.dodgerblue,
61 nTiles=64)]
62 )
63
64 oCSD = mbs.CreateCartesianSpringDamper(bodyNumbers=[ground, mass0],
65 localPosition0=[L,0,0],
66 localPosition1=[0,0,0],
67 stiffness = [spring]*3,
68 damping=[damper]*3)
69
70 #mass point
71 nMass1 = mbs.AddNode(Node1D(referenceCoordinates=[0], initialCoordinates=[0.2]))
72 mass1 = mbs.AddObject(Mass1D(nodeNumber = nMass1,
73 physicsMass=mass*0.5, referencePosition=[0,L,0],
74 visualization=VMass1D(graphicsData=[graphics.Sphere(radius=0.15,
75 color=graphics.color.green,
76 nTiles=64)])))
77
78
79 #ground node
80 nGround = mbs.AddNode(NodePointGround()) #position is irrelevant
81
82 #marker for ground (=fixed):
83 mncGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
84
85 #marker for springDamper for first (x-)coordinate:
86 mncMass1 =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nMass1, coordinate = 0))
87
88 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
89 if case in ['SpringDamper', 'BreakSpring', 'StopMass']:
90 #Spring-Damper between two marker coordinates
91 #this is what we like to replicate with case ReplicateSpringDamper
92 mbs.AddObject(CoordinateSpringDamper(markerNumbers = [mncGround, mncMass1],
93 stiffness = spring,
94 damping = damper,
95 visualization=VCoordinateSpringDamper(show=False)))
96
97 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
98 if case == 'BreakSpring':
99 #use pre-step function to break spring at certain time
100 #this is ideal for a PreStepUserFunction, because we like to decide at the beginning of step
101 springActive = True
102 def PreStepUserFunction(mbs, t):
103 global springActive, oCSD
104 if t>1 and springActive:
105 #print('break spring')
106 springActive = False
107 mbs.SetObjectParameter(oCSD, 'stiffness', [0,0,0])
108 mbs.SetObjectParameter(oCSD, 'damping', [0,0,0])
109 mbs.SetObjectParameter(oCSD, 'Vshow', False) #this only works in visualization during simulation
110 return True
111
112 mbs.SetPreStepUserFunction(PreStepUserFunction)
113
114 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
115 elif case == 'StopMass':
116 #use post-step function to stop motion of mass at certain time instant
117 #we set the system coordinates of mass0 back to start of step position
118 #this is more like a hack but may be needed in such a way
119 #the steps are computed with motion, but motion is reset just before storing/writing results
120 mbs.Assemble() #to compute LTG-references
121 indMass0 = mbs.systemData.GetObjectLTGODE2(mass0) #LTG-mapping, to know indices in system vector
122 def PostStepUserFunction(mbs, t):
123 global mass0, indMass0
124 if t>0.5:
125 qSys = mbs.systemData.GetODE2Coordinates(configuration=exu.ConfigurationType.Current)
126 qSysStart = mbs.systemData.GetODE2Coordinates(configuration=exu.ConfigurationType.StartOfStep)
127 qSys[indMass0] = qSysStart[indMass0]
128 mbs.systemData.SetODE2Coordinates(qSys, configuration=exu.ConfigurationType.Current)
129 return True
130
131 mbs.SetPostStepUserFunction(PostStepUserFunction)
132
133 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
134 elif case == 'ReplicateSpringDamper':
135 #use system user function in each step to replicate the spring-damper;
136 #don't do such, if not needed: CoordinateSpringDamper or ObjectGenericODE2 would be much more efficient
137 #for loads depending on node coordinates, also consider systemData.AddODE2LoadDependencies, which may be easier!!!
138 #as an example, we consider a load that is changed according to the spring-damper behavior
139 #as this has to be computed in every iteration to be resolved by the Newton method, we have to use the NewtonResidual function
140 #further, we need to adjust jacobian entries, while otherwise, it will not converge or lead to much Newton iterations
141
142 load = mbs.AddLoad(LoadCoordinate(markerNumber=mncMass1, load=0))
143
144 mbs.Assemble() #to compute LTG-references
145 def PreNewtonResidualUserFunction(mbs, t, newtonIt, discontinuousIt):
146 global nMass1, spring, damper
147 u = mbs.GetNodeOutput(nMass1, exu.OutputVariableType.Coordinates) #here, we only need to displacement for the spring
148 v = mbs.GetNodeOutput(nMass1, exu.OutputVariableType.Coordinates_t)
149 mbs.SetLoadParameter(load, 'load', (-spring * u - damper * v)) #at RHS, spring and damper forces have negative signs
150 # print('it'+str(newtonIt)+':',u)
151
152 mbs.SetPreNewtonResidualUserFunction(PreNewtonResidualUserFunction)
153
154 indNode1 = mbs.GetNodeODE2Index(nMass1) #global index of node1 in system vector
155
156 #add terms for RHS-jacobian, then total number of Newton iterations is 2000, otherwise 5511 for 2 seconds
157 #NOTE: currently, this is not called for initial accelerations, so you should set computeInitialAccelerations=False
158 def SystemJacobianUserFunction(mbs, t, factorODE2, factorODE2_t, factorODE1):
159 val = -factorODE2 * spring - factorODE2_t * damper #negative because of RHS
160 mc = exu.MatrixContainer()
161 mc.SetWithSparseMatrix([[indNode1, indNode1, val]],3,3)
162 return mc
163
164 mbs.SetSystemJacobianUserFunction(SystemJacobianUserFunction)
165
166 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
167
168 sPos0 = mbs.AddSensor(SensorBody(bodyNumber = mass0, storeInternal=True,
169 outputVariableType=exu.OutputVariableType.Displacement))
170 sPos1 = mbs.AddSensor(SensorBody(bodyNumber = mass1, storeInternal=True,
171 outputVariableType=exu.OutputVariableType.Displacement))
172
173 mbs.Assemble()
174
175 tEnd = 2
176 stepSize = 0.001
177
178 simulationSettings = exu.SimulationSettings()
179 #simulationSettings.solutionSettings.solutionWritePeriod = 2e-3 #output interval
180 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
181 simulationSettings.timeIntegration.endTime = tEnd
182 simulationSettings.solutionSettings.solutionInformation = 'CASE: '+case
183 simulationSettings.displayStatistics = True
184 simulationSettings.timeIntegration.verboseMode = 1
185
186 simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations=False
187
188 SC.visualizationSettings.general.drawWorldBasis = True
189
190
191 if useGraphics:
192 exu.StartRenderer() #start graphics visualization
193 mbs.WaitForUserToContinue() #wait for pressing SPACE bar to continue
194
195 #start solver:
196 mbs.SolveDynamic(simulationSettings)
197
198 if useGraphics:
199 SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit
200 exu.StopRenderer() #safely close rendering window!
201
202 #evaluate final (=current) output values
203 u0 = mbs.GetSensorValues(sPos0)
204 u1 = mbs.GetSensorValues(sPos1)
205 exu.Print('displacement=',u0[0], u0[1], u1[0])
206
207 # for cases 'SpringDamper','ReplicateSpringDamper':
208 # displacement= 0.0959942160568712 -0.016318265957667506 -0.001118733080513271
209 # for case 'BreakSpring':
210 # displacement= 7.840727993539164 -4.554866482139987 -0.001118733080513271
211 # for case 'StopMass':
212 # displacement= 0.3099518908156229 0.24354903131196867 -0.001118733080513271
213
214 result += u0[0]+u0[1]+u1[0]
215
216exu.Print('mainSystemUserFunctionsTest solution=', result)
217exudynTestGlobals.testResult = result
218
219
220if useGraphics:
221 mbs.SolutionViewer()
222
223#+++++++++++++++++++++++++++++++++++++++++++++++++++++
224
225if useGraphics:
226 mbs.PlotSensor(sPos0, components=[0], closeAll=True)
227 mbs.PlotSensor(sPos1, components=[0], newFigure=False)