objectFFRFTest.py
You can view and download this file on Github: objectFFRFTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test for ObjectFFRF with C++ implementation user function for reduced order equations of motion
5# NOTE: this is a development file, with lots of unstructured code; just kept for consistency!
6#
7# Author: Johannes Gerstmayr
8# Date: 2020-05-13
9#
10# 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.
11#
12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
13
14import exudyn as exu
15from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
16import exudyn.graphics as graphics
17from exudyn.FEM import *
18
19testIsActive = exu.sys.get('testIsActive', False)
20
21SC = exu.SystemContainer()
22mbs = SC.AddSystem()
23
24import numpy as np
25
26#==> DO NOT USE THE APPROACH of this test, but use the FEMinterface; also use the more efficient ObjectFFRFreducedOrder !!!!
27
28nodeDrawSize = 0.0025
29
30testMode = 1 #0=MarkerGeneric, 1=MarkerSuperElement
31modeNames = ['FFRF_MG','FFRF']
32
33#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34#Use FEMinterface to import FEM model and create FFRFreducedOrder object
35fem = FEMinterface()
36inputFileName = 'testData/rotorDiscTest' #runTestSuite.py is at another directory
37
38nodes=fem.ImportFromAbaqusInputFile(inputFileName+'.inp', typeName='Instance', name='rotor-1')
39elements = np.array(fem.elements[0]['Hex8'])
40
41fem.ReadMassMatrixFromAbaqus(inputFileName+'MASS1.mtx')
42fem.ReadStiffnessMatrixFromAbaqus(inputFileName+'STIF1.mtx')
43fem.ScaleStiffnessMatrix(1e-2) #for larger deformations, stiffness is reduced to 1%
44
45massMatrix = fem.GetMassMatrix(sparse=False)
46stiffnessMatrix = fem.GetStiffnessMatrix(sparse=False)
47
48#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
49
50exu.Print("nodes size=", nodes.shape)
51exu.Print("elements size=", elements.shape)
52
53minZ = min(nodes[:,2])
54maxZ = max(nodes[:,2])
55midZ = 0.5*(minZ+maxZ)
56
57#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
58
59#nLeft = (78-1)
60#nRight = (77-1)
61#nMid = (43-1)
62nLeft = -1
63nRight = -1
64nMid = -1
65nForce = -1 #40; node where fore is attached
66forceZ = 0.1
67
68#nTopRight = 12 #JG, excitation
69nForce = 12 #JG, excitation
70nTopMid = 9 #alternative: 103; JG, fixed node "rotation"
71
72unbalance = 0.1
73massMatrix[nTopMid*3+0,nTopMid*3+0] += unbalance
74massMatrix[nTopMid*3+1,nTopMid*3+1] += unbalance
75massMatrix[nTopMid*3+2,nTopMid*3+2] += unbalance
76
77#find output nodes:
78for i in range(len(nodes)):
79 n = nodes[i]
80 if abs(n[2] - minZ) < 1e-6 and abs(n[1]) < 1e-6 and abs(n[0]) < 1e-6:
81 nLeft = i
82 if abs(n[2] - maxZ) < 1e-6 and abs(n[1]) < 1e-6 and abs(n[0]) < 1e-6:
83 nRight = i
84 if abs(n[2] - midZ) < 1e-6 and abs(n[1]) < 1e-6 and abs(n[0]) < 1e-6:
85 nMid = i
86 #if abs(n[2] - forceZ) < 1e-6 and abs(n[1]) < 1e-6 and abs(n[0]) < 1e-6:
87 # nForce = i
88
89
90#exu.Print("nLeft=", nLeft, ", nRight=", nRight, ", nMid=", nMid, ", nForce=", nForce)
91
92#posX=0.15 #+/- x coordinate of nodes
93posLeft = nodes[nLeft]
94posRight = nodes[nRight]
95
96nNodes = len(nodes)
97nODE2 = nNodes*3
98exu.Print("nNodes=", nNodes, ", nODE2=", nODE2)
99
100#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
101calcEig = True
102if calcEig:
103 from scipy.linalg import solve, eigh, eig #eigh for symmetric matrices, positive definite
104
105 [eigvals, eigvecs] = eigh(stiffnessMatrix, massMatrix) #this gives omega^2 ... squared eigen frequencies (rad/s)
106 ev = np.sort(a=abs(eigvals))
107
108 listEig = []
109 for i in range(18):
110 listEig += [np.sqrt(ev[i])/(2*np.pi)]
111 exu.Print("eigenvalues =", listEig)
112#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
113
114#eigenvalues of constrained system:
115#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
116calcEigConstrained = False
117if calcEigConstrained:
118
119 constrainedCoordinates = [0]*nODE2
120 constrainedCoordinates[nLeft*3+0] = 1 #X
121 constrainedCoordinates[nLeft*3+1] = 1 #Y
122 constrainedCoordinates[nLeft*3+2] = 1 #Z
123 constrainedCoordinates[nRight*3+1] = 1 #Y
124 constrainedCoordinates[nRight*3+2] = 1 #Z
125
126 nConstrained = sum(constrainedCoordinates)
127 indexList = []
128 cnt = 0
129 for i in range(nODE2):
130 if constrainedCoordinates[i] == 0:
131 indexList+=[i]
132
133 nODE2C = nODE2-nConstrained
134 massMatrixC = np.zeros((nODE2C,nODE2C))
135 stiffnessMatrixC = np.zeros((nODE2C,nODE2C))
136
137 for i in range(nODE2C):
138 for j in range(nODE2C):
139 massMatrixC[i,j] = massMatrix[indexList[i],indexList[j]]
140 stiffnessMatrixC[i,j] = stiffnessMatrix[indexList[i],indexList[j]]
141
142 from scipy.linalg import solve, eigh, eig #eigh for symmetric matrices, positive definite
143
144 [eigvals, eigvecs] = eigh(stiffnessMatrixC, massMatrixC) #this gives omega^2 ... squared eigen frequencies (rad/s)
145 ev = np.sort(a=abs(eigvals))
146
147 listEig = []
148 for i in range(18):
149 listEig += [np.sqrt(ev[i])/(2*np.pi)]
150 exu.Print("eigenvalues of constrained system (Hz)=", listEig)
151
152
153#compute (3 x 3*n) skew matrix from (3*n) vector
154def ComputeSkewMatrix(v):
155 n = int(len(v)/3) #number of nodes
156 sm = np.zeros((3*n,3))
157
158 for i in range(n):
159 off = 3*i
160 x=v[off+0]
161 y=v[off+1]
162 z=v[off+2]
163 mLoc = np.array([[0,-z,y],[z,0,-x],[-y,x,0]])
164 sm[off:off+3,:] = mLoc[:,:]
165
166 return sm
167 #Y0=np.array([[0,0,0],[0,0,1],[0,-1,0]])
168 #Y1=np.array([[0,0,-1],[0,0,0],[1,0,0]])
169 #Y2=np.array([[0,1,0],[-1,0,0],[0,0,0]])
170
171
172#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
173
174nNodesFFRF = nNodes
175nODE2FFRF = nNodes*3
176
177#the following is only working for useFFRFobject = True; with False, it represents an old mode, deactivated with newer ObjectGenericODE2!
178useFFRFobject = True #uses ObjectFFRF instead of ObjectGenericODE2 ==> mesh nodes are indexed from 0 .. n_meshNodes-1
179decFFRFobject = 0 #adapt node numbers if useFFRFobject=True
180if useFFRFobject: decFFRFobject = 1
181
182useFFRF = True
183if useFFRF:
184 p0 = [0,0,midZ*0] #reference position
185 v0 = [0,0,0] #initial translational velocity
186 omega0 = [0,0,50*2*pi] #arbitrary initial angular velocity
187 ep0 = np.array(eulerParameters0) #no rotation
188 ep_t0 = AngularVelocity2EulerParameters_t(omega0, ep0)
189 #adjust mass and stiffness matrices
190 nODE2rigid = len(p0)+len(ep0)
191 nODE2rot = len(ep0) #dimension of rotation parameters
192 dim3D = len(p0) #dimension of position
193 nODE2FFRF = nODE2rigid + nODE2
194 nNodesFFRF = nNodes+1
195
196 Knew = np.zeros((nODE2FFRF,nODE2FFRF))
197 Mnew = np.zeros((nODE2FFRF,nODE2FFRF))
198
199 FillInSubMatrix(stiffnessMatrix, Knew, nODE2rigid, nODE2rigid)
200 FillInSubMatrix(massMatrix, Mnew, nODE2rigid, nODE2rigid)
201
202 Dnew = 2e-4*Knew #add little bit of damping
203 fNew = np.zeros(nODE2FFRF)
204
205 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++
206 #FFRF constant matrices:
207 unit3D = np.eye(3)
208 Phit = np.kron(np.ones(nNodes),unit3D).T
209 PhitTM = Phit.T @ massMatrix
210
211 Mtt = Phit.T @ massMatrix @ Phit
212 Mnew[0:3,0:3] = Mtt
213 totalMass = Mtt[0,0] #must be diagonal matrix with mass in diagonal
214
215 xRef = nodes.flatten() #node reference values in single vector (can be added then to q[7:])
216 xRefTilde = ComputeSkewMatrix(xRef) #rfTilde without q
217
218 inertiaLocal = xRefTilde.T @ massMatrix @ xRefTilde
219 if False:
220 exu.Print("Phit=", Phit[0:6,:])
221 exu.Print("PhitTM=", PhitTM[0:3,0:6])
222 exu.Print("xRef=", xRef[0:6])
223 exu.Print("xRefTilde=", xRefTilde[0:6,:])
224
225 exu.Print("python inertiaLocal=", inertiaLocal)
226 exu.Print("python totalMass=", totalMass)
227 exu.Print("python Mtt=", Mtt)
228
229 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++
230 #compute gravity term
231 g=np.array([0,-9.81,0]) #gravity vector
232 fGravRigid = list(totalMass*g)+[0,0,0,0]
233 #fGrav = np.array(fGravRigid + list((massMatrix @ Phit) @ g) ) #only local vector, without rotation
234 #exu.Print("fGrav=",fGrav)
235 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++
236
237
238
239# mbs.Reset()
240#background
241rect = [-0.3,-0.1,0.3,0.1] #xmin,ymin,xmax,ymax
242background = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[rect[0],rect[1],0, rect[2],rect[1],0, rect[2],rect[3],0, rect[0],rect[3],0, rect[0],rect[1],0]} #background
243oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
244mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=p0))
245#exu.Print("goundMarker=", mGround)
246
247nodeList = []
248nRB = -1
249
250if useFFRF:
251 nRB = mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=p0+list(ep0),
252 initialVelocities=v0+list(ep_t0)))
253 nodeList += [nRB]
254
255
256 #adjust node numbers:
257 #in all cases, triglist is same; elements = elements + 1 - decFFRFobject #increase node numbers, because of FFRFnode
258
259 #boundary nodes not adjusted for old constraints in
260 nLeft += 1
261 nRight += 1
262 nMid += 1
263 nForce += 1
264 nTopMid += 1
265
266for node in nodes:
267 n3 = mbs.AddNode(Point(referenceCoordinates = list(node), visualization=VNodePoint(show = not useFFRF))) #not useFFRF)))
268 nodeList += [n3]
269#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
270
271
272
273
274#exu.Print("nForce=", nForce)
275
276#conventional user function:
277def UFforce(mbs, t, itemIndex, q, q_t):
278 force = np.zeros(nODE2FFRF)
279 Avec = mbs.GetNodeOutput(nRB, exu.OutputVariableType.RotationMatrix)
280 A = Avec.reshape((3,3))
281
282 #implementation for Euler Parameters (Glocal_t*theta_t=0)
283 ep = np.array(q[dim3D:nODE2rigid]) + ep0 #add reference values, q are only the change w.r.t. reference values!
284 G = EulerParameters2GLocal(ep)
285
286 cF_t = np.array(q_t[nODE2rigid:]) #velocities of flexible coordinates
287
288 rF = xRef + np.array(q[nODE2rigid:]) #nodal position
289
290 omega3D = G @ np.array(q_t[dim3D:nODE2rigid])
291 omega3Dtilde = Skew(omega3D)
292 omega = np.array(list(omega3D)*nNodes)
293 omegaTilde = np.kron(np.eye(nNodes),omega3Dtilde)
294
295 #squared angul. vel. matrix:
296 omega3Dtilde2 = Skew(omega3D) @ Skew(omega3D)
297 omegaTilde2 = np.kron(np.eye(nNodes),omega3Dtilde2)
298
299 if True: #for rotordynamics, we assume rigid body motion with constant ang. vel.
300 #these 2 terms are computationally costly:
301 rfTilde = ComputeSkewMatrix(rF) #rfTilde
302 cF_tTilde = ComputeSkewMatrix(cF_t)
303
304 fTrans = A @ (omega3Dtilde @ PhitTM @ rfTilde @ omega3D + 2*PhitTM @ cF_tTilde @ omega3D)
305 force[0:dim3D] = fTrans
306
307 fRot = -G.T@(omega3Dtilde @ rfTilde.T @ massMatrix @ rfTilde @ omega3D +
308 2*rfTilde.T @ massMatrix @ cF_tTilde @ omega3D)
309 force[dim3D:nODE2rigid] = fRot
310
311 fFlex = -massMatrix @ (omegaTilde2 @ rF + 2*(omegaTilde @ cF_t))
312 force[nODE2rigid:] = fFlex
313
314 #add gravity:
315 if False:
316 fGrav = np.array(fGravRigid + list(PhitTM.T @ (A.T @ g)) ) #only local vector, without rotation
317 force += fGrav
318
319
320 return force
321
322#ffrf mass matrix:
323def UFmassGenericODE2(mbs, t, itemIndex, q, q_t):
324 Avec = mbs.GetNodeOutput(nRB, exu.OutputVariableType.RotationMatrix)
325 A = Avec.reshape((3,3))
326 ep = q[dim3D:nODE2rigid] + ep0 #add reference values, q are only the change w.r.t. reference values!
327 G = EulerParameters2GLocal(ep)
328
329 rF = xRef + q[nODE2rigid:] #nodal position
330 rfTilde = ComputeSkewMatrix(rF) #rfTilde
331
332 #Mtr:
333 Mtr = -A @ PhitTM @ rfTilde @ G
334 Mnew[0:dim3D, dim3D:dim3D+nODE2rot] = Mtr
335 Mnew[dim3D:dim3D+nODE2rot, 0:dim3D] = Mtr.T
336 #Mtf:
337 Mtf = A @ PhitTM
338 Mnew[0:dim3D, nODE2rigid:] = Mtf
339 Mnew[nODE2rigid:, 0:dim3D] = Mtf.T
340 #Mrf:
341 Mrf = -G.T @ rfTilde.T @ massMatrix
342 Mnew[dim3D:dim3D+nODE2rot, nODE2rigid:] = Mrf
343 Mnew[nODE2rigid:, dim3D:dim3D+nODE2rot] = Mrf.T
344 #Mrr:
345 Mnew[dim3D:dim3D+nODE2rot, dim3D:dim3D+nODE2rot] = -Mrf @ rfTilde @ G #G.T @ rfTilde.T @ massMatrix @ rfTilde @ G
346
347 #exu.Print(np.linalg.norm(rF))
348 #omega3D = G @ q_t[dim3D:nODE2rigid]
349 #exu.Print(omega3D)
350
351 #exu.Print("Mtt Mtr Mtf=",Mnew[0:3,0:10].round(5))
352 #exu.Print("Mrr=",Mnew[3:7,3:7].round(5))
353 #exu.Print("Mff=",Mnew[7:10,7:13].round(5))
354 #Mnew[:,:] = 0 #for testing
355 return Mnew
356
357
358#convert elements to triangles for drawing:
359trigList = []
360for element in elements:
361 trigList += ConvertHexToTrigs(element)
362trigList = np.array(trigList)
363#exu.Print("trig list=", trigList)
364#exu.Print("trig list size=", trigList.shape)
365
366stiffnessMatrixFF = exu.MatrixContainer()
367stiffnessMatrixFF.SetWithDenseMatrix(stiffnessMatrix,useDenseMatrix=False)
368massMatrixFF = exu.MatrixContainer()
369massMatrixFF.SetWithDenseMatrix(massMatrix,useDenseMatrix=False)
370emptyMC = exu.MatrixContainer()
371
372#add generic body for FFRF-Object:
373if useFFRFobject:
374 oGenericODE2 = mbs.AddObject(ObjectFFRF(nodeNumbers = nodeList,
375 #massMatrixFF=Mnew,
376 stiffnessMatrixFF=stiffnessMatrixFF,
377 #dampingMatrixFF=Dnew,
378 massMatrixFF=massMatrixFF,
379 #dampingMatrixFF=emptyMC,
380 #forceVector=fNew, #now is a global vector!
381 #forceUserFunction=UFforce,
382 #computeFFRFterms=True,
383 #massMatrixUserFunction=UFmassGenericODE2,
384 visualization=VObjectFFRF(triangleMesh = trigList,
385 color=graphics.color.lightred,
386 showNodes = True)))
387else:
388 oGenericODE2 = mbs.AddObject(ObjectGenericODE2(nodeNumbers = nodeList,
389 massMatrix=Mnew,
390 stiffnessMatrix=Knew,
391 dampingMatrix=Dnew,
392 forceVector=fNew, forceUserFunction=UFforce,
393 useFirstNodeAsReferenceFrame=True, #does not exist anymore
394 massMatrixUserFunction=UFmassGenericODE2,
395 visualization=VObjectGenericODE2(triangleMesh = trigList,
396 color=graphics.color.lightred,
397 showNodes = True)))
398
399if nODE2rot == 4: #for euler parameters --> add body to constrain EP
400 epsMass = 1e-3#needed, if not all ffrf terms are included
401 #add rigid body to node for Euler Parameter constraint:
402 nReferenceFrame = mbs.AddObject(ObjectRigidBody(nodeNumber=nRB, mass=epsMass, inertia=[epsMass,epsMass,epsMass,0,0,0]))
403
404mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nRB))
405#exu.Print("rigidNodeMarker=", mRB)
406
407
408#mbs.AddLoad(Torque(markerNumber=mRB, loadVector=[0,0,100*2*pi])) #add drive for reference frame
409
410if False: #OPTIONAL: lock rigid body motion of reference frame (for tests):
411 mbs.AddObject(GenericJoint(markerNumbers=[mGround, mRB], constrainedAxes=[1,1,1, 1,1,0]))
412
413
414
415#ground point:
416nGroundLeft = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,minZ], visualization = VNodePointGround(show=False))) #ground node for coordinate constraint
417nGroundRight = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,maxZ], visualization = VNodePointGround(show=False))) #ground node for coordinate constraint
418
419mGroundLeft = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGroundLeft, coordinate=0)) #Ground node ==> no action
420mGroundRight = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGroundRight, coordinate=0)) #Ground node ==> no action
421
422mGroundPosLeft = mbs.AddMarker(MarkerNodePosition(nodeNumber = nGroundLeft)) #Ground node ==> no action
423mGroundPosRight = mbs.AddMarker(MarkerNodePosition(nodeNumber = nGroundRight)) #Ground node ==> no action
424
425#exu.Print("ground Node/Coordinate Markers =", mGroundLeft, "...", mGroundPosRight)
426
427#++++++++++++++++++++++++++++++++++++++++++
428#find nodes at left and right surface:
429nodeListLeft = []
430nodeListRight = []
431
432for i in range(len(nodes)):
433 n = nodes[i]
434 if abs(n[2] - minZ) < 1e-6:
435 nodeListLeft += [i+useFFRF] #add 1 for rigid body node, which is first node in GenericODE2 object
436 elif abs(n[2] - maxZ) < 1e-6:
437 nodeListRight += [i+useFFRF]
438
439#exu.Print("nodeListLeft =",nodeListLeft)
440#exu.Print("nodeListRight =",nodeListRight)
441
442lenLeft = len(nodeListLeft)
443lenRight = len(nodeListRight)
444weightsLeft = np.array((1./lenLeft)*np.ones(lenLeft))
445weightsRight = np.array((1./lenRight)*np.ones(lenRight))
446
447#exu.Print("nodeLeft =",nLeft)
448#exu.Print("nodeRight =",nRight)
449
450#lock FFRF reference frame:
451for i in range(3):
452 mLeft = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nLeft, coordinate=i))
453# exu.Print("mLeftCoord=", mLeft)
454
455 mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundLeft,mLeft]))
456 if i != 2: #exclude double constraint in z-direction (axis)
457 mRight = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRight, coordinate=i))
458# exu.Print("mRightCoord=", mRight)
459 mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundRight,mRight]))
460
461#lock rotation (also needed in FFRF):
462mTopRight = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nTopMid, coordinate=0)) #x-coordinate of node with y-max
463#exu.Print("mTopRight=", mTopRight)
464mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundRight,mTopRight]))
465
466addSupports = True
467if addSupports:
468 k = 2e8
469 d = k*0.01
470
471 useSpringDamper = True
472
473 if testMode == 0:
474 raise ValueError('does not exist any more')
475 # mLeft = mbs.AddMarker(MarkerGenericBodyPosition(bodyNumber=oGenericODE2,
476 # nodeNumbers=nodeListLeft,
477 # weightingFactors=weightsLeft,
478 # useFirstNodeAsReferenceFrame=useFFRF))
479 # mRight = mbs.AddMarker(MarkerGenericBodyPosition(bodyNumber=oGenericODE2,
480 # nodeNumbers=nodeListRight,
481 # weightingFactors=weightsRight,
482 # useFirstNodeAsReferenceFrame=useFFRF))
483 else:
484 mLeft = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=oGenericODE2,
485 meshNodeNumbers=np.array(nodeListLeft)-1, #these are the meshNodeNumbers
486 weightingFactors=weightsLeft))
487 mRight = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=oGenericODE2,
488 meshNodeNumbers=np.array(nodeListRight)-1, #these are the meshNodeNumbers
489 weightingFactors=weightsRight))
490
491 if useSpringDamper:
492 oSJleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mLeft, mGroundPosLeft],
493 stiffness=[k,k,k], damping=[d,d,d]))
494 oSJright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mRight,mGroundPosRight],
495 stiffness=[k,k,0], damping=[d,d,d]))
496 else:
497 oSJleft = mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosLeft,mLeft], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
498 oSJright= mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosRight,mRight], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
499
500
501fileDir = 'solution/'
502sDisp=mbs.AddSensor(SensorSuperElement(bodyNumber=oGenericODE2, meshNodeNumber=nMid-1, #meshnode is -1
503 storeInternal=True,#fileName=fileDir+'nMidDisplacement'+modeNames[testMode]+'test.txt',
504 outputVariableType = exu.OutputVariableType.Displacement))
505
506
507#exu.Print(mbs)
508mbs.Assemble()
509
510#exu.Print("ltg GenericODE2 left =", mbs.systemData.GetObjectLTGODE2(oSJleft))
511#exu.Print("ltg GenericODE2 right=", mbs.systemData.GetObjectLTGODE2(oSJright))
512
513simulationSettings = exu.SimulationSettings()
514
515SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
516SC.visualizationSettings.nodes.drawNodesAsPoint = False
517SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
518
519SC.visualizationSettings.nodes.show = True
520SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
521SC.visualizationSettings.nodes.basisSize = 0.12
522SC.visualizationSettings.bodies.deformationScaleFactor = 10
523
524SC.visualizationSettings.view0.scene.showFaceEdges = True
525SC.visualizationSettings.view0.scene.showFaces = True
526
527SC.visualizationSettings.sensors.show = True
528SC.visualizationSettings.sensors.drawSimplified = False
529SC.visualizationSettings.sensors.defaultSize = 0.01
530SC.visualizationSettings.markers.drawSimplified = False
531SC.visualizationSettings.markers.show = True
532SC.visualizationSettings.markers.defaultSize = 0.01
533
534SC.visualizationSettings.loads.drawSimplified = False
535
536SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
537SC.visualizationSettings.contour.outputVariableComponent = 2 #z-component
538
539simulationSettings.solution.file.information = modeNames[testMode]
540simulationSettings.solution.file.write=False
541
542h=1e-4
543tEnd = 0.001
544simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
545simulationSettings.timeIntegration.endTime = tEnd
546simulationSettings.solution.file.writePeriod = h
547simulationSettings.timeIntegration.verboseMode = 1
548simulationSettings.timeIntegration.newton.useModifiedNewton = True
549#simulationSettings.timeIntegration.newton.maxModifiedNewtonIterations = 10
550#simulationSettings.timeIntegration.newton.modifiedNewtonJacUpdatePerStep = True #this improves the FFRF simulation slightly
551
552simulationSettings.solution.sensors.writePeriod = h
553
554simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5 #SHOULD work with 0.9 as well
555#simulationSettings.show.statistics = True
556#simulationSettings.show.computationTime = True
557
558#create animation:
559#simulationSettings.solution.recordImagesInterval = 0.0002
560#SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
561
562if not testIsActive:
563 SC.renderer.Start()
564 if 'lastRenderState' in vars():
565 SC.renderer.SetState(lastRenderState) #load last model view
566
567 SC.renderer.DoIdleTasks() #press space to continue
568
569mbs.SolveDynamic(simulationSettings)
570
571data = mbs.GetSensorStoredData(sDisp)
572#data = np.loadtxt(fileDir+'nMidDisplacement'+modeNames[testMode]+'test.txt', comments='#', delimiter=',')
573result = abs(data).sum()
574#pos = mbs.GetObjectOutputBody(objFFRF['oFFRFreducedOrder'],exu.OutputVariableType.Position, localPosition=[0,0,0])
575exu.Print('solution of ObjectFFRF=',result)
576
577exu.sys['testResult'] = result#0.006460010812070858
578
579
580if not testIsActive:
581 SC.renderer.DoIdleTasks()
582 SC.renderer.Stop() #safely close rendering window!
583 lastRenderState = SC.renderer.GetState() #store model view for next simulation
584
585##++++++++++++++++++++++++++++++++++++++++++++++q+++++++
586#plot results
587cList=['r-','g-','b-','k-','c-','r:','g:','b:','k:','c:']
588if not testIsActive:
589
590
591 mbs.PlotSensor(sDisp, components=0, closeAll=True)