← Back to team overview

yade-users team mailing list archive

[Question #661215]: Record the rotation of polyhedra by VTK and plot the i vs rot

 

New question #661215 on Yade:
https://answers.launchpad.net/yade/+question/661215

Dear everybody
I am going on simulating a process of compact ballast by a plate, I want to record the rotation of one ballast of the sample, and make the plot of i vs the rotation ,but I got nothing, the process was simply show below:

>from yade import qt,plot,utils,polyhedra_utils,ymport,export,pack,timing
>from yade import *
>import numpy

>## materials ##
>global gravel
>global steel
vglobal rubber
>gravel = PolyhedraMat()
>...omit

>## objects ##
># loading plate name 'laplt'
>global ldplt
>global ldpltheight
>#thickness of plate is 20mm
>ldpltheight=0.95
>ldplt=polyhedra_utils.polyhedra(steel,v=((-0.15,0,ldpltheight),(-0.075*sqrt(3),-0.075,ldpltheight),(-0.075,->0.075*sqrt(3),ldpltheight),...omit)
>O.bodies.append(ldplt) 

>O.bodies.append(utils.wall(0,axis=2,sense=1, material = steel)) 

>O.bodies.append(geom.facetCylinder(center=(0.0,0.0,0.15),radius=0.152,height=0.3,orientation=Quaternion((0,0,1),0),segmentsNumber=12,wallMask=4,angleRange=(0.0,2*pi),closeGap=False,radiusTopInner=0.0,radiusBottomInner=0.0,material=steel)) #membrane of samples, radius=0.152>radius of loading plate 0.15 for 2mm

>O.bodies.append(geom.facetBunker((0,0,0.3),0.40,0.304,hBunker=0.3,hOutput=0.05,hPipe=0.05,orientation=Quaternion((0,0,1),0),segmentsNumber=12,wallMask=4,angleRange=(0,2*pi),closeGap=False,material=steel)) #bunker for generate samples

>#global ballast_up
>ballast_up=polyhedra_utils.fillBox((-0.15,-0.15,0.4), (0.15,0.15,0.7),gravel,sizemin=[0.025,0.025,0.025],sizemax=[0.05,0.05,0.05],ratio=[1,1,1],seed=4,mask=1)
>#global ballast_low
>ballast_low=polyhedra_utils.fillBox((-0.11,-0.11,0.0), (0.11,0.11,0.4),gravel,sizemin=[0.025,0.025,0.025],sizemax=[0.05,0.05,0.05],ratio=[1,1,1],seed=4,mask=1)


>## engines ##
>O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb()]),
	InteractionLoop(
		[Ig2_Wall_Polyhedra_PolyhedraGeom(), Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
		[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys(), Ip2_FrictMat_PolyhedraMat_FrictPhys()],
		[Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()],   # contact law -- apply forces #Bo1_Cylinder_Aabb()
	),
	NewtonIntegrator(damping=0.3,gravity=(0,0,-9.81)),
	PolyhedraSplitter(iterPeriod=1,label='Splitter'), # allow splittering ballast
	PyRunner(command="vtkExporter.exportPolyhedra(what=[('num','b.id'),('rotate','b.state.ori')])",realPeriod=10),
       ### What' s the problem of this VTK line???
	PyRunner(command='checkUnbalancedI()',realPeriod=1,label='checker')
]

>vtkExporter = export.VTKExporter('postfalling/Jan-')

>O.dt=0.1*polyhedra_utils.PWaveTimeStep()

## selfdefine pyrunner attention the format of the language
>def checkUnbalancedI():
	print "iter %d, time elapsed %f,  time step %.5e, unbalanced forces = %.5f"%(O.iter, O.realtime, O.dt, utils.unbalancedForce())
	if O.iter<6000: return 
	if unbalancedForce()>0.05: return
	O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=50)]+[PyRunner(command='checkersieve()',iterPeriod=500)]
	# next time, do not call this function anymore, but the next one (unloadPlate) instead
	checker.command='cyclicloadingI()'


>def addPlotData():
	Fz=265+135*sin((15/pi)*O.time)
	#Fz=14350+10250*sin((15/pi)*O.time) 
	#Fz=-O.forces.f(ldplt.id)[2]
	W=-ldplt.state.pos[2]+ldplt.state.refPos[2]
	BalRota=O.bodies[45].state.rot()                 ##### What's wrong with this line?
	plot.addData(Fz=Fz,w=W,br=BalRota,i=O.iter,unbalanced=unbalancedForce(),**O.energy) 

>def cyclicloadingI():
	yade.timing.reset()
	O.forces.setPermF(O.bodies[0].id,(0,0,-265-135*sin((15/pi)*O.time))) # 130N to 400N

>def checkersieve():
	polyhedra_utils.SieveCurve()
	polyhedra_utils.SizeRatio()

#from yade import qt
qt.Controller()
V = qt.View()

#plot
plot.plots={'i':('unbalanced',),'w':('Fz',),'i':('br',)}
plot.plot()


## save and load ## 
O.wait()
#O.save('ballastCompact.yade.gz')
O.run(10000)
O.wait()




-- 
You received this question notification because your team yade-users is
an answer contact for Yade.