# yade-users team mailing list archive

## Re: [Question #695154]: Polyhedra contact instability (potential Bug?)

```Question #695154 on Yade changed:

Description changed to:

During the simulation of polyhedra deposition in the box, I encountered instability problems (even for small increments, some grains behave like popcorn). I have made some animations to observe this and came up with the conclusion that it may be related to switching the contact conditions.  I have read the paper of Eliáš (which is full of clever ideas by the way), and I think that the issue may be the calculation of the direction of the normal force.

I conducted a simulation (see the script at the bottom of the post) that
seems to confirm that hypothesis. I have two polyhedra, a cube and a
prism (half of the cube) that are shifted. I slide one polyhedron to
make a contact starting from one side of the cube. As I move the
polyhedron, the intersection becomes symmetrical, and later it is
"closer" to the other (perpendicular) wall. I would expect the forces to
change gradually, but there is a sudden switch between vertical and
horizontal forces.

Everything below is speculation since I am not a C++ developer.

I have found the part of the code responsible for normal force
calculation (function FindNormal in [2]). It utilizes a function, that
has a known flaw [3], so it is a possible suspect for me.

I wish, I could be able to verify this by myself, but it is not my
league (yet ;) ). I am a more Python person.

[1] Eliáš, J. (2014). Simulation of railway ballast using crushable polyhedral particles. Powder Technology, 264, 458-465.
[3] https://github.com/CGAL/cgal/issues/2676

---------------------------------------------
# Yade version 2020.01a on Ubuntu 20.04

######################################## CONSTANTS
d = 0.1

######################################## FUNCTIONS

f0 = O.forces.f(0)[0]
f1 = O.forces.f(0)[1]
f2 = O.forces.f(0)[2]
t0 = O.forces.t(0)[0]
t1 = O.forces.t(0)[1]
t2 = O.forces.t(0)[2]

########################## MATERIALS

polyMat = PolyhedraMat()

########################## BODIES
p1 = polyhedra_utils.polyhedra(polyMat, v= [(0,0,0),(0,d,0),(d,d,0),(d,0,0),(0,0,d),(0,d,d),(d,d,d),(d,0,d)], fixed = True)
p1.state.pos = (0,0,0)

p2 = polyhedra_utils.polyhedra(polyMat, v=
[(0,0,0),(0,d,0),(d/2,d,0),(d/2,0,0),(0,0,d),(0,d,d),(d/2,d,d),(d/2,0,d)],
fixed = True)

p2.state.pos = (0,0.9*d,d)

O.bodies.append((p1,p2))

########################## ENGINES

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Polyhedra_Aabb()]),
InteractionLoop(
[Ig2_Polyhedra_Polyhedra_PolyhedraGeom()],
[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()],
[Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()]
),
NewtonIntegrator(),
]

O.dt= 0.001*polyhedra_utils.PWaveTimeStep()

plot.plots={'t':(('f0','r'),('f1','g'),('f2','b'),('t0','r--'),('t1','g--'),('t2','b--')),}
plot.plot(subPlots =False)

try:
qt.Controller()
v = qt.View()
v.ortho = True
v.viewDir = (-1,0,0)
v.showEntireScene()
except:
pass

p2.state.vel = (0,0,-1)
O.run(20000)

--