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 #only import if it does not conflict
18import numpy as np
19from exudyn.FEM import *
20import time
21
22useGraphics = True #without test
23#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
24#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
25try: #only if called from test suite
26 from modelUnitTests import exudynTestGlobals #for globally storing test results
27 useGraphics = exudynTestGlobals.useGraphics
28except:
29 class ExudynTestGlobals:
30 pass
31 exudynTestGlobals = ExudynTestGlobals()
32#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
33
34
35SC = exu.SystemContainer()
36mbs = SC.AddSystem()
37
38fileName = 'testData/netgenTestMesh' #for load/save of FEM data
39
40#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
41#netgen/meshing part:
42fem = FEMinterface()
43#standard:
44a = 0.1 #half height/width of beam
45b = a
46h = 2*a #much too coarse, only for testing!!!
47L = 1 #Length of beam
48nModes = 4
49
50rho = 1000
51Emodulus=1e6 #smaller to see some deformation
52nu=0.3
53meshCreated = False
54meshOrder = 1 #use order 2 for higher accuracy, but more unknowns
55
56hasNGsolve = False
57
58try:
59 #netgen/ngsolve: https://github.com/NGSolve/ngsolve/releases
60 import ngsolve as ngs
61 from netgen.meshing import *
62 from netgen.csg import *
63 hasNGsolve = True
64except:
65 exu.Print('NGsolve not installed; trying to load mesh')
66
67#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
68if hasNGsolve:
69
70 geo = CSGeometry()
71
72 block = OrthoBrick(Pnt(0,-a,-b),Pnt(L,a,b))
73 geo.Add(block)
74
75 mesh = ngs.Mesh( geo.GenerateMesh(maxh=h))
76 mesh.Curve(1)
77
78 if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
79 # import netgen
80 import netgen.gui
81 ngs.Draw (mesh)
82 for i in range(10000000):
83 netgen.Redraw() #this makes the window interactive
84 time.sleep(0.05)
85
86 meshCreated = True
87 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
88 #Use fem to import FEM model and create FFRFreducedOrder object
89 [bfM, bfK, fes] = fem.ImportMeshFromNGsolve(mesh, density=rho,
90 youngsModulus=Emodulus, poissonsRatio=nu,
91 meshOrder=meshOrder)
92 #if file does not exist, create it - otherwise don't change it!
93 #if you want to replace it, delete the old file!
94 try:
95 fem.LoadFromFile(fileName, mode='PKL')
96 except Exception as e:
97 exu.Print(f'\nNGsolveCMStest: LoadFromFile(...) failed; {e}\n')
98 fem.SaveToFile(fileName, mode='PKL')
99
100#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
101#compute Hurty-Craig-Bampton modes
102try:
103 fem.LoadFromFile(fileName, mode='PKL')
104except Exception as e:
105 exu.Print(f'\nNGsolveCMStest: LoadFromFile(...) failed; {e}\n')
106 raise ValueError('NGsolveCMStest: mesh file not found!')
107
108
109pLeft = [0,-a,-b]
110pRight = [L,-a,-b]
111nTip = fem.GetNodeAtPoint(pRight) #tip node (do not use midpoint, as this may not be a mesh node ...)
112
113nodesLeftPlane = fem.GetNodesInPlane(pLeft, [-1,0,0])
114weightsLeftPlane = fem.GetNodeWeightsFromSurfaceAreas(nodesLeftPlane)
115
116boundaryList = [nodesLeftPlane]
117
118fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
119 nEigenModes=nModes,
120 useSparseSolver=True,
121 computationMode = HCBstaticModeSelection.RBE2)
122
123
124#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
125#compute stress modes for postprocessing (inaccurate for coarse meshes, just for visualization):
126mat = KirchhoffMaterial(Emodulus, nu, rho)
127varType = exu.OutputVariableType.StressLocal
128start_time = time.time()
129fem.ComputePostProcessingModes(material=mat, outputVariableType=varType)
130SC.visualizationSettings.contour.reduceRange=False
131SC.visualizationSettings.contour.outputVariable = varType
132SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
133
134#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
135#now we test load/save for different formats!
136fileName2 = fileName+'2'
137try:
138 fem.SaveToFile(fileName2+'2.npz')
139 fem.LoadFromFile(fileName2+'2', mode='NPZ')
140except:
141 #this should always work!
142 raise ValueError('NGsolveCMStest: load/save with NPZ mode failed!')
143
144try:
145 fem.SaveToFile(fileName2+'.pkl', mode='PKL')
146 fem.LoadFromFile(fileName2, mode='PKL')
147except:
148 #this should always work!
149 raise ValueError('NGsolveCMStest: load/save with PKL mode failed!')
150
151hasHDF5 = False
152try:
153 import h5py
154 hasHDF5 = True
155except:
156 exu.Print('NGsolveCMStest: import of h5py failed; to test, requires to pip install h5py')
157
158try:
159 fem.SaveToFile(fileName2, mode='HDF5')
160 fem.LoadFromFile(fileName2+'.hdf5')
161except:
162 if hasHDF5:
163 raise ValueError('NGsolveCMStest: load/save with HDF5 mode failed even though h5py is installed!')
164 #pass
165
166#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
167cms = ObjectFFRFreducedOrderInterface(fem)
168
169objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
170 initialVelocity=[0,0,0],
171 initialAngularVelocity=[0,0,0],
172 gravity=[0,-9.81,0],
173 color=[0.1,0.9,0.1,1.],
174 )
175
176
177#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
178#add markers and joints
179nodeDrawSize = 0.0025 #for joint drawing
180
181
182mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
183oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
184
185if True:
186 leftMidPoint = [0,0,0]
187
188 mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=leftMidPoint))
189
190 mLeft = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
191 meshNodeNumbers=np.array(nodesLeftPlane), #these are the meshNodeNumbers
192 weightingFactors=weightsLeftPlane))
193 mbs.AddObject(GenericJoint(markerNumbers=[mGround, mLeft],
194 constrainedAxes = [1,1,1,1,1,1*0],
195 visualization=VGenericJoint(axesRadius=0.1*a, axesLength=0.1*a)))
196
197#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
198sensTipDispl = mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'],
199 meshNodeNumber=nTip, #meshnode number!
200 storeInternal=True,
201 outputVariableType = exu.OutputVariableType.Displacement))
202
203mbs.Assemble()
204
205simulationSettings = exu.SimulationSettings()
206
207SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
208SC.visualizationSettings.nodes.drawNodesAsPoint = False
209SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
210
211SC.visualizationSettings.nodes.show = False
212SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
213SC.visualizationSettings.nodes.basisSize = 0.12
214
215SC.visualizationSettings.sensors.show = True
216SC.visualizationSettings.sensors.drawSimplified = False
217SC.visualizationSettings.sensors.defaultSize = 0.01
218SC.visualizationSettings.markers.drawSimplified = False
219SC.visualizationSettings.markers.show = False
220SC.visualizationSettings.markers.defaultSize = 0.01
221
222SC.visualizationSettings.loads.drawSimplified = False
223
224
225h=1e-3
226tEnd = 0.1
227if useGraphics:
228 tEnd = 4
229
230simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
231simulationSettings.timeIntegration.endTime = tEnd
232simulationSettings.solutionSettings.writeSolutionToFile = False
233simulationSettings.timeIntegration.verboseMode = 1
234#simulationSettings.timeIntegration.verboseModeFile = 3
235simulationSettings.timeIntegration.newton.useModifiedNewton = True
236
237simulationSettings.solutionSettings.sensorsWritePeriod = h
238
239simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
240simulationSettings.displayComputationTime = True
241
242#create animation:
243SC.visualizationSettings.window.renderWindowSize=[1920,1080]
244SC.visualizationSettings.openGL.multiSampling = 4
245
246if useGraphics:
247 SC.visualizationSettings.general.autoFitScene=False
248
249 SC.renderer.Start()
250 if 'renderState' in exu.sys: SC.renderer.SetState(exu.sys['renderState']) #load last model view
251
252 SC.renderer.DoIdleTasks() #press space to continue
253
254mbs.SolveDynamic(simulationSettings=simulationSettings)
255
256uTip = mbs.GetSensorValues(sensTipDispl)
257exu.Print("nModes=", nModes, ", tip displacement=", uTip)
258
259
260if useGraphics:
261 SC.renderer.DoIdleTasks()
262 SC.renderer.Stop() #safely close rendering window!
263
264result = np.linalg.norm(uTip)
265exu.Print('solution of NGsolveCMStest=',result)
266
267#%%++++++++++++++++++++++++++++++++++++++++++++++++++++
268
269exudynTestGlobals.testError = result - (0.06953224923173523 )
270exudynTestGlobals.testResult = result
271
272
273#mbs.SolutionViewer()