← Back to team overview

yade-users team mailing list archive

[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.