solutionViewerTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for AddRevoluteJoint utility function
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2021-07-01
  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
 13import exudyn as exu
 14from exudyn.itemInterface import *
 15from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 16import exudyn.graphics as graphics #only import if it does not conflict
 17
 18from math import sin, cos, pi
 19import numpy as np
 20
 21SC = exu.SystemContainer()
 22mbs = SC.AddSystem()
 23
 24
 25#background
 26color = [0.1,0.1,0.8,1]
 27L = 0.4 #length of bodies
 28d = 0.1 #diameter of bodies
 29
 30oGround=mbs.AddObject(ObjectGround(referencePosition= [-0.5*L,0,0]))
 31mPosLast = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oGround, localPosition=[0,0,0]))
 32A0 = np.eye(3)
 33Alast = A0 #previous marker
 34bodyLast = oGround
 35
 36Alist=[]
 37axisList=[]
 38A=RotationMatrixX(0)
 39
 40nBodies = 100
 41for i in range(nBodies):
 42    delta = 0.01*pi
 43    #A = RotationMatrixZ(delta)
 44    v0= A@[0,0,1]
 45    Alist+=[A]
 46    axisList+=[v0]
 47    A = A @ RotationMatrixX(delta)@RotationMatrixZ(2*delta)
 48
 49
 50p0 = [0.,0.,0] #reference position
 51vLoc = np.array([L,0,0]) #last to next joint
 52g = [0,0,9.81]
 53#g = [0,9.81,0]
 54
 55#create a chain of bodies:
 56for i in range(nBodies):
 57    #print("Build Object", i)
 58    inertia = InertiaCuboid(density=1000, sideLengths=[L,d,d])
 59    p0 += Alist[i] @ (0.5*vLoc)
 60    #p0 += (0.5*vLoc)
 61
 62    ep0 = eulerParameters0 #no rotation
 63    graphicsBody = graphics.Brick([0,0,0], [0.96*L,d,d], graphics.color.steelblue)
 64    oRB = mbs.CreateRigidBody(inertia=inertia,
 65                              referencePosition=p0,
 66                              referenceRotationMatrix=Alist[i],
 67                              gravity=g,
 68                              graphicsDataList=[graphicsBody])
 69    nRB= mbs.GetObject(oRB)['nodeNumber']
 70
 71    body0 = bodyLast
 72    body1 = oRB
 73    # point = mbs.GetObjectOutputBody(oRB,exu.OutputVariableType.Position,
 74    #                                 localPosition=[-0.5*L,0,0],
 75    #                                 configuration=exu.ConfigurationType.Reference)
 76    #axis = [0,0,1]
 77    axis = axisList[i]
 78    mbs.CreateRevoluteJoint(bodyNumbers=[body0, body1], position=[0.5*L,0,0],
 79                            axis=Alast.T@axis, useGlobalFrame=False,
 80                            axisRadius=0.6*d, axisLength=1.2*d)
 81    # mbs.CreateRevoluteJoint(bodyNumbers=[body0, body1], position=point,
 82    #                         axis=axis, useGlobalFrame=True,
 83    #                         axisRadius=0.6*d, axisLength=1.2*d)
 84
 85    bodyLast = oRB
 86
 87    p0 += Alist[i] @ (0.5*vLoc)
 88    #p0 += (0.5*vLoc)
 89    Alast = Alist[i]
 90
 91#mbs.AddLoad(LoadForceVector(markerNumber=mPosLast, loadVector=[0,0,20]))
 92
 93mbs.Assemble()
 94
 95simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 96
 97tEnd = 1
 98h=0.0005  #use small step size to detext contact switching
 99
100simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
101simulationSettings.timeIntegration.endTime = tEnd
102simulationSettings.solutionSettings.solutionWritePeriod = 0.005
103simulationSettings.solutionSettings.sensorsWritePeriod = 0.01
104#simulationSettings.timeIntegration.simulateInRealtime = True
105simulationSettings.timeIntegration.realtimeFactor = 0.5
106simulationSettings.timeIntegration.verboseMode = 1
107
108simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
109simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations=True
110simulationSettings.timeIntegration.newton.useModifiedNewton = True
111#simulationSettings.timeIntegration.newton.modifiedNewtonJacUpdatePerStep = True
112simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
113# simulationSettings.parallel.numberOfThreads=4
114
115SC.visualizationSettings.nodes.show = True
116SC.visualizationSettings.nodes.drawNodesAsPoint  = False
117SC.visualizationSettings.nodes.showBasis = True
118SC.visualizationSettings.nodes.basisSize = 0.015
119SC.visualizationSettings.connectors.showJointAxes = True
120
121#for snapshot:
122SC.visualizationSettings.openGL.multiSampling=4
123SC.visualizationSettings.openGL.lineWidth=2
124SC.visualizationSettings.window.renderWindowSize = [800,600]
125SC.visualizationSettings.general.drawCoordinateSystem=False
126SC.visualizationSettings.general.drawWorldBasis=True
127# SC.visualizationSettings.general.useMultiThreadedRendering = False
128SC.visualizationSettings.general.autoFitScene = False #use loaded render state
129useGraphics = True
130if useGraphics:
131    simulationSettings.displayComputationTime = True
132    simulationSettings.displayStatistics = True
133    exu.StartRenderer()
134    if 'renderState' in exu.sys:
135        SC.SetRenderState(exu.sys[ 'renderState' ])
136    #mbs.WaitForUserToContinue()
137else:
138    simulationSettings.solutionSettings.writeSolutionToFile = False
139
140#mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.TrapezoidalIndex2)
141mbs.SolveDynamic(simulationSettings, showHints=True)
142
143if True: #use this to reload the solution and use SolutionViewer
144    #sol = LoadSolutionFile('coordinatesSolution.txt')
145
146    mbs.SolutionViewer() #can also be entered in IPython ...
147
148
149u0 = mbs.GetNodeOutput(nRB, exu.OutputVariableType.Displacement)
150rot0 = mbs.GetNodeOutput(nRB, exu.OutputVariableType.Rotation)
151exu.Print('u0=',u0,', rot0=', rot0)
152
153result = (abs(u0)+abs(rot0)).sum()
154exu.Print('solution of addRevoluteJoint=',result)
155
156
157
158#%%+++++++++++++++++++++++++++++
159if useGraphics:
160    SC.WaitForRenderEngineStopFlag()
161    exu.StopRenderer() #safely close rendering window!