symbolicUserFunctionMasses.py
You can view and download this file on Github: symbolicUserFunctionMasses.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: car with wheels modeled by ObjectConnectorRollingDiscPenalty
5# formulation is still under development and needs more testing
6#
7# Author: Johannes Gerstmayr
8# Date: 2023-11-28
9#
10# 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.
11#
12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
13import sys
14sys.exudynFast = True
15
16import exudyn as exu
17esym = exu.symbolic
18from exudyn.utilities import *
19import numpy as np
20#from math import pi, sin, cos
21
22SC = exu.SystemContainer()
23mbs = SC.AddSystem()
24
25
26#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
27useSymbolicUF = True
28if useSymbolicUF:
29 mysym = esym #np or esym
30 myabs = esym.abs
31 myReal = esym.Real
32else: #regular mode with numpy / Python functions
33 mysym = np #np or esym
34 myabs = abs
35 myReal = float
36
37def springForceUserFunction(mbs, t, itemNumber,
38 deltaL, deltaL_t, stiffness,
39 damping, force):
40 #f0 = damping*deltaL_t + stiffness*deltaL + force
41 f0 = (damping*deltaL_t + stiffness*deltaL)*(mysym.sign(-deltaL+myReal(0.005))+1.)*0.5
42 #f0 = (damping*deltaL_t + stiffness*deltaL)*esym.IfThenElse(deltaL > 0.01,myReal(0),myReal(1))
43 return f0
44
45#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
46
47
48Lx = 2
49Ly= 0.2
50Lz = 0.2
51aFact = 2
52
53ff = 1 #use 2 for finer model
54nx = 25*ff
55ny = 4*ff
56# nx = 15
57# ny = 3
58nz = 1*ny
59
60g = [0,-9.81,0]
61m = 1 / (nx*ny*4)
62stiffness = 200
63damping = 0.5e-3*stiffness
64
65drawSize = 0.2/ff
66showSprings = False
67
68bList = np.zeros((nx,ny,nz),dtype=exu.ObjectIndex)
69cList = []
70for i in range(nx):
71 for j in range(ny):
72 for k in range(nz):
73 x = i/nx
74
75 fact = 3*aFact**(-x*6)+1
76 aLy = Ly*fact
77 aLz = Lz #*fact
78 p = np.array([Lx*x, -0.5*aLy+aLy*j/(ny-1), -0.5*aLz+aLz*k/(nz-1)])
79 if i == 0:
80 b = mbs.AddObject(ObjectGround(referencePosition=p))
81 else:
82 b = mbs.CreateMassPoint(referencePosition=p, physicsMass=m, gravity=g, drawSize = drawSize)
83
84 # print('b=',b,': i',i,': j',j,': k',k)
85 bList[i,j,k] = b
86
87
88for i in range(nx):
89 for j in range(ny):
90 for k in range(nz):
91 if i > 0:
92 rr = 1. - np.random.rand()*0.05
93 cList += [mbs.CreateSpringDamper( bodyList=[bList[i,j,k], bList[i-1,j,k]], stiffness=stiffness*rr, damping=damping,drawSize = 0, show=showSprings)]
94 if j>0:
95 cList += [mbs.CreateSpringDamper(bodyList=[bList[i,j,k], bList[i-1,j-1,k]], stiffness=stiffness*rr, damping=damping,drawSize = 0, show=showSprings)]
96 cList += [mbs.CreateSpringDamper(bodyList=[bList[i,j,k], bList[i ,j-1,k]], stiffness=stiffness*rr, damping=damping,drawSize = 0, show=showSprings)]
97 if k>0:
98 cList += [mbs.CreateSpringDamper(bodyList=[bList[i,j,k], bList[i,j-1,k-1]], stiffness=stiffness*rr, damping=damping,drawSize = 0, show=showSprings)]
99 if j<ny-1:
100 cList += [mbs.CreateSpringDamper(bodyList=[bList[i,j,k], bList[i-1,j+1,k]], stiffness=stiffness*rr, damping=damping,drawSize = 0, show=showSprings)]
101 if k>0:
102 cList += [mbs.CreateSpringDamper(bodyList=[bList[i,j,k], bList[i-1,j,k-1]], stiffness=stiffness*rr, damping=damping,drawSize = 0, show=showSprings)]
103 cList += [mbs.CreateSpringDamper(bodyList=[bList[i,j,k], bList[i ,j,k-1]], stiffness=stiffness*rr, damping=damping,drawSize = 0, show=showSprings)]
104 if k<nz-1:
105 cList += [mbs.CreateSpringDamper(bodyList=[bList[i,j,k], bList[i-1,j,k+1]], stiffness=stiffness*rr, damping=damping,drawSize = 0, show=showSprings)]
106 if j>0:
107 cList += [mbs.CreateSpringDamper(bodyList=[bList[i,j,k], bList[i,j-1,k+1]], stiffness=stiffness*rr, damping=damping,drawSize = 0, show=showSprings)]
108
109
110# CAUTION: this only works, if every object has its own user function!!!
111if useSymbolicUF: #do assemble before adding user functions => then, they will be parallelized
112 mbs.Assemble()
113
114listUF = []
115if True:
116 isNew=1
117 for cc in cList:
118 if useSymbolicUF:
119 #create separate user function for each spring-damper (multithreading)!
120 symbolicFunc = CreateSymbolicUserFunction(mbs, springForceUserFunction, 'springForceUserFunction', cc)
121 mbs.SetObjectParameter(cc, 'springForceUserFunction', symbolicFunc)
122 listUF += [symbolicFunc] #store, such that they are not deleted!!!
123 else:
124 mbs.SetObjectParameter(cc, 'springForceUserFunction', springForceUserFunction)
125
126
127#assemble and solve system for default parameters
128if not useSymbolicUF: #do assemble before adding user functions => then, they will be parallelized
129 mbs.Assemble()
130
131endTime = 10
132stepSize = 1e-4
133if ff==1: stepSize = 5e-4
134if ff>2: stepSize = 0.5e-5
135
136# endTime = 0.01
137simulationSettings = exu.SimulationSettings()
138simulationSettings.solutionSettings.solutionWritePeriod = 0.01
139# simulationSettings.solutionSettings.writeSolutionToFile = False
140
141simulationSettings.timeIntegration.numberOfSteps = int(endTime/stepSize)
142simulationSettings.timeIntegration.endTime = endTime
143simulationSettings.timeIntegration.newton.useModifiedNewton = True
144simulationSettings.parallel.numberOfThreads = 8
145simulationSettings.timeIntegration.verboseMode = 1
146simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
147#simulationSettings.displayComputationTime = True
148
149SC.visualizationSettings.nodes.drawNodesAsPoint = False
150SC.visualizationSettings.nodes.tiling = 6
151SC.visualizationSettings.loads.show = False
152
153SC.visualizationSettings.contour.outputVariableComponent = -1
154SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
155SC.visualizationSettings.window.renderWindowSize = [2000,1600]
156SC.visualizationSettings.openGL.multiSampling = 4
157
158import time
159exu.StartRenderer()
160ts = time.time()
161print('start simulation')
162mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.RK44)
163print('finished: ', time.time()-ts, 'seconds')
164mbs.WaitForUserToContinue()
165exu.StopRenderer()
166# mbs.SolutionViewer()
167#i7-1390, boost
168#results for ExplicitMidpoint, 100 steps, 6272 nodes, RK44, exudynFast=True:
169#8 threads:
170#C++: 2.62s
171#Symbolic user function:4.74s
172#1 thread:
173#C++: 9.83s
174#Symbolic user function:11.04s
175#Python user function: 81.8s => speedup = 38
176
177#i9:
178#results for ExplicitMidpoint, 100 steps, 6272 nodes, RK44, exudynFast=True:
179#8 threads:
180#C++: 3.27s
181#Symbolic user function:6.16s
182#Python user function: 62.72s => speedup = 20.57