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 * #includes itemInterface and rigidBodyUtilities
 15from math import sin, cos
 16import exudyn.graphics as graphics
 17import numpy as np
 18
 19testIsActive = exu.sys.get('testIsActive', False)
 20
 21SC = exu.SystemContainer()
 22mbs = SC.AddSystem()
 23
 24#create an environment for mini example
 25
 26gravity = [0,0,-9.81]
 27
 28planeL = 4
 29p0 = [0.5*planeL,0,0]
 30
 31nPhi = 80  #spheres around circumference
 32nThick = 1 #spheres along cylinder axis (+1)
 33
 34rCyl = 0.25
 35tCyl = 0.1
 36rMarker = 0.02 #radius of spheres for contact points
 37
 38sRad = 0.02
 39k = 1e4 
 40d = 0.0005*k
 41stepSize = 1e-3
 42
 43frictionCoeff = 0.2
 44ss = 10
 45
 46markerList = []
 47radiusList = []
 48gDataList = []
 49
 50
 51gContact = mbs.AddGeneralContact()
 52#gContact.verboseMode = 1
 53gContact.SetFrictionPairings(frictionCoeff*np.eye(1))
 54gContact.SetSearchTreeCellSize(numberOfCells=[ss,ss,ss])
 55# gContact.SetSearchTreeBox(pMin=np.array([-0.5*planeL,-0.5*planeL,-0.1]), 
 56#                           pMax=np.array([0.5*planeL,0.5*planeL,0.1]))
 57#print('treesize=',ssx*ssx*ssy)
 58# gContact.sphereSphereContact = False
 59
 60#%% ground
 61gFloor = graphics.CheckerBoard(p0,size=planeL, nTiles=10)
 62gDataList = [gFloor]
 63
 64
 65nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0] ))
 66mGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround))
 67
 68[meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(gFloor)
 69#[meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
 70gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mGround, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
 71    pointList=meshPoints,  triangleList=meshTrigs)
 72
 73
 74#rigid body containing sphere markers:
 75inertia = InertiaCylinder(1000, tCyl, rCyl, 1)
 76# print(inertia)
 77
 78omegaY = 12
 79bCyl=mbs.CreateRigidBody(referencePosition=[0,0,rCyl*1.0],
 80                    initialVelocity=[omegaY*rCyl,0,0],
 81                    initialAngularVelocity=[0,omegaY*0.5,0],
 82                    initialRotationMatrix=RotationMatrixX(0.1),
 83                    inertia=inertia,
 84                    gravity = gravity,
 85                    nodeType = exu.NodeType.RotationRotationVector,
 86                    graphicsDataList=[graphics.Cylinder(pAxis=[0,-0.5*tCyl,0],vAxis=[0,tCyl,0], radius=rCyl, 
 87                                                           color=[0.3,0.3,0.3,1],
 88                                                           alternatingColor=graphics.color.lightgrey,nTiles=64)]
 89                    )
 90nCyl = mbs.GetObject(bCyl)['nodeNumber']
 91#print(mbs.GetNode(nCyl))
 92sPos = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True, 
 93                                outputVariableType=exu.OutputVariableType.Position))
 94sRot = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True, 
 95                                outputVariableType=exu.OutputVariableType.Rotation))
 96sVel = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True, 
 97                                outputVariableType=exu.OutputVariableType.Velocity))
 98sOmega = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True, 
 99                                outputVariableType=exu.OutputVariableType.AngularVelocity))
100
101for phiI in range(nPhi):
102    phi = phiI/nPhi*2*pi
103    for j in range(nThick+1):
104        rCylMod = rCyl-rMarker
105        #compute local coordinates for markers
106        y = (j/nThick-0.5)*tCyl
107        x = rCylMod*sin(phi)
108        z = rCylMod*cos(phi)
109    
110        m = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bCyl, localPosition=[x,y,z]))
111        gContact.AddSphereWithMarker(m, radius=rMarker, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
112    
113
114
115#put ground here, such that it is transparent in background
116oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
117                                    visualization=VObjectGround(graphicsData=gDataList)))
118
119
120mbs.Assemble()
121
122items=gContact.GetItemsInBox(pMin=[-4,-4,0], pMax=[4,4,20])
123#print('n spheres=',len(items['MarkerBasedSpheres'])) 
124
125
126tEnd = 2
127#tEnd = h*100
128simulationSettings = exu.SimulationSettings()
129simulationSettings.solution.file.write = False
130#simulationSettings.show.computationTime = True
131#simulationSettings.show.statistics = True
132simulationSettings.timeIntegration.verboseMode = 1
133#simulationSettings.parallel.numberOfThreads = 4
134
135SC.visualizationSettings.general.graphicsUpdateInterval=0.02
136SC.visualizationSettings.view0.scene.drawCoordinateSystem=True
137SC.visualizationSettings.loads.show=False
138SC.visualizationSettings.bodies.show=True
139SC.visualizationSettings.markers.show=False
140
141SC.visualizationSettings.nodes.show=True
142SC.visualizationSettings.nodes.drawNodesAsPoint = False
143SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
144SC.visualizationSettings.nodes.tiling = 4
145
146SC.visualizationSettings.view0.window.renderWindowSize=[2000,1200]
147SC.visualizationSettings.openGL.multiSampling = 4
148#improved OpenGL rendering
149
150SC.visualizationSettings.contact.showSpheres = True
151
152SC.visualizationSettings.general.autoFitScene = False
153
154if not testIsActive:
155    SC.renderer.Start()
156    SC.renderer.RestoreSavedState()
157    # SC.renderer.DoIdleTasks()
158
159
160simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
161simulationSettings.timeIntegration.endTime = tEnd
162mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
163
164if not testIsActive:
165    #SC.renderer.DoIdleTasks()
166    SC.renderer.Stop() #safely close rendering window!
167
168#%%+++++++++++++++++++
169q = mbs.GetSensorValues(sPos)
170q += mbs.GetSensorValues(sVel)
171q += mbs.GetSensorValues(sOmega)
172q += mbs.GetSensorValues(sRot)
173#print('q=', q)
174
175u = np.linalg.norm(q)
176exu.Print('solution of generalContactCylinderTest =',u)
177
178exu.sys['testResult'] = u