abaqusImportTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for FEMinterface with different ABAQUS files;
  5#           also test if store/load for FEMinterface and ObjectFFRFreducedOrderInterface works
  6#
  7# Author:   Johannes Gerstmayr
  8# Date:     2020-05-13
  9#
 10# 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.
 11#
 12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 13
 14import exudyn as exu
 15from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 16import exudyn.graphics as graphics #only import if it does not conflict
 17from exudyn.FEM import *
 18
 19useGraphics = True #without test
 20#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 21#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 22try: #only if called from test suite
 23    from modelUnitTests import exudynTestGlobals #for globally storing test results
 24    useGraphics = exudynTestGlobals.useGraphics
 25except:
 26    class ExudynTestGlobals:
 27        pass
 28    exudynTestGlobals = ExudynTestGlobals()
 29#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 30
 31useGraphics=False
 32
 33import numpy as np
 34import time
 35import sys
 36
 37#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 38#Use FEMinterface to import FEM model and create FFRFreducedOrder object
 39fileDir = 'testData/abaqus/block'
 40result = 0
 41
 42elements = ['C3D4','C3D10','C3D8','C3D20','C3D20R']
 43#elements = ['C3D8']
 44
 45for element in elements:
 46    SC = exu.SystemContainer()
 47    mbs = SC.AddSystem()
 48
 49    fem = FEMinterface()
 50    inputFileName = fileDir+element
 51    nodes=fem.ImportFromAbaqusInputFile(inputFileName+'.inp', typeName='Instance', name='rotor-1')
 52
 53    fem.ReadMassMatrixFromAbaqus(inputFileName+'_MASS1.mtx')
 54    fem.ReadStiffnessMatrixFromAbaqus(inputFileName+'_STIF1.mtx')
 55    if True:
 56        fn = 'solution/testFEM'
 57        #test save and import
 58        if np.__version__ <= '2.0': #load save does not work yet!
 59            fem.SaveToFile(fn)
 60            fem = FEMinterface()
 61            fem.LoadFromFile(fn)
 62        else: #for newer numpy version, use pickle (PKL) or HDF5
 63            fem.SaveToFile(fn,mode='PKL')
 64            fem = FEMinterface()
 65            fem.LoadFromFile(fn,mode='PKL')
 66
 67
 68        if False:
 69            exu.Print('size of nodes:', sys.getsizeof(np.array(fem.nodes['Position'])) )
 70            exu.Print('size of elements:', sys.getsizeof(fem.elements[0]) )
 71            exu.Print('size of massMatrix:', sys.getsizeof(fem.massMatrix) )
 72            exu.Print('size of stiffnessMatrix:', sys.getsizeof(fem.stiffnessMatrix) )
 73            exu.Print('size of modeBasis:', sys.getsizeof(fem.modeBasis) )
 74            #print('size of postProcessingModes:', sys.getsizeof(fem.postProcessingModes['matrix']) )
 75            exu.Print('===================')
 76
 77    nModes = 5 #use 2,5,6,7 or 9 but not 8, as in case that symmetric modes are swapped (in Hex case), solution is completely different
 78
 79    pLeft = [0,0,0]
 80    pLeftMid = [0,0.5,0.5]
 81    pRight = [4,0,0]
 82    nTip = fem.GetNodeAtPoint(pRight) #tip node (do not use midpoint, as this may not be a mesh node ...)
 83
 84    nodesLeftPlane = fem.GetNodesInPlane(pLeft, [-1,0,0])
 85    weightsLeftPlane = fem.GetNodeWeightsFromSurfaceAreas(nodesLeftPlane)
 86
 87    nodesRightPlane = fem.GetNodesInPlane(pRight, [-1,0,0])
 88    weightsRightPlane = fem.GetNodeWeightsFromSurfaceAreas(nodesRightPlane)
 89
 90    boundaryList = [nodesLeftPlane]
 91
 92    if useGraphics:
 93        exu.Print("nNodes=",fem.NumberOfNodes())
 94        exu.Print("compute HCB modes... ")
 95    start_time = time.time()
 96    fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
 97                                  nEigenModes=nModes,
 98                                  useSparseSolver=False, #sparse solver gives non-repeatable results ...
 99                                  computationMode = HCBstaticModeSelection.RBE2)
100
101    if useGraphics:
102        exu.Print("HCB modes needed %.3f seconds" % (time.time() - start_time))
103
104    cms = ObjectFFRFreducedOrderInterface(fem)
105    if True: #try save/load
106        fn = 'solution/testCMS'
107        cms.SaveToFile(fn)
108        cms = ObjectFFRFreducedOrderInterface()
109        cms.LoadFromFile(fn)
110
111        if False: #check size of objects
112            for key in cms.__dict__:
113                exu.Print('size of ', key, ':', sys.getsizeof(cms.__dict__[key]) )
114
115    objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
116                                            initialVelocity=[0,0,0],
117                                            initialAngularVelocity=[0,0,0],
118                                            gravity=[0,-9.81,0],
119                                            color=[0.1,0.9,0.1,1.])
120
121    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
122    #add markers and joints
123    nodeDrawSize = 0.0025 #for joint drawing
124
125    mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
126    oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
127
128    mGroundPosLeft = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=pLeftMid))
129    mLeft = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
130                                                  meshNodeNumbers=np.array(nodesLeftPlane), #these are the meshNodeNumbers
131                                                  weightingFactors=weightsLeftPlane))
132
133
134
135    mbs.AddObject(GenericJoint(markerNumbers=[mGroundPosLeft, mLeft], constrainedAxes=[1,1,1, 1,1,1]))
136
137    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
138    sDisp=mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'], meshNodeNumber=nTip, #meshnode number!
139                             storeInternal=True,
140                             outputVariableType = exu.OutputVariableType.DisplacementLocal))
141
142    mbs.Assemble()
143
144    simulationSettings = exu.SimulationSettings()
145
146    # SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
147    SC.visualizationSettings.nodes.drawNodesAsPoint = False
148    SC.visualizationSettings.bodies.deformationScaleFactor = 1e4 #use this factor to scale the deformation of modes
149
150    SC.visualizationSettings.loads.drawSimplified = False
151
152    SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
153    SC.visualizationSettings.contour.outputVariableComponent = 1 #y-component
154
155    # simulationSettings.solutionSettings.solutionInformation = "ObjectFFRFreducedOrder test"
156
157    h=1e-4
158    tEnd = 5e-3 #at almost max. of deflection under gravity
159
160    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
161    simulationSettings.timeIntegration.endTime = tEnd
162    simulationSettings.solutionSettings.solutionWritePeriod = h
163    simulationSettings.timeIntegration.verboseMode = useGraphics
164    #simulationSettings.timeIntegration.verboseModeFile = 3
165    simulationSettings.timeIntegration.newton.useModifiedNewton = True
166
167    simulationSettings.solutionSettings.sensorsWritePeriod = h
168    simulationSettings.solutionSettings.coordinatesSolutionFileName = "solution/coordinatesSolutionCMStest.txt"
169    simulationSettings.solutionSettings.writeSolutionToFile=False
170
171    simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5 #SHOULD work with 0.9 as well
172
173    if useGraphics:
174        exu.StartRenderer()
175        if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
176
177        mbs.WaitForUserToContinue() #press space to continue
178
179    mbs.SolveDynamic(simulationSettings)
180
181    # data = np.loadtxt(fileDir+'nMidDisplacementCMS'+str(nModes)+'Test.txt', comments='#', delimiter=',')
182    data = mbs.GetSensorStoredData(sDisp)
183    exu.Print('u-tip for '+element+' = ', data[-1,1:], ', nNodes=',fem.NumberOfNodes())
184    result += abs(data[-1,1:]).sum()
185
186    if useGraphics:
187        SC.WaitForRenderEngineStopFlag()
188        exu.StopRenderer() #safely close rendering window!
189        lastRenderState = SC.GetRenderState() #store model view for next simulation
190
191exu.Print('solution of abaqusImportTest=',result)
192
193exudynTestGlobals.testError = (result - (0.0005885208722206333))
194exudynTestGlobals.testResult = result
195
196#for small meshes in TestModels:
197# u-tip for C3D4 =  [-1.39753280e-05 -8.83250776e-05  9.86454888e-07] , nNodes= 214
198# u-tip for C3D10 =  [-1.73664007e-05 -1.04237155e-04  1.86678663e-10] , nNodes= 257
199# u-tip for C3D8 =  [-1.70545446e-05 -1.03215074e-04  6.31348275e-08] , nNodes= 176
200# u-tip for C3D20 =  [-1.72124324e-05 -1.04326514e-04 -9.10211089e-11] , nNodes= 171
201# u-tip for C3D20R =  [-1.72278715e-05 -1.04863540e-04  5.83371018e-08] , nNodes= 171
202
203#for larger files (see Experimental folder):
204# u-tip for C3D4 =  [-1.39753280e-05 -8.83250776e-05  9.86454888e-07] , nNodes= 214
205# u-tip for C3D10 =  [-1.74686596e-05 -1.05689104e-04  5.05186882e-08] , nNodes= 1332
206# u-tip for C3D8 =  [-1.72150782e-05 -1.03377335e-04  4.39361119e-08] , nNodes= 176
207# u-tip for C3D20 =  [-1.74978842e-05 -1.05476601e-04 -1.35045610e-08] , nNodes= 600
208# u-tip for C3D20R =  [-1.72957461e-05 -1.05412212e-04  5.01834245e-08] , nNodes= 600