bricardMechanism.py
You can view and download this file on Github: bricardMechanism.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test models for mainSystemExtensions; tests for functions except PlotSensor or SolutionViewer,
5# which are tested already in other functions or cannot be tested with test suite;
6# all tests are self-contained and are included as examples for docu
7#
8# Author: Johannes Gerstmayr
9# Date: 2023-05-19
10#
11# 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.
12#
13#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
14
15import exudyn as exu
16from exudyn.utilities import * #includes itemInterface, graphicsDataUtilities and rigidBodyUtilities
17import numpy as np
18
19useGraphics = True #without test
20#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
21#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
22try: #only if called from test suite
23 from modelUnitTests import exudynTestGlobals #for globally storing test results
24 useGraphics = exudynTestGlobals.useGraphics
25except:
26 class ExudynTestGlobals:
27 pass
28 exudynTestGlobals = ExudynTestGlobals()
29#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
30
31endTime = 1 + 0*useGraphics*4 #test with only 1 second
32
33#%%++++++++++++++++++++++++++++++++++++++++++++++++++++
34#create rigid body with revolute joint:
35SC = exu.SystemContainer()
36mbs = SC.AddSystem()
37
38L = 1 #overall distance
39r = 0.1 #width of cubic bodies, rectangular cross section
40
41inertia = InertiaCuboid(density=100, sideLengths=[L,r,r])
42inertia = inertia.Translated([L*0.5,0,0])
43# exu.Print('m=',inertia.Mass())
44# exu.Print('COM=',inertia.COM())
45
46rotX = [1,0,0]
47rotY = [0,1,0]
48rotZ = [0,0,1]
49
50listP = [
51 [ 0, 0, 0],
52 [ L, 0, 0],
53 [ L,-L, 0],
54
55 [ L,-L, L],
56 [ 0,-L, L],
57 [ 0, 0, L],
58 ]
59rotList = [
60 np.eye(3),
61 RotationMatrixZ(-0.5*pi),
62 RotationMatrixY(-0.5*pi),
63
64 RotationMatrixZ(pi),
65 RotationMatrixZ(0.5*pi),
66 ]
67
68listRotAxes = [rotY,rotZ,rotX, rotY,rotZ,rotX,]
69
70gGround = [GraphicsDataCheckerBoard(point=[0,-2.1,0],normal=[0,1,0],size=6)]
71
72oGround = mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=gGround)))
73listBodies = [oGround]
74
75gData = [GraphicsDataOrthoCubePoint(centerPoint=[0.5*L,0,0],
76 size=[L,r,r], color=color4steelblue)]
77gData +=[GraphicsDataBasis(length = 0.25)]
78
79for i, p in enumerate(listP):
80 if i != 5:
81 b0 = mbs.CreateRigidBody(inertia = inertia,
82 referencePosition = p,
83 referenceRotationMatrix=rotList[i],
84 gravity = [0,-9.81,0],
85 graphicsDataList = gData)
86 else:
87 b0 = oGround
88
89 if False: #True works less good
90 mbs.CreateRevoluteJoint(bodyNumbers=[listBodies[-1], b0],
91 position=p,
92 axis=listRotAxes[i],
93 axisRadius=r, axisLength=1.1*r)
94 else:
95 #using one GenericJoint works slightly better in full Newton case than pure revolute joints
96 if i != 5:
97 mbs.CreateRevoluteJoint(bodyNumbers=[listBodies[-1], b0],
98 position=p,
99 axis=listRotAxes[i],
100 axisRadius=r, axisLength=1.1*r)
101 else:
102 mbs.CreateGenericJoint(bodyNumbers=[listBodies[-1], b0],
103 position=p,
104 constrainedAxes=[1,1,1, 0,1,1],
105 axesRadius=r, axesLength=1.1*r)
106
107 # # as this mechanism contains a redundant constraint and the standard solver cannot cope with that
108 # # we have to use a flexible joint instead
109 # rbd=mbs.CreateRigidBodySpringDamper(bodyOrNodeList=[listBodies[-1], b0],
110 # localPosition0=[ L,0,0],
111 # localPosition1=[ 0,0,L],
112 # stiffness=1e6*np.diag([1,1,1,0,1,1]),
113 # drawSize=r)
114
115 listBodies += [b0]
116
117mbs.Assemble()
118simulationSettings = exu.SimulationSettings() #takes currently set values or default values
119simulationSettings.solutionSettings.solutionWritePeriod = 0.02
120simulationSettings.solutionSettings.writeSolutionToFile = useGraphics
121simulationSettings.timeIntegration.numberOfSteps = 1000
122simulationSettings.timeIntegration.endTime = endTime
123simulationSettings.timeIntegration.verboseMode = 1
124simulationSettings.displayComputationTime = True
125simulationSettings.displayStatistics = True
126
127simulationSettings.timeIntegration.newton.useModifiedNewton = True
128simulationSettings.timeIntegration.newton.relativeTolerance = 1e-6
129
130#the dense solver can treat redundant constraints if according flags turned on
131simulationSettings.linearSolverType = exu.LinearSolverType.EigenDense
132simulationSettings.linearSolverSettings.ignoreSingularJacobian = True
133# simulationSettings.linearSolverSettings.pivotThreshold = 1e-10
134
135#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
136#simulation times for system size 65, last joint=RigidBodySpringDamper!:
137 # useModifiedNewton = False
138 # 10000 steps
139 # endTime=5
140 # EXUdense: tCPU=5.67 / factorization=55.1% / NewtonInc= 1.59% / factTime=3.124
141 # EigenDense / PartPivLU: tCPU=3.57 / factorization=34.7% / NewtonInc= 1.92% / factTime=1.239
142 # EigenDense / FullPivLU: tCPU=5.43 / factorization=55.6% / NewtonInc= 4.83% / factTime=3.019
143#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
144
145
146# SC.visualizationSettings.general.drawWorldBasis = True
147SC.visualizationSettings.openGL.shadow = 0.3
148SC.visualizationSettings.openGL.light0position = [2,12,3,0]
149SC.visualizationSettings.openGL.multiSampling = 4
150
151SC.visualizationSettings.general.autoFitScene = False #prevent from autozoom
152
153if useGraphics:
154 exu.StartRenderer()
155 if 'renderState' in exu.sys:
156 SC.SetRenderState(exu.sys['renderState'])
157 mbs.WaitForUserToContinue()
158
159dof=mbs.ComputeSystemDegreeOfFreedom()
160exu.Print('dof',dof)
161[eigenValues,x] = mbs.ComputeODE2Eigenvalues()
162exu.Print('eigenvalues=',eigenValues)
163
164mbs.SolveDynamic(simulationSettings = simulationSettings)
165
166if useGraphics:
167 exu.StopRenderer()
168
169if False:
170 #%%++++
171 mbs.SolutionViewer()
172
173#%%++++
174testError = np.linalg.norm(mbs.systemData.GetODE2Coordinates())
175testError += dof['degreeOfFreedom'] + dof['redundantConstraints'] + eigenValues[0]
176exu.Print('solution of bricardMechanism test=',testError)
177
178#%%++++++++++++++++++++++++++++++++++++++++++++++++++++
179
180exudynTestGlobals.testError = testError - (4.172189649307425) #2023-06-12: 4.172189649307425
181exudynTestGlobals.testResult = testError