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()