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