New question #683384 on Yade: https://answers.launchpad.net/yade/+question/683384
Hi all, In my case, I want to simulate the progressive loss of fine particles in the soil assembly, the code for this case is attached as follows. After running 1000 steps, we can see the fines indeed get lost from the visualization. However, when using command "len(O.bodies)" in the terminal, we get the particle number remains unchanged. i.e. n=400. How can I monitor and output the instaneous particle number during simulation? Thanks very much! Best, Zheng ################# the code as follows #################### from yade import pack,plot,qt,export import matplotlib.pyplot as plt import numpy as np import random #O.materials.append(FrictMat(young=6.e8,poisson=.8,frictionAngle=.0)) O.materials.append(CohFrictMat(young=1e9,poisson=.8,frictionAngle=.0,normalCohesion=1e6,shearCohesion=1e6,momentRotationLaw=False,etaRoll=.1,label='spheres')) sigmaIso=-1e5 sp = pack.SpherePack() size =0.24 sp.makeCloud(minCorner=(0,0,.05),maxCorner=(size,size,.05),rMean=.005,rRelFuzz=.4,num=400,periodic=True,seed=1) # initial 400 #sp.makeCloud(minCorner=(0,0,.05),maxCorner=(size,size,.05),psdSizes=[0.006,0.0068,0.0072,0.0084,0.0092,0.01,0.0108,0.0116,0.0124,0.0132,0.014], psdCumm=[0,0.1,0.2, 0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0],periodic=True,seed=1,distributeMass=False)#num=400, # if minCorner[k]=maxCorner[k] for one coordinate, it means 2D case; sp.toSimulation() O.cell.hSize = Matrix3(size,0,0, 0,size,0, 0,0,.1) # used for periodic boundaries. #print len(O.bodies) for p in O.bodies: p.state.blockedDOFs = 'zXY' # a sphere can be made to move only in xy plane p.state.mass = 2650 * 0.1 * pi * p.shape.radius**2 # 0.1 = thickness of cylindrical particle inertia = 0.5 * p.state.mass * p.shape.radius**2 p.state.inertia = (.5*inertia,.5*inertia,inertia) O.dt = utils.PWaveTimeStep() print O.dt O.engines = [ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom6D()], [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()], [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(always_use_moment_law=False)] ), PeriTriaxController( dynCell=True, goal=(sigmaIso,sigmaIso,0), stressMask=3, relStressTol=.001, maxUnbalanced=.001, maxStrainRate=(.5,.5,.0), doneHook='shear()',# function to run when finished label='biax' ), NewtonIntegrator(damping=.1), PyRunner(command='delByNum()',iterPeriod=100), PyRunner(command='addPlotData()',iterPeriod=100), #PyRunner(command='deleteFines()',iterPeriod=100), #PyRunner(command='delBelowPerc()',iterPeriod=10000), #PyRunner(command='print ') ] plot.live=True plot.plots={'iter':('sxx','syy'),'iter_':('exx','eyy'),'iter':('poros')} def addPlotData(): plot.addData( iter=O.iter,iter_=O.iter, sxx=biax.stress[0],syy=biax.stress[1], exx=biax.strain[0],eyy=biax.strain[1], Z=avgNumInteractions(), Zm=avgNumInteractions(skipFree=True), poro=porosity(), unbalanced=utils.unbalancedForce(), t=O.time ) plot.saveDataTxt('results',vars=('t','exx','eyy','sxx','syy','Z','Zm','poro')) """ # delete fines by percent def delBelowPerc(): bodyRadius=[] for b in O.bodies: if b.shape.radius<=0.005: #if isinstance(b.shape,Sphere): bodyRadius.append([b,b.shape.radius]) bodyRadius.sort(key=lambda x:x[1]) if len(bodyRadius)>1: maxRad=np.percentile(bodyRadius,10) for b in bodyRadius: if b[0].shape.radius<=maxRad: O.bodies.erase(b[0].id) """ def delByNum(): setContactFriction(0.5) bodyRadius=[] for b in O.bodies: if b.shape.radius<=0.005: # define fine particles #if isinstance(b.shape,Sphere): bodyRadius.append([b,b.shape.radius]) bodyRadius.sort(key=lambda x:x[1]) l=len(bodyRadius) # how many fine particles perc=0.1 delNum=int(l*perc) # how many number to be deleted list=random.sample(bodyRadius,delNum) # list to be deleted for b in list: O.bodies.erase(b[0].id) print 'delete fines' #biax.doneHook='shear()' """ def deleteFines(): #new added here by Zheng 0422 bodiesToBeDeleted=[] for b in O.bodies: if b.shape.radius<=0.004: bodiesToBeDeleted.append(b) for b in bodiesToBeDeleted: O.bodies.erase(b.id) #new added here by Zheng 0422 def term0(): print getStress() biax.goal=(-10,-10,0) biax.doneHook='term()' def term1(): # delete a determined percent of fines after consolidation and then reconsolidation bodyRadius=[] for b in O.bodies: if b.shape.radius<=0.005: # define fine particles #if isinstance(b.shape,Sphere): bodyRadius.append([b,b.shape.radius]) bodyRadius.sort(key=lambda x:x[1]) l=len(bodyRadius) # how many fine particles perc=0.3 delNum=int(l*perc) # how many number to be deleted list=random.sample(bodyRadius,delNum) # list to be deleted for b in list: O.bodies.erase(b[0].id) # output PSD psd = utils.psd(bins=10,mass=True,mask=-1) print psd[0],psd[1] # lists of bins' sizes; cumulative percent of material plt.semilogx(psd[0],psd[1]) #plt.show() plt.savefig('psd.png') fpsd=file('psd.dat','a') fpsd.write('psd[0] psd[1]\n') fpsd.write(str(psd[0])+' '+str(psd[1])+'\n') fpsd.close() biax.doneHook='term()' """ def coh(): O.engines[2].physDispatcher.functors[0].setCohesionNow=True print 'add cohesion' biax.doneHook='shear()' # for biaxial shear below: def shear(): print getStress() #setContactFriction(0.5) # set the current cell configuration to be the reference one O.cell.trsf=Matrix3.Identity biax.goal=(sigmaIso,-0.1,0) biax.stressMask=1 # strain rate along y-axis is 0.01, use a larger x-axis strain rate to better maintian stresses biax.maxStrainRate=(0.5,0.01,0) biax.relStressTol=.01 biax.maxUnbalanced=.01 print 'shearing' biax.doneHook='biaxFinished()' def biaxFinished(): print 'Biaxial Test Finished' O.pause() """ def term(): O.engines = O.engines[:3]+O.engines[4:] print len(O.bodies) print getStress() print O.cell.hSize setContactFriction(0.5) O.cell.trsf=Matrix3.Identity O.cell.velGrad=Matrix3.Zero for p in O.bodies: p.state.vel = Vector3.Zero p.state.angVel = Vector3.Zero p.state.refPos = p.state.pos p.state.refOri = p.state.ori O.save('0.yade.gz') O.pause() """ O.run(1000) O.wait() -- 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