Extract Normals and Shears/tangentials forces on a ballMill

Hello,

please read [1] and provide a MWE.

> I want to extract the normal and/or tangential forces impacted by the balls on the tumbling mill

please be more specific. Forces also with acting points? Just forces (no torques)? ...

> using two laws:

in two different simulations (then there would be no difference) or you want to distinguish contributions from the two laws in one simulation?
(would have been clear from the MWE :-)

extracting forces is easy, just please let us know your requirements in more details.

thanks for info
cheers
Jan

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

Hi Jan,

Thank you for answering my question, and for documentation on How to ask questions.

I'm sorry for my unclear question. I am a new user of Yade; I will make an effort to clarify my questions next time.

Now I will give clarification, based on your reactions to my question:

1. Please be more specific. Forces also with acting points? Just forces (no torques)? ...

 - Precisely, i want to extract the generalized forces (Forces and torques) with their points of contact. (on the tumbling)

2. In two different simulations (then there would be no difference) or you want to distinguish contributions from the two laws in one simulation?

- Yes, I need to distinguish contributions of the both law, compare them, and see which one is more effective than the other.

I hope i have clarified my question with the answers above:

thanks

Gauth

> - Yes, I need to distinguish contributions of the both law, compare them, and see which one is more effective than the other.

do you have
- both laws in one simulation
or
- you run separate simulations, each completely with one law
?

you run separate simulations, each completely with one law ?

No, I'll run the simulations separately and compare the results in an a post-processing program.

just to be sure :-)
1st simulation, entirely with Law2_ScGeom_MindlinPhys_Mindlin
2nd simulation, entirely Law2_ScGeom_ViscElPhys_Basic
is that correct?
thanks
Jan

Exactly!

Yes!

The 1st simulation, entirely with Law2_ScGeom_MindlinPhys_Mindlin
and
The 2st Simulation, entirely with Law2_ScGeom_ViscElPhy_Basic

separatly.

The specific use depends on "extract"/"record" (print, save, save format, ...), when (periodically, just at the end of simulation, ...), the simulation itself (mill is made of facets? boxes? ...?) etc. etc. (please provide a MWE [1])

But the core is to loop over all interactions, skip those you are not interested in and get the desired info from the interactions.

###
for i in O.interactions:
   if isNotInteresting(i):
      continue
   fn = i.phys.normalForce
   fs = i.phys.shearForce
   cp = i.geom.contactPoint
   # store, print, save, ..........
###

cheers
Jan

 Hi Jan
Thank you for your last response. To be clear I complete somes informations to allow you to understand my project.
I'd like to controls the simulation of 0.5 m diameter and 0.3 m length laboratory ball mill. The aim is to extract and save at each time step of simulations the generalized forces(forces and torques), provided by interactions between facets and spheres, their contact point, and the overlap at each contact point. The shell was designed using CAD software and its mesh size file was obtained using GMSH software. and can be imported using ballmill.mesh in GMSH format.
I need to save periodically the informations , The mill is made of facets, and i need to save in the Txt format, which can be used with other softwares.

Here is my code
# DATA
# Geometry of the mill and Filling of the balls
length = 0.3 # [m] mill length
rMean = 0.75/2.0 # [m] profile mean radius (only for filling here)
rBall = 0.0075 # [m] ball radius
nBalls = 20000 # [-] number of balls

# Simulation
omega = 3.1416 # [rad/s] angular rotation speed
tEnd = 10 # [s] total simulation time
O.dt = 5.0e-6 # [s] fixed time step
tRec = 2.0e-2 # [s] recording interval

# imports of the modules necessary for reading the mesh

from yade import ymport,export
import math,os

if os.path.exists('output'):
  print 'Delete old "output" directory.'
  exit()
else:
  os.makedirs('output')

print '\nBall Mill Simulation\n'

# Define the material of the mill and the balls.

millMat = O.materials.append(FrictMat(young=600.0e6,poisson=0.6,density=2.6e3,
                                  frictionAngle=radians(26),label='Friction'))

ballMat = O.materials.append(FrictMat(young=600.0e6,poisson=0.6,density=2.6e3,
                                  frictionAngle=radians(26),label='Friction'))

# Definition of geometry by importing to the simulation the mill mesh in the GMSH format ".mesh".

mill = O.bodies.append(ymport.gmsh('mill.mesh',material=millMat,color=(1,1,1),
                                    scale=1e-3))
walls = O.bodies.append([
          utils.wall((0,0,0),1,material=millMat,color=(1,1,1)),
          utils.wall((0,length,0),1,material=millMat,color=(1,1,1))
        ])

#The mill is filled by insering balls in horizontal planes which are
#progressively moving upwards until all balls have been added to the
#simulation.
#

# sweep gaps and initial coordinates
gap1 = 0.1*rBall # small gap to have no overlap
gap2 = 4*rBall # large gap to have no overlap (due to liner profile)
z = -rMean + rBall + gap2 # vertical coordinate
x = -math.sqrt((rMean-rBall)**2 - z**2) + gap2 # horizontal coordinate
y = rBall + gap1 # axial coordinate

for i in range(nBalls):
  if y > length - rBall - gap1:
    y = rBall + gap1
    x = x + 2*rBall + gap1
  if x > math.sqrt((rMean-rBall)**2 - z**2) - gap2:
    z = z + 2*rBall + gap1
    x = -math.sqrt((rMean-rBall)**2 - z**2) + gap2
  O.bodies.append(utils.sphere(center=(x,y,z),radius=rBall,\
                               material=ballMat,color=(1,1,1)))
  y = y + 2.0*rBall + gap1

print 'Number of balls: '+str(nBalls)+' [-]'

# the simulation is prepared by assembling the engines which
# will be successively called during the simulation loop.
#

rpm = omega*30/math.pi # [rpm] rotation speed

print 'Rotation speed: '+str(rpm)+' [rpm]'

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()] # The Hertz mindlin law or
  ),
  RotationEngine(ids=mill,rotationAxis=[0,1,0],rotateAroundZero=True,
                 angularVelocity=omega),
  NewtonIntegrator(damping=0,gravity=[0,0,-9.81]),

  VTKRecorder(fileName=os.getcwd()+'/output/',
              recorders=['spheres','facets','velocity'],
              iterPeriod=int(math.ceil(tRec/O.dt)))
]
The simulation is run until reaching its end time.

O.run(int(math.ceil(tEnd/O.dt))+1)

 What i need to do is, to add the code which can extract and save:

# - Generalized forces (forces and torque): normal and shear forces.
# - The contact points
# - The overlap corresponding of each contact point
# This can be done for each time-step of simulation
# The result can be saved in the files for post-processing

Hi

- normal and shear forces and contact points can be accessed like in previous answer #8

- interaction torques, additional to the ones caused by the contact forces, in the case of Law2_ScGeom_MindlinPhys_Mindlin seem to be at i.phys.momentBend and .momentTwist [0]. I'm not sure there exists interaction torques in Law2_ScGeom_ViscElPhys_Basic

- overlap is at [1], and is accessed with e.g. i.geom.penetrationDepth where i is an interaction of your simulation

- to repeat the measure and export at each time step, you may wanna look at PyRunner engine [2] and include one in your O.engines

- export in text files can be done following https://yade-dem.org/doc/user.html#exporting if you like gnuplot, or using plot.saveDataTxt()

- The whole doc § of https://yade-dem.org/doc/user.html#tracking-variables (and all User's manual doc) should be very helpful to complete such tasks

[0] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.MindlinPhys.momentBend and below
[1] https://yade-dem.org/doc/yade.wrapper.html?highlight=penetrationdepth#yade.wrapper.ScGeom.penetrationDepth
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.PyRunner
[3] https://yade-dem.org/doc/yade.plot.html?highlight=savedatatxt#yade.plot.saveDataTxt

Hello,

thanks for the code, but next time please try to create a MWE, where M stands for "minimal" - a few facets, no external mesh, just a few balls, no VTKRecoder... (sometimes, but very seldom, the full original code is needed).
Just focusing on your problem, which is saving same specific data.

> to extract and save at each time step

of course it is possible, but I would consider a longer period

> i need to save in the Txt format, which can be used with other softwares.

please be more specific what "Txt format" and "other softwares" are.
**I personally** prefer to save data in some "structured" format like JSON or YAML. Creating "plain text" from them is easy, the other way not.

I think the plot module is not suitable for this kind of exporting. I would save it "manually", see below.

### not tested, I am only willing to test it on a MWE :-)
import json
...
def exportForces():
   data = [] # list of data for each interaction
   for i in O.interactions:
      b1,b2 = [O.bodies[i] for i in (i.id1,i.id2)]
      # skip interactions not containing facets
      if not (isinstance(b1.shape,Facet) or isinstance(b2.shape,Facet)):
         continue
      # gather desired values
      fn = i.phys.normalForce
      fs = i.phys.
      cp = i.geom.contactPoint
      pd = i.geom.penetrationDepth
      # add it to the data collection
      d = dict(normalForce=fn,shearForce=fs,contactPoint=cp,penetrationDepth=pd,id1=i.id1,id2=i.id2)
     data.push(d)
   # save
   fName = "data{:06d}.json".format(O.iter)
   with open(fName,"w") as f:
      json.dump(data,f,indent=2)
...
O.engines = [
   ...
   PyRunner(iterPeriod=1,command="exportForces()"),
]
...
###

cheers
Jan

Hi Jan, and Jerome,

Thank'you for your advices. It's been very helpful. However, I followed the recording method proposed by Jan, the program runs, but i can't get my data recorded, and i get the following error, some times after launching the simulation.

AttributeError Traceback (most recent call last)
/usr/bin/yade in <module>()

/usr/bin/yade in exportForces()
    119 # add it to the data collection
    120 d = dict(normalForce=fn,shearForce=fs,contactPoint=cp,penetrationDepth=pd,id1=i.id1,id2=i.id2)
--> 121 data.push(d)
    122 # save
    123 fName = "data{:06d}.json".format(O.iter)

AttributeError: 'list' object has no attribute 'push'

# this is the MWE i used:

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()]
  ),
  RotationEngine(ids=mill,rotationAxis=[0,1,0],rotateAroundZero=True,angularVelocity=omega),
  NewtonIntegrator(damping=0,gravity=[0,0,-9.81]),
  VTKRecorder(fileName=os.getcwd()+'/output/', recorders=['spheres','facets','velocity'],iterPeriod=int(math.ceil(tRec/O.dt))),
  PyRunner(command='exportForces()',iterPeriod=int(math.ceil(tRec/O.dt))),

]

import json

def exportForces():
   data = [] # list of data for each interaction
   for i in O.interactions:
      b1,b2 = [O.bodies[i] for i in (i.id1,i.id2)]
      # skip interactions not containing facets
      if not (isinstance(b1.shape,Facet) or isinstance(b2.shape,Facet)):
         continue
      # gather desired values
      fn = i.phys.normalForce
      fs = i.phys.shearForce
      cp = i.geom.contactPoint
      pd = i.geom.penetrationDepth
      # add it to the data collection
      d = dict(normalForce=fn,shearForce=fs,contactPoint=cp,penetrationDepth=pd,id1=i.id1,id2=i.id2)
   data.push(d)
   # save
   fName = "data{:06d}.json".format(O.iter)
   with open(fName,"w") as f:
      json.dump(data,f,indent=2)

# The simulation is run until reaching its end time.

O.run(int(math.ceil(tEnd/O.dt))+1

data.push(d)
should have been
data.append(d)

Hi Jan,

Sorry, I replace data.push(d) y data.append(d), however, the data file thoes'nt appear any where.

But i'm trying to use this code to save just fn separately
def exportForces():
fn = 0 # initialization of the normal force
  for i in O.interactions:
    # detection of Force by sphere on facet
    if isinstance(O.bodies[i.id1].shape,Sphere) \
        and isinstance(O.bodies[i.id2].shape,Facet):
      fn = fn + i.phys.normalForce
  out=open(os.getcwd()+'/output/normalForce','a')
  out.write('%g\n'%(fn))
  out.close()
O.engines=[
.......
PyRunner(command='exportForces()',iterPeriod=int(math.ceil(tRec/O.dt)))

......]
The normalForce file appear in the folder of simulation, but the value of fn extracted at any time step is always 0.

What do you think?

i don't know if there's a command missing in the code, that would activate the force?

Hello,

please provide a complete MWE [1]. There might be maaany reasons, it seems to be necessary to be able to try the code to solve the problem.

> PyRunner(command='exportForces()',iterPeriod=int(math.ceil(tRec/O.dt))),

just use iterPeriod=1 if you want the export each iteration (which I understood so far).

> The normalForce file appear in the folder of simulation, but the value of fn extracted at any time step is always 0.
> i don't know if there's a command missing in the code, that would activate the force?

there is no special command (if your simulation works).
But is you do not export each iteration, it is possible that the force just is zero..

> def exportForces():
> fn = 0 # initialization of the normal force
> ...
> fn = fn + i.phys.normalForce

this should be error, fn=0 is int and you cannot add i.phys.normalForce (Vector3) to it.
And as it is given here is syntax error..

waiting for MWE :-)
cheers
Jan

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

Hi Jan,

Thank you for the response. This is a MWE i propose for a better understanding of my problem.
As it is not possible to provide the mill mesh file in gmsh format, The mill is simulated here as a cylinder of the same size, rotating around z-axis. The balls are added to the mill using sphere pack, (p=pack.SpherePack());
However the behavior is not quite what I expected, especially since the mill is not equipped with a lifter to make the balls circulate in the mill as it is physically done. However there are interactions with the mill that normally have to produce normal and tangential forces, That i can't save actually.

I also noticed that the increased rigidity, some balls pass through the walls of the mill, instead of remaining inside.*
Moreover, I expected that at the beginning of the simulation, packages of balls could fall along the y-axis, to disintegrate afterwards. But i notice that it already disintegrates before falling, and go in the whole direction of mill, despite the acceleration of gravity is defined according to the y-axis.*
*
Finaly find the same recording problem mentioned above.

This is the MWE.

# Material
rho = 7640.0 # [kg/m^3] density (of the balls)
Eyoung = 900.0e9 # [N/m] equivalent normal stiffness
vp = 0.6 # [N/m] equivalent tangential stiffness
en = 0.3 # [-] normal coefficient of restitution (sphere/sphere)
mu = 2 # [-] friction coefficient
rMean = 0.8/2.0
length = 0.4
rBall = 0.0075 # [m] ball radius
sphRadFuzz = 0.8

# Simulation
pcv = 70.0 # [%] percentage of critical velocity
tEnd = 10.0 # [s] total simulation time
O.dt = 1.0e-5 # [s] fixed time step
tRec = 1.0e-2 # [s] recording interval
#material of cylinder
cylMat = O.materials.append(FrictMat(poisson=vp,young=Eyoung,density=rho,frictionAngle=radians(mu)))

#material of spheres
ballMat = O.materials.append(FrictMat(poisson=vp,young=Eyoung,density=rho,frictionAngle=radians(mu)))

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))
#walls = O.bodies.append([
# utils.wall((0,0,length/2),2,material=Mat1,color=(1,1,1),sense=+1),
# utils.wall((0,0,-length/2),2,material=Mat1,color=(1,1,1),sense=+1)
# ])
# 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():
  # Mechanical power
  data = [ ] # List of data for each interaction
  for i in O.interactions:

    # Force by sphere on facet
    if isinstance(O.bodies[i.id1].shape,Sphere) \
        and isinstance(O.bodies[i.id2].shape,Facet):
      fn = i.phys.normalForce
      fs = i.phys.shearForce
      cp = i.geom.contactPoint
      pd = i.geom.penetrationDepth
      d = dict(normalForce=fn,shearForce=fs,contactPoint=cp,penetrationDepth=pd, id1=i.id1,id2=i.id2)
      data.append(d)
      print(data)

      # Save data
      dataSim = "data{:06d}.json".format(O.iter)
      with open(DataSim,"w") as f:
        json.dump(data,f,indent=2)

O.dt = 1e-6
###
O.run(int(math.ceil(tEnd/O.dt))+1)

Thank's

Gauthier

> I also noticed that the increased rigidity, some balls pass through the walls of the mill, instead of remaining inside.
> O.dt = 1e-6

use time step according to the critical value, e.g. a fraction of PWaveTimeStep [1]

> Moreover, I expected that at the beginning of the simulation, packages of balls could fall along the y-axis, to disintegrate afterwards. But i notice that it already disintegrates before falling, and go in the whole direction of mill, despite the acceleration of gravity is defined according to the y-axis.

> for i in (0,1): # red and blue spheres
> sp=pack.SpherePack();
> ...
> sp.makeCloud(...)
> O.bodies.append([utils.sphere(s[0],s[1],...) for s in sp])

Here you use two SpherePack, their makeCloud does not know anything about each other and create overlapping cloud, resulting in "disintegration".
Use just one SpherePack:
###
sp=pack.SpherePack() # (!) one instance before the for loop
for i in (0,1): # red and blue spheres
 ...
 sp.makeCloud(...)
 O.bodies.append([utils.sphere(s[0],s[1],...) for s in sp])
###

> there are interactions with the mill that normally have to produce normal and tangential forces, That i can't save actually.

how do you know?
there should be != there are
I have tested your code and simply there were no interactions (for relatively long time for testing), so no forces were saved..
You can try something like "print len([i for i in O.interactions])" in PyRunner to be sure.
Always try to get numerical evidence for such statements.

========
To your exportForces function:

> if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Facet):

I think the order is different, that Facet should be id1. anyway, do not depend on the order of interaction bodies

> with open(DataSim,"w") as f:

DataSim is not defined (dataSim is)

the saving code has wrong indentation, it should be on the "main" level, not inside "for i in O.interactions" cycle

======
The MWE was not real MWE (e.g. materials is irrelevant here, the amount of spheres is too high for playing and testing, some comments were not comments but just commented unused code ...).
Next time please really try to isolate the specific problem, in this case saving sphere-facet interactions

below modified and reduced MWE:
###
import json

rMean = 0.8/2.0
length = 0.4
rBall = 0.075

# Simulation
pcv = 70.0 # [%] percentage of critical velocity
tEnd = 10.0 # [s] total simulation time
O.dt = 1.0e-5 # [s] fixed time step
tRec = 1.0e-2 # [s] recording interval

ids=O.bodies.append(geom.facetCylinder((0,0,0),0.8,0.4,segmentsNumber=64, wallMask=7, angleRange=None,
                   radiusTopInner=-1, radiusBottomInner=-1))
# 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)
sp=pack.SpherePack();
for i in (0,1): # red and blue spheres
 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,0,999999,False)
 O.bodies.append([utils.sphere(s[0],s[1],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()]
    ),
   RotationEngine(ids=ids,angularVelocity=1000,rotateAroundZero=True,rotationAxis=(0,0,1)),
   NewtonIntegrator(damping=0,gravity=[0,0,-100]),
   PyRunner(command='exportForces()',iterPeriod=1)
]
def exportForces():
  # Mechanical power
  data = [ ] # List of data for each interaction
  for i in O.interactions:

    # Force by sphere on facet
    if isinstance(O.bodies[i.id1].shape,Facet) or isinstance(O.bodies[i.id2].shape,Facet):
      fn = tuple(i.phys.normalForce)
      fs = tuple(i.phys.shearForce)
      cp = tuple(i.geom.contactPoint)
      pd = i.geom.penetrationDepth
      d = dict(normalForce=fn,shearForce=fs,contactPoint=cp,penetrationDepth=pd, id1=i.id1,id2=i.id2)
      data.append(d)
      print(data)

  # Save data
  dataSim = "data{:06d}.json".format(O.iter)
  with open(dataSim,"w") as f:
    json.dump(data,f,indent=2)

O.dt = .5*PWaveTimeStep()
O.run(200,True)
###

cheers
Jan

[1] https://yade-dem.org/doc/yade.utils.html#yade._utils.PWaveTimeStep

Hi Jan, and Jêrome

 Thank you for your contributions; This ones gives satisfaction of my project.

regard,

Gauthier

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

Dear All,

With this piece of code, I can calculate for each interaction the normal impact force between a ball and the cylinder at each time step. The result is such that, by recording the time vector with time = O.time, the code registers the time of interaction. This is a good thing. However if I have ten interactions at the instant 0.0502s for instance, O.time, gives me a vector with 0.0502s repeated 10 times. This makes it difficult for me to calculate the spectrum of my signal with such a time vector.
What I would like to improve in my code is :
Rather than recording for each interaction the normal force and the time of interaction, I would rather record for each time step, the number of interactions, as well as the total normal force due to all its interactions in a time t.
 That way, I would have an easier signal to interpret.

Someone could help me to improve my code?

#######
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()]
  ),

  NewtonIntegrator(damping=0,gravity=[0,0,-9.81]),

  RotationEngine(ids=liner,rotationAxis=[0,-1,0],rotateAroundZero=True,
                 angularVelocity=omega),
  PyRunner(command='saveForces()',iterPeriod=1)
]

def saveForces():

  for i in O.interactions:

    if isinstance(O.bodies[i.id1].shape,Facet) or isinstance(O.bodies[i.id2].shape,Facet):
      fn = tuple((i.phys.normalForce))
      time = O.time

    # Save Total force
      out=open(os.getcwd()+'/output/normalForce','a')
      out.write('%s\n'%str(fn))
      out.close()

O.dt = .95*PWaveTimeStep()
O.run(int(math.ceil(tEnd/O.dt))+1)

Hello,

> I can calculate for each interaction the normal impact force between a ball and the cylinder at each time step.
> I would rather record for each time step, the number of interactions, as well as the total normal force due to all its interactions in a time t.

if you have data for all interaction, what is the problem to get overall data?

cheers
Jan

The problem is that I'm not sure that a simple somation for each interaction force will give me the result of the total forces of the interactions during each iteration (time step).

When I use the sum function, instead of getting the total sum of interactions, it calculates the sum of the components (x,y,z) in a vector; this is not at all what I'm looking for.
I don't know if there is a way to insert a loop over the itérations (time step), in which oui can extract the resultant of the interaction forces for each iteration.

You have access to the interactions and forces both in Python and in the saved files, you can do whatever you want with them.. getting number of them as well as the summation seems like trivial task..

Yes, simple summation of the forces will give you the result of the total force, which is the sum of the (x,y,z) components..
If it is not what you are looking for, then please describe more in detail what you are looking for.

Also next time please provide the code, there might be some misunderstanding in the words.

in Python:
###
# extract relevant interactions
inrs = [i for i in O.interactions if isinstance(O.bodies[i.id1].shape,Facet) or isinstance(O.bodies[i.id2].shape,Facet)]
# get data
numberOfInteractions = len(intrs)
totalNormalForce = sum((i.phys.normalForce for i in intrs), Vector3.Zero)
totalNormalForceMagnitude = sum(i.phys.normalForce.norm() for i in intrs)
###

cheers
Jan