NGsolveCraigBampton.py
You can view and download this file on Github: NGsolveCraigBampton.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test for Hurty-Craig-Bampton modes using a simple flexible pendulum meshed with Netgen
5#
6# Author: Johannes Gerstmayr
7# Date: 2021-04-20
8# Update: 2022-07-11: runs now with pip installed ngsolve on Python 3.10
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
14
15import exudyn as exu
16from exudyn.itemInterface import *
17from exudyn.utilities import *
18from exudyn.FEM import *
19from exudyn.graphicsDataUtilities import *
20
21SC = exu.SystemContainer()
22mbs = SC.AddSystem()
23
24import numpy as np
25import time
26#import timeit
27
28import exudyn.basicUtilities as eb
29import exudyn.rigidBodyUtilities as rb
30import exudyn.utilities as eu
31
32import numpy as np
33
34useGraphics = True
35fileName = 'testData/netgenBrick' #for load/save of FEM data
36
37
38
39#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
40#netgen/meshing part:
41fem = FEMinterface()
42#standard:
43a = 0.025 #height/width of beam
44b = a
45h = 0.5*a
46L = 1 #Length of beam
47nModes = 8
48
49# #coarse1:
50# a = 0.025 #height/width of beam
51# b = a
52# h = a
53# L = 1 #Length of beam
54# nModes = 2
55
56#plate:
57# a = 0.025 #height/width of beam
58# b = 0.4
59# L = 1 #Length of beam
60# h = 0.6*a
61# nModes = 40
62
63#for saving:
64# h = 1.25*a
65
66
67rho = 1000
68Emodulus=1e7*10
69nu=0.3
70meshCreated = False
71
72#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
73if True: #needs netgen/ngsolve to be installed with pip install; to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
74
75 import ngsolve as ngs
76 import netgen
77 from netgen.meshing import *
78
79 from netgen.geom2d import unit_square
80 #import netgen.libngpy as libng
81 from netgen.csg import *
82
83 geo = CSGeometry()
84
85 block = OrthoBrick(Pnt(0,-a,-b),Pnt(L,a,b))
86 geo.Add(block)
87
88 #Draw (geo)
89
90 mesh = ngs.Mesh( geo.GenerateMesh(maxh=h))
91 mesh.Curve(1)
92
93 if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
94 # import netgen
95 import netgen.gui
96 ngs.Draw (mesh)
97 for i in range(10000000):
98 netgen.Redraw() #this makes the window interactive
99 time.sleep(0.05)
100
101 meshCreated = True
102 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
103 #Use fem to import FEM model and create FFRFreducedOrder object
104 fem.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
105 if h == 1.25*a:
106 fem.SaveToFile(fileName)
107
108#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
109#compute Hurty-Craig-Bampton modes
110if not meshCreated: #now import mesh as mechanical model to EXUDYN
111 fem.LoadFromFile(fileName)
112
113if True:
114 pLeft = [0,-a,-b]
115 pRight = [L,-a,-b]
116 nTip = fem.GetNodeAtPoint(pRight) #tip node (do not use midpoint, as this may not be a mesh node ...)
117 #print("nMid=",nMid)
118 nodesLeftPlane = fem.GetNodesInPlane(pLeft, [-1,0,0])
119 lenNodesLeftPlane = len(nodesLeftPlane)
120 weightsLeftPlane = np.array((1./lenNodesLeftPlane)*np.ones(lenNodesLeftPlane))
121
122 nodesRightPlane = fem.GetNodesInPlane(pRight, [-1,0,0])
123 lenNodesRightPlane = len(nodesRightPlane)
124 weightsRightPlane = np.array((1./lenNodesRightPlane)*np.ones(lenNodesRightPlane))
125
126 #boundaryList = [nodesLeftPlane, nodesRightPlane] #second boudary (right plane) not needed ...
127 boundaryList = [nodesLeftPlane]
128
129 print("nNodes=",fem.NumberOfNodes())
130
131 print("compute HCB modes... ")
132 start_time = time.time()
133 fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
134 nEigenModes=nModes,
135 useSparseSolver=True,
136 computationMode = HCBstaticModeSelection.RBE2)
137 print("HCB modes needed %.3f seconds" % (time.time() - start_time))
138
139 #alternatives:
140 #fem.ComputeEigenModesWithBoundaryNodes(boundaryNodes=nodesLeftPlane, nEigenModes=nModes, useSparseSolver=False)
141 #fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
142 #print("eigen freq.=", fem.GetEigenFrequenciesHz())
143
144
145 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
146 #compute stress modes for postprocessing (inaccurate for coarse meshes, just for visualization):
147 if False:
148 mat = KirchhoffMaterial(Emodulus, nu, rho)
149 varType = exu.OutputVariableType.StressLocal
150 #varType = exu.OutputVariableType.StrainLocal
151 print("ComputePostProcessingModes ... (may take a while)")
152 start_time = time.time()
153 fem.ComputePostProcessingModes(material=mat,
154 outputVariableType=varType)
155 print(" ... needed %.3f seconds" % (time.time() - start_time))
156 SC.visualizationSettings.contour.reduceRange=False
157 SC.visualizationSettings.contour.outputVariable = varType
158 SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
159 else:
160 SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
161 SC.visualizationSettings.contour.outputVariableComponent = 1
162
163 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
164 print("create CMS element ...")
165 cms = ObjectFFRFreducedOrderInterface(fem)
166
167 objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
168 initialVelocity=[0,0,0],
169 initialAngularVelocity=[0,0,0],
170 gravity=[0,-9.81,0],
171 color=[0.1,0.9,0.1,1.],
172 )
173
174
175 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
176 #animate modes
177 if False:
178 from exudyn.interactive import AnimateModes
179 mbs.Assemble()
180 SC.visualizationSettings.nodes.show = False
181 SC.visualizationSettings.openGL.showFaceEdges = True
182 SC.visualizationSettings.openGL.multiSampling=4
183 #SC.visualizationSettings.window.renderWindowSize = [1600,1080]
184 # SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
185 # SC.visualizationSettings.contour.outputVariableComponent = 0 #component
186
187
188 #%%+++++++++++++++++++++++++++++++++++++++
189 #animate modes of ObjectFFRFreducedOrder (only needs generic node containing modal coordinates)
190 SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved
191
192 nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
193 AnimateModes(SC, mbs, nodeNumber)
194 import sys
195 sys.exit()
196
197
198 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
199 #add markers and joints
200 nodeDrawSize = 0.0025 #for joint drawing
201
202
203 mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
204 oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
205
206 if True:
207 nodesLeftPlane = fem.GetNodesInPlane(pLeft, [-1,0,0])
208 lenNodesLeftPlane = len(nodesLeftPlane)
209 weightsLeftPlane = np.array((1./lenNodesLeftPlane)*np.ones(lenNodesLeftPlane))
210 leftMidPoint = [0,0,0]
211 #print("nodes in plane=",nodesLeftPlane)
212
213 mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=leftMidPoint))
214
215 mLeft = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
216 meshNodeNumbers=np.array(nodesLeftPlane), #these are the meshNodeNumbers
217 weightingFactors=weightsLeftPlane))
218 mbs.AddObject(GenericJoint(markerNumbers=[mGround, mLeft],
219 constrainedAxes = [1,1,1,1,1,1*0],
220 visualization=VGenericJoint(axesRadius=0.1*a, axesLength=0.1*a)))
221
222 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
223 fileDir = 'solution/'
224 sensTipDispl = mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'], meshNodeNumber=nTip, #meshnode number!
225 fileName=fileDir+'nMidDisplacementCMS'+str(nModes)+'Test.txt',
226 outputVariableType = exu.OutputVariableType.Displacement))
227
228 # mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=[GraphicsDataOrthoCubeLines(0, 0, 0, 1, 1, 1)])))
229 # mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=[GraphicsDataOrthoCubeLines(0.2, 0.2, 0.2, 0.8, 0.8, 0.8)], color=color4red)))
230 mbs.Assemble()
231
232 simulationSettings = exu.SimulationSettings()
233
234 SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
235 SC.visualizationSettings.nodes.drawNodesAsPoint = False
236 SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
237
238 SC.visualizationSettings.nodes.show = False
239 SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
240 SC.visualizationSettings.nodes.basisSize = 0.12
241 SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
242
243 SC.visualizationSettings.openGL.showFaceEdges = True
244 SC.visualizationSettings.openGL.showFaces = True
245
246 SC.visualizationSettings.sensors.show = True
247 SC.visualizationSettings.sensors.drawSimplified = False
248 SC.visualizationSettings.sensors.defaultSize = 0.01
249 SC.visualizationSettings.markers.drawSimplified = False
250 SC.visualizationSettings.markers.show = False
251 SC.visualizationSettings.markers.defaultSize = 0.01
252
253 SC.visualizationSettings.loads.drawSimplified = False
254
255
256 simulationSettings.solutionSettings.solutionInformation = "ObjectFFRFreducedOrder test"
257
258 h=1e-3
259 tEnd = 4
260
261 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
262 simulationSettings.timeIntegration.endTime = tEnd
263 simulationSettings.solutionSettings.writeSolutionToFile = False
264 simulationSettings.timeIntegration.verboseMode = 1
265 #simulationSettings.timeIntegration.verboseModeFile = 3
266 simulationSettings.timeIntegration.newton.useModifiedNewton = True
267
268 simulationSettings.solutionSettings.sensorsWritePeriod = h
269
270 simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
271 #simulationSettings.displayStatistics = True
272 simulationSettings.displayComputationTime = True
273
274 #create animation:
275 # simulationSettings.solutionSettings.recordImagesInterval = 0.005
276 # SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
277 SC.visualizationSettings.window.renderWindowSize=[1920,1080]
278 SC.visualizationSettings.openGL.multiSampling = 4
279
280 useGraphics=True
281 if True:
282 if useGraphics:
283 SC.visualizationSettings.general.autoFitScene=False
284
285 exu.StartRenderer()
286 if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
287
288 mbs.WaitForUserToContinue() #press space to continue
289
290 if True:
291 mbs.SolveDynamic(#solverType=exu.DynamicSolverType.TrapezoidalIndex2,
292 simulationSettings=simulationSettings)
293 else:
294 mbs.SolveStatic(simulationSettings=simulationSettings)
295
296 uTip = mbs.GetSensorValues(sensTipDispl)[1]
297 print("nModes=", nModes, ", tip displacement=", uTip)
298
299 if False:
300 # SC.visualizationSettings.exportImages.saveImageFileName = "images/test"
301 SC.visualizationSettings.exportImages.saveImageFormat = "TXT"
302 SC.visualizationSettings.exportImages.saveImageAsTextTriangles=True
303 SC.RedrawAndSaveImage() #uses default filename
304
305 from exudyn.plot import LoadImage, PlotImage
306 data = LoadImage('images/frame00000.txt', trianglesAsLines=True)
307 #PlotImage(data)
308 PlotImage(data, HT=HomogeneousTransformation(RotationMatrixZ(0.*pi)@RotationMatrixX(0.*pi), [0,0,0]), lineWidths=0.5, lineStyles='-',
309 triangleEdgeColors='b', triangleEdgeWidths=0.1, title='', closeAll=True, plot3D=True)
310 # PlotImage(data, HT=HomogeneousTransformation(RotationMatrixZ(0.5*pi)@RotationMatrixX(0.5*pi), [0,0,0]), lineWidths=0.5, title='', closeAll=True, fileName='images/test.pdf')
311
312
313 if useGraphics:
314 SC.WaitForRenderEngineStopFlag()
315 exu.StopRenderer() #safely close rendering window!
316
317
318
319
320#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
321#convergence of static tip-displacement with free-free eigenmodes:
322# nModes= 2 , tip displacement= -0.020705764289813425
323# nModes= 4 , tip displacement= -0.028232935031474123
324# nModes= 8 , tip displacement= -0.03462347289851485
325# nModes= 12 , tip displacement= -0.03744041447559639
326# nModes= 20 , tip displacement= -0.0421200606030284
327# nModes= 32 , tip displacement= -0.045122957364252446
328# nModes= 50 , tip displacement= -0.04711202668188235
329# nModes= 80 , tip displacement= -0.049164046183158706
330# nModes= 120 , tip displacement= -0.050065649361566426
331# nModes= 200 , tip displacement= -0.05054314003738750
332#with correct boundary conditions:
333#nModes= 20 , tip displacement= -0.05254089450183475
334#with Hurty-Craig-Bampton RBE2 boundaries:
335#nModes= 2 , tip displacement= -0.05254074496775043
336#nModes= 8 , tip displacement= -0.05254074496762861