generalContactCylinderTrigsTest.py

You can view and download this file on Github: generalContactCylinderTrigsTest.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  model with GeneralContact of a cylinder modelled by triangles rolling on a bed of spherical markers
  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
 15import exudyn.graphics as graphics #only import if it does not conflict
 16import numpy as np
 17
 18useGraphics = True #without test
 19#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 20#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 21try: #only if called from test suite
 22    from modelUnitTests import exudynTestGlobals #for globally storing test results
 23    useGraphics = exudynTestGlobals.useGraphics
 24except:
 25    class ExudynTestGlobals:
 26        pass
 27    exudynTestGlobals = ExudynTestGlobals()
 28#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 29
 30SC = exu.SystemContainer()
 31mbs = SC.AddSystem()
 32
 33#create an environment for mini example
 34
 35gravity = [0,0,-9.81]
 36
 37p0 = np.array([0,0,0])
 38planeL = 0.3
 39xFact = 5
 40nPlane = 6 #spheres along x and y
 41rMarker = 0.02 #radius of spheres for contact points
 42
 43rCyl = 0.25
 44tCyl = 0.25
 45
 46sRad = 0.02
 47k = 2e4
 48d = 0.0005*k
 49stepSize = 1e-3
 50
 51frictionCoeff = 0.1
 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 #does not influence results (only performance)!
 67
 68
 69#%% ground
 70gFloor = graphics.CheckerBoard(p0+[planeL*xFact*0.5,0,0],size=planeL*xFact, size2=planeL, nTiles=8)
 71oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
 72                                    visualization=VObjectGround(graphicsData=[gFloor])))
 73
 74
 75# nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0] ))
 76# mGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround))
 77
 78gCyl = graphics.Cylinder(pAxis=[0,-0.5*tCyl,0],vAxis=[0,tCyl,0], radius=rCyl,
 79                                       color=[0.3,0.3,0.3,1],
 80                                       alternatingColor=graphics.color.lightgrey,nTiles=64)
 81
 82
 83#rigid body containing sphere markers:
 84omegaY = 10
 85bCyl=mbs.CreateRigidBody(referencePosition=[0,0,rCyl*1.0],
 86                    initialVelocity=[omegaY*rCyl,0,0],
 87                    #initialAngularVelocity=[0,omegaY,0],
 88                    # referenceRotationMatrix=RotationMatrixX(0.1),
 89                    inertia=InertiaCylinder(200, tCyl, rCyl, 1),
 90                    gravity = gravity,
 91                    nodeType = exu.NodeType.RotationRotationVector,
 92                    graphicsDataList=[gCyl],
 93                    )
 94mCyl = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bCyl, localPosition=[0,0,0]))
 95sPos = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
 96                                outputVariableType=exu.OutputVariableType.Position))
 97sRot = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
 98                                outputVariableType=exu.OutputVariableType.Rotation))
 99sVel = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
100                                outputVariableType=exu.OutputVariableType.Velocity))
101sOmega = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
102                                outputVariableType=exu.OutputVariableType.AngularVelocity))
103
104[meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(gCyl)
105
106
107gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mCyl,
108                                    contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
109    pointList=meshPoints,  triangleList=meshTrigs)
110
111
112
113for ix in range(nPlane*xFact+1):
114    for iy in range(nPlane+1):
115        x = (ix/nPlane-0.0)*planeL
116        y = (iy/nPlane-0.5)*planeL
117        z = -rMarker
118
119        m = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[x,y,z]))
120        gContact.AddSphereWithMarker(m, radius=rMarker, contactStiffness=k, contactDamping=d,
121                                     frictionMaterialIndex=0)
122
123
124
125
126
127mbs.Assemble()
128
129items=gContact.GetItemsInBox(pMin=[-4,-4,0], pMax=[4,4,20])
130# print('n spheres=',len(items['MarkerBasedSpheres']))
131
132
133tEnd = 1
134simulationSettings = exu.SimulationSettings()
135simulationSettings.solutionSettings.writeSolutionToFile = False
136simulationSettings.timeIntegration.verboseMode = 1
137#simulationSettings.parallel.numberOfThreads = 4
138
139SC.visualizationSettings.general.graphicsUpdateInterval=0.02
140SC.visualizationSettings.general.drawCoordinateSystem=True
141SC.visualizationSettings.loads.show=False
142SC.visualizationSettings.bodies.show=True
143SC.visualizationSettings.markers.show=False
144
145SC.visualizationSettings.nodes.show=True
146SC.visualizationSettings.nodes.drawNodesAsPoint = False
147SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
148SC.visualizationSettings.nodes.tiling = 4
149
150SC.visualizationSettings.window.renderWindowSize=[1900,1200]
151SC.visualizationSettings.openGL.multiSampling = 4
152#improved OpenGL rendering
153
154SC.visualizationSettings.contact.showSpheres = True
155SC.visualizationSettings.contact.tilingSpheres = 4
156
157
158if useGraphics:
159    exu.StartRenderer()
160    if 'renderState' in exu.sys:
161        SC.SetRenderState(exu.sys['renderState'])
162    # mbs.WaitForUserToContinue()
163
164
165simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
166simulationSettings.timeIntegration.endTime = tEnd
167mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
168
169if useGraphics:
170    #SC.WaitForRenderEngineStopFlag()
171    exu.StopRenderer() #safely close rendering window!
172
173#%%+++++++++++++++++++
174q = mbs.GetSensorValues(sPos)
175q += mbs.GetSensorValues(sVel)
176q += mbs.GetSensorValues(sOmega)
177q += mbs.GetSensorValues(sRot)
178#print('q=', q)
179
180u = NormL2(q)
181exu.Print('solution of generalContactCylinderTest =',u)
182
183exudynTestGlobals.testError = u - (5.486908430912642)
184exudynTestGlobals.testResult = u