# Contact forces of Polyhedra

i am new user in yade. i am using
Using python version: 3.5.2 (default, Apr 16 2020, 17:47:17)
[GCC 5.4.0 20160609]
Ubunto 16.04

*******************************************************************************************************
i have 2 problems:

1. How to Compute the contact forces between 2 polyhedras [i choose].
2. How to display and save these outcome results

hiba
******************************************************************************
#Script

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

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[48].dynamic=False

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

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

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
]

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 *************************************************************************************

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
 Revision history for this message Jan Stránský (honzik) said on 2020-07-09: #1

Hello,

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?

> 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

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

 Revision history for this message khalifeh (hibakh) said on 2020-07-09: #2

[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 numpy import arange
import numpy as np
import random

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

 Revision history for this message Jan Stránský (honzik) said on 2020-07-09: #3

> I have minimize my script

> 1. In one it means both the equations used for computation and how to get this values

Theory: [3], Elias2013, Elias2014
Get: see below

> yes the outcome results means the contact forces

fn = i.phys.normalForce
fs = i.phys.shearForce

"i" could be "for i in O.interactions" or "i=O.interactions[id1,id2]"

> and how to print them.

use standard Python print:
print(fn,fs)

cheers
Jan

 Revision history for this message khalifeh (hibakh) said on 2020-07-10: #4

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