Yade

ptonparticle2 search exceeded 50000 iterations! step:0 0 0

Asked by Liu Zhiruo

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!

Part of my code is displayed here:

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
        color=[0,0.4,0.5]
        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
        color=[0,0.4,0.5]
        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)
        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
        color=[0,0.4,0.5]
        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
        color=[0,0.4,0.5]
        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]
        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

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Vasileios Angelidakis
Solved:
Last query:
Last reply:
Revision history for this message
Best Vasileios Angelidakis (vsangelidakis) said : 		#1

Hi Zhiruo,

First of all, in the future, please try to send Minimal Working Scripts, as described in [1].

Despite this, I can give you some pointers:

1. The main issue with this script, is that you have defined the particles "g" with huge overlaps, while they are dynamic bodies (i.e. fixed=False). So, the contact detection algorithm cannot handle this. In the DEM, we usually deal with small overlaps among the particles. That is, I would advise you to generate the granular material of the bed in elevated positions where the particles are non-overlapping, and let them deposit under gravity inside the box. Then, you can have the ellipsoid impact on them.

2. The second issue lies in the fact that you have 22258 bodies in your simulation, all of which interact with each other. The bodies with mask=3 forming the box are not dynamic, and thus you do not need to perform primary contact detection among them, while you want them to interact with the granular material (g) and the ellipsoid (both of which have a mask=1 by default ). You can stop detecting contact where you don't need it, by assigning a bitmask: avoidSelfInteractionMask=2 inside the collider (line 78 of your script), like this:

InsertionSortCollider([PotentialBlock2AABB()],verletDist=0.00,avoidSelfInteractionMask=2),

You can read the documentation [2] to see how this bitmask works. In essence, using what I suggest, you detect contact between particles with masks 1-1 and 1-3, but you do not detect contact between particles with masks 3-3.

3. Also, I would advise you use a non-zero verletDist, as a fraction of the mean particle dimension, but this is not a critical problem right now.

4. You use a large timestep: The function O.PWaveTimeStep() only works for spheres (or its equivalent, for polyhedra). I would advise you estimate the timestep manually, using the contact stiffness and the mass of the smallest dynamic body (ellipsoid or granular material).

5. To simulate the box, you don't need all these overlapping particles. Would advise you to use one cuboid to simulate each face.

6. Also, for the Potential Blocks, you no longer need to assign the volume and inertia manually. You can have a look in [3].

7. Last, for this number of bodies, I don't think the PotentialBlockVTKRecorder will manage to cope. Use instead the new feature export.exportPotentialBlocks [4] inside a PyRunner statement, if you need a VTK output.

Hope this helps.
All the best,
Vasileios

[1] https://www.yade-dem.org/wiki/Howtoask
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Collider
[3] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/PotentialBlocks/cubePBscaled.py#L70
[4] https://gitlab.com/yade-dev/trunk/-/blob/master/py/export.py#L902

Revision history for this message
Liu Zhiruo (luckyliu) said : 		#2

Hi Vasileios,
Thanks for your advise. It really work when I edit my code according to the pointer one given by you above.
However, when I want to perfect my code following your pointers two and five,the same error occured again! Even though I have tried many times, it still fails.

Part of my critical code is below:

# Boundary Plates
thickness = 40/10.0
wire=False
highlight=False
kPP = 0.0
rPP = 0.2
#
Igeom = array([0.08406805308160123, 0.1994387425401491, 0.2037312986056906]) #----------------------------------------------------------------------------------------------------------------------------------------------- #
# Bottom plate
bb= Body()
bb.mask=3
color=[0,0,1]
aPP = [1,-1,0,0,0,0]
bPP = [0,0,1,-1,0,0]
cPP = [0,0,0,0,1,-1]
dPP = [100.0-rPP , 100.0-rPP, 100-rPP, 100-rPP, thickness-rPP, thickness-rPP]
minmaxAabb = 1.05*Vector3( dPP[0], dPP[2], dPP[4] )
bb.shape=PotentialBlock(k=kPP, r=rPP, R=0.0, a=aPP, b=bPP, c=cPP, d=dPP, id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=minmaxAabb, maxAabb=minmaxAabb, fixedNormal=True, boundaryNormal=Vector3(0,0,-1))
utils._commonBodySetup(bb,bb.shape.volume,Igeom,
material='rockfall_material',pos= [0,0,0], fixed=True)
#bbb.state.ori = bbb.shape.orientation
O.bodies.append(bb)

# Lateral plate A
bA = Body()
bA.mask=3
color=[0,0.5,1]
aPP = [1,-1,0,0,0,0]
bPP = [0,0,1,-1,0,0]
cPP = [0,0,0,0,1,-1]
dPP = [100-rPP , 100-rPP, thickness-rPP, thickness-rPP, 0.5*50-rPP, 0.5*50-rPP]
minmaxAabb = 1.05*Vector3( dPP[0], dPP[2], dPP[4] )
bA.shape=PotentialBlock(k=kPP, r=rPP, R=0.0, a=aPP, b=bPP, c=cPP, d=dPP, id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=minmaxAabb, maxAabb=minmaxAabb, fixedNormal=True, boundaryNormal=Vector3(0,-1,0))
utils._commonBodySetup(bA,bA.shape.volume,Igeom,
material='rockfall_material',pos= [0,100+thickness,25-thickness], fixed=True)
#bA.state.ori = bA.shape.orientation
O.bodies.append(bA)

# ----------------------------------------------------------------------------------------------------------------------------------------------- #
# Lateral plate B
bB = Body()
bB.mask=3
color=[0,0.5,1]
aPP = [1,-1,0,0,0,0]
bPP = [0,0,1,-1,0,0]
cPP = [0,0,0,0,1,-1]
dPP = [100-rPP, 100-rPP, thickness-rPP, thickness-rPP, 0.5*50-rPP, 0.5*50-rPP]
minmaxAabb = 1.05*Vector3( dPP[0], dPP[2], dPP[4] )
bB.shape=PotentialBlock(k=kPP, r=rPP, R=0.0, a=aPP, b=bPP, c=cPP, d=dPP, id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=minmaxAabb, maxAabb=minmaxAabb, fixedNormal=True, boundaryNormal=Vector3(0,1,0))
utils._commonBodySetup(bB,bB.shape.volume,Igeom,
material='rockfall_material',pos= [0,-100-thickness,25-thickness], fixed=True)
#bB.state.ori = bB.shape.orientation
O.bodies.append(bB)

# ----------------------------------------------------------------------------------------------------------------------------------------------- #
# Lateral plate C
bC = Body()
bC.mask=3
color=[0,0.5,1]
aPP = [1,-1,0,0,0,0]
bPP = [0,0,1,-1,0,0]
cPP = [0,0,0,0,1,-1]
dPP = [thickness-rPP, thickness-rPP, 100-rPP, 100-rPP, 0.5*50-rPP, 0.5*50-rPP]
minmaxAabb = 1.05*Vector3( dPP[0], dPP[2], dPP[4] )
bC.shape=PotentialBlock(k=kPP, r=rPP, R=0.0, a=aPP, b=bPP, c=cPP, d=dPP, id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=minmaxAabb, maxAabb=minmaxAabb, fixedNormal=True, boundaryNormal=Vector3(1,0,0))
utils._commonBodySetup(bC,bC.shape.volume,Igeom,
material='rockfall_material',pos= [-100-thickness,0,25-thickness], fixed=True)
#bC.state.ori = bC.shape.orientation
O.bodies.append(bC)

# ----------------------------------------------------------------------------------------------------------------------------------------------- #
# Lateral plate D
bD = Body()
bD.mask=3
color=[0,0.5,1]
aPP = [1,-1,0,0,0,0]
bPP = [0,0,1,-1,0,0]
cPP = [0,0,0,0,1,-1]
dPP = [thickness-rPP, thickness-rPP, 100-rPP, 100-rPP, 0.5*50-rPP, 0.5*50-rPP]
minmaxAabb = 1.05*Vector3( dPP[0], dPP[2], dPP[4] )
bD.shape=PotentialBlock(k=kPP, r=rPP, R=0.0, a=aPP, b=bPP, c=cPP, d=dPP, id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=minmaxAabb, maxAabb=minmaxAabb, fixedNormal=True, boundaryNormal=Vector3(-1,0,0))
utils._commonBodySetup(bD,bD.shape.volume,Igeom,
material='rockfall_material',pos= [100+thickness,0,25-thickness], fixed=True)
#bD.state.ori = bD.shape.orientation
O.bodies.append(bD)

#granular

sp=pack.SpherePack()
distanceToCentre = 0.5
mn,mx=Vector3(90,90,50),Vector3(-90,-90,80)
sp.makeCloud(mn,mx,0.5,0.5,1000,False)

r=0.01
for s in sp:
 bb=Body()
 bb.mask=1
 bb.aspherical=True
 wire=False
 color=Vector3(random.random(),random.random(),random.random())
 highlight=False
 bb.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1],
d=[distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r], isBoundary=False, color=color,
wire=wire, highlight=highlight, minAabb=sqrt(3)*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), maxAabb=sqrt(3)*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), AabbMinMax=True, fixedNormal=False, id=len(O.bodies))
 utils._commonBodySetup(bb, bb.shape.volume, Igeom, material='rockfall_material',pos=s[0], fixed=False)
 b.state.pos = s[0] #s[0] stores center
 b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
 O.bodies.append(bb)

Revision history for this message
Vasileios Angelidakis (vsangelidakis) said : 		#3

Hi Zhiruo,

The error "ptonparticle2 search exceeded 50000 iterations! step:0 0 0" usually occurs when the overlap between two particles is too large. I cannot replicate the error in your latest message (#2), but would advise you to check that the particles you generate with makeCloud are not overlapping initially and so, generate them inside non-overlapping spheres. I.e, for a cube with an edge of 5, the radius of the smallest bounding sphere is 5/2*sqrt(3).

I see that in #2 you generated a larger box with dimensions around 100*100*4. Can you please try to generate a box with the dimensions of the original question?

Giving you answers on two different problems in the same thread will confuse us both :)

Cheers,
Vasileios

Revision history for this message
Liu Zhiruo (luckyliu) said : 		#4

Hi Vasileios,
Thanks for your suggestion again!
I download the newest edition of YADE on github which surprisingly solved my problem. So I guess there may be some different functions of different edition of YADE.
But I really learnt a lot about YADE because of your careful help which also inspired me to continue to learn it deeply!

Revision history for this message
Liu Zhiruo (luckyliu) said : 		#5

Thanks Vasileios Angelidakis, that solved my question.