Measuring Force,Trapdoor displacement and Surface displacement

Hello,

1,2,3)
to get force on one body, use
O.forces.f(id)

to get sum from several bodies, use
sum((O.forces.f(id) for id in ids), Vector3.Zero)
or
f = Vector3.Zero
for i in ids:
   f += O.forces.f(i)

for the values of left, right and middle, you only have to pass correct ids of corresponding facets
in the specific code example:
middle: 0,1
left: 2,3
right: 4,5

4)
displacement at given time step:
trapdoor.state.displ()
where trapdoor is one of the facets of the "middle square" (e.g. O.bodies[0])
to record the values over time, you can use e.g. [1]

5)
you have access to the positions and orientations of each body and its vertices
##
for b in O.bodies:
   if not isinstance(b.shape,Polyhedra): continue
   pos = b.state.pos # center of polyhedron
   ori = b.state.ori
   vs = b.shape.v
   vs = [pos + ori*v for v in vs] # vertices of polyhedron in global coordinate system
##
you can store the values at the start of simulation and compare it with state during simulation.
then it is up to your definition of surface displacement

cheers
Jan

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

thank you

I tried to put all together but got a lot of errors. I am new to python and yade now still learning but have to complete this project for my these, I will be very grateful if you can help me identify and correct the errors, please. Thanks

###########
######################################################################
# INPUTS
######################################################################
gravity = 100

# box dimensions
widthl = .3
widthr = .3
widthc = .3
height = .3
thick = .1
deep = -.2

# size of grains
sizeMin = 40e-3
sizeMax = 60e-3

frictionAngle = .5
young = 1e8 # stiffness

dt = 1e-3 # time step

nGravityDeposition = 250 # how long to run initial gravity deposition
nCycles = 3 # how many jumps to run afterwards
nStepsBetweenCycles = 200 # number of time steps between jumps
dspl = 20e-3

# how much larger the initial make cloud box should be
fillBoxHFactor = 3
######################################################################
from yade import polyhedra_utils

width = widthl+widthc+widthr

# mat, engines, ...
mat = PolyhedraMat(young=young,poisson=10,frictionAngle=frictionAngle)
O.materials.append(mat)
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
      [Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
      [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()],
      [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()],
   ),
   NewtonIntegrator(damping=0.4,gravity=(0,0,-gravity)),
 PyRunner(iterPeriod=1,command='checker()'),
]
O.dt = dt
def checker():
 for i in range(nCycles):
  ii = nGravityDeposition+i*nStepsBetweenCycles
  if O.iter == ii:
   moveBottom()
  if O.iter == ii+1:
   stopBottom()
 if O.iter == nGravityDeposition+nCycles*nStepsBetweenCycles:
  O.pause()
def moveBottom():
 v = dspl / O.dt
 for b in movables:
  b.state.vel = (0,0,-v)
def stopBottom():
 for b in movables:
  b.state.vel = (0,0,0)

# box
p000 = Vector3(0,0,0)
p100 = Vector3(widthl,0,0)
p200 = Vector3(widthl+widthc,0,0)
p300 = Vector3(widthl+widthc+widthr,0,0)
pxs = (p000,p100,p200,p300)
p001,p101,p201,p301 = [p+Vector3(0,0,height) for p in pxs]
p010,p110,p210,p310 = [p+Vector3(0,thick,0) for p in pxs]
p011,p111,p211,p311 = [p+Vector3(0,thick,height) for p in pxs]
p00b,p10b,p20b,p30b = [p+Vector3(0,0,deep) for p in pxs]
p01b,p11b,p21b,p31b = [p+Vector3(0,thick,deep) for p in pxs]
def rect(vs,**kw):
 v1,v2,v3,v4 = vs
 return [
  facet((v1,v2,v3),**kw),
  facet((v1,v3,v4),**kw),
 ]
movables = rect((p100,p200,p210,p110)) # bottom center
rects = (
 (p000,p100,p110,p010), # bottom left
 (p200,p300,p310,p210), # bottom left
 (p000,p010,p011,p001), # left
 (p300,p310,p311,p301), # right
 (p000,p100,p101,p001), # front left
 (p100,p200,p201,p101), # front center
 (p200,p300,p301,p201), # front right
 (p010,p110,p111,p011), # back left
 (p110,p210,p211,p111), # back center
 (p210,p310,p311,p211), # back right
 (p100,p200,p20b,p10b), # front center below
 (p110,p210,p21b,p11b), # back center below
 (p100,p110,p11b,p10b), # left below
 (p200,p210,p21b,p20b), # right below
)
rects = movables + sum((rect(r) for r in rects),[])
O.bodies.append(rects)

# gravel
polyhedra_utils.fillBox((0,0,0),(width,thick,fillBoxHFactor*height),mat,sizemin=3*[sizeMin],sizemax=3*[sizeMax],seed=1)
###########
#to record the values

O.forces.f(id)
sum((O.forces.f(id) for id in ids), Vector3.Zero)
#
#f = Vector3.Zero
#for i in ids:
# f += O.forces.f(i)

middle= 0,1
left = 2,3
right = 4,5

trapdoor.state.displ()
#O.bodies[1].state.displ()

##
for b in O.bodies:
   if not isinstance(b.shape,Polyhedra): continue
   pos = b.state.pos # center of polyhedron
   ori = b.state.ori
   vs = b.shape.v
   vs = [pos + ori*v for v in vs] # vertices of polyhedron in global coordinate system
##

please post your code here and also the errors you get
thanks
Jan

This is the code please :

###########
######################################################################
# INPUTS
######################################################################
gravity = 100

# box dimensions
widthl = .3
widthr = .3
widthc = .3
height = .3
thick = .1
deep = -.2

# size of grains
sizeMin = 40e-3
sizeMax = 60e-3

frictionAngle = .5
young = 1e8 # stiffness

dt = 1e-3 # time step

nGravityDeposition = 250 # how long to run initial gravity deposition
nCycles = 3 # how many jumps to run afterwards
nStepsBetweenCycles = 200 # number of time steps between jumps
dspl = 20e-3

# how much larger the initial make cloud box should be
fillBoxHFactor = 3
######################################################################
from yade import polyhedra_utils

width = widthl+widthc+widthr

# mat, engines, ...
mat = PolyhedraMat(young=young,poisson=10,frictionAngle=frictionAngle)
O.materials.append(mat)
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
      [Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
      [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()],
      [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()],
   ),
   NewtonIntegrator(damping=0.4,gravity=(0,0,-gravity)),
 PyRunner(iterPeriod=1,command='checker()'),
]
O.dt = dt
def checker():
 for i in range(nCycles):
  ii = nGravityDeposition+i*nStepsBetweenCycles
  if O.iter == ii:
   moveBottom()
  if O.iter == ii+1:
   stopBottom()
 if O.iter == nGravityDeposition+nCycles*nStepsBetweenCycles:
  O.pause()
def moveBottom():
 v = dspl / O.dt
 for b in movables:
  b.state.vel = (0,0,-v)
def stopBottom():
 for b in movables:
  b.state.vel = (0,0,0)

# box
p000 = Vector3(0,0,0)
p100 = Vector3(widthl,0,0)
p200 = Vector3(widthl+widthc,0,0)
p300 = Vector3(widthl+widthc+widthr,0,0)
pxs = (p000,p100,p200,p300)
p001,p101,p201,p301 = [p+Vector3(0,0,height) for p in pxs]
p010,p110,p210,p310 = [p+Vector3(0,thick,0) for p in pxs]
p011,p111,p211,p311 = [p+Vector3(0,thick,height) for p in pxs]
p00b,p10b,p20b,p30b = [p+Vector3(0,0,deep) for p in pxs]
p01b,p11b,p21b,p31b = [p+Vector3(0,thick,deep) for p in pxs]
def rect(vs,**kw):
 v1,v2,v3,v4 = vs
 return [
  facet((v1,v2,v3),**kw),
  facet((v1,v3,v4),**kw),
 ]
movables = rect((p100,p200,p210,p110)) # bottom center
rects = (
 (p000,p100,p110,p010), # bottom left
 (p200,p300,p310,p210), # bottom left
 (p000,p010,p011,p001), # left
 (p300,p310,p311,p301), # right
 (p000,p100,p101,p001), # front left
 (p100,p200,p201,p101), # front center
 (p200,p300,p301,p201), # front right
 (p010,p110,p111,p011), # back left
 (p110,p210,p211,p111), # back center
 (p210,p310,p311,p211), # back right
 (p100,p200,p20b,p10b), # front center below
 (p110,p210,p21b,p11b), # back center below
 (p100,p110,p11b,p10b), # left below
 (p200,p210,p21b,p20b), # right below
)
rects = movables + sum((rect(r) for r in rects),[])
O.bodies.append(rects)

# gravel
polyhedra_utils.fillBox((0,0,0),(width,thick,fillBoxHFactor*height),mat,sizemin=3*[sizeMin],sizemax=3*[sizeMax],seed=1)
###########
#to record the values

#O.forces.f(id)
#sum((O.forces.f(id) for id in ids), Vector3.Zero)
#
#f = Vector3.Zero
#for i in ids:
# f += O.forces.f(i)

middle = 0,1
lef= 2,3
right = 4,5

#trapdoor.state.displ()
O.bodies[1].state.displ()

##
for b in O.bodies:
   if not isinstance(b.shape,Polyhedra): continue
   pos = b.state.pos # center of polyhedron
   ori = b.state.ori
   vs = b.shape.v
   vs = [pos + ori*v for v in vs] # vertices of polyhedron in global coordinate system
##

And this is the error:

wise@wise-CM6330-CM6630-CM6730-CM6830-M11AA-8:~/wise$ yade nanawise.py
Welcome to Yade 2016.06a
TCP python prompt on localhost:9000, auth cookie `ssaduc'
XMLRPC info provider on http://localhost:21000
Running script nanawise.py
generated 136 polyhedrons
[[ ^L clears screen, ^U kills line. F12 controller, F11 3d view (use h-key for showing help), F10 both, F9 generator, F8 plot. ]]

Yade [1]: ERROR /build/yade-iRsyM_/yade-2016.06a/pkg/common/InsertionSortCollider.cpp:240 action: verletDist is set to 0 because no spheres were found. It will result in suboptimal performances, consider setting a positive verletDist in your script.
import yade.plot; yade.plot.plot();
^[[19~^[[19~^[[19~^[[19~---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
/usr/bin/yade in <module>()
----> 1 import yade.plot; yade.plot.plot();

/usr/lib/x86_64-linux-gnu/yade/py/yade/plot.py in plot(noShow, subPlots)
    593 global currLineRefs
    594 figs=set([l.line.axes.get_figure() for l in currLineRefs])
--> 595 if not hasattr(list(figs)[0],'show') and not noShow:
    596 import warnings
    597 warnings.warn('plot.plot not showing figure (matplotlib using headless backend?)')

IndexError: list index out of range

Yade [2]: import yade.plot; yade.plot.plot();
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
/usr/bin/yade in <module>()
----> 1 import yade.plot; yade.plot.plot();

/usr/lib/x86_64-linux-gnu/yade/py/yade/plot.py in plot(noShow, subPlots)
    593 global currLineRefs
    594 figs=set([l.line.axes.get_figure() for l in currLineRefs])
--> 595 if not hasattr(list(figs)[0],'show') and not noShow:
    596 import warnings
    597 warnings.warn('plot.plot not showing figure (matplotlib using headless backend?)')

IndexError: list index out of range

I meant complete code..

> Yade [1]: ERROR /build/yade-iRsyM_/yade-2016.06a/pkg/common/InsertionSortCollider.cpp:240 action: verletDist is set to 0 because no spheres were found. It will result in suboptimal performances, consider setting a positive verletDist in your script.
import yade.plot; yade.plot.plot();

this is actually no error

The other IndexErrors: I think you should define plot.plots [1,2] before plot.plot(), will check tomorrow with computer.

Jan

[1] https://yade-dem.org/doc/user.html#plotting-variables
[2] https://yade-dem.org/doc/yade.plot.html#yade.plot.plots

THIS IS THE COMPLETE CODE:

###########
######################################################################
# INPUTS
######################################################################
gravity = 100

# box dimensions
widthl = .3
widthr = .3
widthc = .3
height = .3
thick = .1
deep = -.2

# size of grains
sizeMin = 40e-3
sizeMax = 60e-3

frictionAngle = .5
young = 1e8 # stiffness

dt = 1e-3 # time step

nGravityDeposition = 250 # how long to run initial gravity deposition
nCycles = 3 # how many jumps to run afterwards
nStepsBetweenCycles = 200 # number of time steps between jumps
dspl = 20e-3

# how much larger the initial make cloud box should be
fillBoxHFactor = 3
######################################################################
from yade import polyhedra_utils

width = widthl+widthc+widthr

# mat, engines, ...
mat = PolyhedraMat(young=young,poisson=10,frictionAngle=frictionAngle)
O.materials.append(mat)
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
      [Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
      [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()],
      [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()],
   ),
   NewtonIntegrator(damping=0.4,gravity=(0,0,-gravity)),
 PyRunner(iterPeriod=1,command='checker()'),
]
O.dt = dt
def checker():
 for i in range(nCycles):
  ii = nGravityDeposition+i*nStepsBetweenCycles
  if O.iter == ii:
   moveBottom()
  if O.iter == ii+1:
   stopBottom()
 if O.iter == nGravityDeposition+nCycles*nStepsBetweenCycles:
  O.pause()
def moveBottom():
 v = dspl / O.dt
 for b in movables:
  b.state.vel = (0,0,-v)
def stopBottom():
 for b in movables:
  b.state.vel = (0,0,0)

# box
p000 = Vector3(0,0,0)
p100 = Vector3(widthl,0,0)
p200 = Vector3(widthl+widthc,0,0)
p300 = Vector3(widthl+widthc+widthr,0,0)
pxs = (p000,p100,p200,p300)
p001,p101,p201,p301 = [p+Vector3(0,0,height) for p in pxs]
p010,p110,p210,p310 = [p+Vector3(0,thick,0) for p in pxs]
p011,p111,p211,p311 = [p+Vector3(0,thick,height) for p in pxs]
p00b,p10b,p20b,p30b = [p+Vector3(0,0,deep) for p in pxs]
p01b,p11b,p21b,p31b = [p+Vector3(0,thick,deep) for p in pxs]
def rect(vs,**kw):
 v1,v2,v3,v4 = vs
 return [
  facet((v1,v2,v3),**kw),
  facet((v1,v3,v4),**kw),
 ]
movables = rect((p100,p200,p210,p110)) # bottom center
rects = (
 (p000,p100,p110,p010), # bottom left
 (p200,p300,p310,p210), # bottom left
 (p000,p010,p011,p001), # left
 (p300,p310,p311,p301), # right
 (p000,p100,p101,p001), # front left
 (p100,p200,p201,p101), # front center
 (p200,p300,p301,p201), # front right
 (p010,p110,p111,p011), # back left
 (p110,p210,p211,p111), # back center
 (p210,p310,p311,p211), # back right
 (p100,p200,p20b,p10b), # front center below
 (p110,p210,p21b,p11b), # back center below
 (p100,p110,p11b,p10b), # left below
 (p200,p210,p21b,p20b), # right below
)
rects = movables + sum((rect(r) for r in rects),[])
O.bodies.append(rects)

# gravel
polyhedra_utils.fillBox((0,0,0),(width,thick,fillBoxHFactor*height),mat,sizemin=3*[sizeMin],sizemax=3*[sizeMax],seed=1)
###########
#to record the values

#O.forces.f(id)
#sum((O.forces.f(id) for id in ids), Vector3.Zero)
#
#f = Vector3.Zero
#for i in ids:
# f += O.forces.f(i)

middle = 0,1
lef= 2,3
right = 4,5

#trapdoor.state.displ()
O.bodies[1].state.displ()

##
for b in O.bodies:
   if not isinstance(b.shape,Polyhedra): continue
   pos = b.state.pos # center of polyhedron
   ori = b.state.ori
   vs = b.shape.v
   vs = [pos + ori*v for v in vs] # vertices of polyhedron in global coordinate system
##

Hello,
what is written in the script directly (not in function etc) is executed and never called again (like O.bodies[1].state.displ()). To be able to plot the data, first you need to coleect the data during simulation:

#############
from yade import plot # put such import at the beginning of the file

...

O.engines = [
   ...
   PyRunner(iterPeriod=10,command='plotAddData()'),
]

...

def computeSurfaceDisplacement():
   vertices = []
   for b in O.bodies:
      if not isinstance(b.shape,Polyhedra): continue
      pos = b.state.pos # center of polyhedron
      ori = b.state.ori
      vs = b.shape.v
      vs = [pos + ori*v for v in vs] # vertices of polyhedron in global coordinate system
      vertices.extend(vs)
   maxZ = max(v[2] for v in vertices) # you would need something much more sophisticated, just to give you an idea
   return maxZ
##

def plotAddData():
   sd = computeSurfaceDisplacement()
   plot.addData(
      i = O.iter,
      tDspl = O.bodies[1].state.displ(),
      fLeft = sum(O.forces.f(id)[2] for id in (2,3)),
      surfDspl = sd,
   )

plot.plots = {'i':'surfDspl', 'i ':'fLeft', 'i ':'tDspl'}
##############

for proper surfaceDisplacement, you would probably need some record of "initial" surface after gravity deposition. You can do it in checker function:
#############
def checker():
   if O.iter == nGravityDeposition:
      recordInitialSurface()
   ...

def recordInitialSurface():
   global initSurf
   ... # some code here to "save" what the surface is
#############

cheers
Jan

hi,

I hope you are well, you helped me with the code below : now I got a little problem, when I run it the middle ( the trapdoor does not drop) can you please help me modify it, such the the trapdoor moves at a rate of 1.0mm per min for 100mm( thus after deposition)

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

rom yade import pack
from math import *
from yade import plot
from yade import qt

gravity = 100

#box dimensions
widthl = .3
widthr = .3
widthc = .3
height = .3
thick = .1
deep = -.2

#size of grains
sizeMin = 40e-3
sizeMax = 60e-3

frictionAngle = .5
young = 1e8 #stiffness
poisson = 10

dt = 1e-3 #time step
nGravityDeposition = 250 #how long to run initial gravity deposition
nCycles = 3 #how many jumps to run afterwards
nStepsBetweenCycles = 200 #number of times steps between jumps
dspl = 20e-3

#how much larger the initial make cloud box should be
fillBoxHFactor = 3

from yade import polyhedra_utils

width = widthl + widthc + widthr

#mat
mat = PolyhedraMat(young=young,poisson=10,frictionAngle=frictionAngle)
O.materials.append(mat)

def computeSurfaceDisplacement():
 vertices = []
 for b in O.bodies:
   if not isinstance(b.shape,Polyhedra): continue
   pos = b.state.pos
   ori = b.state.ori
   vs = b.shape.v
   vs = [pos + ori*v for v in vs]
   vertices.extend(vs)
   maxZ = max(v[2] for v in vertices)
   return maxZ

def plotAddData():
 sd = computeSurfaceDisplacement()
 plot.addData(
  i = O.iter,
  tDspl = O.bodies[1].state.displ().norm(),
  fLeft = sum(O.forces.f(id)[2] for id in (2,3)),
  fRight = sum(O.forces.f(id)[2] for id in (4,5)),
  surfDspl = sd,
  )
 #plot.plots = {'i ':'fRight'}
 #plot.plots = {'i ':'fRight'}
 plot.plots = {'i ':'surfDspl'}
 #plot.plots = {'i':'tDspl'}
 #plot.plots = {'i':'surfDspl','i':'fLeft','i':'tDspl'}
#################################

#engines
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
      [Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
      [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()],
      [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()],
   ),
   NewtonIntegrator(damping=0.4,gravity=(0,0,-gravity)),
 #commnet out the original PyRunner, replace with new one
 #PyRunner(iterPeriod=1,command='checker()'),
 PyRunner(iterPeriod = 10,command='plotAddData()'),
]

O.dt = dt

def moveBottom():
 v = dspl / O.dt
 for b in movables:
  b.state.vel = (0,0,-v)
def stopBottom():
 for b in movables:
  b.state.vel = (0,0,0)
def checker():
 for i in range(nCycles):
  ii = nGravityDeposition+i*nStepsBetweenCycles
  if O.iter == ii:
   moveBottom()
  if O.iter == ii+1:
   stopBottom()
 if O.iter == nGravityDeposition+nCycles*nStepsBetweenCycles:
  O.pause()

# box
p000 = Vector3(0,0,0)
p100 = Vector3(widthl,0,0)
p200 = Vector3(widthl+widthc,0,0)
p300 = Vector3(widthl+widthc+widthr,0,0)
pxs = (p000,p100,p200,p300)
p001,p101,p201,p301 = [p+Vector3(0,0,height) for p in pxs]
p010,p110,p210,p310 = [p+Vector3(0,thick,0) for p in pxs]
p011,p111,p211,p311 = [p+Vector3(0,thick,height) for p in pxs]
p00b,p10b,p20b,p30b = [p+Vector3(0,0,deep) for p in pxs]
p01b,p11b,p21b,p31b = [p+Vector3(0,thick,deep) for p in pxs]

def rect(vs,**kw):
 v1,v2,v3,v4 = vs
 return [
  facet((v1,v2,v3),**kw),
  facet((v1,v3,v4),**kw),
 ]

movables = rect((p100,p200,p210,p110)) # bottom center
rects = (
 (p000,p100,p110,p010), # bottom left
 (p200,p300,p310,p210), # bottom left
 (p000,p010,p011,p001), # left
 (p300,p310,p311,p301), # right
 (p000,p100,p101,p001), # front left
 (p100,p200,p201,p101), # front center
 (p200,p300,p301,p201), # front right
 (p010,p110,p111,p011), # back left
 (p110,p210,p211,p111), # back center
 (p210,p310,p311,p211), # back right
 (p100,p200,p20b,p10b), # front center below
 (p110,p210,p21b,p11b), # back center below
 (p100,p110,p11b,p10b), # left below
 (p200,p210,p21b,p20b), # right below
)
rects = movables + sum((rect(r) for r in rects),[])
O.bodies.append(rects)

# gravel
polyhedra_utils.fillBox((0,0,0),(width,thick,fillBoxHFactor*height),mat,sizemin=3*[sizeMin],sizemax=3*[sizeMax],seed=1)
###########

Hello,
the trapdoor motion is done by moveBottom function, which is called in checker() function, which is not used in your script. Just uncomment the line
 PyRunner(iterPeriod=1,command='checker()'),
and it works.
How much and how often is defined by 'nStepsBetweenCycles' and 'dspl' at the beginning of the script (it jumps by dspl every nStepsBetweenCycles time steps, i.e. each nStepsBetweenCycles*O.dt seconds. You can adjust those values to meet your requirements.
cheers
Jan