Contact forces of Polyhedra

Hello,

please read [1] and try to provide a MWE.

W = working.
I have tried your code with:
O.bodies[48].dynamic=False
IndexError: Body id out of range.

M=minimal.
In this case two (or three if you want to have more than one interaction) polyhedrons is enough.
Furthermore they can be much simpler, e.g. tetrahedrons or cubes or something like that should be enough to demonstrate the problem (which does not seem to be related to polyhedron shape).

> 1. How to Compute the contact forces between 2 polyhedras [i choose].

What does "compute" mean?
- the theory what equations are used for force computation?
- how to get the value already computed by Yade?

> 2. How to display and save these outcome results

What does "display" mean? print the values? graphical output?
How you want to use the saved value?
What does "these outcome results" mean? the values of contact forces?

cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.InteractionContainer

PS: a few notes to the script (and its relation to MWE)

> from __future__ import print_function
> from builtins import range
> from numpy import arange
> import numpy as np
> import random

print, range, arange, np and random are not used in the script

> for i in O.interactions:
> if not i.isReal: continue

Interactions in "for i in O.interactions" are always real [2]

> waitIfBatch()

also probably not needed in the MWE

Thank you for your quick answer:
[I have minimize my script and forgot to correct the id number of O.bodies[48].dynamic=False. (19 instead of 48)

 1. In one it means both the equations used for computation and how to get this values
 2. yes the outcome results means the contact forces and how to print them.

My corrected script is:
*************************************************************************************
from __future__ import print_function
from builtins import range
from yade import qt
from numpy import arange
import numpy as np
import random
from yade import polyhedra_utils
from yade import plot

cg1=0.23 #center of gravity of block range 1 level 1 (m)
cg2 =0.69 #center of gravity of block range 2 level 2 (m)

m = PolyhedraMat() #for the straw material
m.density = 120 #kg/m3
m.young = 3e7 #Pa (N/m2) 3e7
m.frictionAngle = 0.78
O.materials.append(m)

n = PolyhedraMat() #for the floor/roof material made from wood
n.density = 740 #kg/m3
n.young = 8.5e8 #Pa 8.5e8
n.frictionAngle = 0.65
O.materials.append(n)

# dimensions of the rock blocks ****************************************************
width = .36 # m
length = .72 # m
height = .46 # m

# COMPOSITION OF RANGE 1 level 1 ***************************************************

#vertices1 = [[0,0,0],[72,0,0],[0,36,0],[72,36,0],[0,0,46],[72,0,46],[0,36,46],[72,36,46]]
vertices1 = [[0,0,0],[length,0,0],[0,width,0],[length,width,0],[0,0,height],[length,0,height],[0,width,height],[length,width,height]]
b1 = polyhedra_utils.polyhedra(m,v=vertices1)
b1.state.pos = (length*0.5,width*0.5,cg1)
O.bodies.append(b1)

#vertices2 = [[72,0,0],[144,0,0],[72,36,0],[144,36,0],[72,0,46],[144,0,46],[72,36,46],[144,36,46]]
vertices2 = [[length,0,0],[length*2,0,0],[length,width,0],[length*2,width,0],[length,0,height],[length*2,0,height],[length,width,height],[length*2,width,height]]
b2 = polyhedra_utils.polyhedra(m,v=vertices2)
b2.state.pos = (length*1.5,width*0.5,cg1)
O.bodies.append(b2)

#vertices3 = [[144,0,0],[216,0,0],[144,36,0],[216,36,0],[144,0,46],[216,0,46],[144,36,46],[216,36,46]]
vertices3 = [[length*2,0,0],[length*3,0,0],[length*2,width,0],[length*3,width,0],[length*2,0,height],[length*3,0,height],[length*2,width,height],[length*3,width,height]]
b3 = polyhedra_utils.polyhedra(m,v=vertices3)
b3.state.pos = (length*2.5,width*0.5,cg1)
O.bodies.append(b3)

#vertices4 = [[180,36,0],[216,36,0],[180,108,0],[216,108,0],[180,36,46],[216,36,46],[180,108,46],[216,108,46]]
vertices4 = [[length*2.5,width,0],[length*3,width,0],[length*2.5,width*3,0],[length*3,width*3,0],[length*2.5,width,height],[length*3,width,height],[length*2.5,width*3,height],[length*3,width*3,height]]
b4 = polyhedra_utils.polyhedra(m,v=vertices4)
b4.state.pos = (length*2.75,width*2,cg1)
O.bodies.append(b4)

#vertices5 = [[180,108,0],[216,108,0],[180,180,0],[180,180,0],[180,108,46],[216,108,46],[180,180,46],[216,180,46]]
vertices5 = [[length*2.5,width*3,0],[length*3,width*3,0],[length*2.5,width*5,0],[length*3,width*5,0],[length*2.5,width*3,height],[length*3,width*3,height],[length*2.5,width*5,height],[length*3,width*5,height]]
b5 = polyhedra_utils.polyhedra(m,v=vertices5)
b5.state.pos = (length*2.75,width*4,cg1)
O.bodies.append(b5)

#vertices6 = [[144,180,0],[216,180,0],[144,216,0],[216,216,0],[144,180,46],[216,180,46],[144,216,46],[216,216,46]]
vertices6 = [[length*2,width*5,0],[length*3,width*5,0],[length*2,width*6,0],[length*3,width*6,0],[length*2,width*5,height],[length*3,width*5,height],[length*2,width*6,height],[length*3,width*6,height]]
b6 = polyhedra_utils.polyhedra(m,v=vertices6)
b6.state.pos = (length*2.5,width*5.5,cg1)
O.bodies.append(b6)

#vertices7 = [[72,180,0],[144,180,0],[72,216,0],[144,216,0],[72,180,46],[144,180,46],[72,216,46],[144,216,46]]
vertices7 = [[length,width*5,0],[length*2,width*5,0],[length,width*6,0],[length*2,width*6,0],[length,width*5,height],[length*2,width*5,height],[length,width*6,height],[length*2,width*6,height]]
b7 = polyhedra_utils.polyhedra(m,v=vertices7)
b7.state.pos = (length*1.5,width*5.5,cg1)
O.bodies.append(b7)

#vertices8 = [[0,180,0],[72,180,0],[0,216,0],[72,216,0],[0,180,46],[72,180,46],[0,216,46],[72,216,46]]
vertices8 = [[0,width*5,0],[length,width*5,0],[0,width*6,0],[length,width*6,0],[0,width*5,height],[length,width*5,height],[0,width*6,height],[length,width*6,height]]
b8 = polyhedra_utils.polyhedra(m,v=vertices8)
b8.state.pos = (length*.5,width*5.5,cg1)
O.bodies.append(b8)

#vertices9 = [[0,144,0],[36,144,0],[0,180,0],[36,180,0],[0,144,46],[36,144,46],[0,180,46],[36,180,46]]
vertices9 = [[0,width*4,0],[length*.5,width*4,0],[0,width*5,0],[length*.5,width*5,0],[0,width*4,height],[length*.5,width*4,height],[0,width*5,height],[length*.5,width*5,height]]
b9 = polyhedra_utils.polyhedra(m,v=vertices9)
b9.state.pos = (length*0.25,width*4.5,cg1)
O.bodies.append(b9)

#vertices10 = [[0,36,0],[36,36,0],[0,72,0],[36,72,0],[0,36,46],[36,36,46],[0,72,46],[36,72,46]]
vertices10 = [[0,width,0],[length*.5,width,0],[0,width*2,0],[length*.5,width*2,0],[0,width,height],[length*.5,width,height],[0,width*2,height],[length*.5,width*2,height]]
b10 = polyhedra_utils.polyhedra(m,v=vertices10)
b10.state.pos = (length*0.25,width*1.5,cg1)
O.bodies.append(b10)

# COMPOSITION OF RANGE 2 level 2 *******************************************************
#vertices11 = [[0,0,46],[36,0,46],[0,72,46],[36,72,46],[0,0,92],[36,0,92],[0,72,92],[36,72,92]]
vertices11 = [[0,0,height],[length*.5,0,height],[0,width*2,height],[length*.5,width*2,height],[0,0,height*2],[length*.5,0,height*2],[0,width*2,height*2],[length*.5,width*2,height*2]]
b11 = polyhedra_utils.polyhedra(m,v=vertices11)
b11.state.pos = (length*.25,width,cg2)
O.bodies.append(b11)

#vertices12 = [[36,0,46],[108,0,46],[36,36,46],[108,36,46],[36,0,92],[108,0,92],[36,36,92],[108,36,92]]
vertices12 = [[length*.5,0,height],[length*1.5,0,height],[length*.5,width,height],[length*1.5,width,height],[length*.5,0,height*2],[length*1.5,0,height*2],[length*.5,width,height*2],[length*1.5,width,height*2]]
b12 = polyhedra_utils.polyhedra(m,v=vertices12)
b12.state.pos = (length*1,width*.5,cg2)
O.bodies.append(b12)

#vertices13 = [[108,0,46],[180,0,46],[108,36,46],[180,36,46],[108,0,92],[180,0,92],[108,36,92],[180,36,92]]
vertices13 = [[length*1.5,0,height],[length*2.5,0,height],[length*1.5,width,height],[length*2.5,width,height],[length*1.5,0,height*2],[length*2.5,0,height*2],[length*1.5,width,height*2],[length*2.5,width,height*2]]
b13 = polyhedra_utils.polyhedra(m,v=vertices13)
b13.state.pos = (length*2,width*.5,cg2)
O.bodies.append(b13)

#vertices14 = [[180,0,46],[216,0,46],[180,72,46],[216,72,46],[180,0,92],[216,0,92],[180,72,92],[216,72,92]]
vertices14 = [[length*2.5,0,height],[length*3,0,height],[length*2.5,width*2,height],[length*3,width*2,height],[length*2.5,0,height*2],[length*3,0,height*2],[length*2.5,width*2,height*2],[length*3,width*2,height*2]]
b14 = polyhedra_utils.polyhedra(m,v=vertices14)
b14.state.pos = (length*2.75,width*1,cg2)
O.bodies.append(b14)

#vertices15 = [[180,72,46],[216,72,46],[180,144,46],[216,144,46],[180,72,92],[216,72,92],[180,144,92],[216,144,92]]
vertices15 = [[length*2.5,width*2,height],[length*3,width*2,height],[length*2.5,width*4,height],[length*3,width*4,height],[length*2.5,width*2,height*2],[length*3,width*2,height*2],[length*2.5,width*4,height*2],[length*3,width*4,height*2]]
b15 = polyhedra_utils.polyhedra(m,v=vertices15)
b15.state.pos = (length*2.75,width*3,cg2)
O.bodies.append(b15)

#vertices16 = [[180,144,46],[216,144,46],[180,216,46],[216,216,46],[180,144,92],[216,144,92],[180,216,92],[216,216,92]]
vertices16 = [[length*2.5,width*4,height],[length*3,width*4,height],[length*2.5,width*6,height],[length*3,width*6,height],[length*2.5,width*4,height*2],[length*3,width*4,height*2],[length*2.5,width*6,height*2],[length*3,width*6,height*2]]
b16 = polyhedra_utils.polyhedra(m,v=vertices16)
b16.state.pos = (length*2.75,width*5,cg2)
O.bodies.append(b16)

#vertices17 = [[108,180,46],[180,180,46],[108,216,46],[180,216,46],[108,180,92],[180,180,92],[108,216,92],[180,216,92]]
vertices17 = [[length*1.5,0,height],[length*2.5,0,height],[length*1.5,width,height],[length*2.5,width,height],[length*1.5,0,height*2],[length*2.5,0,height*2],[length*1.5,width,height*2],[length*2.5,width,height*2]]
b17 = polyhedra_utils.polyhedra(m,v=vertices17)
b17.state.pos = (length*2,width*5.5,cg2)
O.bodies.append(b17)

#vertices18 = [[36,180,46],[108,180,46],[36,216,46],[108,216,46],[36,180,92],[108,180,92],[36,216,92],[108,216,92]]
vertices18 = [[length*0.5,width*5,height],[length*1.5,width*5,height],[length*0.5,width*6,height],[length*1.5,width*6,height],[length*0.5,width*5,height*2],[length*1.5,width*5,height*2],[length*0.5,width*6,height*2],[length*1.5,width*6,height*2]]
b18 = polyhedra_utils.polyhedra(m,v=vertices18)
b18.state.pos = (length,width*5.5,cg2)
O.bodies.append(b18)

#vertices19 = [[0,144,46],[36,144,46],[0,216,46],[36,216,46],[0,144,92],[36,144,92],[0,216,92],[36,216,92]]
vertices19 = [[0,width*4,height],[length*.5,width*4,height],[0,width*6,height],[length*.5,width*6,height],[0,width*4,height*2],[length*.5,width*4,height*2],[0,width*6,height*2],[length*.5,width*6,height*2]]
b19 = polyhedra_utils.polyhedra(m,v=vertices19)
b19.state.pos = (length*.25,width*5,cg2)
O.bodies.append(b19)

# floor **********************************************************
# floor dimensions ************************************************
lengthf= 4 # m
widthf =4 # m
heightf=-0.1 # m
#vertices_floor = [[0,0,0],[400,0,0],[0,400,0],[400,400,0],[0,0,10],[400,0,10],[0,400,10],[400,400,10]]

vertices_floor = [[0,0,0],[lengthf,0,0],[0,widthf,0],[lengthf,widthf,0],[0,0,heightf],[lengthf,0,heightf],[0,widthf,heightf],[lengthf,widthf,heightf]]
floor = polyhedra_utils.polyhedra(n,v=vertices_floor, color=(.6,.4,.3))
floor.state.pos = (1.08,1.08,heightf*0.5)
O.bodies.append(floor)
O.bodies[19].dynamic=False

qt.Controller()
V = qt.View()

# ENGINE ********************************************************************************

def checkUnbalanced():
   print ('iter %d, time elapsed %f, unbalanced forces = %f'%(O.iter, O.realtime, utils.unbalancedForce()))
   #if unbalancedForce()<.05: # to ensure stability
    # O.pause()
     # plot.saveDataTxt('strawbale-static-state.txt.bz2')

for i in O.interactions:
  if not i.isReal: continue
  cp = i.geom.contactPoint
  fn = i.phys.normalForce
  fs = i.phys.shearForce

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb()]),
   InteractionLoop(
      [Ig2_Wall_Polyhedra_PolyhedraGeom(), Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
      [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()], # collision "physics" geometry
      [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()] # contact law -- apply forces
   ),
   NewtonIntegrator(damping=0.2,gravity=(0,0,-9.81)), #apply gravity force to particles. Damping, numerical dissipation of energy
#call the check unbalanced function every 2 sec, the label cretaes an automatic variable referring to this engine
   PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
   ]

O.dt=0.00000164
O.run()

# when running with yade-batch, the script must not finish until the simulation is done fully
# this command will wait for that (has no influence in the non-batch mode)
waitIfBatch()
# end

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