Re: [Yade-users] [Question #289662]: ZeroDivisionError: float division by zero
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
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
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 > isinstance(O.bodies[i.id2].shape,Box)): > E+=0.5*(i.phys.
[Yade-users] [Question #289662]: ZeroDivisionError: float division by zero
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 triax.goal3=confine