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
 86        import netgen.gui
 87        ngs.Draw (mesh)
 88        for i in range(10000000):
 89            netgen.Redraw() #this makes the window interactive
 90            time.sleep(0.05)
 91
 92    meshCreated = True
 93    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 94    #Use fem to import FEM model and create FFRFreducedOrder object
 95    [bfM, bfK, fes] = fem.ImportMeshFromNGsolve(mesh, density=rho,
 96                                                youngsModulus=Emodulus, poissonsRatio=nu,
 97                                                meshOrder=meshOrder)
 98    if h == 1.25*a:
 99        fem.SaveToFile(fileName)
100
101#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
102#compute Hurty-Craig-Bampton modes
103if not meshCreated: #now import mesh as mechanical model to EXUDYN
104    fem.LoadFromFile(fileName)
105
106if True:
107    pLeft = [0,-a,-b]
108    pRight = [L,-a,-b]
109    nTip = fem.GetNodeAtPoint(pRight) #tip node (do not use midpoint, as this may not be a mesh node ...)
110    #print("nMid=",nMid)
111    nodesLeftPlane = fem.GetNodesInPlane(pLeft, [-1,0,0])
112    # lenNodesLeftPlane = len(nodesLeftPlane)
113    # weightsLeftPlane = np.array((1./lenNodesLeftPlane)*np.ones(lenNodesLeftPlane))
114    weightsLeftPlane = fem.GetNodeWeightsFromSurfaceAreas(nodesLeftPlane)
115
116    nodesRightPlane = fem.GetNodesInPlane(pRight, [-1,0,0])
117    # lenNodesRightPlane = len(nodesRightPlane)
118    # weightsRightPlane = np.array((1./lenNodesRightPlane)*np.ones(lenNodesRightPlane))
119    weightsRightPlane = fem.GetNodeWeightsFromSurfaceAreas(nodesRightPlane)
120
121    #boundaryList = [nodesLeftPlane, nodesRightPlane] #second boudary (right plane) not needed ...
122    boundaryList = [nodesLeftPlane]
123
124    print("nNodes=",fem.NumberOfNodes())
125
126    print("compute HCB modes... ")
127    start_time = time.time()
128    fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
129                                  nEigenModes=nModes,
130                                  useSparseSolver=True,
131                                  computationMode = HCBstaticModeSelection.RBE2)
132    print("HCB modes needed %.3f seconds" % (time.time() - start_time))
133
134    #alternatives:
135    #fem.ComputeEigenModesWithBoundaryNodes(boundaryNodes=nodesLeftPlane, nEigenModes=nModes, useSparseSolver=False)
136    #fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
137    #print("eigen freq.=", fem.GetEigenFrequenciesHz())
138
139
140    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
141    #compute stress modes for postprocessing (inaccurate for coarse meshes, just for visualization):
142    if False:
143        mat = KirchhoffMaterial(Emodulus, nu, rho)
144        varType = exu.OutputVariableType.StressLocal
145        #varType = exu.OutputVariableType.StrainLocal
146        print("ComputePostProcessingModes ... (may take a while)")
147        start_time = time.time()
148        if True: #faster with ngsolve; requires fes
149            fem.ComputePostProcessingModesNGsolve(fes, material=mat,
150                                           outputVariableType=varType)
151        else:
152            fem.ComputePostProcessingModes(material=mat,
153                                            outputVariableType=varType)
154        print("   ... needed %.3f seconds" % (time.time() - start_time))
155        SC.visualizationSettings.contour.reduceRange=False
156        SC.visualizationSettings.contour.outputVariable = varType
157        SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
158    else:
159        SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
160        SC.visualizationSettings.contour.outputVariableComponent = 1
161
162    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
163    print("create CMS element ...")
164    cms = ObjectFFRFreducedOrderInterface(fem)
165
166    objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
167                                                  initialVelocity=[0,0,0],
168                                                  initialAngularVelocity=[0,0,0],
169                                                  gravity=[0,-9.81,0],
170                                                  color=[0.1,0.9,0.1,1.],
171                                                  )
172
173    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
174    #animate modes
175    if False:
176        from exudyn.interactive import AnimateModes
177        mbs.Assemble()
178        SC.visualizationSettings.nodes.show = False
179        SC.visualizationSettings.openGL.showFaceEdges = True
180        SC.visualizationSettings.openGL.multiSampling=4
181        #SC.visualizationSettings.window.renderWindowSize = [1600,1080]
182        # SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
183        # SC.visualizationSettings.contour.outputVariableComponent = 0 #component
184
185
186        #%%+++++++++++++++++++++++++++++++++++++++
187        #animate modes of ObjectFFRFreducedOrder (only needs generic node containing modal coordinates)
188        SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved
189
190        nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
191        AnimateModes(SC, mbs, nodeNumber)
192        import sys
193        sys.exit()
194
195
196    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
197    #add markers and joints
198    nodeDrawSize = 0.0025 #for joint drawing
199
200
201    mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
202    oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
203
204    if True:
205        leftMidPoint = [0,0,0]
206
207        mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=leftMidPoint))
208
209        mLeft = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
210                                                      meshNodeNumbers=np.array(nodesLeftPlane), #these are the meshNodeNumbers
211                                                      weightingFactors=weightsLeftPlane))
212        mbs.AddObject(GenericJoint(markerNumbers=[mGround, mLeft],
213                                   constrainedAxes = [1,1,1,1,1,1*0],
214                                   visualization=VGenericJoint(axesRadius=0.1*a, axesLength=0.1*a)))
215
216    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
217    fileDir = 'solution/'
218    sensTipDispl = mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'], meshNodeNumber=nTip, #meshnode number!
219                             fileName=fileDir+'nMidDisplacementCMS'+str(nModes)+'Test.txt',
220                             outputVariableType = exu.OutputVariableType.Displacement))
221
222    mbs.Assemble()
223
224    simulationSettings = exu.SimulationSettings()
225
226    SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
227    SC.visualizationSettings.nodes.drawNodesAsPoint = False
228    SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
229
230    SC.visualizationSettings.nodes.show = False
231    SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
232    SC.visualizationSettings.nodes.basisSize = 0.12
233    SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
234
235    SC.visualizationSettings.openGL.showFaceEdges = True
236    SC.visualizationSettings.openGL.showFaces = True
237
238    SC.visualizationSettings.sensors.show = True
239    SC.visualizationSettings.sensors.drawSimplified = False
240    SC.visualizationSettings.sensors.defaultSize = 0.01
241    SC.visualizationSettings.markers.drawSimplified = False
242    SC.visualizationSettings.markers.show = False
243    SC.visualizationSettings.markers.defaultSize = 0.01
244
245    SC.visualizationSettings.loads.drawSimplified = False
246
247
248    simulationSettings.solutionSettings.solutionInformation = "ObjectFFRFreducedOrder test"
249
250    h=1e-3
251    tEnd = 4
252
253    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
254    simulationSettings.timeIntegration.endTime = tEnd
255    # simulationSettings.solutionSettings.writeSolutionToFile = False
256    simulationSettings.timeIntegration.verboseMode = 1
257    #simulationSettings.timeIntegration.verboseModeFile = 3
258    simulationSettings.timeIntegration.newton.useModifiedNewton = True
259
260    simulationSettings.solutionSettings.sensorsWritePeriod = h
261
262    simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
263    #simulationSettings.displayStatistics = True
264    simulationSettings.displayComputationTime = True
265
266    #create animation:
267    # simulationSettings.solutionSettings.recordImagesInterval = 0.005
268    # SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
269    SC.visualizationSettings.window.renderWindowSize=[1920,1080]
270    SC.visualizationSettings.openGL.multiSampling = 4
271
272    useGraphics=True
273    if True:
274        if useGraphics:
275            SC.visualizationSettings.general.autoFitScene=False
276
277            exu.StartRenderer()
278            if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
279
280            mbs.WaitForUserToContinue() #press space to continue
281
282        if True:
283            mbs.SolveDynamic(#solverType=exu.DynamicSolverType.TrapezoidalIndex2,
284                              simulationSettings=simulationSettings)
285        else:
286            mbs.SolveStatic(simulationSettings=simulationSettings)
287
288        uTip = mbs.GetSensorValues(sensTipDispl)[1]
289        print("nModes=", nModes, ", tip displacement=", uTip)
290
291        if False:
292            # SC.visualizationSettings.exportImages.saveImageFileName = "images/test"
293            SC.visualizationSettings.exportImages.saveImageFormat = "TXT"
294            SC.visualizationSettings.exportImages.saveImageAsTextTriangles=True
295            SC.RedrawAndSaveImage() #uses default filename
296
297            from exudyn.plot import LoadImage, PlotImage
298            data = LoadImage('images/frame00000.txt', trianglesAsLines=True)
299            #PlotImage(data)
300            PlotImage(data, HT=HomogeneousTransformation(RotationMatrixZ(0.*pi)@RotationMatrixX(0.*pi), [0,0,0]), lineWidths=0.5, lineStyles='-',
301                      triangleEdgeColors='b', triangleEdgeWidths=0.1, title='', closeAll=True, plot3D=True)
302            # 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')
303
304
305        if useGraphics:
306            SC.WaitForRenderEngineStopFlag()
307            exu.StopRenderer() #safely close rendering window!
308
309
310        #mbs.SolutionViewer()
311
312
313#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
314#convergence of static tip-displacement with free-free eigenmodes:
315# nModes= 2 , tip displacement= -0.020705764289813425
316# nModes= 4 , tip displacement= -0.028232935031474123
317# nModes= 8 , tip displacement= -0.03462347289851485
318# nModes= 12 , tip displacement= -0.03744041447559639
319# nModes= 20 , tip displacement= -0.0421200606030284
320# nModes= 32 , tip displacement= -0.045122957364252446
321# nModes= 50 , tip displacement= -0.04711202668188235
322# nModes= 80 , tip displacement= -0.049164046183158706
323# nModes= 120 , tip displacement= -0.050065649361566426
324# nModes= 200 , tip displacement= -0.05054314003738750
325#with correct boundary conditions:
326#nModes= 20 , tip displacement= -0.05254089450183475
327#with Hurty-Craig-Bampton RBE2 boundaries:
328#nModes= 2 , tip displacement= -0.05254074496775043
329#nModes= 8 , tip displacement= -0.05254074496762861