Yade

Highest clump position in clump cloud

Asked by Sergio

HIGHEST CLUMP IN CLOUD

Hello Yade mates,

I have made 3 clump clouds, at different heights each one. I'd like to know the position of the highest clump of each cloud at each time step.

To achieve that, I've made a Pyrunner command, as you can read below. I intended to do it by writting different loops.
However, when I run the simulation for a while, and check the Ymax values, their values result the same as prior to the simulation run, they don't vary.
I suppose the problem is related to the first Ymax and HBmax values assigned, but I don't know how to solve it.

Thank you for your advices
Sergio

from yade import pack, qt, plot, ymport, export

#////// MATERIALS \\\\\\
MatBal=O.materials.append(FrictMat(
 young=27e9,poisson=0.30,density=5400,
  frictionAngle=0.7,label='balasto'))

#////// BOUNDARIES (bx - htot - bz) \\\\\\

O.bodies.append(utils.geom.facetBox(

 (0,0,0),(1.5,0.8,1.5),wallMask=55))

#////// SPHERE PACK CREATION \\\\\\\\\\\\\\\\

c1=pack.SpherePack([((0.007,0,0),0.022),((-0.011,0,0),0.0205),
 ((0,0.012,0),0.024)

 ])

c2=pack.SpherePack([((0,0,0),0.0122),

 ((0,0.005,0),0.0122),((0.006,0.0075,0),0.0121),

  ((0,0.00875,0),0.0118),
 ])

#////// PROPERTIES \\\\\\

###### Layer 1 ######

capa1=pack.SpherePack()

capa1.makeClumpCloud(

 (-0.1,-0.2,-0.1),(0.1,0,0.2),

 [c1,c2],periodic=False)

capa1.toSimulation(material=MatBal,color=(0.3,0.9,1))

len1=len(O.bodies)

HBmax=len(capa1)-1

Ymax1=O.bodies[HBmax].state.pos[1]

##### Layer 2 ########

capa2=pack.SpherePack()
capa2.makeClumpCloud(

 (-0.1,-0.5,-0.1),(0.1,-0.3,0.2),

 [c1,c2],periodic=False)

capa2.toSimulation(material=MatBal,color=(0.9,0.2,1))
len2=len(O.bodies)
HBmax2=len2
Ymax2=O.bodies[HBmax2-1].state.pos[1]

##### Layer 3 #########

capa3=pack.SpherePack()
capa3.makeClumpCloud(

 (-0.1,-0.8,-0.1),(0.1,-0.6,0.2),

 [c1,c2],periodic=False)
capa3.toSimulation(material=MatBal,color=(0.3,0.2,0.1))

len3=len(O.bodies)

HBmax3=len3
Ymax3=O.bodies[HBmax3-1].state.pos[1]

#///////////// Layer heights \\\\\\\\\\\\\\\

def alturasMaximasCapas():

 HBmax= len1-1

 Ymax1=O.bodies[HBmax].state.pos[1]

 HBmax2=len2-1

 Ymax2=O.bodies[HBmax2].state.pos[1]

 HBmax3=len3-1

 Ymax3=O.bodies[HBmax3].state.pos[1]

# Layer heigth 1

 for ii in range(0,len(capa1)):

  if O.bodies[ii].isClump==True:

   if Ymax1<O.bodies[ii].state.pos[1]:

    Ymax1=O.bodies[ii].state.pos[1]

    HBmax=ii

# Layer heigth 2

 for jj in range(len(capa1),len(capa1)+len(capa2)-1):

  if O.bodies[jj].isClump==True:

   if Ymax2<O.bodies[jj].state.pos[1]:

    Ymax2=O.bodies[jj].state.pos[1]

    HBmax2=jj

# Layer heigth 3

 for zz in range(len(capa2),len(capa1)+len(capa2)+len(capa3)-1):

  if O.bodies[zz].isClump==True:

   if Ymax3<O.bodies[zz].state.pos[1]:

    Ymax3=O.bodies[zz].state.pos[1]

    HBmax3=zz

# ////////////////////// Engines \\\\\\\\\\\\\\\\

O.engines=[

 ForceResetter(),

 InsertionSortCollider(

  [Bo1_Sphere_Aabb(),

   Bo1_Facet_Aabb(),

   Bo1_Box_Aabb(),

   Bo1_Wall_Aabb()]),

 InteractionLoop(

  [Ig2_Sphere_Sphere_ScGeom(),

   Ig2_Facet_Sphere_ScGeom(),

   Ig2_Wall_Sphere_ScGeom(),

   Ig2_Box_Sphere_ScGeom()],

  [Ip2_FrictMat_FrictMat_FrictPhys()],

  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),

 PyRunner(command='alturasMaximasCapas()',

   iterPeriod=1,label='checker'),

 NewtonIntegrator(damping=0.5,gravity=[0,-9.81,0],label='newton')

]

factor=1
O.dt=factor*utils.PWaveTimeStep() #time step

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Sergio
Solved:
Last query:
Last reply:
Revision history for this message
Christian Jakob (jakob-ifgt) said : 		#1

hi sergio,

it is a bit tricky, but possible. one idea:

- give your clump clouds different color
- loop over clumps
- find maximum

example:

for b in O.bodies:
   if b.color= ???:
     maxHeight1 = max(maxHeight1,b.state.pos[2])
   if b.color= ???:
     maxHeight2 = max(maxHeight2,b.state.pos[2])
...

regards,

christian

p.s. i know, that i should use launchpad, but i can not access it
right now. blocked in china?!

Revision history for this message
Sergio (serrodam4) said : 		#2

Hi Christian,

Thank you for your comment. The loops taking in care the colors seems a very good idea!

I have tried to implement it, but the results are not as expected. I have write the following function (called by a Pyrunner):

def alturasMaximasCapas():
 for b in O.bodies:
  if b.isClump==True and b.shape.color==(0.3,0.9,1):
   maxHeight1=max(maxHeight1,b.state.pos[1])
   print maxHeight1
  elif b.isClump==True and b.shape.color==(0.9,0.2,1):
   maxHeight1=max(maxHeight1,b.state.pos[1])
   print maxHeight2
  elif b.isClump==True and b.shape.color==(0.3,0.2,0.1):
   maxHeight1=max(maxHeight1,b.state.pos[1])
   print maxHeight3

But I have no output about any maxHeightX. If I try to call maxHeight1 in terminal, it returns 'name 'maxHeight1' is not defined'
If I give a initial value out of the def function, i.e. maxHeight1=-1000 to ensure it has to change, the terminal returns the initial value without any change.

I don't know where the problem could be. It is the first time I use a Python based program, so I have no knowledge of this programming language. :(

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#3

"global" is the keyword.
Bruno

a=0
def f():
    a=1
    print a

f()
>> 1
 print a
>> 0

a=0
def f():
    global a
    a=1
    print a

f()
>> 1
 print a
>> 1

Revision history for this message
Sergio (serrodam4) said : 		#4

Thank you for the comment Bruno!

So, I've tried to modify the code as follows:

maxHeight1=-1000
maxHeight2=-1000
maxHeight3=-1000
#///////////// Layer heights \\\\\\\\\\\\\\\
def altCapas():
 global maxHeight1
 global maxHeight2
 global maxHeight3
 for b in O.bodies:
  if b.shape.color==(0.3,0.9,1):
   maxHeight1=max(maxHeight1,b.state.pos[1])
   print 'Maximum height layer 1:',maxHeight1
  elif b.shape.color==(0.9,0.2,1):
   maxHeight2=max(maxHeight2,b.state.pos[1])
   print 'Maximum height layer 2:',maxHeight2
  elif b.shape.color==(0.3,0.2,0.1):
   maxHeight3=max(maxHeight3,b.state.pos[1])
   print 'Maximum height layer 3:',maxHeight3

Now, the function works. BUT there is a problem. Since I want the highest clump position of each layer, this function only works properly if the clumps move up. If they move down, as in my case, the highest position in each layer results in first iter. (As iter 1 pos > iter 2 pos). It would be needed something else to achieve the highest clump position in each time step. I'll keep trying.

Thanks again for the help!!!

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#5

Yes, because the max should be reste each time. Actually, it seems to me that you simply want to return values.
Why not return them (then you can skip "global")?

def altCapas():
   ...
   return maxHeight1,maxHeight2,maxHeight3

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#6

reste <= reseted

Revision history for this message
Sergio (serrodam4) said : 		#7

Thank you very much again, Bruno!

Finally, I could make a function that solves my problem:

def altCapas():
 global maxHeight1
 global maxHeight2
 global maxHeight3
 for b in O.bodies:
  if b.shape.color==(0.3,0.9,1):
   maxHeight1=max(maxHeight1,b.state.pos[1])
  elif b.shape.color==(0.9,0.2,1):
   maxHeight2=max(maxHeight2,b.state.pos[1])
# print 'Maximum height layer 2:',maxHeight2
# print 'Position y: ', b.state.pos[1]
  elif b.shape.color==(0.3,0.2,0.1):
   maxHeight3=max(maxHeight3,b.state.pos[1])
 print 'Maximum height layer 1:',maxHeight1
 print 'Maximum height layer 2:',maxHeight2
 print 'Maximum height layer 3:',maxHeight3
 print 'The number of particles that exceed the limit is: ', numOfPart

 maxHeight1=-500
 maxHeight2=-500
 maxHeight3=-500

Now it works. I tried to use the 'return', but I did'n achieved a good results. Thank you again Christian and Bruno for your help ;)