bicycleIftommBenchmark.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  bicycle Iftomm benchmark example
  5#           https://www.iftomm-multibody.org/benchmark/problem/Uncontrolled_bicycle/
  6#
  7# Author:   Johannes Gerstmayr
  8# Date:     2021-06-22
  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
 14
 15import exudyn as exu
 16from exudyn.itemInterface import *
 17from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 18import exudyn.graphics as graphics #only import if it does not conflict
 19from exudyn.graphicsDataUtilities import *
 20
 21from math import sin, cos, pi
 22import numpy as np
 23
 24SC = exu.SystemContainer()
 25mbs = SC.AddSystem()
 26
 27
 28#%%++++++++++++++++++++++++++++++++++++++++++++++++
 29#coordinate system according to IFToMM:
 30#note: here, wheels are rotated as local wheel axis=x, z points upwards in EXUDYN model
 31#                            ox P2
 32#                        oooo  o
 33#                    oooo       o
 34#      +++       oooo            o +++
 35#     +   +   ooo                 o   +
 36#    +   oooo                    + o   +
 37#    +  xP1+                     +  xP3+
 38#    +     +                     +     +
 39#     +   +                       +   +
 40#      +++                         +++
 41#       x---------------------------x----------------------> x
 42#       |         <- w ->
 43#       v  z
 44
 45
 46#parameters
 47sZ = -1         #switch z coordinate compared to IFToMM description
 48w = 1.02        #wheel base (distance of wheel centers)
 49c = 0.08        #trail
 50lam = pi/10     #steer axis tilt (rad)
 51g = [0,0,9.81*sZ]     #gravity in m/s^2
 52
 53#rear wheel R:
 54rR = 0.3            #rear wheel radius
 55mR = 2              #rear wheel mass
 56IRxx = 0.0603       #rear wheel inertia xx = zz ; but in EXUDYN, x=rotation axis!
 57IRyy = 0.12         #rear wheel inertia yy (around wheel axis)
 58inertiaR = RigidBodyInertia(mass=mR, inertiaTensor=np.array([[IRyy,0,0],[0,IRxx,0],[0,0,IRxx]]))
 59
 60#front wheel F:
 61rF = 0.35           #rear wheel radius
 62mF = 3              #rear wheel mass
 63IFxx = 0.1405       #rear wheel inertia xx = zz ; but in EXUDYN, x=rotation axis!
 64IFyy = 0.28         #rear wheel inertia yy (around wheel axis)
 65inertiaF = RigidBodyInertia(mass=mF, inertiaTensor=np.array([[IFyy,0,0],[0,IFxx,0],[0,0,IFxx]]))
 66
 67#rear body B:
 68xB = 0.3            #COM x
 69zB = -0.9*sZ        #COM z
 70bCOM = np.array([xB, 0, zB])
 71mB = 85             #rear body (and driver) mass
 72zz=-1
 73inertiaB = RigidBodyInertia(mass=mB,
 74                            inertiaTensor=np.array([[9.2,0,2.4*zz],[0,11,0],[2.4*zz,0,2.8]]),
 75                            # inertiaTensor=np.diag([1,1,1]),
 76                            com=np.zeros(3)) #reference position = COM for this body
 77                            # com=bCOM)
 78
 79#front Handlebar H:
 80xH = 0.9            #COM x
 81zH = -0.7*sZ        #COM z
 82hCOM = np.array([xH, 0, zH])
 83mH = 4              #handle bar mass
 84inertiaH = RigidBodyInertia(mass=mH,
 85                            inertiaTensor=np.array([[0.05892, 0, -0.00756*zz],[0,0.06,0],[-0.00756*zz, 0, 0.00708]]),
 86                            # inertiaTensor=np.diag([1,1,1]),
 87                            com=np.zeros(3)) #reference position = COM for this body
 88                            # com=hCOM)
 89
 90#geometrical parameters for joints
 91P1 = np.array([0,0,-0.3*sZ])
 92P2 = np.array([0.82188470506, 0, -0.85595086466*sZ])
 93P3 = np.array([w, 0, -0.35*sZ])
 94
 95
 96#stable velocity limits according to linear theory:
 97vMin = 4.29238253634111
 98vMax = 6.02426201538837
 99
100maneuver = 'M1'
101if maneuver == 'M1':
102    vX0 = 4.                #initial forward velocity in x-direction
103    omegaX0 = 0.05          #initial roll velocity around x axis
104elif maneuver == 'M2':
105    vX0 = 4.6               #initial forward velocity in x-direction
106    omegaX0 = 0.5           #initial roll velocity around x axis
107elif maneuver == 'M3':
108    vX0 = 8                 #initial forward velocity in x-direction
109    omegaX0 = 0.05          #initial roll velocity around x axis
110
111omegaRy0 = -vX0/rR*sZ   #initial angular velocity of rear wheel
112omegaFy0 = -vX0/rF*sZ   #initial angular velocity of front wheel
113
114#%%++++++++++++++++++++++++++++++++++++++++++++++++
115#visualization:
116dY = 0.02
117#graphicsFrame = graphics.Brick(centerPoint=[0,0,0],size=[dFoot*1.1,0.7*rFoot,0.7*rFoot], color=graphics.color.lightred)
118graphicsR = graphics.Cylinder(pAxis=[-1*dY,0,0], vAxis=[dY*2,0,0], radius=rR, color=graphics.color.steelblue, nTiles=4)
119graphicsF = graphics.Cylinder(pAxis=[-1*dY,0,0], vAxis=[dY*2,0,0], radius=rF, color=graphics.color.steelblue, nTiles=4)
120graphicsB = graphics.Cylinder(pAxis=P1-bCOM, vAxis=P2-P1, radius=dY*1.5, color=graphics.color.lightred)
121graphicsB2 = graphics.Sphere(point=[0,0,0], radius=3*dY, color=graphics.color.lightgrey)
122graphicsH = graphics.Cylinder(pAxis=P3-hCOM, vAxis=P2-P3, radius=dY*1.3, color=graphics.color.lightgreen)
123
124#option to track motion of bicycle
125#old version; does not work any more to set state during graphics user function => use trackMarker
126if False:
127    #add user function to track bicycle frame
128    def UFgraphics(mbs, objectNum):
129        n = mbs.variables['nTrackNode']
130        p = mbs.GetNodeOutput(n,exu.OutputVariableType.Position,
131                              configuration=exu.ConfigurationType.Visualization)
132        rs=SC.renderer.GetState()
133        A = np.array(rs['modelRotation'])
134        p = A.T @ p
135        rs['centerPoint']=[p[0],p[1],p[2]]
136        SC.renderer.SetState(rs)
137        return []
138
139    #add object with graphics user function
140    oGround2 = mbs.AddObject(ObjectGround(visualization=
141                                          VObjectGround(graphicsDataUserFunction=UFgraphics)))
142#add rigid bodies
143#rear wheel
144resultR = mbs.CreateRigidBody(
145    referencePosition = P1,
146    referenceRotationMatrix = RotationMatrixZ(np.pi*0.5),
147    initialVelocity = [vX0, omegaX0*P1[2]*sZ, 0],
148    initialAngularVelocity = [omegaX0, omegaRy0, 0],
149    inertia = inertiaR,
150    gravity = g,
151    graphicsDataList = [graphicsR],
152    returnDict = True)
153
154nR, bR = resultR['nodeNumber'], resultR['bodyNumber']
155
156mbs.variables['nTrackNode'] = nR #node to be tracked
157
158#front wheel
159resultF = mbs.CreateRigidBody(
160    referencePosition = P3,
161    referenceRotationMatrix = RotationMatrixZ(pi*0.5),
162    initialVelocity = [vX0, omegaX0*P3[2]*sZ, 0],
163    initialAngularVelocity = [omegaX0, omegaFy0, 0],
164    inertia = inertiaF,
165    gravity = g,
166    nodeType = exu.NodeType.RotationEulerParameters,
167    graphicsDataList = [graphicsF],
168    returnDict = True
169)
170nF, bF = resultF['nodeNumber'], resultF['bodyNumber']
171
172#read body
173resultB = mbs.CreateRigidBody(
174    referencePosition = bCOM,
175    initialVelocity = [vX0, omegaX0*bCOM[2]*sZ, 0],
176    initialAngularVelocity = [omegaX0, 0, 0],
177    inertia = inertiaB,
178    gravity = g,
179    nodeType = exu.NodeType.RotationEulerParameters,
180    graphicsDataList = [graphicsB, graphicsB2],
181    returnDict = True
182)
183nB, bB = resultB['nodeNumber'], resultB['bodyNumber']
184
185#handle
186resultH = mbs.CreateRigidBody(
187    referencePosition = hCOM,
188    initialVelocity = [vX0, omegaX0*hCOM[2]*sZ, 0],
189    initialAngularVelocity = [omegaX0, 0, 0],
190    inertia = inertiaH,
191    gravity = g,
192    nodeType = exu.NodeType.RotationEulerParameters,
193    graphicsDataList = [graphicsH],
194    returnDict = True
195)
196nH, bH = resultH['nodeNumber'], resultH['bodyNumber']
197
198
199#%%++++++++++++++++++++++++++++++++++++++++++++++++
200#ground body and marker
201gGround = graphics.CheckerBoard(point=[0,0,0], size=50, nTiles=64)
202oGround = mbs.CreateGround(graphicsDataList=[gGround])
203markerGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,0]))
204
205markerR = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bR, localPosition=[0,0,0]))
206markerF = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bF, localPosition=[0,0,0]))
207markerB1 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bB, localPosition=P1-bCOM))
208
209sMarkerR = mbs.AddSensor(SensorMarker(markerNumber=markerR, outputVariableType=exu.OutputVariableType.Position))
210sMarkerB1= mbs.AddSensor(SensorMarker(markerNumber=markerB1,outputVariableType=exu.OutputVariableType.Position))
211
212#%%++++++++++++++++++++++++++++++++++++++++++++++++
213#add joints:
214useJoints = True
215if useJoints:
216    oJointRW = mbs.CreateRevoluteJoint(bodyNumbers=[bR, bB], position=P1, axis=[0,1,0],
217                            axisRadius=0.5*dY, axisLength=5*dY)
218    oJointFW = mbs.CreateRevoluteJoint(bodyNumbers=[bF, bH], position=P3, axis=[0,1,0],
219                            axisRadius=0.5*dY, axisLength=5*dY)
220    oJointSteer = mbs.CreateRevoluteJoint(bodyNumbers=[bB, bH],
221                                          position=P2-bCOM, useGlobalFrame=False,
222                                          axis=RotationMatrixY(-lam) @ [0,0,1],
223                                          axisRadius=0.5*dY, axisLength=5*dY)
224#%%++++++++++++++++++++++++++++++++++++++++++++++++
225#add 'rolling disc' for wheels:
226cStiffness = 5e4*10 #spring stiffness: 50N==>F/k = u = 0.001m (penetration)
227cDamping = cStiffness*0.05*0.1 #think on a one-mass spring damper
228nGenericR = mbs.AddNode(NodeGenericData(initialCoordinates=[0,0,0], numberOfDataCoordinates=3))
229if False:
230    oRollingR=mbs.AddObject(ObjectConnectorRollingDiscPenalty(markerNumbers=[markerGround, markerR],
231                                                              nodeNumber = nGenericR,
232                                                              discRadius=rR,
233                                                              planeNormal=[0,0,1],
234                                                              dryFriction=[0.8,0.8],
235                                                              dryFrictionProportionalZone=1e-2,
236                                                              rollingFrictionViscous=0.,
237                                                              contactStiffness=cStiffness,
238                                                              contactDamping=cDamping,
239                                                              #activeConnector = False, #set to false to deactivated
240                                                              visualization=VObjectConnectorRollingDiscPenalty(show=True,
241                                                                                                               discWidth=dY, color=graphics.color.blue)))
242
243    nGenericF = mbs.AddNode(NodeGenericData(initialCoordinates=[0,0,0], numberOfDataCoordinates=3))
244    oRollingF=mbs.AddObject(ObjectConnectorRollingDiscPenalty(markerNumbers=[markerGround, markerF],
245                                                              nodeNumber = nGenericF,
246                                                              discRadius=rF,
247                                                              planeNormal=[0,0,1],
248                                                              dryFriction=[0.8,0.8],
249                                                              dryFrictionProportionalZone=1e-2,
250                                                              rollingFrictionViscous=0.,
251                                                              contactStiffness=cStiffness,
252                                                              contactDamping=cDamping,
253                                                              #activeConnector = False, #set to false to deactivated
254                                                              visualization=VObjectConnectorRollingDiscPenalty(show=True, discWidth=dY, color=graphics.color.blue)))
255else:
256    if True:
257        oRollingR=mbs.AddObject(ObjectJointRollingDisc(markerNumbers=[markerGround, markerR],
258                                                        discRadius=rR,
259                                                        visualization=VObjectJointRollingDisc(show=True, discWidth=dY, color=graphics.color.blue)))
260
261        oRollingF=mbs.AddObject(ObjectJointRollingDisc(markerNumbers=[markerGround, markerF],
262                                                        discRadius=rF,
263                                                        visualization=VObjectJointRollingDisc(show=True, discWidth=dY, color=graphics.color.blue)))
264
265
266
267#%%++++++++++++++++++++++++++++++++++++++++++++++++
268#add sensors
269addSensors = True
270if addSensors:
271    sForwardVel = mbs.AddSensor(SensorBody(bodyNumber=bB, fileName='solution/bicycleBvelLocal.txt',
272                                          localPosition=P1-bCOM,
273                                          outputVariableType=exu.OutputVariableType.VelocityLocal))
274
275    sBAngVelLocal = mbs.AddSensor(SensorBody(bodyNumber=bB, fileName='solution/bicycleBangVelLocal.txt',
276                                                outputVariableType=exu.OutputVariableType.AngularVelocityLocal))
277    sBrot = mbs.AddSensor(SensorBody(bodyNumber=bB, fileName='solution/bicycleBrot.txt',
278                                                outputVariableType=exu.OutputVariableType.Rotation))
279
280
281    bodies = [bB, bH, bR, bF]
282    massBodies = [mB, mH, mR, mF]
283    inertiaBodies = [inertiaB.inertiaTensor,
284                     inertiaH.inertiaTensor,
285                     inertiaR.inertiaTensor,
286                     inertiaF.inertiaTensor]
287
288    nBodies = len(bodies)
289    sList = []
290    for b in bodies:
291        sPosCOM = mbs.AddSensor(SensorBody(bodyNumber=b, writeToFile=False,
292                                                    outputVariableType=exu.OutputVariableType.Position))
293        sVelCOM = mbs.AddSensor(SensorBody(bodyNumber=b, writeToFile=False,
294                                                    outputVariableType=exu.OutputVariableType.Velocity))
295        sAngVelLocal = mbs.AddSensor(SensorBody(bodyNumber=b, writeToFile=False,
296                                                    outputVariableType=exu.OutputVariableType.AngularVelocityLocal))
297
298        sList += [sPosCOM,sVelCOM,sAngVelLocal]
299
300    if useJoints:
301        sSteerAngle = mbs.AddSensor(SensorObject(objectNumber=oJointSteer, fileName='solution/bicycleSteerAngle.txt',
302                                                    outputVariableType=exu.OutputVariableType.Rotation))
303        sSteerVel = mbs.AddSensor(SensorObject(objectNumber=oJointSteer, fileName='solution/bicycleSteerVelocity.txt',
304                                                    outputVariableType=exu.OutputVariableType.AngularVelocityLocal))
305
306
307    #create user joint for kinetic and potential energy
308    def UFsensorEnergy(mbs, t, sensorNumbers, factors, configuration):
309        E = 0
310        P = 0
311        for i in range(nBodies):
312            pos = mbs.GetSensorValues(sensorNumbers[i*3+0])
313            vel = mbs.GetSensorValues(sensorNumbers[i*3+1]) #vel
314            omega = mbs.GetSensorValues(sensorNumbers[i*3+2]) #ang vel local
315            E += 0.5 * NormL2(vel)**2 * massBodies[i]
316            E += 0.5 * np.array(omega) @ inertiaBodies[i] @ omega
317
318            P -= np.dot(g,pos)*massBodies[i]
319        return [P, E, E+P] #return potential, kinetic and total mechanical energy
320
321    sEnergy = mbs.AddSensor(SensorUserFunction(sensorNumbers=sList, #factors=[180/pi],
322                                             fileName='solution/sensorKineticPotentialEnergy.txt',
323                                             sensorUserFunction=UFsensorEnergy))
324
325    def UFsensorResults(mbs, t, sensorNumbers, factors, configuration):
326        energy =        mbs.GetSensorValues(sensorNumbers[0])
327        forwardVel =    mbs.GetSensorValues(sensorNumbers[1])
328        angVelLocalB =  mbs.GetSensorValues(sensorNumbers[2])
329        rotB =          mbs.GetSensorValues(sensorNumbers[3])
330        steerAngle =    mbs.GetSensorValues(sensorNumbers[4])
331        steerVel =      mbs.GetSensorValues(sensorNumbers[5])
332        return [rotB[0], angVelLocalB[0], forwardVel[0], energy[0], energy[1], energy[2], -steerAngle[2], -steerVel[2]] #return kinetic, potential and total mechanical energy
333
334    # 1=roll angle, 2=roll angular velocity, 3=forward speed, 4=potential energy, 5=kinetic energy, 6=mechanical energy, 7=steer angle, and 8=steer velocity
335    sResults = mbs.AddSensor(SensorUserFunction(sensorNumbers=[sEnergy,sForwardVel,sBAngVelLocal,sBrot, sSteerAngle, sSteerVel],
336                                             fileName='solution/sensorResults'+maneuver+'.txt',
337                                             sensorUserFunction=UFsensorResults))
338
339#%%++++++++++++++++++++++++++++++++++++++++++++++++
340#simulate:
341mbs.Assemble()
342
343pR = mbs.GetSensorValues(sMarkerR)
344pB1 = mbs.GetSensorValues(sMarkerB1)
345print("pR=",pR)
346print("pB1=",pB1)
347simulationSettings = exu.SimulationSettings() #takes currently set values or default values
348
349tEnd = 5
350h=0.001  #use small step size to detext contact switching
351
352simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
353simulationSettings.timeIntegration.endTime = tEnd
354simulationSettings.solutionSettings.writeSolutionToFile= False #set False for CPU performance measurement
355simulationSettings.solutionSettings.sensorsWritePeriod = 0.01
356simulationSettings.solutionSettings.outputPrecision = 16
357
358simulationSettings.timeIntegration.verboseMode = 1
359#simulationSettings.linearSolverSettings.ignoreSingularJacobian = True
360
361# simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
362# simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
363simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.7
364simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations=True
365simulationSettings.timeIntegration.newton.useModifiedNewton = True
366
367#SC.visualizationSettings.general.useMultiThreadedRendering = False
368SC.visualizationSettings.nodes.show = True
369SC.visualizationSettings.nodes.drawNodesAsPoint  = False
370SC.visualizationSettings.nodes.showBasis = True
371SC.visualizationSettings.nodes.basisSize = 0.015
372
373SC.visualizationSettings.view0.camera.trackMarker = markerR #rear wheel
374
375if False: #record animation frames:
376    SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
377    SC.visualizationSettings.view0.window.renderWindowSize=[1600,1024]
378    SC.visualizationSettings.openGL.multiSampling = 4
379    simulationSettings.solutionSettings.recordImagesInterval = 0.02
380
381SC.visualizationSettings.general.autoFitScene = False #use loaded render state
382useGraphics = True
383if useGraphics:
384    SC.renderer.Start()
385    if 'renderState' in exu.sys:
386        SC.renderer.SetState(exu.sys[ 'renderState' ])
387    SC.renderer.DoIdleTasks()
388
389mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.TrapezoidalIndex2)
390#mbs.SolveDynamic(simulationSettings, showHints=True)
391
392
393#%%+++++++++++++++++++++++++++++
394if useGraphics:
395    SC.renderer.DoIdleTasks()
396    SC.renderer.Stop() #safely close rendering window!
397
398#%%++++++++++++++++++++++++++++++++++++++++++++++q+++++++
399if addSensors:
400    #plot results
401
402
403
404    import matplotlib.pyplot as plt
405    import matplotlib.ticker as ticker
406    plt.close('all')
407
408
409    # mbs.PlotSensor(sensorNumbers=[sBpos,sBpos,sBpos], components=[0,1,2])
410    #plt.figure('lateral position')
411    #mbs.PlotSensor(sensorNumbers=[sBpos], components=[1])
412
413    plt.figure('forward velocity')
414    mbs.PlotSensor(sensorNumbers=[sForwardVel], components=[0])
415    # mbs.PlotSensor(sensorNumbers=[sBvelLocal,sBvelLocal,sBvelLocal], components=[0,1,2])
416
417    # plt.figure('local ang velocities')
418    # mbs.PlotSensor(sensorNumbers=[sBAngVelLocal,sBAngVelLocal,sBAngVelLocal], components=[0,1,2])
419    # if False:
420    #     import matplotlib.pyplot as plt
421    #     import matplotlib.ticker as ticker
422
423    # 1=roll angle, 2=roll angular velocity, 3=forward speed, 4=potential energy, 5=kinetic energy, 6=mechanical energy, 7=steer angle, and 8=steer velocity
424    data = np.loadtxt('solution/uncontrolledBicycleGonzalez.txt')#, comments='#', delimiter='')
425    plt.plot(data[:,0], data[:,9], 'b:',label='')
426
427    data2 = np.loadtxt('solution/uncontrolledBicycleSanjurjo.txt')#, comments='#', delimiter='')
428    plt.plot(data2[:,0], data2[:,3+8], 'g:',label='')
429
430    plt.figure('steer vel')
431    mbs.PlotSensor(sensorNumbers=[sSteerVel], components=[2])
432    plt.plot(data2[:,0], -data2[:,8+8], 'g:',label='')
433
434    plt.figure('steer ang')
435    mbs.PlotSensor(sensorNumbers=[sSteerAngle], components=[2])
436    plt.plot(data2[:,0], -data2[:,7+8], 'g:',label='')
437
438    plt.figure('roll')
439    mbs.PlotSensor(sensorNumbers=[sBrot], components=[0])
440    plt.plot(data2[:,0], data2[:,1+8], 'g:',label='')
441
442    plt.figure('roll ang vel')
443    mbs.PlotSensor(sensorNumbers=[sBAngVelLocal], components=[0])
444    plt.plot(data2[:,0], data2[:,2+8], 'g:',label='')
445
446
447    plt.figure('potential energy')
448    mbs.PlotSensor(sensorNumbers=[sEnergy], components=[0])
449    plt.plot(data2[:,0], data2[:,4+8], 'g:',label='')
450
451    plt.figure('kinetic energy')
452    mbs.PlotSensor(sensorNumbers=[sEnergy], components=[1])
453    plt.plot(data2[:,0], data2[:,5+8], 'g:',label='')
454
455    plt.figure('total energy')
456    mbs.PlotSensor(sensorNumbers=[sEnergy], components=[2])
457
458    dataE = np.loadtxt('solution/sensorKineticPotentialEnergy.txt', comments='#', delimiter=',')
459    performance = 100*(max(dataE[:,3]) - min(dataE[:,3])) / dataE[0,3]
460    print("performance = ", performance, "(must by < 1e-3)")
461
462    #CPU performance with 20 seconds simulation time, maneuver 2
463    #performance = 0.000915 < 0.001: max h=0.014; CPU time = 0.596 seconds on Intel Core i9
464    #reference solution computed with:
465    #  performance = 6.423e-06: max h=0.001; CPU time = 5.5935 seconds on Intel Core i9
466
467
468#%%+++++++++++++++++
469#merge result files for IFToMM
470if True:
471    dataM1 = np.loadtxt('solution/sensorResultsM1.txt', comments='#', delimiter=',')
472    dataM2 = np.loadtxt('solution/sensorResultsM2.txt', comments='#', delimiter=',')
473    dataM3 = np.loadtxt('solution/sensorResultsM3.txt', comments='#', delimiter=',')
474
475    data = np.hstack((dataM1,dataM2[:,1:],dataM3[:,1:]))
476    np.savetxt('solution/bicycleResultsIFToMM.txt', data, fmt='%1.15e')
477
478# benchmark results:
479# 6.423e-06
480# 5.5935
481# Intel(R) Core(TM) i9-7940X CPU @ 3.10GHz
482# Simulated using C++/Python library EXUDYN, see https://github.com/jgerstmayr/EXUDYN.
483# Solved using implicit trapezoidal rule (energy conserving) with Index 2 constraint reduction, using redundant coordinate formulation. Rigid bodies are modeled with Euler parameters.
484# With a larger step size of 0.014 we still obtain a performance <0.001, but have CPU time of 0.596 seconds.