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