SpringDamperMassUserFunction.py
You can view and download this file on Github: SpringDamperMassUserFunction.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: The 2D movement of a point mass system is simulated.
5# As compared to a similar example, here it uses the itemInterface.py and
6# it uses user functions for springs and dampers
7#
8# Author: Johannes Gerstmayr
9# Date: 2019-12-04
10# Update: 2023-12-08 (symbolic user function)
11#
12# 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.
13#
14#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15
16import exudyn as exu
17from exudyn.itemInterface import *
18from exudyn.utilities import *
19
20import numpy as np
21
22SC = exu.SystemContainer()
23mbs = SC.AddSystem()
24
25useSymbolicUserFunction = True
26
27#defines relative displacement, relative velocity, stiffness k, damping d, and additional spring force f0
28def springForce(mbs, t, itemIndex, u, v, k, d, f0):
29 return u*k+v*d
30
31sqrt2 = 2**0.5
32nBodies = 24 #24; 240*4 #480 for Eigen factorization test
33nBodies2 = 3 #6; 30*4 #60 for Eigen factorization test
34
35coList = []
36
37for j in range(nBodies2):
38# body = mbs.AddObject({'objectType': 'Ground', 'referencePosition': [0,j,0]})
39# mbs.AddMarker({'markerType': 'BodyPosition', 'bodyNumber': body, 'localPosition': [0.0, 0.0, 0.0], 'bodyFixed': False})
40 body = mbs.AddObject(ObjectGround(referencePosition=[0,j,0]))
41 mbs.AddMarker(MarkerBodyPosition(bodyNumber = body, localPosition = [0.0, 0.0, 0.0]))
42
43 for i in range(nBodies-1):
44 #2D:
45 node = mbs.AddNode(NodePoint2D(referenceCoordinates=[i+1, j], initialCoordinates=[0, 0]))
46 body = mbs.AddObject(MassPoint2D(physicsMass=10, nodeNumber=node))
47 mBody = mbs.AddMarker(MarkerBodyPosition(bodyNumber=body, localPosition=[0,0,0]))
48 #dynamic/explicit:
49 #mbs.AddLoad(LoadForceVector(markerNumber = mBody, loadVector = [0, -0.025*100, 0]))
50 #static:
51 mbs.AddLoad(LoadForceVector(markerNumber = mBody, loadVector = [0, -0.025, 0]))
52
53#add spring-dampers:
54for j in range(nBodies2-1):
55 for i in range(nBodies-1):
56 coList += [mbs.AddObject(ObjectConnectorSpringDamper(markerNumbers=[j*nBodies + i,j*nBodies + i+1], stiffness=4000,
57 damping=10, force=0, referenceLength=1, springForceUserFunction = springForce))]
58 coList += [mbs.AddObject(ObjectConnectorSpringDamper(markerNumbers=[j*nBodies + i,(j+1)*nBodies + i], stiffness=4000,
59 damping=10, force=0, referenceLength=1, springForceUserFunction = springForce))]
60 coList += [mbs.AddObject(ObjectConnectorSpringDamper(markerNumbers=[j*nBodies + i,(j+1)*nBodies + i+1], stiffness=4000,
61 damping=10, force=0, referenceLength=sqrt2, springForceUserFunction = springForce))] #diagonal elements
62
63for i in range(nBodies-1):
64 j = nBodies2-1
65 coList += [mbs.AddObject(ObjectConnectorSpringDamper(markerNumbers=[j*nBodies + i,j*nBodies + i+1], stiffness=4000,
66 damping=10, force=0, referenceLength=1, springForceUserFunction = springForce))]
67for j in range(nBodies2-1):
68 i = nBodies-1
69 coList += [mbs.AddObject(ObjectConnectorSpringDamper(markerNumbers=[j*nBodies + i,(j+1)*nBodies + i], stiffness=4000,
70 damping=10, force=0, referenceLength=1, springForceUserFunction = springForce))]
71
72#now set symbolic user functions
73if useSymbolicUserFunction:
74 #create symbolic version of Python user function (only works for limited kinds of functions)
75 #we have to keep this handle, do not overwrite:
76 symbolicFunc = CreateSymbolicUserFunction(mbs, springForce, 'springForceUserFunction', coList[0])
77 for co in coList:
78 #now inject symbolic user function into object, directly done inside C++ (no Pybind overhead):
79 #symbolicFunc.TransferUserFunction2Item(mbs, co, 'springForceUserFunction')
80 mbs.SetObjectParameter(co, 'springForceUserFunction', symbolicFunc)
81
82
83#optional:
84nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[-0.5,0,0])) #ground node for coordinate constraint
85mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
86mNC1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = 1, coordinate=1))
87##add constraint for testing (does not work in explicit computation):
88#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mNC1]))
89
90mbs.Assemble()
91print(mbs)
92
93useGraphics = True
94if useGraphics:
95 exu.StartRenderer()
96
97simulationSettings = exu.SimulationSettings()
98simulationSettings.timeIntegration.numberOfSteps = 100*200
99simulationSettings.timeIntegration.endTime = 1*200
100simulationSettings.solutionSettings.writeSolutionToFile = False
101simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
102simulationSettings.timeIntegration.verboseMode = 1
103simulationSettings.displayStatistics = True
104simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
105#simulationSettings.linearSolverType = exu.LinearSolverType.EXUdense
106simulationSettings.displayComputationTime = True
107
108SC.visualizationSettings.nodes.show = True
109SC.visualizationSettings.bodies.show = False
110SC.visualizationSettings.loads.show = False
111SC.visualizationSettings.markers.show = False
112SC.visualizationSettings.nodes.defaultSize = 0.2*2
113
114SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
115SC.visualizationSettings.contour.outputVariableComponent = 0 #y-component
116
117mbs.SolveDynamic(simulationSettings, solverType = exudyn.DynamicSolverType.ExplicitMidpoint)
118#u = mbs.GetNodeOutput(nBodies-2, exu.OutputVariableType.Position) #tip node
119#print('dynamic tip displacement (y)=', u[1]) #dense: -11.085967426937412, sparse:-11.085967426937431
120
121simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-7
122simulationSettings.staticSolver.newton.relativeTolerance = 1e-6*1e5 # make this large for linear computation
123simulationSettings.staticSolver.newton.absoluteTolerance = 1e-1
124simulationSettings.staticSolver.verboseMode = 1
125
126mbs.SolveStatic(simulationSettings)
127
128u = mbs.GetNodeOutput(nBodies-2, exu.OutputVariableType.Position) #tip node
129print('static tip displacement (y)=', u[1])
130staticError = u[1]-(-0.44056224799446486)
131
132if useGraphics:
133 SC.WaitForRenderEngineStopFlag()
134 exu.StopRenderer()