sparseMatrixSpringDamperTest.py

You can view and download this file on Github: sparseMatrixSpringDamperTest.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Spring damper model for test of sparse matrix and sparse solver
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2019-11-20
  8#
  9# 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.
 10#
 11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 12
 13import exudyn as exu           #c++ bibliothek, liest Dictionaries
 14from exudyn.itemInterface import *     # conversion of data to exudyn dictionaries C interface
 15
 16useGraphics = True #without test
 17#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 18#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 19try: #only if called from test suite
 20    from modelUnitTests import exudynTestGlobals #for globally storing test results
 21    useGraphics = exudynTestGlobals.useGraphics
 22except:
 23    class ExudynTestGlobals:
 24        pass
 25    exudynTestGlobals = ExudynTestGlobals()
 26#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 27
 28SC = exu.SystemContainer()
 29mbs = SC.AddSystem()
 30
 31
 32sqrt2 = 2**0.5
 33nBodies = 72
 34nBodies2 = 6
 35
 36for j in range(nBodies2):
 37    body = mbs.AddObject({'objectType': 'Ground', 'referencePosition': [0,j,0]})
 38    mbs.AddMarker({'markerType': 'BodyPosition',  'bodyNumber': body,  'localPosition': [0.0, 0.0, 0.0], 'bodyFixed': False})
 39    for i in range(nBodies-1):
 40        #2D:
 41        node = mbs.AddNode(NodePoint2D(referenceCoordinates=[i+1, j], initialCoordinates=[0, 0]))
 42        body = mbs.AddObject(MassPoint2D(physicsMass=10, nodeNumber=node))
 43        mBody = mbs.AddMarker(MarkerBodyPosition(bodyNumber=body, localPosition=[0,0,0]))
 44        mbs.AddLoad({'loadType': 'ForceVector',  'markerNumber': mBody,  'loadVector': [0, -0.025, 0]})
 45
 46#add spring-dampers:
 47for j in range(nBodies2-1):
 48    for i in range(nBodies-1):
 49        mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': 4000, 'damping': 10, 'force': 0,
 50                        'referenceLength':1, 'markerNumbers': [j*nBodies + i,j*nBodies + i+1]})
 51        mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': 4000, 'damping': 10, 'force': 0,
 52                        'referenceLength':1, 'markerNumbers': [j*nBodies + i,(j+1)*nBodies + i]})
 53        #diagonal spring: l*sqrt(2)
 54        mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': 4000, 'damping': 10, 'force': 0,
 55                        'referenceLength':sqrt2, 'markerNumbers': [j*nBodies + i,(j+1)*nBodies + i+1]})
 56
 57for i in range(nBodies-1):
 58    j = nBodies2-1
 59    mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': 4000, 'damping': 10, 'force': 0,
 60                    'referenceLength':1, 'markerNumbers': [j*nBodies + i,j*nBodies + i+1]})
 61for j in range(nBodies2-1):
 62    i = nBodies-1
 63    mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': 4000, 'damping': 10, 'force': 0,
 64                    'referenceLength':1, 'markerNumbers': [j*nBodies + i,(j+1)*nBodies + i]})
 65
 66#add constraints for testing:
 67nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[-0.5,0,0])) #ground node for coordinate constraint
 68mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 69
 70mNC1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = 1, coordinate=1))
 71mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mNC1]))
 72
 73mbs.Assemble()
 74#exu.Print(mbs)
 75
 76#useGraphics = True
 77if useGraphics:
 78    exu.StartRenderer()
 79
 80simulationSettings = exu.SimulationSettings()
 81#simulationSettings.displayStatistics = True
 82
 83SC.visualizationSettings.nodes.show = False
 84SC.visualizationSettings.loads.show = False
 85SC.visualizationSettings.markers.show = False
 86#SC.visualizationSettings.nodes.defaultSize = 0.05
 87
 88
 89simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-5*0.01
 90#simulationSettings.staticSolver.newton.relativeTolerance = 1e-6#*1e5
 91#simulationSettings.staticSolver.newton.absoluteTolerance = 1e-1
 92#simulationSettings.staticSolver.numberOfLoadSteps = 10
 93#simulationSettings.staticSolver.loadStepGeometric = True
 94simulationSettings.staticSolver.verboseMode = 2
 95
 96#dense solver:
 97simulationSettings.linearSolverType = exu.LinearSolverType.EXUdense
 98mbs.SolveStatic(simulationSettings)
 99
100u = mbs.GetNodeOutput(nBodies-2, exu.OutputVariableType.Position) #tip node
101exu.Print('static tip displacement (y)=', u[1])
102exudynTestGlobals.testError = u[1]-(-6.779862983765133) #72 x 6 bodies; CPUtime surface: 0.55 seconds
103exudynTestGlobals.testResult = u[1]
104
105#sparse solver:
106simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
107mbs.SolveStatic(simulationSettings)
108
109u = mbs.GetNodeOutput(nBodies-2, exu.OutputVariableType.Position) #tip node
110exu.Print('static tip displacement (y)=', u[1])
111
112#factor 1e-2: 32bit version shows 2.1e-12 error
113exudynTestGlobals.testError = 1e-2*(u[1]-(-6.779862983766792)) #72 x 6 bodies; CPUtime surface: 0.029 seconds
114exudynTestGlobals.testResult = 1e-2*u[1]
115
116if useGraphics:
117    SC.WaitForRenderEngineStopFlag()
118    exu.StopRenderer()