After deletion of particles, why the output of particle number remains unchanged

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=10000),
   #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:
    O.bodies.erase(b[0].id)

   # output PSD
   psd = utils.psd(bins=10,mass=True,mask=-1)
   print psd[0],psd[1] # lists of bins' sizes; cumulative percent of material
   plt.semilogx(psd[0],psd[1])
   #plt.show()
   plt.savefig('psd.png')

   fpsd=file('psd.dat','a')
   fpsd.write('psd[0] psd[1]\n')
   fpsd.write(str(psd[0])+' '+str(psd[1])+'\n')
   fpsd.close()

   biax.doneHook='term()'
"""

def coh():
   O.engines[2].physDispatcher.functors[0].setCohesionNow=True
   print 'add cohesion'
   biax.doneHook='shear()'

# for biaxial shear below:
def shear():
   print getStress()
   #setContactFriction(0.5)
   # set the current cell configuration to be the reference one
   O.cell.trsf=Matrix3.Identity
   biax.goal=(sigmaIso,-0.1,0)
   biax.stressMask=1
   # strain rate along y-axis is 0.01, use a larger x-axis strain rate to better maintian stresses
   biax.maxStrainRate=(0.5,0.01,0)
   biax.relStressTol=.01
   biax.maxUnbalanced=.01
   print 'shearing'
   biax.doneHook='biaxFinished()'

def biaxFinished():
   print 'Biaxial Test Finished'
   O.pause()

"""
def term():
   O.engines = O.engines[:3]+O.engines[4:]
   print len(O.bodies)
   print getStress()
   print O.cell.hSize
   setContactFriction(0.5)
   O.cell.trsf=Matrix3.Identity
   O.cell.velGrad=Matrix3.Zero
   for p in O.bodies:
      p.state.vel = Vector3.Zero
      p.state.angVel = Vector3.Zero
      p.state.refPos = p.state.pos
      p.state.refOri = p.state.ori
   O.save('0.yade.gz')
   O.pause()
"""

O.run(1000)
O.wait()