genericODE1duplicateNodeTest.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  The numerical Jacobian of ODE1 coordinates that an object addresses twice (#1424): an
 5#           ObjectGenericODE1 on the same NodeGenericODE1 twice, q_t = -q written as two halves, against the
 6#           same equations on the node once. The Jacobian is computed numerically for ObjectGenericODE1;
 7#           with each coordinate differentiated once, the linear problem converges in one Newton step per
 8#           time step in both models (a coordinate differentiated twice doubled its column: four steps).
 9#
10# Author:   Johannes Gerstmayr
11# Date:     2026-09-30
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
19
20import numpy as np
21
22testIsActive = exu.sys.get('testIsActive', False)
23
24testResult = 0
25for twice in [False, True]:
26    SC = exu.SystemContainer(); mbs = SC.AddSystem()
27    node = mbs.AddNode(NodeGenericODE1(referenceCoordinates=[0, 0], initialCoordinates=[1, 0.5], numberOfODE1Coordinates=2))
28    n = 4 if twice else 2
29    mbs.AddObject(ObjectGenericODE1(nodeNumbers=[node, node] if twice else [node],
30                                    systemMatrix=-np.eye(n)*(0.5 if twice else 1), rhsVector=[0]*n))
31    mbs.Assemble()
32    simulationSettings = exu.SimulationSettings()
33    simulationSettings.timeIntegration.numberOfSteps = 100
34    simulationSettings.timeIntegration.endTime = 1
35    simulationSettings.timeIntegration.newton.useModifiedNewton = False
36    simulationSettings.solution.file.write = False
37    solver = exu.MainSolverImplicitSecondOrder()
38    solver.SolveSystem(mbs, simulationSettings)
39    q = mbs.GetNodeOutput(node, exu.OutputVariableType.Coordinates)
40    exu.Print('node addressed', 'twice' if twice else 'once', ': q =', q, ', Newton steps', solver.it.newtonStepsCount)
41    testResult += sum(q) + solver.it.newtonStepsCount/100
42
43exu.Print('solution of genericODE1duplicateNodeTest=', testResult)
44exu.sys['testResult'] = testResult