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 * #includes itemInterface and rigidBodyUtilities
 18import exudyn.graphics as graphics
 19from exudyn.FEM import *
 20
 21SC = exu.SystemContainer()
 22mbs = SC.AddSystem()
 23
 24import numpy as np
 25import time
 26
 27
 28#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 29#netgen/meshing part:
 30fem = FEMinterface()
 31
 32nModes = 16
 33
 34#steel:
 35rho = 7850
 36nu=0.3
 37Emodulus=1e8#use some very soft material to visualize deformations
 38
 39#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 40if True: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
 41    import sys
 42    import ngsolve as ngs
 43    import netgen
 44    from netgen.meshing import *
 45
 46    print('load stl file...')
 47    import netgen.stl as nstl
 48    #load STL file; needs to be closed (no holes) and consistent!
 49    #               and may not have defects (may require some processing of STL files!)
 50    geom = nstl.STLGeometry('../Examples/testData/gyro.stl') #Peter's gyro
 51
 52    maxh=0.01
 53    mesh = ngs.Mesh( geom.GenerateMesh(maxh=maxh))
 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.gui #this starts netgen gui; Press button "Visual" and activate "Auto-redraw after (sec)"; Then select "Mesh"
 59
 60
 61    # sys.exit()
 62    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 63    #Use fem to import FEM model and create FFRFreducedOrder object
 64    [bfM, bfK, fes] = fem.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus, 
 65                                                poissonsRatio=nu, meshOrder=1)
 66
 67
 68#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 69#compute Hurty-Craig-Bampton modes
 70if True: #now import mesh as mechanical model to EXUDYN
 71    print("nNodes=",fem.NumberOfNodes())
 72
 73
 74    cyl=np.array([0,0,0])
 75    rCyl = 0.011/2
 76    nodesOnCyl = fem.GetNodesOnCylinder(cyl-[0,0.01,0], cyl+[0,0.01,0], radius=rCyl, tolerance=0.001)
 77    # #print("boundary nodes bolt=", nodesOnBolt)
 78    nodesOnCylWeights = fem.GetNodeWeightsFromSurfaceAreas(nodesOnCyl)
 79    pMid = fem.GetNodePositionsMean(nodesOnCyl)
 80    print('cyl midpoint=', pMid)
 81
 82
 83    #boundaryList = [nodesOnBolt, nodesOnBolt, nodesOnBushing] #for visualization, use first interface twice
 84    boundaryList = [nodesOnCyl] 
 85            
 86    print("compute HCB modes... (may take some seconds)")
 87    fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList, 
 88                                  nEigenModes=nModes, 
 89                                  useSparseSolver=True,
 90                                  computationMode = HCBstaticModeSelection.RBE2)
 91    
 92    print("eigen freq.=", fem.GetEigenFrequenciesHz())
 93
 94    #draw cylinder to see geometry of hole    
 95    # gGround = [graphics.Cylinder([0,0,0],[0,0.02,0], radius=0.011/2, color=graphics.color.dodgerblue, nTiles=128)]
 96    # oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData=gGround)))
 97    
 98        
 99    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
100    #compute stress modes for postprocessing (inaccurate for coarse meshes, just for visualization):
101    if True:
102        mat = KirchhoffMaterial(Emodulus, nu, rho)
103        varType = exu.OutputVariableType.StressLocal
104        print("ComputePostProcessingModes ... (may take a while)")
105        start_time = time.time()
106        fem.ComputePostProcessingModesNGsolve(fes, material=mat,
107                                       outputVariableType=varType)
108        SC.visualizationSettings.contour.reduceRange=False
109        SC.visualizationSettings.contour.outputVariable = varType
110        SC.visualizationSettings.contour.outputVariableComponent = -1 #norm
111    else:
112        varType = exu.OutputVariableType.DisplacementLocal
113        SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
114        SC.visualizationSettings.contour.outputVariableComponent = 0
115    
116    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
117    print("create CMS element ...")
118    cms = ObjectFFRFreducedOrderInterface(fem)
119    
120    objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0], 
121                                                  initialVelocity=[0,0,0], 
122                                                  initialAngularVelocity=[0,0,0],
123                                                  color=[0.9,0.9,0.9,1.],
124                                                  gravity=[0,0,-9.81] 
125                                                  )
126    
127    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
128    #add markers and joints
129    nodeDrawSize = 0.0005 #for joint drawing
130
131    #add constraint for cylinder
132    if True:
133        
134        oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
135
136        altApproach = True
137        mCyl = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'], 
138                                                      meshNodeNumbers=np.array(nodesOnCyl), #these are the meshNodeNumbers
139                                                      weightingFactors=nodesOnCylWeights))
140
141        #due to meshing effects and weighting, the center point is not exactly at [0,1.5,0] as intended ...
142        pm0 = mbs.GetMarkerOutput(mCyl, exu.OutputVariableType.Position,exu.ConfigurationType.Reference)
143        print('marker0 ref pos=', pm0)
144
145        mGroundCyl = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, 
146                                                    localPosition=pm0,
147                                                    visualization=VMarkerBodyRigid(show=True)))
148        mbs.AddObject(GenericJoint(markerNumbers=[mGroundCyl, mCyl], 
149                                    constrainedAxes = [1]*6,
150                                    visualization=VGenericJoint(show=False, axesRadius=0.01, axesLength=0.01)))
151
152    
153    if False: #activate to animate modes
154        from exudyn.interactive import AnimateModes
155        mbs.Assemble()
156        SC.visualizationSettings.nodes.show = False
157        SC.visualizationSettings.view0.scene.showFaceEdges = True
158        SC.visualizationSettings.openGL.multiSampling=4
159        SC.visualizationSettings.openGL.lineWidth=2
160        SC.visualizationSettings.view0.window.renderWindowSize = [1600,1080]
161        SC.visualizationSettings.contour.advanced.showColorBar = False
162        SC.visualizationSettings.view0.window.globalFontSize = 16
163        
164        #%%+++++++++++++++++++++++++++++++++++++++
165        #animate modes of ObjectFFRFreducedOrder (only needs generic node containing modal coordinates)
166        SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved
167        
168        nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
169        AnimateModes(SC, mbs, nodeNumber, period=0.1, showTime=False, renderWindowText='Hurty-Craig-Bampton: 2 x 6 static modes and 8 eigenmodes\n',
170                     runOnStart=True)
171        # import sys
172        # sys.exit()
173        
174    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
175    #animate modes
176    SC.visualizationSettings.markers.show = True
177    SC.visualizationSettings.markers.defaultSize=nodeDrawSize
178    SC.visualizationSettings.markers.drawSimplified = False
179
180    SC.visualizationSettings.loads.show = False
181
182    SC.visualizationSettings.openGL.multiSampling=4
183    SC.visualizationSettings.openGL.lineWidth=2
184    
185    mbs.Assemble()
186    
187    simulationSettings = exu.SimulationSettings()
188    
189    SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
190    SC.visualizationSettings.nodes.drawNodesAsPoint = False
191    SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
192    
193    SC.visualizationSettings.nodes.show = False
194    SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
195    SC.visualizationSettings.nodes.basisSize = 0.12
196    SC.visualizationSettings.bodies.deformationScaleFactor = 100 #use this factor to scale the deformation of modes
197    
198    SC.visualizationSettings.view0.scene.showFaceEdges = True
199    SC.visualizationSettings.view0.scene.showFaces = True
200    
201    SC.visualizationSettings.sensors.show = True
202    SC.visualizationSettings.sensors.drawSimplified = False
203    SC.visualizationSettings.sensors.defaultSize = 0.01
204    
205    h=2e-5 #make small to see some oscillations
206    tEnd = 0.5
207    
208    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
209    simulationSettings.timeIntegration.endTime = tEnd
210    simulationSettings.solution.file.write = False
211    simulationSettings.timeIntegration.verboseMode = 1
212    simulationSettings.timeIntegration.realtime.active = True
213    simulationSettings.timeIntegration.realtime.factor = 0.01
214    simulationSettings.timeIntegration.newton.useModifiedNewton = True
215    
216    simulationSettings.solution.sensors.writePeriod = h
217    
218    simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
219    #simulationSettings.show.statistics = True
220    simulationSettings.show.computationTime = True
221    
222    SC.visualizationSettings.view0.window.renderWindowSize=[1920,1080]
223    SC.visualizationSettings.openGL.multiSampling = 4
224
225    SC.visualizationSettings.general.autoFitScene=False
226
227    SC.renderer.Start()
228    SC.renderer.RestoreSavedState() #load last model view
229
230    SC.renderer.DoIdleTasks() #press space to continue
231
232    mbs.SolveDynamic(simulationSettings=simulationSettings)
233    
234    if varType == exu.OutputVariableType.StressLocal:
235        mises = CMSObjectComputeNorm(mbs, 0, exu.OutputVariableType.StressLocal, 'Mises')
236        print('max von-Mises stress=',mises)
237    
238    SC.renderer.DoIdleTasks()
239    SC.renderer.Stop() #safely close rendering window!
240    
241    # mbs.PlotSensor(sensorNumbers=[sensBushingVel], components=[1])