solutionViewerTest.py

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

  1# -*- coding: utf-8 -*-
  2#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  3# This is an EXUDYN example
  4#
  5# Details:  Test for AddRevoluteJoint utility function
  6#
  7# Author:   Johannes Gerstmayr
  8# Date:     2021-07-01
  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.itemInterface import *
 16from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 17import exudyn.graphics as graphics #only import if it does not conflict
 18
 19from math import sin, cos, pi
 20import numpy as np
 21
 22SC = exu.SystemContainer()
 23mbs = SC.AddSystem()
 24
 25
 26#background
 27color = [0.1,0.1,0.8,1]
 28L = 0.4 #length of bodies
 29d = 0.1 #diameter of bodies
 30
 31oGround=mbs.AddObject(ObjectGround(referencePosition= [-0.5*L,0,0]))
 32mPosLast = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oGround, localPosition=[0,0,0]))
 33A0 = np.eye(3)
 34Alast = A0 #previous marker
 35bodyLast = oGround
 36
 37Alist=[]
 38axisList=[]
 39A=RotationMatrixX(0)
 40
 41nBodies = 100
 42for i in range(nBodies):
 43    delta = 0.01*pi
 44    #A = RotationMatrixZ(delta)
 45    v0= A@[0,0,1]
 46    Alist+=[A]
 47    axisList+=[v0]
 48    A = A @ RotationMatrixX(delta)@RotationMatrixZ(2*delta)
 49
 50
 51p0 = [0.,0.,0] #reference position
 52vLoc = np.array([L,0,0]) #last to next joint
 53g = [0,0,9.81]
 54#g = [0,9.81,0]
 55
 56#create a chain of bodies:
 57for i in range(nBodies):
 58    #print("Build Object", i)
 59    inertia = InertiaCuboid(density=1000, sideLengths=[L,d,d])
 60    p0 += Alist[i] @ (0.5*vLoc)
 61    #p0 += (0.5*vLoc)
 62
 63    ep0 = eulerParameters0 #no rotation
 64    graphicsBody = graphics.Brick([0,0,0], [0.96*L,d,d], graphics.color.steelblue)
 65    oRB = mbs.CreateRigidBody(inertia=inertia,
 66                              referencePosition=p0,
 67                              referenceRotationMatrix=Alist[i],
 68                              gravity=g,
 69                              graphicsDataList=[graphicsBody])
 70    nRB= mbs.GetObject(oRB)['nodeNumber']
 71
 72    body0 = bodyLast
 73    body1 = oRB
 74    # point = mbs.GetObjectOutputBody(oRB,exu.OutputVariableType.Position,
 75    #                                 localPosition=[-0.5*L,0,0],
 76    #                                 configuration=exu.ConfigurationType.Reference)
 77    #axis = [0,0,1]
 78    axis = axisList[i]
 79    mbs.CreateRevoluteJoint(bodyNumbers=[body0, body1], position=[0.5*L,0,0],
 80                            axis=Alast.T@axis, useGlobalFrame=False,
 81                            axisRadius=0.6*d, axisLength=1.2*d)
 82    # mbs.CreateRevoluteJoint(bodyNumbers=[body0, body1], position=point,
 83    #                         axis=axis, useGlobalFrame=True,
 84    #                         axisRadius=0.6*d, axisLength=1.2*d)
 85
 86    bodyLast = oRB
 87
 88    p0 += Alist[i] @ (0.5*vLoc)
 89    #p0 += (0.5*vLoc)
 90    Alast = Alist[i]
 91
 92#mbs.AddLoad(LoadForceVector(markerNumber=mPosLast, loadVector=[0,0,20]))
 93
 94mbs.Assemble()
 95
 96simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 97
 98tEnd = 1
 99h=0.0005  #use small step size to detext contact switching
100
101simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
102simulationSettings.timeIntegration.endTime = tEnd
103simulationSettings.solutionSettings.solutionWritePeriod = 0.005
104simulationSettings.solutionSettings.sensorsWritePeriod = 0.01
105#simulationSettings.timeIntegration.simulateInRealtime = True
106simulationSettings.timeIntegration.realtimeFactor = 0.5
107simulationSettings.timeIntegration.verboseMode = 1
108
109simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
110simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations=True
111simulationSettings.timeIntegration.newton.useModifiedNewton = True
112#simulationSettings.timeIntegration.newton.modifiedNewtonJacUpdatePerStep = True
113simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
114# simulationSettings.parallel.numberOfThreads=4
115
116SC.visualizationSettings.nodes.show = True
117SC.visualizationSettings.nodes.drawNodesAsPoint  = False
118SC.visualizationSettings.nodes.showBasis = True
119SC.visualizationSettings.nodes.basisSize = 0.015
120SC.visualizationSettings.connectors.showJointAxes = True
121
122#for snapshot:
123SC.visualizationSettings.openGL.multiSampling=4
124SC.visualizationSettings.openGL.lineWidth=2
125SC.visualizationSettings.view0.window.renderWindowSize = [800,600]
126SC.visualizationSettings.view0.scene.drawCoordinateSystem=False
127SC.visualizationSettings.view0.scene.drawWorldBasis=True
128# SC.visualizationSettings.general.useMultiThreadedRendering = False
129SC.visualizationSettings.general.autoFitScene = False #use loaded render state
130
131
132#test UTF-8 characters:
133text = 'Demo UTF-8 text:ΓΔΘΛΞΠΣΦΨΩ\nαβγδεζηθικλμνξοπρστυφχψωϕϵ\n'
134text+= 'x₀₁₂₃₄₅₆₇₈₉x⁰¹²³⁴⁵⁶⁷⁸⁹\n∂∫♥√≈∞🙂😒°×·\nüöäÜÖÄßéèáàاÿ~'
135
136SC.visualizationSettings.general.renderWindowString = text
137SC.visualizationSettings.view0.window.globalFontSize = 14 #to see special characters
138useGraphics = True
139if useGraphics:
140    simulationSettings.displayComputationTime = True
141    simulationSettings.displayStatistics = True
142    SC.renderer.Start()
143    if 'renderState' in exu.sys:
144        SC.renderer.SetState(exu.sys[ 'renderState' ])
145    #SC.renderer.DoIdleTasks()
146else:
147    simulationSettings.solutionSettings.writeSolutionToFile = False
148
149#mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.TrapezoidalIndex2)
150mbs.SolveDynamic(simulationSettings, showHints=True)
151
152if True: #use this to reload the solution and use SolutionViewer
153    #sol = LoadSolutionFile('coordinatesSolution.txt')
154
155    mbs.SolutionViewer() #can also be entered in IPython ...
156
157
158u0 = mbs.GetNodeOutput(nRB, exu.OutputVariableType.Displacement)
159rot0 = mbs.GetNodeOutput(nRB, exu.OutputVariableType.Rotation)
160exu.Print('u0=',u0,', rot0=', rot0)
161
162result = (abs(u0)+abs(rot0)).sum()
163exu.Print('solution of addRevoluteJoint=',result)
164
165
166
167#%%+++++++++++++++++++++++++++++
168if useGraphics:
169    SC.renderer.DoIdleTasks()
170    SC.renderer.Stop() #safely close rendering window!