createSphereQuadContact2.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for CreateSphereQuadContact, combining sphere-sphere and 2 sphere-triangle contacts
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2025-06-29
  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 *
 15from math import sin, cos
 16import exudyn.graphics as graphics
 17import numpy as np
 18
 19testIsActive = exu.sys.get('testIsActive', False)
 20#testIsActive = True
 21
 22testSolution = 0
 23
 24SC = exu.SystemContainer()
 25mbs = SC.AddSystem()
 26
 27
 28radius=0.1
 29mass = 0.2                              #mass in kg
 30contactStiffness = 2e4*10               #stiffness of spring-damper in N/m
 31contactDamping = 1e-4*contactStiffness  #damping constant in N/(m/s); as we have always contact, we only need some damping for initial effects
 32dynamicFriction = 0.2
 33
 34isExplicitSolver = False
 35tEnd = 1     #end time of simulation
 36if not testIsActive:
 37    tEnd = 10
 38
 39stepSize = 5e-4 #*10
 40
 41g = 9.81
 42
 43size = 1
 44
 45oGround = mbs.CreateGround(graphicsDataList=[graphics.CheckerBoard(point=[0,0,0],size=2*size,
 46                                                                   color=graphics.color.lightgrey[0:3]+[graphics.material.indexChrome],
 47                                                                   alternatingColor=graphics.color.lightgrey2[0:3]+[graphics.material.indexChrome],
 48                                                                   ), ])
 49mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround))
 50mGround2 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,-radius]))
 51
 52#for explicit solver, we need Lie group node:
 53nodeType=exu.NodeType.RotationRotationVector if isExplicitSolver else exu.NodeType.RotationEulerParameters
 54
 55
 56nSpheres = 5
 57listMasses=[]
 58sPosList=[]
 59
 60for iPhi in range(nSpheres):
 61    
 62    phi = iPhi / nSpheres * (2*pi)
 63    x = 2*radius*cos(phi)
 64    y = 2*radius*sin(phi)
 65    
 66    vv = 0.5*2
 67    vx = -vv*sin(phi)
 68    vy = vv*cos(phi)
 69
 70    oMass = mbs.CreateRigidBody(referencePosition=[x,y,radius],
 71                                initialVelocity=[vx,vy,0],
 72                                initialAngularVelocity=[0,0,0],
 73                                nodeType=nodeType,
 74                                inertia=InertiaSphere(mass=mass, radius=radius),
 75                                gravity = [0,0,-g],
 76                                graphicsDataList=[graphics.Sphere(radius=radius, color=graphics.colorList[iPhi], 
 77                                                                  nTiles=32)],
 78                                )
 79    listMasses.append(oMass)
 80    mMass = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass))
 81
 82    if False:
 83        #create contact between each sphere:
 84        for oMass2 in listMasses[:-1]:
 85            restitutionCoefficient = 0.75
 86            impactModel = 2
 87            oSSC = mbs.CreateSphereSphereContact(itemNumbers=[oMass, oMass2],
 88                                                 spheresRadii=[radius, radius],
 89                                                 contactStiffness = contactStiffness,
 90                                                 dynamicFriction = dynamicFriction,
 91                                                 impactModel = impactModel,
 92                                                 restitutionCoefficient = restitutionCoefficient,
 93                                                 show=True,
 94                                                 )
 95
 96    quadPoints = exu.Vector3DList([[-size,-size,0],[size,-size,0],[size,size,0],[-size,size,0]])
 97    oSSC = mbs.CreateSphereQuadContact(itemNumbers=[oMass, oGround],
 98                                       quadPoints=quadPoints,
 99                                       includeEdges=15, #all edges
100                                       sphereRadius=radius,
101                                       contactStiffness = contactStiffness,
102                                       contactDamping = contactDamping,
103                                       dynamicFriction=dynamicFriction,
104                                       show = True
105                                       )
106    sPos=mbs.AddSensor(SensorBody(bodyNumber=oMass, storeInternal=True,
107                                  outputVariableType=exu.OutputVariableType.Position))
108    sPosList.append(sPos)
109
110
111sizeTop = 8*radius
112zTop = 0.5*radius
113oBodyTop = mbs.CreateRigidBody(referencePosition=[0,0,2*radius+0.5*zTop],
114                               initialVelocity=[0,0,0],
115                               initialAngularVelocity=[0,0,0], #not initialized and will follow initial velocity of spheres
116                               nodeType=nodeType,
117                               inertia=InertiaCuboid(0.01*2800, sideLengths=[sizeTop,sizeTop,zTop]),
118                               gravity = [0,0,-g],
119                               color = graphics.color.red[0:3]+[0.3],
120                            )
121sBodyAngVel = mbs.AddSensor(SensorBody(bodyNumber=oBodyTop, 
122                                    storeInternal=True, 
123                                    outputVariableType=exu.OutputVariableType.AngularVelocity))
124# mbs.CreateGenericJoint(itemNumbers=[oGround, oBodyTop], position = [0,0,2*radius+0.5*zTop])
125
126
127for oMass in listMasses:
128    quadPoints = exu.Vector3DList([[-0.5*sizeTop,-0.5*sizeTop,-0.5*zTop],
129                                   [0.5*sizeTop,-0.5*sizeTop,-0.5*zTop],
130                                   [0.5*sizeTop,0.5*sizeTop,-0.5*zTop],
131                                   [-0.5*sizeTop,0.5*sizeTop,-0.5*zTop]])
132    oSSC = mbs.CreateSphereQuadContact(itemNumbers=[oMass,oBodyTop],
133                                       quadPoints=quadPoints,
134                                       includeEdges=15, #all edges
135                                       sphereRadius=radius*1,
136                                       contactStiffness = contactStiffness,
137                                       contactDamping = contactDamping,
138                                       dynamicFriction=dynamicFriction,
139                                       show = True
140                                       )
141
142
143#exu.Print(mbs)
144mbs.Assemble()
145
146simulationSettings = exu.SimulationSettings()
147simulationSettings.solution.file.write = (not testIsActive) #only the SolutionViewer reads it (#2492)
148simulationSettings.solution.file.writePeriod = 0.005
149simulationSettings.solution.sensors.writePeriod = 0.001  #output interval
150simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
151simulationSettings.timeIntegration.endTime = tEnd
152#simulationSettings.timeIntegration.realtime.active = True
153simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-6
154simulationSettings.timeIntegration.newton.relativeTolerance = 1e-6
155
156simulationSettings.timeIntegration.stepInformation = 3 #remove flag 64 which shows step reduction warnings
157
158
159simulationSettings.timeIntegration.newton.useModifiedNewton = True
160simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
161
162simulationSettings.show.statistics = True
163simulationSettings.timeIntegration.verboseMode = 1
164SC.visualizationSettings.view0.scene.drawCoordinateSystem = False
165SC.visualizationSettings.general.showSolverInformation = False
166SC.visualizationSettings.connectors.showContact = False
167SC.visualizationSettings.nodes.showBasis = True
168SC.visualizationSettings.nodes.basisSize = 1.5*radius
169#SC.visualizationSettings.openGL.advanced.depthSorting = True
170
171if not testIsActive:
172    SC.renderer.Start()              #start graphics visualization
173    SC.renderer.DoIdleTasks()    #wait for pressing SPACE bar to continue
174
175mbs.SolveDynamic(simulationSettings)
176
177if not testIsActive:
178    SC.renderer.Stop()               #safely close rendering window!
179
180    if False:
181        mbs.PlotSensor(sPosList, components=[2]*len(sPosList))
182    
183#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
184ode2 = mbs.systemData.GetODE2Coordinates()
185
186testSolution += 0.1*np.linalg.norm(ode2)
187
188
189#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
190exu.Print('solution of createSphereQuadContact2=',testSolution) 
191exu.sys['testResult'] = testSolution
192#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
193
194mbs.PlotSensor(sBodyAngVel,components=2)
195
196
197if not testIsActive and False:
198    mbs.SolutionViewer()