NGsolveCMStest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for Hurty-Craig-Bampton modes using a simple flexible pendulum meshed with Netgen
  5#
  6# Author:   Johannes Gerstmayr 
  7# Date:     2021-04-20
  8# Update:   2022-07-11: runs now with pip installed ngsolve on Python 3.10
  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
 18import numpy as np
 19from exudyn.FEM import *
 20import os   #the mesh file is checked for existence, not loaded speculatively (#2469)
 21import time
 22
 23testIsActive = exu.sys.get('testIsActive', False)
 24
 25
 26SC = exu.SystemContainer()
 27mbs = SC.AddSystem()
 28
 29#for load/save of FEM data. The tracked file is an .npz of plain numpy arrays since 2026-09-17
 30#(#2469): the .pkl it replaces carried exudyn C++ types, so loading it under a second module
 31#(exudynCPPfast) raised 'type "Real" is already registered' and the model then regenerated the
 32#mesh - changing its own reference value. The .npz was CONVERTED from that .pkl, so the mesh and
 33#the reference value are unchanged. (The PKL/HDF5 round trip further below still tests those
 34#formats; it writes into the output directory, not here.)
 35fileName = 'testData/netgenTestMesh'
 36
 37#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 38#netgen/meshing part:
 39fem = FEMinterface()
 40#standard:
 41a = 0.1     #half height/width of beam
 42b = a
 43h = 2*a     #much too coarse, only for testing!!! 
 44L = 1       #Length of beam
 45nModes = 4
 46
 47rho = 1000
 48Emodulus=1e6 #smaller to see some deformation
 49nu=0.3
 50meshCreated = False
 51meshOrder = 1 #use order 2 for higher accuracy, but more unknowns
 52
 53hasNGsolve = False
 54
 55try:
 56    #netgen/ngsolve: https://github.com/NGSolve/ngsolve/releases
 57    import ngsolve as ngs
 58    from netgen.meshing import *
 59    from netgen.csg import *
 60    hasNGsolve = True
 61except:
 62    exu.Print('NGsolve not installed; trying to load mesh')
 63
 64#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 65if hasNGsolve: 
 66    
 67    geo = CSGeometry()
 68    
 69    block = OrthoBrick(Pnt(0,-a,-b),Pnt(L,a,b))
 70    geo.Add(block)
 71    
 72    mesh = ngs.Mesh( geo.GenerateMesh(maxh=h))
 73    mesh.Curve(1)
 74
 75    if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
 76        # import netgen
 77        import netgen.gui
 78        ngs.Draw (mesh)
 79        for i in range(10000000):
 80            netgen.Redraw() #this makes the window interactive
 81            time.sleep(0.05)
 82
 83    meshCreated = True
 84    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 85    #Use fem to import FEM model and create FFRFreducedOrder object
 86    [bfM, bfK, fes] = fem.ImportMeshFromNGsolve(mesh, density=rho, 
 87                                                youngsModulus=Emodulus, poissonsRatio=nu,
 88                                                meshOrder=meshOrder)
 89    #if file does not exist, create it - otherwise don't change it!
 90    #if you want to replace it, delete the old file!
 91    #
 92    #DECIDE ON THE FILE, NOT ON AN EXCEPTION (#2469): this used to save whenever the LOAD failed,
 93    #and a load fails for reasons that have nothing to do with the file being absent - under a
 94    #second C++ module it raises 'type "Real" is already registered', and a new ngsolve version
 95    #would do the same. The mesh file is TRACKED and the reference value belongs to it, so an
 96    #overwrite silently moved the result by 2.4e-8 and left the working tree dirty.
 97    if not os.path.isfile(fileName + '.npz'):
 98        exu.Print('\nNGsolveCMStest: ' + fileName + '.npz does not exist; creating it from the '
 99                  'mesh just generated\n')
100        fem.SaveToFile(fileName, mode='NPZ')
101
102#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
103#compute Hurty-Craig-Bampton modes
104#the file is the reference: it is never rewritten here, so a failure to load has to be reported
105#with its reason rather than repaired silently (#2469)
106if not os.path.isfile(fileName + '.npz'):
107    raise ValueError('NGsolveCMStest: ' + fileName + '.npz not found and ngsolve is not '
108                     'available to create it')
109try:
110    fem.LoadFromFile(fileName, mode='NPZ')
111except Exception as e:
112    exu.Print(f'\nNGsolveCMStest: LoadFromFile(...) failed; {e}\n')
113    raise ValueError('NGsolveCMStest: ' + fileName + '.npz exists but cannot be loaded (see the '
114                     'message above); the file is tracked and holds the reference solution, so it '
115                     'is NOT regenerated. Delete it deliberately to have it rebuilt.')
116
117    
118pLeft = [0,-a,-b]
119pRight = [L,-a,-b]
120nTip = fem.GetNodeAtPoint(pRight) #tip node (do not use midpoint, as this may not be a mesh node ...)
121
122nodesLeftPlane = fem.GetNodesInPlane(pLeft, [-1,0,0])
123weightsLeftPlane = fem.GetNodeWeightsFromSurfaceAreas(nodesLeftPlane)
124
125boundaryList = [nodesLeftPlane] 
126
127fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList, 
128                              nEigenModes=nModes, 
129                              useSparseSolver=True,
130                              computationMode = HCBstaticModeSelection.RBE2)
131
132
133#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
134#compute stress modes for postprocessing (inaccurate for coarse meshes, just for visualization):
135mat = KirchhoffMaterial(Emodulus, nu, rho)
136varType = exu.OutputVariableType.StressLocal
137start_time = time.time()
138fem.ComputePostProcessingModes(material=mat, outputVariableType=varType)
139SC.visualizationSettings.contour.reduceRange=False
140SC.visualizationSettings.contour.outputVariable = varType
141SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
142
143#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
144#now we test load/save for different formats!
145#the files of this round trip are OUTPUT of the test: write them into
146#exudyn.config.outputDirectory, which the test suite sets per model, so that two models cannot
147#write the same file (#2418); the mesh loaded further above is INPUT and is read as given
148from exudyn.basicUtilities import OutputFilePath
149fileName2 = OutputFilePath(fileName+'2', 'NGsolveCMStest')
150try:
151    fem.SaveToFile(fileName2+'2.npz')
152    fem.LoadFromFile(fileName2+'2', mode='NPZ')
153except:
154    #this should always work!
155    raise ValueError('NGsolveCMStest: load/save with NPZ mode failed!')
156
157try:
158    fem.SaveToFile(fileName2+'.pkl', mode='PKL')
159    fem.LoadFromFile(fileName2, mode='PKL')
160except:
161    #this should always work!
162    raise ValueError('NGsolveCMStest: load/save with PKL mode failed!')
163
164hasHDF5 = False
165try:
166    import h5py
167    hasHDF5 = True
168except:
169    exu.Print('NGsolveCMStest: import of h5py failed; to test, requires to pip install h5py')
170    
171try:
172    fem.SaveToFile(fileName2, mode='HDF5')
173    fem.LoadFromFile(fileName2+'.hdf5')
174except:
175    if hasHDF5:
176        raise ValueError('NGsolveCMStest: load/save with HDF5 mode failed even though h5py is installed!')
177        #pass
178    
179#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
180cms = ObjectFFRFreducedOrderInterface(fem)
181
182objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0], 
183                                              initialVelocity=[0,0,0], 
184                                              initialAngularVelocity=[0,0,0],
185                                              gravity=[0,-9.81,0],
186                                              color=[0.1,0.9,0.1,1.],
187                                              )
188
189
190#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
191#add markers and joints
192nodeDrawSize = 0.0025 #for joint drawing
193
194
195mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
196oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
197
198if True:
199    leftMidPoint = [0,0,0]
200    
201    mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=leftMidPoint))
202
203    mLeft = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'], 
204                                                  meshNodeNumbers=np.array(nodesLeftPlane), #these are the meshNodeNumbers
205                                                  weightingFactors=weightsLeftPlane))
206    mbs.AddObject(GenericJoint(markerNumbers=[mGround, mLeft], 
207                               constrainedAxes = [1,1,1,1,1,1*0],
208                               visualization=VGenericJoint(axesRadius=0.1*a, axesLength=0.1*a)))
209
210#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
211sensTipDispl = mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'], 
212                                                meshNodeNumber=nTip, #meshnode number!
213                                                storeInternal=True,
214                                                outputVariableType = exu.OutputVariableType.Displacement))
215    
216mbs.Assemble()
217
218simulationSettings = exu.SimulationSettings()
219
220SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
221SC.visualizationSettings.nodes.drawNodesAsPoint = False
222SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
223
224SC.visualizationSettings.nodes.show = False
225SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
226SC.visualizationSettings.nodes.basisSize = 0.12
227   
228SC.visualizationSettings.sensors.show = True
229SC.visualizationSettings.sensors.drawSimplified = False
230SC.visualizationSettings.sensors.defaultSize = 0.01
231SC.visualizationSettings.markers.drawSimplified = False
232SC.visualizationSettings.markers.show = False
233SC.visualizationSettings.markers.defaultSize = 0.01
234
235SC.visualizationSettings.loads.drawSimplified = False
236
237
238h=1e-3
239tEnd = 0.1
240if not testIsActive:
241    tEnd = 4
242
243simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
244simulationSettings.timeIntegration.endTime = tEnd
245simulationSettings.solution.file.write = False
246simulationSettings.timeIntegration.verboseMode = 1
247#simulationSettings.timeIntegration.verboseModeFile = 3
248simulationSettings.timeIntegration.newton.useModifiedNewton = True
249
250simulationSettings.solution.sensors.writePeriod = h
251
252simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
253simulationSettings.show.computationTime = True
254
255#create animation:
256SC.visualizationSettings.view0.window.renderWindowSize=[1920,1080]
257SC.visualizationSettings.openGL.multiSampling = 4
258
259if not testIsActive:
260    SC.visualizationSettings.general.autoFitScene=False
261
262    SC.renderer.Start()
263    SC.renderer.RestoreSavedState() #load last model view
264
265    SC.renderer.DoIdleTasks() #press space to continue
266
267mbs.SolveDynamic(simulationSettings=simulationSettings)
268    
269uTip = mbs.GetSensorValues(sensTipDispl)
270exu.Print("nModes=", nModes, ", tip displacement=", uTip)
271
272    
273if not testIsActive:
274    SC.renderer.DoIdleTasks()
275    SC.renderer.Stop() #safely close rendering window!
276
277result = np.linalg.norm(uTip)
278exu.Print('solution of NGsolveCMStest=',result)
279
280#%%++++++++++++++++++++++++++++++++++++++++++++++++++++
281
282exu.sys['testResult'] = result
283
284
285#mbs.SolutionViewer()