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. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp