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