# polyhedrals after collide to facet, throwing away strongly

Hi everyone
I run this code for some polyhedrals, they fall down with gravity, but when polyhedrals collide to facet in bottom, all of them move up with high speed and spread in space. While they should stay on the floor for a little while.

note: Most of polyhedrals have sharp angle, does this have effect on this problem?

I would appreciate if you can help me.

ubuntu : 18.04

too many input data.
my codes:

import os
import math
import pylab
import matplotlib; matplotlib.rc('axes',grid=True)
from matplotlib import pyplot
import numpy as np
from numpy import *
from yade import export as expt

Dolomite = PolyhedraMat()
Dolomite.density = 2870 #kg/m^3
Dolomite.young = 1e7 #Pa
Dolomite.poisson = 0.2

Bound = PolyhedraMat()
Bound.density = 2870 #kg/m^3
Bound.young = 1e300 #Pa
Bound.poisson = 0.2

p1=((3.20999113285824,4.56682535905544,0.25),(3.20999113285824,4.92812419373700,0.25),(3.34925246942789,5.13990997478261,0.25),(3.48096338587174,5.13990997478261,0.25),(3.48096338587174,4.97471882295321,0.25),(3.38909183943153,4.70539658408990,0.25),(3.20999113285824,4.56682535905544,0.25),(3.20999113285824,4.56682535905544,-0.25),(3.20999113285824,4.92812419373700,-0.25),(3.34925246942789,5.13990997478261,-0.25),(3.48096338587174,5.13990997478261,-0.25),(3.48096338587174,4.97471882295321,-0.25),(3.38909183943153,4.70539658408990,-0.25),(3.20999113285824,4.56682535905544,-0.25))
p2=((3.60265201447845,7.35103416471193,0.25),(3.61010182622721,7.35103416471193,0.25),(3.71828876180858,7.48918606800436,0.25),(3.74275101349167,7.67253715703928,0.25),(3.74275101349167,7.82185097460873,0.25),(3.67394962688404,7.75228435593977,0.25),(3.60265201447845,7.35103416471193,0.25),(3.60265201447845,7.35103416471193,-0.25),(3.61010182622721,7.35103416471193,-0.25),(3.71828876180858,7.48918606800436,-0.25),(3.74275101349167,7.67253715703928,-0.25),(3.74275101349167,7.82185097460873,-0.25),(3.67394962688404,7.75228435593977,-0.25),(3.60265201447845,7.35103416471193,-0.25))

O.bodies.append(polyhedra_utils.polyhedra(Dolomite,v=p1,fixed=False, color=(0.1,0.5,0.2)))
O.bodies.append(polyhedra_utils.polyhedra(Dolomite,v=p2,fixed=False, color=(0.1,0.5,0.2)))

v1=((-2.4993,-0.01,2),(-2,-0.01,2),(-2.4993,-0.01,-2))
v2=((-2,-0.01,2),(-2,-0.01,-2),(-2.4993,-0.01,-2))
v3=((-1.99,-0.01,2),(-1.99,-0.01,-2),(-1,-0.01,2))
v4=((-1,-0.01,2),(-1,-0.01,-2),(-1.99,-0.01,-2))
v5=((-0.99,-0.01,2),(-0.99,-0.01,-2),(1,-0.01,2))
v6=((1,-0.01,2),(1,-0.01,-2),(-0.99,-0.01,-2))
v7=((1.01,-0.01,2),(1.01,-0.01,-2),(2,-0.01,2))
v8=((2,-0.01,2),(2,-0.01,-2),(1.01,-0.01,-2))
v9=((2.01,-0.01,2),(2.01,-0.01,-2),(3,-0.01,2))
v10=((3,-0.01,2),(3,-0.01,-2),(2.01,-0.01,-2))
v11=((3.01,-0.01,2),(3.01,-0.01,-2),(4.58,-0.01,2))
v12=((4.58,-0.01,2),(4.58,-0.01,-2),(3.01,-0.01,-2))
v13=((4.6162,0.0059,2),(5.25,0.93,2),(4.6162,0.0059,-2))
v14=((5.25,0.93,2),(5.25,0.93,-2),(4.6162,0.0059,-2))
v15=((5.25,0.94,2),(5.25,0.94,-2),(6.42,2.57,2))
v16=((6.42,2.57,2),(6.42,2.57,-2),(5.25,0.94,-2))
v17=((6.42,2.58,2),(6.42,2.58,-2),(6.7906,3.01,2))
v18=((6.7906,3.01,2),(6.7906,3.01,-2),(6.42,2.58,-2))

O.bodies.append(utils.facet(v1,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v2,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v3,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v4,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v5,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v6,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v7,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v8,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v9,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v10,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v11,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v12,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v13,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v14,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v15,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v16,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v17,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v18,color=(1,0,1), material=Bound))

for b in O.bodies:
if b.material is Dolomite:
b.state.blockedDOFs='zXY'

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Facet_Aabb()],verletDist=0.05),
InteractionLoop(
[Ig2_Polyhedra_Polyhedra_PolyhedraGeom(),Ig2_Facet_Polyhedra_PolyhedraGeom()],
[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()],
[Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()]
),
NewtonIntegrator(gravity=(0,-9.81,0),damping=0.2),
]

utils.calm()

O.dt=10e-6
O.saveTmp()
O.run()

## Question information

Language:
English Edit question
Status:
Open
For:
Assignee:
No assignee Edit question
Last query:
2020-03-21
2020-03-19
 Jan Stránský (honzik) said on 2020-02-26: #1

Hi,
this is a typical symptom of too large time step (especially if decreasing stiffness helps), try to decrease it.
cheers
Jan

 Mahdeyeh (mahdiye.sky) said on 2020-02-26: #2

Thank you a lot for all your helps Jan Stransky

>> this is a typical symptom of too large time step
Do you think if I decrease time step, maybe this problem is solved?
for more particles, I set time step:10e-8.

>> note: Most of polyhedrals have sharp angle, does this have effect on this problem?

 Mahdeyeh (mahdiye.sky) said on 2020-02-29: #3

Actually I think contact model in polyhedrals has problem!As the two polyhedrals collide very slowly, their velocity and kinetic energy increase abnormally!!

>> this is a typical symptom of too large time step (especially if decreasing stiffness helps)
Unfortunately decreasing stiffness or time step did help at all. :(

Your help in this would be appreciated

 Jan Stránský (honzik) said on 2020-03-02: #4

> note: Most of polyhedrals have sharp angle, does this have effect on this problem?

no

> decreasing young module until 10e3 can help the situation
> Unfortunately decreasing stiffness or time step did help at all. :(

please be consistent. So does it help or not?

> for more particles, I set time step:10e-8.

Why? the value of time step should not depend on number of particles..

> Actually I think contact model in polyhedrals has problem!As the two polyhedrals collide very slowly, their velocity and kinetic energy increase abnormally!!

can you please provide a Minimal Working Example (e.g. two polhedrons) reproducing this? thanks

cheers
Jan

 Mahdeyeh (mahdiye.sky) said on 2020-03-02: #5

> please be consistent. So does it help or not?

("decreasing young module until 10e3 can help the situation"): this was a mistake and decreasing stiffness or time step did help at all.
If stiffness reaches less than 10e5, then particles are passed through facet.

> Why? the value of time step should not depend on number of particles..
Appropriate time step for simulation with all particles is 10e-8 to avoid exploding particles, in this test (for saving time) 10e-5 was selected for time step because particles number are small and they do not explode, whereas I checked, changing time step does not change situation (does not solve this problem).

 Jan Stránský (honzik) said on 2020-03-02: #7

Hello,

I have tried the script provided in OP, resulting in:
NameError: name 'ListVer2' is not defined
could you please provide a test script which we can try? thanks

cheers
Jan

 Mahdeyeh (mahdiye.sky) said on 2020-03-02: #8

> can you please provide a Minimal Working Example (e.g. two polhedrons) reproducing this?

Here is small data and input for my code, please copy in file:

*LWPOLYLINE
4.54097400033104,6.49184245861868
4.65891804001884,6.58769165827740
4.51425920185592,6.93490024048091
4.46361293495481,6.82361276486310
4.54097400033104,6.49184245861868
*LWPOLYLINE
3.38239216165079,4.01327469219457
3.43869334720801,4.27025377921128
3.31240868369778,4.27786275211975
3.18944272805255,4.21917403263651
3.17292105915418,4.00731315957292
3.38239216165079,4.01327469219457
*LWPOLYLINE
3.07472367219926,4.22146522059046
3.00951705252640,4.00710000000000
3.15872105915418,4.00710000000000
3.16994885128804,4.22031962661348
3.07472367219926,4.22146522059046
*LWPOLYLINE
3.19496730804472,4.54686225863838
3.10768290255748,4.56486506811422
3.08101342018444,4.41477399856679
3.07472367219926,4.23195436841659
3.17133336059188,4.25583338934358
3.19496730804472,4.54686225863838
*LWPOLYLINE
3.71267645811818,4.17335865185637
3.80626391854923,4.37990041469040
3.85709007230427,4.37990041469040
3.81834238444732,4.15093746143616
3.71267645811818,4.17335865185637
*LWPOLYLINE
3.86833315936705,7.18115411676252
3.75805619368548,6.85052980253022
3.96022603149273,6.83823476145616
3.98624257576362,7.12695113182557
3.86833315936705,7.18115411676252
*LWPOLYLINE
3.75705101349167,7.68673972685043
3.75705101349167,8.04757108283743
3.89622020287160,8.25921672832907
4.02762326650518,8.25921672832907
4.02762326650518,8.09425875012654
3.93578182400268,7.82502476121894
3.75705101349167,7.68673972685043
*LWPOLYLINE
4.29367658088641,7.54996433412856
4.25127773179602,7.75135995952029
4.17124793309971,7.75135995952029
4.11338686757211,7.59716107043463
4.21414905751508,7.54996433412856
4.29367658088641,7.54996433412856
*LWPOLYLINE
3.60265201447845,7.35103416471193
3.61010182622721,7.35103416471193
3.71828876180858,7.48918606800436
3.74275101349167,7.67253715703928
3.74275101349167,7.82185097460873
3.67394962688404,7.75228435593977
3.60265201447845,7.35103416471193
*LWPOLYLINE
3.75665696623051,7.66810154888463
3.94362927435906,7.81276307636047
4.02303237072992,7.79446112714408
3.89479310068892,7.42404447036142
3.73297934687869,7.49095982637576
3.75665696623051,7.66810154888463

 Jan Stránský (honzik) said on 2020-03-02: #9

Thanks for the working script. I have tried it. With O.dt=10e-5, I see clearly what was described.
Again, it is a typical symptom of too large time step (making the simulation unstable) and I don't believe it is not overcome with decreasing the time step to a **sufficient** value

> Unfortunately decreasing ... time step did help at all. :(

what values have you tried?

> Actually I think contact model in polyhedrals has problem!As the two polyhedrals collide very slowly, their velocity and kinetic energy increase abnormally!!

Can you please provide a script (two bodies as is described) reproducing this behavior?
(even the number of bodies is decreased, still the simulation is far from minimal to show and study such phenomenons in detail)

cheers
Jan

 Mahdeyeh (mahdiye.sky) said on 2020-03-03: #10

> what values have you tried?
O.dt=10e-8

> Can you please provide a script (two bodies as is described) reproducing this behavior?

import os
import math
import pylab
import matplotlib; matplotlib.rc('axes',grid=True)
from matplotlib import pyplot
import numpy as np
from numpy import *
from yade import export as expt

Dolomite = PolyhedraMat()
Dolomite.density = 2870 #kg/m^3
Dolomite.young = 24.36e9 #Pa
Dolomite.poisson = 0.2

Bound = PolyhedraMat()
Bound.density = 2870 #kg/m^3
Bound.young = 60e10 #Pa
Bound.poisson = 0.2

p1=((4.54097400033104,6.49184245861868,0.25),(4.65891804001884,6.58769165827740,0.25),(4.51425920185592,6.93490024048091,0.25),(4.46361293495481,6.82361276486310,0.25),(4.54097400033104,6.49184245861868,0.25),(4.54097400033104,6.49184245861868,-0.25),(4.65891804001884,6.58769165827740,-0.25),(4.51425920185592,6.93490024048091,-0.25),(4.46361293495481,6.82361276486310,-0.25),(4.54097400033104,6.49184245861868,-0.25))
p2=((3.60265201447845,7.35103416471193,0.25),(3.61010182622721,7.35103416471193,0.25),(3.71828876180858,7.48918606800436,0.25),(3.74275101349167,7.67253715703928,0.25),(3.74275101349167,7.82185097460873,0.25),(3.67394962688404,7.75228435593977,0.25),(3.60265201447845,7.35103416471193,0.25),(3.60265201447845,7.35103416471193,-0.25),(3.61010182622721,7.35103416471193,-0.25),(3.71828876180858,7.48918606800436,-0.25),(3.74275101349167,7.67253715703928,-0.25),(3.74275101349167,7.82185097460873,-0.25),(3.67394962688404,7.75228435593977,-0.25),(3.60265201447845,7.35103416471193,-0.25))

O.bodies.append(polyhedra_utils.polyhedra(Dolomite,v=p1,fixed=False, color=(0.1,0.5,0.2)))
O.bodies.append(polyhedra_utils.polyhedra(Dolomite,v=p2,fixed=False, color=(0.1,0.5,0.2)))

v1=((-2.4993,-0.01,2),(-2,-0.01,2),(-2.4993,-0.01,-2))
v2=((-2,-0.01,2),(-2,-0.01,-2),(-2.4993,-0.01,-2))
v3=((-1.99,-0.01,2),(-1.99,-0.01,-2),(-1,-0.01,2))
v4=((-1,-0.01,2),(-1,-0.01,-2),(-1.99,-0.01,-2))
v5=((-0.99,-0.01,2),(-0.99,-0.01,-2),(1,-0.01,2))
v6=((1,-0.01,2),(1,-0.01,-2),(-0.99,-0.01,-2))
v7=((1.01,-0.01,2),(1.01,-0.01,-2),(2,-0.01,2))
v8=((2,-0.01,2),(2,-0.01,-2),(1.01,-0.01,-2))
v9=((2.01,-0.01,2),(2.01,-0.01,-2),(3,-0.01,2))
v10=((3,-0.01,2),(3,-0.01,-2),(2.01,-0.01,-2))
v11=((3.01,-0.01,2),(3.01,-0.01,-2),(4.58,-0.01,2))
v12=((4.58,-0.01,2),(4.58,-0.01,-2),(3.01,-0.01,-2))
v13=((4.6162,0.0059,2),(5.25,0.93,2),(4.6162,0.0059,-2))
v14=((5.25,0.93,2),(5.25,0.93,-2),(4.6162,0.0059,-2))
v15=((5.25,0.94,2),(5.25,0.94,-2),(6.42,2.57,2))
v16=((6.42,2.57,2),(6.42,2.57,-2),(5.25,0.94,-2))
v17=((6.42,2.58,2),(6.42,2.58,-2),(6.7906,3.01,2))
v18=((6.7906,3.01,2),(6.7906,3.01,-2),(6.42,2.58,-2))
v19=((6.7906,3.0162,2),(6.7939,4,2),(6.7906,3.0162,-2))
v20=((6.7906,3.0162,-2),(6.7939,4,-2),(6.7939,4,2))
v21=((6.7953,4.0079,2),(7.02,5.6,2),(6.7953,4.0079,-2))
v22=((7.02,5.6,2),(7.02,5.6,-2),(6.7953,4.0079,-2))
v23=((7.02,5.61,2),(7.02,5.61,-2),(7.3,7.01,2))
v24=((7.3,7.01,2),(7.3,7.01,-2),(7.02,5.61,-2))
v25=((7.3,7.02,2),(7.3,7.02,-2),(7.54,8.51,2))
v26=((7.54,8.51,2),(7.54,8.51,-2),(7.3,7.02,-2))
v27=((7.54,8.52,2),(7.54,8.52,-2),(7.8,10,2))
v28=((7.8,10,2),(7.8,10,-2),(7.54,8.52,-2))
v29=((7.8,10.01,2),(7.8,10.01,-2),(8.14,12,2))
v30=((8.14,12,2),(8.14,12,-2),(7.8,10.01,-2))
v31=((8.14,12.01,2),(8.14,12.01,-2),(8.5,14,2))
v32=((8.5,14,2),(8.5,14,-2),(8.14,12,-2))
v33=((8.5,14.01,2),(8.5,14.01,-2),(8.75,15.5,2))
v34=((8.75,15.5,2),(8.75,15.5,-2),(8.5,14.01,-2))
v35=((8.75,15.51,2),(8.75,15.51,-2),(9.0615,16.9981,2))
v36=((9.0615,16.9981,2),(9.0615,16.9981,-2),(8.75,15.51,-2))
v37=((9.0621,17.0187,2),(9.09,19,2),(9.0621,17.087,-2))
v38=((9.09,19,2),(9.09,19,-2),(9.062,17.0187,-2))
v39=((9.09,19.01,2),(9.09,19.01,-2),(9.2006,20.9992,2))
v40=((9.2006,20.9992,2),(9.2006,20.9992,-2),(9.09,19.01,-2))
v41=((2.2009,21.0924,2),(2.2009,21.0924,-2),(2.22,19.01,2))
v42=((2.22,19.01,2),(2.22,19.01,-2),(2.2009,21.0924,-2))
v43=((2.22,19,2),(2.22,19,-2),(2.1538,17.043,2))
v44=((2.1538,17.043,2),(2.1538,17.043,-2),(2.22,19,-2))
v45=((2.0524,17.0104,2),(2.0524,17.0104,-2),(1.76,15.01,2))
v46=((1.76,15.01,2),(1.76,15.01,-2),(2.0524,17.0104,-2))
v47=((1.76,15,2),(1.76,15,-2),(1.42,13.01,2))
v48=((1.42,13.01,2),(1.42,13.01,-2),(1.76,15,-2))
v49=((1.42,13,2),(1.42,13,-2),(1.16,11.49,2))
v50=((1.16,11.49,2),(1.16,11.49,-2),(1.42,13,-2))
v51=((1.16,11.48,2),(1.16,11.48,-2),(0.8699,10.0876,2))
v52=((0.8699,10.0876,2),(0.8699,10.0876,-2),(1.16,11.48,-2))
v53=((1.1059,9.9916,2),(1.1059,9.9916,-2),(2.5,9.9927,2))
v54=((2.5,9.9927,2),(2.5,9.9927,-2),(1.1059,9.9916,-2))
v55=((2.51,9.9927,2),(2.51,9.9927,-2),(4.0026,9.9927,2))
v56=((2.51,9.9927,-2),(4.0026,9.9927,-2),(4.0026,9.9927,2))
v57=((4,9.94,2),(4,9.94,-2),(3.87,9,2))
v58=((3.87,9,2),(3.87,9,-2),(4,9.94,-2))
v59=((3.87,8.99,2),(3.87,8.99,-2),(3.69,8,2))
v60=((3.69,8,2),(3.69,8,-2),(3.87,8.99,-2))
v61=((3.69,7.99,2),(3.69,7.99,-2),(3.51,7,2))
v62=((3.51,7,2),(3.51,7,-2),(3.69,7.99,-2))
v63=((3.51,6.99,2),(3.51,6.99,-2),(3.34,5.99,2))
v64=((3.34,5.99,2),(3.34,5.99,-2),(3.51,6.99,-2))
v65=((3.34,5.98,2),(3.34,5.98,-2),(3.16,4.98,2))
v66=((3.16,4.98,2),(3.16,4.98,-2),(3.34,5.98,-2))
v67=((3.16,4.97,2),(3.16,4.97,-2),(2.95,4.0001,2))
v68=((2.95,4.0001,2),(2.95,4.0001,-2),(3.16,4.97,-2))
v69=((0,4,2),(0,4,-2),(2.9,3.9,2))
v70=((0,4,-2),(2.9,3.9,2),(2.9,3.9,-2))

O.bodies.append(utils.facet(v1,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v2,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v3,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v4,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v5,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v6,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v7,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v8,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v9,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v10,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v11,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v12,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v13,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v14,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v15,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v16,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v17,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v18,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v19,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v20,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v21,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v22,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v23,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v24,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v25,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v26,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v27,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v28,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v29,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v30,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v31,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v32,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v33,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v34,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v35,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v36,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v37,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v38,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v39,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v40,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v41,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v42,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v43,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v44,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v45,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v46,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v47,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v48,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v49,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v50,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v51,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v52,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v53,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v54,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v55,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v56,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v57,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v58,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v59,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v60,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v61,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v62,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v63,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v64,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v65,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v66,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v67,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v68,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v69,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v70,color=(1,0,1), material=Bound))

for b in O.bodies:
if b.material is Dolomite:
b.state.blockedDOFs='zXY'

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Facet_Aabb()],verletDist=0.05),
InteractionLoop(
[Ig2_Polyhedra_Polyhedra_PolyhedraGeom(),Ig2_Facet_Polyhedra_PolyhedraGeom()],
[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()],
[Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()]
),
NewtonIntegrator(gravity=(0,-9.81,0),damping=0.2),
]

utils.calm()

O.dt=10e-5
O.saveTmp()
O.run()

 Jan Stránský (honzik) said on 2020-03-03: #11

>> what values have you tried?
>O.dt=10e-8

what about 1e-10? 1e-12? 1e-15? 1e-20?

> If stiffness reaches less than 10e5, then particles are passed through facet.

what are the striffness values of both polyhedra and facets?
You can decrease stiffness of polyhedrons and increase stiffness of facets artificially (e.g. 1e300, the harmonic average in contact law makes it ok)

> Can you please provide a script (two bodies as is described) reproducing this behavior?

thanks, but it is not really a good script:
- there is a lot of not needed facets
- the two polyhedrons may not interact at all. E.g. with
Dolomite.young = 5e7
Bound.young = 1e300
O.dt=10e-5
the simulation behaves quite OK (because the polyhedrons does not interact)

for some other values, the interaction polyhedrons interact and "bounces excessively".

> it is a typical symptom of too large time step (making the simulation unstable) and I don't believe it is not overcome with decreasing the time step to a **sufficient** value

this still holds. It seems O.dt=1e-8 is still too large.

cheers
Jan

 Mahdeyeh (mahdiye.sky) said on 2020-03-03: #12

I really appreciate your patience Jan.

> Can you please provide a script (two bodies as is described) reproducing this behavior?
> but it is not really a good script

I corrected the codes in question and two polyhedrons have interact (With using your suggestion for young modules (Dolomite.young = 5e7 Bound.young = 1e300) and O.dt=10e-6). But when there are more polygons, problem appears.

> You can decrease stiffness of polyhedrons and increase stiffness of facets artificially (e.g. 1e300, the harmonic average in contact law makes it ok)

You know, this works for two polyhedrons! When I checked with more polyhedrons, the problem was still there.

> what about 1e-10? 1e-12? 1e-15? 1e-20?
In fact, I have no idea because if problem is really solved with these time steps, simulation will be impossible for me because my original code have 920 particles and it takes a long time.

 Jan Stránský (honzik) said on 2020-03-03: #13

Yes, DEM is expensive.. there are several option you can try:
- improve current polyhedrons implementation
- try different shapes, like PFacet, PotentialParticle or PotentialBlock
- try totally different method (some event-driven particle simulation
- try different software

In any of the cases, I have no idea how much you can gain..

cheers
Jan

 Mahdeyeh (mahdiye.sky) said on 2020-03-04: #14

> - try different shapes, like PFacet, PotentialParticle or PotentialBlock

Thanks a lot Jan , You right.

I want to try clump, hope to be successful.

Again Thanks for all your helps.

 Jan Stránský (honzik) said on 2020-03-04: #15

I would like to emphasize that I really have no idea how much (if any) performance you can gain.
On the other hand, to get some basic feeling, it should not be too difficult and time consuming to try e.g. clumped PFacets (although I have no personal experience with them).
cheers
Jan

 Mahdeyeh (mahdiye.sky) said on 2020-03-04: #16

> I would like to emphasize that I really have no idea how much (if any) performance you can gain.
> On the other hand, to get some basic feeling, it should not be too difficult and time consuming to try e.g. clumped PFacets

Yes Jan, I understand what you mean. I'll try it.

Best regards,
Mahdeyeh

 Launchpad Janitor (janitor) said on 2020-03-19: #17

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

 Mahdeyeh (mahdiye.sky) said on 2020-03-21: #18

Dear Jan Stransky

> - try totally different method (some event-driven particle simulation

is there event-driven particle simulation in the YADE right now?

 Jan Stránský (honzik) said on 2020-03-21: #19

no..
cheers
Jan

 Mahdeyeh (mahdiye.sky) said on 2020-03-21: #20

Thanks Jan Stransky.

I appreciate for introducing if you know any software with the ability of event-driven particle simulation.

Best regards,
Mahdeyeh