Particle size distribution curve after particle breaking

Hello,

I guess you want something like this:

def getRadiusArray():
    radii=[]
    for b in O.bodies:
 if isinstance(b.shape,Sphere):
         radii.append(b.shape.radius)

getRadiusArray()
values, base = np.histogram(data, bins=40)
cumulative = np.array(np.cumsum(values),dtype=float)
perc = cumulative/max(cumulative)*100
plt.semilogx(base[:-1], cumulative, c='blue')

Cheers,

Robert

**noticed it was missing some stuff**

def getRadiusArray():
....

    return radii

radii = getRadiusArray()
values,base=np.histogram(radii,bins=40)

Hi Robert:

Thank you for your reply !! But I meet some problems:
I added the code you provided to mine, there seems to be some conflict between plt.show() and plot.plot(), Or my understanding is wrong. Help me plz..

This is my whole code:

readParamsFromTable(rMean=.05,rRelFuzz=.3,maxLoad=1e6,minLoad=1e4)
from yade.params.table import *

from yade import pack, plot
O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=rMean)
sp.toSimulation()

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),
   PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
   PyRunner(command='replace()',iterPeriod=100,label ='rep'),
   #PyRunner(command='getRadiusArray()',iterPeriod=100,label ='psd'),

]
O.dt=.5*PWaveTimeStep()
##############################################################
def selectBodiesToCrush():
    stresses = bodyStressTensors() # e.g. based on overall particle stress
    ret = []
    for i,stress in enumerate(stresses):
        if stress[0,0] > 50: # Particle breakage criterion
           ret.append(i)
    return ret

def createChildren(b):
    p = b.state.pos
    r = b.shape.radius
    # Particle replacement mode
    s1 = sphere(p+(.5*r,0,0),.5*r)
    s2 = sphere(p-(.5*r,0,0),.5*r)
    return s1,s2

def replace():
    bodiesToDelete = selectBodiesToCrush()
    for i in bodiesToDelete:
        if isinstance(O.bodies[i].shape,Sphere): # Limitation of the smallest particles that can break
            b = O.bodies[i]
            children = createChildren(b)
            O.bodies.erase(i)
            O.bodies.append(children)

##############################################################
#import matplotlib.pyplot as plt
#import numpy as np
#def getRadiusArray():
# radii=[]
# for b in O.bodies:
# if isinstance(b.shape,Sphere):
# radii.append(b.shape.radius)
# return radii
#
#radii=getRadiusArray()
#values, base = np.histogram(radii, bins=40)
#cumulative = np.array(np.cumsum(values),dtype=float)
#perc = cumulative/max(cumulative)*100
#plt.semilogx(base[:-1], cumulative, c='blue')
#plt.show()
##############################################################

def checkUnbalanced():
   if O.iter<5000: return
   if unbalancedForce()>.1: return
   O.bodies.append(wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1))
   global plate
   plate=O.bodies[-1]
   plate.state.vel=(0,0,-.1)
   O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=200)]
   checker.command='unloadPlate()'

def unloadPlate():
   if abs(O.forces.f(plate.id)[2])>maxLoad:
      plate.state.vel*=-1
      checker.command='stopUnloading()'

def stopUnloading():
   if abs(O.forces.f(plate.id)[2])<minLoad:
      plot.saveDataTxt(O.tags['d.id']+'.txt')
      O.pause()

def addPlotData():
   if not isinstance(O.bodies[-1].shape,Wall):
      plot.addData(); return
   Fz=O.forces.f(plate.id)[2]
   plot.addData(Fz=Fz,w=plate.state.pos[2]-plate.state.refPos[2],unbalanced=unbalancedForce(),i=O.iter)

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

O.run()
#from yade import utils
#binsSizes, binsProc, binsSumCum= utils.psd(bins=5, mass=True)

waitIfBatch()

Thanks Robert Caulk, that solved my question.

I'm sorry, I disagree with your choice to edit the thread after you received help from the community. It may come as a surprise to you, but others will want to use this thread in the future. Therefore, I am going to go ahead and re-post the response that you removed:

>>>>>>>>>>>>>>>>>>>>>
Hi Robert:

Thank you for your reply !! But I meet some problems:
I added the code you provided to mine, there seems to be some conflict between plt.show() and plot.plot(), Or my understanding is wrong. Help me plz..

This is my whole code:

readParamsFromTable(rMean=.05,rRelFuzz=.3,maxLoad=1e6,minLoad=1e4)
from yade.params.table import *

from yade import pack, plot
O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=rMean)
sp.toSimulation()

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),
   PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
   PyRunner(command='replace()',iterPeriod=100,label ='rep'),
   #PyRunner(command='getRadiusArray()',iterPeriod=100,label ='psd'),

]
O.dt=.5*PWaveTimeStep()
##############################################################
def selectBodiesToCrush():
    stresses = bodyStressTensors() # e.g. based on overall particle stress
    ret = []
    for i,stress in enumerate(stresses):
        if stress[0,0] > 50: # Particle breakage criterion
           ret.append(i)
    return ret

def createChildren(b):
    p = b.state.pos
    r = b.shape.radius
    # Particle replacement mode
    s1 = sphere(p+(.5*r,0,0),.5*r)
    s2 = sphere(p-(.5*r,0,0),.5*r)
    return s1,s2

def replace():
    bodiesToDelete = selectBodiesToCrush()
    for i in bodiesToDelete:
        if isinstance(O.bodies[i].shape,Sphere): # Limitation of the smallest particles that can break
            b = O.bodies[i]
            children = createChildren(b)
            O.bodies.erase(i)
            O.bodies.append(children)

##############################################################
#import matplotlib.pyplot as plt
#import numpy as np
#def getRadiusArray():
# radii=[]
# for b in O.bodies:
# if isinstance(b.shape,Sphere):
# radii.append(b.shape.radius)
# return radii
#
#radii=getRadiusArray()
#values, base = np.histogram(radii, bins=40)
#cumulative = np.array(np.cumsum(values),dtype=float)
#perc = cumulative/max(cumulative)*100
#plt.semilogx(base[:-1], cumulative, c='blue')
#plt.show()
##############################################################

def checkUnbalanced():
   if O.iter<5000: return
   if unbalancedForce()>.1: return
   O.bodies.append(wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1))
   global plate
   plate=O.bodies[-1]
   plate.state.vel=(0,0,-.1)
   O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=200)]
   checker.command='unloadPlate()'

def unloadPlate():
   if abs(O.forces.f(plate.id)[2])>maxLoad:
      plate.state.vel*=-1
      checker.command='stopUnloading()'

def stopUnloading():
   if abs(O.forces.f(plate.id)[2])<minLoad:
      plot.saveDataTxt(O.tags['d.id']+'.txt')
      O.pause()

def addPlotData():
   if not isinstance(O.bodies[-1].shape,Wall):
      plot.addData(); return
   Fz=O.forces.f(plate.id)[2]
   plot.addData(Fz=Fz,w=plate.state.pos[2]-plate.state.refPos[2],unbalanced=unbalancedForce(),i=O.iter)

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

O.run()
#from yade import utils
#binsSizes, binsProc, binsSumCum= utils.psd(bins=5, mass=True)

waitIfBatch()
>>>>>>>>>>>>>>>>>>>>>

Yes, please, yangjunjie, don't re-write the history of discussions.
I'm also curious why the first part of the question was removed: "This is my last question about particle crushing".
It seems last question on crushing ends up as first question on the additivity of python codes. ;)

Please think twice _before_ posting, not after posting.

Bruno

Sorry, that's my fault.

Turns out, Yade might have a nice built in function for this[1]. I haven't used it myself, but it could be worth trying.

[1]https://yade-dem.org/doc/yade.utils.html?highlight=utils.psd#yade.utils.psd