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

Asked by zheng on 2019-08-27

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()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Robert Caulk
Solved:
2019-09-17
Last query:
2019-09-17
Last reply:
2019-08-28
Robert Caulk (rcaulk) said : #1

Yikes, that is the longest way to show:

O.bodies.append(sphere(Vector3(0,0,0),1))
O.run(1,1)
len(O.bodies)
O.bodies.erase(0)
O.run(1,1)
len(O.bodies)

That I have ever seen.

It is because the body id can be reused if the particle is reinserted.

Best Robert Caulk (rcaulk) said : #2

If you iterate on the bodies, it will not count erased bodies:

count = 0
for b in O.bodies:
    count+=1

zheng (tigerhz) said : #3

Thanks Robert Caulk, that solved my question.