← Back to team overview

yade-users team mailing list archive

  1. yade-users team
  2. Mailing list archive
  3. Message #17381

[Question #669734]: polyhedra and facet interaction :error /InsertionSortcollider.cpp:61 insertionSortParallel

  Thread Previous • Date Previous • • Thread Next 
New question #669734 on Yade:
https://answers.launchpad.net/yade/+question/669734

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]	


-- 
You received this question notification because your team yade-users is
an answer contact for Yade.
  Thread Previous • Date Previous • • Thread Next