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