New question #683384 on Yade:
https://answers.launchpad.net/yade/+question/683384
Hi all,
In my case, I want to simulate the progressive loss of fine particles in the
soil assembly, the code for this case is attached as follows. After running
1000 steps, we can see the fines indeed get lost from the visualization.
However, when using command "len(O.bodies)" in the terminal, we get the
particle number remains unchanged. i.e. n=400. How can I monitor and output the
instaneous particle number during simulation? Thanks very much!
Best,
Zheng
# the code as follows
from yade import pack,plot,qt,export
import matplotlib.pyplot as plt
import numpy as np
import random
#O.materials.append(FrictMat(young=6.e8,poisson=.8,frictionAngle=.0))
O.materials.append(CohFrictMat(young=1e9,poisson=.8,frictionAngle=.0,normalCohesion=1e6,shearCohesion=1e6,momentRotationLaw=False,etaRoll=.1,label='spheres'))
sigmaIso=-1e5
sp = pack.SpherePack()
size =0.24
sp.makeCloud(minCorner=(0,0,.05),maxCorner=(size,size,.05),rMean=.005,rRelFuzz=.4,num=400,periodic=True,seed=1)
# initial 400
#sp.makeCloud(minCorner=(0,0,.05),maxCorner=(size,size,.05),psdSizes=[0.006,0.0068,0.0072,0.0084,0.0092,0.01,0.0108,0.0116,0.0124,0.0132,0.014],
psdCumm=[0,0.1,0.2,
0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0],periodic=True,seed=1,distributeMass=False)#num=400,
# if minCorner[k]=maxCorner[k] for one coordinate, it means 2D case;
sp.toSimulation()
O.cell.hSize = Matrix3(size,0,0, 0,size,0, 0,0,.1) # used for periodic
boundaries.
#print len(O.bodies)
for p in O.bodies:
p.state.blockedDOFs = 'zXY' # a sphere can be made to move only in xy plane
p.state.mass = 2650 * 0.1 * pi * p.shape.radius**2 # 0.1 = thickness of
cylindrical particle
inertia = 0.5 * p.state.mass * p.shape.radius**2
p.state.inertia = (.5*inertia,.5*inertia,inertia)
O.dt = utils.PWaveTimeStep()
print O.dt
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D()],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(always_use_moment_law=False)]
),
PeriTriaxController(
dynCell=True,
goal=(sigmaIso,sigmaIso,0),
stressMask=3,
relStressTol=.001,
maxUnbalanced=.001,
maxStrainRate=(.5,.5,.0),
doneHook='shear()',# function to run when finished
label='biax'
),
NewtonIntegrator(damping=.1),
PyRunner(command='delByNum()',iterPeriod=100),
PyRunner(command='addPlotData()',iterPeriod=100),
#PyRunner(command='deleteFines()',iterPeriod=100),
#PyRunner(command='delBelowPerc()',iterPeriod=1),
#PyRunner(command='print ')
]
plot.live=True
plot.plots={'iter':('sxx','syy'),'iter_':('exx','eyy'),'iter':('poros')}
def addPlotData():
plot.addData(
iter=O.iter,iter_=O.iter,
sxx=biax.stress[0],syy=biax.stress[1],
exx=biax.strain[0],eyy=biax.strain[1],
Z=avgNumInteractions(),
Zm=avgNumInteractions(skipFree=True),
poro=porosity(),
unbalanced=utils.unbalancedForce(),
t=O.time
)
plot.saveDataTxt('results',vars=('t','exx','eyy','sxx','syy','Z','Zm','poro'))
"""
# delete fines by percent
def delBelowPerc():
bodyRadius=[]
for b in O.bodies:
if b.shape.radius<=0.005:
#if isinstance(b.shape,Sphere):
bodyRadius.append([b,b.shape.radius])
bodyRadius.sort(key=lambda x:x[1])
if len(bodyRadius)>1:
maxRad=np.percentile(bodyRadius,10)
for b in bodyRadius:
if b[0].shape.radius<=maxRad:
O.bodies.erase(b[0].id)
"""
def delByNum():
setContactFriction(0.5)
bodyRadius=[]
for b in O.bodies:
if b.shape.radius<=0.005: # define fine particles
#if isinstance(b.shape,Sphere):
bodyRadius.append([b,b.shape.radius])
bodyRadius.sort(key=lambda x:x[1])
l=len(bodyRadius) # how many fine particles
perc=0.1
delNum=int(l*perc) # how many number to be deleted
list=random.sample(bodyRadius,delNum) # list to be deleted
for b in list:
O.bodies.erase(b[0].id)
print 'delete fines'
#biax.doneHook='shear()'
"""
def deleteFines():
#new added here by Zheng 0422
bodiesToBeDeleted=[]
for b in O.bodies:
if b.shape.radius<=0.004:
bodiesToBeDeleted.append(b)
for b in bodiesToBeDeleted:
O.bodies.erase(b.id)
#new added here by Zheng 0422
def term0():
print getStress()
biax.goal=(-10,-10,0)
biax.doneHook='term()'
def term1(): # delete a determined percent of fines after consolidation and
then reconsolidation
bodyRadius=[]
for b in O.bodies:
if b.shape.radius<=0.005: # define fine particles
#if isinstance(b.shape,Sphere):
bodyRadius.append([b,b.shape.radius])
bodyRadius.sort(key=lambda x:x[1])
l=len(bodyRadius) # how many fine particles
perc=0.3
delNum=int(l*perc) # how many number to be deleted
list=random.sample(bodyRadius,delNum) # list to be deleted
for b in list:
