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