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'

When	Who	Action	New State
2016-03-31	Jan Stránský	Answer	Answered
2016-03-31	Yor1	Comment	Answered
2016-04-12	Jérôme Duriez	Answer	Answered
2016-05-02	Yor1	Confirm	Solved

Yade

History of question #289662