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