generalContactImplicit1.py

You can view and download this file on Github: generalContactImplicit1.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.8      #dry friction
 41m = 0.025     #mass
 42k = 1e3       #(linear) normal contact stiffness
 43d = 2*1e-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
 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
 85omega0 = -0.*np.array([5,1.,0.])
 86pRef = [-0.1*L,-0.1*L,0.5*L+r-0*m*gFact/k]
 87RBinertia = InertiaSphere(m, r)
 88oMass = mbs.CreateRigidBody(inertia=RBinertia, 
 89                            nodeType=exu.NodeType.RotationRotationVector,
 90                            referencePosition=pRef,
 91                            #referenceRotationMatrix=RotationMatrixX(0.),
 92                            initialVelocity=-np.cross([0,0,r],omega0),
 93                            initialAngularVelocity=omega0, 
 94                            gravity=g, 
 95                            graphicsDataList=gList,
 96                            )
 97
 98nMass = mbs.GetObject(oMass)['nodeNumber']
 99mbs.SetNodeParameter(nMass, 'Vshow', True) #node not shown with CreateRigidBody ...
100
101nNode0 = nMass
102mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
103
104gContact.AddSphereWithMarker(mNode, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
105if not testIsActive:
106    sNode0 = mbs.AddSensor(SensorNode(nodeNumber=nNode0, storeInternal=True, #fileName='solution/contactNode0.txt',
107                                      outputVariableType=exu.OutputVariableType.Displacement))
108    vNode0 = mbs.AddSensor(SensorNode(nodeNumber=nNode0, storeInternal=True, #fileName='solution/contactNode0Vel.txt',
109                                      outputVariableType=exu.OutputVariableType.Velocity))
110    omegaNode0 = mbs.AddSensor(SensorNode(nodeNumber=nNode0, storeInternal=True, #fileName='solution/contactNode0Vel.txt',
111                                      outputVariableType=exu.OutputVariableType.AngularVelocity))
112evalNodes += [nMass] 
113
114
115#add as last because of transparency
116oGround = mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=gDataList)))
117
118#%%+++++++++++++++++++++++++++++++++
119mbs.Assemble()
120
121tEnd = 2 #tEnd = 0.8 for test suite
122h= 0.0002  #h= 0.0002 for test suite
123# h*=0.1
124# tEnd*=3
125simulationSettings = exu.SimulationSettings()
126#simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
127simulationSettings.solution.file.write = False
128if not testIsActive:
129    simulationSettings.solution.file.writePeriod = 0.005
130    simulationSettings.solution.file.write = True
131    simulationSettings.solution.file.name = 'solution/coordinatesSolution.txt'
132else:
133    simulationSettings.solution.file.export.accelerations = False
134    simulationSettings.solution.file.export.velocities = False
135    
136simulationSettings.solution.sensors.writePeriod = h*10
137simulationSettings.solution.precision = 8 #make files smaller
138# simulationSettings.show.computationTime = True
139# simulationSettings.show.globalTimers = True
140#simulationSettings.show.statistics = True
141simulationSettings.timeIntegration.verboseMode = 1
142simulationSettings.parallel.numberOfThreads = 1
143
144simulationSettings.timeIntegration.newton.numericalDifferentiation.forODE2 = False
145simulationSettings.timeIntegration.newton.useModifiedNewton = False
146
147SC.visualizationSettings.general.graphicsUpdateInterval=0.05
148# SC.visualizationSettings.view0.scene.drawWorldBasis = True
149SC.visualizationSettings.general.circleTiling=200
150SC.visualizationSettings.view0.scene.drawCoordinateSystem=True
151SC.visualizationSettings.loads.show=False
152SC.visualizationSettings.bodies.show=True
153SC.visualizationSettings.markers.show=False
154
155SC.visualizationSettings.nodes.show=True
156SC.visualizationSettings.nodes.showBasis = True
157SC.visualizationSettings.nodes.basisSize = 1
158SC.visualizationSettings.nodes.drawNodesAsPoint = False
159SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
160SC.visualizationSettings.nodes.tiling = 4
161SC.visualizationSettings.openGL.drawFaceNormals = False
162
163SC.visualizationSettings.openGL.multiSampling = 4
164SC.visualizationSettings.openGL.light0.shadow = 0.25
165SC.visualizationSettings.openGL.light0.position = [-3,3,10,0]
166
167if not testIsActive:
168    SC.visualizationSettings.general.autoFitScene = False
169    SC.renderer.Start()
170    SC.renderer.RestoreSavedState()
171    SC.renderer.DoIdleTasks()
172
173simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
174simulationSettings.timeIntegration.endTime = tEnd
175simulationSettings.timeIntegration.explicit.computeEndOfStepAccelerations = False #increase performance, accelerations less accurate
176mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
177# mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ODE23)
178
179#compute error:
180uSum=0
181for node in evalNodes:
182    u = mbs.GetNodeOutput(node, exu.OutputVariableType.Coordinates)
183    exu.Print('coords node'+str(node)+' =',u)
184    for c in u:
185        uSum += abs(c) #add up all coordinates for comparison
186
187
188exu.Print('solution of generalContactImplicit1=',uSum)
189
190exu.sys['testResult'] = uSum
191
192    
193if not testIsActive:
194    SC.renderer.DoIdleTasks()
195
196    if True:
197        SC.visualizationSettings.general.autoFitScene = False
198        SC.visualizationSettings.general.graphicsUpdateInterval=0.02
199        
200        mbs.SolutionViewer(timeout=0.01)
201
202    SC.renderer.Stop() #safely close rendering window!
203
204if not testIsActive:
205    mbs.PlotSensor([], closeAll=True)
206    mbs.PlotSensor([sNode0]*3, [0,1,2], figureName='sphere position')
207    mbs.PlotSensor([omegaNode0]*3, [0,1,2], figureName='sphere angular vel')