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