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