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