yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #22677
[Question #689771]: ptonparticle2 search exceeded 50000 iterations! step:0 0 0
New question #689771 on Yade:
https://answers.launchpad.net/yade/+question/689771
I want to use the Potential Blocks to build a open box in the simulation,but when I run the code,the system just told me this:"ptonparticle2 search exceeded 50000 iterations! step:0 0 0
" .
And this sentence is being repeated all the time and the time step still stops at zero.
This question only occurs when I use the PB module.
I sincerely hope that anyone can come to my rescue! Thanks a lot!
My code is here:
from yade import pack,qt
import math
import os
import errno
from numpy import array
from decimal import *
import sys
global orien
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Clear the directory where the VTK results will be written of all files (but not subdirectories).
# The directory is cleared, so files from previous runs do not interfere with files from new runs.
# If the directory does not exist, it is created.
import os, shutil
folder = './vtk_rockfall'
try:
os.mkdir(folder[2:])
except:
for the_file in os.listdir(folder):
file_path = os.path.join(folder, the_file)
try:
if os.path.isfile(file_path):
os.unlink(file_path)
except Exception as e:
print(e)
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
aa=array([-0.73629877, -0.73629877, -0.73629877, -0.73629877, -0.53666249, -0.53666249, -0.53666249, -0.53666249, -0.30648893, -0.30648893, -0.30648893, -0.30648893, -0.19906207, -0.19906207, -0.19906207, -0.19906207, -0.16970517, -0.16970517, -0.16970517, -0.16970517, -0.11430817, -0.11430817, -0.11430817, -0.11430817, -0.09863695, -0.09863695, -0.09863695, -0.09863695, -0.04037159, -0.04037159, -0.04037159, -0.04037159, 0.04037159, 0.04037159, 0.04037159, 0.04037159, 0.09863695, 0.09863695, 0.09863695, 0.09863695, 0.11430817, 0.11430817, 0.11430817, 0.11430817, 0.16970517, 0.16970517, 0.16970517, 0.16970517, 0.19906207, 0.19906207, 0.19906207, 0.19906207, 0.30648893, 0.30648893, 0.30648893, 0.30648893, 0.53666249, 0.53666249, 0.53666249, 0.53666249, 0.73629877, 0.73629877, 0.73629877, 0.73629877])
bb=array([-0.29291786, -0.29291786, 0.29291786, 0.29291786, -0.71717283, -0.71717283, 0.71717283, 0.71717283, -0.91738619, -0.91738619, 0.91738619, 0.91738619, -0.19118596, -0.19118596, 0.19118596, 0.19118596, -0.54751162, -0.54751162, 0.54751162, 0.54751162, -0.82601961, -0.82601961, 0.82601961, 0.82601961, -0.99176290, -0.99176290, 0.99176290, 0.99176290, -0.97998568, -0.97998568, 0.97998568, 0.97998568, -0.97998568, -0.97998568, 0.97998568, 0.97998568, -0.99176290, -0.99176290, 0.99176290, 0.99176290, -0.82601961, -0.82601961, 0.82601961, 0.82601961, -0.54751162, -0.54751162, 0.54751162, 0.54751162, -0.19118596, -0.19118596, 0.19118596, 0.19118596, -0.91738619, -0.91738619, 0.91738619, 0.91738619, -0.71717283, -0.71717283, 0.71717283, 0.71717283, -0.29291786, -0.29291786, 0.29291786, 0.29291786])
cc=array([-0.60996987, 0.60996987, -0.60996987, 0.60996987, -0.44458577, 0.44458577, -0.44458577, 0.44458577, -0.25390374, 0.25390374, -0.25390374, 0.25390374, -0.96115671, 0.96115671, -0.96115671, 0.96115671, -0.81940904, 0.81940904, -0.81940904, 0.81940904, -0.55192866, 0.55192866, -0.55192866, 0.55192866, -0.08171353, 0.08171353, -0.08171353, 0.08171353, -0.19493127, 0.19493127, -0.19493127, 0.19493127, -0.19493127, 0.19493127, -0.19493127, 0.19493127, -0.08171353, 0.08171353, -0.08171353, 0.08171353, -0.55192866, 0.55192866, -0.55192866, 0.55192866, -0.81940904, 0.81940904, -0.81940904, 0.81940904, -0.96115671, 0.96115671, -0.96115671, 0.96115671, -0.25390374, 0.25390374, -0.25390374, 0.25390374, -0.44458577, 0.44458577, -0.44458577, 0.44458577, -0.60996987, 0.60996987, -0.60996987, 0.60996987])
dd=array([0.73629877, 0.73629877, 0.73629877, 0.73629877, 0.63303657, 0.63303657, 0.63303657, 0.63303657, 0.54106540, 0.54106540, 0.54106540, 0.54106540, 0.48057836, 0.48057836, 0.48057836, 0.48057836, 0.48327944, 0.48327944, 0.48327944, 0.48327944, 0.48717828, 0.48717828, 0.48717828, 0.48717828, 0.49588145, 0.49588145, 0.49588145, 0.49588145, 0.48999284, 0.48999284, 0.48999284, 0.48999284, 0.48999284, 0.48999284, 0.48999284, 0.48999284, 0.49588145, 0.49588145, 0.49588145, 0.49588145, 0.48717828, 0.48717828, 0.48717828, 0.48717828, 0.48327944, 0.48327944, 0.48327944, 0.48327944, 0.48057836, 0.48057836, 0.48057836, 0.48057836, 0.54106540, 0.54106540, 0.54106540, 0.54106540, 0.63303657, 0.63303657, 0.63303657, 0.63303657, 0.73629877, 0.73629877, 0.73629877, 0.73629877])
Igeom = array([0.08406805308160123, 0.1994387425401491, 0.2037312986056906])
Volume = 0.9069255090831569
minAabbRotatedValue = array([1.000000, 0.500000, 0.500000])
maxAabbRotatedValue = array([1.000000, 0.500000, 0.500000])
R_bound = 1.0000000000000000
# Parameters table with DEFAULT values, used if the script is NOT run in batch mode
utils.readParamsFromTable(
#material parameters
Kn=2.0e9, # (N/m) Normal stiffness
Ks=6.0e8, # (N/m) Shear stiffness
newtonDamping =0.0, # Fix it. No damping proportional to acceleration
viscousDamping=0.1, # 10% damping proportional to the relative velocity of the particles in contact
frictionAngle=30.0, # (deg) Friction angle
density=2650.0, # (kg/m^3) Density. Standard nominal value for many rocks, e.g. granite. We can change it to 2000 kg/m3 or anything else you want
#kinematic parameters
velMagnitude=7.0, # (m/s) Magnitude of the initial (translational) velocity
impactAngle=60.0, # (deg) The anngle between the velocity and the plane
angVelMagnitude_X=3, #Assume the same angular velocity around 3 axises
angVelMagnitude_Y=3,
angVelMagnitude_Z=3,
oriX=0, # random.random(), Fix the orientation as Professor's advice
oriY=0, # random.random(),
oriZ=0, # random.random(),
oriAngle=0, # random.random(),
#case ID
ID=1,
noTableOk=True # use default values if not run in batch
)
from yade.params import table
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Engines
O.engines=[
ForceResetter(),
InsertionSortCollider([PotentialBlock2AABB()],verletDist=0.00),
InteractionLoop(
[Ig2_PB_PB_ScGeom()],
[Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=table.Kn, ks_i=table.Ks, Knormal=table.Kn, Kshear=table.Ks, viscousDamping=table.viscousDamping)],
[Law2_SCG_KnKsPBPhys_KnKsPBLaw(label='law',neverErase=False)]
),
NewtonIntegrator(damping=table.newtonDamping, exactAsphericalRot=True, gravity=[0,0,0],label='newton'), #No gravity
PotentialBlockVTKRecorder(fileName=folder+'/rockfall_', iterPeriod=500, twoDimension=False, sampleX=60,sampleY=60,sampleZ=60, maxDimension=10, label='vtkRecorder')
#VTKRecorder(fileName='3d-vtk-',recorders=['all'],iterPeriod=1000),
#qt.SnapshotEngine(fileBase='3d-',iterPeriod=10000,label='snapshot')
]
#O.engines[1].avoidSelfInteractionMask=3 #Avoid contact detection for my case. But for your case this line may be deleted.
O.engines=O.engines+[PyRunner(iterPeriod=2,command='falling()')]+[PyRunner(iterPeriod=2,command='record()')]
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# The "young" and "poisson" parameters are not used. The normal and shear stiffnesses of all contacts are pre-defined within the IPhys functor
# The "frictionAngle" is used to calculate the maximum shear force before contacts slip
# The "density" is used to calculate the "mass" and "inertia tensor" of the particles
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(table.frictionAngle),density=table.density,label='rockfall_material'))
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Rockfall particle body definition
#from rockfall_particle_PB import * # Import particle characteristics from the Matalb-generated file
count = 0 #used for particle id
b=Body()
color=[0,0.5,1.0]
r=0.2*R_bound # Fix it
b.aspherical=True
b.shape=PotentialBlock(k=0.0, r=r, R=R_bound,
a=aa, b=bb, c=cc, d=dd-r,
isBoundary=False,
color=color,
AabbMinMax=True,fixedNormal=False,id=count)
geomInert=Igeom
utils._commonBodySetup(b,Volume,geomInert, material='rockfall_material',pos=[0,0,0], fixed=False)
b.volume=Volume
b.state.pos = [0,3.5*R_bound,R_bound*3]
b.state.ori = Quaternion((table.oriX,table.oriY,table.oriZ),(table.oriAngle))
O.bodies.append(b)
count = count+1
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Boundary plate body definition. For your case, you may not make the bodie fixed.
edge=R_bound/6
#O.bodies.append(geom.facetBox((-5*edge+50*edge/2.0,-5*edge+50*edge/2.0,0),(48*edge/2.0,48*edge/2.0,0.5),wallMask=31))
for x in range(1,100):
for y in range(1,100):
bbb=Body()
bbb.mask=3 #Avoid contact detection because I am simulating a fixed layer of particles. But you may concern their motion so comment this line seems ok.
color=[0,0.4,0.5]
# It controls the size of the plate particle
r=0.2*edge
bbb.shape=PotentialBlock(k=0.0, r=r, R=edge,
a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1],
d=[edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r],
isBoundary=True,
color=color,
minAabb = 1.00*Vector3(edge/2.0,edge/2.0,edge/2.0),
maxAabb = 1.00*Vector3(edge/2.0,edge/2.0,edge/2.0),
AabbMinMax=True,fixedNormal=False,id=count)
V=edge**3.0
geomInert=(1./6.)*V*edge**2.0
utils._commonBodySetup(bbb,V,Vector3(geomInert,geomInert,geomInert), material='rockfall_material',pos=[0,0,0], fixed=True) #You may make fixe=False and make them deposit in an open box
bbb.volume=V
bbb.state.pos = [-5*edge+x*edge/2.0,-5*edge+y*edge/2.0,0]
bbb.state.ori = Quaternion((1,1,0),radians(45))
O.bodies.append(bbb)
count = count+1
for y in range(1,100):
for z in range(1,12):
ccc=Body()
ccc.mask=3 #Avoid contact detection because I am simulating a fixed layer of particles. But you may concern their motion so comment this line seems ok.
color=[0,0.4,0.5]
# It controls the size of the plate particle
r=0.2*edge
ccc.shape=PotentialBlock(k=0.0, r=r, R=edge,
a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1],
d=[edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r],
#isBoundary=True,
color=color,
minAabb = 1.00*Vector3(edge/2.0,edge/2.0,edge/2.0),
maxAabb = 1.00*Vector3(edge/2.0,edge/2.0,edge/2.0),
AabbMinMax=True,fixedNormal=False,id=count)
V=edge**3.0
geomInert=(1./6.)*V*edge**2.0
utils._commonBodySetup(ccc,V,Vector3(geomInert,geomInert,geomInert), material='rockfall_material',pos=[0,0,0], fixed=True) #You may make fixe=False and make them deposit in an open box
ccc.volume=V
ccc.state.pos = [-5*edge+1*edge/2.0,-5*edge+y*edge/2.0,0+z*edge/2.0]
ccc.state.ori = Quaternion((1,1,0),radians(45))
O.bodies.append(ccc)
count = count+1
for y in range(1,100):
for z in range(1,12):
ddd=Body()
ddd.mask=3 #Avoid contact detection because I am simulating a fixed layer of particles. But you may concern their motion so comment this line seems ok.
color=[0,0.4,0.5]
# It controls the size of the plate particle
r=0.2*edge
ddd.shape=PotentialBlock(k=0.0, r=r, R=edge,
a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1],
d=[edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r],
#isBoundary=True,
color=color,
minAabb = 1.00*Vector3(edge/2.0,edge/2.0,edge/2.0),
maxAabb = 1.00*Vector3(edge/2.0,edge/2.0,edge/2.0),
AabbMinMax=True,fixedNormal=False,id=count)
V=edge**3.0
geomInert=(1./6.)*V*edge**2.0
utils._commonBodySetup(ddd,V,Vector3(geomInert,geomInert,geomInert), material='rockfall_material',pos=[0,0,0], fixed=True) #You may make fixe=False and make them deposit in an open box
ddd.volume=V
ddd.state.pos = [-5*edge+99*edge/2.0,-5*edge+y*edge/2.0,0+z*edge/2.0]
ddd.state.ori = Quaternion((1,1,0),radians(45))
O.bodies.append(ddd)
count = count+1
for x in range(1,100):
for z in range(1,12):
eee=Body()
eee.mask=3 #Avoid contact detection because I am simulating a fixed layer of particles. But you may concern their motion so comment this line seems ok.
color=[0,0.4,0.5]
# It controls the size of the plate particle
r=0.2*edge
eee.shape=PotentialBlock(k=0.0, r=r, R=edge,
a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1],
d=[edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r],
#isBoundary=True,
color=color,
minAabb = 1.00*Vector3(edge/2.0,edge/2.0,edge/2.0),
maxAabb = 1.00*Vector3(edge/2.0,edge/2.0,edge/2.0),
AabbMinMax=True,fixedNormal=False,id=count)
V=edge**3.0
geomInert=(1./6.)*V*edge**2.0
utils._commonBodySetup(eee,V,Vector3(geomInert,geomInert,geomInert), material='rockfall_material',pos=[0,0,0], fixed=True) #You may make fixe=False and make them deposit in an open box
eee.volume=V
eee.state.pos = [-5*edge+x*edge/2.0,-5*edge+1*edge/2.0,0+z*edge/2.0]
eee.state.ori = Quaternion((1,1,0),radians(45))
O.bodies.append(eee)
count = count+1
for x in range(1,100):
for z in range(1,12):
fff=Body()
fff.mask=3 #Avoid contact detection because I am simulating a fixed layer of particles. But you may concern their motion so comment this line seems ok.
color=[0,0.4,0.5]
# It controls the size of the plate particle
r=0.2*edge
fff.shape=PotentialBlock(k=0.0, r=r, R=edge,
a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1],
d=[edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r],
#isBoundary=True,
color=color,
minAabb = 1.00*Vector3(edge/2.0,edge/2.0,edge/2.0),
maxAabb = 1.00*Vector3(edge/2.0,edge/2.0,edge/2.0),
AabbMinMax=True,fixedNormal=False,id=count)
V=edge**3.0
geomInert=(1./6.)*V*edge**2.0
utils._commonBodySetup(fff,V,Vector3(geomInert,geomInert,geomInert), material='rockfall_material',pos=[0,0,0], fixed=True) #You may make fixe=False and make them deposit in an open box
fff.volume=V
fff.state.pos = [-5*edge+x*edge/2.0,-5*edge+99*edge/2.0,0+z*edge/2.0]
fff.state.ori = Quaternion((1,1,0),radians(45))
O.bodies.append(fff)
count = count+1
for x in range(5,95):
for y in range(5,95):
g=Body()
color=[0,0,0.3]
# It controls the size of the plate particle
r=0.1*random.random()
g.shape=PotentialBlock(k=0.0, r=r, R=edge,
a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1],
d=[edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r,edge/2.0-r],
#isBoundary=True,
color=color,
#minAabb = 1.00*Vector3(edge/2.0,edge/2.0,edge/2.0),
#maxAabb = 1.00*Vector3(edge/2.0,edge/2.0,edge/2.0),
AabbMinMax=True,fixedNormal=False,id=count)
V=edge**3.0
geomInert=(1./6.)*V*edge**2.0
utils._commonBodySetup(g,V,Vector3(geomInert,geomInert,geomInert), material='rockfall_material',pos=[0,0,0], fixed=False) #You may make fixe=False and make them deposit in an open box
g.volume=V
g.state.pos = [-5*edge+x*edge/2.0,-5*edge+y*edge/2.0,0.5]
g.state.ori = Quaternion((1,1,0),random.random())
O.bodies.append(g)
count = count+1
def activateGravity():
O.engines[3].gravity=Vector3(0,0,-9.81)
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Initial (translational) velocity
velMagnitude = table.velMagnitude # m/s
velDir=Vector3(tan(radians(90-table.impactAngle)),0,-1)
velDir.normalize()
velVector=velMagnitude*velDir
b.state.vel=velVector # Apply initial translational velocity on the rockfall particle
# Initial angular velocity
desiredAngVel=Vector3(table.angVelMagnitude_X,table.angVelMagnitude_Y,table.angVelMagnitude_Z)
test=numpy.diag(b.state.inertia)*(b.state.ori.conjugate()*desiredAngVel)
test1=Vector3(test[0][0],test[1][1],test[2][2])
b.state.angMom = b.state.ori*test1
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Time step.We need to ensure we use a small-enough timestep.
O.dt=.05*PWaveTimeStep()
#O.dt = 0.1*sqrt(0.3*O.bodies[0].state.mass/table.Ks)
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Function "activateGravity" to activate gravity upon impact. You may activate it in the IP functor right at the beginning of the simulation
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
previousVel=0
v0=velVector # Record velocities before first impact
firstImpact=True
isReboundVelocity=False
def falling():
global isReboundVelocity
global firstImpact
currentVel=b.state.vel[2]
if O.interactions.countReal()>0:
if firstImpact==True:
activateGravity()
firstImpact=False
if firstImpact==False:
if currentVel>0:
if O.interactions.countReal()==0:
O.engines[3].gravity=Vector3(0,0,0)
isReboundVelocity=True
def record():
global isReboundVelocity
global v0
#The logic is that when the impact starts the gravity is activated because we want a desired impact velocity. When the
#impact is complete we remove the gravity and record corrresponding data/here you can write any post-procesing or data recording commands
if O.interactions.countReal()==0:
if isReboundVelocity==True:
if b.state.pos[2]>1.2*R_bound:
v1=b.state.vel
w1=b.state.angVel
i1=b.state.inertia
ID1=table.ID
CoR_Tangential = math.sqrt(v1[0]**2+v1[1]**2)/math.sqrt(v0[0]**2+v0[1]**2)
CoR_Normal = -v1[2]/v0[2]
CoR_Overall = sqrt(v1[0]**2+v1[1]**2+v1[2]**2)/sqrt(v0[0]**2+v0[1]**2+v0[2]**2)
DispersionAngle=math.degrees(atan(v1[1]/v1[0]))
DispersionAngle_z=math.degrees(asin(abs(v1[2])/sqrt(v1[0]**2+v1[1]**2+v1[2]**2)))
TransEnergy_X=0.5*b.state.mass*v1[0]**2
TransEnergy_Y=0.5*b.state.mass*v1[1]**2
TransEnergy_Z=0.5*b.state.mass*v1[2]**2
TE=0.5*b.state.mass*(v1[0]**2+v1[1]**2+v1[2]**2)
TE_Original=0.5*b.state.mass*(v0[0]**2+v0[1]**2+v0[2]**2)
CoR_angVel_X=b.state.angVel[0]/table.angVelMagnitude_X
CoR_angVel_Y=b.state.angVel[1]/table.angVelMagnitude_Y
CoR_angVel_Z=b.state.angVel[2]/table.angVelMagnitude_Z
#f=open('./RESULTS.txt','a+')
#f.write('{}\t{}\t{}\t{}\n'.format(ID1,str(Decimal(CoR_Tangential)),str(Decimal(CoR_Normal)),str(Decimal(abs(CoR_Overall)))))
#f.close()
O.pause() # Pause scene
#makeVideo(snapshot.snapshots,'3d.mpeg',fps=10,bps=10000)
O.saveTmp()
from yade import qt
qt.Controller()
v=qt.View()
import yade.timing
O.timingEnabled = True
yade.timing.reset()
--
You received this question notification because your team yade-users is
an answer contact for Yade.
