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