Periodic random packing of spheres in a cube

Hello,

pack.randomPeriPack() function could be useful for you. It is described here:
https://yade-dem.org/doc/yade.pack.html?highlight=randomperipack#yade.pack.randomPeriPack

Hi!

Yes, i tried that, but am not sure how to export the radius and position of the spheres. This is what i tried in command terminal:

 yade.pack.randomPeriPack(2, [ 20,20,20], 0.3, memoizeDb=None, noPrint = False)
 -> [1]: <yade._packSpheres.SpherePack at 0x7f2ebe257368>

So it says it returned SpherePack, how should i export the position and radius of all the spheres?
-Shailesh

Indeed this returned SpherePack is not directly used in your simulation. But, once you're used to, it is not so difficult to do so, thanks to "toSimulation() function (see https://yade-dem.org/doc/yade.pack.html#yade._packSpheres.SpherePack.toSimulation)

Finally all may fit in one line:

pack.randomPeriPack(..........).toSimulation()

Yes! I did that, in the controller and primary view window, i can see the packing generated. How to export the co-ordinate position and radius of all the spheres in the packing including those on the edges? I am not sure how to apply this command:

yade.export.textExt(konda1,format='x_y_z_r',comment=" ",mask=-1,attrs=[])
error: Attribute error

Is there something am missing in the format?

Hi Shailesh,

>
> yade.export.textExt(konda1,format='x_y_z_r',comment=" ",mask=-1,attrs=[])
> error: Attribute error
>

In the case of similar errors, please always give us more details, i.e.
full error message and also some info about your variables (like what
'konda1' is). It will help us solve your problem more quickly :-)

Concerning periodic boundary conditions, you may also call it periodic
contact detection. There is a nice section in 'DEM formulation' section in
docs [1]. All particles are in the simulation physically just once, but
collisions are detected also with periodic images of other particles. Thus
there is nothing like edge particles and you would need some manual
postprocessing if you need such information.

Concerning SpherePack, you can use toList or other approaches to extract
info into a list and then do whatever you want with the values, see [2].
Just to complete the toSimulation answer.

Cheers
Jan

[1]
https://yade-dem.org/doc/formulation.html#collision-detection-in-periodic-cell
[2] https://yade-dem.org/doc/yade.pack.html#yade._packSpheres.SpherePack
Dne 7.10.2014 5:51 "skonda2" <email address hidden>
napsal(a):

> Question #255437 on Yade changed:
> https://answers.launchpad.net/yade/+question/255437
>
> Status: Answered => Open
>
> skonda2 is still having a problem:
> Yes! I did that, in the controller and primary view window, i can see
> the packing generated. How to export the co-ordinate position and radius
> of all the spheres in the packing including those on the edges? I am not
> sure how to apply this command:
>
> yade.export.textExt(konda1,format='x_y_z_r',comment=" ",mask=-1,attrs=[])
> error: Attribute error
>
> Is there something am missing in the format?
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Ok, I did not get that you wanted to export in a text file position and radii of your existing sample.. So

- if you want only position and radii, export.text() function is enough, you do not need export.textExt().
And by the way, all arguments, with "=" sign in the documentation (as "format=.." , "comment=..") are optional arguments, with their default values appearing after the equal. So it is useless to define them, if you keep the default values

- whatever the function you use, the first argument has to be a string. So, unless you defined previously "konda1" as a string before (and in this case, I do not see where the error would come from), you have to type

export.text('konda1') (note the quotes that define this argument as a string !)

Hi!

Jan,
Documentation for Periodic packing was very helpful, i did go through it completely.
Jerome,
I implemented what you suggested.

So i wrote a script of what am trying to do:

from yade import pack, export
yade.pack.randomPeriPack(2, [ 25,25,25], 0.3, noPrint =True).toSimulation()
yade.qt.Controller(), yade.qt.View()
yade.export.text("konda1")

What i got is a reasonable packing structure with around 108 spheres and also i was able to export the data in the file "konda1". The problem is, i do not see any spheres on the boundaries the periodic images in the primary view window.
However, when i check the periodic cell in the inspect window or click the play button, am able to see the (extra) spheres on the boundaries with their periodic images on the opposite sides (wire frame view for spheres near boundaries). That is the sort of packing am trying to get.
I am not sure how to export the position and radii for these set of spheres.? I feel am very close to the final objective which is get to periodic dense packing of non-overlapping spheres in a cube.

-Shailesh

Hi Shailesh

> What i got is a reasonable packing structure with around 108 spheres and
> also i was able to export the data in the file "konda1". The problem is, i
> do not see any spheres on the boundaries the periodic images in the primary
> view window.

However, when i check the periodic cell in the inspect window or click the
> play button, am able to see the (extra) spheres on the boundaries with
> their periodic images on the opposite sides (wire frame view for spheres
> near boundaries). That is the sort of packing am trying to get.
> I am not sure how to export the position and radii for these set of
> spheres.? I feel am very close to the final objective which is get to
> periodic dense packing of non-overlapping spheres in a cube.
>

the visualization of "edge particles" is some kind of postprocessing. The
particles, whose bounding boxes are close to the edge (I am not sure about
actual implementation) of periodic cell are rendered also as the periodic
images.

If you want to export this, you can try something like (example of images
in x direction):

from yade import pack
pack.randomPeriPack(1,(10,10,10)).toSimulation()
O.step() # to initialize bounding boxes
xImages = [] # list of "edge" particles
for b in O.bodies:
xmin,xmax = O.cell.wrap(b.bound.min)[0], O.cell.wrap(b.bound.max)[0] # wrap
is important as physically the particles can be anywhere
if xmin > xmax: # this means that bounding box crosses pariodic cell. You
can define various different conditions..
xImages.append(b)
pr = [[O.cell.wrap(b.state.pos),b.shape.radius] for b in xImages] # list or
[pos,radius] of "edge" particles. Now pos is always inside the cell
for i,(pos,r) in enumerate(pr):
shift = Vector3(O.cell.size[0],0,0) * (1 if pos[0]<.5*O.cell.size[0] else
-1) # determine shift to create a periodic image. +x direction if particles
is near -x face and vice versa
pr[i][0] += shift
# saves images into a file
f = open("/tmp/ximages.dat","w")
f.write("x y z r\n")
f.writelines("%g %g %g %g\n"%(pos[0],pos[1],pos[2],r) for pos,r in pr)
f.close()

I don't know how well you know Python, so if you found something unclear
(code itself or its meaning), just ask :-)

cheers
Jan

Hi Jan,

I tried to run the script, but it generated empty output file with just the first line as " x y z "
I have attached the script that i ran, is there something wrong in that?

# Script for exporting data of all the spheres, including the one on the edges and periodic images
from yade import pack
pack.randomPeriPack(1,(10,10,10)).toSimulation()
O.step()
xImages = []
for b in O.bodies:
 xmin,xmax = O.cell.wrap(b.bound.min)[0], O.cell.wrap(b.bound.max)[0]

if xmin>xmax:
 xImages.append(b)

pr = [[O.cell.wrap(b.state.pos),b.shape.radius] for b in xImages]

for i, (pos,r) in enumerate(pr):
 shift = Vector3(O.cell.size[0],0,0)*(1 if pos[0]<.5*O.cell.size[0] else -1)
 pr[i][0] += shift

f = open("/home/konda/ximages_1.dat","w")
f.write("x y z r\n")
f.writelines("%g %g %g %g\n" % (pos[0],pos[1],pos[2],r) for pos, r in pr)
f.close()

So, " for b in O.bodies" does this also consider the bodies on the edges and their periodic images too? Because before intializing the box, it had say 72 bodies and after initialising )O.step) there should be more bodies including those on the edges. So i tried len(O.bodies), which came out as 72. I was thinking, may be in the for loop it does not consider any bodies on the edges, so there is no output. Let me kno what you think of my understanding.

-Shailesh

Hi Shailesh,

>
> I tried to run the script, but it generated empty output file with just
> the first line as " x y z "
> I have attached the script that i ran, is there something wrong in that?
>

the indentation is wrong. The indents are extremely important in Python :-)
please try to use the indents I did in my previous email. If you have still
problems, I can send you the script directly (outside this forum, as
attaching files is not possible).

>
> So, " for b in O.bodies" does this also consider the bodies on the edges
> and their periodic images too? Because before intializing the box, it had
> say 72 bodies and after initialising )O.step) there should be more bodies
> including those on the edges. So i tried len(O.bodies), which came out as
> 72. I was thinking, may be in the for loop it does not consider any bodies
> on the edges, so there is no output. Let me kno what you think of my
> understanding.
>
>
O.bodies has always the same length. Let me remind you once more that there
are *no* "bodies on the edges" and every particle is present only once in
the simulation. It might not be easy to understand at first, especially if
you have some background in FEM-like periodic simulations (I remember my
problems). The periodicity is stored "merely" in the interactions,
remembering that the second particle of the interaction is xyz periods form
the first one. Nothing else. So, O.bodies does not contain periodic images.

You can see some "ghost" particles in 3D window, but it is only a
postprocessing and visualization. The program determines which particles
are sufficiently close to the boundary and draw also their periodic images.
But from simulation point of view, number of particles is always the same.

HTH
Jan

Hi Jan,

There was no indentations in your previous reply, i guess in this forum it is difficult to do it. If possible, you can possibly mail me the script at "<email address hidden>".
I did run the script with proper indentations and got some output. So, i believe the pos and radius that is generated is of the spheres that are on the edges with their periodic images.
I generated script for y direction and have attached the script below, do you think there is anything wrong?

#script for edge particles in y-direction
from yade import pack
pack.randomPeriPack(1,(10,10,10)).toSimulation()
O.step() # to initialize bounding boxes
yImages = [] # list of "edge" particles
for b in O.bodies:
 ymin,ymax = O.cell.wrap(b.bound.min)[1], O.cell.wrap(b.bound.max)[1] # wrap is important as physically the particles can be anywhere
 if ymin > ymax: # this means that bounding box crosses pariodic cell. You can define various different conditions..
  yImages.append(b)
pr = [[O.cell.wrap(b.state.pos),b.shape.radius] for b in yImages] # list or [pos,radius] of "edge" particles. Now pos is always inside the cell

for i,(pos,r) in enumerate(pr):
 shift = Vector3(0,O.cell.size[1],0) * (1 if pos[1]<.5*O.cell.size[1] else -1) # determine shift to create a periodic image. +y... direction if particles is near -y face and vice versa
 pr[i][0] += shift

# saves images into a file
f = open("/home/konda/yimages_1.dat","w")
f.write("x y z r\n")
f.writelines("%g %g %g %g\n"%(pos[0],pos[1],pos[2],r) for pos,r in pr)
f.close()

Similarly, i generated script for z direction. So now, i have position and radius of all spheres near the boundaries in x,y,z directions. Now, i also exported the pos n radius of spheres which are already in simulation. So in total, there are pos and radius of all the "spheres in simulation" and "ghost" spheres.
I tried to visualize the whole packing, and found that there are some overlapping spheres, may there is something wrong with the program am using for y and z directions. Let me know, what you think.

I guess am beginning to understand the periodic conditions, being from FEM background i was under the impression that periodic random packing would mean it would generate packing including those close to the edges with their periodic images. I was also working on molecular dynamics based hard sphere packing, which includes those spheres on the edges and the copy of them. Thanks for your detailed explanation.

So, in "yade" post- processing we see those "ghost" particles which is the sort of packing am looking for and very similar to what we have it in FEM and molecular dynamics simulation. Just to put in perspective, am trying to carry out high level CFD calculations on that sort of packing structures. So for now, my objective is to get packing structure (including "ghost" spheres) in yade and export there position and radius.

-Shailesh

Hi Shailesh,

> There was no indentations in your previous reply, i guess in this forum it
> is difficult to do it. If possible, you can possibly mail me the script at "
> <email address hidden>".
> I did run the script with proper indentations and got some output. So, i
> believe the pos and radius that is generated is of the spheres that are on
> the edges with their periodic images.
> I generated script for y direction and have attached the script below, do
> you think there is anything wrong?
>

I believe it is correct. Manipulating with positions is done only in term
of periodic cell size, so it should be ok..

> I tried to visualize the whole packing, and found that there are some
> overlapping spheres, may there is something wrong with the program am using
> for y and z directions. Let me know, what you think.
>

how much do they overlap? Could you compare your vizualization with that in
yade? Both cases should look the same..

cheers
Jan

Hi Jan,

I was making an error while running program. I figured that out. Now, i do not see any overlapping spheres, but the packing visualization is not the same.
So, what i do right now is, export all the x,y,z directions periodic image spheres ("ghost" spheres) from the program that you provided and then i export all the spheres from the simulation, then i see the visualization.
I think, when i export the spheres from the simulation, there are already some spheres which are on the edges of the periodic cell (position beyond or on periodic cell) DO you think that could be one reason?
The periodic cell size is 7.86 for the above simulation, so when i expor the spheres from simulation, do you think i should only consider spheres which are completely inside the periodic cell?, since we already have spheres from our program.

-Shailesh

Hi!

I forgot to add script in my last comment, so the script am using is below:

from yade import pack
pack.randomPeriPack(1,(10,10,10)).toSimulation()
O.step() # to initialize bounding boxes
xImages = [] # list of "edge" particles

for b in O.bodies:
 xmin,xmax = O.cell.wrap(b.bound.min)[0], O.cell.wrap(b.bound.max)[0] # wrap is important as physically the particles can be anywhere

 if xmin > xmax: # this means that bounding box crosses periodic cell. You can define various different conditions..
  xImages.append(b)

pr = [[O.cell.wrap(b.state.pos),b.shape.radius] for b in xImages] # list or [pos,radius] of "edge" particles. Now pos is always inside the cell

for i,(pos,r) in enumerate(pr):
 shift = Vector3(O.cell.size[0],0,0) * (1 if pos[0]<0.5*O.cell.size[0] else -1) # determine shift to create a periodic image +x... direction if particles is near -x face and vice versa
 pr[i][0] += shift

# saves images into a file
f = open("/home/konda/imp/1_1_x_1.dat","w")
f.write("x y z r\n")
f.writelines("%g %g %g %g\n"%(pos[0],pos[1],pos[2],r) for pos,r in pr)
f.close()

yImages = []

for b in O.bodies:
 ymin,ymax = O.cell.wrap(b.bound.min)[1], O.cell.wrap(b.bound.max)[1]
 if ymin > ymax:
  yImages.append(b)

pr1 = [[O.cell.wrap(b.state.pos),b.shape.radius] for b in yImages]

for i,(pos,r) in enumerate(pr1):
 shift = Vector3(0,O.cell.size[1],0) * (1 if pos[1]<0.5*O.cell.size[1] else -1)
 pr1[i][0] += shift

# saves images into a file
f = open("/home/konda/imp/1_1_y_1.dat","w")
f.write("x y z r\n")
f.writelines("%g %g %g %g\n"%(pos[0],pos[1],pos[2],r) for pos,r in pr1)
f.close()

zImages = []

for b in O.bodies:
 zmin,zmax = O.cell.wrap(b.bound.min)[2], O.cell.wrap(b.bound.max)[2]
 if zmin > zmax:
  zImages.append(b)

pr2 = [[O.cell.wrap(b.state.pos),b.shape.radius] for b in zImages]

for i,(pos,r) in enumerate(pr2):
 shift = Vector3(0,0,O.cell.size[2]) * (1 if pos[2]<0.5*O.cell.size[2] else -1)
 pr2[i][0] += shift

# saves images into a file
f = open("/home/konda/imp/1_1_z_1.dat","w")
f.write("x y z r\n")
f.writelines("%g %g %g %g\n"%(pos[0],pos[1],pos[2],r) for pos,r in pr2)
f.close()

So this program generates three files, for the periodic spheres in three directions. Is there anything wrong with script am using?

-SHailesh

Hi

2014-10-13 18:01 GMT+02:00 skonda2 <email address hidden>:

> Question #255437 on Yade changed:
> https://answers.launchpad.net/yade/+question/255437
>
> Status: Answered => Open
>
> skonda2 is still having a problem:
> Hi Jan,
>
> I was making an error while running program. I figured that out. Now, i do
> not see any overlapping spheres, but the packing visualization is not the
> same.
>

do you compare the same view and the same run (since every realization is
different)?

> So, what i do right now is, export all the x,y,z directions periodic image
> spheres ("ghost" spheres) from the program that you provided and then i
> export all the spheres from the simulation, then i see the visualization.
> I think, when i export the spheres from the simulation, there are already
> some spheres which are on the edges of the periodic cell (position beyond
> or on periodic cell) DO you think that could be one reason?
>

It is also a feature of periodic cell in Yade. The spheres may have
arbitrary position :-) specifically in the case of your current approach,
all spheres starts "inside" periodic cell, i.e. with zero period shift, but
during compaction stage (in randomPeriPack), it may happen that spheres
move to another ("neigboring") period. But the yade visualization of them
are always inside the cell. So to export "normal" particles, you should do:

f = open('normal.dat')
f.write('x y z r\n')
for b in bodies:
   p,r = O.cell.wrap(b.state.pos),b.shape.radius
   f.write('%g %g %g %g\n'%(p[0],p[1],p[2],r))
f.close()

instead of the export function you use now. It is another tricky step,
being difficult for beginners :-)

> The periodic cell size is 7.86 for the above simulation, so when i expor
> the spheres from simulation, do you think i should only consider spheres
> which are completely inside the periodic cell?, since we already have
> spheres from our program.
>

no, export all spheres from O.bodies ("normal") + ghosts, see above.

cheers
Jan

>
>
>
> I forgot to add script in my last comment, so the script am using is
> below:
>

I will have a look later, sorry..

>
> So this program generates three files, for the periodic spheres in three
> directions. Is there anything wrong with script am using?
>

Initially I showed you only x direction as an example, I think you can then
merge all directions into one single file. It is also possible, that you
would have to specially treat corners :-)

I know that the beginning of working with periodic yade is not easy and
needs some time, so do not feel upset not understanding everything
imediatelly :-)

cheers
Jan

PS: a mind training: generate your periodic packing, and than run:

import random
ri = random.randint

calm() # sets zero velocity to all particles
for b in O.bodies:
   x,y,z = ri(-5,5), ri(-5,5), ri(-5,5) # some random integers
   shift = Vector3(O.cell.size[0]*x, O.cell.size[1]*y, O.cell.size[2]*z) #
shift of particles to be applied, as a random multiplication of cell size
   vel = shift/O.dt # transfer shift to velocity as Yade don't like direct
position assignment
   b.state.vel = vel
O.step() # apply new positions
calm() # sets zero velocities

now your packing should look almost the same as before these steps, but the
actual positions of particles are complete mess for the first sight (if you
print or export it)

Hi Jan,

I guess exporting "normal" spheres did the trick. So, i will post the final script that am using for periodic packing, which can be useful for any future yade users in this forum.

from yade import pack
pack.randomPeriPack(3.5,(45,45,45), rRelFuzz = 0.5).toSimulation()
O.step() # to initialize bounding boxes
xImages = [] # list of "edge" particles

for b in O.bodies:
 xmin,xmax = O.cell.wrap(b.bound.min)[0], O.cell.wrap(b.bound.max)[0] # wrap is important as physically the particles can be anywhere

 if xmin > xmax: # this means that bounding box crosses periodic cell. You can define various different conditions..
  xImages.append(b)

pr = [[O.cell.wrap(b.state.pos),b.shape.radius] for b in xImages] # list or [pos,radius] of "edge" particles. Now pos is always inside the cell

for i,(pos,r) in enumerate(pr):
 shift = Vector3(O.cell.size[0],0,0) * (1 if pos[0]<0.5*O.cell.size[0] else -1) # determine shift to create a periodic image +x... direction if particles is near -x face and vice versa
 pr[i][0] += shift

# saves images into a file
f = open("/home/konda/imp/Extra_imp/one_packing/1_1_x.dat","w")
f.write("x y z r\n")
f.writelines("%g %g %g %g\n"%(pos[0],pos[1],pos[2],r) for pos,r in pr)
f.close()

yImages = []

for b in O.bodies:
 ymin,ymax = O.cell.wrap(b.bound.min)[1], O.cell.wrap(b.bound.max)[1]
 if ymin > ymax:
  yImages.append(b)

pr1 = [[O.cell.wrap(b.state.pos),b.shape.radius] for b in yImages]

for i,(pos,r) in enumerate(pr1):
 shift = Vector3(0,O.cell.size[1],0) * (1 if pos[1]<0.5*O.cell.size[1] else -1)
 pr1[i][0] += shift

# saves images into a file
f = open("/home/konda/imp/Extra_imp/one_packing/1_1_y.dat","w")
f.write("x y z r\n")
f.writelines("%g %g %g %g\n"%(pos[0],pos[1],pos[2],r) for pos,r in pr1)
f.close()

zImages = []

for b in O.bodies:
 zmin,zmax = O.cell.wrap(b.bound.min)[2], O.cell.wrap(b.bound.max)[2]
 if zmin > zmax:
  zImages.append(b)

pr2 = [[O.cell.wrap(b.state.pos),b.shape.radius] for b in zImages]

for i,(pos,r) in enumerate(pr2):
 shift = Vector3(0,0,O.cell.size[2]) * (1 if pos[2]<0.5*O.cell.size[2] else -1)
 pr2[i][0] += shift

# saves images into a file
f = open("/home/konda/imp/Extra_imp/one_packing/1_1_z.dat","w")
f.write("x y z r\n")
f.writelines("%g %g %g %g\n"%(pos[0],pos[1],pos[2],r) for pos,r in pr2)
f.close()

# saving normal spheres
f = open("/home/konda/imp/Extra_imp/one_packing/normal.dat", "w")
f.write('x y z r\n')
for b in O.bodies:
 p,r = O.cell.wrap(b.state.pos),b.shape.radius
 f.write('%g %g %g %g\n' %(p[0],p[1],p[2],r))
f.close()

So, this script produces periodic packing of spheres in x,y,z directions. I tested the visualization, it is very much similar to what we observe in yade after exporting the position. There is something wrong with the corner spheres. I am not sure how to export positions of corner spheres and there periodic images? It is very hard to understand, how to incorporate those in my script?

Thanks Jan Stránský, that solved my question.