objectFFRFreducedOrderNetgen.py

You can view and download this file on Github: objectFFRFreducedOrderNetgen.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
 16import exudyn as exu
 17from exudyn.itemInterface import *
 18from exudyn.utilities import *
 19from exudyn.FEM import *
 20from exudyn.graphicsDataUtilities import *
 21
 22SC = exu.SystemContainer()
 23mbs = SC.AddSystem()
 24
 25import numpy as np
 26
 27import timeit
 28
 29import exudyn.basicUtilities as eb
 30import exudyn.rigidBodyUtilities as rb
 31import exudyn.utilities as eu
 32
 33import numpy as np
 34
 35useGraphics = True
 36fileName = 'testData/netgenBrick' #for load/save of FEM data
 37#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 38#netgen/meshing part:
 39femInterface = FEMinterface()
 40a = 0.025 #height/width of beam
 41L = 1     #Length of beam
 42
 43if False: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
 44
 45    from ngsolve import *
 46    from netgen.geom2d import unit_square
 47    import netgen.libngpy as libng
 48    from netgen.csg import *
 49
 50    geo = CSGeometry()
 51
 52    block = OrthoBrick(Pnt(0,-a,-a),Pnt(L,a,a))
 53    geo.Add(block)
 54
 55    #Draw (geo)
 56
 57    mesh = Mesh( geo.GenerateMesh(maxh=a))
 58    mesh.Curve(1)
 59
 60    if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
 61        import netgen.gui
 62        Draw (mesh)
 63        netgen.Redraw()
 64
 65    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 66    #Use FEMinterface to import FEM model and create FFRFreducedOrder object
 67    femInterface.ImportMeshFromNGsolve(mesh, density=1000, youngsModulus=5e6, poissonsRatio=0.3)
 68    femInterface.SaveToFile(fileName)
 69
 70if True: #now import mesh as mechanical model to EXUDYN
 71    femInterface.LoadFromFile(fileName)
 72
 73    nModes = 12
 74    femInterface.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
 75    #print("eigen freq.=", femInterface.GetEigenFrequenciesHz())
 76
 77    cms = ObjectFFRFreducedOrderInterface(femInterface)
 78
 79    #user functions should be defined outside of class:
 80    def UFmassFFRFreducedOrder(mbs, t, itemIndex, qReduced, qReduced_t):
 81        return cms.UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
 82
 83    def UFforceFFRFreducedOrder(mbs, t, itemIndex, qReduced, qReduced_t):
 84        return cms.UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
 85
 86    objFFRF = cms.AddObjectFFRFreducedOrderWithUserFunctions(exu, mbs, positionRef=[0,0,0], eulerParametersRef=eulerParameters0,
 87                                                  initialVelocity=[0,0,0], initialAngularVelocity=[0,0,0],
 88                                                  gravity = [0,-9.81,0],
 89                                                  UFforce=UFforceFFRFreducedOrder, UFmassMatrix=UFmassFFRFreducedOrder,
 90                                                  color=[0.1,0.9,0.1,1.])
 91
 92    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 93    #add markers and joints
 94    nodeDrawSize = 0.0025 #for joint drawing
 95
 96    pLeft = [0,-a,-a]
 97    pRight = [0,-a,a]
 98    nTip = femInterface.GetNodeAtPoint([L,-a,-a]) #tip node
 99    #print("nMid=",nMid)
100
101    mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
102    oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
103
104    mGroundPosLeft = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pLeft))
105    mGroundPosRight = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pRight))
106
107    #++++++++++++++++++++++++++++++++++++++++++
108    #find nodes at left and right surface:
109    #nodeListLeft = femInterface.GetNodesInPlane(pLeft, [0,0,1])
110    #nodeListRight = femInterface.GetNodesInPlane(pRight, [0,0,1])
111    nodeListLeft = [femInterface.GetNodeAtPoint(pLeft)]
112    nodeListRight = [femInterface.GetNodeAtPoint(pRight)]
113
114
115    lenLeft = len(nodeListLeft)
116    lenRight = len(nodeListRight)
117    weightsLeft = np.array((1./lenLeft)*np.ones(lenLeft))
118    weightsRight = np.array((1./lenRight)*np.ones(lenRight))
119
120    addSupports = True
121    if addSupports:
122        k = 10e8     #joint stiffness
123        d = k*0.01  #joint damping
124
125        useSpringDamper = True
126
127        mLeft = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=objFFRF['oFFRFreducedOrder'],
128                                                        meshNodeNumbers=np.array(nodeListLeft), #these are the meshNodeNumbers
129                                                        weightingFactors=weightsLeft))
130        mRight = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=objFFRF['oFFRFreducedOrder'],
131                                                        meshNodeNumbers=np.array(nodeListRight), #these are the meshNodeNumbers
132                                                        weightingFactors=weightsRight))
133        if useSpringDamper:
134            oSJleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mLeft, mGroundPosLeft],
135                                                stiffness=[k,k,k], damping=[d,d,d]))
136            oSJright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mRight,mGroundPosRight],
137                                                stiffness=[k,k,0], damping=[d,d,d]))
138        else:
139            oSJleft = mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosLeft,mLeft], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
140            oSJright= mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosRight,mRight], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
141
142
143    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
144    fileDir = 'solution/'
145    mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'], meshNodeNumber=nTip, #meshnode number!
146                             fileName=fileDir+'nMidDisplacementCMS'+str(nModes)+'Test.txt',
147                             outputVariableType = exu.OutputVariableType.Displacement))
148
149    mbs.Assemble()
150
151    simulationSettings = exu.SimulationSettings()
152
153    SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
154    SC.visualizationSettings.nodes.drawNodesAsPoint = False
155    SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
156
157    SC.visualizationSettings.nodes.show = False
158    SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
159    SC.visualizationSettings.nodes.basisSize = 0.12
160    SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
161
162    SC.visualizationSettings.openGL.showFaceEdges = True
163    SC.visualizationSettings.openGL.showFaces = True
164
165    SC.visualizationSettings.sensors.show = True
166    SC.visualizationSettings.sensors.drawSimplified = False
167    SC.visualizationSettings.sensors.defaultSize = 0.01
168    SC.visualizationSettings.markers.drawSimplified = False
169    SC.visualizationSettings.markers.show = True
170    SC.visualizationSettings.markers.defaultSize = 0.01
171
172    SC.visualizationSettings.loads.drawSimplified = False
173
174    SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
175    SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
176
177    simulationSettings.solutionSettings.solutionInformation = "ObjectFFRFreducedOrder test"
178
179    h=5e-4
180    tEnd = 3
181
182    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
183    simulationSettings.timeIntegration.endTime = tEnd
184    simulationSettings.solutionSettings.solutionWritePeriod = h
185    simulationSettings.timeIntegration.verboseMode = 1
186    #simulationSettings.timeIntegration.verboseModeFile = 3
187    simulationSettings.timeIntegration.newton.useModifiedNewton = True
188
189    simulationSettings.solutionSettings.sensorsWritePeriod = h
190    simulationSettings.solutionSettings.coordinatesSolutionFileName = "solution/coordinatesSolutionCMStest.txt"
191
192    simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5 #SHOULD work with 0.9 as well
193    #simulationSettings.displayStatistics = True
194    #simulationSettings.displayComputationTime = True
195
196    #create animation:
197    #simulationSettings.solutionSettings.recordImagesInterval = 0.0002
198    #SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
199
200    if True:
201        if useGraphics:
202            exu.StartRenderer()
203            if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
204
205            mbs.WaitForUserToContinue() #press space to continue
206
207        mbs.SolveDynamic(solverType=exu.DynamicSolverType.TrapezoidalIndex2,
208                         simulationSettings=simulationSettings)
209
210
211
212        if useGraphics:
213            SC.WaitForRenderEngineStopFlag()
214            exu.StopRenderer() #safely close rendering window!
215            lastRenderState = SC.GetRenderState() #store model view for next simulation