objectFFRFTest2.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for ObjectFFRF with C++ implementation user function for reduced order equations of motion
  5# NOTE: this is a development file, with lots of unstructured code; just kept for consistency!
  6#
  7# Author:   Johannes Gerstmayr
  8# Date:     2020-05-13
  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 *
 16from exudyn.FEM import *
 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#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 31
 32SC = exu.SystemContainer()
 33mbs = SC.AddSystem()
 34
 35
 36
 37#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 38#Use FEMinterface to import FEM model and create FFRFreducedOrder object
 39fem = FEMinterface()
 40inputFileName = 'testData/rotorDiscTest' #runTestSuite.py is at another directory
 41#if useGraphics:
 42#    inputFileName = 'testData/rotorDiscTest'        #if executed in current directory
 43
 44nodes=fem.ImportFromAbaqusInputFile(inputFileName+'.inp', typeName='Instance', name='rotor-1')
 45
 46fem.ReadMassMatrixFromAbaqus(inputFileName+'MASS1.mtx')
 47fem.ReadStiffnessMatrixFromAbaqus(inputFileName+'STIF1.mtx')
 48fem.ScaleStiffnessMatrix(1e-2) #for larger deformations, stiffness is reduced to 1%
 49
 50nodeNumberUnbalance = 9  #on disc, max y-value
 51unbalance = 0.1
 52fem.AddNodeMass(nodeNumberUnbalance, unbalance)
 53#print(fem.GetMassMatrix()[8*3:11*3,:])
 54
 55nModes = 8
 56fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
 57#print("eigen freq.=", fem.GetEigenFrequenciesHz())
 58
 59ffrf = ObjectFFRFinterface(fem)
 60
 61##user functions should be defined outside of class:
 62#def UFmassFFRF(t, qReduced, qReduced_t):
 63#    return cms.UFmassFFRF(exu, mbs, t, qReduced, qReduced_t)
 64#
 65#def UFforceFFRF(t, qReduced, qReduced_t):
 66#    return cms.UFforceFFRF(exu, mbs, t, qReduced, qReduced_t)
 67
 68objFFRF = ffrf.AddObjectFFRF(exu, mbs, positionRef=[0,0,0], eulerParametersRef=eulerParameters0,
 69                             initialVelocity=[0,0,0], initialAngularVelocity=[0,0,50*2*pi],
 70                             gravity = [0,-0*9.81,0],
 71                             #UFforce=UFforceFFRFreducedOrder, UFmassMatrix=UFmassFFRFreducedOrder,
 72                             color=[0.1,0.9,0.1,1.])
 73
 74#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 75#add markers and joints
 76nodeDrawSize = 0.0025 #for joint drawing
 77
 78pLeft = [0,0,0]
 79pRight = [0,0,0.5]
 80nMid = fem.GetNodeAtPoint([0,0,0.25])
 81#print("nMid=",nMid)
 82
 83mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
 84oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
 85
 86mGroundPosLeft = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pLeft))
 87mGroundPosRight = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pRight))
 88
 89#torque on reference frame:
 90#mbs.AddLoad(Torque(markerNumber=mRB, loadVector=[0,0,100*2*pi]))
 91
 92
 93#++++++++++++++++++++++++++++++++++++++++++
 94#find nodes at left and right surface:
 95nodeListLeft = fem.GetNodesInPlane(pLeft, [0,0,1])
 96nodeListRight = fem.GetNodesInPlane(pRight, [0,0,1])
 97#nLeft = fem.GetNodeAtPoint(pLeft)
 98#nRight = fem.GetNodeAtPoint(pRight)
 99
100
101lenLeft = len(nodeListLeft)
102lenRight = len(nodeListRight)
103weightsLeft = np.array((1./lenLeft)*np.ones(lenLeft))
104weightsRight = np.array((1./lenRight)*np.ones(lenRight))
105
106addSupports = True
107if addSupports:
108    k = 2e8     #joint stiffness
109    d = k*0.01  #joint damping
110
111    useSpringDamper = True
112
113    mLeft = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=objFFRF['oFFRF'],
114                                                    meshNodeNumbers=np.array(nodeListLeft), #these are the meshNodeNumbers
115                                                    weightingFactors=weightsLeft))
116    mRight = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=objFFRF['oFFRF'],
117                                                    meshNodeNumbers=np.array(nodeListRight), #these are the meshNodeNumbers
118                                                    weightingFactors=weightsRight))
119    if useSpringDamper:
120        oSJleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mLeft, mGroundPosLeft],
121                                            stiffness=[k,k,k], damping=[d,d,d]))
122        oSJright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mRight,mGroundPosRight],
123                                            stiffness=[k,k,0], damping=[d,d,d]))
124    else:
125        oSJleft = mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosLeft,mLeft], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
126        oSJright= mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosRight,mRight], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
127
128
129fileDir = 'solution/'
130#keep files, as they are checked in the .git repo:
131sDisp=mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRF'], meshNodeNumber=nMid, #meshnode number!
132                         storeInternal=True,#fileName=fileDir+'nMidDisplacementFFRFtest.txt',
133                         outputVariableType = exu.OutputVariableType.Displacement))
134
135sAngVel=mbs.AddSensor(SensorNode(nodeNumber=objFFRF['nRigidBody'],
136                         storeInternal=True,#fileName=fileDir+'nRigidBodyAngVelFFRFtest.txt',
137                         outputVariableType = exu.OutputVariableType.AngularVelocity))
138
139mbs.Assemble()
140
141simulationSettings = exu.SimulationSettings()
142
143SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
144SC.visualizationSettings.nodes.drawNodesAsPoint = False
145SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
146
147SC.visualizationSettings.nodes.show = True
148SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
149SC.visualizationSettings.nodes.basisSize = 0.12
150SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
151
152SC.visualizationSettings.openGL.showFaceEdges = True
153SC.visualizationSettings.openGL.showFaces = True
154
155SC.visualizationSettings.sensors.show = True
156SC.visualizationSettings.sensors.drawSimplified = False
157SC.visualizationSettings.sensors.defaultSize = 0.01
158SC.visualizationSettings.markers.drawSimplified = False
159SC.visualizationSettings.markers.show = True
160SC.visualizationSettings.markers.defaultSize = 0.01
161
162SC.visualizationSettings.loads.drawSimplified = False
163
164SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
165SC.visualizationSettings.contour.outputVariableComponent = 1 #y-component
166
167simulationSettings.solutionSettings.solutionInformation = "ObjectFFRF test"
168simulationSettings.solutionSettings.writeSolutionToFile=False
169
170h=1e-4
171tEnd = 0.0025
172if useGraphics:
173    tEnd = 0.0025
174
175simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
176simulationSettings.timeIntegration.endTime = tEnd
177simulationSettings.solutionSettings.solutionWritePeriod = h
178simulationSettings.timeIntegration.verboseMode = 1
179#simulationSettings.timeIntegration.verboseModeFile = 3
180simulationSettings.timeIntegration.newton.useModifiedNewton = True
181
182simulationSettings.solutionSettings.sensorsWritePeriod = h
183simulationSettings.solutionSettings.coordinatesSolutionFileName = "solution/coordinatesSolutionFFRFtest.txt"
184simulationSettings.solutionSettings.writeSolutionToFile=False
185
186simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5 #SHOULD work with 0.9 as well
187#simulationSettings.displayStatistics = True
188#simulationSettings.displayComputationTime = True
189
190#create animation:
191#simulationSettings.solutionSettings.recordImagesInterval = 0.0002
192#SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
193
194if useGraphics:
195    exu.StartRenderer()
196    if 'lastRenderState' in vars():
197        SC.SetRenderState(lastRenderState) #load last model view
198
199    mbs.WaitForUserToContinue() #press space to continue
200
201mbs.SolveDynamic(simulationSettings)
202
203
204#data = np.loadtxt(fileDir+'nMidDisplacementFFRFtest.txt', comments='#', delimiter=',')
205data = mbs.GetSensorStoredData(sDisp)
206result = abs(data).sum()
207#pos = mbs.GetObjectOutputBody(objFFRF['oFFRFreducedOrder'],exu.OutputVariableType.Position, localPosition=[0,0,0])
208exu.Print('solution of ObjectFFRFtest2=',result)
209
210exudynTestGlobals.testError = result - (0.03552188069017914) #2022-02-20: changed to internal sensor storage; 2020-05-26 (tEnd=0.0025, h=1e-4): 0.03553746369388042
211exudynTestGlobals.testResult = result
212
213if useGraphics:
214    SC.WaitForRenderEngineStopFlag()
215    exu.StopRenderer() #safely close rendering window!
216    lastRenderState = SC.GetRenderState() #store model view for next simulation
217
218##++++++++++++++++++++++++++++++++++++++++++++++q+++++++
219#plot results
220if useGraphics:
221
222
223    mbs.PlotSensor([fileDir+'nMidDisplacementCMS8.txt',sDisp], components=1, closeAll=True)
224
225    # import matplotlib.pyplot as plt
226    # import matplotlib.ticker as ticker
227    # cList=['r-','g-','b-','k-','c-','r:','g:','b:','k:','c:']
228
229    # data = np.loadtxt(fileDir+'nMidDisplacementCMS8.txt', comments='#', delimiter=',') #new result from this file
230    # plt.plot(data[:,0], data[:,2], cList[1],label='uMid,CMS8') #numerical solution, 1 == x-direction
231
232    # data = np.loadtxt(fileDir+'nMidDisplacementFFRFtest.txt', comments='#', delimiter=',')
233    # plt.plot(data[:,0], data[:,2], cList[2],label='uMid,FFRF') #numerical solution, 1 == x-direction
234
235    # ax=plt.gca() # get current axes
236    # ax.grid(True, 'major', 'both')
237    # ax.xaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
238    # ax.yaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
239    # plt.tight_layout()
240    # plt.legend()
241    # plt.show()