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