sphereTriangleTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for SphereTrigContact
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2025-06-14
  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 exudyn.graphics as graphics
 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
 30
 31SC = exu.SystemContainer()
 32mbs = SC.AddSystem()
 33
 34radius=0.1
 35mass = 1.6                              #mass in kg
 36contactStiffness = 1e5                  #stiffness of spring-damper in N/m
 37contactDamping = 0*5e-4*contactStiffness  #damping constant in N/(m/s)
 38dynamicFriction = 0.2
 39restitutionCoefficient = 0.5
 40impactModel = 2
 41
 42tEnd = 0.5     #end time of simulation
 43stepSize = 2e-4 #*10
 44g = 10
 45
 46a = 0.25
 47
 48
 49oGround = mbs.CreateGround(graphicsDataList=[#graphics.CheckerBoard(point=[0,0,-0.001],size=4),
 50                                             graphics.Quad([[-a,-a,0],[2*a,-a,0],[2*a,2*a,0],[-a,2*a,0]],
 51                                                           color=graphics.color.dodgerblue)
 52                                             ])
 53mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround))
 54mGround2 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,-radius]))
 55
 56listMasses=[]
 57sPosList=[]
 58nMasses=6
 59for i in range(nMasses):
 60    #ground node
 61    vx = sin(i/nMasses*2*pi)
 62    vy = cos(i/nMasses*2*pi)
 63    oMass = mbs.CreateRigidBody(referencePosition=[2*radius*vx+0.2,2*radius*vy+0.2,radius*1.1],
 64                                    initialVelocity=[1.2*vx,1.2*vy,0],
 65                                    # initialAngularVelocity=[5*pi*2,0,pi*3*0],
 66                                    inertia=InertiaSphere(mass=mass, radius=radius),
 67                                    gravity = [0,0,-g],
 68                                    graphicsDataList=[graphics.Sphere(radius=radius,
 69                                                                      color=graphics.color.orange, nTiles=32)],
 70                                    )
 71    listMasses.append(oMass)
 72    mMass = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass))
 73
 74    trianglePoints0 = exu.Vector3DList([[-a,-a,0],[2*a,-a,0],[-a,2*a,0]])
 75    trianglePoints1 = exu.Vector3DList([[2*a,-a,0],[2*a,2*a,0],[-a,2*a,0]])
 76    includeEdgesList = [5,3] #watch difference to includeEdgesList = [0,0] with g=100!
 77    trigList = [trianglePoints0,trianglePoints1]
 78    for k, trianglePoints in enumerate(trigList):
 79        nData1 = mbs.AddNode(NodeGenericData(initialCoordinates=[0.1,0,0,0],
 80                                            numberOfDataCoordinates=4))
 81        oSSC = mbs.AddObject(ObjectContactSphereTriangle(markerNumbers=[mMass, mGround],
 82                                                        nodeNumber=nData1,
 83                                                        trianglePoints=trianglePoints,
 84                                                        includeEdges=includeEdgesList[k],
 85                                                        radiusSphere=radius,
 86                                                        contactStiffness = contactStiffness,
 87                                                        dynamicFriction=dynamicFriction,
 88                                                        contactDamping = contactDamping,
 89                                                        impactModel = impactModel,
 90                                                        # frictionProportionalZone = 0.001,
 91                                                        restitutionCoefficient = restitutionCoefficient,
 92                                                        #minimumImpactVelocity = 1e-2,
 93                                                        visualization=VObjectContactSphereSphere(show=True),
 94                                                        ))
 95
 96    sPos=mbs.AddSensor(SensorBody(bodyNumber=oMass, storeInternal=True,
 97                                    outputVariableType=exu.OutputVariableType.Position))
 98    sPosList.append(sPos)
 99
100mbs.Assemble()
101
102simulationSettings = exu.SimulationSettings()
103simulationSettings.solutionSettings.writeSolutionToFile = True
104simulationSettings.solutionSettings.solutionWritePeriod = 0.005
105simulationSettings.solutionSettings.sensorsWritePeriod = stepSize  #output interval
106simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
107simulationSettings.timeIntegration.endTime = tEnd
108#simulationSettings.timeIntegration.simulateInRealtime = True
109simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-8
110simulationSettings.timeIntegration.newton.relativeTolerance = 1e-6
111
112simulationSettings.timeIntegration.newton.useModifiedNewton = True
113simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
114
115simulationSettings.displayStatistics = True
116simulationSettings.timeIntegration.verboseMode = 1
117
118SC.visualizationSettings.openGL.multiSampling = 1
119SC.visualizationSettings.openGL.shadow = 0.2
120SC.visualizationSettings.openGL.depthSorting = True
121
122SC.visualizationSettings.window.renderWindowSize=[1600,1200]
123SC.visualizationSettings.nodes.showBasis = True
124SC.visualizationSettings.nodes.basisSize = radius*1.5
125
126if useGraphics:
127    SC.renderer.Start()              #start graphics visualization
128    SC.renderer.DoIdleTasks()    #wait for pressing SPACE bar to continue
129
130mbs.SolveDynamic(simulationSettings,
131                 solverType=exu.DynamicSolverType.TrapezoidalIndex2)
132
133if useGraphics:
134    SC.renderer.DoIdleTasks()
135    SC.renderer.Stop()               #safely close rendering window!
136
137    mbs.PlotSensor(sPosList, components=[2]*len(sPosList))
138
139    mbs.SolutionViewer()
140
141
142#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
143ode2 = mbs.systemData.GetODE2Coordinates()
144
145u = np.linalg.norm(ode2)
146exu.Print('solution of sphereTriangleTest=',u)
147
148exudynTestGlobals.testResult = u
149#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++