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