MathematicalPendulumTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  A point mass on a massless rod of length L under gravity, modelled twice: once with a
  5#           `DistanceConstraint` and once with a stiff `SpringDamper`. Both are integrated over 2 s
  6#           and compared with their own reference; the result is the sum of the two errors.
  7#           The model compares against a reference value written into it, so its
  8#           result is that difference and its reference solution is 0 (#2632).
  9#
 10# Author:   Johannes Gerstmayr
 11# Date:     2019-11-01, reworked 2026-09-24
 12#
 13# 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.
 14#
 15#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 16
 17import exudyn as exu
 18from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 19import exudyn.graphics as graphics
 20
 21testIsActive = exu.sys.get('testIsActive', False)
 22exu.sys['testTolerance'] = 4e-13 #the tolerance RunAllModelUnitTests used for these ten
 23
 24SC = exu.SystemContainer()
 25mbs = SC.AddSystem()
 26
 27L = 0.8 #distance
 28nodeList = [0,0]
 29
 30for k in range(2):
 31    n1=mbs.AddNode(Point(referenceCoordinates = [L,0,0], initialCoordinates = [0,0,0]))
 32    nodeList[k] = n1
 33
 34    mass = 2.5
 35    g = 9.81
 36
 37    #add mass points and ground object:
 38    objectGround = mbs.AddObject(ObjectGround(referencePosition = [0,0,0]))
 39    massPoint = mbs.AddObject(MassPoint(mass = mass, nodeNumber = n1))
 40
 41    #marker for constraint / springDamper
 42    groundMarker = mbs.AddMarker(MarkerBodyPosition(bodyNumber = objectGround, localPosition= [0, 0, 0]))
 43    bodyMarker = mbs.AddMarker(MarkerBodyPosition(bodyNumber = massPoint, localPosition= [0, 0, 0]))
 44
 45    if k==0:
 46        mbs.AddObject(DistanceConstraint(distance = L, markerNumbers = [groundMarker,bodyMarker]))
 47    else:
 48        k = 40000 #spring stiffness
 49        d = 200  #damping coefficient
 50        mbs.AddObject(SpringDamper(stiffness = k, damping = d, force = 0, referenceLength = L, 
 51                                   markerNumbers = [groundMarker,bodyMarker]))
 52
 53    #add loads:
 54    mbs.AddLoad({'loadType': 'ForceVector',  'markerNumber': bodyMarker,  'loadVector': [0, -mass*g, 0]}) 
 55
 56#exu.Print(mbs)
 57mbs.Assemble()
 58if not testIsActive: 
 59    SC.renderer.Start()
 60
 61simulationSettings = exu.SimulationSettings()
 62
 63simulationSettings.solution.file.write=False
 64simulationSettings.timeIntegration.numberOfSteps = 1000
 65simulationSettings.timeIntegration.endTime = 2
 66
 67simulationSettings.timeIntegration.generalizedAlpha.useNewmark = False #better convergence with index2 and correct initial accelerations
 68simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = False
 69
 70#simulationSettings.timeIntegration.newton.relativeTolerance = 1e-6 #SHOULD work with standard values ...
 71simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
 72simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = True
 73simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6 
 74#simulationSettings.show.statistics = False
 75
 76mbs.SolveDynamic(simulationSettings)
 77if not testIsActive: 
 78    SC.renderer.DoIdleTasks()
 79    SC.renderer.Stop() #safely close rendering window!
 80
 81u1 = mbs.GetNodeOutput(nodeList[0], exu.OutputVariableType.Position) #tip node
 82exu.Print('solution mathematicalPendulum Constraint=',u1[1])
 83
 84u2 = mbs.GetNodeOutput(nodeList[1], exu.OutputVariableType.Position) #tip node
 85exu.Print('solution mathematicalPendulum SpringDamper=',u2[1])
 86if True:
 87    errorConstraint= u1[1] -   (-0.06808284314701757)   #2021-02-04: -0.06808284314701757
 88    errorSpringDamper= u2[1] - (-0.07148800507819238)   #2021-02-04: -0.07148800507819238
 89else:
 90    errorConstraint= u1[1] - (-0.0714264053422459)  #2019-12-26: -0.0714264053422459; 15.12.2019: -0.07242503089584812; before 15.12.2019: (-0.0724256565815142) #-(-0.7823882479152345) with endtime=10 and 10000 steps
 91    errorSpringDamper= u2[1] - (-0.07477852383438113) #2019-12-26: -0.07477852383438113; 15.12.2019: (-0.07579968609949864); before 15.12.2019: (-0.07579967194309412)  #-(-0.7738842923525007)
 92
 93testResult = abs(errorConstraint)+abs(errorSpringDamper)
 94exu.Print('solution of MathematicalPendulumTest=', testResult)
 95exu.sys['testResult'] = testResult
 96
 97#solution converged for 6 digits with index2 + accelerations initial conditions:
 98#100000
 99#solution mathematicalPendulum Constraint= -0.068039063184649
100#solution mathematicalPendulum SpringDamper= -0.071442303558492