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 : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
