yade-users team mailing list archive

[De zhang <question661262@xxxxxxxxxxxxxxxxxxxxx>]

Question #661262 on Yade changed:
https://answers.launchpad.net/yade/+question/661262

Description changed to:
I made a simulation of polyhedra compaction model and I want to record the rotation angle of polyhedra,
I tried :
> def addPlotData():
>   Fz=100+50*sin((15/pi)*O.time)
>   w=-plate.state.pos[2]+plate.state.refPos[2]
>   s=O.bodies[45].state.angVel
>   plot.addData(Fz=Fz,w=w,s=s)

I also tried 'state.rot' and 'state.ori' , but can't get any information of rotate angle of polyhedra? How to deal with this? 
The whole process is fllowing:

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

## materials ##
global gravel
global steel
global rubber
gravel = PolyhedraMat()
gravel.density = 2600 #kg/m^3 
gravel.young = 1E7 #Pa
gravel.poisson = 20000/1E7
gravel.frictionAngle = 0.5 #rad
gravel.strength = 1E9 # Pa crushable 1000MPa

steel = PolyhedraMat()
steel.density = 7850 #kg/m^3 
steel.young = 10*gravel.young #inital steel was 10*gravel.young
steel.poisson = gravel.poisson
steel.frictionAngle = 0.4 #rad

rubber = PolyhedraMat()
rubber.density = 1000 #kg/m^3 
rubber.young = gravel.young/10
rubber.poisson = gravel.poisson
rubber.frictionAngle = 0.7 #rad

## 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),(0,-0.15,ldpltheight),(0.075,-0.075*sqrt(3),ldpltheight),(0.075*sqrt(3),-0.075,ldpltheight),(0.15,0,ldpltheight),(0.075*sqrt(3),0.075,ldpltheight),(0.075,0.075*sqrt(3),ldpltheight),(0,0.15,ldpltheight),(-0.075,0.075*sqrt(3),ldpltheight),(-0.075*sqrt(3),0.075,ldpltheight),(-0.15,0,ldpltheight+0.02),(-0.075*sqrt(3),-0.075,ldpltheight+0.02),(-0.075,-0.075*sqrt(3),ldpltheight+0.02),(0,-0.15,ldpltheight+0.02),(0.075,-0.075*sqrt(3),ldpltheight+0.02),(0.075*sqrt(3),-0.075,ldpltheight+0.02),(0.15,0,ldpltheight+0.02),(0.075*sqrt(3),0.075,ldpltheight+0.02),(0.075,0.075*sqrt(3),ldpltheight+0.02),(0,0.15,ldpltheight+0.02),(-0.075,0.075*sqrt(3),ldpltheight+0.022),(-0.075*sqrt(3),0.075,ldpltheight+0.02)),fixed=False, color=(0.75,0.65,0.65))
O.bodies.append(ldplt) # area is 0.082 m^2 & weight is 7850*0.082*0.02=12.874kg

O.bodies.append(utils.wall(0,axis=2,sense=1, material = steel)) #bottom
of samples, wall objects are 'fixed' by default, i.e. not subject to
forces

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=[('n','b.state.ori[1]')])",realPeriod=10),
	PyRunner(command='checkUnbalancedI()',realPeriod=1,label='checker')
	#PyRunner(command='VTKExporter.exportFacet(what=[('n','b.id'),,('n','b.state.pos')])',realPeriod=50)
]

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]
	s=O.bodies[45].state.angVel
	plot.addData(Fz=Fz,w=W,s=s,i=O.iter,unbalanced=unbalancedForce(),**O.energy) #w=ldplt.state.pos[2]-ldplt.state.refPos[2]

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()
V.screenSize = (550,450)
V.sceneRadius = 1
V.eyePosition = (0.7,0.5,0.1)
V.upVector = (0,0,1)
V.lookAt = (0.15,0.15,0.1)

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


## save and load ## y
#O.run(30000)
O.wait()
#O.save('Jantest.yade.gz')
O.run(10000)
O.wait()

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