Hello everyone,
I try to model the 3-points bending test with several number of particle in a beam 240x80x30 mm.
To calculate the elastic energy i use this formula:
E+=0.5*(i.phys.normalForce.squaredNorm()/i.phys.kn + i.phys.shearForce.squaredNorm()/i.phys.ks)
But when i execute my code i have this error message: "ZeroDivisionError: float division by zero".
This is my code:
from yade import ymport, utils, pack, export
from yade import plot
from pylab import *
import math
PACKING='X240Y80Z30_2k'
OUT=PACKING+'_flexion_3_points'
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
intR=1.54726
DENS=4000
YOUNG=65e9
FRICT=10
ALPHA=0.4
TENS=8e6
COH=160e6
sphereMat = JCFpmMat(type=1,density=DENS,young=YOUNG,poisson = ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
wallMat = JCFpmMat(type=0,density=DENS,young=YOUNG,frictionAngle=radians(0))
for mat in (sphereMat,wallMat):
O.materials.append(mat) # then wallMat will be used if material is not specified
O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0)))
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=X/15.
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()
piston3 = O.bodies.append(geom.facetCylinder(center=(xinf+X/12.,yinf-r+0.1*Rmean,Z/2.),radius=r,height=Z,orientation=Quaternion((1, 0, 0), 0),segmentsNumber=20,wire=False,material=wallMat)) # bas gauche
piston2 = O.bodies.append(geom.facetCylinder(center=(xinf+11*X/12.,yinf-r+0.1*Rmean,Z/2.),radius=r,height=Z,orientation=Quaternion((1, 0, 0), 0),segmentsNumber=20,wire=False,material=wallMat)) # bas droite
piston1 = O.bodies.append(geom.facetCylinder(center=(0.5*X,Y+0.9*r,Z/2.),radius=r,height=Z,dynamic=False,orientation=Quaternion((1, 0, 0), 0),segmentsNumber=20,wire=False,material=wallMat)) # haut
p0=O.bodies[piston1[0]].state.pos[1]
for i in range(0,len(piston1)):
O.bodies[piston1[i]].state.vel[1]=-0.001
beam=O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))
import gts
alpha = 90 #angle of the crack with horizontal (pi/4.)
c = 0.375*Y # length of the crack
ptA = gts.Vertex( X/2., c, 8./7.*Z)
ptB = gts.Vertex( X/2., 0, 8./7.*Z)
ptApr = gts.Vertex(X/2., c, -1./7.*Z)
ptBpr = gts.Vertex(X/2., 0, -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')
for o in O.bodies:
if isinstance(o.shape,Sphere):
o.shape.color=(0.7,0.5,0.3)
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb'),Bo1_Facet_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Facet_Sphere_ScGeom()],
[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.4, defaultDt=0.1*utils.PWaveTimeStep()),
NewtonIntegrator(gravity=(0,-9.81,0),damping=DAMP,label="newton"),
PyRunner(iterPeriod=int(saveData),initRun=True,command='recorder()',label='data'),
VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['facets','spheres','bstresses','jcfpm','cracks'],Key=OUT,label='vtk')
]
g=[0,-9.81,0]
tensCks=shearCks=cks=cks0=0
def recorder():
Ub=0
E=0
global tensCks, shearCks
tensCks=0
shearCks=0
for o in O.bodies:
if isinstance(o.shape,Sphere):
tensCks+=o.state.tensBreak
shearCks+=o.state.shearBreak
Ub+=sum([b.state.mass*abs(g[1])*abs(b.state.displ()[1]) for b in O.bodies])
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):
E+=0.5*(i.phys.normalForce.squaredNorm()/i.phys.kn + i.phys.shearForce.squaredNorm()/i.phys.ks)
yade.plot.addData({ 't':O.time
,'i':O.iter
,'f1':utils.sumForces(piston1,Vector3(0,1,0))
,'f2':utils.sumForces(piston2,Vector3(0,-1,0))
,'f3':utils.sumForces(piston3,Vector3(0,-1,0))
,'displacement':p0-O.bodies[piston1[0]].state.pos[1]
,'E':E
,'Wp':utils.sumForces(piston1,Vector3(0,1,0))*(p0-O.bodies[piston1[0]].state.pos[1])
,'Ub':Ub
,'tc':0.5*tensCks,'sc':0.5*shearCks,'unbF':utils.unbalancedForce()
}
)
yade.plot.saveDataTxt(OUT)
And this is the code "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)
O.dt=0.001*utils.PWaveTimeStep()
############################ 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
for b in O.bodies:
if isinstance(b.shape,Facet):
O.bodies.erase(b.id)
O.resetTime()
O.interactions.clear()
print '\n IdentificationSpheresOnJoint executed ! Spheres onJoint (and so on...) detected, facets deleted\n\n'
