compareFullModifiedNewton.py
You can view and download this file on Github: compareFullModifiedNewton.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: pendulum example: check difference of full and modified newton
5#
6# Author: Johannes Gerstmayr
7# Date: 2020-06-19
8#
9# 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.
10#
11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12
13import exudyn as exu
14from exudyn.utilities import *
15
16import numpy as np
17import matplotlib.pyplot as plt
18import matplotlib.ticker as ticker
19
20useGraphics = True #without test
21#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
22#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
23try: #only if called from test suite
24 from modelUnitTests import exudynTestGlobals #for globally storing test results
25 useGraphics = exudynTestGlobals.useGraphics
26except:
27 class ExudynTestGlobals:
28 pass
29 exudynTestGlobals = ExudynTestGlobals()
30#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31
32SC = exu.SystemContainer()
33mbs = SC.AddSystem()
34
35
36#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
37#rigid pendulum:
38rect = [-2,-2,2,2] #xmin,ymin,xmax,ymax
39background = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[rect[0],rect[1],0, rect[2],rect[1],0, rect[2],rect[3],0, rect[0],rect[3],0, rect[0],rect[1],0]} #background
40oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
41a = 0.5 #x-dim of pendulum
42b = 0.05 #y-dim of pendulum
43massRigid = 12
44inertiaRigid = massRigid/12*(2*a)**2
45g = 9.81 # gravity
46
47graphics2 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[-a,-b,0, a,-b,0, a,b,0, -a,b,0, -a,-b,0]} #background
48omega0 = 5*0 #inconsistent initial conditions lead to integration problems!
49nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[a,0,0], initialVelocities=[0,omega0*a,omega0]));
50oRigid = mbs.AddObject(RigidBody2D(physicsMass=massRigid, physicsInertia=inertiaRigid,nodeNumber=nRigid,visualization=VObjectRigidBody2D(graphicsData= [graphics2])))
51
52mR1 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[-a,0.,0.])) #support point
53mR2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ a,0.,0.])) #end point
54
55mG0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=[0,0,0.]))
56mbs.AddObject(RevoluteJoint2D(markerNumbers=[mG0,mR1]))
57
58mbs.AddLoad(Force(markerNumber = mR2, loadVector = [0, -massRigid*g, 0]))
59
60mbs.Assemble()
61#exu.Print(mbs)
62
63simulationSettings = exu.SimulationSettings() #takes currently set values or default values
64
65simulationSettings.timeIntegration.numberOfSteps = 100 #100
66simulationSettings.timeIntegration.endTime = 2
67simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8
68simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-4
69simulationSettings.timeIntegration.verboseMode = 1
70
71simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
72#simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
73#simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
74simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
75simulationSettings.displayStatistics = True
76simulationSettings.solutionSettings.solutionWritePeriod = 1e-4
77
78simulationSettings.timeIntegration.newton.useModifiedNewton = True
79simulationSettings.solutionSettings.coordinatesSolutionFileName = "solution/modifiedNewton.txt"
80mbs.SolveDynamic(simulationSettings)#, experimentalNewSolver=False)
81
82simulationSettings.timeIntegration.newton.useModifiedNewton = False
83simulationSettings.solutionSettings.coordinatesSolutionFileName = "solution/fullNewton.txt"
84mbs.SolveDynamic(simulationSettings)#, experimentalNewSolver=False)
85
86
87#%%*****************************************
88#post processing for mass point system
89
90dataM = np.loadtxt('solution/modifiedNewton.txt', comments='#', delimiter=',')
91dataF = np.loadtxt('solution/fullNewton.txt', comments='#', delimiter=',')
92
93# exu.Print("solFullNewton =",sum(abs(dataF[:,5])))
94# exu.Print("solModifiedNewton =",sum(abs(dataM[:,5])))
95
96u=sum(abs(dataM[:,5]-dataF[:,5]))
97exu.Print("compareFullModifiedNewton u=",u)
98
99exudynTestGlobals.testError = u - (0.0001583478719999567 ) #2020-12-18: 0.0001583478719999567
100exudynTestGlobals.testResult = u
101
102import os
103os.remove('solution/modifiedNewton.txt')
104os.remove('solution/fullNewton.txt')
105
106if useGraphics:
107 # plt.plot(dataM[:,0], dataM[:,3+2], 'b-') #plot column i over column 0 (time)
108 # plt.plot(dataF[:,0], dataF[:,3+2], 'r-') #plot column i over column 0 (time)
109 plt.plot(dataF[:,0], dataF[:,5]-dataM[:,5], 'r-')
110
111 ax=plt.gca() # get current axes
112 ax.grid(True, 'major', 'both')
113 ax.xaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
114 ax.yaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
115 plt.tight_layout()
116 plt.show()