generalContactImplicit2.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  test implicit integration of contact with spheres and triangles, focus on normal contact
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2024-06-02
  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
 16
 17
 18import numpy as np
 19
 20testIsActive = exu.sys.get('testIsActive', False)
 21
 22SC = exu.SystemContainer()
 23mbs = SC.AddSystem()
 24
 25doImplicit = True
 26#%%+++++++++++++++++++++++++++++++++
 27#sphere-sphere with coordinate constraints, prestressed; fixed torque on one side, linear increasing torque on other side
 28#sphere on ground, rolling
 29#cube on ground, sliding (f=[f, f*mu*t, 0]), tangential force changing
 30#cube on ground with initial velocity
 31#cube-cube contact (meshed)
 32
 33
 34L = 1 #surrounding
 35a = 0.1 #base dimention of objects
 36r = 0.5*a #radius
 37t = 0.25*a #thickness
 38
 39#contact coefficients:
 40mu = 0.1      #dry friction
 41m = 0.025     #mass
 42k = 1e3       #(linear) normal contact stiffness
 43d = 2e-4*k  #(linear) contact damping
 44gFact = 10
 45g = [0,0,-gFact]
 46
 47gContact = mbs.AddGeneralContact()
 48gContact.verboseMode = 1
 49#gContact.sphereSphereContact = False
 50gContact.frictionProportionalZone = 1e-3*0.1
 51#gContact.excludeDuplicatedTrigSphereContactPoints = False
 52fricMat = mu*np.eye(1)
 53gContact.SetFrictionPairings(fricMat)
 54gContact.SetSearchTreeCellSize(numberOfCells=[4,4,4])
 55
 56#%% ground
 57p0 = np.array([0,0,-0.5*t])
 58color4wall = [0.9,0.9,0.7,0.5]
 59addNormals = False
 60gFloor = graphics.Brick(p0,[L,L,t],graphics.color.steelblue,addNormals)
 61
 62gDataList = [gFloor]
 63
 64
 65nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0] ))
 66mGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround))
 67mGroundC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
 68
 69[meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(gFloor)
 70#[meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
 71# [meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
 72gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mGround, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
 73    pointList=meshPoints,  triangleList=meshTrigs)
 74
 75if True: #looses color
 76    gFloor = graphics.FromPointsAndTrigs(meshPoints, meshTrigs, color=color4wall) #show refined mesh
 77    gDataList = [gFloor]
 78
 79evalNodes = [] #collect nodes that are evaluated for test
 80#%%++++++++++++++++++++++++++++++++++++++++++++
 81#free rolling sphere:
 82gList = [graphics.Sphere(point=[0,0,0], radius=r, color= graphics.color.red, nTiles=24)]
 83gList += [graphics.Basis(length=2*r)]
 84
 85pRef = [-0.4*L,0*L,r - m*gFact/(0.5*k)] #need to use series of 2 k's!
 86
 87RBinertia = InertiaSphere(m, r)
 88oMass = mbs.CreateRigidBody(inertia=RBinertia, 
 89                            #nodeType=exu.NodeType.RotationRotationVector,
 90                            nodeType=exu.NodeType.RotationRxyz,
 91                            referencePosition=pRef,
 92                            #referenceRotationMatrix=RotXYZ2RotationMatrix([1,2,3]),
 93                            initialVelocity=[1,0,0],
 94                            gravity=g, 
 95                            graphicsDataList=gList,
 96                            )
 97
 98nMass = mbs.GetObject(oMass)['nodeNumber']
 99mbs.SetNodeParameter(nMass, 'Vshow', True) #node not shown with CreateRigidBody ...
100mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
101
102gContact.AddSphereWithMarker(mNode, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
103
104N = m*gFact
105F = mu*N
106exu.Print('contact force N =',N)
107exu.Print('friction force F=',F)
108exu.Print('acceleration    =',F/m)
109
110evalNodes += [nMass] 
111
112#add as last because of transparency
113oGround = mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=gDataList)))
114
115#only implicit:
116# mbs.CreateGenericJoint(itemNumbers=[oGround, oMass], position=pRef,
117#                        constrainedAxes=[0,0,0, 1,1,1])
118
119mRot0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=3))
120mRot1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=4))
121mRot2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=5))
122
123ccR0 = mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundC, mRot0]))
124ccR1 = mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundC, mRot1]))
125ccR2 = mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundC, mRot2]))
126
127if not testIsActive:
128    sNode0 = mbs.AddSensor(SensorNode(nodeNumber=nMass, storeInternal=True, 
129                                      outputVariableType=exu.OutputVariableType.Displacement))
130    vNode0 = mbs.AddSensor(SensorNode(nodeNumber=nMass, storeInternal=True, 
131                                      outputVariableType=exu.OutputVariableType.Velocity))
132    # omegaNode0 = mbs.AddSensor(SensorNode(nodeNumber=nMass, storeInternal=True, 
133    #                                   outputVariableType=exu.OutputVariableType.AngularVelocity))
134    
135    #lagrange multiplier; only available with implicit integration
136    sRot1 = mbs.AddSensor(SensorObject(objectNumber=ccR1, storeInternal=True, 
137                                      outputVariableType=exu.OutputVariableType.Force))
138
139#%%+++++++++++++++++++++++++++++++++
140mbs.Assemble()
141
142tEnd = 2 #tEnd = 0.8 for test suite
143h= 0.0002  #h= 0.0002 for test suite
144# h*=0.1
145# tEnd*=3
146simulationSettings = exu.SimulationSettings()
147#simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
148simulationSettings.solution.file.write = False
149if not testIsActive:
150    simulationSettings.solution.file.writePeriod = 0.005
151    simulationSettings.solution.file.write = True
152    simulationSettings.solution.file.name = 'solution/coordinatesSolution.txt'
153else:
154    simulationSettings.solution.file.export.accelerations = False
155    simulationSettings.solution.file.export.velocities = False
156    
157simulationSettings.solution.sensors.writePeriod = h*10
158simulationSettings.solution.precision = 8 #make files smaller
159# simulationSettings.show.computationTime = True
160# simulationSettings.show.globalTimers = True
161#simulationSettings.show.statistics = True
162simulationSettings.timeIntegration.verboseMode = 1
163simulationSettings.parallel.numberOfThreads = 1
164
165simulationSettings.timeIntegration.newton.numericalDifferentiation.forODE2 = False
166simulationSettings.timeIntegration.newton.useModifiedNewton = True
167
168SC.visualizationSettings.general.graphicsUpdateInterval=0.05
169# SC.visualizationSettings.view0.scene.drawWorldBasis = True
170SC.visualizationSettings.general.circleTiling=200
171SC.visualizationSettings.view0.scene.drawCoordinateSystem=True
172SC.visualizationSettings.loads.show=False
173SC.visualizationSettings.bodies.show=True
174SC.visualizationSettings.markers.show=False
175
176SC.visualizationSettings.nodes.show=True
177SC.visualizationSettings.nodes.showBasis = True
178SC.visualizationSettings.nodes.basisSize = 1
179SC.visualizationSettings.nodes.drawNodesAsPoint = False
180SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
181SC.visualizationSettings.nodes.tiling = 4
182SC.visualizationSettings.openGL.drawFaceNormals = False
183
184SC.visualizationSettings.openGL.multiSampling = 4
185SC.visualizationSettings.openGL.light0.shadow = 0.25
186SC.visualizationSettings.openGL.light0.position = [-3,3,10,0]
187
188if not testIsActive:
189    SC.visualizationSettings.general.autoFitScene = False
190    SC.renderer.Start()
191    SC.renderer.RestoreSavedState()
192    # SC.renderer.DoIdleTasks()
193
194simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
195simulationSettings.timeIntegration.endTime = tEnd
196simulationSettings.timeIntegration.explicit.computeEndOfStepAccelerations = False #increase performance, accelerations less accurate
197#mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
198mbs.SolveDynamic(simulationSettings)
199
200#compute error:
201uSum=0
202for node in evalNodes:
203    u = mbs.GetNodeOutput(node, exu.OutputVariableType.Coordinates)
204    exu.Print('coords node'+str(node)+' =',u)
205    for c in u:
206        uSum += abs(c) #add up all coordinates for comparison
207
208
209exu.Print('solution of generalContactImplicit1=',uSum)
210
211exu.sys['testResult'] = uSum
212
213    
214if not testIsActive:
215    # SC.renderer.DoIdleTasks()
216
217    if True:
218        SC.visualizationSettings.general.autoFitScene = False
219        SC.visualizationSettings.general.graphicsUpdateInterval=0.02
220        
221        mbs.SolutionViewer(timeout=0.01)
222
223    SC.renderer.Stop() #safely close rendering window!
224
225if not testIsActive:
226    mbs.PlotSensor([], closeAll=True)
227    mbs.PlotSensor([sNode0]*3, [0,1,2], figureName='sphere position')
228    mbs.PlotSensor([vNode0]*3, [0,1,2], figureName='sphere velocity')
229    mbs.PlotSensor([sRot1], [0], figureName='Torque y') #only available in implicit integrators!