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