solutionViewerMultipleSimulations.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for multiple static solutions merged into one solution file
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2022-12-20
  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 *
 15from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 16import exudyn.graphics as graphics #only import if it does not conflict
 17
 18SC = exu.SystemContainer()
 19mbs = SC.AddSystem()
 20
 21#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 22#set up simple ANCF model
 23#background
 24rect = [-0.5,-2,2.5,0.5] #xmin,ymin,xmax,ymax
 25background = {'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
 26oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
 27
 28#cable:
 29
 30L=2                    # length of ANCF element in m
 31E=2.07e11             # Young's modulus of ANCF element in N/m^2
 32rho=7800               # density of ANCF element in kg/m^3
 33b=0.1                  # width of rectangular ANCF element in m
 34h=0.1                  # height of rectangular ANCF element in m
 35A=b*h                  # cross sectional area of ANCF element in m^2
 36I=b*h**3/12            # second moment of area of ANCF element in m^4
 37f=2*3*E*I/L**2           # tip load applied to ANCF element in N
 38
 39nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 40mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 41
 42#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 43#generate ANCF beams with utilities function
 44cableTemplate = Cable2D(#physicsLength = L / nElements, #set in GenerateStraightLineANCFCable2D(...)
 45                        physicsMassPerLength = rho*A,
 46                        physicsBendingStiffness = E*I,
 47                        physicsAxialStiffness = E*A,
 48                        useReducedOrderIntegration = 1,
 49                        #nodeNumbers = [0, 0], #will be filled in GenerateStraightLineANCFCable2D(...)
 50                        )
 51
 52positionOfNode0 = [0, 0, 0] # starting point of line
 53positionOfNode1 = [L, 0, 0] # end point of line
 54numberOfElements = 16
 55
 56#alternative to mbs.AddObject(Cable2D(...)) with nodes:
 57ancf=GenerateStraightLineANCFCable2D(mbs,
 58                positionOfNode0, positionOfNode1,
 59                numberOfElements,
 60                cableTemplate, #this defines the beam element properties
 61                massProportionalLoad = [0,-9.81*0,0], #optionally add gravity
 62                fixedConstraintsNode0 = [1,1,0,1], #add constraints for pos and rot (r'_y)
 63                fixedConstraintsNode1 = [0,0,0,0])
 64mANCFLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=ancf[0][-1])) #ancf[0][-1] = last node
 65nLoad = mbs.AddLoad(Force(markerNumber = mANCFLast, loadVector = [f*0, -f, 0])) #will be changed in load steps
 66
 67
 68#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 69mbs.Assemble()
 70# print(mbs)
 71simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 72
 73simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/coordinatesSolution.txt'
 74simulationSettings.solutionSettings.writeSolutionToFile = True
 75simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/1000
 76simulationSettings.displayComputationTime = False
 77#simulationSettings.displayStatistics = True
 78#simulationSettings.displayComputationTime = True
 79
 80SC.visualizationSettings.nodes.defaultSize = 0.01
 81
 82simulationSettings.solutionSettings.solutionInformation = "Cantilever"
 83simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
 84
 85simulationSettings.staticSolver.verboseMode = 0
 86simulationSettings.staticSolver.newton.newtonResidualMode = 1
 87
 88#adapt these settings for better solution file with multiple simulations:
 89#**************************************************
 90simulationSettings.solutionSettings.appendToFile = False
 91simulationSettings.solutionSettings.writeFileFooter = False #never write footer as it would be seen between the solution steps
 92#**************************************************
 93
 94useGraphics=False
 95if useGraphics:
 96    exu.StartRenderer()
 97    mbs.WaitForUserToContinue()
 98
 99nLoadSteps = 25 #this is the number of individual computations; could also be done with staticSolver.numberOfLoadSteps
100                #  but here, we want to show how to do multiple steps merged into one solution file
101for loadSteps in range(nLoadSteps):
102    #loadValue = f**((loadSteps+1)/nLoadSteps) #geometric increment of loads
103    loadValue = 2*f*(loadSteps+1)/(nLoadSteps)
104
105    mbs.SetLoadParameter(nLoad, 'loadVector', [0, -loadValue,0])
106    #print('load vector=' + str(mbs.GetLoadParameter(nLoad, 'loadVector')) )
107
108    simulationSettings.staticSolver.loadStepStart = loadSteps
109    # simulationSettings.staticSolver.numberOfLoadSteps = 5
110    mbs.SolveStatic(simulationSettings, updateInitialValues=True)
111
112    #**************************************************
113    #after first STEP, add this:
114    simulationSettings.solutionSettings.writeInitialValues = False #to avoid duplication of output times (start/end)
115    simulationSettings.solutionSettings.writeFileHeader = False
116    simulationSettings.solutionSettings.appendToFile = True
117    #**************************************************
118
119    sol = mbs.systemData.GetODE2Coordinates()
120
121    n = len(sol)
122    print('load=',loadValue, ', tip: x='+str(sol[n-4])+', y='+str(sol[n-3]))
123
124if useGraphics:
125    SC.WaitForRenderEngineStopFlag()
126    exu.StopRenderer() #safely close rendering window!
127
128if True:
129    #%%
130
131    t=LoadSolutionFile('solution/coordinatesSolution.txt', verbose=False, safeMode=True)
132    mbs.SolutionViewer(solution=t)