contactCurveExample.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Simple test for ObjectContactCurveCircles
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2025-05-10
  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#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 12import exudyn as exu
 13from exudyn.utilities import *
 14from math import sin, cos
 15import exudyn.graphics as graphics
 16import numpy as np
 17
 18testIsActive = exu.sys.get('testIsActive', False)
 19
 20SC = exu.SystemContainer()
 21mbs = SC.AddSystem()
 22
 23L = 2
 24a = 0.1 
 25b = 0.1
 26rPin = a*0.75
 27
 28oGround = mbs.CreateGround(graphicsDataList=[graphics.CheckerBoard(point=[0,0,-b], size=3*L)])
 29
 30contactStiffness = 2e5*10
 31contactDamping = 5e2*10
 32
 33inertiaBrick=InertiaCuboid(1000, sideLengths=[L,a,b])
 34
 35#Example 1: arm moving in given curve:
 36if True:
 37    pOff = np.array([-L,0,0])
 38    inertiaBrick1 = inertiaBrick.Translated([-0.5*L,0,0])
 39    lever = mbs.CreateRigidBody(
 40                                referencePosition=pOff+[1*L+1*rPin,-0.5*L+rPin,0],
 41                                inertia=inertiaBrick1,
 42                                gravity = [0,-g,0],
 43                                graphicsDataList=[graphics.Brick(centerPoint=[-0.5*L,0,0], size=[L,a,b], color=graphics.color.dodgerblue),
 44                                                  graphics.Basis(inertiaBrick1.COM()),
 45                                                  graphics.Sphere(point=[-L,0,0], radius=rPin, nTiles=16, color=graphics.color.darkgrey)],
 46                                #create2D=True, #not possible here, as COM must be 0 with 2D
 47                                )
 48    
 49    mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=lever, localPosition=[-L,0,0]))
 50    
 51    #create 2D points for contact curve (linear segments)    
 52    pList = []
 53    n = 64
 54    
 55    for i in range(2*n+2):
 56        phi = -1*pi+i/(n-1)*1*pi
 57        r = L*0.5
 58        if i >= n:
 59            phi = 0.*pi-phi
 60            r = 0.5*L-2*rPin
 61        x = r*sin(phi)
 62        y = -r*cos(phi)
 63        if i == n:
 64            x = 0.25*L
 65            y = -0.5*L
 66        if i == n+1:
 67            x = 0.25*L
 68            y = -0.5*L+2*rPin
 69            
 70        pList.append([x,y])
 71    
 72    
 73    #transform points into contact segments:
 74    nPoints = len(pList)
 75    segmentsData = np.zeros((nPoints,4))
 76    segmentsData[:,0:2] = pList
 77    segmentsData[:,2:4] = np.roll(pList,2) #roll is element wise on rows and columns, therefore 2>shift one row
 78    
 79    segmentsData = segmentsData[1:,:] #we do not like to close curve, so remove first segment
 80    nSeg = len(segmentsData)
 81            
 82    mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=pOff+[0,0,0]))
 83    nGenericData = mbs.AddNode(NodeGenericData(initialCoordinates=[-1,0,0]*nSeg,
 84                                               numberOfDataCoordinates=3*nSeg))
 85
 86    mbs.AddObject(ObjectContactCurveCircles(markerNumbers=[mGround]+[mBody], 
 87                                            nodeNumber=nGenericData,
 88                                            circlesRadii=[rPin], 
 89                                            segmentsData=exu.MatrixContainer(segmentsData), 
 90                                            contactStiffness=contactStiffness, contactDamping=contactDamping,
 91                                            visualization=VObjectContactCurveCircles(show=True, color=graphics.color.blue)
 92                                            ))
 93
 94    mbs.CreateCoordinateConstraint(itemNumbers=[lever,None],coordinates=[1,None])
 95    mbs.CreateCoordinateConstraint(itemNumbers=[lever,None],coordinates=[2,None])
 96
 97    def UFoffset_t(mbs, t, itemNumber, lOffset): #time derivative of UFoffset
 98        return 2.2*L*sin(2*pi*t)
 99    
100    def UFoffset(mbs, t, itemNumber, lOffset): #time derivative of UFoffset
101        return -0.7*(cos(2*pi*t)-1)
102
103    mbs.CreateCoordinateConstraint(itemNumbers=[lever,None],coordinates=[0,None],
104                                   #velocityLevel=True,
105                                   offsetUserFunction_t=UFoffset_t,
106                                   offsetUserFunction=UFoffset,
107                                   )
108
109#Example 2: pendulum in bucket:
110if True:
111    pOff = np.array([L*1.5,0,0])
112
113    lever = mbs.CreateRigidBody(referencePosition=pOff+[-0.9*L,rPin,0],
114                                inertia=inertiaBrick,
115                                gravity = [0,-g,0],
116                                create2D=True,
117                                graphicsDataList=[graphics.Brick(size=[L,a,b], color=graphics.color.orange),
118                                                  graphics.Sphere(point=[0.5*L,0,0], radius=rPin, nTiles=16, color=graphics.color.darkgrey)],
119                                )
120    
121    mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=lever, localPosition=[0.5*L,0,0]))
122
123    #create 2D points for contact curve (linear segments)    
124    pList = []
125    
126    pList.append([-1,0.5])
127    pList.append([-1,0.])
128    pList.append([1,0.])
129    pList.append([1,0.5])
130
131    #transform points into contact segments:
132    nPoints = len(pList)
133    segmentsData = np.zeros((nPoints,4))
134    segmentsData[:,0:2] = pList
135    segmentsData[:,2:4] = np.roll(pList,2) #roll is element wise on rows and columns, therefore 2>shift one row
136    
137    segmentsData = segmentsData[1:,:]
138    nSeg = len(segmentsData)
139            
140    mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=pOff+[0,0,0]))
141    nGenericData = mbs.AddNode(NodeGenericData(initialCoordinates=[-1,0,0]*nSeg,
142                                               numberOfDataCoordinates=3*nSeg))
143
144    mbs.AddObject(ObjectContactCurveCircles(markerNumbers=[mGround]+[mBody], 
145                                            nodeNumber=nGenericData,
146                                            circlesRadii=[rPin], 
147                                            segmentsData=exu.MatrixContainer(segmentsData), 
148                                            contactStiffness=contactStiffness, contactDamping=contactDamping,
149                                            visualization=VObjectContactCurveCircles(show=True, color=graphics.color.blue)
150                                            ))
151
152
153mbs.Assemble()
154
155stepSize=0.001
156tEnd = 1
157simulationSettings = exu.SimulationSettings()
158simulationSettings.solution.file.write = (not testIsActive)
159simulationSettings.solution.file.writePeriod = 0.01
160simulationSettings.solution.sensors.writePeriod = stepSize  #output interval
161simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
162simulationSettings.timeIntegration.endTime = tEnd
163# simulationSettings.timeIntegration.realtime.active = True
164#simulationSettings.timeIntegration.realtime.factor = 0.5
165# simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-2
166# simulationSettings.timeIntegration.discontinuous.useRecommendedStepSize = False
167
168#simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
169simulationSettings.timeIntegration.newton.useModifiedNewton = True
170#simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
171
172simulationSettings.timeIntegration.verboseMode = 1
173
174SC.visualizationSettings.general.graphicsUpdateInterval = 0.02
175SC.visualizationSettings.view0.window.renderWindowSize=[1600,2000]
176SC.visualizationSettings.openGL.multiSampling=4
177#SC.visualizationSettings.view0.scene.facesTransparent=True
178SC.visualizationSettings.openGL.light0.shadow=0.3*(not testIsActive)
179SC.visualizationSettings.loads.show = False
180SC.visualizationSettings.connectors.showContact = True
181
182
183mbs.SolveDynamic(simulationSettings)
184
185#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
186ode2 = mbs.systemData.GetODE2Coordinates()
187
188u = 0.1*np.linalg.norm(ode2)
189exu.Print('solution of contactCurveExample=',u) 
190
191exu.sys['testResult'] = u
192#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
193
194if not testIsActive:
195    mbs.SolutionViewer()