yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #12337
Re: [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 <question289662@xxxxxxxxxxxxxxxxxxxxx>:
> 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=2000000 # maximum number of iteration (to be adjusted)
> saveVTK=10000 # 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=confinement
> triax.max_vel=1
>
> O.run(iterMax)
>
>
> and this is the code names identifBis.py:
>
> interactionRadius=intR;
> O.engines=[
>
> ForceResetter(),
>
> InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius,label='is2aabb'),Bo1_Facet_Aabb()]),
> InteractionLoop(
>
> [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='ss2d3dg'),Ig2_Facet_Sphere_ScGeom()],
>
> [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
>
> [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=True,label='interactionLaw')]
> ),
> NewtonIntegrator(damping=1)
>
> ]
>
> ############################ timestep + opening yade windows
> O.dt=0.001*utils.PWaveTimeStep()
>
> # from yade import qt
> # v=qt.Controller()
> # v=qt.View()
>
> ############################ Identification spheres on joint
> #### color set for particles on joint
> jointcolor1=(0,1,0)
> jointcolor2=(1,0,0)
> jointcolor3=(0,0,1)
> jointcolor4=(1,1,1)
> jointcolor5=(0,0,0)
>
> #### first step-> find spheres on facet
> O.step();
>
> for i in O.interactions:
> ##if not i.isReal : continue
> ### Rk: facet are only stored in id1
> if isinstance(O.bodies[i.id1].shape,Facet) and
> isinstance(O.bodies[i.id2].shape,Sphere):
> vertices=O.bodies[i.id1].shape.vertices
> normalRef=vertices[0].cross(vertices[1]) # defines the normal to
> the facet normalRef
> nRef=normalRef/(normalRef.norm()) ## normalizes normalRef
> normalFacetSphere=i.geom.normal # geom.normal is oriented from id1
> to id2 -> normalFacetSphere from facet (i.id1) to sphere (i.id2)
>
> if O.bodies[i.id2].state.onJoint==False : ## particles has not yet
> been identified as belonging to a joint plane
> O.bodies[i.id2].state.onJoint=True
> O.bodies[i.id2].state.joint=1
> O.bodies[i.id2].shape.color=jointcolor1
> if nRef.dot(normalFacetSphere)>=0 :
> O.bodies[i.id2].state.jointNormal1=nRef
> elif nRef.dot(normalFacetSphere)<0 :
> O.bodies[i.id2].state.jointNormal1=-nRef
> elif O.bodies[i.id2].state.onJoint==True : ## particles has
> already been identified as belonging to, at least, 1 facet
> if O.bodies[i.id2].state.joint==1 and
> ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05) : ##
> particles has already been identified as belonging to only 1 facet
> O.bodies[i.id2].state.joint=2
> O.bodies[i.id2].shape.color=jointcolor2
> if nRef.dot(normalFacetSphere)>=0 :
> O.bodies[i.id2].state.jointNormal2=nRef
> elif nRef.dot(normalFacetSphere)<0 :
> O.bodies[i.id2].state.jointNormal2=-nRef
> elif O.bodies[i.id2].state.joint==2 and
> ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05) and
> ((O.bodies[i.id2].state.jointNormal2.cross(nRef)).norm()>0.05): ##
> particles has already been identified as belonging to more than 1 facet
> O.bodies[i.id2].state.joint=3
> O.bodies[i.id2].shape.color=jointcolor3
> if nRef.dot(normalFacetSphere)>=0 :
> O.bodies[i.id2].state.jointNormal3=nRef
> elif nRef.dot(normalFacetSphere)<0 :
> O.bodies[i.id2].state.jointNormal3=-nRef
> elif O.bodies[i.id2].state.joint==3 and
> ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05) and
> ((O.bodies[i.id2].state.jointNormal2.cross(nRef)).norm()>0.05) and
> ((O.bodies[i.id2].state.jointNormal3.cross(nRef)).norm()>0.05):
> O.bodies[i.id2].state.joint=4
> O.bodies[i.id2].shape.color=jointcolor5
>
> #### second step -> find spheres interacting with spheres on facet (could
> be executed in the same timestep as step 1?)
> for j in O.interactions:
> #if not i.isReal : continue
> ## Rk: facet are only stored in id1
> if isinstance(O.bodies[j.id1].shape,Facet) and
> isinstance(O.bodies[j.id2].shape,Sphere):
> vertices=O.bodies[j.id1].shape.vertices
> normalRef=vertices[0].cross(vertices[1]) # defines the normal to
> the facet normalRef
> nRef=normalRef/(normalRef.norm()) ## normalizes normalRef
> if ((O.bodies[j.id2].state.jointNormal1.cross(nRef)).norm()<0.05) :
> jointNormalRef=O.bodies[j.id2].state.jointNormal1
> elif
> ((O.bodies[j.id2].state.jointNormal2.cross(nRef)).norm()<0.05) :
> jointNormalRef=O.bodies[j.id2].state.jointNormal2
> elif
> ((O.bodies[j.id2].state.jointNormal3.cross(nRef)).norm()<0.05) :
> jointNormalRef=O.bodies[j.id2].state.jointNormal3
> else : continue
> facetCenter=O.bodies[j.id1].state.pos
> #### seek for each sphere interacting with the identified sphere
> j.id2
> for n in O.interactions.withBody(j.id2) :
> if isinstance(O.bodies[n.id1].shape,Sphere) and
> isinstance(O.bodies[n.id2].shape,Sphere):
> if j.id2==n.id1: # the sphere that was detected on facet
> (that is, j.id2) is id1 of interaction n
> sphOnF=n.id1
> othSph=n.id2
> elif j.id2==n.id2: # here, this sphere that was detected
> on facet (that is, j.id2) is id2 of interaction n
> sphOnF=n.id2
> othSph=n.id1
> facetSphereDir=(O.bodies[othSph].state.pos-facetCenter)
> if O.bodies[othSph].state.onJoint==True :
> if O.bodies[othSph].state.joint==3 and
> ((O.bodies[othSph].state.jointNormal1.cross(jointNormalRef)).norm()>0.05)
> and
> ((O.bodies[othSph].state.jointNormal2.cross(jointNormalRef)).norm()>0.05)
> and
> ((O.bodies[othSph].state.jointNormal3.cross(jointNormalRef)).norm()>0.05):
> O.bodies[othSph].state.joint=4
> O.bodies[othSph].shape.color=jointcolor5
> elif O.bodies[othSph].state.joint==2 and
> ((O.bodies[othSph].state.jointNormal1.cross(jointNormalRef)).norm()>0.05)
> and
> ((O.bodies[othSph].state.jointNormal2.cross(jointNormalRef)).norm()>0.05):
> O.bodies[othSph].state.joint=3
> if facetSphereDir.dot(jointNormalRef)>=0:
>
> O.bodies[othSph].state.jointNormal3=jointNormalRef
> elif facetSphereDir.dot(jointNormalRef)<0:
>
> O.bodies[othSph].state.jointNormal3=-jointNormalRef
> elif O.bodies[othSph].state.joint==1 and
> ((O.bodies[othSph].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) :
> O.bodies[othSph].state.joint=2
> if facetSphereDir.dot(jointNormalRef)>=0:
>
> O.bodies[othSph].state.jointNormal2=jointNormalRef
> elif facetSphereDir.dot(jointNormalRef)<0:
>
> O.bodies[othSph].state.jointNormal2=-jointNormalRef
> elif O.bodies[othSph].state.onJoint==False :
> O.bodies[othSph].state.onJoint=True
> O.bodies[othSph].state.joint=1
> O.bodies[othSph].shape.color=jointcolor4
> if facetSphereDir.dot(jointNormalRef)>=0:
> O.bodies[othSph].state.jointNormal1=jointNormalRef
> elif facetSphereDir.dot(jointNormalRef)<0:
> O.bodies[othSph].state.jointNormal1=-jointNormalRef
>
>
> O.resetTime()
> O.interactions.clear()
> print '\n IdentificationSpheresOnJoint executed ! Spheres onJoint (and so
> on...) detected, facets deleted\n\n'
>
>
>
>
>
>
> --
> 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@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>
--
You received this question notification because your team yade-users is
an answer contact for Yade.
