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