yade-users team mailing list archive

Thread
Date
Re: [Question #235138]: how to define friction at the boundaries in triaxial test

To: yade-users@xxxxxxxxxxxxxxxxxxx
From: stephen <question235138@xxxxxxxxxxxxxxxxxxxxx>
Date: Thu, 05 Sep 2013 10:31:35 -0000
Reply-to: question235138@xxxxxxxxxxxxxxxxxxxxx
Sender: bounces@xxxxxxxxxxxxx
Question #235138 on Yade changed:
https://answers.launchpad.net/yade/+question/235138

    Status: Answered => Open

stephen is still having a problem:
hi christian,
thanks for your observation,
i have made the necessary corrections and get the best output that correspond to exertly what i was expecting.
how ever, i have some additional development during the simulation. an error named " FATAL /build/buildd/yade-0.70.0/core/ThreadRunner.cpp:31 run: Exception occured: 
Body #37 has velocity==NaN! "
occured and this automatically terminates the simulation. 

do you know anything that could be the cause of this problem and how to fix it?.
i have checked the similar post here but the problem still persist.

about setting the material for sphere, i used this command to do that

 O.materials[idSand].young=1e6.

thanks for your great assistance....am greatful

best regards
stephen

you may also wish to go through the recent script shown below

#file=open("data1.txt","w")
#file.close()

#stra=open("strain.txt","w")
#stra.write(str(0)+"\n") #because we need to add the zero strain. #look at later, may be we should also add zero force???
#stra.close()

devi=open("deviatoric.txt","w") #deviatoric stresses file
#devi.write(str(0)+"\n") #adding zero deviatoric stress
devi.close()

strain=open("strain.txt","w") 
#strain.write(str(0)+"\n") #adding zero strain
strain.close() 

sigma11=open("sigma11.txt","w") # Axial stress file
#strain.write(str(0)+"\n") #adding zero strain
sigma11.close() 

shear11=open("shear11.txt","w")# shear stress file. shear stress = ((sigma1-sigma2)/2)
#strain.write(str(0)+"\n") #adding zero strain
shear11.close() 

epsellion=("epsellion.txt","w")
#epsellion.close()

file=open("sigma1.txt","w")
file1=open("sigma2.txt","w")
file.close()
file1.close()

hop=open("leng.txt","w")
#strain.write(str(0)+"\n") #adding zero strain
hop.close() 


from yade import pack, plot

young_modulus=1e6 #young modulus of the spheres [Pa] (1e6 kPa=1e9 Pa)
poison_ratio=0.25 #poisson's ratios of the spheres
fa=0.5 #friction angle of the spheres [rad]  0.2 radian=11.46 degree, 0.5 radian= 28.6 degree, 0.4 rad=22.9 degreee, 0.8 rad=45.8 degree
fb=0.5 #friction at the boundaries in 
number=1400


##create materials for spheres and plates

idSand=O.materials.append(FrictMat(young=young_modulus,poisson=poison_ratio,frictionAngle=fa,label='idSand'))
#defining sphere properties

#SphereMat =
FrictMat(young=young_modulus,poisson=poisson_ratio,density=density,frictionAngle=0.5)

BoundaryMat =
FrictMat(young=young_modulus,poisson=poison_ratio,frictionAngle=fb)

O.materials.append(BoundaryMat)

O.bodies.append(geom.facetBox((0.035,0.035,0.035),(0.035,0.035,0.035),wallMask=15, material = BoundaryMat, dynamic=1))
O.bodies.append(utils.geom.facetBox((0.035,0.035,0.035),(0.035,0.035,0.035),wallMask=16, material = BoundaryMat,dynamic=0))#bottom of the sample, which is not allowed to move. because of that dynamic=0 [O means "False"]


O.bodies.append(utils.geom.facetBox((0.035,0.035,0.105),(0.035,0.035,0.035),wallMask=15, material = BoundaryMat))#!!!! in order to control the particles!!! after the end of the stabilization, we will remove it


psdSizes,psdCumm=[0.00525,0.00603,0.0105,0.01206,0.02099],[0,0.1,0.5,0.6,1.0] #d50=10.5mm, cu=2


mincorner,maxcorner=Vector3(0.001,0.001,0.001),Vector3(0.07,0.07,0.14)  #min and max coordinates of the wall


sp=pack.SpherePack()   #burda pack modulunun Spherepack modulunu cagiriyo, bi altinda da o cagirdinin,piramit gibi lan iste!
sp.makeCloud(mincorner,maxcorner,psdSizes=psdSizes,num=number,psdCumm=psdCumm,distributeMass=1)  
#sp.makeCloud(mincorner,maxcorner,rMean=0.005,rRelFuzz=0.3)
Gl1_Sphere.stripes=True  #showing the spheres in two colors

sp.toSimulation()#color=(0,0,1))

unb=utils.unbalancedForce()

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]), 
   InteractionLoop(                                             
      
      [Ig2_Sphere_Sphere_L3Geom(),Ig2_Facet_Sphere_L3Geom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_L3Geom_FrictPhys_ElPerfPl()]
   ),
   GravityEngine(gravity=(0,0,-9.81)),
   NewtonIntegrator(damping=0.8), #after the first phase we are going to change it to 0.5     
   PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
     
   
]
   
O.dt=0.3*utils.PWaveTimeStep()

#### D E P O S I T I O N ####

print '### Deposion is Started ###'


def checkUnbalanced():

  print 'iteration number:',O.iter,'unbalanced force:', unbalancedForce()
  if O.iter<150000: return 
   # the rest will be run only if unbalanced is < .1 (stabilized packing)
  if utils.unbalancedForce()>0.1: return
   
  if utils.unbalancedForce()<0.1: O.bodies.append(geom.facetBox((0.035,0.035,0.035),(0.035,0.035,0.035),wallMask=32,material = BoundaryMat, dynamic=1))
 
     
  

  checker.command='erasingwalls()'
  print "Depostion is Completed and the Top Part of the Membran is Added!"

def erasingwalls():
   print "### Additional Walls are Erased! ###"
   O.bodies.erase(10)
   O.bodies.erase(11)
   O.bodies.erase(12)
   O.bodies.erase(13)
   O.bodies.erase(14)
   O.bodies.erase(15)
   O.bodies.erase(16)
   O.bodies.erase(17)
   
   checker.command='load()'


print "loading is starting!"      


aim=100000  #PASCAL!!

vel=0.003     #applyied velocity
#0.003 calisii


#### C O N F I N I N G   P R E S S U R E ####
 
def load():   
  
  O.materials[idSand].young=1e6
  O.engines=O.engines+[PyRunner(command='writingDatas()',iterPeriod=100000)] 
  print ('CONFINING PRESSURE')
  print 'stress at the top :',(abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))
  
  print 'stress at the front:',(abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))
   

  print 'stress at the back:',(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))
  
  print 'stress at the right:',(abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))
  
  print 'stress at the left:',(abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))  
  

    
  if (abs((O.forces.f(O.bodies[number+18].id)[2])+(O.forces.f(O.bodies[number+19].id)[2])))/((abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1]))))<aim:
      
      O.bodies[(number+18)].state.vel=(0,0,(-vel)) #top
      O.bodies[(number+19)].state.vel=(0,0,(-vel)) 
 
  if (abs((O.forces.f(O.bodies[number+18].id)[2])+(O.forces.f(O.bodies[number+19].id)[2])))/((abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1]))))>aim:
      O.bodies[(number+18)].state.vel=(0,0,vel/2) #top
      O.bodies[(number+19)].state.vel=(0,0,vel/2)


  if (abs((O.forces.f(O.bodies[2].id)[0])+(O.forces.f(O.bodies[3].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))<aim:
      O.bodies[2].state.vel=(-vel,0,0) #front
      O.bodies[3].state.vel=(-vel,0,0)
  if (abs((O.forces.f(O.bodies[2].id)[0])+(O.forces.f(O.bodies[3].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))>aim:
      O.bodies[2].state.vel=(vel/2,0,0) #front
      O.bodies[3].state.vel=(vel/2,0,0)


  if (abs((O.forces.f(O.bodies[0].id)[0])+(O.forces.f(O.bodies[1].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))<aim:
      O.bodies[0].state.vel=(vel,0,0) #behind
      O.bodies[1].state.vel=(vel,0,0)
  if (abs((O.forces.f(O.bodies[0].id)[0])+(O.forces.f(O.bodies[1].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))>aim:
      O.bodies[0].state.vel=(-vel/2,0,0) #behind
      O.bodies[1].state.vel=(-vel/2,0,0)


  if (abs((O.forces.f(O.bodies[6].id)[1])+(O.forces.f(O.bodies[7].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
      O.bodies[6].state.vel=(0,-vel,0) #right
      O.bodies[7].state.vel=(0,-vel,0)
  if (abs((O.forces.f(O.bodies[6].id)[1])+(O.forces.f(O.bodies[7].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
      O.bodies[6].state.vel=(0,(vel/2),0) #right
      O.bodies[7].state.vel=(0,(vel/2),0)


  if (abs((O.forces.f(O.bodies[4].id)[1])+(O.forces.f(O.bodies[5].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
      O.bodies[4].state.vel=(0,vel,0) #left
      O.bodies[5].state.vel=(0,vel,0)    
  if (abs((O.forces.f(O.bodies[4].id)[1])+(O.forces.f(O.bodies[5].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
      O.bodies[4].state.vel=(0,-vel/2,0) #left
      O.bodies[5].state.vel=(0,-vel/2,0)

      
  if (abs((O.forces.f(O.bodies[number+18].id)[2])+(O.forces.f(O.bodies[number+19].id)[2])))/((abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1]))))>aim and (abs((O.forces.f(O.bodies[2].id)[0])+(O.forces.f(O.bodies[3].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))>aim and (abs((O.forces.f(O.bodies[0].id)[0])+(O.forces.f(O.bodies[1].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))>aim and (abs((O.forces.f(O.bodies[6].id)[1])+(O.forces.f(O.bodies[7].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim and (abs((O.forces.f(O.bodies[4].id)[1])+(O.forces.f(O.bodies[5].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
     
     print "loading is stopped"
     uzun=(abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))
     print uzun
     hip=open("leng.txt","a")
     hip.write(str(uzun)+"\n")
     hip.close() 
     checker.command='axial()'
     print "### Axial Loading is Starting ###"
    


#### A X I A L   L O A D I N G ####

def axial():
  O.materials[idSand].young=1e6
  #O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=200)]

  
  O.engines=O.engines+[PyRunner(command='writingLength()',iterPeriod=150000)] 


  print 'stress at the top :',(abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1]))),'iteration number:',O.iter
  
  print 'stress at the front:',(abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))
   

  print 'stress at the back:',(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))
  
  print 'stress at the right:',(abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))
  
  print 'stress at the left:',(abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))  
      
  

  if (abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))<4*aim:
      O.bodies[number+18].state.vel=(0,0,-2*(vel)) #top
      O.bodies[number+19].state.vel=(0,0,-2*(vel)) 
  

  if (abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
      O.bodies[2].state.vel=((-vel/4),0,0) #front
      O.bodies[3].state.vel=((-vel/4),0,0)

  if (abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
      O.bodies[2].state.vel=((vel/4),0,0) #front
      O.bodies[3].state.vel=((vel/4),0,0)
      O.bodies[number+18].state.vel=(0,0,-(vel)) #top
      O.bodies[number+19].state.vel=(0,0,-(vel))
  


  if (abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
      O.bodies[0].state.vel=((vel/4),0,0) #behind
      O.bodies[1].state.vel=((vel/4),0,0)

  if (abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
      O.bodies[0].state.vel=((-vel/4),0,0) #behind
      O.bodies[1].state.vel=((-vel/4),0,0)
      O.bodies[number+18].state.vel=(0,0,-(vel)) #top
      O.bodies[number+19].state.vel=(0,0,-(vel))
 


  if (abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
      O.bodies[6].state.vel=(0,(-vel/4),0) #right
      O.bodies[7].state.vel=(0,(-vel/4),0)
  
  if (abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
      O.bodies[6].state.vel=(0,(vel/4),0) #right
      O.bodies[7].state.vel=(0,(vel/4),0)
      O.bodies[number+18].state.vel=(0,0,-(vel)) #top
      O.bodies[number+19].state.vel=(0,0,-(vel))
  


  if (abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
      O.bodies[4].state.vel=(0,(vel/4),0) #left
      O.bodies[5].state.vel=(0,(vel/4),0)    
  
  if (abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
      O.bodies[4].state.vel=(0,(-vel/4),0) #left
      O.bodies[5].state.vel=(0,(-vel/4),0)
      O.bodies[number+18].state.vel=(0,0,-(vel)) #top
      O.bodies[number+19].state.vel=(0,0,-(vel))
  

      
  if (abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))>4*aim: 
     
     print "### Triaxial Test is Completed ###"
     print "###     ENDE GUT ALLES GUT     ###"
     print "###   Thank you JAMAL          ###"

     checker.command='strain()'
     
       
      
 
def writingDatas():
  
  

x=(abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))

y=(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))


  file=open("sigma1.txt","a")
  file1=open("sigma2.txt","a")

  file.write(str(x)+"\n")
  file1.write(str(y)+"\n")
  #+"\n")
  file.close()
  file1.close()

for b in O.bodies: b.shape.color=utils.scalarOnColorScale(b.state.rot().norm(),0,pi/2.)
  # rot() gives rotation vector between reference and current position
    
  
  

def writingLength():
  
  leng=open("leng.txt","r")
  uzunluk=leng.readlines()
  
  first=(float(uzunluk[0]))
  leng.close()
  
    
 
  if (abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))<4*aim:

s=((0.07)-abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))/(0.07)

E=((0.07)-abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1])))/(0.07)


    shr=((((abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))-(((abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))/4)))/1000)/2
    sig=((((abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))+(((abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))/4)))/1000)
    


    y=(((abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))-(((abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))/4)))/1000
    
    devi=open("deviatoric.txt","a")
    devi.write(str(y)+"\n")
    devi.close()

    strain=open("strain.txt","a")
    strain.write(str(s)+"\n")
    strain.close() 
    
    epsellion=open("epsellion.txt","a")
    epsellion.write(str(E)+"\n")
    epsellion.close()
    

    sigma11=open("sigma11.txt","a")
    sigma11.write(str(sig)+"\n")
    sigma11.close()
    
    shear11=open("shear11.txt","a")
    shear11.write(str(shr)+"\n")
    shear11.close()
    
    
    
    print ("strain:",s,"  ","epsellion:",E,"   ","deviatoric:",y) 


def strain():
  
 # length=open("length.txt","r")
  #numbers=length.readlines()
  
  
  #for i in range(0,(len(numbers)-1)):
   #"stra" defined as the file of strains calculated for each iteration.the calculation starts with the end of axial loading.
   # stra=open("strain.txt","a")
    #strain=(abs(float(numbers[i]))-(float(numbers[i+1])))/(float(numbers[i]))
    #stra.write(str(strain)+"\n")
    #stra.close()   	
 
  #print ("### calculation of the strain completed ###")
  print ("_____________________________________________________________________")
  O.pause()   #this stops yade!
   
O.run()

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.