geometricallyExactBeamJacobianTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  The Jacobian of ObjectBeamGeometricallyExact against a numerical Jacobian (#1100, #1550): one
  5#           element in a bent and rotating configuration, with Euler parameter nodes and with Tait-Bryan
  6#           nodes. The stiffness part is compared on the tangent space of the Euler parameters (the part
  7#           normal to the unit sphere is fixed by the norm constraint), the velocity part - the quadratic
  8#           velocity vector - with a thick cross section. Then a cantilever with Tait-Bryan nodes under a
  9#           large tip force and torque, solved statically with the analytic Jacobian. The result is the
 10#           tip position plus the number of Jacobian comparisons that pass.
 11#
 12# Author:   Johannes Gerstmayr
 13# Date:     2026-09-29
 14#
 15# 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.
 16#
 17#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 18
 19import exudyn as exu
 20from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 21
 22import numpy as np
 23
 24testIsActive = exu.sys.get('testIsActive', False)
 25
 26E = 1e8; rho = 1000; h = 0.002; b = 0.01; nu = 0.3
 27A = b*h; Izz = b*h**3/12; Iyy = h*b**3/12; G = E/(2*(1+nu))
 28
 29
 30def Section(inertiaFactor=1):
 31    section = exu.BeamSection()
 32    section.stiffnessMatrix = np.diag([E*A, G*A, G*A, G*(Iyy+Izz), E*Iyy, E*Izz])
 33    section.inertia = inertiaFactor*np.diag([rho*(Iyy+Izz), rho*Iyy, rho*Izz])
 34    section.massPerLength = rho*A
 35    return section
 36
 37
 38def Node(nodeType, p, R, omega):
 39    if nodeType == 'EP':
 40        ep = RotationMatrix2EulerParameters(R)
 41        return NodeRigidBodyEP(referenceCoordinates=list(p)+list(ep),
 42                               initialVelocities=[0, 0, 0]+list(AngularVelocity2EulerParameters_t(omega, ep)))
 43    rot = RotationMatrix2RotXYZ(R)
 44    return NodeRigidBodyRxyz(referenceCoordinates=list(p)+list(rot),
 45                             initialVelocities=[0, 0, 0]+list(AngularVelocity2RotXYZ_t(omega, rot)))
 46
 47
 48def Jacobian(nodeType, numerical, velocity, inertiaFactor):
 49    """the ODE2 Jacobian of one bent, rotating element: of the coordinates, or of the velocities"""
 50    SC = exu.SystemContainer(); mbs = SC.AddSystem()
 51    L = 0.5
 52    n0 = mbs.AddNode(Node(nodeType, [0, 0, 0], RotationMatrixY(0.1) @ RotationMatrixZ(0.2), [1, 2, 3]))
 53    n1 = mbs.AddNode(Node(nodeType, [0.98*L, 0.05, -0.03], RotationMatrixZ(0.5) @ RotationMatrixX(0.3), [-2, 1, 4]))
 54    mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[n0, n1], length=L, sectionData=Section(inertiaFactor)))
 55    mbs.Assemble()
 56    s = exu.SimulationSettings()
 57    s.timeIntegration.newton.numericalDifferentiation.forODE2 = numerical
 58    s.timeIntegration.newton.numericalDifferentiation.relativeEpsilon = 1e-6 if velocity else 1e-7
 59    s.linearSolver.solverType = exu.LinearSolverType.EXUdense
 60    solver = exu.MainSolverImplicitSecondOrder()
 61    solver.InitializeSolver(mbs, s)
 62    n = mbs.systemData.ODE2Size()
 63    solver.ComputeJacobianODE2RHS(mbs, scalarFactor_ODE2=0 if velocity else 1, scalarFactor_ODE2_t=1 if velocity else 0,
 64                                  scalarFactor_ODE1=0)
 65    J = solver.GetSystemJacobian()[:n, :n]
 66    #tangent space: remove the direction of the Euler parameters themselves
 67    T = np.eye(n)
 68    if nodeType == 'EP':
 69        for (node, i0) in [(n0, 3), (n1, 10)]:
 70            ep = np.array(mbs.GetNodeParameter(node, 'referenceCoordinates')[3:])
 71            T[i0:i0+4, i0:i0+4] -= np.outer(ep, ep)
 72    return T @ J @ T
 73
 74
 75passed = 0
 76for nodeType in ['EP', 'Rxyz']:
 77    for (velocity, inertiaFactor) in [(False, 1), (True, 1e6)]:
 78        JAnalytic = Jacobian(nodeType, False, velocity, inertiaFactor)
 79        JNumerical = Jacobian(nodeType, True, velocity, inertiaFactor)
 80        difference = np.linalg.norm(JAnalytic-JNumerical)/np.linalg.norm(JNumerical)
 81        exu.Print(nodeType, 'velocities' if velocity else 'coordinates', ': relative difference analytic - numerical Jacobian', difference)
 82        passed += difference < 1e-5
 83
 84
 85#a cantilever with Tait-Bryan nodes, large tip force and torque, static
 86SC = exu.SystemContainer(); mbs = SC.AddSystem()
 87L = 1; nElements = 8; lElement = L/nElements
 88nodes = [mbs.AddNode(NodeRigidBodyRxyz(referenceCoordinates=[i*lElement, 0, 0, 0, 0, 0])) for i in range(nElements+1)]
 89for i in range(nElements):
 90    mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[nodes[i], nodes[i+1]], length=lElement, sectionData=Section()))
 91mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=mbs.AddNode(NodePointGround()), coordinate=0))
 92for i in range(6):
 93    mbs.AddObject(ObjectConnectorCoordinate(markerNumbers=[mGround, mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nodes[0], coordinate=i))]))
 94mTip = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodes[-1]))
 95mbs.AddLoad(LoadForceVector(markerNumber=mTip, loadVector=[0, -2*E*Izz/L**2, 0.01*E*Iyy/L**2]))
 96mbs.AddLoad(LoadTorqueVector(markerNumber=mTip, loadVector=[0.3*G*(Iyy+Izz)/L, 0, 0]))
 97mbs.Assemble()
 98s = exu.SimulationSettings()
 99s.linearSolver.solverType = exu.LinearSolverType.EigenSparse
100s.staticSolver.numberOfLoadSteps = 5
101s.staticSolver.newton.relativeTolerance = 1e-10
102solver = exu.MainSolverStatic()
103solver.SolveSystem(mbs, s)
104pTip = mbs.GetNodeOutput(nodes[-1], exu.OutputVariableType.Position)
105exu.Print('cantilever: tip', pTip.round(8), ', Newton steps', solver.it.newtonStepsCount)
106
107testResult = passed + sum(pTip)
108exu.Print('solution of geometricallyExactBeamJacobianTest=', testResult)
109exu.sys['testResult'] = testResult