manualExplicitIntegrator.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  ANCF Cable2D cantilever bent with manual explicit integrator
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2020-01-08
  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.itemInterface import *
 15
 16useGraphics = True #without test
 17#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 18#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 19try: #only if called from test suite
 20    from modelUnitTests import exudynTestGlobals #for globally storing test results
 21    useGraphics = exudynTestGlobals.useGraphics
 22except:
 23    class ExudynTestGlobals:
 24        pass
 25    exudynTestGlobals = ExudynTestGlobals()
 26#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 27
 28SC = exu.SystemContainer()
 29mbs = SC.AddSystem()
 30
 31
 32exu.Print("\n\n++++++++++++++++++++++++++\nStart EXUDYN version "+exu.GetVersionString()+"\n")
 33
 34#background
 35rect = [-2,-2,2,2] #xmin,ymin,xmax,ymax
 36background0 = {'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
 37background1 = {'type':'Circle', 'radius': 0.1, 'position': [-1.5,0,0]}
 38background2 = {'type':'Text', 'position': [-1,-1,0], 'text':'Example with text\nin two lines:.=!'} #background
 39oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background0, background1, background2])))
 40
 41#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 42#cable:
 43mypi = 3.141592653589793
 44
 45L=2.                   # length of ANCF element in m
 46#L=mypi                 # length of ANCF element in m
 47E=2.07e11*1e-5              # Young's modulus of ANCF element in N/m^2
 48rho=7800               # density of ANCF element in kg/m^3
 49b=0.1                  # width of rectangular ANCF element in m
 50h=0.1                  # height of rectangular ANCF element in m
 51A=b*h                  # cross sectional area of ANCF element in m^2
 52I=b*h**3/12            # second moment of area of ANCF element in m^4
 53f=3*E*I/L**2           # tip load applied to ANCF element in N
 54
 55exu.Print("load f="+str(f))
 56exu.Print("EI="+str(E*I))
 57
 58nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 59mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 60
 61cableList=[]
 62
 63
 64
 65nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
 66nElements = 4 #for tests use nElements = 4
 67lElem = L / nElements
 68nLast = 0
 69for i in range(nElements):
 70    nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
 71    elem=mbs.AddObject(Cable2D(physicsLength=lElem, physicsMassPerLength=rho*A,
 72                               physicsBendingStiffness=E*I, physicsAxialStiffness=E*A*0.1,
 73                               nodeNumbers=[int(nc0)+i,int(nc0)+i+1], useReducedOrderIntegration=True))
 74    cableList+=[elem]
 75
 76mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=0))
 77mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=1))
 78mANCF2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=3))
 79
 80#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
 81#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
 82#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2]))
 83
 84mANCFLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=nLast)) #force
 85mbs.AddLoad(Force(markerNumber = mANCFLast, loadVector = [0, -1000, 0])) #will be changed in load steps
 86#mANCFrigid = mbs.AddMarker(MarkerBodyRigid(bodyNumber=elem, localPosition=[lElem,0,0])) #local position L = beam tip
 87#mbs.AddLoad(Torque(markerNumber = mANCFrigid, loadVector = [0, 0, E*I*0.25*mypi]))
 88#mANCFnode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nLast)) #local position L = beam tip
 89#mbs.AddLoad(Torque(markerNumber = mANCFnode, loadVector = [0, 0, 0.4*E*I*0.25*mypi]))
 90#mbs.AddLoad(Force(markerNumber = mANCFnode, loadVector = [0, 0.4*E*I*0.25*mypi,0]))
 91
 92
 93
 94mbs.Assemble()
 95#exu.Print(mbs)
 96
 97simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 98
 99
100#SC.visualizationSettings.bodies.showNumbers = False
101SC.visualizationSettings.nodes.defaultSize = 0.025
102dSize=0.01
103SC.visualizationSettings.bodies.defaultSize = [dSize, dSize, dSize]
104
105#simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-9
106simulationSettings.staticSolver.verboseMode = 1
107simulationSettings.staticSolver.verboseModeFile = 2
108simulationSettings.solutionSettings.solverInformationFileName = 'solution/solverInformation.txt'
109
110#simulationSettings.staticSolver.newton.absoluteTolerance = 1e-8
111simulationSettings.staticSolver.newton.relativeTolerance = 1e-6 #1e-5 works for 64 elements
112simulationSettings.staticSolver.newton.maxIterations = 20 #50 for bending into circle
113
114if useGraphics: #only start graphics once, but after background is set
115    exu.StartRenderer()
116
117simulationSettings.staticSolver.numberOfLoadSteps = 10
118simulationSettings.staticSolver.adaptiveStep = True
119
120import numpy as np
121
122testRefVal = 0
123#compute eigenvalues manually:
124calcEig = True
125if calcEig:
126    from scipy.linalg import solve, eigh, eig #eigh for symmetric matrices, positive definite
127
128    staticSolver = exu.MainSolverStatic()
129    #staticSolver.SolveSystem(mbs, simulationSettings)
130
131    staticSolver.InitializeSolver(mbs, simulationSettings)
132
133    staticSolver.ComputeMassMatrix(mbs)
134    m = staticSolver.GetSystemMassMatrix()
135    #exu.Print("m =",m)
136
137    staticSolver.ComputeJacobianODE2RHS(mbs, scalarFactor_ODE2=-1, scalarFactor_ODE2_t=0)
138    staticSolver.ComputeJacobianAE(mbs)
139    K = staticSolver.GetSystemJacobian()
140    #exu.Print("K =",K)
141    nODE2 = staticSolver.GetODE2size()
142
143
144    K2 = K[0:nODE2,0:nODE2]
145
146    [eigvals, eigvecs] = eigh(K2, m) #this gives omega^2 ... squared eigen frequencies (rad/s)
147    ev = np.sort(a=abs(eigvals))
148    #exu.Print("ev =",ev)
149    if (len(ev) >= 7):
150        f6 = np.sqrt(abs(ev[6]))/(2*np.pi)
151        exu.Print("ev=", f6)
152        testRefVal += f6 #first bending eigenmode
153
154    staticSolver.FinalizeSolver(mbs, simulationSettings)
155
156#++++++++++++++++++++++++++++++++++++++++++++++++++
157#TEST
158def UserFunctionInitializeStep(mainSolver, mainSys, sims):
159    #exu.Print("t=", mainSolver.it.currentTime)
160    mainSolver.UpdateCurrentTime(mainSys, sims)
161    mainSys.systemData.SetTime(mainSolver.it.currentTime);
162    return True
163
164#test for explicit integrator:
165def UserFunctionNewton(mainSolver, mainSys, sims):
166
167    nODE2 = mainSolver.GetODE2size()
168    nAE = mainSolver.GetAEsize()
169    #nSys = nODE2+nAE
170    #print("u=", mainSys.systemData.GetODE2Coordinates())
171    dynamicSolver.ComputeODE2RHS(mbs)
172    res = dynamicSolver.GetSystemResidual()
173    Fode2 = res[0:nODE2]
174    #print("res=", Fode2)
175
176    dynamicSolver.ComputeMassMatrix(mbs)
177    M = dynamicSolver.GetSystemMassMatrix()
178    a = np.linalg.solve(M,Fode2) #acceleration
179
180    h = dynamicSolver.it.currentStepSize
181
182    u0 = mainSys.systemData.GetODE2Coordinates()
183    v0 = mainSys.systemData.GetODE2Coordinates_t()
184
185    mainSys.systemData.SetODE2Coordinates(u0+h*v0)
186    mainSys.systemData.SetODE2Coordinates_t(v0+h*a)
187
188    return True
189
190dynamicSolver = exu.MainSolverImplicitSecondOrder()
191
192simulationSettings.timeIntegration.numberOfSteps = 5000 #1000 steps for test suite/error
193simulationSettings.timeIntegration.endTime = 0.05              #1s for test suite / error
194
195simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
196#simulationSettings.displayComputationTime = True
197simulationSettings.timeIntegration.verboseMode = 1
198
199#dynamicSolver.SetUserFunctionInitializeStep(mbs, UserFunctionInitializeStep)
200dynamicSolver.SetUserFunctionNewton(mbs, UserFunctionNewton)
201
202dynamicSolver.SolveSystem(mbs, simulationSettings)
203#mbs.SolveDynamic(simulationSettings)
204
205uy=mbs.GetNodeOutput(nLast,exu.OutputVariableType.Position)[1] #y-coordinate of tip
206exu.Print("uy=", uy)
207exu.Print("testResult=", testRefVal + uy)
208exudynTestGlobals.testError = testRefVal + uy - (2.280183538481952-0.2204849087896498) #2020-01-16: 2.280183538481952-0.2204849087896498
209exudynTestGlobals.testResult = testRefVal + uy
210
211if useGraphics: #only start graphics once, but after background is set
212    SC.WaitForRenderEngineStopFlag()
213    exu.StopRenderer() #safely close rendering window!