genericODE2matrixFormatsTest.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  ObjectGenericODE2 with its mass, stiffness and damping matrices given dense and sparse, with the
 5#           dense and the sparse linear solver (#777): a chain of four masses on springs, fixed at one end and
 6#           pulled at the other, as one GenericODE2 on four 1D coordinates. All four combinations give the same
 7#           motion, and the static deflection is the one of the springs in series.
 8#
 9# Author:   Johannes Gerstmayr
10# Date:     2026-10-05
11#
12# 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.
13#
14#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15
16import exudyn as exu
17from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
18import numpy as np
19
20testIsActive = exu.sys.get('testIsActive', False)
21
22n = 4
23(mass, k, d, force) = (0.5, 1e3, 2., 10.)
24M = mass*np.eye(n)
25K = np.zeros((n,n))
26for i in range(n): #spring i between mass i-1 (or the wall) and mass i
27    K[i,i] += k
28    if i > 0:
29        K[i-1,i-1] += k; K[i-1,i] -= k; K[i,i-1] -= k
30D = d/k*K
31
32def Container(matrix, sparse):
33    container = exu.MatrixContainer()
34    if sparse:
35        triplets = [[i, j, matrix[i,j]] for i in range(n) for j in range(n) if matrix[i,j] != 0]
36        container.SetWithSparseMatrix(np.array(triplets), numberOfRows=n, numberOfColumns=n, useDenseMatrix=False)
37    else:
38        container.SetWithDenseMatrix(matrix, useDenseMatrix=True)
39    return container
40
41def Run(sparseMatrices, linearSolverType, static=False):
42    SC = exu.SystemContainer()
43    mbs = SC.AddSystem()
44    nodes = [mbs.AddNode(Node1D(referenceCoordinates=[0])) for i in range(n)]
45    mbs.AddObject(ObjectGenericODE2(nodeNumbers=nodes, massMatrix=Container(M, sparseMatrices),
46                                    stiffnessMatrix=Container(K, sparseMatrices), dampingMatrix=Container(D, sparseMatrices),
47                                    forceVector=[0]*(n-1)+[force]))
48    mbs.Assemble()
49    settings = exu.SimulationSettings()
50    settings.linearSolver.solverType = linearSolverType
51    settings.solution.file.write = False
52    settings.timeIntegration.verboseMode = 0
53    settings.staticSolver.verboseMode = 0
54    settings.timeIntegration.numberOfSteps = 200
55    settings.timeIntegration.endTime = 0.5
56    if static:
57        mbs.SolveStatic(settings)
58    else:
59        mbs.SolveDynamic(settings)
60    return mbs.systemData.GetODE2Coordinates()
61
62results = []
63for sparseMatrices in [False, True]:
64    for linearSolverType in [exu.LinearSolverType.EXUdense, exu.LinearSolverType.EigenSparse]:
65        results.append(Run(sparseMatrices, linearSolverType))
66        exu.Print('sparse matrices', sparseMatrices, ',', linearSolverType, ':', results[-1])
67difference = max(np.linalg.norm(result - results[0]) for result in results)
68static = Run(True, exu.LinearSolverType.EigenSparse, static=True)
69staticExpected = force/k*np.arange(1, n+1) #springs in series: each spring carries the force
70exu.Print('largest difference of the combinations', difference, ', static', static, ', expected', staticExpected)
71
72u = np.sum(results[0]) + int(difference < 1e-12) + int(np.linalg.norm(static - staticExpected) < 1e-8) #within the Newton tolerance
73exu.Print('solution of genericODE2matrixFormatsTest=', u)
74
75exu.sys['testResult'] = u