Re: [Yade-users] [Question #289662]: ZeroDivisionError: float division by zero

[Yor1]

Question #289662 on Yade changed:
https://answers.launchpad.net/yade/+question/289662

Status: Answered => Solved

Yor1 confirmed that the question is solved:
Hi Jerome

Thank you for your advice.
I will try to put the essential part of the script without comments :)

Jabrane.

-- 
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

Re: [Yade-users] [Question #289662]: ZeroDivisionError: float division by zero

[Jérôme Duriez]

Question #289662 on Yade changed:
https://answers.launchpad.net/yade/+question/289662

Jérôme Duriez proposed the following answer:
Hi Jabrane,

The error says "Division by 0" thus it means Yade has at one point to
divide by 0 and obviously it can not handle that. The error also says
this division by 0 happens when such operation is made:

E+=0.5*(i.phys.normalForce.squaredNorm()/i.phys.kn +
i.phys.shearForce.squaredNorm()/i.phys.ks)

=> Yade encountered either one zero i.phys.kn, or one zero i.phys.ks.

The solution is for you to modify your mechanical parameters to avoid this, or 
to modify the use of your elastic energy formula. An interaction with zero 
stifness has certainly zero force hence zero elastic energy. No need to try 
"0/0" to obtain 0.
You could check all this asking Yade to give you more details when it faces a 
zero kn or ks (which interaction ? between who and who ?...)


As a general remark, your problem has been without solutions for 2 weeks now, 
which might have been disappointing to you. If you would like to improve you 
chances getting help, please do your best to improve a lot your question 
compared with this one, in the "minimal working script" direction : honestly, I 
think very few Yade users test the script of others when they are more than few 
lines long. In any case, here with your scripts and this 20 000 lines long 
.sphere files, it seems noone did it. 

It would have been straightforward for you to remove all the commented
lines from the script, or everything what concerns colors, or also to
state directly that "identifBis.py" is taken directly from the trunk
version (or to state the differences, if any).


It would much much more entice people to give a close look to your problem (and 
make it much easier).

Jerome

-- 
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

Re: [Yade-users] [Question #289662]: ZeroDivisionError: float division by zero

[Jan Stránský]

Question #289662 on Yade changed:
https://answers.launchpad.net/yade/+question/289662

Status: Open => Answered

Jan Stránský proposed the following answer:
Hi Jabrane,

please provide also the file X1Y2Z1_20k.spheres such that we can test it
ourselves.

Also the error says exactly on which line the error occurred, could you
please send the line alone in next message (such that we do not need to go
through the whole file)? Or just verify that the error is really at the
energy computation line and not related to some other divisions in the
scripts.

Thanks
Jan


2016-03-31 11:43 GMT+02:00 Yor1 :

> New question #289662 on Yade:
> https://answers.launchpad.net/yade/+question/289662
>
> Hello,
>
> I try to simulate triaxial compression test on fractured sample in order
> to determine the energy balance.
> To calculate the elastic energy i use this formula:
> Elastic energy=Fn^2/kn + Fs^2/ks with Fn, Fs normal and shear forces and
> kn, ks normal and shear stiffnesses.
>
> When i run the code i have this error message:
> ZeroDivisionError: float division by zero.
>
> Best regards.
> Jabrane.
>
> This is the code:
>
> from yade import ymport, utils , plot
> import math
>
> PACKING='X1Y2Z1_20k'
> OUT=PACKING+'_compressionTest'
>
>
> DAMP=0.4 # numerical damping
> saveData=100 # data record interval
> iterMax=200 # maximum number of iteration (to be adjusted)
> saveVTK=1 # Vtk files record interval
>
>
> confinement=-1e6
> #uniaxial_stress=-1e6
> #delta_stress=-1e6
> #stress_max=-200e6
> strainRate=-0.02
>
>
> intR=1.455
> DENS=4000
> YOUNG=65e9
> FRICT=10
> ALPHA=0.4
> TENS=8e6
> COH=160e6
>
>
> def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson =
> ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
> def wallMat(): return
> JCFpmMat(type=0,density=DENS,young=YOUNG,frictionAngle=radians(0))
>
>
>
> O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))
>
>
> dim=utils.aabbExtrema()
> xinf=dim[0][0]
> xsup=dim[1][0]
> X=xsup-xinf
> yinf=dim[0][1]
> ysup=dim[1][1]
> Y=ysup-yinf
> zinf=dim[0][2]
> zsup=dim[1][2]
> Z=zsup-zinf
>
> R=0
> Rmax=0
> numSpheres=0.
> for o in O.bodies:
>  if isinstance(o.shape,Sphere):
>numSpheres+=1
>R+=o.shape.radius
>if o.shape.radius>Rmax:
>  Rmax=o.shape.radius
> Rmean=R/numSpheres
>
> O.reset()
>
>
> mn,mx=Vector3(xinf+0.1*Rmean,yinf+0.1*Rmean,zinf+0.1*Rmean),Vector3(xsup-0.1*Rmean,ysup-0.1*Rmean,zsup-0.1*Rmean)
>
> walls=utils.aabbWalls(oversizeFactor=1.5,extrema=(mn,mx),thickness=min(X,Y,Z)/100.,material=wallMat)
> wallIds=O.bodies.append(walls)
>
>
> O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))
>
> import gts
>
> c=X/4
> alpha=math.pi/4
>
> ptA =  gts.Vertex( X/2. + c/2.*cos(alpha), Y/2. + c/2.*sin(alpha), 8./7.*Z)
> ptB =  gts.Vertex( X/2. - c/2.*cos(alpha), Y/2. - c/2.*sin(alpha), 8./7.*Z)
> ptApr = gts.Vertex(X/2. + c/2.*cos(alpha), Y/2. + c/2.*sin(alpha),
> -1./7.*Z)
> ptBpr = gts.Vertex(X/2. - c/2.*cos(alpha), Y/2. - c/2.*sin(alpha),
> -1./7.*Z)
> e1 = gts.Edge(ptA,ptB)
> e2 = gts.Edge(ptA,ptApr)
> e3 = gts.Edge(ptApr,ptB)
> f1 = gts.Face(e1,e2,e3)
>
> e4 = gts.Edge(ptB,ptBpr)
> e5 = gts.Edge(ptBpr,ptApr)
> f2 = gts.Face(e4,e5,e3)
>
> s1 = gts.Surface()
> s1.add(f1)
> s1.add(f2)
> facet = gtsSurface2Facets(s1,wire = False,material=wallMat)
> O.bodies.append(facet)
>
> execfile('identifBis.py')
>
> O.engines=[
> ForceResetter(),
>
> InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
> InteractionLoop(
>
> [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Box_Sphere_ScGeom()],
>
> [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
>
> [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
> ),
> TriaxialStressController(internalCompaction=False,label='triax'),
>
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.4,
> defaultDt=0.1*utils.PWaveTimeStep()),
> NewtonIntegrator(damping=DAMP,label="newton"),
>
> PyRunner(iterPeriod=int(saveData),initRun=True,command='recorder()',label='data'),
>
> VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks'],Key=OUT,label='vtk')
> ]
>
>
> tensCks=shearCks=cks=cks0=0
>
> def recorder():
>
> E=0
>
> global tensCks, shearCks, e10,e20,e30
> tensCks=0
> shearCks=0
> e10=0
> e20=0
> e30=0
>
> for i in O.interactions:
>  if not i.isReal : continue
>  if (isinstance(O.bodies[i.id1].shape,Sphere) and
> isinstance(O.bodies[i.id2].shape,Sphere)) or
> (isinstance(O.bodies[i.id1].shape,Box) and
> isinstance(O.bodies[i.id2].shape,Sphere)) or
> (isinstance(O.bodies[i.id1].shape,Sphere) and
> 

[Yade-users] [Question #289662]: ZeroDivisionError: float division by zero

[Yor1]

New question #289662 on Yade:
https://answers.launchpad.net/yade/+question/289662

Hello,

I try to simulate triaxial compression test on fractured sample in order to 
determine the energy balance.
To calculate the elastic energy i use this formula:
Elastic energy=Fn^2/kn + Fs^2/ks with Fn, Fs normal and shear forces and kn, ks 
normal and shear stiffnesses.

When i run the code i have this error message:
ZeroDivisionError: float division by zero.

Best regards.
Jabrane.

This is the code:

from yade import ymport, utils , plot
import math

PACKING='X1Y2Z1_20k'
OUT=PACKING+'_compressionTest'


DAMP=0.4 # numerical damping
saveData=100 # data record interval
iterMax=200 # maximum number of iteration (to be adjusted)
saveVTK=1 # Vtk files record interval


confinement=-1e6
#uniaxial_stress=-1e6
#delta_stress=-1e6
#stress_max=-200e6
strainRate=-0.02


intR=1.455
DENS=4000
YOUNG=65e9 
FRICT=10
ALPHA=0.4
TENS=8e6 
COH=160e6


def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson = 
ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
def wallMat(): return 
JCFpmMat(type=0,density=DENS,young=YOUNG,frictionAngle=radians(0))


O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))


dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf

R=0
Rmax=0
numSpheres=0.
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   numSpheres+=1
   R+=o.shape.radius
   if o.shape.radius>Rmax:
 Rmax=o.shape.radius
Rmean=R/numSpheres

O.reset() 

mn,mx=Vector3(xinf+0.1*Rmean,yinf+0.1*Rmean,zinf+0.1*Rmean),Vector3(xsup-0.1*Rmean,ysup-0.1*Rmean,zsup-0.1*Rmean)
walls=utils.aabbWalls(oversizeFactor=1.5,extrema=(mn,mx),thickness=min(X,Y,Z)/100.,material=wallMat)
wallIds=O.bodies.append(walls)

O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))

import gts

c=X/4
alpha=math.pi/4

ptA =  gts.Vertex( X/2. + c/2.*cos(alpha), Y/2. + c/2.*sin(alpha), 8./7.*Z)
ptB =  gts.Vertex( X/2. - c/2.*cos(alpha), Y/2. - c/2.*sin(alpha), 8./7.*Z)
ptApr = gts.Vertex(X/2. + c/2.*cos(alpha), Y/2. + c/2.*sin(alpha), -1./7.*Z)
ptBpr = gts.Vertex(X/2. - c/2.*cos(alpha), Y/2. - c/2.*sin(alpha), -1./7.*Z)
e1 = gts.Edge(ptA,ptB)
e2 = gts.Edge(ptA,ptApr)
e3 = gts.Edge(ptApr,ptB)
f1 = gts.Face(e1,e2,e3)

e4 = gts.Edge(ptB,ptBpr)
e5 = gts.Edge(ptBpr,ptApr)
f2 = gts.Face(e4,e5,e3)

s1 = gts.Surface()
s1.add(f1)
s1.add(f2)
facet = gtsSurface2Facets(s1,wire = False,material=wallMat)
O.bodies.append(facet)

execfile('identifBis.py')

O.engines=[
ForceResetter(),

InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
InteractionLoop(

[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Box_Sphere_ScGeom()],

[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],

[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
),
TriaxialStressController(internalCompaction=False,label='triax'),

GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.4,
 defaultDt=0.1*utils.PWaveTimeStep()),
NewtonIntegrator(damping=DAMP,label="newton"),

PyRunner(iterPeriod=int(saveData),initRun=True,command='recorder()',label='data'),

VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks'],Key=OUT,label='vtk')
]


tensCks=shearCks=cks=cks0=0

def recorder():
   
E=0

global tensCks, shearCks, e10,e20,e30
tensCks=0
shearCks=0
e10=0
e20=0
e30=0 

for i in O.interactions:
 if not i.isReal : continue
 if (isinstance(O.bodies[i.id1].shape,Sphere) and 
isinstance(O.bodies[i.id2].shape,Sphere)) or 
(isinstance(O.bodies[i.id1].shape,Box) and 
isinstance(O.bodies[i.id2].shape,Sphere)) or 
(isinstance(O.bodies[i.id1].shape,Sphere) and 
isinstance(O.bodies[i.id2].shape,Box)):
E+=0.5*(i.phys.normalForce.squaredNorm()/i.phys.kn + 
i.phys.shearForce.squaredNorm()/i.phys.ks)

  
triax.stressMask=7
triax.goal1=confinement
triax.goal2=confinement
triax.goal3=confinement
triax.max_vel=0.01

while 1:
  if confinement==0: 
O.run(1000,True) # to stabilize the system
break
  O.run(100,True)
  unb=unbalancedForce()
  
meanS=abs(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
  print 'unbalanced force:',unb,' mean stress: ',meanS
  if unb<0.005 and abs(meanS-abs(confinement))/abs(confinement)<0.001:
O.run(1000,True) # to stabilize the system
e10=triax.strain[0]
e20=triax.strain[1]
e30=triax.strain[2]
break

triax.stressMask=5
triax.goal1=confinement
triax.goal2=strainRate