ObjectFFRFconvergenceTestHinge.py
You can view and download this file on Github: ObjectFFRFconvergenceTestHinge.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test for Hurty-Craig-Bampton modes using a simple flexible pendulum meshed with Netgen
5#
6# Author: Johannes Gerstmayr
7# Date: 2021-04-20
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
13
14import exudyn as exu
15#a local run must not open a plot window when windows are suppressed (#2477)
16if exu.special.userInterface.suppressPlots: import matplotlib; matplotlib.use('Agg')
17from exudyn.itemInterface import *
18from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
19import exudyn.graphics as graphics
20from exudyn.FEM import *
21import time
22
23SC = exu.SystemContainer()
24mbs = SC.AddSystem()
25
26import numpy as np
27
28#import timeit
29
30import exudyn.basicUtilities as eb
31import exudyn.rigidBodyUtilities as rb
32import exudyn.utilities as eu
33
34import numpy as np
35
36useGraphics = True
37fileName = OutputFilePath('solution/netgenHinge', 'ObjectFFRFconvergenceTestHinge') #generated output, not input: goes to the ignored solution/ (#2491)
38
39#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
40#netgen/meshing part:
41fem = FEMinterface()
42
43#geometrical parameters:
44L = 0.4 #Length of plate (X)
45w = 0.04 #width of plate (Y)
46h = 0.02 #height of plate (Z)
47d = 0.03 #diameter of bolt
48D = d*2 #diameter of bushing
49b = 0.05 #length of bolt
50nModes = 32 #128
51meshH = 0.01 #0.01 is default, 0.002 gives 100000 nodes and is fairly converged
52#meshH = 0.0014 #203443 nodes, takes 1540 seconds for eigenmode computation (free-free) and 753 seconds for postprocessing on i9
53
54#steel:
55rho = 7850
56Emodulus=2.1e11
57nu=0.3
58
59#test high flexibility
60Emodulus=2e8
61# nModes = 32
62
63
64#helper function for cylinder with netgen
65def CSGcylinder(p0,p1,r):
66 v = (np.array(p1) - p0)
67 v = Normalize(v)
68 cyl = Cylinder(Pnt(p0[0],p0[1],p0[2]), Pnt(p1[0],p1[1],p1[2]),
69 r) * Plane(Pnt(p0[0],p0[1],p0[2]), Vec(-v[0],-v[1],-v[2])) * Plane(Pnt(p1[0],p1[1],p1[2]), Vec(v[0],v[1],v[2]))
70 return cyl
71
72mBushing = None
73meshCreated = False
74
75#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
76if True: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
77
78 import ngsolve as ngs
79 import netgen
80 from netgen.meshing import *
81
82 from netgen.geom2d import unit_square
83 #import netgen.libngpy as libng
84 from netgen.csg import *
85
86 geo = CSGeometry()
87
88 #plate
89 block = OrthoBrick(Pnt(0, 0, -0.5*h),Pnt(L, w, 0.5*h))
90
91 #bolt
92 bolt0 = CSGcylinder(p0=[0,w,0], p1=[0,0,0], r=1.6*h)
93 bolt = CSGcylinder(p0=[0,0.5*w,0], p1=[0,-b,0], r=0.5*d)
94
95 #bushing
96 bushing = (CSGcylinder(p0=[L,w,0], p1=[L,-b,0], r=0.5*D) -
97 CSGcylinder(p0=[L,0,0], p1=[L,-b*1.1,0], r=0.5*d))
98
99 geo.Add(block+bolt0+bolt+bushing)
100
101 mesh = ngs.Mesh( geo.GenerateMesh(maxh=meshH))
102 mesh.Curve(1)
103
104 if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
105 # import netgen
106 import netgen.gui
107 ngs.Draw(mesh)
108 netgen.Redraw()
109
110 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
111 #Use fem to import FEM model and create FFRFreducedOrder object
112 fem.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
113 meshCreated = True
114 if (meshH==0.01): fem.SaveToFile(fileName)
115
116#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
117#compute Hurty-Craig-Bampton modes
118if True: #now import mesh as mechanical model to EXUDYN
119 if not meshCreated: fem.LoadFromFile(fileName)
120
121 boltP1=[0,0,0]
122 boltP2=[0,-b,0]
123 nodesOnBolt = fem.GetNodesOnCylinder(boltP1, boltP2, radius=0.5*d)
124 #print("boundary nodes bolt=", nodesOnBolt)
125 nodesOnBoltLen = len(nodesOnBolt)
126 nodesOnBoltWeights = np.array((1./nodesOnBoltLen)*np.ones(nodesOnBoltLen))
127
128 bushingP1=[L,0,0]
129 bushingP2=[L,-b,0]
130 nodesOnBushing = fem.GetNodesOnCylinder(bushingP1, bushingP2, radius=0.5*d)
131 #print("boundary nodes bushing=", nodesOnBushing)
132 nodesOnBushingLen = len(nodesOnBushing)
133 nodesOnBushingWeights = np.array((1./nodesOnBushingLen)*np.ones(nodesOnBushingLen))
134
135 print("nNodes=",fem.NumberOfNodes())
136
137 strMode = ''
138 if False: #pure eigenmodes
139 print("compute eigen modes... ")
140 start_time = time.time()
141 fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
142 print("eigen modes computation needed %.3f seconds" % (time.time() - start_time))
143 print("eigen freq.=", fem.GetEigenFrequenciesHz())
144
145 elif False:
146 strMode = 'HCB'
147 #boundaryList = [nodesOnBolt, nodesOnBolt, nodesOnBushing] #for visualization, use first interface twice
148 boundaryList = [nodesOnBolt, nodesOnBushing]
149
150 print("compute HCB modes... ")
151 start_time = time.time()
152 fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
153 nEigenModes=nModes,
154 useSparseSolver=True,
155 computationMode = HCBstaticModeSelection.RBE2)
156
157 print("eigen freq.=", fem.GetEigenFrequenciesHz())
158 print("HCB modes needed %.3f seconds" % (time.time() - start_time))
159 else:
160 strMode = 'HCBsingle'
161 #boundaryList = [nodesOnBolt, nodesOnBolt, nodesOnBushing] #for visualization, use first interface twice
162 boundaryList = [nodesOnBolt]
163
164 print("compute HCB single modes... ")
165 start_time = time.time()
166 fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
167 nEigenModes=nModes,
168 useSparseSolver=True,
169 computationMode = HCBstaticModeSelection.RBE2)
170
171 print("eigen freq.=", fem.GetEigenFrequenciesHz())
172 print("HCB modes needed %.3f seconds" % (time.time() - start_time))
173
174
175
176 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
177 #compute stress modes for postprocessing (inaccurate for coarse meshes, just for visualization):
178 if False:
179 mat = KirchhoffMaterial(Emodulus, nu, rho)
180 varType = exu.OutputVariableType.StressLocal
181 #varType = exu.OutputVariableType.StrainLocal
182 print("ComputePostProcessingModes ... (may take a while)")
183 start_time = time.time()
184 fem.ComputePostProcessingModes(material=mat,
185 outputVariableType=varType)
186 print(" ... needed %.3f seconds" % (time.time() - start_time))
187 SC.visualizationSettings.contour.reduceRange=False
188 SC.visualizationSettings.contour.outputVariable = varType
189 SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
190 else:
191 SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
192 SC.visualizationSettings.contour.outputVariableComponent = -1
193
194 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
195 print("create CMS element ...")
196 cms = ObjectFFRFreducedOrderInterface(fem)
197
198 objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
199 initialVelocity=[0,0,0],
200 initialAngularVelocity=[0,0,0],
201 color=[0.9,0.9,0.9,1.],
202 )
203
204 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
205 #add markers and joints
206 nodeDrawSize = 0.0025 #for joint drawing
207
208
209 #mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
210
211 if False:
212 boltMidPoint = 0.5*(np.array(boltP1)+boltP2)
213
214 oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
215
216 altApproach = True
217 mBolt = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
218 meshNodeNumbers=np.array(nodesOnBolt), #these are the meshNodeNumbers
219 #referencePosition=boltMidPoint,
220 useAlternativeApproach=altApproach,
221 weightingFactors=nodesOnBoltWeights))
222 bushingMidPoint = 0.5*(np.array(bushingP1)+bushingP2)
223
224 #add marker for visualization of boundary nodes
225 mBushing = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
226 meshNodeNumbers=np.array(nodesOnBushing), #these are the meshNodeNumbers
227 #referencePosition=bushingMidPoint,
228 useAlternativeApproach=altApproach,
229 weightingFactors=nodesOnBushingWeights))
230
231 lockedAxes=[1,1,1,1,1*0,1]
232 if True:
233
234 mGroundBolt = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround,
235 localPosition=boltMidPoint,
236 visualization=VMarkerBodyRigid(show=True)))
237 mbs.AddObject(GenericJoint(markerNumbers=[mGroundBolt, mBolt],
238 constrainedAxes = lockedAxes,
239 visualization=VGenericJoint(show=False, axesRadius=0.1*b, axesLength=0.1*b)))
240
241 else:
242
243 mGroundBushing = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=bushingMidPoint))
244 mbs.AddObject(GenericJoint(markerNumbers=[mGroundBushing, mBushing],
245 constrainedAxes = lockedAxes,
246 visualization=VGenericJoint(axesRadius=0.1*b, axesLength=0.1*b)))
247
248
249 if False:
250 cms = ObjectFFRFreducedOrderInterface(fem)
251
252 objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
253 initialVelocity=[990,990,990],
254 initialAngularVelocity=[0,0,0],
255 color=[0.9,0.9,0.9,1.],
256 )
257
258 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
259 #animate modes
260 SC.visualizationSettings.markers.show = True
261 SC.visualizationSettings.markers.defaultSize=0.0075
262 SC.visualizationSettings.markers.drawSimplified = False
263
264 SC.visualizationSettings.loads.show = False
265 SC.visualizationSettings.loads.drawSimplified = False
266 SC.visualizationSettings.loads.defaultSize=0.1
267 SC.visualizationSettings.loads.defaultRadius = 0.002
268
269 SC.visualizationSettings.openGL.multiSampling=4
270 SC.visualizationSettings.openGL.lineWidth=2
271
272 if False: #activate to animate modes
273 from exudyn.interactive import AnimateModes
274 mbs.Assemble()
275 SC.visualizationSettings.nodes.show = False
276 SC.visualizationSettings.view0.scene.showFaceEdges = True
277 SC.visualizationSettings.openGL.multiSampling=4
278 SC.visualizationSettings.openGL.lineWidth=2
279 SC.visualizationSettings.view0.window.renderWindowSize = [1600,1080]
280 SC.visualizationSettings.contour.advanced.showColorBar = False
281 SC.visualizationSettings.view0.window.globalFontSize = 16
282
283 #%%+++++++++++++++++++++++++++++++++++++++
284 #animate modes of ObjectFFRFreducedOrder (only needs generic node containing modal coordinates)
285 SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved
286
287 nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
288 AnimateModes(SC, mbs, nodeNumber, period=0.1, showTime=False, renderWindowText='Hurty-Craig-Bampton: 2 x 6 static modes and 8 eigenmodes\n')
289 # import sys
290 # sys.exit()
291
292 #add gravity (not necessary if user functions used)
293 oFFRF = objFFRF['oFFRFreducedOrder']
294 mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber=oFFRF))
295 mbs.AddLoad(LoadMassProportional(markerNumber=mBody, loadVector= [0,0,-9.81*0]))
296
297
298 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
299 if mBushing != None:
300 fileDir = 'solution/'
301 # sensBolt = mbs.AddSensor(SensorMarker(markerNumber=mBolt,
302 # fileName=fileDir+'hingePartBoltPos'+str(nModes)+strMode+'.txt',
303 # outputVariableType = exu.OutputVariableType.Position))
304 # sensBushing= mbs.AddSensor(SensorMarker(markerNumber=mBushing,
305 # fileName=fileDir+'hingePartBushingPos'+str(nModes)+strMode+'.txt',
306 # outputVariableType = exu.OutputVariableType.Position))
307 sensBushingVel= mbs.AddSensor(SensorMarker(markerNumber=mBushing,
308 fileName=fileDir+'hingePartBushingVel'+str(nModes)+strMode+'.txt',
309 outputVariableType = exu.OutputVariableType.Velocity))
310 sensBushing= mbs.AddSensor(SensorMarker(markerNumber=mBushing,
311 fileName=fileDir+'hingePartBushing'+str(nModes)+strMode+'.txt',
312 outputVariableType = exu.OutputVariableType.Position))
313
314 mbs.Assemble()
315
316 simulationSettings = exu.SimulationSettings()
317
318 SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
319 SC.visualizationSettings.nodes.drawNodesAsPoint = False
320 SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
321
322 SC.visualizationSettings.nodes.show = False
323 SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
324 SC.visualizationSettings.nodes.basisSize = 0.12
325 SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
326
327 SC.visualizationSettings.view0.scene.showFaceEdges = True
328 SC.visualizationSettings.view0.scene.showFaces = True
329
330 SC.visualizationSettings.sensors.show = True
331 SC.visualizationSettings.sensors.drawSimplified = False
332 SC.visualizationSettings.sensors.defaultSize = 0.01
333
334
335 simulationSettings.solution.file.information = "CMStutorial "+str(nModes)+" "+strMode+"modes"
336
337 h=0.25e-3
338 tEnd = 1
339
340 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
341 simulationSettings.timeIntegration.endTime = tEnd
342 simulationSettings.solution.file.write = False
343 simulationSettings.timeIntegration.verboseMode = 1
344 #simulationSettings.timeIntegration.verboseModeFile = 3
345 simulationSettings.timeIntegration.newton.useModifiedNewton = True
346
347 simulationSettings.solution.sensors.writePeriod = h
348
349 simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
350 #simulationSettings.show.statistics = True
351 simulationSettings.show.computationTime = True
352
353 #create animation:
354 # simulationSettings.solution.recordImagesInterval = 0.005
355 # SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
356 SC.visualizationSettings.view0.window.renderWindowSize=[1920,1080]
357 SC.visualizationSettings.openGL.multiSampling = 4
358
359 if True:
360 if useGraphics:
361 SC.visualizationSettings.general.autoFitScene=False
362
363 SC.renderer.Start()
364 SC.renderer.RestoreSavedState() #load last model view
365
366 SC.renderer.DoIdleTasks() #press space to continue
367
368 if True:
369 # mbs.SolveDynamic(solverType=exu.DynamicSolverType.TrapezoidalIndex2,
370 # simulationSettings=simulationSettings)
371 mbs.SolveDynamic(simulationSettings=simulationSettings)
372 else:
373 mbs.SolveStatic(simulationSettings=simulationSettings)
374
375
376 if useGraphics:
377 SC.renderer.DoIdleTasks()
378 SC.renderer.Stop() #safely close rendering window!
379
380 if mBushing != None:
381 uTip = mbs.GetSensorValues(sensBushing)
382 print("nModes="+strMode, nModes, ", bushing position=", uTip)
383 if False:
384
385 mbs.PlotSensor(sensorNumbers=[sensBushingVel], components=[1])
386
387#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
388if False:
389 import matplotlib.pyplot as plt
390 import matplotlib.ticker as ticker
391 import exudyn as exu
392 from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
393 import exudyn.graphics as graphics
394 CC = PlotLineCode
395 comp = 3 #1=x, 2=y, ...
396 var = 'Vel'
397 # data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'2.txt'), comments='#', delimiter=',')
398 # plt.plot(data[:,0], data[:,comp], CC(7), label='2 eigenmodes')
399 # data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'4.txt'), comments='#', delimiter=',')
400 # plt.plot(data[:,0], data[:,comp], CC(8), label='4 eigenmodes')
401 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'8.txt'), comments='#', delimiter=',')
402 plt.plot(data[:,0], data[:,comp], CC(8), label='8 eigenmodes')
403 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'16.txt'), comments='#', delimiter=',')
404 plt.plot(data[:,0], data[:,comp], CC(9), label='16 eigenmodes')
405 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'32.txt'), comments='#', delimiter=',')
406 plt.plot(data[:,0], data[:,comp], CC(10), label='32 eigenmodes')
407 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'64.txt'), comments='#', delimiter=',')
408 plt.plot(data[:,0], data[:,comp], CC(11), label='64 eigenmodes')
409 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'64.txt'), comments='#', delimiter=',')
410 plt.plot(data[:,0], data[:,comp], CC(12), label='64 eigenmodes')
411 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'128.txt'), comments='#', delimiter=',')
412 plt.plot(data[:,0], data[:,comp], CC(13), label='128 eigenmodes')
413
414 # data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'2HCB.txt'), comments='#', delimiter=',')
415 # plt.plot(data[:,0], data[:,comp], CC(1), label='HCB + 2 eigenmodes')
416 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'4HCB.txt'), comments='#', delimiter=',')
417 plt.plot(data[:,0], data[:,comp], CC(2), label='HCB2 + 4 eigenmodes')
418 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'8HCB.txt'), comments='#', delimiter=',')
419 plt.plot(data[:,0], data[:,comp], CC(3), label='HCB2 + 8 eigenmodes')
420 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'16HCB.txt'), comments='#', delimiter=',')
421 plt.plot(data[:,0], data[:,comp], CC(4), label='HCB2 + 16 eigenmodes')
422 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'32HCB.txt'), comments='#', delimiter=',')
423 plt.plot(data[:,0], data[:,comp], CC(5), label='HCB2 + 32 eigenmodes')
424 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'64HCB.txt'), comments='#', delimiter=',')
425 plt.plot(data[:,0], data[:,comp], CC(6), label='HCB2 + 64 eigenmodes')
426 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'128HCB.txt'), comments='#', delimiter=',')
427 plt.plot(data[:,0], data[:,comp], CC(7), label='HCB2 + 128 eigenmodes')
428
429 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'2HCBsingle.txt'), comments='#', delimiter=',')
430 plt.plot(data[:,0], data[:,comp], CC(14), label='HCB1 + 2 eigenmodes')
431 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'4HCBsingle.txt'), comments='#', delimiter=',')
432 plt.plot(data[:,0], data[:,comp], CC(15), label='HCB1 + 4 eigenmodes')
433 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'8HCBsingle.txt'), comments='#', delimiter=',')
434 plt.plot(data[:,0], data[:,comp], CC(16), label='HCB1 + 8 eigenmodes')
435 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'16HCBsingle.txt'), comments='#', delimiter=',')
436 plt.plot(data[:,0], data[:,comp], CC(17), label='HCB1 + 16 eigenmodes')
437 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'32HCBsingle.txt'), comments='#', delimiter=',')
438 plt.plot(data[:,0], data[:,comp], CC(18), label='HCB1 + 32 eigenmodes')
439 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'64HCBsingle.txt'), comments='#', delimiter=',')
440 plt.plot(data[:,0], data[:,comp], CC(19), label='HCB1 + 64 eigenmodes')
441 data = np.loadtxt(OutputFilePath('solution/hingePartBushing'+var+'128HCBsingle.txt'), comments='#', delimiter=',')
442 plt.plot(data[:,0], data[:,comp], CC(20), label='HCB1 + 128 eigenmodes')
443
444
445 ax=plt.gca() # get current axes
446 ax.grid(True, 'major', 'both')
447 ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
448 ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
449 #
450 plt.xlabel("time (s)")
451 plt.ylabel("y-component of tip velocity of hinge (m)")
452 plt.legend() #show labels as legend
453 plt.tight_layout()
454 plt.show()
455
456#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
457if True:
458 varList = ['','HCB','HCBsingle']
459 for var in varList:
460 for i in range(6):
461 n = 4*2**i
462 filename = 'solution/hingePartBushingVel'+str(n)+var+'.txt'
463 #print(filename)
464 data = np.loadtxt(filename, comments='#', delimiter=',')
465 s = var + " eigenmodes"
466 print("solution with "+str(n)+" "+s+" = ",data[-1,1],", ",data[-1,2],", ",data[-1,3],sep="")
467
468#++++++++++++++++++++++
469#(x,y,z-position results for h=0.25e-3, tEnd = 1:
470# solution with 4 eigenmodes = -0.1218716941, -0.02212539352, -0.3826646827
471# solution with 8 eigenmodes = -0.1246493313, -0.02134551124, -0.3817672439
472# solution with 16 eigenmodes = -0.125718746, -0.02220973667, -0.3817761998
473# solution with 32 eigenmodes = -0.1227923675, -0.02232804332, -0.3826703705
474# solution with 64 eigenmodes = -0.1211624347, -0.02256801385, -0.3830241186
475# solution with 128 eigenmodes = -0.1211098342, -0.02258891649, -0.3830239774
476
477# solution with 4 HCB eigenmodes = -0.137803822, -0.02140481771, -0.377325894
478# solution with 8 HCB eigenmodes = -0.09278682737, -0.02088216306, -0.3910735225
479# solution with 16 HCB eigenmodes = -0.1006048749, -0.0210529449, -0.3890880585
480# solution with 32 HCB eigenmodes = -0.1418260115, -0.02137465745, -0.3755985975
481# solution with 64 HCB eigenmodes = -0.1261576272, -0.02133676138, -0.3811615539
482# solution with 128 HCB eigenmodes = -0.1249497117, -0.02134015915, -0.381582143
483
484# solution with 4 HCBsingle eigenmodes = -0.1236432594, -0.02127703127, -0.3822381713
485# solution with 8 HCBsingle eigenmodes = -0.1553884175, -0.02144366871, -0.3712096711
486# solution with 16 HCBsingle eigenmodes = -0.1096747619, -0.02127260753, -0.3871797944
487# solution with 32 HCBsingle eigenmodes = -0.130126813, -0.02149842833, -0.3807721171
488# solution with 64 HCBsingle eigenmodes = -0.1261109915, -0.02147756767, -0.3821287225
489# solution with 128 HCBsingle eigenmodes = -0.1269092416, -0.02148461514, -0.3818634658
490
491#NOTE: main differences due to different initial conditions (USE offset, bad convergence of HCB modes for gravity, etc.)
492#(x,y,z-velocity results for h=0.25e-3, tEnd = 1:
493# solution with 4 eigenmodes = 2.798215342, 0.0123889876, -0.894408541
494# solution with 8 eigenmodes = 2.753795922, 0.001046355507, -1.033353889
495# solution with 16 eigenmodes = 2.862677224, 0.05041922189, -0.70615996
496# solution with 32 eigenmodes = 2.886092992, 0.04990608422, -0.783893511
497# solution with 64 eigenmodes = 2.82897851, -0.02284196211, -0.9656913985
498# solution with 128 eigenmodes = 2.839233628, 0.001567636751, -0.9556805815
499#
500# solution with 4 HCB eigenmodes = 2.841690471, 0.02171168723, -0.8530592818
501# solution with 8 HCB eigenmodes = 2.96737056, -0.01208003067, -0.6819585453
502# solution with 16 HCB eigenmodes = 2.919615786, -0.01640113107, -0.7205707584
503# solution with 32 HCB eigenmodes = 2.803855522, 0.01284070602, -0.9694702614
504# solution with 64 HCB eigenmodes = 2.86587674, 0.01787123237, -0.8448990047
505# solution with 128 HCB eigenmodes = 2.87133748, 0.03213267314, -0.8176578849
506#
507# solution with 4 HCBsingle eigenmodes = 2.790998662, 0.007480706365, -0.9071953092
508# solution with 8 HCBsingle eigenmodes = 2.71735531, 0.005031127492, -1.102723094
509# solution with 16 HCBsingle eigenmodes = 2.889954015, -0.005524615368, -0.8508318815
510# solution with 32 HCBsingle eigenmodes = 2.856518668, 0.03496577193, -0.8353875884
511# solution with 64 HCBsingle eigenmodes = 2.867595936, 0.03403208487, -0.8067800302
512# solution with 128 HCBsingle eigenmodes = 2.865221368, 0.03422539291, -0.8118038999