# calculate the mass of the particles

Asked by Estefany Carmona Alvarez on 2020-02-18

I have generated a code for the simulation of the sand free fall.

This simulation consists of dropping the spheres that simulate the sand, so that they are deposited inside two cylindrical vessels. As a result of this simulation I need to obtain the weight or mass of the particles that remain within each cylinder.

In the Yade documentation that I have consulted, I saw that the state.mass function fulfills this need, but since I do not have enough domain in this programming language, I needed to know if anyone has any example of how to use this code or if they have more information ...

The code I am using is as follows:
postscript: the code comments are in Spanish (Latin)

#Altura de caida del material
hc=0

#EQUIPO

x=0
y=0
z=0

f=0.45
a=0.50
h=1

#estrutura del equipo de pluviacion

#Mallas
#O.bodies.append(utils.wall((0,0,(.2+hc)),axis=2,color=(1,0.945,0.078),sense=-1))

#(0,0,(.2+hc)),axis=2,color(1,0.945,0.078),sense=0))

#O.bodies.append(utils.wall((0,0,.2),axis=2,sense=0))

#PROBETA

#Altura de la probeta
height=0.1

#Probeta 1
xp1=0.1
yp1=0
zp1=1

#Probeta 2
xp2=0.1
yp2=0
zp2=1

#esquina inferior
xi=0.225
yi=0.225
zi=0.2

#esquina superior
xs=0.225
ys=0.25
zs=0.5

p=0.903

#tamices
#los diametrso y los porcentajes acumulados siempre se deben ubicar de menor a mayor
#por cada t debe haber un p

#Diametro de la abertura de tamices psdSizes
t1=0.008
t2=0.015
t3=0.030
t4=0.085
t5=0.200
t6=0.475

#Porcentajes retenidos del material psdCumm
p1=0
p2=0.002
p3=0.111
p4=0.889
p5=0.9994
p6=1

#Material
#angulo de friccion
#Modulo de Young
young=(1e8) #Pa
#Poisson
poisson=0.5
#contenidod de vacios
#label="SndG"
))

#Algoritmo YADE para el cuerpo de las esferas

sp=pack.SpherePack()

sp.makeCloud((xi,-yi,(zi+hc)), (-xs,ys,(zs+hc)), porosity=porosidad, psdSizes=[t1,t2,t3,t4,t5,t6], psdCumm=[p1,p2,p3,p4,p5,p6], distributeMass=False)

sp.toSimulation(material=snd,color=(0.882,0.866,0.839))

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),
Bo1_Facet_Aabb(),
Bo1_Cylinder_Aabb(),
Bo1_Wall_Aabb()
]),
InteractionLoop(
# handle sphere+sphere and facet+sphere collisions
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
#GravityEngine(gravity=(0,0,-9.81)),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
PyRunner()
]
O.dt=PWaveTimeStep()

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
2020-02-24
2020-02-24
 Jan Stránský (honzik) said on 2020-02-18: #1

Hello,

I have tried your code, but the spheres are larger than the cylinders and therefore no spheres are deposited inside any cylinder and their mass is zero.
Please update the code such that it results in something meaningful.

Otherwise the task is easy, something like:
###
cyl1 = pack.inCylinder(...) # , predicate for a cylinder
mass = 0.0
for b in O.bodies:
if cyl1(b.state.pos): # b is inside cyl1
mass += b.state.mass
###

or a one-liner :-)
mass = sum(b.state.mass for b in O.bodies if cyl1(b.state.pos))

cheers
Jan

 Estefany Carmona Alvarez (jecarmona91) said on 2020-02-19: #2

1. In this code is the actual size of the particles, which if they are deposited in the cylinders ...

#Altura de caida del material
hc=0

#EQUIPO

x=0
y=0
z=0

f=0.45
a=0.50
h=1

#estrutura del equipo de pluviacion

#Mallas
#O.bodies.append(utils.wall((0,0,(.2+hc)),axis=2,color=(1,0.945,0.078),sense=-1))

#(0,0,(.2+hc)),axis=2,color(1,0.945,0.078),sense=0))

#O.bodies.append(utils.wall((0,0,.2),axis=2,sense=0))

#PROBETA

#Altura de la probeta
height=0.1

#Probeta 1
xp1=0.1
yp1=0
zp1=1

#Probeta 2
xp2=0.1
yp2=0
zp2=1

#esquina inferior
xi=0.225
yi=0.225
zi=0.2

#esquina superior
xs=0.225
ys=0.25
zs=0.5

p=0.903

#tamices
#los diametrso y los porcentajes acumulados siempre se deben ubicar de menor a mayor
#por cada t debe haber un p

#Diametro de la abertura de tamices psdSizes
t1=0.008
t2=0.015
t3=0.030
t4=0.085
t5=0.200
t6=0.475

#Porcentajes retenidos del material psdCumm
p1=0
p2=0.002
p3=0.111
p4=0.889
p5=0.9994
p6=1

#Material
#angulo de friccion
#Modulo de Young
young=(1e8) #Pa
#Poisson
poisson=0.5
#contenidod de vacios
#label="SndG"
))

#Algoritmo YADE para el cuerpo de las esferas

sp=pack.SpherePack()

sp.makeCloud((xi,-yi,(zi+hc)), (-xs,ys,(zs+hc)), porosity=porosidad, psdSizes=[t1,t2,t3,t4,t5,t6], psdCumm=[p1,p2,p3,p4,p5,p6], distributeMass=False)

sp.toSimulation(material=snd,color=(0.882,0.866,0.839))

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),
Bo1_Facet_Aabb(),
Bo1_Cylinder_Aabb(),
Bo1_Wall_Aabb()
]),
InteractionLoop(
# handle sphere+sphere and facet+sphere collisions
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
#GravityEngine(gravity=(0,0,-9.81)),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
PyRunner()
]
O.dt=PWaveTimeStep()

###
cyl1 = pack.inCylinder(...) # , predicate for a cylinder
mass = 0.0
for b in O.bodies:
if cyl1(b.state.pos): # b is inside cyl1
mass += b.state.mass
###

or a one-liner :-)
mass = sum(b.state.mass for b in O.bodies if cyl1(b.state.pos))

I am new to Yade and in the available documents there is not much information about it.

Thanks for your help Jan !!

 Jan Stránský (honzik) said on 2020-02-19: #3

> I am new to Yade and in the available documents there is not much information about it.

Yade uses Python as an user interface. In python, you have access to almost any data in the simulation (tutorial  and user's manual  is full of this).
So you have access to positions of all particles, then it is easy to check if they are in a cylinder or not and sum mass of those which are. There are several options e.g.
###
d = body.state.pos - center # vector between cylinder center and particle center
d = 0 # only xy components are relevant
distance = d.norm() # horizontal distance
return distance < radius # particle is inside, if distance < raidus
def masses():
mass1 = mass2 = 0.0 # intially zero
center1 = Vector3(xp1,yp1,0) # center of 1st cylinder
center2 = Vector3(-xp2,yp2,0) # center of 2nd cylinder
for b in O.bodies: # test all particles
if not isinstance(b.shape,Sphere): continue # skip non-spherical
if isBodyInCylinder(b,center1,radius): # if a particle is inside 1st cylinder ...
mass1 += b.state.mass # ... increase mass1 by the particle mass
if isBodyInCylinder(b,center2,radius): # if a particle is inside 2nd cylinder ...
mass2 += b.state.mass # ... increase mass2 by the particle mass
return mass1,mass2 # return computed values
###

cheers
Jan

 Estefany Carmona Alvarez (jecarmona91) said on 2020-02-20: #4

> I included the code you indicated, but now I don't know where I can check the mass of the particles or if it is necessary to include some other function to be able to visualize this value.

#Altura de caida del material
hc=0

#EQUIPO

x=0
y=0
z=0

f=0.45
a=0.50
h=1

#estrutura del equipo de pluviacion

#Mallas
#O.bodies.append(utils.wall((0,0,(.2+hc)),axis=2,color=(1,0.945,0.078),sense=-1))

#(0,0,(.2+hc)),axis=2,color(1,0.945,0.078),sense=0))

#O.bodies.append(utils.wall((0,0,.2),axis=2,sense=0))

#PROBETA

#Altura de la probeta
height=0.1

#Probeta 1
xp1=0.1
yp1=0
zp1=1

#Probeta 2
xp2=0.1
yp2=0
zp2=1

#esquina inferior
xi=0.225
yi=0.225
zi=0.2

#esquina superior
xs=0.225
ys=0.25
zs=0.5

p=0.903

#tamices
#los diametrso y los porcentajes acumulados siempre se deben ubicar de menor a mayor
#por cada t debe haber un p

#Diametro de la abertura de tamices psdSizes
t1=0.00008
t2=0.00015
t3=0.00030
t4=0.00085
t5=0.00200
t6=0.00475

#Porcentajes retenidos del material psdCumm
p1=0
p2=0.002
p3=0.111
p4=0.889
p5=0.9994
p6=1

#Material
#angulo de friccion
#Modulo de Young
young=(1e8) #Pa
#Poisson
poisson=0.5
#contenidod de vacios
#label="SndG"
))

#Algoritmo YADE para el cuerpo de las esferas

sp=pack.SpherePack()

sp.makeCloud((xi,-yi,(zi+hc)), (-xs,ys,(zs+hc)), porosity=porosidad, psdSizes=[t1,t2,t3,t4,t5,t6], psdCumm=[p1,p2,p3,p4,p5,p6], distributeMass=False)

sp.toSimulation(material=snd,color=(0.882,0.866,0.839))

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),
Bo1_Facet_Aabb(),
Bo1_Cylinder_Aabb(),
Bo1_Wall_Aabb()
]),
InteractionLoop(
# handle sphere+sphere and facet+sphere collisions
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
#GravityEngine(gravity=(0,0,-9.81)),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
PyRunner()
]
O.dt=PWaveTimeStep()

d=body.state.pos-center
d=0
distance=d.norm()

def masses():
mass1=mass2=0.0
center1=(xp1,yp1,0)
center2=(-xp2,yp2,0)
for b in O.bodies:
if not isinstance(b.shape,Sphere): continue
mass1+=b.state.mass
mass2+=b.state.mass
return mass1,mass2

> In the previous answer, you mentioned some supporting documents. You could attach them to me to consult them. They would be helpful to understand some functions that I need and for which there is not much information on the web.

thanks!!!

 Jan Stránský (honzik) said on 2020-02-20: #5

> but now I don't know where I can check the mass of the particles or if it is necessary to include some other function to be able to visualize this value.

run your simulation until a point when it is considered finished. Stop it. Call
masses()
in yade interactive terminal. The returned value are the masses in the two cylinders. You can print it or do whatever you want:
##
mass1,mass2 = masses()
print(mass1,mass2)
totalMass = mass1 + mass2
...
##

> In the previous answer, you mentioned some supporting documents. You could attach them to me to consult them. They would be helpful to understand some functions that I need

sure, just forgot to add it last time, see below

> for which there is not much information on the web.

I don't understand this..
"for b in O.bodies:" is surely documented as well as that b.state.mass is the mass of the body and b.state.pos its center.
Definitely you can easily find information on the web how to do a sum or how to check if a point is inside a cylinder..
BUT, Yade is open source and you are free (and welcome actually) to add such example or improve the documentation (or its part) if you find it not good! :-)

cheers
Jan

 Estefany Carmona Alvarez (jecarmona91) said on 2020-02-20: #6

Jan !! Help...

It seemed an easy function to use, but despite defining the function that stops the simulation in relation to the reduction of force, the simulation continues without showing me the results.

I included the following code in the simulation

####
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),
Bo1_Facet_Aabb(),
Bo1_Cylinder_Aabb(),
Bo1_Wall_Aabb()
]),
InteractionLoop(
# handle sphere+sphere and facet+sphere collisions
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
#GravityEngine(gravity=(0,0,-9.81)),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
PyRunner(command="checkUnbalanced()"),
PyRunner(command='isBodyInCylinder()'),
PyRunner(command='masses()')
]
O.dt=PWaveTimeStep()

def checkUnbalanced():
if unbalancedForce<0.05: O.pause()

d=body.state.pos-center
d=0
distance=d.norm()

def masses():
mass1=mass2=0.0
center1=(xp1,yp1,0)
center2=(-xp2,yp2,0)
for b in O.bodies:
print(mass1,mass2)
if not isinstance(b.shape,Sphere): continue
mass1+=b.state.mass
mass2+=b.state.mass
return mass1,mass2

###
According to what I understand, this function causes the simulation to stop when the equilibrium force of the particles is less than 0.05, but although the particles visually no longer interact, the simulation continues and in the script there is no result of mass.

Thanks for the documents, they were very necessary ... I hope that with the development of this project I can contribute something very interesting to the Yade community ...

It seems easy to use this software, but it has become an excellent professional challenge ...

thanks Jan!

 Jan Stránský (honzik) said on 2020-02-21: #7

> PyRunner(command="checkUnbalanced()"),
> PyRunner(command='isBodyInCylinder()'),
> PyRunner(command='masses()')

use exatly one of periods (iterPeriod , virtPeriod  or realPeriod ) for PyRunner, otherwise it is not run.

cheers
Jan

 Estefany Carmona Alvarez (jecarmona91) said on 2020-02-21: #8

Jan!
With the commands you indicated, I can now see that the simulation stops for the indicated condition.

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),
Bo1_Facet_Aabb(),
Bo1_Cylinder_Aabb(),
Bo1_Wall_Aabb()
]),
InteractionLoop(
# handle sphere+sphere and facet+sphere collisions
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
#GravityEngine(gravity=(0,0,-9.81)),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
PyRunner(command='checkUnbalanced()',realPeriod=1),
PyRunner(command='isBodyInCylinder()',realPeriod=1),
PyRunner(command='masses()',realPeriod=1)
]
O.dt=PWaveTimeStep()

def checkUnbalanced():
if unbalancedForce()<0.05: O.pause()

d=body.state.pos-center
d=0
distance=d.norm()

def masses():
mass1=mass2=0.0
center1=(xp1,yp1,zp1)
center2=(-xp2,yp2,zp2)
for b in O.bodies:

if not isinstance(b.shape,Sphere): continue
mass1+=b.state.mass
mass2+=b.state.mass
return mass1,mass2
print(mass1,mass2)
>>>
1. But when executing the isBodyInCylinder function the script shows the following alert:

###
TypeError: isBodyInCylinder () takes exactly 3 arguments (0 give)
###

I understand that the arguments requested are the body (b), the center (center1 or center2) and the radius. I understand that these arguments are defined, so I don't know what the error might be.

>>>
2. If it is ok to add a PyRunner for each command created
###
PyRunner (command = 'checkUnbalanced ()', realPeriod = 1),
PyRunner (command = 'isBodyInCylinder ()', realPeriod = 1),
PyRunner (command = 'masses ()', realPeriod = 1)
###

>>>>>
3. Is there any reference value for the real Period or how can I ensure that the period indicated for the function is sufficient.

Thanks Jan!!

 Bruno Chareyre (bruno-chareyre) said on 2020-02-21: #9

> TypeError: isBodyInCylinder () takes exactly 3 arguments (0 give)

Dear Estefany,
Please keep in mind point 5 in .
Although Jan may keep replying the questions because he is super-nice, it's now clearly not about computing total mass.
It may be beneficial for you to step back maybe, and study python functionalities and error messages a bit, before diving again into Yade.
It's very confusing for someone to not distinguish what's python and what's Yade, and I smell it could be your case. :)
Regards
Bruno

Side note: consider indenting with tabs or four spaces (at least), most text editors can do that for you. It's much more readable.

 Jan Stránský (honzik) said on 2020-02-24: #10

Hello,

1.
as the error says, isBodyInCylinder () takes exactly 3 arguments, but you gave 0 arguments..
3 arguments = 3 values inside "calling" brackets.

solution:
isBodyInCylinder is just an auxiliary function for masses() function, just do not call it from PyRunner (delete that PyRunner from the O.engines)

2.
I guess yes for this case. If to use iterPeriod, virtPeriod or realPeriod depends always on specific context.

3.
see above, it is always different. If a specific value is sufficient, it is easy to test.

cheers
Jan

 Estefany Carmona Alvarez (jecarmona91) said on 2020-02-24: #11

Jan !!! you cannot imagine how grateful I am with your help, it is very valuable and important to me.
If there is any way I can contribute to your help, please let me know.

With your contribution I corrected the code and this was the result:

#Altura de caida del material
hc=0

#EQUIPO

x=0
y=0
z=0

f=0.45
a=0.50
h=1

#estrutura del equipo de pluviacion

#Mallas
#O.bodies.append(utils.wall((0,0,(.2+hc)),axis=2,color=(1,0.945,0.078),sense=-1))

#(0,0,(.2+hc)),axis=2,color(1,0.945,0.078),sense=0))

#O.bodies.append(utils.wall((0,0,.2),axis=2,sense=0))

#PROBETA

#Altura de la probeta
height=0.1

#Probeta 1
xp1=0.1
yp1=0
zp1=1

#Probeta 2
xp2=0.1
yp2=0
zp2=1

#esquina inferior
xi=0.225
yi=0.225
zi=0.2

#esquina superior
xs=0.225
ys=0.25
zs=0.5

p=0.903

#tamices
#los diametrso y los porcentajes acumulados siempre se deben ubicar de menor a mayor
#por cada t debe haber un p

#Diametro de la abertura de tamices psdSizes
t1=0.00008
t2=0.00015
t3=0.00030
t4=0.00085
t5=0.00200
t6=0.00475

#Porcentajes retenidos del material psdCumm
p1=0
p2=0.002
p3=0.111
p4=0.889
p5=0.9994
p6=1

#Material
#angulo de friccion
#Modulo de Young
young=(1e8) #Pa
#Poisson
poisson=0.5
#contenidod de vacios
#label="SndG"
))

#Algoritmo YADE para el cuerpo de las esferas

sp=pack.SpherePack()

sp.makeCloud((xi,-yi,(zi+hc)), (-xs,ys,(zs+hc)), porosity=porosidad, psdSizes=[t1,t2,t3,t4,t5,t6], psdCumm=[p1,p2,p3,p4,p5,p6], distributeMass=False)

sp.toSimulation(material=snd,color=(0.882,0.866,0.839))

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),
Bo1_Facet_Aabb(),
Bo1_Cylinder_Aabb(),
Bo1_Wall_Aabb()
]),
InteractionLoop(
# handle sphere+sphere and facet+sphere collisions
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
#GravityEngine(gravity=(0,0,-9.81)),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
PyRunner(command='checkUnbalanced()',realPeriod=1),
PyRunner(command='masses()',realPeriod=1)
]
O.dt=PWaveTimeStep()

def checkUnbalanced():
if unbalancedForce()<0.05: O.pause()

d=body.state.pos-center
d=0
distance=d.norm()

def masses():
mass1=mass2=0
center1=(xp1,yp1,zp1)
center2=(-xp2,yp2,zp2)
for b in O.bodies:

if not isinstance(b.shape,Sphere): continue
mass1+=b.state.mass

mass2+=b.state.mass
return mass1, mass2
mass1,mass2=masses()
print"masa de cada cilindro (g):",(mass1, mass2)

totalMass=mass1+mass2
print"masa total de los cilindros (g):",(totalMass)

masa=sum([b.state.mass for b in O.bodies])
print"masa del material de simulacion (g):",(masa)

Now I only have a small question, the value that the script shows me before visualizing the simulation, corresponds to the value of the mass of the particles that remain inside each cylinder? is that when compared with the tax experiment they are much inferior.

 Jan Stránský (honzik) said on 2020-02-24: #12

> the value ... corresponds to the value of the mass of the particles that remain inside each cylinder?

yes

> is that when compared with the tax experiment they are much inferior.

yes, it is a simulation, not an experiment :-)

cheers
Jan

 Chareyre (bruno-chareyre-9) said on 2020-02-24: #13

Is there O.run() somewhere, or you are speaking of mass at iteration 0?
B

Le lun. 24 févr. 2020 19:48, Jan Stránský <
<email address hidden>> a écrit :

> Question #688837 on Yade changed:
>
> Jan Stránský posted a new comment:
> > the value ... corresponds to the value of the mass of the particles
> that remain inside each cylinder?
>
> yes
>
> > is that when compared with the tax experiment they are much inferior.
>
> yes, it is a simulation, not an experiment :-)
>
> cheers
> Jan
>
> --
>
> _______________________________________________