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