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:  A chain of 100 rigid bodies, each turned against the previous one, linked by revolute joints and
  6#           falling under gravity; the frames of the chain are homogeneous transformations (exu.HT): the
  7#           bodies take theirs as referenceHT, the joint markers as localHT. The solution is shown again
  8#           with the SolutionViewer.
  9#
 10# Author:   Johannes Gerstmayr 
 11# Date:     2021-07-01
 12#
 13# 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.
 14#
 15#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 16
 17import exudyn as exu
 18from exudyn.itemInterface import *
 19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 20import exudyn.graphics as graphics
 21
 22from math import sin, cos, pi
 23import numpy as np
 24
 25SC = exu.SystemContainer()
 26mbs = SC.AddSystem()
 27
 28
 29#background
 30color = [0.1,0.1,0.8,1]
 31L = 0.4 #length of bodies
 32d = 0.1 #diameter of bodies
 33
 34oGround=mbs.AddObject(ObjectGround(referenceHT=exu.HT(translation=[-0.5*L,0,0])))
 35nBodies = 100
 36delta = 0.01*pi
 37g = [0,0,9.81]
 38
 39#the frames of the chain are homogeneous transformations (exu.HT): the joint frame of each body is the joint
 40#frame of the previous one, moved along the body by L and turned by Hturn; the body sits halfway in between
 41Hhalf = exu.HT(translation=[0.5*L,0,0])                    #from a joint to the body center, and on to the next joint
 42Hturn = exu.HT().SetRotationX(delta) * exu.HT().SetRotationZ(2*delta) #the turn from one body to the next
 43Hjoint = exu.HT()                                          #the first joint at the origin; its z-axis is the joint axis
 44bodyLast = oGround
 45HturnLast = exu.HT()                                       #the ground is not turned against the first body
 46
 47for i in range(nBodies):
 48    inertia = InertiaCuboid(density=1000, sideLengths=[L,d,d])
 49    graphicsBody = graphics.Brick([0,0,0], [0.96*L,d,d], graphics.color.steelblue)
 50    oRB = mbs.CreateRigidBody(inertia=inertia,
 51                              referenceHT=Hjoint*Hhalf,
 52                              gravity=g,
 53                              graphicsDataList=[graphicsBody])
 54    nRB= mbs.GetObject(oRB)['nodeNumber']
 55
 56    #the joint frame in each body is the localHT of its marker: at the end of the previous body, turned to this
 57    #body, and at the start of this body
 58    mLast = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bodyLast, localHT=Hhalf*HturnLast))
 59    mThis = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRB, localHT=Hhalf.Inverse()))
 60    mbs.AddObject(RevoluteJointZ(markerNumbers=[mLast, mThis],
 61                                 visualization=VRevoluteJointZ(axisRadius=0.6*d, axisLength=1.2*d)))
 62
 63    bodyLast = oRB
 64    HturnLast = Hturn
 65    Hjoint = Hjoint*Hhalf*Hhalf*Hturn                      #the next joint frame
 66
 67
 68mbs.Assemble()
 69
 70simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 71
 72tEnd = 1
 73h=0.0005  #use small step size to detext contact switching
 74
 75simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
 76simulationSettings.timeIntegration.endTime = tEnd
 77simulationSettings.solution.file.writePeriod = 0.005
 78simulationSettings.solution.sensors.writePeriod = 0.01
 79#simulationSettings.timeIntegration.realtime.active = True
 80simulationSettings.timeIntegration.realtime.factor = 0.5
 81simulationSettings.timeIntegration.verboseMode = 1
 82
 83simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
 84simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations=True
 85simulationSettings.timeIntegration.newton.useModifiedNewton = True
 86#simulationSettings.timeIntegration.newton.modifiedNewtonJacUpdatePerStep = True
 87simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
 88# simulationSettings.parallel.numberOfThreads=4
 89
 90SC.visualizationSettings.nodes.show = True
 91SC.visualizationSettings.nodes.drawNodesAsPoint  = False
 92SC.visualizationSettings.nodes.showBasis = True
 93SC.visualizationSettings.nodes.basisSize = 0.015
 94SC.visualizationSettings.connectors.showJointAxes = True
 95
 96#for snapshot:
 97SC.visualizationSettings.openGL.multiSampling=4
 98SC.visualizationSettings.openGL.lineWidth=2
 99SC.visualizationSettings.view0.window.renderWindowSize = [800,600]
100SC.visualizationSettings.view0.scene.drawCoordinateSystem=False
101SC.visualizationSettings.view0.scene.drawWorldBasis=True
102# SC.visualizationSettings.general.useMultiThreadedRendering = False
103SC.visualizationSettings.general.autoFitScene = False #use loaded render state
104
105
106#test UTF-8 characters:
107text = 'Demo UTF-8 text:ΓΔΘΛΞΠΣΦΨΩ\nαβγδεζηθικλμνξοπρστυφχψωϕϵ\n'
108text+= 'x₀₁₂₃₄₅₆₇₈₉x⁰¹²³⁴⁵⁶⁷⁸⁹\n∂∫♥√≈∞🙂😒°×·\nüöäÜÖÄßéèáàاÿ~'
109
110SC.visualizationSettings.general.renderWindowString = text
111SC.visualizationSettings.view0.window.globalFontSize = 14 #to see special characters
112useGraphics = True
113if useGraphics:
114    simulationSettings.show.computationTime = True
115    simulationSettings.show.statistics = True
116    SC.renderer.Start()
117    SC.renderer.RestoreSavedState()
118    #SC.renderer.DoIdleTasks()
119else:
120    simulationSettings.solution.file.write = False
121
122#mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.TrapezoidalIndex2)
123mbs.SolveDynamic(simulationSettings, showHints=True)
124
125if True: #use this to reload the solution and use SolutionViewer
126    #sol = LoadSolutionFile(OutputFilePath('coordinatesSolution.txt'))
127    
128    mbs.SolutionViewer() #can also be entered in IPython ...
129
130
131u0 = mbs.GetNodeOutput(nRB, exu.OutputVariableType.Displacement)
132rot0 = mbs.GetNodeOutput(nRB, exu.OutputVariableType.Rotation)
133exu.Print('u0=',u0,', rot0=', rot0)
134
135result = (abs(u0)+abs(rot0)).sum()
136exu.Print('solution of addRevoluteJoint=',result)
137
138
139
140#%%+++++++++++++++++++++++++++++
141if useGraphics:
142    SC.renderer.DoIdleTasks()
143    SC.renderer.Stop() #safely close rendering window!