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
 16testIsActive = exu.sys.get('testIsActive', False)
 17
 18SC = exu.SystemContainer()
 19mbs = SC.AddSystem()
 20
 21
 22sqrt2 = 2**0.5
 23nBodies = 72
 24nBodies2 = 6
 25
 26for j in range(nBodies2): 
 27    body = mbs.AddObject({'objectType': 'Ground', 'referencePosition': [0,j,0]})
 28    mbs.AddMarker({'markerType': 'BodyPosition',  'bodyNumber': body,  'localPosition': [0.0, 0.0, 0.0], 'bodyFixed': False})
 29    for i in range(nBodies-1): 
 30        #2D:
 31        node = mbs.AddNode(NodePoint2D(referenceCoordinates=[i+1, j], initialCoordinates=[0, 0]))
 32        body = mbs.AddObject(MassPoint2D(mass=10, nodeNumber=node))
 33        mBody = mbs.AddMarker(MarkerBodyPosition(bodyNumber=body, localPosition=[0,0,0]))
 34        mbs.AddLoad({'loadType': 'ForceVector',  'markerNumber': mBody,  'loadVector': [0, -0.025, 0]})
 35
 36#add spring-dampers:
 37for j in range(nBodies2-1): 
 38    for i in range(nBodies-1): 
 39        mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': 4000, 'damping': 10, 'force': 0,
 40                        'referenceLength':1, 'markerNumbers': [j*nBodies + i,j*nBodies + i+1]})
 41        mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': 4000, 'damping': 10, 'force': 0,
 42                        'referenceLength':1, 'markerNumbers': [j*nBodies + i,(j+1)*nBodies + i]})
 43        #diagonal spring: l*sqrt(2)
 44        mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': 4000, 'damping': 10, 'force': 0,
 45                        'referenceLength':sqrt2, 'markerNumbers': [j*nBodies + i,(j+1)*nBodies + i+1]})
 46
 47for i in range(nBodies-1): 
 48    j = nBodies2-1
 49    mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': 4000, 'damping': 10, 'force': 0,
 50                    'referenceLength':1, 'markerNumbers': [j*nBodies + i,j*nBodies + i+1]})
 51for j in range(nBodies2-1): 
 52    i = nBodies-1
 53    mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': 4000, 'damping': 10, 'force': 0,
 54                    'referenceLength':1, 'markerNumbers': [j*nBodies + i,(j+1)*nBodies + i]})
 55
 56#add constraints for testing:
 57nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[-0.5,0,0])) #ground node for coordinate constraint
 58mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 59
 60mNC1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = 1, coordinate=1))
 61mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mNC1]))
 62
 63mbs.Assemble()
 64#exu.Print(mbs)
 65
 66#testIsActive = False
 67if not testIsActive: 
 68    SC.renderer.Start()
 69
 70simulationSettings = exu.SimulationSettings()
 71#simulationSettings.show.statistics = True
 72
 73SC.visualizationSettings.nodes.show = False
 74SC.visualizationSettings.loads.show = False
 75SC.visualizationSettings.markers.show = False
 76#SC.visualizationSettings.nodes.defaultSize = 0.05
 77
 78
 79simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-5*0.01
 80#simulationSettings.staticSolver.newton.relativeTolerance = 1e-6#*1e5
 81#simulationSettings.staticSolver.newton.absoluteTolerance = 1e-1
 82#simulationSettings.staticSolver.numberOfLoadSteps = 10
 83#simulationSettings.staticSolver.loadStepGeometric = True
 84simulationSettings.staticSolver.verboseMode = 2
 85
 86#dense solver:
 87simulationSettings.linearSolver.solverType = exu.LinearSolverType.EXUdense
 88mbs.SolveStatic(simulationSettings)
 89
 90u = mbs.GetNodeOutput(nBodies-2, exu.OutputVariableType.Position) #tip node
 91exu.Print('static tip displacement (y)=', u[1])
 92exu.sys['testResult'] = u[1]
 93
 94#sparse solver:
 95simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
 96mbs.SolveStatic(simulationSettings)
 97
 98u = mbs.GetNodeOutput(nBodies-2, exu.OutputVariableType.Position) #tip node
 99exu.Print('static tip displacement (y)=', u[1])
100
101#factor 1e-2: 32bit version shows 2.1e-12 error
102exu.sys['testResult'] = 1e-2*u[1]
103
104if not testIsActive: 
105    SC.renderer.DoIdleTasks()
106    SC.renderer.Stop()