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