createFunctionsTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test model for MainSystem Create functions; the model itself is not very meaningful, but it uses many Create functions
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2025-05-11
  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
 17import numpy as np
 18
 19testIsActive = exu.sys.get('testIsActive', False)
 20
 21
 22#set up new multibody system to work with
 23SC = exu.SystemContainer()
 24mbs = SC.AddSystem()
 25
 26graphicsGround = graphics.CheckerBoard(point=[0,2,-1], size=8)
 27oGround = mbs.CreateGround(referencePosition=[0,0,0], graphicsDataList=[graphicsGround])
 28
 29oMass = mbs.CreateMassPoint(mass=5, referencePosition=[1,0,0], 
 30                            initialDisplacement=[0,0,0],
 31                            initialVelocity=[0,0.5,0],
 32                            gravity=[0,-9.81,0],
 33                            drawSize=0.3, color=graphics.color.orange,
 34                            )
 35
 36oSpringDamper = mbs.CreateSpringDamper(itemNumbers=[oGround, oMass], #[body0,body1]
 37                                      localPosition0=[0,0,0], #locally on body0
 38                                      localPosition1=[0,0,0], #locally on body1
 39                                      referenceLength=None, #usually set to None (default) => takes the (relaxed) length in reference configuration
 40                                      stiffness=1e2, damping=1)
 41
 42oDistance = mbs.CreateDistanceConstraint(itemNumbers=[oGround, oMass], #[body0,body1]
 43                                         localPosition0 = [0,0,0], #locally on body0
 44                                         localPosition1 = [0,0,0], #locally on body1
 45                                         distance=None)
 46
 47#%%++++++++++++++++++++++++++++++++++++++++++++++++++
 48
 49mass = 10 
 50outerRadius = 0.45
 51length = 2
 52volume = np.pi * outerRadius**2 * length
 53inertiaCylinder = InertiaCylinder(density=mass/volume, length=length, outerRadius=outerRadius, 
 54                                  axis=0, innerRadius=0.7*outerRadius)
 55
 56oBody = mbs.CreateRigidBody(inertia=inertiaCylinder,
 57                            referencePosition=[3,0,0],
 58                            initialVelocity=[0,0,0],
 59                            initialAngularVelocity=[2*np.pi,0,0],
 60                            gravity=[0,-9.81,0],
 61                            #graphicsDataList=[graphicsCyl]
 62                            )
 63
 64
 65rbX = 2
 66rbY = 0.5
 67rbZ = 0.2
 68
 69#create inertia instance for cuboid (brick) shape
 70mass=10 #kg
 71volume=rbX*rbY*rbZ
 72inertiaCube2 = InertiaCuboid(density=mass/volume, sideLengths=[rbX,rbY,rbZ])
 73inertiaCube2 = inertiaCube2.Translated([0.5,0,0])
 74
 75oBody2 = mbs.CreateRigidBody(inertia = inertiaCube2,
 76                             referencePosition = [4,-0.5*rbX,0], 
 77                             referenceRotationMatrix = RotationVector2RotationMatrix([0,0,-0.5*np.pi]),
 78                             initialAngularVelocity = [0,0,0.2*2*np.pi],
 79                             initialVelocity = [0.2*np.pi*0.5*rbX,0,0],
 80                             gravity = [0,-9.81,0],
 81                             color=graphics.color.blue,
 82                             #graphicsDataList=[graphicsCube, graphics.Basis(origin=inertiaCube2.COM())],
 83                             )
 84
 85loadMassPoint = mbs.CreateForce(itemNumber=oMass, loadVector=[10,0,0])
 86loadRigidBody = mbs.CreateForce(itemNumber=oBody, localPosition=[0,0,0], loadVector=[10,0,0])
 87
 88def UFforce(mbs, t, loadVector):
 89    return (10+5*np.sin(t*10*2*np.pi))*np.array([0,5,0])
 90
 91mbs.CreateForce(itemNumber=oBody,
 92                localPosition=[0,1.2,0.5], #position at body
 93                loadVectorUserFunction=UFforce)
 94
 95def UserFunctionTorque(mbs, t, loadVector):
 96    return np.cos(t*2*np.pi)*np.array(loadVector)
 97
 98#add torque with user function
 99mbs.CreateTorque(itemNumber=oBody, loadVector=[5,0,0], 
100                 loadVectorUserFunction=UserFunctionTorque)
101
102#+++++++++
103inertiaSphere = InertiaSphere(mass=2, radius=0.25)
104gSphere = inertiaSphere.GetGraphics(graphics.color.red)
105
106oMass1 = mbs.CreateMassPoint(mass=2, referencePosition=[1,0,0],
107                             gravity=[0,-9.81,0],
108                             # drawSize=0.25, color=graphics.color.red,
109                             graphicsDataList=[gSphere]
110                             )
111
112#create spherical joint between ground and mass point; could also be applied to two bodies; possible bodies: mass point or rigid body
113mbs.CreateSphericalJoint(itemNumbers=[oGround, oMass1],
114                         position=[1,0,0], #global position of joint (in reference configuration)
115                         constrainedAxes=[1,0,0]) #x,y,z directions: 1 = fixed, 0 = free motion
116
117#constrain Y and Z coordinate of mass point to move z=10*y
118mbs.CreateCoordinateConstraint(itemNumbers=[oMass1, oMass1], 
119                               coordinates=[2,1],
120                               factor1=10)
121
122mbs.CreatePrismaticJoint(itemNumbers=[oGround, oBody], 
123                         position=[3,0,0], #global position of joint
124                         axis=[1,0,0], #global axis of joint, can move in global x-direction
125                         useGlobalFrame=True) #use local coordinates for joint definition
126
127mbs.CreateRevoluteJoint(itemNumbers=[oBody, oBody2], 
128                        position=[4,-0.5,0], #global position of joint
129                        axis=[0,0,1], #rotation along global z-axis
130                        useGlobalFrame=True)
131
132mbs.CreateTorsionalSpringDamper(itemNumbers=[oBody, oBody2],
133                                position=[4,-0.5,0], #global position of spring-damper
134                                axis=[0,0,1],        #global rotation axis
135                                stiffness=1000,
136                                damping=20,
137                                useGlobalFrame=True)
138
139mass = 5 
140mu = 0.05 #static friction
141rDisc = 0.5
142length = 0.1
143volume = np.pi * rDisc**2 * length
144inertiaCylinder = InertiaCylinder(density=mass/volume, length=length, outerRadius=rDisc, axis=0)
145
146oDisc = mbs.CreateRigidBody(inertia=inertiaCylinder,
147                            nodeType=exu.NodeType.RotationRxyz, #this allows to constrain Z-rotation
148                            referencePosition=[0,2,rDisc-1], #reference position x/y/z of COM
149                            initialVelocity=[0,2*np.pi*rDisc,0],   
150                            initialAngularVelocity=[-2*np.pi,0.2,0],
151                            gravity=[0,0,-9.81])
152
153mbs.CreateRollingDiscPenalty(bodyNumbers=[oGround, oDisc], discRadius = rDisc, 
154                             axisPosition=[0,0,0], axisVector=[1,0,0],
155                             planePosition = [0,0,-1], planeNormal = [0,0,1],
156                             contactStiffness = 1e4, contactDamping = 2e2, 
157                             dryFriction = [mu,mu], color=graphics.color.red)
158
159#constrain Z and Y-rotation of disc, so it rotates only around X
160mbs.CreateCoordinateConstraint(itemNumbers=[oDisc, None], coordinates=[5,None])
161mbs.CreateCoordinateConstraint(itemNumbers=[oDisc, oGround], coordinates=[4,None])#also works with ground
162
163mass = 5 
164rDisc = 0.5
165length = 0.1
166volume = np.pi * rDisc**2 * length
167inertiaCylinder = InertiaCylinder(density=mass/volume, length=length, outerRadius=rDisc, axis=0)
168
169#create a free rigid body using the defined inertia properties and applies gravity (y-direction).
170oDisc1 = mbs.CreateRigidBody(inertia=inertiaCylinder,
171                            referencePosition=[1,2,rDisc-1], #reference position x/y/z of COM
172                            initialVelocity=[0,2*np.pi*rDisc,0],   
173                            initialAngularVelocity=[-2*np.pi,0.2,0],
174                            gravity=[0,0,-9.81])
175
176mbs.CreateForce(itemNumber=oDisc1, loadVector=[10,0,0])
177
178#create a 'rolling' joint between flat ground defined by plane, lying on oGround, and rigid body given as oDisc:
179mbs.CreateRollingDisc(bodyNumbers=[oGround, oDisc1], discRadius = rDisc, 
180                      axisPosition=[0,0,0], axisVector=[1,0,0], #relative to oDisc
181                      planePosition = [0,0,-1], planeNormal = [0,0,1], #defines plane
182                      constrainedAxes=[1,1,1] #constrain 3 axes: lateral motion, forward motion and normal contact
183                      )
184
185#constrain the two discs to have same X-motion
186mbs.CreateCoordinateConstraint(itemNumbers=[oDisc,oDisc1],
187                               velocityLevel=True,
188                               coordinates=[0,0])
189
190mbs.CreateCartesianSpringDamper(
191    itemNumbers=[oGround, oBody2], #[body0,body1]
192    localPosition0=[4,-0.5*rbX,0], #for body0
193    localPosition1=[0,0,0],        #for body1
194    stiffness = [100,10,10], #x,y,z stiffness
195    damping = [5,2,2], #x,y,z damping
196)
197
198#%%++++++++++++++++++++++++++++++++++++++++++++++
199mbs.Assemble()
200
201endTime = 0.5 #simulation time in seconds
202if not testIsActive:
203    endTime = 2.5
204
205stepSize = 0.002 #should be small enough to achieve sufficient accuracy
206
207#some simulation parameters:
208simulationSettings = exu.SimulationSettings() #takes currently set values or default values
209simulationSettings.timeIntegration.numberOfSteps = int(endTime/stepSize)
210simulationSettings.timeIntegration.endTime = endTime
211simulationSettings.timeIntegration.verboseMode = 1
212simulationSettings.timeIntegration.newton.useModifiedNewton = True
213
214#for redundant constraints, use these two settings
215simulationSettings.linearSolver.ignoreSingularJacobian=True
216simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenDense #use EigenSparse for larger systems alternatively
217#simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
218# simulationSettings.timeIntegration.newton.useModifiedNewton = True
219SC.visualizationSettings.openGL.multiSampling = 4
220
221mbs.SolveDynamic(simulationSettings)
222
223
224#%%
225#TAG: visualization
226#DIFFICULTY: 10000
227SC.visualizationSettings.nodes.drawNodesAsPoint = False
228SC.visualizationSettings.nodes.showBasis = True
229
230if not testIsActive:
231    mbs.SolutionViewer()
232
233#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
234ode2 = mbs.systemData.GetODE2Coordinates()
235
236u = 0.01*np.linalg.norm(ode2)
237exu.Print('solution of createFunctionsTest=',u) 
238
239exu.sys['testResult'] = u
240#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++