initialAccelerationsTest.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  The initial accelerations of the generalized-alpha method, timeIntegration.generalizedAlpha.
 5#           computeInitialAccelerations (#2235): a pendulum - a mass point on a distance constraint of length
 6#           L - starts horizontally with velocity v perpendicular to the rod. Its initial acceleration is the
 7#           centripetal v^2/L towards the pivot plus gravity; the term of the initial velocity comes from the
 8#           derivative of the constraint Jacobian, and the multiplier is the rod force m v^2/L. Without
 9#           computeInitialAccelerations the accelerations are zero.
10#
11# Author:   Johannes Gerstmayr
12# Date:     2026-10-05
13#
14# 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.
15#
16#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
17
18import exudyn as exu
19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
20import numpy as np
21
22testIsActive = exu.sys.get('testIsActive', False)
23
24L = 0.8
25v = 1.5
26g = 9.81
27
28def InitialAccelerations(computeInitialAccelerations, linearSolverType=exu.LinearSolverType.EXUdense):
29    SC = exu.SystemContainer()
30    mbs = SC.AddSystem()
31    oGround = mbs.CreateGround()
32    oMass = mbs.CreateMassPoint(referencePosition=[L,0,0], initialVelocity=[0,v,0], mass=2, gravity=[0,-g,0])
33    mbs.CreateDistanceConstraint(itemNumbers=[oGround, oMass])
34    mbs.Assemble()
35    settings = exu.SimulationSettings()
36    settings.linearSolver.solverType = linearSolverType
37    settings.solution.file.write = False
38    settings.timeIntegration.verboseMode = 0
39    settings.timeIntegration.generalizedAlpha.computeInitialAccelerations = computeInitialAccelerations
40    solver = exu.MainSolverImplicitSecondOrder()
41    solver.InitializeSolver(mbs, settings)
42    return (mbs.systemData.GetODE2Coordinates_tt(), mbs.systemData.GetAECoordinates())
43
44expected = np.array([-v**2/L, -g, 0])
45u = 0
46for linearSolverType in [exu.LinearSolverType.EXUdense, exu.LinearSolverType.EigenSparse]:
47    (a, multipliers) = InitialAccelerations(True, linearSolverType)
48    exu.Print(linearSolverType, ': initial accelerations', a, ', expected', expected, ', multiplier', multipliers)
49    u += int(np.linalg.norm(a - expected) < 1e-8) + int(abs(abs(multipliers[0]) - 2*v**2/L) < 1e-8) + a[0] + 0.1*a[1]
50(a, multipliers) = InitialAccelerations(False)
51exu.Print('without computeInitialAccelerations:', a)
52u += int(np.linalg.norm(a) == 0)
53
54exu.Print('solution of initialAccelerationsTest=', u)
55
56exu.sys['testResult'] = u