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
 15from exudyn.itemInterface import *
 16from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 17import exudyn.graphics as graphics #only import if it does not conflict
 18from exudyn.FEM import *
 19from exudyn.graphicsDataUtilities import *
 20import time
 21
 22SC = exu.SystemContainer()
 23mbs = SC.AddSystem()
 24
 25import numpy as np
 26
 27#import timeit
 28
 29import exudyn.basicUtilities as eb
 30import exudyn.rigidBodyUtilities as rb
 31import exudyn.utilities as eu
 32
 33import numpy as np
 34
 35useGraphics = True
 36fileName = 'testData/netgenHinge' #for load/save of FEM data
 37
 38#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 39#netgen/meshing part:
 40fem = FEMinterface()
 41
 42#geometrical parameters:
 43L = 0.4  #Length of plate (X)
 44w = 0.04 #width of plate  (Y)
 45h = 0.02 #height of plate (Z)
 46d = 0.03    #diameter of bolt
 47D = d*2 #diameter of bushing
 48b = 0.05 #length of bolt
 49nModes = 32 #128
 50meshH = 0.01 #0.01 is default, 0.002 gives 100000 nodes and is fairly converged
 51#meshH = 0.0014 #203443 nodes, takes 1540 seconds for eigenmode computation (free-free) and 753 seconds for postprocessing on i9
 52
 53#steel:
 54rho = 7850
 55Emodulus=2.1e11
 56nu=0.3
 57
 58#test high flexibility
 59Emodulus=2e8
 60# nModes = 32
 61
 62
 63#helper function for cylinder with netgen
 64def CSGcylinder(p0,p1,r):
 65    v = VSub(p1,p0)
 66    v = Normalize(v)
 67    cyl = Cylinder(Pnt(p0[0],p0[1],p0[2]), Pnt(p1[0],p1[1],p1[2]),
 68                   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]))
 69    return cyl
 70
 71mBushing = None
 72meshCreated = False
 73
 74#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 75if True: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
 76
 77    import ngsolve as ngs
 78    import netgen
 79    from netgen.meshing import *
 80
 81    from netgen.geom2d import unit_square
 82    #import netgen.libngpy as libng
 83    from netgen.csg import *
 84
 85    geo = CSGeometry()
 86
 87    #plate
 88    block = OrthoBrick(Pnt(0, 0, -0.5*h),Pnt(L, w, 0.5*h))
 89
 90    #bolt
 91    bolt0 = CSGcylinder(p0=[0,w,0], p1=[0,0,0], r=1.6*h)
 92    bolt = CSGcylinder(p0=[0,0.5*w,0], p1=[0,-b,0], r=0.5*d)
 93
 94    #bushing
 95    bushing = (CSGcylinder(p0=[L,w,0], p1=[L,-b,0], r=0.5*D) -
 96               CSGcylinder(p0=[L,0,0], p1=[L,-b*1.1,0], r=0.5*d))
 97
 98    geo.Add(block+bolt0+bolt+bushing)
 99
100    mesh = ngs.Mesh( geo.GenerateMesh(maxh=meshH))
101    mesh.Curve(1)
102
103    if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
104        # import netgen
105        import netgen.gui
106        ngs.Draw(mesh)
107        netgen.Redraw()
108
109    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
110    #Use fem to import FEM model and create FFRFreducedOrder object
111    fem.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
112    meshCreated = True
113    if (meshH==0.01): fem.SaveToFile(fileName)
114
115#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
116#compute Hurty-Craig-Bampton modes
117if True: #now import mesh as mechanical model to EXUDYN
118    if not meshCreated: fem.LoadFromFile(fileName)
119
120    boltP1=[0,0,0]
121    boltP2=[0,-b,0]
122    nodesOnBolt = fem.GetNodesOnCylinder(boltP1, boltP2, radius=0.5*d)
123    #print("boundary nodes bolt=", nodesOnBolt)
124    nodesOnBoltLen = len(nodesOnBolt)
125    nodesOnBoltWeights = np.array((1./nodesOnBoltLen)*np.ones(nodesOnBoltLen))
126
127    bushingP1=[L,0,0]
128    bushingP2=[L,-b,0]
129    nodesOnBushing = fem.GetNodesOnCylinder(bushingP1, bushingP2, radius=0.5*d)
130    #print("boundary nodes bushing=", nodesOnBushing)
131    nodesOnBushingLen = len(nodesOnBushing)
132    nodesOnBushingWeights = np.array((1./nodesOnBushingLen)*np.ones(nodesOnBushingLen))
133
134    print("nNodes=",fem.NumberOfNodes())
135
136    strMode = ''
137    if False: #pure eigenmodes
138        print("compute eigen modes... ")
139        start_time = time.time()
140        fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
141        print("eigen modes computation needed %.3f seconds" % (time.time() - start_time))
142        print("eigen freq.=", fem.GetEigenFrequenciesHz())
143
144    elif False:
145        strMode = 'HCB'
146        #boundaryList = [nodesOnBolt, nodesOnBolt, nodesOnBushing] #for visualization, use first interface twice
147        boundaryList = [nodesOnBolt, nodesOnBushing]
148
149        print("compute HCB modes... ")
150        start_time = time.time()
151        fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
152                                      nEigenModes=nModes,
153                                      useSparseSolver=True,
154                                      computationMode = HCBstaticModeSelection.RBE2)
155
156        print("eigen freq.=", fem.GetEigenFrequenciesHz())
157        print("HCB modes needed %.3f seconds" % (time.time() - start_time))
158    else:
159        strMode = 'HCBsingle'
160        #boundaryList = [nodesOnBolt, nodesOnBolt, nodesOnBushing] #for visualization, use first interface twice
161        boundaryList = [nodesOnBolt]
162
163        print("compute HCB single modes... ")
164        start_time = time.time()
165        fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
166                                      nEigenModes=nModes,
167                                      useSparseSolver=True,
168                                      computationMode = HCBstaticModeSelection.RBE2)
169
170        print("eigen freq.=", fem.GetEigenFrequenciesHz())
171        print("HCB modes needed %.3f seconds" % (time.time() - start_time))
172
173
174
175    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
176    #compute stress modes for postprocessing (inaccurate for coarse meshes, just for visualization):
177    if False:
178        mat = KirchhoffMaterial(Emodulus, nu, rho)
179        varType = exu.OutputVariableType.StressLocal
180        #varType = exu.OutputVariableType.StrainLocal
181        print("ComputePostProcessingModes ... (may take a while)")
182        start_time = time.time()
183        fem.ComputePostProcessingModes(material=mat,
184                                       outputVariableType=varType)
185        print("   ... needed %.3f seconds" % (time.time() - start_time))
186        SC.visualizationSettings.contour.reduceRange=False
187        SC.visualizationSettings.contour.outputVariable = varType
188        SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
189    else:
190        SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
191        SC.visualizationSettings.contour.outputVariableComponent = -1
192
193    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
194    print("create CMS element ...")
195    cms = ObjectFFRFreducedOrderInterface(fem)
196
197    objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
198                                                  initialVelocity=[0,0,0],
199                                                  initialAngularVelocity=[0,0,0],
200                                                  color=[0.9,0.9,0.9,1.],
201                                                  )
202
203    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
204    #add markers and joints
205    nodeDrawSize = 0.0025 #for joint drawing
206
207
208    #mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
209
210    if False:
211        boltMidPoint = 0.5*(np.array(boltP1)+boltP2)
212
213        oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
214
215        altApproach = True
216        mBolt = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
217                                                      meshNodeNumbers=np.array(nodesOnBolt), #these are the meshNodeNumbers
218                                                      #referencePosition=boltMidPoint,
219                                                      useAlternativeApproach=altApproach,
220                                                      weightingFactors=nodesOnBoltWeights))
221        bushingMidPoint = 0.5*(np.array(bushingP1)+bushingP2)
222
223        #add marker for visualization of boundary nodes
224        mBushing = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
225                                                      meshNodeNumbers=np.array(nodesOnBushing), #these are the meshNodeNumbers
226                                                      #referencePosition=bushingMidPoint,
227                                                      useAlternativeApproach=altApproach,
228                                                      weightingFactors=nodesOnBushingWeights))
229
230        lockedAxes=[1,1,1,1,1*0,1]
231        if True:
232
233            mGroundBolt = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround,
234                                                        localPosition=boltMidPoint,
235                                                        visualization=VMarkerBodyRigid(show=True)))
236            mbs.AddObject(GenericJoint(markerNumbers=[mGroundBolt, mBolt],
237                                        constrainedAxes = lockedAxes,
238                                        visualization=VGenericJoint(show=False, axesRadius=0.1*b, axesLength=0.1*b)))
239
240        else:
241
242            mGroundBushing = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=bushingMidPoint))
243            mbs.AddObject(GenericJoint(markerNumbers=[mGroundBushing, mBushing],
244                                        constrainedAxes = lockedAxes,
245                                        visualization=VGenericJoint(axesRadius=0.1*b, axesLength=0.1*b)))
246
247
248    if False:
249        cms = ObjectFFRFreducedOrderInterface(fem)
250
251        objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
252                                                      initialVelocity=[990,990,990],
253                                                      initialAngularVelocity=[0,0,0],
254                                                      color=[0.9,0.9,0.9,1.],
255                                                      )
256
257    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
258    #animate modes
259    SC.visualizationSettings.markers.show = True
260    SC.visualizationSettings.markers.defaultSize=0.0075
261    SC.visualizationSettings.markers.drawSimplified = False
262
263    SC.visualizationSettings.loads.show = False
264    SC.visualizationSettings.loads.drawSimplified = False
265    SC.visualizationSettings.loads.defaultSize=0.1
266    SC.visualizationSettings.loads.defaultRadius = 0.002
267
268    SC.visualizationSettings.openGL.multiSampling=4
269    SC.visualizationSettings.openGL.lineWidth=2
270
271    if False: #activate to animate modes
272        from exudyn.interactive import AnimateModes
273        mbs.Assemble()
274        SC.visualizationSettings.nodes.show = False
275        SC.visualizationSettings.openGL.showFaceEdges = True
276        SC.visualizationSettings.openGL.multiSampling=4
277        SC.visualizationSettings.openGL.lineWidth=2
278        SC.visualizationSettings.window.renderWindowSize = [1600,1080]
279        SC.visualizationSettings.contour.showColorBar = False
280        SC.visualizationSettings.general.textSize = 16
281
282        #%%+++++++++++++++++++++++++++++++++++++++
283        #animate modes of ObjectFFRFreducedOrder (only needs generic node containing modal coordinates)
284        SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved
285
286        nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
287        AnimateModes(SC, mbs, nodeNumber, period=0.1, showTime=False, renderWindowText='Hurty-Craig-Bampton: 2 x 6 static modes and 8 eigenmodes\n')
288        # import sys
289        # sys.exit()
290
291    #add gravity (not necessary if user functions used)
292    oFFRF = objFFRF['oFFRFreducedOrder']
293    mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber=oFFRF))
294    mbs.AddLoad(LoadMassProportional(markerNumber=mBody, loadVector= [0,0,-9.81*0]))
295
296
297    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
298    if mBushing != None:
299        fileDir = 'solution/'
300        # sensBolt = mbs.AddSensor(SensorMarker(markerNumber=mBolt,
301        #                                       fileName=fileDir+'hingePartBoltPos'+str(nModes)+strMode+'.txt',
302        #                                       outputVariableType = exu.OutputVariableType.Position))
303        # sensBushing= mbs.AddSensor(SensorMarker(markerNumber=mBushing,
304        #                                       fileName=fileDir+'hingePartBushingPos'+str(nModes)+strMode+'.txt',
305        #                                       outputVariableType = exu.OutputVariableType.Position))
306        sensBushingVel= mbs.AddSensor(SensorMarker(markerNumber=mBushing,
307                                              fileName=fileDir+'hingePartBushingVel'+str(nModes)+strMode+'.txt',
308                                              outputVariableType = exu.OutputVariableType.Velocity))
309        sensBushing= mbs.AddSensor(SensorMarker(markerNumber=mBushing,
310                                              fileName=fileDir+'hingePartBushing'+str(nModes)+strMode+'.txt',
311                                              outputVariableType = exu.OutputVariableType.Position))
312
313    mbs.Assemble()
314
315    simulationSettings = exu.SimulationSettings()
316
317    SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
318    SC.visualizationSettings.nodes.drawNodesAsPoint = False
319    SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
320
321    SC.visualizationSettings.nodes.show = False
322    SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
323    SC.visualizationSettings.nodes.basisSize = 0.12
324    SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
325
326    SC.visualizationSettings.openGL.showFaceEdges = True
327    SC.visualizationSettings.openGL.showFaces = True
328
329    SC.visualizationSettings.sensors.show = True
330    SC.visualizationSettings.sensors.drawSimplified = False
331    SC.visualizationSettings.sensors.defaultSize = 0.01
332
333
334    simulationSettings.solutionSettings.solutionInformation = "CMStutorial "+str(nModes)+" "+strMode+"modes"
335
336    h=0.25e-3
337    tEnd = 1
338
339    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
340    simulationSettings.timeIntegration.endTime = tEnd
341    simulationSettings.solutionSettings.writeSolutionToFile = False
342    simulationSettings.timeIntegration.verboseMode = 1
343    #simulationSettings.timeIntegration.verboseModeFile = 3
344    simulationSettings.timeIntegration.newton.useModifiedNewton = True
345
346    simulationSettings.solutionSettings.sensorsWritePeriod = h
347
348    simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
349    #simulationSettings.displayStatistics = True
350    simulationSettings.displayComputationTime = True
351
352    #create animation:
353    # simulationSettings.solutionSettings.recordImagesInterval = 0.005
354    # SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
355    SC.visualizationSettings.window.renderWindowSize=[1920,1080]
356    SC.visualizationSettings.openGL.multiSampling = 4
357
358    if True:
359        if useGraphics:
360            SC.visualizationSettings.general.autoFitScene=False
361
362            exu.StartRenderer()
363            if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
364
365            mbs.WaitForUserToContinue() #press space to continue
366
367        #SC.RedrawAndSaveImage()
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.WaitForRenderEngineStopFlag()
378            exu.StopRenderer() #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 True:
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
393import exudyn.graphics as graphics #only import if it does not conflict
394    CC = PlotLineCode
395    comp = 3 #1=x, 2=y, ...
396    var = 'Vel'
397    # data = np.loadtxt('solution/hingePartBushing'+var+'2.txt', comments='#', delimiter=',')
398    # plt.plot(data[:,0], data[:,comp], CC(7), label='2 eigenmodes')
399    # data = np.loadtxt('solution/hingePartBushing'+var+'4.txt', comments='#', delimiter=',')
400    # plt.plot(data[:,0], data[:,comp], CC(8), label='4 eigenmodes')
401    data = np.loadtxt('solution/hingePartBushing'+var+'8.txt', comments='#', delimiter=',')
402    plt.plot(data[:,0], data[:,comp], CC(8), label='8 eigenmodes')
403    data = np.loadtxt('solution/hingePartBushing'+var+'16.txt', comments='#', delimiter=',')
404    plt.plot(data[:,0], data[:,comp], CC(9), label='16 eigenmodes')
405    data = np.loadtxt('solution/hingePartBushing'+var+'32.txt', comments='#', delimiter=',')
406    plt.plot(data[:,0], data[:,comp], CC(10), label='32 eigenmodes')
407    data = np.loadtxt('solution/hingePartBushing'+var+'64.txt', comments='#', delimiter=',')
408    plt.plot(data[:,0], data[:,comp], CC(11), label='64 eigenmodes')
409    data = np.loadtxt('solution/hingePartBushing'+var+'64.txt', comments='#', delimiter=',')
410    plt.plot(data[:,0], data[:,comp], CC(12), label='64 eigenmodes')
411    data = np.loadtxt('solution/hingePartBushing'+var+'128.txt', comments='#', delimiter=',')
412    plt.plot(data[:,0], data[:,comp], CC(13), label='128 eigenmodes')
413
414    # data = np.loadtxt('solution/hingePartBushing'+var+'2HCB.txt', comments='#', delimiter=',')
415    # plt.plot(data[:,0], data[:,comp], CC(1), label='HCB + 2 eigenmodes')
416    data = np.loadtxt('solution/hingePartBushing'+var+'4HCB.txt', comments='#', delimiter=',')
417    plt.plot(data[:,0], data[:,comp], CC(2), label='HCB2 + 4 eigenmodes')
418    data = np.loadtxt('solution/hingePartBushing'+var+'8HCB.txt', comments='#', delimiter=',')
419    plt.plot(data[:,0], data[:,comp], CC(3), label='HCB2 + 8 eigenmodes')
420    data = np.loadtxt('solution/hingePartBushing'+var+'16HCB.txt', comments='#', delimiter=',')
421    plt.plot(data[:,0], data[:,comp], CC(4), label='HCB2 + 16 eigenmodes')
422    data = np.loadtxt('solution/hingePartBushing'+var+'32HCB.txt', comments='#', delimiter=',')
423    plt.plot(data[:,0], data[:,comp], CC(5), label='HCB2 + 32 eigenmodes')
424    data = np.loadtxt('solution/hingePartBushing'+var+'64HCB.txt', comments='#', delimiter=',')
425    plt.plot(data[:,0], data[:,comp], CC(6), label='HCB2 + 64 eigenmodes')
426    data = np.loadtxt('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('solution/hingePartBushing'+var+'2HCBsingle.txt', comments='#', delimiter=',')
430    plt.plot(data[:,0], data[:,comp], CC(14), label='HCB1 + 2 eigenmodes')
431    data = np.loadtxt('solution/hingePartBushing'+var+'4HCBsingle.txt', comments='#', delimiter=',')
432    plt.plot(data[:,0], data[:,comp], CC(15), label='HCB1 + 4 eigenmodes')
433    data = np.loadtxt('solution/hingePartBushing'+var+'8HCBsingle.txt', comments='#', delimiter=',')
434    plt.plot(data[:,0], data[:,comp], CC(16), label='HCB1 + 8 eigenmodes')
435    data = np.loadtxt('solution/hingePartBushing'+var+'16HCBsingle.txt', comments='#', delimiter=',')
436    plt.plot(data[:,0], data[:,comp], CC(17), label='HCB1 + 16 eigenmodes')
437    data = np.loadtxt('solution/hingePartBushing'+var+'32HCBsingle.txt', comments='#', delimiter=',')
438    plt.plot(data[:,0], data[:,comp], CC(18), label='HCB1 + 32 eigenmodes')
439    data = np.loadtxt('solution/hingePartBushing'+var+'64HCBsingle.txt', comments='#', delimiter=',')
440    plt.plot(data[:,0], data[:,comp], CC(19), label='HCB1 + 64 eigenmodes')
441    data = np.loadtxt('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