NGsolveFFRF.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for FFRF with Netgen/NGsolve, using OpenCascade (OCC) geometry;
  5#           this example also makes use of boundary conditions (BC) and multi-domain materials
  6#
  7# Author:   Johannes Gerstmayr 
  8# Date:     2025-06-12
  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
 17from math import sin
 18import exudyn.graphics as graphics
 19from exudyn.FEM import *
 20
 21SC = exu.SystemContainer()
 22mbs = SC.AddSystem()
 23
 24import time
 25
 26import numpy as np
 27
 28useGraphics = True
 29fileName = OutputFilePath('solution/netgenFFRF2', 'NGsolveFFRF') #generated output, not input: goes to the ignored solution/ (#2491)
 30
 31
 32#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 33#netgen/meshing part:
 34fem = FEMinterface()
 35#standard:
 36a = 0.12 #height/width of beam
 37b = a*0.75
 38meshSize = 0.5*a
 39L = 1     #Length of beam
 40nModes = 16
 41
 42
 43meshCreated = False
 44meshOrder = 1 #use order 2 for higher accuracy, but more unknowns
 45
 46#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 47if True: 
 48    #needs netgen/ngsolve to be installed with pip install; to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
 49
 50    import ngsolve as ngs
 51    from netgen.occ import *
 52
 53    #fillets lead to smaller mesh size
 54    addFillets = True
 55
 56    #we create a simple mechanism with a container and some filling    
 57    box0 = Box((0,-0.5*a,-0.5*b), (L,0.5*a,0.5*b))
 58    box1 = Box((L-b*0.5,-0.5*a,-0.5*b), (L+b*0.5,0.5*a,0.5*b+L))
 59    box2 = Box((0,-0.5*a,-0.5*b+L), (L,0.5*a,0.5*b+L))
 60    cyl0 = Cylinder(Pnt(a,0,-b), Z, r=0.44*a, h=2*b)
 61    cyl2 = Cylinder(Pnt(a,0,-b+L), Z, r=0.44*a, h=2*b)
 62    box = box0+box1+box2
 63    
 64    if addFillets:
 65        box = box.MakeFillet (box.edges, a*0.1)
 66        cyl0 = cyl0.MakeChamfer(cyl0.edges, a*0.05)
 67        cyl2 = cyl2.MakeChamfer(cyl2.edges, a*0.05)
 68    
 69    cyl0.faces.Min(Z).name='cyl0_minZ'
 70    cyl0.faces.Max(Z).name='cyl0_maxZ'
 71    cyl2.faces.Min(Z).name='cyl2_minZ'
 72    cyl2.faces.Max(Z).name='cyl2_maxZ'
 73
 74    
 75    t=0.04
 76    cylContainer = Cylinder(Pnt(L-0.5*b+0.5*t,0,0.5*L), -X, r=0.4*L, h=1.2*L) 
 77                    
 78    cylContainerInner = Cylinder(Pnt(L-0.5*b-0.5*t,0,0.5*L), -X, r=0.4*L-t, h=1.2*L)
 79    cylContainer = (cylContainer-cylContainerInner)
 80    if addFillets:
 81        cylContainer = cylContainer.MakeFillet(cylContainer.edges, 0.4*t)
 82    
 83    #filling (with different density and stiffness)
 84    cylContent = Cylinder(Pnt(L-0.5*b-0.5*t,0,0.5*L), -X, r=0.4*L-t*0.95, h=0.5*L)
 85    cylContent.name = 'filling'
 86    frame = box+cyl0+cyl2+cylContainer
 87    geo = OCCGeometry(Glue((frame,cylContent)) ) #to keep different materials
 88    
 89    mesh = ngs.Mesh( geo.GenerateMesh(maxh=meshSize,
 90                                      curvaturesafety=0.5*1, #0.5 is coarse, 2 is quite fine 
 91                                      grading=0.8, #default=0.3
 92                                      ))
 93
 94    print('elements=',mesh.ne)
 95    print('materials=',mesh.GetMaterials())
 96    # mesh.Curve(1)
 97
 98    boundariesList = ['cyl0_minZ','cyl0_maxZ','cyl2_minZ','cyl2_maxZ']
 99
100    #set this to true, if you want to visualize the mesh inside netgen/ngsolve
101    if False:
102        import netgen.gui #this starts netgen gui; Press button "Visual" and activate "Auto-redraw after (sec)"; Then select "Mesh"
103
104    #%%
105    #define materials; note: we make them softer to see some deformations
106    materials = {'default':{'youngsModulus':2.1e11*0.01, 'poissonsRatio':0.3, 'density':7800},
107                  'filling':{'youngsModulus':5e8, 'poissonsRatio':0.4, 'density':2800}, #very soft
108                  }
109
110    meshCreated = True
111    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
112    #Use fem to import FEM model and create FFRFreducedOrder object
113    [bfM, bfK, fes] = fem.ImportMeshFromNGsolve(mesh, 
114                                                materials=materials,
115                                                boundaryNamesList=boundariesList,
116                                                meshOrder=meshOrder)
117    fem.SaveToFile(fileName)
118
119#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
120#compute Hurty-Craig-Bampton modes
121if not meshCreated: #now import mesh as mechanical model to EXUDYN
122    fem.LoadFromFile(fileName)
123
124#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
125#create Exudyn FFRF model from NGsolve mesh
126
127boundaryNodesList = []
128boundaryWeightsList = []
129for nodeSet in fem.nodeSets:
130    boundaryNodesList.append(nodeSet['NodeNumbers'])
131    boundaryWeightsList.append(nodeSet['NodeWeights'])
132
133
134print("nNodes=",fem.NumberOfNodes())
135
136print("compute HCB modes... ")
137start_time = time.time()
138fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryNodesList, 
139                                  nEigenModes=nModes, 
140                                  useSparseSolver=True,
141                                  computationMode = HCBstaticModeSelection.RBE2)
142print("HCB modes needed %.3f seconds" % (time.time() - start_time))
143
144
145
146
147#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
148#compute stress modes for postprocessing (inaccurate for coarse meshes, just for visualization):
149if True:
150    #mat = KirchhoffMaterial( **materials['default'] ) #same material for all domains
151    mat = KirchhoffMaterial(materials=materials, fes=fes) #use multi-domain materials
152    varType = exu.OutputVariableType.StressLocal
153    #for strain-computation:
154    # mat = 0
155    # varType = exu.OutputVariableType.StrainLocal
156    print("ComputePostProcessingModes ... (may take a while)")
157    start_time = time.time()
158    fem.ComputePostProcessingModesNGsolve(fes, material=mat,
159                                   outputVariableType=varType)
160    print("   ... needed %.3f seconds" % (time.time() - start_time))
161    SC.visualizationSettings.contour.reduceRange=True
162    SC.visualizationSettings.contour.outputVariable = varType
163    SC.visualizationSettings.contour.outputVariableComponent = -1 #norm
164else:
165    SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
166    SC.visualizationSettings.contour.outputVariableComponent = -1 
167
168#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
169print("create CMS element ...")
170cms = ObjectFFRFreducedOrderInterface(fem)
171
172objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0], 
173                                        initialVelocity=[0,0,0], 
174                                        initialAngularVelocity=[0,0,0],
175                                        stiffnessProportionalDamping=1e-2,
176                                        gravity=[9.81,0,0],
177                                        color=[0.1,0.9,0.1,1.],
178                                        )
179
180boundaryMarkerList = []
181for i, nodeSet in enumerate(boundaryNodesList):
182    weights = boundaryWeightsList[i]
183    marker = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'], 
184                                                   meshNodeNumbers=np.array(nodeSet), #these are the meshNodeNumbers
185                                                   weightingFactors=weights))
186    boundaryMarkerList.append(marker)
187
188
189SC.visualizationSettings.markers.drawSimplified = False
190SC.visualizationSettings.markers.show = True
191SC.visualizationSettings.markers.defaultSize = 0.01
192SC.visualizationSettings.view0.window.renderWindowSize = [1200,1000]
193SC.visualizationSettings.openGL.multiSampling=4
194
195#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
196#animate modes
197if False:
198    from exudyn.interactive import AnimateModes
199    mbs.Assemble()
200    SC.visualizationSettings.nodes.show = False
201    SC.visualizationSettings.view0.scene.showFaceEdges = True
202    
203    
204    #%%+++++++++++++++++++++++++++++++++++++++
205    #animate modes of ObjectFFRFreducedOrder (only needs generic node containing modal coordinates)
206    SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved
207    
208    nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
209    AnimateModes(SC, mbs, nodeNumber)
210    import sys
211    sys.exit()
212
213
214#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
215#add markers and joints
216nodeDrawSize = 0.0025 #for joint drawing
217
218
219oGround = mbs.CreateGround(referencePosition= [0,0,0],
220                           graphicsDataList=[graphics.CheckerBoard(point=[L*1.2,0,0.5*L],normal=[-1,0,0],size=4*L,
221                                                                   materialIndex=graphics.material.indexChrome)])
222
223for i, nodeSet in enumerate(boundaryNodesList):
224    markerFFRF = boundaryMarkerList[i]
225    markerGround = GetOtherMarker(mbs, bodyNumber=oGround, existingMarker=markerFFRF, show=False)
226    #print('pos',i,'=',mbs.GetMarkerOutput(markerGround,variableType=exu.OutputVariableType.Position,
227    #                                      configuration=exu.ConfigurationType.Reference))
228    mbs.AddObject(GenericJoint(markerNumbers=[markerGround, markerFFRF], 
229                                constrainedAxes = [1,1,1*(i==0), 0,0,0],
230                                visualization=VGenericJoint(axesRadius=0.1*a, axesLength=0.1*a)))
231
232torque = mbs.AddLoad(Torque(markerNumber=boundaryMarkerList[0],
233                            loadVector=[0,0,0]))
234
235#%%++++++++++++++++++++++++++++++
236#apply torque (only on one side):
237def PreStepUserFunction(mbs, t):
238    M = 7000*sin(pi*(t-1)/4) if (t < 5 and t>=1) else 0
239    mbs.SetLoadParameter(loadNumber=torque, parameterName='loadVector',value=[0,0,M])
240    return True
241
242mbs.SetPreStepUserFunction(PreStepUserFunction)    
243
244#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
245mbs.Assemble()
246
247simulationSettings = exu.SimulationSettings()
248
249SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
250SC.visualizationSettings.nodes.drawNodesAsPoint = False
251SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
252
253SC.visualizationSettings.nodes.show = False
254SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
255SC.visualizationSettings.nodes.basisSize = 0.12
256SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
257
258SC.visualizationSettings.view0.scene.showFaceEdges = False
259SC.visualizationSettings.view0.scene.showFaces = True
260SC.visualizationSettings.openGL.lineWidth = 2
261SC.visualizationSettings.openGL.light0.position = [-8,-5,1,1]
262SC.visualizationSettings.openGL.light0.shadow = 0.25
263
264SC.visualizationSettings.sensors.show = True
265SC.visualizationSettings.sensors.drawSimplified = False
266SC.visualizationSettings.sensors.defaultSize = 0.01
267
268SC.visualizationSettings.loads.drawSimplified = False
269SC.visualizationSettings.view0.scene.drawWorldBasis = True
270SC.visualizationSettings.view0.scene.worldBasisSize = 0.5
271
272
273stepSize=4e-3
274tEnd = 10
275
276simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
277simulationSettings.timeIntegration.endTime = tEnd
278# simulationSettings.solution.file.write = False
279simulationSettings.solution.file.writePeriod = 0.02
280simulationSettings.timeIntegration.verboseMode = 1
281simulationSettings.timeIntegration.newton.useModifiedNewton = True
282simulationSettings.solution.sensors.writePeriod = stepSize
283simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
284simulationSettings.show.computationTime = True
285
286if True: #set False if you do not like to do raytracing
287    #raytracing options
288    SC.visualizationSettings.general.useMultiThreadedRendering = False
289    SC.visualizationSettings.openGL.multiSampling = 1
290    SC.visualizationSettings.openGL.light1.enable = False
291    SC.visualizationSettings.raytracer.numberOfThreads = 16 #number of threads!
292    SC.visualizationSettings.view0.camera.useRaytracer = False #set True for raytracing
293    SC.visualizationSettings.openGL.lightModelAmbient = [0.5,0.5,0.5,1]
294    SC.visualizationSettings.raytracer.advanced.backgroundColorReflections = [0.3,0.3,0.3,1]
295    #SC.visualizationSettings.raytracer.advanced.searchTreeFactor = 8
296    SC.visualizationSettings.raytracer.imageSizeFactor=3 #for faster rendering
297    SC.visualizationSettings.contour.alphaTransparency = graphics.material.indexShiny
298    SC.visualizationSettings.general.graphicsUpdateInterval = 0.1 #for record images, avoid intermediate redraw
299    SC.visualizationSettings.exportImages.saveImageTimeOut = 200000
300    SC.visualizationSettings.raytracer.keepWindowActive= True
301
302    #adjust the material for the contour color
303    mat0=SC.renderer.materials.Get(5)
304    mat0.alpha=0.5
305    mat0.ior=1.15
306    mat0.reflectivity=0.15
307    SC.renderer.materials.Set(5,mat0)
308
309useGraphics=True
310if useGraphics:
311    SC.visualizationSettings.general.autoFitScene=False
312
313    SC.renderer.Start()
314    SC.renderer.RestoreSavedState() #load last model view
315
316    SC.renderer.DoIdleTasks() #press space to continue
317
318mbs.SolveDynamic(simulationSettings=simulationSettings)
319    
320    
321if useGraphics:
322    SC.renderer.DoIdleTasks()
323    SC.renderer.Stop() #safely close rendering window!
324
325
326mbs.SolutionViewer()