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