Yade

polyhedra and facet interaction :error /InsertionSortcollider.cpp:61 insertionSortParallel

Asked by De zhang

I made a simulation of Polyhedra particles compaction. The boundary walls were formed by facet assemebly. I used utils.wall to simulate boundary before, but I found that the force chain of polyhdra and walls were all concentrated to the center of walls. Using utils.facet and polyhedra, I used "InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb()])", but after run(), it showed that:
 "ERROR /data/trunk/pkg/common/InsertionSortCollider.cpp:61 insertionSortParallel: Parallel insertion sort needs verletDist>0"
Though, it could keep runing, but I want to know what's the problem of this error message?
The following is my procedure:
####################################################################
from yade import pack,qt,plot,utils,polyhedra_utils,ymport,export,pack,timing
from yade import *
import numpy
from pprint import pprint
import random
from random import uniform
#from random import randint
import math
from math import *

##################################
#material:ballast and loadingplate
global gravel
global steel
gravel = PolyhedraMat()
gravel.density = 2600 #kg/m^3
gravel.young = 2.0E9 #5.5E8# 1e9
gravel.poisson = 0.21 # real 0.21
gravel.frictionAngle = 0.83 #rad radians(48) // change for rad math.radians(31)
#gravel.strength = 200.0e0 # Pa crushable 1000MPa
#gravel.IsSplitable = True
#The rock real young 5.5E10 (5.5e8 also acceptable)
steel = PolyhedraMat()
steel.density = 7850 #kg/m^3
steel.young = 8.0E9 #inital steel was 10*gravel.young
steel.poisson = 0.3
steel.frictionAngle = 0.55 #rad radians(31)
##next
#################################################################

##ID=4
# Area of the confining Wall
global A1,A2,A3,A4
global LoadPos,IniLoadPos,plateF,IniTime,forceA
global WallStress,ConfStress,ConfDevi,MoveVel,AxiDevi,MoveAxial,AreaPlate
#unit:m^2
IniTime=0
plateF=0.1 #Unit:kPa
LoadPos=0.6
IniLoadPos=LoadPos # (link to Area of Walls)
forceA=200 # Unit:kPa,P=N/A;N=P*0.0615*1000;A=0.0615
A1=LoadPos*0.3
A2=LoadPos*0.3
A3=LoadPos*0.3
A4=LoadPos*0.3
WallStress=0 # Unit:kPa
ConfStress=70 # Unit:kPa==> 70kPa---> 240kPa
ConfDevi=0
AxiDevi=0
MoveVel=0
MoveAxial=0
AreaPlate=0.09
# Id of different substances
global NumLoad_1,NumLoad_2,NumEndBall,StepNum,NumEnd
NumLoad_1=1
NumLoad_2=1
NumEndBall=1
StepNum=1
NumEnd=1
global xratio,yratio,zratio,NumContact,WallContact
xratio=1
yratio=1
zratio=1
NumContact=4
WallContact=1
# Position and other parameters record
######################parameters
global theta,thega,Vol
theta=0
thega=0
Vol=0

global EndPosDel,EndPos1,EndPos2
EndPosDel=0.1
EndPos1=0.5
EndPos2=0.6
##
global iternum
iternum=0
#bottom_wall=utils.wall(0,axis=2,sense=1,material=steel,mask=1)
#O.bodies.append(bottom_wall)
#bottom_wall.state.blockedDOFs='xyzXYZ'
for i in range(1,13):
 for j in range(1,13):
  O.bodies.append(utils.facet(vertices=((-0.18+0.03*(i-1),-0.18+0.03*(j-1),0),(-0.18+0.03*i,-0.18+0.03*(j-1),0),(-0.18+0.03*i,-0.18+0.03*j,0)),dynamic=None,fixed=True,wire=True,highlight=False,noBound=False,material=steel,mask=1,chain=-1))
  O.bodies.append(utils.facet(vertices=((-0.18+0.03*(i-1),-0.18+0.03*(j-1),0),(-0.18+0.03*(i-1),-0.18+0.03*j,0),(-0.18+0.03*i,-0.18+0.03*j,0)),dynamic=None,fixed=True,wire=True,highlight=False,noBound=False,material=steel,mask=1,chain=-1))
#Bottom Walls
#ID=12*12*2=288==>0-287
for i in range(1,13):
 for j in range(1,21):
  O.bodies.append(utils.facet(vertices=((-0.18+0.03*(i-1),-0.15,0.03*(j-1)),(-0.18+0.03*i,-0.15,0.03*(j-1)),(-0.18+0.03*i,-0.15,0.03*j)),dynamic=None,fixed=True,wire=True,highlight=False,noBound=False,material=steel,mask=1,chain=-1))
  O.bodies.append(utils.facet(vertices=((-0.18+0.03*(i-1),-0.15,0.03*(j-1)),(-0.18+0.03*(i-1),-0.15,0.03*j),(-0.18+0.03*i,-0.15,0.03*j)),dynamic=None,fixed=True,wire=True,highlight=False,noBound=False,material=steel,mask=1,chain=-1))
#Front Walls
#ID=12*20*2=480==>288-767
for i in range(1,13):
 for j in range(1,21):
  O.bodies.append(utils.facet(vertices=((-0.18+0.03*(i-1),0.15,0.03*(j-1)),(-0.18+0.03*i,0.15,0.03*(j-1)),(-0.18+0.03*i,0.15,0.03*j)),dynamic=None,fixed=True,wire=True,highlight=False,noBound=False,material=steel,mask=1,chain=-1))
  O.bodies.append(utils.facet(vertices=((-0.18+0.03*(i-1),0.15,0.03*(j-1)),(-0.18+0.03*(i-1),0.15,0.03*j),(-0.18+0.03*i,0.15,0.03*j)),dynamic=None,fixed=True,wire=True,highlight=False,noBound=False,material=steel,mask=1,chain=-1))
#Behind Walls
#ID=12*20*2=480==>768-1247
for i in range(1,13):
 for j in range(1,21):
  O.bodies.append(utils.facet(vertices=((-0.15,-0.18+0.03*(i-1),0.03*(j-1)),(-0.15,-0.18+0.03*i,0.03*(j-1)),(-0.15,-0.18+0.03*i,0.03*j)),dynamic=None,fixed=True,wire=True,highlight=False,noBound=False,material=steel,mask=1,chain=-1))
  O.bodies.append(utils.facet(vertices=((-0.15,-0.18+0.03*(i-1),0.03*(j-1)),(-0.15,-0.18+0.03*(i-1),0.03*j),(-0.15,-0.18+0.03*i,0.03*j)),dynamic=None,fixed=True,wire=True,highlight=False,noBound=False,material=steel,mask=1,chain=-1))
#left Walls
#ID=12*20*2=480==>1248-1727
for i in range(1,13):
 for j in range(1,21):
  O.bodies.append(utils.facet(vertices=((0.15,-0.18+0.03*(i-1),0.03*(j-1)),(0.15,-0.18+0.03*i,0.03*(j-1)),(0.15,-0.18+0.03*i,0.03*j)),dynamic=None,fixed=True,wire=True,highlight=False,noBound=False,material=steel,mask=1,chain=-1))
  O.bodies.append(utils.facet(vertices=((0.15,-0.18+0.03*(i-1),0.03*(j-1)),(0.15,-0.18+0.03*(i-1),0.03*j),(0.15,-0.18+0.03*i,0.03*j)),dynamic=None,fixed=True,wire=True,highlight=False,noBound=False,material=steel,mask=1,chain=-1))
#Right Walls
#ID=12*20*2=480==>1728-2207
ballast1=polyhedra_utils.fillBox((-0.15,-0.15,0), (0.15,0.15,0.6),gravel,sizemin=[0.025,0.025,0.025],sizemax=[0.036,0.036,0.036],ratio=[1,1,1],seed=4,mask=1)
#O.bodies.append(ballast1)
NumEndBall=O.bodies[-1].id
O.dt=1.15e-5 #Check it!
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
      [Ig2_Wall_Polyhedra_PolyhedraGeom(),Ig2_Polyhedra_Polyhedra_PolyhedraGeom(),Ig2_Facet_Polyhedra_PolyhedraGeom()],
      [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys(),Ip2_FrictMat_PolyhedraMat_FrictPhys()],
      [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()],
   ),
   NewtonIntegrator(damping=0.6,gravity=(0,0,-9.81),label='newton'),
   PyRunner(command='TraiStep1()',iterPeriod=1,label='checker1'),
   PyRunner(command='TraiStep2()',iterPeriod=1,label='checker2',dead=True),
   PyRunner(command='TraiStep3()',iterPeriod=1,label='checker3',dead=True),
   #PyRunner(command='LoadAxial()',iterPeriod=1,label='loadkeep',dead=True),
   PyRunner(command='AxialLoading()',iterPeriod=1,label='axialload',dead=True),
   PyRunner(command='addPlotData()',iterPeriod=1,label='plotdata',dead=True),
   #PyRunner(command='Confining()',iterPeriod=1,label='keepconf',dead=True),
   PyRunner(command='Monitor_1()',iterPeriod=5000,label='Monitor1',dead=True),
   PyRunner(command='Monitor_2()',iterPeriod=105000,label='Monitor2',dead=True),
]

##Fullfill the box
def TraiStep1():
 global WallStress,ConfStress,ConfDevi,MoveVel,AxiDevi,MoveAxial
 global A1
 global LoadPos,NumLoad_1,NumLoad_2,NumEndBall,IniLoadPos,plateF,IniTime,forceA,StepNum,NumEnd,iternum,AreaPlate
 global xratio,yratio,ztario,NumContact,WallContact
 global EndPosDel,EndPos1,EndPos2
 global theta,thega,Vol
 ######
 if StepNum == 1:
  ##### make circle dormetory
  ##Wall generation
  print("Come on! @ Dez")
  StepNum=StepNum+1
 elif StepNum == 2:
  O.dt=polyhedra_utils.PWaveTimeStep()
  print "2-UnForces = %.5f"%(utils.unbalancedForce())
  if O.iter < 5000: return
  if utils.unbalancedForce() > 1.0: return
  NumEnd=0.35
  ldplt=polyhedra_utils.polyhedra(steel,v=((-0.148,-0.148,NumEnd),(0.148,-0.148,NumEnd),(0.148,0.148,NumEnd),(-0.148,0.148,NumEnd),
   (-0.148,-0.148,NumEnd+0.02),(0.148,-0.148,NumEnd+0.02),(0.148,0.148,NumEnd+0.02),(-0.148,0.148,NumEnd+0.02)),fixed=False, color=(0.75,0.65,0.65))
  O.bodies.append(ldplt)
  NumLoad_1=O.bodies[-1].id
  O.bodies[NumLoad_1].state.blockedDOFs='xyXYZ'
  StepNum=StepNum+1
  #O.pause()
  #-- generation of first plate to make compaction of assembly --#
 elif StepNum == 3:
  #AreaPlate=(O.bodies[4].state.pos[0]-O.bodies[3].state.pos[0])*(O.bodies[4].state.pos[0]-O.bodies[3].state.pos[0])
  #AreaPlate=(O.bodies[1800].state.pos[0]-O.bodies[1300].state.pos[0])*(O.bodies[1800].state.pos[0]-O.bodies[1300].state.pos[0])
  plateF=(O.forces.f(NumLoad_1)[2])/90
  print "3-Load-force= %.5f"%(plateF)
  if plateF < 25: return
  O.bodies[NumLoad_1].state.blockedDOFs='xyzXYZ'
  StepNum=StepNum+1
  #O.pause()
 elif StepNum == 4:
  LoadPos=O.bodies[NumLoad_1].state.pos[2]
  print "4-Load-force= %.5f"%(plateF)
  #AreaPlate=(O.bodies[1800].state.pos[0]-O.bodies[1300].state.pos[0])*(O.bodies[1800].state.pos[0]-O.bodies[1300].state.pos[0])
  plateF=(O.forces.f(NumLoad_1)[2])/90 #P=F/A=F/(0.0615*1000)=F/61.5 Unit:kPa
  if plateF < 160:
   O.bodies[NumLoad_1].state.vel=(0,0,-0.005) #5mm/s
   O.bodies[NumLoad_1].state.blockedDOFs='xyzXYZ'
  else:
   O.bodies[NumLoad_1].state.vel=(0,0,0)
   O.bodies[NumLoad_1].state.blockedDOFs='xyzXYZ'
   StepNum=StepNum+1
   O.save('Step1.yade')
   O.pause()
   ###StepNum=4,NumLoad_1=O.bodies[-1].id,LoadPos=O.bodies[NumLoad_1].state.pos[2]

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
Last reply:
Revision history for this message
Best Jan Stránský (honzik) said : 		#1

Hello,

> but I found that the force chain of polyhdra and walls were all concentrated to the center of walls

you will get the same problem with facets.. also see [1].
With finer facets you improve this postprocessing problem, but increase number of bodies (simulation time) and introduce unneeded negative effects of facets (mainly wrong stiffness for a body in contact with multiple facets).
So this workaround works, but the clean way would be to use walls and update the postprocessing :-)

> "ERROR /data/trunk/pkg/common/InsertionSortCollider.cpp:61 insertionSortParallel: Parallel insertion sort needs verletDist>0"
> Though, it could keep running, but I want to know what's the problem of this error message?

It is an error from point of InsertionSortCollider (not the simulation itself). It is trying to run in parallel, but cannot because verletDist is set to -.5 by default [2], possibly automatically updated according to spheres in the simulation.
I **think** the error means that the collider will not run in parallel.
Either ignore the message or set manually:
InsertionSortCollider(..., verletDist=.05), # or another positive value

cheers
Jan

[1] https://answers.launchpad.net/yade/+question/668264
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.InsertionSortCollider.verletDist
[3] https://github.com/yade/trunk/blob/master/pkg/common/InsertionSortCollider.cpp#L61

Revision history for this message
De zhang (dzlearnyade) said : 		#2

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

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#3

Hi,
For verletDist if you can't rely on automatic settings (because you have no sphere) you could set it to e.g. +0.5*referencePolyhedronSize.
If you don't it will run but it will be slow.
It should actually be a warning more than an error (I suspect "error" is returning immediately, hence blocking the simulation?).

Bruno

Revision history for this message
De zhang (dzlearnyade) said : 		#4

Hi Bruno,
Thank you for your reply!
I tried verletDist with a little positive value (1e-10), it could run, and actually if with the error message for verletDist <0, it didn't stop running the programme.

Best regard!
De Zhang