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!