Extract normal force

Hello,

> I am interested on the extraction of the normal contact forces

What is "the extraction"? printing, saving, ...?
What is "the normal contact forces"? individual force vectors? magnitudes? sum of the vectors? sum of magnitudes? ... ?

> I use i.phys.normalForce.dot(i.geom.normal)

what is "i"?

> the length of the simulation is 4098

"length" is number of time steps?

> by extracting the normal forces

what is "extracting"?
what is "the normal forces"? vectors? magnitudes? ... ?

> I expected to have a force vector of the same length as the number of simulation

what is "the number of simulation"?
what is "a force vector"? I would expect "a force vector" to be a "standard Euclidean vector"..

> I would like to know why I don't get the length of the vector equal to the length of the simulation?
> How can i do to have the length of the vector normal force, with the same length of the total number of iteration?

What is "the vector"?
What is "the length of the vector"?
What is "the length of the simulation" and why would it equal to "the length of the vector"?
What is "the vector normal force"?

From the question, we can hardly even guess, please read [1] and provide more information.

> Can someone help me?

Yes, but first please help us. E.g. by providing a MWE [1], which would (should) have answered most of the "what is ..." questions..

cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

Hi,

Thanks for your response, and you can find here after some response at your questions;

What is "the extraction"? printing, saving, ...?
The term extraction used here means , Saving force
What is "the normal contact forces"? individual force vectors? magnitudes? sum of the vectors? sum of magnitudes? ... ?
I call normal contact force, the magnitudes of the force of contact between ball and mill in the normal direction, saving at each step time; i need to save this force of contact at each time step. To know the evolution of this in time.
what is "the number of simulation"?
the number of simulation is the number of iteration here is 4098
what is "a force vector"? I would expect "a force vector" to be a "standard Euclidean vector"..
I call force vector, or normal force, in my question here, the vector contenning the value of magnitudes force in normal direction saving at each step time.

 the force at each step time. if the length of simulation, or the vector time length is 4098, i expect to have the saving vector force with the same length.

Briefly i want to save the evolution of interaction force in time; reason why i insert a loop for j in O.time in my script; But i can't have the result i need.
I don't know if there is an other way, to save the contact force at each step time, and to plot his evolution in time;

I' m trying to do it in my mwe here after:

# Material
rho = 7640.0 # [kg/m^3] density (of the balls)
kn = 1.0e5 # [N/m] equivalent normal stiffness
kt = kn # [N/m] equivalent tangential stiffness
en = 0.3 # [-] normal coefficient of restitution (sphere/sphere)
mu = 0.75 # [-] friction coefficient
cn = 0.75
ct = cn

# Geometry/Filling
length = 0.4 # [m] axial length of the mill
rMean = 0.772/2.0 # [m] liner profile mean radius
rBall = 0.01 # [m] ball radius
fRatio = 1.0 # [%] filling ratio
sphRadFuzz = 0.8

# Simulation

tEnd = 1.0 # [s] total simulation time
O.dt = 2.44140e-4 # [s] fixed time step
tRec = 1.0e-2 # [s] recording interval

#material of cylinder
cylMat = O.materials.append(ViscElMat(kn=2*kn,ks=2*kt,cn=2*cn,
                                             cs=2*ct,frictionAngle=.1))
#material of spheres
ballMat = O.materials.append(ViscElMat(kn=2*kn,ks=2*kt,cn=2*cn,
                                             cs=2*ct,frictionAngle=.1,density=rho))

ids=O.bodies.append(geom.facetCylinder((0,0,0),0.8,0.4,segmentsNumber=64, wallMask=7, angleRange=None,
                   radiusTopInner=-1, radiusBottomInner=-1,material=cylMat))

# define domains for initial cloud of red and blue spheres
packHt=.8*rMean # size of the area
bboxes=[(Vector3(-.5*length,-.5*packHt,-.5*packHt),Vector3(.5*length,0,.5*packHt)),(Vector3(-.5*length,-.5,-.5*packHt),Vector3(.5*length,.5*packHt,.5*packHt))]
colors=(1,0,0),(0,0,1)
for i in (0,1): # red and blue spheres
 sp=pack.SpherePack();
 bb=bboxes[i];
 vol=(bb[1][0]-bb[0][0])*(bb[1][1]-bb[0][1])*(bb[1][2]-bb[0][2])
 sp.makeCloud(bb[0],bb[1],rBall,sphRadFuzz,int(.25*vol/((4./3)*pi*rBall**3)),False)
 O.bodies.append([utils.sphere(s[0],s[1],material=ballMat,color=colors[i]) for s in sp])

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_MindlinPhys()],
        [Law2_ScGeom_MindlinPhys_Mindlin()]
    ),

   #GravityEngine(gravity=(0,-3e4,0)), # gravity artificially high, to make it faster going ;-)
   RotationEngine(ids=ids,angularVelocity=1000,rotateAroundZero=True,rotationAxis=(0,0,1)),
   NewtonIntegrator(damping=0,gravity=[0,0,-9.81]),
   PyRunner(command='exportForces()',iterPeriod=1)
]
def exportForces():

  # get data

for j in O.time:
  for i in O.interactions:
    if isinstance(O.bodies[i.id1].shape,Sphere) \
        and isinstance(O.bodies[i.id2].shape,Facet):

  totalNormalForce = i.phys.normalForce

  # Save data

  out=open(os.getcwd()+'/output/NormalForceMagnitude','a')
  out.write('%s\n'%str(NormalForceMagnitude))
  out.close()
###
O.run(int(math.ceil(tEnd/O.dt))+1)

Hello,

> I call normal contact force, the magnitudes of the force of contact between ball and mill in the normal direction
> > the force at each step time
> i want to save the evolution of interaction force in time

still this is very confusing.
In the simulation, you have (possibly) many balls and many ball-mill interactions and many forces.
What is "(the) force"?
- force of one specific interaction?
- sum of force vectors?
- sum of force magnitudes?
- ... ?

> I call force vector, or normal force, in my question here, the vector contenning the value of magnitudes force in normal direction saving at each step time.

I strongly recommend to use a different name for the "vector" in this case not to be confused with "ordinary" force vector.
Maybe array? field?

> i insert a loop for j in O.time in my script; But i can't have the result i need.

it results in an error, no surprise you do not have the result.
Fix errors first.

> I don't know if there is an other way, to save the contact force at each step time

There are many ways.
To choose a suitable one, first of all you need to know (or explain it properly in the case you want help from others) what you want. Which does not seem to be the case (see above).

> ... mwe ...

W = working, which is not the case here. Please, fix the indentation to be valid Python code.
Moreover, O.time is a number and you cannot iterate it "for j in O.time"

> totalNormalForce = i.phys.normalForce

**guessing** it should be inside the sphere-facet condition:
- why the equal sign? Shouldn't the value of totalNormalForce be incremented by i.phys.normalForce instead of rewriting?
- you write you want normal force magnitudes. here you assign/sum normal force vector. What exactly do you want?
- you can save the data "manually" as you do here. But you should clear the output file first (the code appends data to whatever already is in the file)
- instead of "manual" saving, Yade "standardly" uses plot module for this purposes [1]

cheers
Jan

[1] https://yade-dem.org/doc/user.html#tracking-variables

Hi, Thank you, This some response,

still this is very confusing.
In the simulation, you have (possibly) many balls and many ball-mill interactions and many forces.
What is "(the) force"?
- force of one specific interaction?
- sum of force vectors?
- sum of force magnitudes?

Exactly you are right. In the simulation, several interactions take place between balls and mill, because of the presence of several balls. That's why I want to save at each iteration (or step time), the sum of forces magnitudes.

I'will provide my mwe soon

This is my MWE;

As I said before, I want to save for each iteration, the sum of force (amplitude) of interactions between balls and mill;
In this MWE, I have rather calculated the forces when there is interaction; But I wish to calculate it rather for each iteration; this will allow me to represent the evolution of its forces in time;

from yade import pack,utils, plot,ymport,qt
#######################################################################################

fr = 0.1;rho=7640.0;length = 4.0;rMean = 0.772/2.0;sphRadFuzz = 0.8;rBall = 0.04 ;

MillMat = O.materials.append(ViscElMat(kn=2*1.0e5,ks=2*1.0e5,frictionAngle=2*fr))
ballMat = O.materials.append(ViscElMat(kn=2*1.0e5,ks=2*1.0e5,frictionAngle=fr,density=rho))

#######################################################################################

# Create Mill from facets
Mill=O.bodies.append(geom.facetCylinder((0.0,0.0,0.5),8.0,height=4.0,segmentsNumber=64,wallMask=7,material=MillMat))
#Walls=O.bodies.append([utils.wall(position=(0,0,0),axis=2,sense=0,material=MillMat),
# utils.wall(position=(0,length,0),axis=2,sense=0,material=MillMat)
# ])
# define domains for initial cloud of red and blue spheres
packHt=.8*rMean # size of the area
bboxes=[(Vector3(-.5*length,-.5*packHt,-.5*packHt),Vector3(.5*length,0,.5*packHt)),(Vector3(-.5*length,-.5,-.5*packHt),Vector3(.5*length,.5*packHt,.5*packHt))]
colors=(1,0,0),(0,0,1)
for i in (0,1):
 sp=pack.SpherePack();
 bb=bboxes[i];
 vol=(bb[1][0]-bb[0][0])*(bb[1][1]-bb[0][1])*(bb[1][2]-bb[0][2])
 sp.makeCloud(bb[0],bb[1],rBall,sphRadFuzz,int(.25*vol/((4./3)*pi*rBall**3)),False)
 O.bodies.append([utils.sphere(s[0],s[1],material=ballMat,color=colors[i]) for s in sp])

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,-9.81,0),damping=0.0),

    RotationEngine(ids=Mill,rotationAxis=[0,0,-1],rotateAroundZero=True,
                 angularVelocity=1),

   ]

O.dt=.5*PWaveTimeStep()

def saveForces():

   for i in O.interactions:

    if isinstance(O.bodies[i.id1].shape,Facet) or isinstance(O.bodies[i.id2].shape,Facet):
       MagnitudeNormalForce = sum((i.phys.normalForce.norm()))
    else :
       MagnitudeNormalForce = 0

   # out=open(os.getcwd()+'/Resultat/MagnitudenormalForce','a')
   # out.write('%s\n'%str(MagnitudenormalForce))
   # out.close()

O.run(100000)
#O.save('Mymodel2.yade')

You defined saveForces, but it is not used anywhere.
You have to use saveForces somehow, most likely form a PyRunner (as you already had in the previous codes).
cheers
Jan

Hi Jan,

This is my updated MWE, as i tell you before, i need the simulation show the plot on the screen, and update while the simulation runs in time, and save the result in txt format.

rom yade import pack,utils, plot,ymport,qt
#######################################################################################

fr = 0.1;rho=7640.0;length = 4.0;rMean = 0.772/2.0;sphRadFuzz = 0.8;rBall = 0.04 ;

MillMat = O.materials.append(ViscElMat(kn=2*1.0e5,ks=2*1.0e5,frictionAngle=2*fr))
ballMat = O.materials.append(ViscElMat(kn=2*1.0e5,ks=2*1.0e5,frictionAngle=fr,density=rho))

#######################################################################################

# Create Mill from facets
Mill=O.bodies.append(geom.facetCylinder((0.0,0.0,0.5),8.0,height=4.0,segmentsNumber=64,wallMask=7,material=MillMat))
#Walls=O.bodies.append([utils.wall(position=(0,0,0),axis=2,sense=0,material=MillMat),
# utils.wall(position=(0,length,0),axis=2,sense=0,material=MillMat)
# ])
# define domains for initial cloud of red and blue spheres
packHt=.8*rMean # size of the area
bboxes=[(Vector3(-.5*length,-.5*packHt,-.5*packHt),Vector3(.5*length,0,.5*packHt)),(Vector3(-.5*length,-.5,-.5*packHt),Vector3(.5*length,.5*packHt,.5*packHt))]
colors=(1,0,0),(0,0,1)
for i in (0,1):
 sp=pack.SpherePack();
 bb=bboxes[i];
 vol=(bb[1][0]-bb[0][0])*(bb[1][1]-bb[0][1])*(bb[1][2]-bb[0][2])
 sp.makeCloud(bb[0],bb[1],rBall,sphRadFuzz,int(.25*vol/((4./3)*pi*rBall**3)),False)
 O.bodies.append([utils.sphere(s[0],s[1],material=ballMat,color=colors[i]) for s in sp])

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,-9.81,0),damping=0.0),

    RotationEngine(ids=Mill,rotationAxis=[0,0,-1],rotateAroundZero=True,
                 angularVelocity=1),
    PyRunner(command='addPlotData()',iterPeriod=100),
   ]

O.dt=.5*PWaveTimeStep()

# collect history of data which will be plotted
def addPlotData():

  for i in O.interactions:

    if isinstance(O.bodies[i.id1].shape,Facet) or isinstance(O.bodies[i.id2].shape,Facet):
       Force = sum((O.forces.f(i).norm()))
    else :
       Force = 0

  plot.addData(i=O.iter,Force=sum(O.forces.f(i).norm()))

plot.plots={'i':('Force')}

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

O.run(100000)
#O.save('Mymodel2.yade')

but when i run the simulation i get this error;

ArgumentError: Python argument types in
    ForceContainer.f(ForceContainer, Interaction)
did not match C++ signature:
    f(yade::pyForceContainer {lvalue}, long id, bool sync=False)
In [1]: Erreur de segmentation (core dumped)

Can i have some correction in my code please?

code 1:
> def exportForces():
> totalNormalForce = i.phys.normalForce

code 2:
> def saveForces():
> MagnitudeNormalForce = sum((i.phys.normalForce.norm()))

code 3:
> def addPlotData():
> Force = sum((O.forces.f(i).norm()))

Every update, you completely changed the namings. Yes, it is just naming and here it is clear, but next time please try to stay consistent.
There is no reason for such renaming and then I would prefer not to rename the stuff to prevent confusion.

Every update, you changed the force computation method, without any comment. Next time please consider giving some comment about the reasons.
Ok, from the discussion it came out that the magnitude should be force.norm() instead of force vector itself, so i.phys.normalForce -> i.phys.normalForce.norm() is reasonable.
But why the sum (inside for loop)?
And why i.phys.normalForce -> O.forces.f ?(!?!?)

> This is my updated MWE
> rom yade import pack,utils, plot,ymport,qt

W=working, starting a script with "rom yade..." surely is not :-)

> for i in O.interactions:
> ... O.forces.f(i)

O.forces.f needs body ID, see documentation [1]. In your code you pass "i", which is Interaction... (this is what the error says)

O.forces.f returns resulting force on the body. It has no relation to individual interaction forces.
(O.forces.f is sum of all interaction forces on the body, but with no backwards influence)

> sum((O.forces.f(i).norm()))

force.norm() is one number.
What you expect from sum(number)?
try it in python. e.g. sum(2)

> for i in O.interactions:
> if ...:
> Force = sum((O.forces.f(i).norm()))
> else :
> Force = 0

this way, you loop over all interactions just to assign Force to be the force from the very last interaction or 0...

A possible solution:
###
def addPlotData():
    Force = 0.0 # initial zero value
    for i in O.interactions:
        if isinstance(O.bodies[i.id1].shape,Facet) or isinstance(O.bodies[i.id2].shape,Facet):
            Force += i.phys.normalForce.norm() # increment the Force
    # or even a one liner
    # Force = sum(i.phys.normalForce.norm() for i in O.interactions if isinstance(O.bodies[i.id1].shape,Facet) or isinstance(O.bodies[i.id2].shape,Facet))
    plot.addData(i=O.iter,Force=Force)
###

You can save the data with plot.saveDataTxt [2], either inside addPlotData, or in some other PyRunner with lower frequency, or after the simulation is finished, or ... whenever you need.

cheers
Jan

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ForceContainer.f
[2] https://yade-dem.org/doc/yade.plot.html#yade.plot.saveDataTxt

Hi, Jan thank you for your response;

What can i do, if i want to save de force on a specific point on the mill?

Can you help me ?

What is the problem? getting the value? saving the value to a file? ... ?
What is "a specific point"?
cheers
Jan

Thank's Jan,
Plaese find the raison why i need to save the forces (Forces&Torques) at the specific point.
During the operation of a ball mill, the drum lifts the balls along one side of the mill until it reaches a falling point. From this position the particles begin to fall independently of the movement of the mill. The balls then fall onto the inner surface of the mill drum, creating impact forces (forces and torques). It is these forces that I wish to save.However, there is a region of the mill where the magnitude of the forces is maximum around the circumference of the mill. Exactly between 180 and 360 degrees in the third and fourth quadrant, if we consider that the mill rotates anti-clockwise. It is in this zone that I wish to record the efforts (forces and torque), and trace the evolution of its efforts in space (angle in degrees) and in time (time in seconds). To know exactly at which angular positions of the drum to place different sensors.

But I don't know yet how I can do it.

Thank's Jan,
Plaese find the raison why i need to save the forces (Forces&Torques) at the specific point.
During the operation of a ball mill, the drum lifts the balls along one side of the mill until it reaches a falling point. From this position the particles begin to fall independently of the movement of the mill. The balls then fall onto the inner surface of the mill drum, creating impact forces (forces and torques). It is these forces that I wish to save.However, there is a region of the mill where the magnitude of the forces is maximum around the circumference of the mill. Exactly between 180 and 360 degrees in the third and fourth quadrant, if we consider that the mill rotates anti-clockwise. It is in this zone that I wish to record the efforts (forces and torque), and trace the evolution of its efforts in space (angular position in degrees) and in time (time in seconds). To know exactly at which angular positions of the drum to place different sensors.

But I don't know yet how I can do it.

Hello,
thanks for more information.
If I got it correctly, there are several approaches, depending on your actual needs. Basically all of them are consists of:
1) getting all forces and contact points (action points of the forces)
2) do some postprocessing (averaging, smoothing, ...)

1) must be done in Yade, 2) may be done in the Yade script or somehow externally from saved data.

I **personally** would save the raw data from Yade and do external postprocessing (probably using Python or some dedicated tool).
This way, you can do the postprocessing from the same data with different parameters / methods.

The postprocessing may be both in space and time sense, summing / averaging / smoothing / ... , really many options, depending what you want..

cheers
Jan