How to calculate the height of particle accumulation?

Hi,

Function aabbExtrema() returns boundary box of particle packing (min and max corner positions).

If you are only interested in dimension, you can use aabbDim().

Both functions work for spherical bodies.

Cheers
Karol

[1] https://yade-dem.org/doc/yade.utils.html?highlight=aabbextrema#yade._utils.aabbExtrema
[2] https://yade-dem.org/doc/yade.utils.html?highlight=aabbdim#yade.utils.aabbDim

And one more comment
Are you sure that you generate particles with this function geom.facetbox ((,), (,), wallmask=)?

You generate bodies (facets), but I wouldn't call them particles by default.

Cheers,
Karol

My code for generating particles is as follows:
O.bodies.append(geom.facetBox((0.05, 3, 4), (0.1, 3, 4), wallMask=63))
sp = pack.SpherePack()
sp.makeCloud(Vector3(0.05,0,0),Vector3(0.05,6,8), rMean=0.05, rRelFuzz=0)
sp.toSimulation()
What I understand is that o.bodies.append (geom. Facetbox ((,), (,), wallmask =) is used to create faces and generate six faces, that is, a cuboid or other shaped objects with a certain size. Then I can fill them with sp.makecloud. I don't know whether this understanding is correct.

Hi,

Yes, first you generate facets, next particles. I only meant that facetbox() function creates facets, not particles (as stated in the first post).

Referring to your question, this is how you can obtain the height of spherical particle packing:
#################
O.bodies.append(geom.facetBox((0.05, 3, 4), (0.1, 3, 4), wallMask=63))
sp = pack.SpherePack()
sp.makeCloud(Vector3(0.05,0,0),Vector3(0.05,6,8), rMean=0.05, rRelFuzz=0)
sp.toSimulation()
O.engines[4].gravity = (0,0,-10)
O.run(3000,True)

(xmin,ymin,zmin),(xmax,ymax,zmax) = aabbExtrema() # you can use aabbDim() instead in those two lines if you only need height
height = zmax-zmin
print("Height is ",height)
################

Cheers,
Karol

For example:

O.bodies.append(geom.facetBox((0.05, 3, 4), (0.1, 3, 4), wallMask=63))
sp = pack.SpherePack()
sp.makeCloud(Vector3(0.05,0,0),Vector3(0.05,6,8), rMean=0.05, rRelFuzz=0)
sp.toSimulation()

(xmin,ymin,zmin),(xmax,ymax,zmax) = aabbDim()
height = zmax-zmin
print("Height is ",height)

Is that okay? Is there something wrong with the code?

aabbDim() gives you dimensions (three values), so it should rather be:

#####
O.bodies.append(geom.facetBox((0.05, 3, 4), (0.1, 3, 4), wallMask=63))
sp = pack.SpherePack()
sp.makeCloud(Vector3(0.05,0,0),Vector3(0.05,6,8), rMean=0.05, rRelFuzz=0)
sp.toSimulation()

(xdim,ydim,zdim) = aabbDim()

print("Height is ",zdim)
####

Cheers,
Karol

I just tried it and I found that the code works, but there is no zdim output in the console, why is that?

Are you sure? What is the console output after running this code?

I tried the code again and there was no output from zdim.
The code is as follows:

from yade import pack,plot
O.bodies.append(geom.facetBox((0.05, 9, 10), (0.1, 9, 10), wallMask=63))
sp = pack.SpherePack()
sp.makeCloud(Vector3(0.05,0,0),Vector3(0.05,18,20), rMean=0.1, rRelFuzz=1/3)
sp.toSimulation()

(xdim,ydim,zdim)= aabbDim()

print("Height is ",zdim)

for b in O.bodies:
   if isinstance(b.shape,Sphere):
       b.state.blockedDOFs='ZxY'
       b.shape.color=(3.,2.,1.)

circleRadius=1.5
circleCenter = Vector3(0.05,6,6)
circleRadius1=1.5
circleCenter1 = Vector3(0.05,12,6)
   #
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
        InteractionLoop(
                # handle sphere+sphere and facet+sphere collisions
                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.4),
        # call the checkUnbalanced function (defined below) every 2 seconds
        PyRunner(command='checkUnbalanced()', realPeriod=2),
        # call the addPlotData function every 200 steps
        PyRunner(command='addPlotData()', iterPeriod=100)
]
O.dt = 0.5 * PWaveTimeStep()

# enable energy tracking; any simulation parts supporting it
# can create and update arbitrary energy types, which can be
# accessed as O.energy['energyName'] subsequently
O.trackEnergy = True

# if the unbalanced forces goes below .05, the packing
# is considered stabilized, therefore we stop collected
# data history and stop
def checkUnbalanced():
 if unbalancedForce() < .05:
  O.pause()
  plot.saveDataTxt('bbb.txt.bz2')
                for b in O.bodies:
      if isinstance(b.shape,Sphere):
                        #b.state.blockedDOFs='zxy'
                        b.state.vel=(0,0,0)
                        b.state.angVel=(0,0,0)
                    d = (b.state.pos - circleCenter).norm() # distance of circleCenter and center of "b"
                    d1 = (b.state.pos - circleCenter1).norm() # distance of circleCenter and center of "b"
                    if d < circleRadius:
                       O.bodies.erase(b.id)
                    if d1 < circleRadius1:
                       O.bodies.erase(b.id)

  O.bodies.append(sphere(center=(0.05,7.378730854,6.1431074487337), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.315569669,6.37882820065594), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.212435565,6.6), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.07246222,6.79990265356116), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.899902654,6.97246222036657), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.7,7.11243556529821), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.478828201,7.21556966910027), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.243107449,7.27873085421709), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6,7.3), radius=.1274, fixed=True,color=(1.,1.,1.)))

  O.bodies.append(sphere(center=(0.05,5.756892551,7.27873085421709), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.521171799,7.21556966910027), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.3,7.11243556529822), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.100097346,6.97246222036657), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.92753778,6.79990265356116), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.787564435,6.6), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.684430331,6.37882820065594), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.621269146,6.1431074487337), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.6,5.9), radius=.1274, fixed=True,color=(1.,1.,1.)))

  O.bodies.append(sphere(center=(0.05,4.621269146,5.6568925512663), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.684430331,5.42117179934406), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.787564435,5.2), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.92753778,5.00009734643885), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.100097346,4.82753777963343), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.3,4.68756443470179), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.521171799,4.58443033089973), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.756892551,4.52126914578291), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6,4.5), radius=.1274, fixed=True,color=(1.,1.,1.)))

  O.bodies.append(sphere(center=(0.05,6.243107449,4.52126914578291), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.478828201,4.58443033089973), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.7,4.68756443470179), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.899902654,4.82753777963343), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.07246222,5.00009734643885), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.212435565,5.2), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.315569669,5.42117179934406), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.378730854,5.6568925512663), radius=.1274, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.4,5.9), radius=.1274, fixed=True,color=(1.,1.,1.)))

# collect history of data which will be plotted
def addPlotData():
 # each item is given a names, by which it can be the unsed in plot.plots
 # the **O.energy converts dictionary-like O.energy to plot.addData arguments
 plot.addData(i=O.iter, unbalanced=unbalancedForce(), **O.energy)

# define how to plot data: 'i' (step number) on the x-axis, unbalanced force
# on the left y-axis, all energies on the right y-axis
# (O.energy.keys is function which will be called to get all defined energies)
# None separates left and right y-axis

plot.plots = {'i': ('unbalanced', None, O.energy.keys)}

# show the plot on the screen, and update while the simulation runs
plot.plot()

O.saveTmp()

# to see it
from yade import qt
qt.Controller()
qt.View()
qt.View()

Yes, I added these two lines of code again, and zdim got the output. What does zdim mean? Is it the highest point? Nadir? Or is it an average?

If you use

####
(xmin,ymin,zmin),(xmax,ymax,zmax) = aabbExtrema()
####

zmax is the highest point.

If you use

####
(xdim,ydim,zdim)= aabbDim()
####

zdim is packing height (the highest point minus the lowest)

Karol

Thanks Karol Brzezinski, that solved my question.