How to detect the coordinates of each particle?

Hello,

> How do I get the coordinates of each particle?

it depends on the definition of "get".

You can print them in a for loop:
for b in O.bodies: print(b.state.pos)

You can get the coordinates e.g. as a list:
coords = [b.state.pos for b in O.bodies]

You can save them as a file:
yade.export.text("some_file_name.txt")

You can ...

> how do I know the z-value of the uppermost particle?

zMax = max(b.state.pos[2] for b in O.bodies)

Cheers
Jan

I tried it and found that the code had no errors, but no z-value output for the uppermost particles.
My code is as follows:

from yade import pack,plot
O.bodies.append(geom.facetBox((0.05, 9, 10), (0.1, 9, 10), wallMask=63))
sp = pack.SpherePack()
sp.makeCloud(Vector3(0.05,0,0),Vector3(0.05,18,20), rMean=0.1, rRelFuzz=1/3)
sp.toSimulation()

(xdim,ydim,zdim)= aabbDim()

print("Height is ",zdim)

for b in O.bodies:
   if isinstance(b.shape,Sphere):
       b.state.blockedDOFs='ZxY'
       b.shape.color=(3.,2.,1.)

circleRadius=1.5
circleCenter = Vector3(0.05,6,6)
circleRadius1=1.5
circleCenter1 = Vector3(0.05,12,6)
   #
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
        InteractionLoop(
                # handle sphere+sphere and facet+sphere collisions
                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.4),
        # call the checkUnbalanced function (defined below) every 2 seconds
        PyRunner(command='checkUnbalanced()', realPeriod=2),
        # call the addPlotData function every 200 steps
        PyRunner(command='addPlotData()', iterPeriod=100)
]
O.dt = 0.5 * PWaveTimeStep()

# enable energy tracking; any simulation parts supporting it
# can create and update arbitrary energy types, which can be
# accessed as O.energy['energyName'] subsequently
O.trackEnergy = True

# if the unbalanced forces goes below .05, the packing
# is considered stabilized, therefore we stop collected
# data history and stop
def checkUnbalanced():
 if unbalancedForce() < .05:
  O.pause()
  plot.saveDataTxt('bbb.txt.bz2')
                for b in O.bodies:
      if isinstance(b.shape,Sphere):
                        #b.state.blockedDOFs='zxy'
                        b.state.vel=(0,0,0)
                        b.state.angVel=(0,0,0)
                    d = (b.state.pos - circleCenter).norm() # distance of circleCenter and center of "b"
                    d1 = (b.state.pos - circleCenter1).norm() # distance of circleCenter and center of "b"
                    if d < circleRadius:
                       O.bodies.erase(b.id)
                    if d1 < circleRadius1:
                       O.bodies.erase(b.id)
      #coords = [b.state.pos for b in O.bodies]
                    #print(b.state.pos)
   #zMax = max(b.state.pos[2] for b in O.bodies)

  O.bodies.append(sphere(center=(0.05,7.38276367683319,6.11900825105632), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.33147912281321,6.33262379212493), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.24740913386372,6.53558669963537), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.13262379212493,6.72289935320946), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.98994949366117,6.88994949366117), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.82289935320946,7.03262379212493), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.63558669963537,7.14740913386372), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.43262379212493,7.23147912281322), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.21900825105632,7.28276367683319), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6,7.3), radius=.1019, fixed=True,color=(1.,1.,1.)))

  O.bodies.append(sphere(center=(0.05,5.78099174894368,7.28276367683319), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.56737620787507,7.23147912281322), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.36441330036463,7.14740913386372), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.17710064679054,7.03262379212493), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.01005050633883,6.88994949366117), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.86737620787507,6.72289935320946), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.75259086613629,6.53558669963537), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.66852087718679,6.33262379212493), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.61723632316681,6.11900825105632), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.6,5.9), radius=.1019, fixed=True,color=(1.,1.,1.)))

  O.bodies.append(sphere(center=(0.05,4.61723632316681,5.68099174894368), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.66852087718678,5.46737620787507), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.75259086613629,5.26441330036464), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,4.86737620787507,5.07710064679054), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.01005050633883,4.91005050633883), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.17710064679054,4.76737620787507), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.36441330036463,4.65259086613629), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.56737620787507,4.56852087718679), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,5.78099174894368,4.51723632316681), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6,4.5), radius=.1019, fixed=True,color=(1.,1.,1.)))

  O.bodies.append(sphere(center=(0.05,6.21900825105632,4.51723632316681), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.43262379212493,4.56852087718679), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.63558669963537,4.65259086613629), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.82289935320946,4.76737620787507), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,6.98994949366117,4.91005050633883), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.13262379212493,5.07710064679054), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.24740913386372,5.26441330036464), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.33147912281321,5.46737620787507), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.38276367683319,5.68099174894368), radius=.1019, fixed=True,color=(1.,1.,1.)))
  O.bodies.append(sphere(center=(0.05,7.4,5.9), radius=.1019, fixed=True,color=(1.,1.,1.)))
  zMax = max(b.state.pos[2] for b in O.bodies)

  (xdim,ydim,zdim)= aabbDim()

  print("Height is ",zdim)

# collect history of data which will be plotted
def addPlotData():
 # each item is given a names, by which it can be the unsed in plot.plots
 # the **O.energy converts dictionary-like O.energy to plot.addData arguments
 plot.addData(i=O.iter, unbalanced=unbalancedForce(), **O.energy)

# define how to plot data: 'i' (step number) on the x-axis, unbalanced force
# on the left y-axis, all energies on the right y-axis
# (O.energy.keys is function which will be called to get all defined energies)
# None separates left and right y-axis

plot.plots = {'i': ('unbalanced', None, O.energy.keys)}

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

O.saveTmp()

# to see it
from yade import qt
qt.Controller()
qt.View()
qt.View()

Cheers

Hello,

> My code is as follows:

please make your code MWE [1]:
W = working. This code throw IndentationError: unexpected indent
M = minimal. For the problem, there is no need of so many O.bodies.append, no qt is needed, plot is not needed, plenty of commented lines are not needed, etc. etc.

> I tried it and found that the code had no errors, but no z-value output for the uppermost particles.
> zMax = max(b.state.pos[2] for b in O.bodies)

what do you mean by "no z-value output"?
zMax is not used anywhere (e.g. not printed), is this the problem?

Cheers
Jan

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

My code is as follows:

from yade import pack,plot
O.bodies.append(geom.facetBox((0.05, 9, 10), (0.1, 9, 10), wallMask=63))
sp = pack.SpherePack()
sp.makeCloud(Vector3(0.05,0,0),Vector3(0.05,18,20), rMean=0.1, rRelFuzz=1/3)
sp.toSimulation()
(xdim,ydim,zdim)= aabbDim()

print("Height is ",zdim)

for b in O.bodies:
   if isinstance(b.shape,Sphere):
       b.state.blockedDOFs='ZxY'
       b.shape.color=(3.,2.,1.)

O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.4),
        PyRunner(command='checkUnbalanced()', realPeriod=2),
        PyRunner(command='addPlotData()', iterPeriod=100)
]
O.dt = 0.5 * PWaveTimeStep()
O.trackEnergy = True

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

      #coords = [b.state.pos for b in O.bodies]
                    #print(b.state.pos)
   #zMax = max(b.state.pos[2] for b in O.bodies)

  zMax = max(b.state.pos[2] for b in O.bodies)
  print("Z ",zMax)
  (xdim,ydim,zdim)= aabbDim()

  print("Height is ",zdim)

def addPlotData():
 plot.addData(i=O.iter, unbalanced=unbalancedForce(), **O.energy)

O.saveTmp()
from yade import qt
qt.Controller()
#qt.View()

>This code throw IndentationError: unexpected indent

This may be due to an indentation error caused by copying over.

>what do you mean by "no z-value output"?

Yes, I found this problem, I did not set the print for output.

But the zMax = max(b.state.pos[2] for b in O.bodies) string of code didn't solve my problem. It can only get the z-value of the highest particle before the run, and I need to get the z-value of the uppermost particles after gravity.

Thanks for the MWE.

> But the zMax = max(b.state.pos[2] for b in O.bodies) string of code didn't solve my problem. It can only get the z-value of the highest particle before the run, and I need to get the z-value of the uppermost particles after gravity.

the command gets the zMax value at the time it is executed.
If you execute it before the run, it is evaluated before the run.
If you execute it "after gravity", it is evaluated "after gravity".
Here in your code, it is evaluated inside checkUnbalanced() function if unbalancedForce() < .05 is fulfilled.

Also, max(b.state.pos[2] for b in O.bodies) takes into account all bodies, including facets.
If you are interested only in spheres, you can modify it as:
zMax = max(b.state.pos[2] for b in O.bodies if isinstance(b.shape,Sphere))

Cheers
Jan

Yes, I added zMax = max(b.state.pos[2] for b in O.body if isinstance(b.shape, Sphere)) to the checkUnbalanced() function, z-values can get the output, but only one zMax output, zMAX is the highest z-value, and I need to output the z-value of each particle in the uppermost layer.

> ... coordinates of each particle?
> ... the z-value of the uppermost particle?
> ... the z-value of each particle in the uppermost layer.

please be consistent.

> in the uppermost layer.

first, define rigorously "the uppermost layer", i.e. how mathematically determine if particle belongs to this layer or not.
Then just:
###
def isUpperMostLayer(body):
    ... # your definition
    return False # or True

zValsUpperMostLayer = [b.state.pos[2] for b in O.bodies if isinstance(b.shape, Sphere) and isUpperMostLayer(b)]
print(zValsUpperMostLayer) # or whatever else with the values
###

Cheers
Jan

I see what you mean. But how can I strictly define particles as belonging to the top layer, I don't seem to have a train of thought.
My first thought was to use coordinates to define the uppermost particles, but the coordinates of these particles are also unknown quantities, and I can't use unknown quantities to solve for unknown quantities.

Cheers

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