netgenSTLtest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Example to import .stl mesh, mesh with Netgen, create FEM model,
  5#           reduced order CMS and simulate under gravity
  6#
  7# Author:   Johannes Gerstmayr
  8# Date:     2023-04-21
  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
 27
 28useGraphics = True
 29
 30#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 31#netgen/meshing part:
 32fem = FEMinterface()
 33
 34nModes = 16
 35
 36#steel:
 37rho = 7850
 38nu=0.3
 39Emodulus=1e8#use some very soft material to visualize deformations
 40
 41#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 42if True: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
 43    import sys
 44    import ngsolve as ngs
 45    import netgen
 46    from netgen.meshing import *
 47
 48    import netgen.stl as nstl
 49    #load STL file; needs to be closed (no holes) and consistent!
 50    #               and may not have defects (may require some processing of STL files!)
 51    geom = nstl.STLGeometry('testData/gyro.stl') #gyro bei Peter Manzl
 52
 53    mesh = ngs.Mesh( geom.GenerateMesh(maxh=0.003))
 54    # mesh.Curve(1) #don't do that!
 55
 56    #set True to see mesh in netgen tool:
 57    if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
 58        # import netgen
 59        import netgen.gui
 60        ngs.Draw(mesh)
 61        for i in range(10000000):
 62            netgen.Redraw() #this makes the netgen window interactive
 63            time.sleep(0.05)
 64
 65
 66    # sys.exit()
 67    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 68    #Use fem to import FEM model and create FFRFreducedOrder object
 69    [bfM, bfK, fes] = fem.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus,
 70                                                poissonsRatio=nu, meshOrder=1)
 71
 72
 73#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 74#compute Hurty-Craig-Bampton modes
 75if True: #now import mesh as mechanical model to EXUDYN
 76    print("nNodes=",fem.NumberOfNodes())
 77
 78
 79    cyl=np.array([0,0,0])
 80    rCyl = 0.011/2
 81    nodesOnCyl = fem.GetNodesOnCylinder(cyl-[0,0.01,0], cyl+[0,0.01,0], radius=rCyl, tolerance=0.001)
 82    # #print("boundary nodes bolt=", nodesOnBolt)
 83    nodesOnCylWeights = fem.GetNodeWeightsFromSurfaceAreas(nodesOnCyl)
 84    pMid = fem.GetNodePositionsMean(nodesOnCyl)
 85    print('cyl midpoint=', pMid)
 86
 87
 88    #boundaryList = [nodesOnBolt, nodesOnBolt, nodesOnBushing] #for visualization, use first interface twice
 89    boundaryList = [nodesOnCyl]
 90
 91    print("compute HCB modes... (may take some seconds)")
 92    fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
 93                                  nEigenModes=nModes,
 94                                  useSparseSolver=True,
 95                                  computationMode = HCBstaticModeSelection.RBE2)
 96
 97    print("eigen freq.=", fem.GetEigenFrequenciesHz())
 98
 99    #draw cylinder to see geometry of hole
100    # gGround = [GraphicsDataCylinder([0,0,0],[0,0.02,0], radius=0.011/2, color=color4dodgerblue, nTiles=128)]
101    # oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData=gGround)))
102
103
104    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
105    #compute stress modes for postprocessing (inaccurate for coarse meshes, just for visualization):
106    if True:
107        mat = KirchhoffMaterial(Emodulus, nu, rho)
108        varType = exu.OutputVariableType.StressLocal
109        print("ComputePostProcessingModes ... (may take a while)")
110        start_time = time.time()
111        fem.ComputePostProcessingModesNGsolve(fes, material=mat,
112                                       outputVariableType=varType)
113        SC.visualizationSettings.contour.reduceRange=False
114        SC.visualizationSettings.contour.outputVariable = varType
115        SC.visualizationSettings.contour.outputVariableComponent = -1 #norm
116    else:
117        varType = exu.OutputVariableType.DisplacementLocal
118        SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
119        SC.visualizationSettings.contour.outputVariableComponent = 0
120
121    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
122    print("create CMS element ...")
123    cms = ObjectFFRFreducedOrderInterface(fem)
124
125    objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
126                                                  initialVelocity=[0,0,0],
127                                                  initialAngularVelocity=[0,0,0],
128                                                  color=[0.9,0.9,0.9,1.],
129                                                  gravity=[0,0,-9.81]
130                                                  )
131
132    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
133    #add markers and joints
134    nodeDrawSize = 0.0005 #for joint drawing
135
136    #add constraint for cylinder
137    if True:
138
139        oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
140
141        altApproach = True
142        mCyl = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
143                                                      meshNodeNumbers=np.array(nodesOnCyl), #these are the meshNodeNumbers
144                                                      weightingFactors=nodesOnCylWeights))
145
146        #due to meshing effects and weighting, the center point is not exactly at [0,1.5,0] as intended ...
147        pm0 = mbs.GetMarkerOutput(mCyl, exu.OutputVariableType.Position,exu.ConfigurationType.Reference)
148        print('marker0 ref pos=', pm0)
149
150        mGroundCyl = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround,
151                                                    localPosition=pm0,
152                                                    visualization=VMarkerBodyRigid(show=True)))
153        mbs.AddObject(GenericJoint(markerNumbers=[mGroundCyl, mCyl],
154                                    constrainedAxes = [1]*6,
155                                    visualization=VGenericJoint(show=False, axesRadius=0.01, axesLength=0.01)))
156
157
158    if False: #activate to animate modes
159        from exudyn.interactive import AnimateModes
160        mbs.Assemble()
161        SC.visualizationSettings.nodes.show = False
162        SC.visualizationSettings.openGL.showFaceEdges = True
163        SC.visualizationSettings.openGL.multiSampling=4
164        SC.visualizationSettings.openGL.lineWidth=2
165        SC.visualizationSettings.window.renderWindowSize = [1600,1080]
166        SC.visualizationSettings.contour.showColorBar = False
167        SC.visualizationSettings.general.textSize = 16
168
169        #%%+++++++++++++++++++++++++++++++++++++++
170        #animate modes of ObjectFFRFreducedOrder (only needs generic node containing modal coordinates)
171        SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved
172
173        nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
174        AnimateModes(SC, mbs, nodeNumber, period=0.1, showTime=False, renderWindowText='Hurty-Craig-Bampton: 2 x 6 static modes and 8 eigenmodes\n',
175                     runOnStart=True)
176        # import sys
177        # sys.exit()
178
179    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
180    #animate modes
181    SC.visualizationSettings.markers.show = True
182    SC.visualizationSettings.markers.defaultSize=nodeDrawSize
183    SC.visualizationSettings.markers.drawSimplified = False
184
185    SC.visualizationSettings.loads.show = False
186
187    SC.visualizationSettings.openGL.multiSampling=4
188    SC.visualizationSettings.openGL.lineWidth=2
189
190    mbs.Assemble()
191
192    simulationSettings = exu.SimulationSettings()
193
194    SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
195    SC.visualizationSettings.nodes.drawNodesAsPoint = False
196    SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
197
198    SC.visualizationSettings.nodes.show = False
199    SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
200    SC.visualizationSettings.nodes.basisSize = 0.12
201    SC.visualizationSettings.bodies.deformationScaleFactor = 100 #use this factor to scale the deformation of modes
202
203    SC.visualizationSettings.openGL.showFaceEdges = True
204    SC.visualizationSettings.openGL.showFaces = True
205
206    SC.visualizationSettings.sensors.show = True
207    SC.visualizationSettings.sensors.drawSimplified = False
208    SC.visualizationSettings.sensors.defaultSize = 0.01
209
210    h=2e-5 #make small to see some oscillations
211    tEnd = 5
212
213    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
214    simulationSettings.timeIntegration.endTime = tEnd
215    simulationSettings.solutionSettings.writeSolutionToFile = True
216    simulationSettings.timeIntegration.verboseMode = 1
217    simulationSettings.timeIntegration.simulateInRealtime = True
218    simulationSettings.timeIntegration.realtimeFactor = 0.01
219    simulationSettings.timeIntegration.newton.useModifiedNewton = True
220
221    simulationSettings.solutionSettings.sensorsWritePeriod = h
222
223    simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
224    #simulationSettings.displayStatistics = True
225    simulationSettings.displayComputationTime = True
226
227    #create animation:
228    # simulationSettings.solutionSettings.recordImagesInterval = 0.005
229    # SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
230    SC.visualizationSettings.window.renderWindowSize=[1920,1080]
231    SC.visualizationSettings.openGL.multiSampling = 4
232
233    useGraphics=True
234    if True:
235        if useGraphics:
236            SC.visualizationSettings.general.autoFitScene=False
237
238            exu.StartRenderer()
239            if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
240
241            mbs.WaitForUserToContinue() #press space to continue
242
243        #SC.RedrawAndSaveImage()
244        if True:
245            # mbs.SolveDynamic(solverType=exu.DynamicSolverType.TrapezoidalIndex2,
246            #                   simulationSettings=simulationSettings)
247            mbs.SolveDynamic(simulationSettings=simulationSettings)
248        else:
249            mbs.SolveStatic(simulationSettings=simulationSettings)
250
251        # uTip = mbs.GetSensorValues(sensTipDispl)[1]
252        # print("nModes=", nModes, ", tip displacement=", uTip)
253
254        if varType == exu.OutputVariableType.StressLocal:
255            mises = CMSObjectComputeNorm(mbs, 0, exu.OutputVariableType.StressLocal, 'Mises')
256            print('max von-Mises stress=',mises)
257
258        if useGraphics:
259            SC.WaitForRenderEngineStopFlag()
260            exu.StopRenderer() #safely close rendering window!
261
262        if False:
263
264            mbs.PlotSensor(sensorNumbers=[sensBushingVel], components=[1])
265
266#%%
267if False:
268    import matplotlib.pyplot as plt
269    import matplotlib.ticker as ticker
270    import exudyn as exu
271    from exudyn.utilities import *
272    CC = PlotLineCode
273    comp = 1 #1=x, 2=y, ...
274    var = ''
275    # data = np.loadtxt('solution/hingePartBushing'+var+'2.txt', comments='#', delimiter=',')
276    # plt.plot(data[:,0], data[:,comp], CC(7), label='2 eigenmodes')
277    # data = np.loadtxt('solution/hingePartBushing'+var+'4.txt', comments='#', delimiter=',')
278    # plt.plot(data[:,0], data[:,comp], CC(8), label='4 eigenmodes')
279    data = np.loadtxt('solution/hingePartBushing'+var+'8.txt', comments='#', delimiter=',')
280    plt.plot(data[:,0], data[:,comp], CC(9), label='8 eigenmodes')
281    data = np.loadtxt('solution/hingePartBushing'+var+'16.txt', comments='#', delimiter=',')
282    plt.plot(data[:,0], data[:,comp], CC(10), label='16 eigenmodes')
283    data = np.loadtxt('solution/hingePartBushing'+var+'32.txt', comments='#', delimiter=',')
284    plt.plot(data[:,0], data[:,comp], CC(11), label='32 eigenmodes')
285
286    data = np.loadtxt('solution/hingePartBushing'+var+'2HCB.txt', comments='#', delimiter=',')
287    plt.plot(data[:,0], data[:,comp], CC(1), label='HCB + 2 eigenmodes')
288    data = np.loadtxt('solution/hingePartBushing'+var+'4HCB.txt', comments='#', delimiter=',')
289    plt.plot(data[:,0], data[:,comp], CC(2), label='HCB + 4 eigenmodes')
290    data = np.loadtxt('solution/hingePartBushing'+var+'8HCB.txt', comments='#', delimiter=',')
291    plt.plot(data[:,0], data[:,comp], CC(3), label='HCB + 8 eigenmodes')
292    data = np.loadtxt('solution/hingePartBushing'+var+'16HCB.txt', comments='#', delimiter=',')
293    plt.plot(data[:,0], data[:,comp], CC(4), label='HCB + 16 eigenmodes')
294    data = np.loadtxt('solution/hingePartBushing'+var+'32HCB.txt', comments='#', delimiter=',')
295    plt.plot(data[:,0], data[:,comp], CC(5), label='HCB + 32 eigenmodes')
296    data = np.loadtxt('solution/hingePartBushing'+var+'64HCB.txt', comments='#', delimiter=',')
297    plt.plot(data[:,0], data[:,comp], CC(6), label='HCB + 64 eigenmodes')
298    data = np.loadtxt('solution/hingePartBushing'+var+'128HCB.txt', comments='#', delimiter=',')
299    plt.plot(data[:,0], data[:,comp], CC(7), label='HCB + 128 eigenmodes')
300
301
302    ax=plt.gca() # get current axes
303    ax.grid(True, 'major', 'both')
304    ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
305    ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
306    #
307    plt.xlabel("time (s)")
308    plt.ylabel("y-component of tip velocity of hinge (m)")
309    plt.legend() #show labels as legend
310    plt.tight_layout()
311    plt.show()