genericODE2test.py

You can view and download this file on Github: genericODE2test.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test model for GenericODE2
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2020-03-09
  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.utilities import * #includes itemInterface and rigidBodyUtilities
 15import exudyn.graphics as graphics #only import if it does not conflict
 16
 17import numpy as np
 18
 19useGraphics = True #without test
 20#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 21#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 22try: #only if called from test suite
 23    from modelUnitTests import exudynTestGlobals #for globally storing test results
 24    useGraphics = exudynTestGlobals.useGraphics
 25except:
 26    class ExudynTestGlobals:
 27        pass
 28    exudynTestGlobals = ExudynTestGlobals()
 29#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 30
 31SC = exu.SystemContainer()
 32mbs = SC.AddSystem()
 33
 34
 35#linear spring-damper
 36L=0.5
 37mass = 1.6
 38k = 4000
 39omega0 = 50 # sqrt(4000/1.6)
 40dRel = 0.05
 41d = dRel * 2 * 80 #80=sqrt(1.6*4000)
 42u0=-0.08
 43v0=1
 44f = 80
 45#static result = f/k = 0.01;
 46x0 = f/k
 47#(exact) dynamic result with u0 and v0: (see bottom of file)
 48
 49#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 50#node for mass point:
 51n1=mbs.AddNode(Point(referenceCoordinates = [L,0,0], initialCoordinates = [u0,0,0], initialVelocities= [v0,0,0]))
 52
 53#add mass points and ground object:
 54objectGround = mbs.AddObject(ObjectGround(referencePosition = [0,0,0]))
 55massPoint1 = mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = n1))
 56
 57#marker for constraint / springDamper
 58groundMarker = mbs.AddMarker(MarkerBodyPosition(bodyNumber = objectGround, localPosition= [0, 0, 0]))
 59bodyMarker = mbs.AddMarker(MarkerBodyPosition(bodyNumber = massPoint1, localPosition= [0, 0, 0]))
 60
 61mbs.AddObject(CartesianSpringDamper(markerNumbers = [groundMarker, bodyMarker], stiffness = [k,k,k], damping = [d,0,0], offset = [L,0,0]))
 62
 63#add loads:
 64mbs.AddLoad(Force(markerNumber = bodyMarker, loadVector = [f, 0, 0]))
 65
 66
 67#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 68#node for mass point:
 69
 70n2=mbs.AddNode(Point(referenceCoordinates = [2*L,0,0], initialCoordinates = [u0,f/k*0.999,0], initialVelocities= [v0,0,0]))
 71
 72M=np.diag([mass,mass,mass])
 73#exu.Print("M =",M)
 74K=np.diag([k,k,k])
 75#exu.Print("K =",K)
 76D=np.diag([d,0,d])
 77#exu.Print("D =",D)
 78fv=np.array([f,f,0])
 79#exu.Print("fv =",fv)
 80
 81fdyn=np.array([0,0,10])
 82
 83def Sweep(t, t1, f0, f1):
 84    k = (f1-f0)/t1
 85    return np.sin(2*np.pi*(f0+k*0.5*t)*t) #take care of factor 0.5 in k*0.5*t, in order to obtain correct frequencies!!!
 86
 87def UFgenericODE2(mbs, t, itemIndex, q, q_t):
 88    #f = np.sin(t*2*np.pi*10)*fdyn
 89    f = Sweep(t,10,1,100)*fdyn
 90    return f
 91    #exu.Print("t =", t, ", f =", f)
 92
 93def UFmassGenericODE2(mbs, t, itemIndex, q, q_t):
 94    return 1*M
 95
 96mbs.AddObject(ObjectGenericODE2(nodeNumbers = [n2], massMatrix=M, stiffnessMatrix=K, dampingMatrix=D, forceVector=fv,
 97                                forceUserFunction=UFgenericODE2, massMatrixUserFunction=UFmassGenericODE2))
 98
 99#exu.Print(mbs)
100mbs.Assemble()
101
102simulationSettings = exu.SimulationSettings()
103
104tEnd = 1
105steps = 2000
106simulationSettings.timeIntegration.numberOfSteps = steps
107simulationSettings.timeIntegration.endTime = tEnd
108simulationSettings.solutionSettings.solutionWritePeriod = tEnd/steps
109simulationSettings.timeIntegration.verboseMode = 1
110#simulationSettings.solutionSettings.solutionWritePeriod = tEnd/steps
111
112simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #SHOULD work with 0.9 as well
113
114mbs.SolveDynamic(simulationSettings)
115
116u1 = mbs.GetNodeOutput(n1, exu.OutputVariableType.Coordinates)
117#exu.Print("u1 =", u1)
118u2 = mbs.GetNodeOutput(n2, exu.OutputVariableType.Coordinates)
119#exu.Print("u2 =", u2)
120
121u=NormL2(u1) + NormL2(u2)
122exu.Print('solution of genericODE2test=',u)
123
124exudynTestGlobals.testError = u - (0.03604546349898683) #2020-04-22: 0.03604546349898683
125exudynTestGlobals.testResult = u