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()