NGsolvePistonEngine.py
You can view and download this file on Github: NGsolvePistonEngine.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: generate a piston engine with finite element mesh
5# created with NGsolve and with variable number of pistons
6#
7# Author: Johannes Gerstmayr
8# Date: 2020-06-12
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
14
15import sys
16import exudyn as exu
17
18from exudyn.itemInterface import *
19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
20import exudyn.graphics as graphics #only import if it does not conflict
21from exudyn.rigidBodyUtilities import *
22from exudyn.FEM import *
23
24import time
25
26import mkl
27mkl.set_num_threads(20)
28
29from ngsolve import *
30from netgen.geom2d import unit_square
31
32import netgen.libngpy as libng
33
34netgenDrawing = False #set true, to show geometry and mesh in NETGEN
35#if netgenDrawing, uncomment the following line and execute in external terminal, not in spyder (see preferences "Run"):
36#import netgen.gui
37
38from netgen.csg import *
39
40import numpy as np
41import timeit
42
43verbose = True
44meshSize = 0.005*2*2 #fast: 0.005*2; standard:0.005; fine: 0.0011: memory limit (96GB) for NGsolve; < 0.0015 makes problems with scipy eigensolver
45meshOrder = 1 #2 for stresses!
46showStresses = True #may take very long for large number of modes/nodes
47
48#++++++++++++++++++++++++++++++++++++
49#helper functions (copied from EXUDYN):
50def RotationMatrixZ(angleRad):
51 return np.array([ [np.cos(angleRad),-np.sin(angleRad), 0],
52 [np.sin(angleRad), np.cos(angleRad), 0],
53 [0, 0, 1] ]);
54
55def VAdd(v0, v1):
56 if len(v0) != len(v1): print("ERROR in VAdd: incompatible vectors!")
57 n = len(v0)
58 v = [0]*n
59 for i in range(n):
60 v[i] = v0[i]+v1[i]
61 return v
62
63def VSub(v0, v1):
64 if len(v0) != len(v1): print("ERROR in VSub: incompatible vectors!")
65 n = len(v0)
66 v = [0]*n
67 for i in range(n):
68 v[i] = v0[i]-v1[i]
69 return v
70
71def NormL2(vector):
72 value = 0
73 for x in vector:
74 value += x**2
75 return value**0.5
76
77def Normalize(v):
78 v2=[0]*len(v)
79
80 fact = NormL2(v)
81 fact = 1./fact
82 for i in range(len(v2)):
83 v2[i]=fact*v[i]
84 return v2
85#++++++++++++++++++++++++++++++++++++
86startTotal = timeit.default_timer()
87#parameters
88
89#crank:
90b1 = 0.012 #width of journal bearing
91r1 = 0.012 #radius of journal bearing
92dk = 0.015 #crank arm width (z)
93bk = 0.032 #crank arm size (y)
94
95l3 = 0.030
96l4 = 0.040
97#l4x= 0.005 #offset of counterweight
98lk = 0.030 #l4*0.5+l3 #crank arm length (x)
99bm = 0.065
100dBevel = dk*0.5
101#shaft:
102r0 = 0.012 #0.012
103d0 = 0.020 #shaft length at left/right support
104d1 = 0.012 #shaft length at intermediate support
105
106#distance rings:
107db = 0.002 #width of distance ring
108rdb0 = r0+db #total radius of distance ring, shaft
109rdb1 = r1+db #total radius of distance ring, crank
110
111#conrod:
112bc = 0.024 #height of conrod
113dc = 0.012 #width of conrod
114lc = 0.080 #length of conrod (axis-axis)
115r1o= r1+0.006 #outer radius of conrod at crank joint
116r2 = 0.008 #radius of piston journal bearing
117r2o= r2+0.006 #outer radius of conrod at piston joint
118
119cylOffZ=0.010 #z-offset of cylinder cut out of conrod
120cylR = 0.008 #radius of cylinder cut out of conrod
121
122angC = 4*np.pi/180
123
124#piston:
125dpb = r2o-0.000 #axis inside piston
126r2p = r2o+0.004 #0.018
127lp = 0.034
128bp = 0.050
129lpAxis = dc+2*db
130lOffCut = 0.011 #offset for cutout of big cylinder
131
132#total length of one segment:
133lTotal = db+dk+db+b1+db+dk+db+d1
134
135#eps
136eps = 5e-4 #added to faces, to avoid CSG-problems
137
138#++++++++++++++++++++++++++++++++++++
139#points
140pLB = [0 ,0,-d0]
141p0B = [0 ,0,0]
142p1B = [0 ,0,db]
143#p2B = [0, 0,db+dk]
144p21B =[lk,0,db+dk]
145p31B = [lk,0,db+dk+db]
146p41B = [lk,0,db+dk+db+b1]
147p51B =[lk,0,db+dk+db+b1+db]
148p6B = [0 ,0,db+dk+db+b1+db+dk]
149p7B = [0 ,0,db+dk+db+b1+db+dk+db]
150p8B = [0 ,0,lTotal]
151
152def CSGcylinder(p0,p1,r):
153 v = VSub(p1,p0)
154 v = Normalize(v)
155 cyl = Cylinder(Pnt(p0[0],p0[1],p0[2]), Pnt(p1[0],p1[1],p1[2]),
156 r) * Plane(Pnt(p0[0],p0[1],p0[2]), Vec(-v[0],-v[1],-v[2])) * Plane(Pnt(p1[0],p1[1],p1[2]), Vec(v[0],v[1],v[2]))
157 return cyl
158
159def CSGcube(pCenter,size):
160 s2 = [0.5*size[0],0.5*size[1],0.5*size[2]]
161 p0 = VSub(pCenter,s2)
162 p1 = VAdd(pCenter,s2)
163 brick = OrthoBrick(Pnt(p0[0],p0[1],p0[2]),Pnt(p1[0],p1[1],p1[2]))
164 return brick
165
166
167#transform points
168def TransformCrank(p, zOff, zRot):
169 p2 = RotationMatrixZ(zRot) @ p
170 pOff=[0,0,zOff]
171 return VAdd(p2,pOff)
172
173#cube only in XY-plane, z infinite
174def CSGcubeXY(pCenter,sizeX,sizeY,ex,ey):
175 #print("pCenter=",pCenter)
176 pl1 = Plane(Pnt(pCenter[0]-0.5*sizeX*ex[0],pCenter[1]-0.5*sizeX*ex[1],0),Vec(-ex[0],-ex[1],-ex[2]))
177 pl2 = Plane(Pnt(pCenter[0]+0.5*sizeX*ex[0],pCenter[1]+0.5*sizeX*ex[1],0),Vec( ex[0], ex[1], ex[2]))
178
179 pl3 = Plane(Pnt(pCenter[0]-0.5*sizeY*ey[0],pCenter[1]-0.5*sizeY*ey[1],0),Vec(-ey[0],-ey[1],-ey[2]))
180 pl4 = Plane(Pnt(pCenter[0]+0.5*sizeY*ey[0],pCenter[1]+0.5*sizeY*ey[1],0),Vec( ey[0], ey[1], ey[2]))
181
182 return pl1*pl2*pl3*pl4
183
184
185#create one crank face at certain z-offset and rotation; side=1: left, side=-1: right
186def GetCrankFace(zOff, zRot, side=1):
187 ex = RotationMatrixZ(zRot) @ [1,0,0]
188 ey = RotationMatrixZ(zRot) @ [0,1,0]
189 #print("zOff=",zOff, "zRot=", zRot, "side=", side,"ex=", ex)
190 pLeft = [0,0,zOff]
191 pRight = [0,0,zOff+dk]
192 pMid = [0,0,zOff+0.5*dk]
193
194 pcLeft=VAdd(pLeft,lk*ex)
195 pcRight=VAdd(pRight,lk*ex)
196 f=0.5**0.5
197 cyl1pl = Plane(Pnt(pcLeft[0],pcLeft[1],pcLeft[2]+0.5*dk-side*dk),Vec(f*ex[0],f*ex[1],f*ex[2]-side*f))
198 cyl1 = Cylinder(Pnt(pcLeft[0],pcLeft[1],pcLeft[2]-1), Pnt(pcRight[0],pcRight[1],pcRight[2]+1), 0.5*bk)*cyl1pl
199
200 #cone2 = Cylinder(Pnt(pcLeft[0],pcLeft[1],pcLeft[2]-1), Pnt(pcRight[0],pcRight[1],pcRight[2]+1), lk+l4)
201 cone2 = Cone(Pnt(pcLeft[0],pcLeft[1],pcLeft[2]-side*dBevel+0.5*dk), Pnt(pcLeft[0],pcLeft[1],pcLeft[2]+side*dBevel+0.5*dk), lk+l4-1.5*dBevel, lk+l4-0.5*dBevel)
202 cube1 = CSGcubeXY(VAdd(pMid,0.49*l3*ex),1.02*l3,bk,ex,ey) #make l3 a little longer, to avoid bad edges
203 cube2 = CSGcubeXY(VAdd(pMid,-0.5*l4*ex),1.0*l4,bm,ex,ey)*cone2
204
205 pc3a = VAdd(pLeft,0.*l3*ex+(0.5*bk+0.4*l3)*ey)
206 cyl3a = Cylinder(Pnt(pc3a[0],pc3a[1],pc3a[2]-1), Pnt(pc3a[0],pc3a[1],pc3a[2]+1), 0.42*l3)
207 pc3b = VAdd(pLeft,0.*l3*ex+(-0.5*bk-0.4*l3)*ey)
208 cyl3b = Cylinder(Pnt(pc3b[0],pc3b[1],pc3b[2]-1), Pnt(pc3b[0],pc3b[1],pc3b[2]+1), 0.42*l3)
209 #cube3a = (CSGcubeXY(VAdd(pMid,0.26*l3*ex+(0.5*bk+0.26*l3)*ey),0.5*l3,0.5*l3,ex,ey)-cyl3a)
210
211 return ((cube1+cube2+cyl1)-(cyl3a+cyl3b))*Plane(Pnt(0,0,pLeft[2]),Vec(0,0,-1))*Plane(Pnt(0,0,pRight[2]),Vec(0,0,1))
212 #return (cube1+cube2+cyl1)*Plane(Pnt(0,0,pLeft[2]),Vec(0,0,-1))*Plane(Pnt(0,0,pRight[2]),Vec(0,0,1))
213
214#generate one crank, rotated around z-axis in radiant
215def GenerateCrank(zOff, zRot):
216 pL = TransformCrank(pLB,zOff, zRot)
217 p0 = TransformCrank(p0B,zOff, zRot)
218 p1 = TransformCrank(p1B,zOff, zRot)
219
220 p21 = TransformCrank(p21B,zOff, zRot)
221 p31 = TransformCrank(p31B,zOff, zRot)
222 p41 = TransformCrank(p41B,zOff, zRot)
223 p51 = TransformCrank(p51B,zOff, zRot)
224
225 p6 = TransformCrank(p6B,zOff, zRot)
226 p7 = TransformCrank(p7B,zOff, zRot)
227 p8 = TransformCrank(p8B,zOff, zRot)
228
229 crank0 = CSGcylinder(pL,[p0[0],p0[1],p0[2]+eps],r0)
230 crank1 = CSGcylinder(p0,[p1[0],p1[1],p1[2]+eps],rdb0)
231
232 #conrod bearing:
233 crank3 = CSGcylinder([p21[0],p21[1],p21[2]-eps],p31,rdb1)
234 crank7 = CSGcylinder(p31,p41,r1)
235 crank8 = CSGcylinder(p41,[p51[0],p51[1],p51[2]+eps],rdb1)
236
237 crank9 = CSGcylinder([p6[0],p6[1],p6[2]-eps],p7,rdb0)
238 crank10 = CSGcylinder([p7[0],p7[1],p7[2]-eps],p8,r0)
239
240 #return crank0+crank1+crank3+crank4+crank5+crank6+crank7+crank8+crank4b+crank5b+crank6b+crank9+crank10
241 if zOff==0:#add first shaft
242 crank1 = crank1+crank0
243 return crank1+GetCrankFace(db+zOff,zRot,1)+crank3+crank7+crank8+GetCrankFace(db+2*db+dk+b1+zOff,zRot,-1)+crank10+crank9
244
245
246geoCrank = CSGeometry()
247
248#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
249#choose configuration for crankshaft:
250#crankConfig = [0] #1-piston
251#crankConfig = [np.pi/2] #1-piston
252#crankConfig = [0,np.pi] #2-piston
253#crankConfig = [0,np.pi*2./3.,2.*np.pi*2./3.] #3-piston
254#crankConfig = [0,np.pi,np.pi,0] #4-piston
255crankConfig = [0,np.pi*2./3.,2.*np.pi*2./3.,2.*np.pi*2./3.,np.pi*2./3.,0] #6-piston
256#crankConfig = crankConfig*2 #12-piston
257
258nPistons = len(crankConfig)
259
260crank = GenerateCrank(0, crankConfig[0])
261zPos = lTotal
262for i in range(len(crankConfig)-1):
263 angle = crankConfig[i+1]
264 crank += GenerateCrank(zPos, angle)
265 zPos += lTotal
266
267# crank = (GenerateCrank(0, 0) + GenerateCrank(lTotal, np.pi*2./3.) + GenerateCrank(2*lTotal, np.pi*2.*2./3.)+
268# GenerateCrank(3*lTotal, np.pi*2.*2./3.) + GenerateCrank(4*lTotal, np.pi*2./3.))
269
270geoCrank.Add(crank)
271
272#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
273#conrod model:
274def GenerateConrod(zOff):
275 ey0 = [0,1,0] #top/bottom face vector of conrod
276 ey1 = [0,-1,0]
277
278 ex0 = [1,0,0] #top/bottom face vector of conrod
279 ex1 = [1,0,0]
280
281 ey0 = RotationMatrixZ(-angC)@ey0
282 ey1 = RotationMatrixZ(angC)@ey1
283 ex0 = RotationMatrixZ(-angC)@ex0
284 ex1 = RotationMatrixZ(angC)@ex1
285
286
287 pl1 = Plane(Pnt(0, 0.5*bc,0),Vec(ey0[0],ey0[1],ey0[2]))
288 pl2 = Plane(Pnt(0,-0.5*bc,0),Vec(ey1[0],ey1[1],ey1[2]))
289
290 pl3 = Plane(Pnt(-0.5*lc,0,0),Vec(-1,0,0))
291 pl4 = Plane(Pnt( 0.5*lc,0,0),Vec( 1,0,0))
292
293 pl5 = Plane(Pnt( 0,0,-0.5*dc+zOff),Vec( 0,0,-1))
294 pl6 = Plane(Pnt( 0,0, 0.5*dc+zOff),Vec( 0,0, 1))
295
296
297 cylC1 = Cylinder(Pnt(-0.5*lc,0,-1), Pnt(-0.5*lc,0,1), r1)
298 #cylC1o = Cylinder(Pnt(-0.5*lc,0,-1), Pnt(-0.5*lc,0,1), r1o)
299 cylC1o = Sphere(Pnt(-0.5*lc,0,zOff), r1o) #in fact is a sphere
300
301 cylC2 = Cylinder(Pnt( 0.5*lc,0,-1), Pnt( 0.5*lc,0,1), r2)
302 #cylC2o = Cylinder(Pnt(0.5*lc,0,-1), Pnt( 0.5*lc,0,1), r2o)
303 cylC2o = Sphere(Pnt(0.5*lc,0,zOff), r2o) #in fact is a sphere
304
305 cylSideA = (Cylinder(Pnt(-0.5*lc+r1o,0,cylOffZ+zOff), Pnt(0.5*lc-r2o,0,cylOffZ+zOff), cylR)*
306 Plane(Pnt(-0.5*lc+r1o-0.002,0,0),Vec(-1,0,0))*
307 Plane(Pnt( 0.5*lc-r2o+0.002,0,0),Vec( 1,0,0)))
308
309 cylSideB = (Cylinder(Pnt(-0.5*lc+r1o,0,-cylOffZ+zOff), Pnt(0.5*lc-r2o,0,-cylOffZ+zOff), cylR)*
310 Plane(Pnt(-0.5*lc+r1o-0.002,0,0),Vec(-1,0,0))*
311 Plane(Pnt( 0.5*lc-r2o+0.002,0,0),Vec( 1,0,0)))
312
313
314 return ((pl1*pl2*pl3*pl4+cylC1o+cylC2o)-cylC1-cylC2)*pl5*pl6-cylSideA-cylSideB
315 #return pl1*pl2*pl3*pl4*pl5*pl6
316
317geoConrod = CSGeometry()
318conrod = GenerateConrod(0)#db+dk+db+0.5*b1
319geoConrod.Add(conrod)
320
321# if netgenDrawing:
322# Draw(geoCrank)
323
324#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
325#conrod model:
326def GeneratePiston(zOff):
327 p0 = [-dpb,0,zOff]
328 p1 = [-dpb+lp,0,zOff]
329 cylPo = CSGcylinder(p0, p1, 0.5*bp) #piston outside
330 cylPaxis= CSGcylinder([0,0,-0.5*lpAxis-eps+zOff], [0,0, 0.5*lpAxis+eps+zOff], r2) #piston axis
331 cylPaxis0= CSGcylinder([0,0,-0.5*lpAxis-eps+zOff], [0,0,-0.5*lpAxis+db+zOff], r2+db) #piston axis
332 cylPaxis1= CSGcylinder([0,0, 0.5*lpAxis-db+zOff], [0,0, 0.5*lpAxis+eps+zOff], r2+db) #piston axis
333 cylPin = CSGcylinder([0,0,-0.5*lpAxis+zOff], [0,0, 0.5*lpAxis+zOff], r2p) #piston inner cutout
334
335 #box = CSGcube([0,0,zOff], [dpb+r2p,2*(r2p),lpAxis])
336 box = CSGcube([-0.5*dpb,0,zOff], [dpb,2*(r2p)-0.002,lpAxis-0.000])
337
338 cylCut = CSGcylinder([-(l4+l3+lOffCut),0,-bp+zOff], [-(l4+l3+lOffCut),0, bp+zOff], l4+l3) #piston inner cutout
339
340 return (cylPo-box-cylCut-cylPin)+cylPaxis+cylPaxis0+cylPaxis1
341
342geoPiston = CSGeometry()
343piston = GeneratePiston(0)#db+dk+db+0.5*b1
344geoPiston.Add(piston)
345
346if verbose: print("Generate meshes ...")
347#do meshing, if geometry is successful
348if True:
349 meshCrank = Mesh( geoCrank.GenerateMesh(maxh=meshSize))
350 meshCrank.Curve(1)
351 if netgenDrawing:
352 Draw(meshCrank)
353
354 if False:
355 #save mesh to file:
356 meshCrank.ngmesh.Export('testData/crankshaft.mesh','Neutral Format')
357
358if True:
359 meshConrod = Mesh( geoConrod.GenerateMesh(maxh=meshSize)) #in videos 0.003
360 meshConrod.Curve(1)
361 if netgenDrawing:
362 Draw(meshConrod)
363 if False:
364 meshConrod.ngmesh.Export('testData/conrod.mesh','Neutral Format')
365 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++
366
367if True:
368 meshPiston = Mesh( geoPiston.GenerateMesh(maxh=meshSize+0.001*0))
369 meshPiston.Curve(1)
370 if netgenDrawing:
371 Draw(meshPiston)
372 if False:
373 meshPiston.ngmesh.Export('testData/piston.mesh','Neutral Format')
374 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++
375
376#here starts the EXUDYN part
377if True:
378 SC = exu.SystemContainer()
379 mbs = SC.AddSystem()
380
381 #crankshaft and piston mechanical parameters:
382 density = 7850
383 youngsModulus = 2.1e11 *1e-1
384 poissonsRatio = 0.3
385 fRotorStart = 20 #initial revolutions per second, only crankshaft
386
387 totalFEcoordinates = 0 #accumulated FE-mesh coordinates
388 #%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
389 #import crankshaft mesh into EXUDYN FEMinterface
390 femCrank = FEMinterface()
391 eigenModesNGsolve=True
392 nModes=8
393
394 [bfM, bfK, fes] = femCrank.ImportMeshFromNGsolve(meshCrank, density, youngsModulus, poissonsRatio,
395 verbose = True, meshOrder = meshOrder)
396 # computeEigenmodes=eigenModesNGsolve, excludeRigidBodyModes = 6,
397 # numberOfModes = nModes, maxEigensolveIterations=20)
398
399 nModes = 20
400 excludeRigidBodyModes = 6
401 if verbose: print("number of coordinates crank =", femCrank.NumberOfCoordinates())
402 if verbose: print("Compute eigenmodes crank ....")
403
404 if not eigenModesNGsolve:
405 startCrank = timeit.default_timer()
406 femCrank.ComputeEigenmodes(nModes, excludeRigidBodyModes = excludeRigidBodyModes, useSparseSolver = True)
407 stopCrank = timeit.default_timer()
408 print("\ncrank eigen analysis time=", stopCrank-startCrank)
409 else:
410 start_time = time.time()
411 femCrank.ComputeEigenmodesNGsolve(bfM, bfK, nModes=nModes,
412 excludeRigidBodyModes=excludeRigidBodyModes, maxEigensolveIterations=20)
413 print("NGsolve mode computation needed %.3f seconds" % (time.time() - start_time))
414
415 totalFEcoordinates+=femCrank.NumberOfCoordinates()
416 print("eigen freq. crank=", femCrank.GetEigenFrequenciesHz()[0:nModes])
417
418 #+++++++++++++++++++++++++++++++++++++++++++++++++++++
419 #compute stress modes:
420 SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
421 mat = KirchhoffMaterial(youngsModulus, poissonsRatio, density)
422 varType = exu.OutputVariableType.DisplacementLocal
423 #varType = exu.OutputVariableType.StrainLocal
424 if showStresses:
425 print("ComputePostProcessingModes femCrank ... ")
426 start_time = time.time()
427 varType = exu.OutputVariableType.StressLocal
428 femCrank.ComputePostProcessingModesNGsolve(fes, material=mat,
429 outputVariableType=varType)
430 print("--- %s seconds ---" % (time.time() - start_time))
431
432 SC.visualizationSettings.contour.outputVariable = varType
433
434 #print("Create CMS object and matrices ....")
435 cmsCrank = ObjectFFRFreducedOrderInterface(femCrank)
436
437 #user functions should be defined outside of class:
438 def UFmassFFRFreducedOrderCrank(mbs, t, itemIndex, qReduced, qReduced_t):
439 return cmsCrank.UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
440
441 def UFforceFFRFreducedOrderCrank(mbs, t, itemIndex, qReduced, qReduced_t):
442 return cmsCrank.UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
443
444 objFFRFcrank = cmsCrank.AddObjectFFRFreducedOrderWithUserFunctions(exu, mbs,
445 positionRef=[0,0,0],
446 eulerParametersRef=eulerParameters0,
447 initialVelocity=[0,0,0], initialAngularVelocity=[0,0,1*fRotorStart*2*pi],
448 gravity = [0,-0*9.81,0],
449 #UFforce=UFforceFFRFreducedOrderCrank,
450 #UFmassMatrix=UFmassFFRFreducedOrderCrank,
451 color=[0.1,0.9,0.1,1.])
452 mbs.SetObjectParameter(objFFRFcrank['oFFRFreducedOrder'],'VshowNodes',False)
453
454
455 if False:#animate eigenmodes of crankshaft
456 from exudyn.interactive import AnimateModes
457 mbs.Assemble()
458
459 SC.visualizationSettings.general.textSize = 16 #30 for cover figure
460 SC.visualizationSettings.general.useGradientBackground = True
461 SC.visualizationSettings.openGL.lineWidth = 2
462 SC.visualizationSettings.openGL.showFaceEdges = True
463 SC.visualizationSettings.openGL.showFaces = True
464 SC.visualizationSettings.openGL.multiSampling = 4
465 SC.visualizationSettings.nodes.show = False
466 SC.visualizationSettings.window.renderWindowSize = [1600,1080]
467
468 SC.visualizationSettings.contour.outputVariableComponent = 0
469
470 SC.visualizationSettings.general.autoFitScene=False
471
472 AnimateModes(SC, mbs, 1, period=0.2)
473 exit()
474
475 #%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
476 #import conrod and piston mesh into EXUDYN FEMinterface and compute eigenmodes
477 nModes = 8
478 excludeRigidBodyModes = 6
479 femConrod = FEMinterface()
480 # femConrod.ImportMeshFromNGsolve(meshConrod, density, youngsModulus, poissonsRatio, verbose = False)
481 [bfM, bfK, fes] = femConrod.ImportMeshFromNGsolve(meshConrod, density, youngsModulus, poissonsRatio,
482 verbose = False, meshOrder = meshOrder)
483 # computeEigenmodes=eigenModesNGsolve, excludeRigidBodyModes = 6,
484 # numberOfModes = nModes, maxEigensolveIterations=20)
485 if verbose: print("number of coordinates conrod =", femConrod.NumberOfCoordinates())
486 if verbose: print("Compute eigenmodes conrod ....")
487
488 if not eigenModesNGsolve:
489 femConrod.ComputeEigenmodes(nModes, excludeRigidBodyModes = excludeRigidBodyModes, useSparseSolver = True)
490 else:
491 femConrod.ComputeEigenmodesNGsolve(bfM, bfK, nModes=nModes, excludeRigidBodyModes=excludeRigidBodyModes)
492
493 totalFEcoordinates+=femConrod.NumberOfCoordinates()
494 if verbose: print("eigen freq. conrod=", femConrod.GetEigenFrequenciesHz()[0:nModes])
495
496 if showStresses:
497 print("ComputePostProcessingModes femConrod ... ")
498 start_time = time.time()
499 femConrod.ComputePostProcessingModesNGsolve(fes, material=mat,
500 outputVariableType=varType)
501 print("--- %s seconds ---" % (time.time() - start_time))
502
503 #%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
504 #import piston mesh into EXUDYN FEMinterface
505 femPiston = FEMinterface()
506 #femPiston.ImportMeshFromNGsolve(meshPiston, density, youngsModulus, poissonsRatio, verbose = False)
507 [bfM, bfK, fes] = femPiston.ImportMeshFromNGsolve(meshPiston, density, youngsModulus, poissonsRatio, verbose = False, meshOrder = meshOrder)
508
509 if verbose: print("number of coordinates piston =", femPiston.NumberOfCoordinates())
510 if verbose: print("Compute eigenmodes piston ....")
511
512 if not eigenModesNGsolve:
513 femPiston.ComputeEigenmodes(nModes, excludeRigidBodyModes = excludeRigidBodyModes, useSparseSolver = True)
514 else:
515 femPiston.ComputeEigenmodesNGsolve(bfM, bfK, nModes=nModes, excludeRigidBodyModes=excludeRigidBodyModes)
516
517 totalFEcoordinates+=femPiston.NumberOfCoordinates()
518 if verbose: print("eigen freq. Piston=", femPiston.GetEigenFrequenciesHz()[0:nModes])
519
520 if showStresses:
521 print("ComputePostProcessingModes femPiston ... ")
522 start_time = time.time()
523 femPiston.ComputePostProcessingModesNGsolve(fes, material=mat,
524 outputVariableType=varType)
525 print("--- %s seconds ---" % (time.time() - start_time))
526
527 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
528 #import multiple conrods and pistons
529
530 #user functions should be defined outside of class:
531 def UFmassFFRFreducedOrderConrod0(mbs, t, itemIndex, qReduced, qReduced_t):
532 return cmsConrodList[0].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
533 def UFmassFFRFreducedOrderConrod1(mbs, t, itemIndex, qReduced, qReduced_t):
534 return cmsConrodList[1].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
535 def UFmassFFRFreducedOrderConrod2(mbs, t, itemIndex, qReduced, qReduced_t):
536 return cmsConrodList[2].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
537 def UFmassFFRFreducedOrderConrod3(mbs, t, itemIndex, qReduced, qReduced_t):
538 return cmsConrodList[3].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
539 def UFmassFFRFreducedOrderConrod4(mbs, t, itemIndex, qReduced, qReduced_t):
540 return cmsConrodList[4].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
541 def UFmassFFRFreducedOrderConrod5(mbs, t, itemIndex, qReduced, qReduced_t):
542 return cmsConrodList[5].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
543
544 def UFforceFFRFreducedOrderConrod0(mbs, t, itemIndex, qReduced, qReduced_t):
545 return cmsConrodList[0].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
546 def UFforceFFRFreducedOrderConrod1(mbs, t, itemIndex, qReduced, qReduced_t):
547 return cmsConrodList[1].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
548 def UFforceFFRFreducedOrderConrod2(mbs, t, itemIndex, qReduced, qReduced_t):
549 return cmsConrodList[2].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
550 def UFforceFFRFreducedOrderConrod3(mbs, t, itemIndex, qReduced, qReduced_t):
551 return cmsConrodList[3].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
552 def UFforceFFRFreducedOrderConrod4(mbs, t, itemIndex, qReduced, qReduced_t):
553 return cmsConrodList[4].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
554 def UFforceFFRFreducedOrderConrod5(mbs, t, itemIndex, qReduced, qReduced_t):
555 return cmsConrodList[5].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
556
557 #user functions should be defined outside of class:
558 def UFmassFFRFreducedOrderPiston0(mbs, t, itemIndex, qReduced, qReduced_t):
559 return cmsPistonList[0].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
560 def UFmassFFRFreducedOrderPiston1(mbs, t, itemIndex, qReduced, qReduced_t):
561 return cmsPistonList[1].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
562 def UFmassFFRFreducedOrderPiston2(mbs, t, itemIndex, qReduced, qReduced_t):
563 return cmsPistonList[2].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
564 def UFmassFFRFreducedOrderPiston3(mbs, t, itemIndex, qReduced, qReduced_t):
565 return cmsPistonList[3].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
566 def UFmassFFRFreducedOrderPiston4(mbs, t, itemIndex, qReduced, qReduced_t):
567 return cmsPistonList[4].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
568 def UFmassFFRFreducedOrderPiston5(mbs, t, itemIndex, qReduced, qReduced_t):
569 return cmsPistonList[5].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
570
571 def UFforceFFRFreducedOrderPiston0(mbs, t, itemIndex, qReduced, qReduced_t):
572 return cmsPistonList[0].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
573 def UFforceFFRFreducedOrderPiston1(mbs, t, itemIndex, qReduced, qReduced_t):
574 return cmsPistonList[1].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
575 def UFforceFFRFreducedOrderPiston2(mbs, t, itemIndex, qReduced, qReduced_t):
576 return cmsPistonList[2].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
577 def UFforceFFRFreducedOrderPiston3(mbs, t, itemIndex, qReduced, qReduced_t):
578 return cmsPistonList[3].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
579 def UFforceFFRFreducedOrderPiston4(mbs, t, itemIndex, qReduced, qReduced_t):
580 return cmsPistonList[4].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
581 def UFforceFFRFreducedOrderPiston5(mbs, t, itemIndex, qReduced, qReduced_t):
582 return cmsPistonList[5].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
583
584 #lists for multiple objects in conrods and pistons:
585 UFmassFFRFreducedOrderConrodList=[UFmassFFRFreducedOrderConrod0,UFmassFFRFreducedOrderConrod1,
586 UFmassFFRFreducedOrderConrod2,UFmassFFRFreducedOrderConrod3,
587 UFmassFFRFreducedOrderConrod4,UFmassFFRFreducedOrderConrod5]
588 UFforceFFRFreducedOrderConrodList=[UFforceFFRFreducedOrderConrod0,UFforceFFRFreducedOrderConrod1,
589 UFforceFFRFreducedOrderConrod2,UFforceFFRFreducedOrderConrod3,
590 UFforceFFRFreducedOrderConrod4,UFforceFFRFreducedOrderConrod5]
591 objFFRFconrodList=[]
592 cmsConrodList=[]
593 UFmassFFRFreducedOrderPistonList=[UFmassFFRFreducedOrderPiston0,UFmassFFRFreducedOrderPiston1,
594 UFmassFFRFreducedOrderPiston2,UFmassFFRFreducedOrderPiston3,
595 UFmassFFRFreducedOrderPiston4,UFmassFFRFreducedOrderPiston5]
596 UFforceFFRFreducedOrderPistonList=[UFforceFFRFreducedOrderPiston0,UFforceFFRFreducedOrderPiston1,
597 UFforceFFRFreducedOrderPiston2,UFforceFFRFreducedOrderPiston3,
598 UFforceFFRFreducedOrderPiston4,UFforceFFRFreducedOrderPiston5]
599 objFFRFpistonList=[]
600 cmsPistonList=[]
601 pkList = []
602 pcList = []
603 ppList = []
604 zOffsetList = []
605 for iCrank in range(len(crankConfig)):
606 zOffset = db+dk+db + lTotal*iCrank #left end of conrod, for multiple conrods in a loop
607 zOffsetList.append(zOffset)
608 #compute crank (pK), conrod (pC) and piston position (pP) for any crank angle:
609 phi = crankConfig[iCrank]
610 pK = np.array([lk*np.cos(phi),lk*np.sin(phi),0])
611 alpha=np.arcsin(pK[1]/lc)
612 pC = pK + np.array([0.5*lc*np.cos(alpha),-0.5*lc*np.sin(alpha),0])
613 pP = pK + np.array([lc*np.cos(alpha),-lc*np.sin(alpha),0])
614 pkList.append(pK)
615 pcList.append(pC)
616 ppList.append(pP)
617 #print("pK=",pK)
618 #print("pC=",pC)
619 #print("pP=",pP)
620
621 eulerParametersInit = RotationMatrix2EulerParameters(RotationMatrixZ(-alpha))
622 #pRef = [lk+0.5*lc,0,zOffset+0.5*b1] #0-degree
623 pRef = pC + [0,0,zOffset+0.5*b1]
624
625 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
626 #import conrod CMS
627 cmsConrod = ObjectFFRFreducedOrderInterface(femConrod)
628 cmsConrodList.append(cmsConrod)
629 objFFRFconrod = cmsConrod.AddObjectFFRFreducedOrderWithUserFunctions(exu, mbs,
630 positionRef=pRef,
631 eulerParametersRef=eulerParametersInit,
632 initialVelocity=[0,0,0],
633 initialAngularVelocity=[0,0,0*fRotorStart*2*pi],
634 gravity = [0,-0*9.81,0],
635 #UFforce=UFforceFFRFreducedOrderConrodList[iCrank],
636 #UFmassMatrix=UFmassFFRFreducedOrderConrodList[iCrank],
637 color=[0.1,0.9,0.1,1.])
638 mbs.SetObjectParameter(objFFRFconrod['oFFRFreducedOrder'],'VshowNodes',False)
639 objFFRFconrodList.append(objFFRFconrod)
640
641 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
642 #import piston CMS
643 cmsPiston = ObjectFFRFreducedOrderInterface(femPiston)
644 cmsPistonList.append(cmsPiston)
645
646 objFFRFpiston = cmsPiston.AddObjectFFRFreducedOrderWithUserFunctions(exu, mbs,
647 positionRef=pP+[0,0,zOffset+0.5*b1],
648 eulerParametersRef=eulerParameters0,
649 initialVelocity=[0,0,0], initialAngularVelocity=[0,0,0*fRotorStart*2*pi],
650 gravity = [0,-0*9.81,0],
651 #UFforce=UFforceFFRFreducedOrderPistonList[iCrank],
652 #UFmassMatrix=UFmassFFRFreducedOrderPistonList[iCrank],
653 color=[0.1,0.9,0.1,1.])
654 mbs.SetObjectParameter(objFFRFpiston['oFFRFreducedOrder'],'VshowNodes',False)
655 objFFRFpistonList.append(objFFRFpiston)
656
657 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
658 if True: #connect bodies:
659 k = 1e6 #joint stiffness
660 d = k*0.002 #joint damping
661 nMarkerPerPiston = 10 #number of markers per crank/conrod/piston part
662
663 genMarkerPos = [[0,0,-d0],[0,0,lTotal*nPistons]]
664 genMarkerR = [r0,r0]
665 genMarkerFEM = [femCrank,femCrank]
666 genMarkerObject = [objFFRFcrank,objFFRFcrank]
667
668 for iCrank in range(len(crankConfig)):
669 genMarkerPos += [pkList[iCrank]+[0,0,zOffsetList[iCrank]],pkList[iCrank]+[0,0,zOffsetList[iCrank]+b1],
670 [-0.5*lc,0,-0.5*dc],[-0.5*lc,0, 0.5*dc],[0.5*lc,0,-0.5*dc],[0.5*lc,0, 0.5*dc],
671 [0,0,-0.5*dc],[0,0,0.5*dc], [-dpb,0,0],[lp-dpb,0,0]]
672 genMarkerR += [r1,r1,
673 r1,r1,r2,r2,
674 r2,r2,0.5*bp,0.5*bp]
675 genMarkerFEM += [femCrank,femCrank,
676 femConrod,femConrod,femConrod,femConrod,
677 femPiston,femPiston,femPiston,femPiston]
678 genMarkerObject += [objFFRFcrank,objFFRFcrank,
679 objFFRFconrodList[iCrank],objFFRFconrodList[iCrank],objFFRFconrodList[iCrank],objFFRFconrodList[iCrank],
680 objFFRFpistonList[iCrank],objFFRFpistonList[iCrank],objFFRFpistonList[iCrank],objFFRFpistonList[iCrank]]
681
682 markerList = []
683 #generate markers for joints:
684 for i in range(len(genMarkerPos)):
685 p = genMarkerPos[i]
686 nodeList=[]
687 if p[2] != 0:
688 nodeList= genMarkerFEM[i].GetNodesOnCircle(p, [0,0,1], genMarkerR[i])
689 else:
690 nodeList= genMarkerFEM[i].GetNodesOnCircle(p, [1,0,0], genMarkerR[i])
691 #print("nodeList"+str(i)+":", nodeList)
692 lenNodeList = len(nodeList)
693 weights = np.array((1./lenNodeList)*np.ones(lenNodeList))
694
695 markerList += [mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=genMarkerObject[i]['oFFRFreducedOrder'],
696 meshNodeNumbers=np.array(nodeList), #these are the meshNodeNumbers
697 weightingFactors=weights))]
698
699 oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
700
701 mGroundPosLeft = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=genMarkerPos[0]))
702 mGroundPosRight = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=genMarkerPos[1]))
703
704
705 #joints for crankshaft/ground
706 oSJleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGroundPosLeft, markerList[0]],
707 stiffness=[k,k,k], damping=[d,d,d]))
708 oSJright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGroundPosRight, markerList[1]],
709 stiffness=[k,k,k], damping=[d,d,d]))
710
711 for iCrank in range(len(crankConfig)):
712 mOff = nMarkerPerPiston*iCrank
713 #joints for crankshaft/conrod:
714 oJointCCleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[markerList[mOff+2], markerList[mOff+4]],
715 stiffness=[k,k,k], damping=[d,d,d]))
716 oJointCCright= mbs.AddObject(CartesianSpringDamper(markerNumbers=[markerList[mOff+3], markerList[mOff+5]],
717 stiffness=[k,k,k], damping=[d,d,d]))
718
719 #joints for conrod/piston:
720 oJointCPleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[markerList[mOff+6], markerList[mOff+8]],
721 stiffness=[k,k,k], damping=[d,d,d]))
722 oJointCPright= mbs.AddObject(CartesianSpringDamper(markerNumbers=[markerList[mOff+7], markerList[mOff+9]],
723 stiffness=[k,k,k], damping=[d,d,d]))
724
725 mGroundPosPiston = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround,
726 localPosition=[ppList[iCrank][0],0,zOffsetList[iCrank]+0.5*b1]))
727 oJointPGleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGroundPosPiston, markerList[mOff+10]],
728 stiffness=[0,k,k], damping=[0,d,d]))
729 oJointPGright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGroundPosPiston, markerList[mOff+11]],
730 stiffness=[0,k,k], damping=[0,d,d]))
731
732
733 stopTotal = timeit.default_timer()
734 print("\ntotal elapsed time=", stopTotal-startTotal)
735 mbs.Assemble()
736
737 #now simulate model in exudyn:
738 #%%+++++++++++++++++++++
739 if True:
740 print("totalFEcoordinates=",totalFEcoordinates)
741
742 simulationSettings = exu.SimulationSettings()
743
744 nodeDrawSize = 0.0005
745 SC.visualizationSettings.general.textSize = 14 #30 for cover figure
746 SC.visualizationSettings.general.useGradientBackground = True
747 SC.visualizationSettings.openGL.lineWidth = 2
748
749 SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
750 SC.visualizationSettings.nodes.drawNodesAsPoint = False
751 SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
752 SC.visualizationSettings.connectors.show = False
753
754 SC.visualizationSettings.nodes.show = False
755 SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
756 SC.visualizationSettings.nodes.basisSize = 0.12
757 SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
758
759 SC.visualizationSettings.openGL.showFaceEdges = True
760 SC.visualizationSettings.openGL.showFaces = True
761 SC.visualizationSettings.openGL.multiSampling = 4
762
763 SC.visualizationSettings.sensors.show = True
764 SC.visualizationSettings.sensors.drawSimplified = False
765 SC.visualizationSettings.sensors.defaultSize = 0.01
766 SC.visualizationSettings.markers.drawSimplified = False
767 SC.visualizationSettings.markers.show = False
768 SC.visualizationSettings.markers.defaultSize = 0.01
769
770 SC.visualizationSettings.loads.drawSimplified = False
771
772 #SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
773 SC.visualizationSettings.contour.outputVariableComponent = -1
774 SC.visualizationSettings.contour.reduceRange = True
775 #SC.visualizationSettings.contour.automaticRange = False
776 #SC.visualizationSettings.contour.maxValue = 3e7
777 # SC.visualizationSettings.contour.minValue = -0.0003
778 # SC.visualizationSettings.contour.maxValue = 0.0003
779
780 simulationSettings.solutionSettings.solutionInformation = "NGsolve/NETGEN engine test"
781
782 h=0.05e-3
783 tEnd = 2
784
785 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
786 simulationSettings.timeIntegration.endTime = tEnd
787 simulationSettings.solutionSettings.solutionWritePeriod = h*10 #writing already costs much time
788 simulationSettings.timeIntegration.verboseMode = 1
789 #simulationSettings.timeIntegration.verboseModeFile = 3
790 simulationSettings.timeIntegration.newton.useModifiedNewton = True
791
792 simulationSettings.solutionSettings.sensorsWritePeriod = h
793 #simulationSettings.solutionSettings.coordinatesSolutionFileName = "solution/coordinatesSolution.txt"
794 simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse #faster, because system size already quite large
795
796 simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5 #SHOULD work with 0.9 as well
797 simulationSettings.displayStatistics = True
798 #simulationSettings.displayComputationTime = True
799 SC.visualizationSettings.general.autoFitScene = False #for reloading of renderState to work
800
801 #create animation:
802 if False:
803 simulationSettings.solutionSettings.recordImagesInterval = 0.001
804 SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
805 SC.visualizationSettings.window.renderWindowSize=[1920,1080]
806
807 exu.StartRenderer()
808 if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
809
810 mbs.WaitForUserToContinue() #press space to continue
811
812 simulate = True #set false to show last stored solution
813 if simulate:
814 mbs.SolveDynamic(simulationSettings)
815 else:
816 SC.visualizationSettings.general.autoFitScene = False
817 sol = LoadSolutionFile('coordinatesSolution.txt')
818 if False: #directly show animation
819 AnimateSolution(mbs, solution=sol, rowIncrement = 1, timeout=0.01,
820 createImages = False, runLoop = True)
821 else: #interact with animation
822
823 mbs.SolutionViewer(sol, rowIncrement=1, timeout=0.02)
824
825
826 if False: #draw with matplotlib, export as pdf
827 SC.visualizationSettings.exportImages.saveImageFormat = "TXT"
828 SC.visualizationSettings.exportImages.saveImageAsTextTriangles=True
829 SC.RedrawAndSaveImage() #uses default filename
830
831 from exudyn.plot import LoadImage, PlotImage
832
833 # plot 2D
834 # data = LoadImage('images/frame00000.txt', trianglesAsLines=True)
835 # PlotImage(data, HT=HomogeneousTransformation(RotationMatrixZ(0.5*pi)@RotationMatrixX(0.5*pi), [0,0,0]),
836 # lineWidths=0.5, lineStyles='-', title='', closeAll=True, plot3D=False,
837 # fileName='images/test.pdf')
838
839 data = LoadImage('images/frame00000.txt', trianglesAsLines=False)
840 PlotImage(data, HT=HomogeneousTransformation(2.5*RotationMatrixZ(0.5*pi)@RotationMatrixY(-0.5*pi), [0,1,0.25]),
841 lineWidths=0.5, lineStyles='-', triangleEdgeColors='black', triangleEdgeWidths=0.25, title='', closeAll=True, plot3D=True,
842 fileName='images/test3D.pdf')
843
844 SC.WaitForRenderEngineStopFlag()
845 exu.StopRenderer() #safely close rendering window!
846 lastRenderState = SC.GetRenderState() #store model view for next simulation