NGsolvePostProcessingStresses.py
You can view and download this file on Github: NGsolvePostProcessingStresses.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test for meshing with NETGEN and import of FEM model;
5# Model is a simple flexible pendulum meshed with tet elements;
6# Note that the model is overly flexible (linearized strain assumption not valid),
7# but it should serve as a demonstration of the FFRFreducedOrder modeling
8#
9# Author: Johannes Gerstmayr
10# Date: 2021-02-05
11#
12# 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.
13#
14#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15
16
17import exudyn as exu
18from exudyn.itemInterface import *
19from exudyn.utilities import *
20from exudyn.FEM import *
21from exudyn.graphicsDataUtilities import *
22
23SC = exu.SystemContainer()
24mbs = SC.AddSystem()
25
26import numpy as np
27
28import timeit
29
30import exudyn.basicUtilities as eb
31import exudyn.rigidBodyUtilities as rb
32import exudyn.utilities as eu
33
34import numpy as np
35
36useGraphics = True
37fileName = 'testData/netgenBrick' #for load/save of FEM data
38
39
40if __name__ == '__main__': #needed to use multiprocessing for mode computation
41
42
43 #+++++++++++++++++++++++++++++++++++++++++++++++++++++
44 #netgen/meshing part:
45 fem = FEMinterface()
46 #standard:
47 a = 0.025 #height/width of beam
48 b = a
49 h = 0.3*a
50 L = 1 #Length of beam
51 nModes = 10
52
53 #plate:
54 # a = 0.025 #height/width of beam
55 # b = 0.4
56 # L = 1 #Length of beam
57 # h = 0.6*a
58 # nModes = 40
59
60 rho = 1000
61 Emodulus=1e7
62 nu=0.3
63 meshCreated = False
64
65 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
66 if True: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
67
68 import ngsolve as ngs
69 from netgen.geom2d import unit_square
70 import netgen.libngpy as libng
71 from netgen.csg import *
72
73 geo = CSGeometry()
74
75 block = OrthoBrick(Pnt(0,-a,-b),Pnt(L,a,b))
76 geo.Add(block)
77
78 #Draw (geo)
79
80 mesh = ngs.Mesh( geo.GenerateMesh(maxh=h))
81 mesh.Curve(1)
82
83 if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
84 import netgen.gui
85 Draw (mesh)
86 netgen.Redraw()
87
88 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
89 #Use fem to import FEM model and create FFRFreducedOrder object
90 fem.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
91 meshCreated = True
92 if (h==a): #save only if it has smaller size
93 fem.SaveToFile(fileName)
94
95 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
96 if not meshCreated: fem.LoadFromFile(fileName)
97
98 print("nNodes=",fem.NumberOfNodes())
99 fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
100 #print("eigen freq.=", fem.GetEigenFrequenciesHz())
101
102 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
103 #compute stress modes:
104 mat = KirchhoffMaterial(Emodulus, nu, rho)
105 varType = exu.OutputVariableType.StressLocal
106 #varType = exu.OutputVariableType.StrainLocal
107 print("ComputePostProcessingModes ... (may take a while)")
108 start_time = time.time()
109 fem.ComputePostProcessingModes(material=mat,
110 outputVariableType=varType,
111 numberOfThreads=5)
112 print("--- %s seconds ---" % (time.time() - start_time))
113
114 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
115 print("create CMS element ...")
116 cms = ObjectFFRFreducedOrderInterface(fem)
117
118 objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
119 initialVelocity=[0,0,0], initialAngularVelocity=[0,0,0],
120 color=[0.1,0.9,0.1,1.])
121
122 # mbs.SetObjectParameter(objectNumber=objFFRF['oFFRFreducedOrder'],
123 # parameterName='outputVariableModeBasis',
124 # value=stressModes)
125
126 # mbs.SetObjectParameter(objectNumber=objFFRF['oFFRFreducedOrder'],
127 # parameterName='outputVariableTypeModeBasis',
128 # value=exu.OutputVariableType.StressLocal) #type=stress modes ...
129
130
131 #add gravity (not necessary if user functions used)
132 oFFRF = objFFRF['oFFRFreducedOrder']
133 mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber=oFFRF))
134 mbs.AddLoad(LoadMassProportional(markerNumber=mBody, loadVector= [0,-9.81,0]))
135
136 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
137 #add markers and joints
138 nodeDrawSize = 0.0025 #for joint drawing
139
140 pLeft = [0,-a,-b]
141 pRight = [0,-a,b]
142 nTip = fem.GetNodeAtPoint([L,-a,-b]) #tip node
143 #print("nMid=",nMid)
144
145 mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
146 oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
147
148 mGroundPosLeft = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pLeft))
149 mGroundPosRight = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pRight))
150
151 #++++++++++++++++++++++++++++++++++++++++++
152 #find nodes at left and right surface:
153 #nodeListLeft = fem.GetNodesInPlane(pLeft, [0,0,1])
154 #nodeListRight = fem.GetNodesInPlane(pRight, [0,0,1])
155 nodeListLeft = [fem.GetNodeAtPoint(pLeft)]
156 nodeListRight = [fem.GetNodeAtPoint(pRight)]
157
158
159 lenLeft = len(nodeListLeft)
160 lenRight = len(nodeListRight)
161 weightsLeft = np.array((1./lenLeft)*np.ones(lenLeft))
162 weightsRight = np.array((1./lenRight)*np.ones(lenRight))
163
164 addSupports = True
165 if addSupports:
166 k = 10e8 #joint stiffness
167 d = k*0.01 #joint damping
168
169 useSpringDamper = True
170
171 mLeft = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=objFFRF['oFFRFreducedOrder'],
172 meshNodeNumbers=np.array(nodeListLeft), #these are the meshNodeNumbers
173 weightingFactors=weightsLeft))
174 mRight = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=objFFRF['oFFRFreducedOrder'],
175 meshNodeNumbers=np.array(nodeListRight), #these are the meshNodeNumbers
176 weightingFactors=weightsRight))
177 if useSpringDamper:
178 oSJleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mLeft, mGroundPosLeft],
179 stiffness=[k,k,k], damping=[d,d,d]))
180 oSJright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mRight,mGroundPosRight],
181 stiffness=[k,k,0], damping=[d,d,d]))
182 else:
183 oSJleft = mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosLeft,mLeft], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
184 oSJright= mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosRight,mRight], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
185
186
187 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
188 fileDir = 'solution/'
189 mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'], meshNodeNumber=nTip, #meshnode number!
190 fileName=fileDir+'nMidDisplacementCMS'+str(nModes)+'Test.txt',
191 outputVariableType = exu.OutputVariableType.Displacement))
192
193 mbs.Assemble()
194
195 simulationSettings = exu.SimulationSettings()
196
197 SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
198 SC.visualizationSettings.nodes.drawNodesAsPoint = False
199 SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
200
201 SC.visualizationSettings.nodes.show = False
202 SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
203 SC.visualizationSettings.nodes.basisSize = 0.12
204 SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
205
206 SC.visualizationSettings.openGL.showFaceEdges = True
207 SC.visualizationSettings.openGL.showFaces = True
208
209 SC.visualizationSettings.sensors.show = True
210 SC.visualizationSettings.sensors.drawSimplified = False
211 SC.visualizationSettings.sensors.defaultSize = 0.01
212 SC.visualizationSettings.markers.drawSimplified = False
213 SC.visualizationSettings.markers.show = False
214 SC.visualizationSettings.markers.defaultSize = 0.01
215
216 SC.visualizationSettings.loads.drawSimplified = False
217
218 # SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
219 # SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
220 SC.visualizationSettings.contour.reduceRange=False
221 SC.visualizationSettings.contour.outputVariable = varType
222 SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
223
224 simulationSettings.solutionSettings.solutionInformation = "ObjectFFRFreducedOrder test"
225
226 h=0.25e-3
227 tEnd = 0.05
228
229 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
230 simulationSettings.timeIntegration.endTime = tEnd
231 simulationSettings.solutionSettings.writeSolutionToFile = False
232 simulationSettings.timeIntegration.verboseMode = 1
233 #simulationSettings.timeIntegration.verboseModeFile = 3
234 simulationSettings.timeIntegration.newton.useModifiedNewton = True
235
236 simulationSettings.solutionSettings.sensorsWritePeriod = h
237
238 simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8 #SHOULD work with 0.9 as well
239 #simulationSettings.displayStatistics = True
240 #simulationSettings.displayComputationTime = True
241
242 #create animation:
243 # simulationSettings.solutionSettings.recordImagesInterval = 0.005
244 # SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
245 SC.visualizationSettings.window.renderWindowSize=[1920,1080]
246 SC.visualizationSettings.openGL.multiSampling = 4
247
248 if True:
249 if useGraphics:
250 SC.visualizationSettings.general.autoFitScene=False
251
252 exu.StartRenderer()
253 if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
254
255 mbs.WaitForUserToContinue() #press space to continue
256
257 mbs.SolveDynamic(solverType=exu.DynamicSolverType.TrapezoidalIndex2,
258 simulationSettings=simulationSettings)
259
260
261
262 if useGraphics:
263 SC.WaitForRenderEngineStopFlag()
264 exu.StopRenderer() #safely close rendering window!