generalContactCylinderTest.py
You can view and download this file on Github: generalContactCylinderTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: model with GeneralContact of a cylinder modelled by spheres rolling on triangle ground mesh
5#
6# Author: Johannes Gerstmayr
7# Date: 2024-03-16
8#
9# 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.
10#
11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12
13import exudyn as exu
14from exudyn.utilities import *
15import numpy as np
16
17useGraphics = True #without test
18#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
19#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
20try: #only if called from test suite
21 from modelUnitTests import exudynTestGlobals #for globally storing test results
22 useGraphics = exudynTestGlobals.useGraphics
23except:
24 class ExudynTestGlobals:
25 pass
26 exudynTestGlobals = ExudynTestGlobals()
27#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
28
29SC = exu.SystemContainer()
30mbs = SC.AddSystem()
31
32#create an environment for mini example
33
34gravity = [0,0,-9.81]
35
36planeL = 4
37p0 = [0.5*planeL,0,0]
38
39nPhi = 80 #spheres around circumference
40nThick = 1 #spheres along cylinder axis (+1)
41
42rCyl = 0.25
43tCyl = 0.1
44rMarker = 0.02 #radius of spheres for contact points
45
46sRad = 0.02
47k = 1e4
48d = 0.0005*k
49stepSize = 1e-3
50
51frictionCoeff = 0.2
52ss = 10
53
54markerList = []
55radiusList = []
56gDataList = []
57
58
59gContact = mbs.AddGeneralContact()
60#gContact.verboseMode = 1
61gContact.SetFrictionPairings(frictionCoeff*np.eye(1))
62gContact.SetSearchTreeCellSize(numberOfCells=[ss,ss,ss])
63# gContact.SetSearchTreeBox(pMin=np.array([-0.5*planeL,-0.5*planeL,-0.1]),
64# pMax=np.array([0.5*planeL,0.5*planeL,0.1]))
65#print('treesize=',ssx*ssx*ssy)
66# gContact.sphereSphereContact = False
67
68#%% ground
69gFloor = GraphicsDataCheckerBoard(p0,size=planeL, nTiles=10)
70gDataList = [gFloor]
71
72
73nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0] ))
74mGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround))
75
76[meshPoints, meshTrigs] = GraphicsData2PointsAndTrigs(gFloor)
77#[meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
78gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mGround, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
79 pointList=meshPoints, triangleList=meshTrigs)
80
81
82#rigid body containing sphere markers:
83inertia = InertiaCylinder(1000, tCyl, rCyl, 1)
84# print(inertia)
85
86omegaY = 12
87bCyl=mbs.CreateRigidBody(referencePosition=[0,0,rCyl*1.0],
88 initialVelocity=[omegaY*rCyl,0,0],
89 initialAngularVelocity=[0,omegaY*0.5,0],
90 initialRotationMatrix=RotationMatrixX(0.1),
91 inertia=inertia,
92 gravity = gravity,
93 nodeType = exu.NodeType.RotationRotationVector,
94 graphicsDataList=[GraphicsDataCylinder(pAxis=[0,-0.5*tCyl,0],vAxis=[0,tCyl,0], radius=rCyl,
95 color=[0.3,0.3,0.3,1],
96 alternatingColor=color4lightgrey,nTiles=64)]
97 )
98nCyl = mbs.GetObject(bCyl)['nodeNumber']
99#print(mbs.GetNode(nCyl))
100sPos = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
101 outputVariableType=exu.OutputVariableType.Position))
102sRot = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
103 outputVariableType=exu.OutputVariableType.Rotation))
104sVel = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
105 outputVariableType=exu.OutputVariableType.Velocity))
106sOmega = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
107 outputVariableType=exu.OutputVariableType.AngularVelocity))
108
109for phiI in range(nPhi):
110 phi = phiI/nPhi*2*pi
111 for j in range(nThick+1):
112 rCylMod = rCyl-rMarker
113 #compute local coordinates for markers
114 y = (j/nThick-0.5)*tCyl
115 x = rCylMod*sin(phi)
116 z = rCylMod*cos(phi)
117
118 m = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bCyl, localPosition=[x,y,z]))
119 gContact.AddSphereWithMarker(m, radius=rMarker, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
120
121
122
123#put ground here, such that it is transparent in background
124oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
125 visualization=VObjectGround(graphicsData=gDataList)))
126
127
128mbs.Assemble()
129
130items=gContact.GetItemsInBox(pMin=[-4,-4,0], pMax=[4,4,20])
131#print('n spheres=',len(items['MarkerBasedSpheres']))
132
133
134tEnd = 2
135#tEnd = h*100
136simulationSettings = exu.SimulationSettings()
137simulationSettings.solutionSettings.writeSolutionToFile = False
138#simulationSettings.displayComputationTime = True
139#simulationSettings.displayStatistics = True
140simulationSettings.timeIntegration.verboseMode = 1
141# simulationSettings.parallel.numberOfThreads = 4
142# simulationSettings.timeIntegration.simulateInRealtime = True
143
144SC.visualizationSettings.general.graphicsUpdateInterval=0.02
145SC.visualizationSettings.general.drawCoordinateSystem=True
146SC.visualizationSettings.loads.show=False
147SC.visualizationSettings.bodies.show=True
148SC.visualizationSettings.markers.show=False
149
150SC.visualizationSettings.nodes.show=True
151SC.visualizationSettings.nodes.drawNodesAsPoint = False
152SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
153SC.visualizationSettings.nodes.tiling = 4
154
155SC.visualizationSettings.window.renderWindowSize=[2000,1200]
156SC.visualizationSettings.openGL.multiSampling = 4
157#improved OpenGL rendering
158
159SC.visualizationSettings.contact.showSpheres = True
160
161SC.visualizationSettings.general.autoFitScene = False
162
163if useGraphics:
164 exu.StartRenderer()
165 if 'renderState' in exu.sys:
166 SC.SetRenderState(exu.sys['renderState'])
167 # mbs.WaitForUserToContinue()
168
169
170simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
171simulationSettings.timeIntegration.endTime = tEnd
172mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
173
174if useGraphics:
175 #SC.WaitForRenderEngineStopFlag()
176 exu.StopRenderer() #safely close rendering window!
177
178#%%+++++++++++++++++++
179q = mbs.GetSensorValues(sPos)
180q += mbs.GetSensorValues(sVel)
181q += mbs.GetSensorValues(sOmega)
182q += mbs.GetSensorValues(sRot)
183#print('q=', q)
184
185u = NormL2(q)
186exu.Print('solution of generalContactCylinderTest =',u)
187
188exudynTestGlobals.testError = u - (12.42377622187738 )
189exudynTestGlobals.testResult = u