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