Question #695935 on Yade changed:
https://answers.launchpad.net/yade/+question/695935
Yuxuan Wen gave more information on the question:
And then I run the second file. The "
Floating point exception (core dumped)" error will show when running this file:
#---------------------------------------------------------------------
## length (m), time (s), mass (1kg), force (N), pressure (Pa)
from yade import pack, plot, qt, export, os
O.load('consolidation.xml')
step0=O.iter
t0=O.time
rate=0.5
velocity=rate*0.1
radius=O.bodies[2].shape.radius
diameter=2*radius
## delete 2 up and low boxes
z1=O.bodies[0].state.pos[2]
z2=O.bodies[1].state.pos[2]
thickness=diameter
upzsurf=max(z1,z2)-0.5*thickness
lowzsurf=min(z1,z2)+0.5*thickness
lx=O.cell.size[0]
ly=O.cell.size[1]
lz=upzsurf-lowzsurf
V=lx*ly*lz
Vs0=O.cell.size[0]*O.cell.size[1]*O.cell.size[2]*(1-porosity())
O.bodies.erase(0,True)
O.bodies.erase(1,True)
## record all the sphere's id
spherelist=[]
msphere=0.0
for i in O.bodies: # bodies 0 and 1 are already deleted
spherelist.append(i.id)
msphere=msphere+i.state.mass
Vs=msphere/2650.0
## clump the walls
uplist=[]
lowlist=[]
clumplist=[]
for i in O.bodies:
if i.state.pos[2] >= (upzsurf-1.5*diameter):
uplist.append(i.id)
if i.state.pos[2] <= (lowzsurf+1.5*diameter):
lowlist.append(i.id)
clumplist=uplist+lowlist
upclump=O.bodies.clump(uplist)
lowclump=O.bodies.clump(lowlist)
upposy0=O.bodies[upclump].state.pos[1]
lowposy0=O.bodies[lowclump].state.pos[1]
## record the data of sphere body that is not in the clump walls
shearspherelist=[]
checklist=[]
for i in spherelist:
if i not in uplist and i not in lowlist:
shearspherelist.append(i)
checklist.append(i)
checklist.append(upclump)
checklist.append(lowclump)
#O.bodies[upclump].dynamic=False
#O.bodies[lowclump].dynamic=False
#O.bodies[upclump].state.vel=(0,velocity,0)
######################### define engines #########################
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
),
PyRunner(command='servo()',iterPeriod=1), # add force must before
NewtonIntegrator(), otherwise it will be ForceResetter()
GlobalStiffnessTimeStepper(),
NewtonIntegrator(gravity=(0,0,0),damping=0.2), # numerical damping,
dissipate Ekine because of quasi-static loading
PyRunner(command='if
O.iter-step0<=10000:addPlotData()',iterPeriod=100),
PyRunner(command='if
O.iter-step0>10000:addPlotData()',iterPeriod=1000),
PyRunner(command='if O.iter-step0<=100000:veldata()',iterPeriod=1000),
PyRunner(command='if O.iter-step0<=10000:fabricdata()',iterPeriod=100),
PyRunner(command='if O.iter-step0>10000:fabricdata()',iterPeriod=1000),
PyRunner(command='finished()',iterPeriod=1000),
]
plot.plots={'t':('sxx','syy','szz'),'t ':('upfx','upfy','upfz'),'t
':('upvx','upvy','upvz'),'t ':('Eunbal')}
plot.plot()
######################### define functions #########################
def servo():
O.bodies[upclump].dynamic=False
O.bodies[lowclump].dynamic=False
O.bodies[upclump].state.vel=(0,velocity,0) # can't blockDOF and add
velocity at the same time when dynamic is turned to False
def addPlotData():
upforcex=0
upforcey=0
upforcez=0
lowforcex=0
lowforcey=0
lowforcez=0
for i in O.interactions:
if not i.isReal: continue
point = i.geom.contactPoint
norm = i.geom.normal
ovp = i.geom.penetrationDepth
Fn = i.phys.normalForce
Fs = i.phys.shearForce
if (i.id1 in uplist and i.id2 not in uplist): # force direction
is from id1 to id2
upforcex=upforcex-Fn[0]-Fs[0]
upforcey=upforcey-Fn[1]-Fs[1]
upforcez=upforcez-Fn[2]-Fs[2]
if (i.id1 not in uplist and i.id2 in uplist):
upforcex=upforcex+Fn[0]+Fs[0]
upforcey=upforcey+Fn[1]+Fs[1]
upforcez=upforcez+Fn[2]+Fs[2]
if (i.id1 in lowlist and i.id2 not in lowlist):
lowforcex=lowforcex-Fn[0]-Fs[0]
lowforcey=lowforcey-Fn[1]-Fs[1]
lowforcez=lowforcez-Fn[2]-Fs[2]
if (i.id1 not in lowlist and i.id2 in lowlist):
lowforcex=lowforcex+Fn[0]+Fs[0]
lowforcey=lowforcey+Fn[1]+Fs[1]
lowforcez=lowforcez+Fn[2]+Fs[2]
n0=0
n1=0
for i in shearspherelist:
if len(O.bodies[i].intrs())==0:
n0=n0+1
if len(O.bodies[i].intrs())==1:
n1=n1+1
plot.addData(
N0=n0,
N1=n1,
step=O.iter,
t=O.time,
phi=Vs/V,
void=(V-Vs)/Vs,
sxx=utils.getStress(volume=V)[0][0],
sxy=utils.getStress(volume=V)[0][1], sxz=utils.getStress(volume=V)[0][2],
syx=utils.getStress(volume=V)[1][0],
syy=utils.getStress(volume=V)[1][1], syz=utils.getStress(volume=V)[1][2],
szx=utils.getStress(volume=V)[2][0],
szy=utils.getStress(volume=V)[2][1], szz=utils.getStress(volume=V)[2][2],
exx=O.cell.size[0], eyy=O.cell.size[1], ezz=O.cell.size[2],
upfx=upforcex, upfy=upforcey, upfz=upforcez,
lowfx=lowforcex, lowfy=lowforcey, lowfz=lowforcez,
upclumpfx=O.forces.f(upclump)[0],
upclumpfy=O.forces.f(upclump)[1], upclumpfz=O.forces.f(upclump)[2],
lowclumpfx=O.forces.f(lowclump)[0],
lowclumpfy=O.forces.f(lowclump)[1], lowclumpfz=O.forces.f(lowclump)[2],
upvx=O.bodies[upclump].state.vel[0],
upvy=O.bodies[upclump].state.vel[1], upvz=O.bodies[upclump].state.vel[2],
lowvx=O.bodies[lowclump].state.vel[0],lowvy=O.bodies[lowclump].state.vel[1],lowvz=O.bodies[lowclump].state.vel[2],
upposx=O.bodies[upclump].state.pos[0],upposy=O.bodies[upclump].state.pos[1],upposz=O.bodies[upclump].state.pos[2],
lowposx=O.bodies[lowclump].state.pos[0],lowposy=O.bodies[lowclump].state.pos[1],lowposz=O.bodies[lowclump].state.pos[2],
Z=avgNumInteractions(), # interactions between clump
particle-walls are also considered in YADE 2018.02b
Zm=avgNumInteractions(skipFree=True),
nshearsphere=len(shearspherelist), nclumpsphere=len(clumplist),
nsphere=len(spherelist),
Ekine=utils.kineticEnergy(),
Etot=utils.O.energy.total,
Eunbal=utils.unbalancedForce(),**O.energy # **O.energy converts
dictionary-like O.energy to plot.addData arguments
)
# overall state variables under different step are saved in one file,
no loop
plot.saveDataTxt(os.getcwd()+'/normalfiles/shearstate',vars=('t','step','Eunbal','phi','void','nshearsphere','nclumpsphere','nsphere'))
plot.saveDataTxt(os.getcwd()+'/normalfiles/shearforce',vars=('t','step','upfx','upfy','upfz','lowfx','lowfy','lowfz','upposy','upvy','lowposy','lowvy','N0','N1'))
#plot.saveDataTxt(os.getcwd()+'/normalfiles/shearflowclump',vars=('t','step','upclumpfx','upclumpfy','upclumpfz','lowclumpfx','lowclumpfy','lowclumpfz','upvx','upvy','upvz','lowvx','lowvy','lowvz'))
def veldata():
filename=os.getcwd()+'/velocitydata/'+'velocity'+str(O.iter)
f = open(filename,'w')
f.write('id vx vy vz pos_x pos_y pos_z ovp\n')
for i in checklist: # data of all bodies/intrs at each step is saved
under a unique file named str(O.iter)
vx=O.bodies[i].state.vel[0]
vy=O.bodies[i].state.vel[1]
vz=O.bodies[i].state.vel[2]
pos_x=O.bodies[i].state.pos[0]
pos_y=O.bodies[i].state.pos[1]
pos_z=O.bodies[i].state.pos[2]
f.write('%-8i %-12g %-12g %-12g %-12g %-12g
%-12g\n'%(i,vx,vy,vz,pos_x,pos_y,pos_z))
f.close
def fabricdata():
filename1=os.getcwd()+'/fabricdata/'+'fabric'+str(O.iter)
f1 = open(filename1,'w')
f1.write('id1 id2 x_cp y_cp z_cp n_x n_y n_z Fn_x Fn_y Fn_z Fs_x Fs_y
Fs_z ovp\n')
for i in O.interactions:
if not i.isReal: continue
point = i.geom.contactPoint
norm = i.geom.normal
ovp = i.geom.penetrationDepth
Fn = i.phys.normalForce
Fs = i.phys.shearForce
if not (i.id1 in clumplist and i.id2 in clumplist): # the
clumplist interactions should be constant
f1.write('%-8i %-8i %-12g %-12g %-12g %-12g %-12g %-12g
%-12g %-12g %-12g %-12g %-12g %-12g
%-12g\n'%(i.id1,i.id2,point[0],point[1],point[2],norm[0],norm[1],norm[2],Fn[0],Fn[1],Fn[2],Fs[0],Fs[1],Fs[2],ovp))
f1.close
filename2=os.getcwd()+'/N0N1data/'+'N0_id'+str(O.iter)
f2 = open(filename2,'w')
f2.write('N0_id\n')
for i in shearspherelist:
if len(O.bodies[i].intrs())==0:
f2.write('%-8i\n'%(i))
f2.close
filename3=os.getcwd()+'/N0N1data/'+'N1_id'+str(O.iter)
f3 = open(filename3,'w')
f3.write('N1_id\n')
for i in shearspherelist:
if len(O.bodies[i].intrs())==1:
f3.write('%-8i\n'%(i))
f3.close
def finished():
unbal=utils.unbalancedForce()
upvely=O.bodies[upclump].state.vel[1]
upposy=O.bodies[upclump].state.pos[1]
strain=(upposy-upposy0)/lz
if (O.time-t0)>(20*ly/velocity):
print ('Shearing Finished')
print ('Unbalanced force='+str(unbal))
print ('Shear strain='+str(strain))
O.pause()
--
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
