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