A bug for using of Pyrunner[]

Hi, can you reproduce the "bug" with a MWE [1]? You can't expect us to look
at 800 lines of code. Also, I can run the script so there is no bug. Try to
be more specific when asking questions.

And finally, looks like you are running a triaxial test. Why not using
TriaxialStressController [2]? For an example see [3].

HTH Klaus

[1] https://yade-dem.org/wiki/Howtoask
[2]
https://yade-dem.org/doc/yade.wrapper.html?highlight=triaxialstresscontroller#yade.wrapper.TriaxialStressController
[3] https://github.com/yade/trunk/tree/master/examples/triax-tutorial

On Sun, May 20, 2018 at 1:50 PM, De zhang <
<email address hidden>> wrote:

> New question #669404 on Yade:
> https://answers.launchpad.net/yade/+question/669404
>
> Hello,
> Following the last question, I found that a bug for using PyRunner[].
> I found that at the begining of O.engines=[...PyRunner[1...],PyRunner[2...,label=2.],],
> the PyRunner[2...,label=2.] was called in the process of PyRunner[1...] by
> using 'label[2].dead=False', the the PyRunner[2...] works continuously,
> but if the PyRunner[2...] was called in the process of PyRunner[1...] by
> using the O.engines=O.engines+ PyRunner[2...], the PyRunner[2...] works
> discontinuously.
> is it a bug for using the PyRunner[]?
> The following is two use of PyRunner.
> ###########################################
> ####1st using of '----O.engines=O.engines+PyRunner[]'---###
> #This simulation for triaxial experiment of ballast which size betweeen
> 30cm~45cm
> #Friction angle for 48 degree
> from yade import pack,qt,plot,utils,polyhedra_utils,ymport,export,pack,
> timing
> from yade import *
> import numpy
> from pprint import pprint
> import random
> from random import uniform
> #from random import randint
> import math
> from math import *
> global gravel,steel
> gravel = FrictMat()
> gravel.density = 2600 #kg/m^3
> gravel.young = 2e9
> gravel.poisson = 0.21 # real 0.21
> gravel.frictionAngle = 0.83 #rad radians(48) // change for rad
> math.radians(31)
> steel = FrictMat()
> steel.density = 7850 #kg/m^3
> steel.young = 10*gravel.young #inital steel was 10*gravel.young
> steel.poisson = 0.3
> steel.frictionAngle = 0.55 #rad radians(31)
> ##
> bottom_wall=utils.wall(0.00,axis=2,sense=1,material=steel)
> O.bodies.append(bottom_wall)
> bottom_wall.state.blockedDOFs='xyzXYZ'
> ###Number for 7 walls
> for i in range(1,8):
> O.bodies.append(utils.facet(vertices=((-0.3,-0.15,0.1*(i-
> 1)),(0.3,-0.15,0.1*(i-1)),(0.3,-0.15,0.1*i)),dynamic=None,
> fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
> O.bodies.append(utils.facet(vertices=((-0.3,-0.15,0.1*(i-
> 1)),(-0.3,-0.15,0.1*i),(0.3,-0.15,0.1*i)),dynamic=None,
> fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
>
> for i in range(1,8):
> O.bodies.append(utils.facet(vertices=((-0.3,0.15,0.1*(i-1)
> ),(0.3,0.15,0.1*(i-1)),(0.3,0.15,0.1*i)),dynamic=None,fixed=
> True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
> O.bodies.append(utils.facet(vertices=((-0.3,0.15,0.1*(i-1)
> ),(-0.3,0.15,0.1*i),(0.3,0.15,0.1*i)),dynamic=None,fixed=
> True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
>
> for i in range(1,8):
> O.bodies.append(utils.facet(vertices=((-0.15,-0.3,0.1*(i-
> 1)),(-0.15,0.3,0.1*(i-1)),(-0.15,0.3,0.1*i)),dynamic=None,
> fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
> O.bodies.append(utils.facet(vertices=((-0.15,-0.3,0.1*(i-
> 1)),(-0.15,-0.3,0.1*i),(-0.15,0.3,0.1*i)),dynamic=None,
> fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
>
> for i in range(1,8):
> O.bodies.append(utils.facet(vertices=((0.15,-0.3,0.1*(i-1)
> ),(0.15,0.3,0.1*(i-1)),(0.15,0.3,0.1*i)),dynamic=None,
> fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
> O.bodies.append(utils.facet(vertices=((0.15,-0.3,0.1*(i-1)
> ),(0.15,-0.3,0.1*i),(0.15,0.3,0.1*i)),dynamic=None,fixed=
> True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
>
> global Wall1Stressx,Wall2Stressx,Wall3Stressx,Wall4Stressx
> global Wall1Stressy,Wall2Stressy,Wall3Stressy,Wall4Stressy
> global Wall1S,Wall2S,Wall3S,Wall4S
> Wall1Stressx=0
> Wall2Stressx=0
> Wall3Stressx=0
> Wall4Stressx=0
> Wall1Stressy=0
> Wall2Stressy=0
> Wall3Stressy=0
> Wall4Stressy=0
> Wall1S=0
> Wall2S=0
> Wall3S=0
> Wall4S=0
> global WallStress,ConfStress,ConfDevi,MoveVel,AxiDevi,MoveAxial #stress
> control
> # Area of the confining Wall
> global A1,A2,A3,A4
> global LoadPos,IniLoadPos,plateF,IniTime,forceA
> global WallStress,ConfStress,ConfDevi,MoveVel,AxiDevi,MoveAxial,AreaPlate
> #unit:m^2
> IniTime=0
> plateF=0 #Unit:kPa
> LoadPos=0.6
> IniLoadPos=LoadPos # (link to Area of
> Walls)
> forceA=200 # Unit:kPa,P=N/A;N=P*0.0615*1000;A=0.0615
> A1=LoadPos*0.3
> A2=LoadPos*0.3
> A3=LoadPos*0.3
> A4=LoadPos*0.3
> WallStress=0 # Unit:kPa
> ConfStress=100 # Unit:kPa
> ConfDevi=0
> AxiDevi=0
> MoveVel=0
> MoveAxial=0
> AreaPlate=0.09
> # Id of different substances
> global NumLoad,NumEndBall,StepNum,NumEnd,xratio,yratio,zratio,NumContact,WallContact
>
> NumLoad=1
> NumEndBall=1
> StepNum=1
> NumEnd=1
> xratio=1
> yratio=1
> zratio=1
> NumContact=4
> WallContact=1
> # Position and other parameters record
> ######################parameters
> sp=pack.SpherePack()
> sp.makeCloud((-0.15,-0.15,0),(0.15,0.15,0.8),rMean=0.016,rRelFuzz=0.25)
> sp.toSimulation(material=gravel)
> NumEndBall=O.bodies[-1].id#Mark Sphere
> global iternum
> iternum=0
> #O.dt=1.0e-6 #Check it!
> O.dt=8.0e-6 #Check it!
>
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_
> Aabb(),Bo1_Facet_Aabb(),Bo1_Sphere_Aabb()]),
> InteractionLoop(
> [Ig2_Wall_Polyhedra_PolyhedraGeom(),Ig2_Polyhedra_
> Polyhedra_PolyhedraGeom(),Ig2_Facet_Polyhedra_PolyhedraGeom(
> ),Ig2_Sphere_Sphere_ScGeom(),Ig2_Sphere_Polyhedra_ScGeom(),
> Ig2_Facet_Sphere_L3Geom(),Ig2_Wall_Sphere_ScGeom()],
> [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys(),
> Ip2_FrictMat_PolyhedraMat_FrictPhys(),Ip2_FrictMat_FrictMat_FrictPhys()],
> [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric(),
> Law2_L3Geom_FrictPhys_ElPerfPl(),Law2_ScGeom_FrictPhys_CundallStrack(),
> Law2_L6Geom_FrictPhys_Linear()],
> ),
> NewtonIntegrator(damping=0.6,gravity=(0,0,-9.81),label='newton'),
> PyRunner(command='TraiStep()',iterPeriod=1,label='checker'),
> PyRunner(command='LoadAxial100kPa()',iterPeriod=
> 1,label='loadkeep100kPa'),
> ]
>
> ##Fullfill the box
> def TraiStep():
> global Wall1Stressx,Wall2Stressx,Wall3Stressx,Wall4Stressx
> global Wall1Stressy,Wall2Stressy,Wall3Stressy,Wall4Stressy
> global Wall1S,Wall2S,Wall3S,Wall4S
> global WallStress,ConfStress,ConfDevi,MoveVel,AxiDevi,MoveAxial
> #stress control,WallStress,ConfStress,ConfDevi,MoveVel
> global A1,A2,A3,A4
> global LoadPos,NumLoad,NumEndBall,IniLoadPos,plateF,IniTime,
> forceA,StepNum,NumEnd,iternum,AreaPlate
> global xratio,yratio,ztario,NumContact,WallContact
> ######
> #Step1=> add the loadingplate
> #Step2=> apply the initial axial force and confing force
> #Step3=> apply the loading force and confining stress
> if StepNum == 1:
> loadkeep100kPa.dead=True
> #loadkeep200kPa.dead=True
> StepNum=StepNum+1
> elif StepNum == 2:
> print "2-unbalanced forces = %.5f"%(utils.unbalancedForce()
> )
> if O.iter < 30000: return
> if utils.unbalancedForce() > 0.01: return
> iternum=O.iter
> m=max([b.state.pos[2]+b.shape.radius for b in O.bodies if
> isinstance(b.shape,Sphere)])
> O.bodies.append(utils.wall(m,
> axis=2,sense=0,material=steel))
> NumLoad=O.bodies[-1].id
> NumEnd=O.bodies[-1].id
> LoadPos=O.bodies[NumLoad].state.pos[2]
> StepNum=StepNum+1
> elif StepNum == 3:
> LoadPos=O.bodies[NumLoad].state.pos[2]
> print "3-Loadplate force= %.5f"%(plateF)
> AreaPlate=(O.bodies[50].state.
> pos[0]-O.bodies[36].state.pos[0])*(O.bodies[50].state.pos[0]
> -O.bodies[36].state.pos[0])
> plateF=(O.forces.f(NumLoad)[2])/(AreaPlate*1000)
> #P=F/A=F/(0.0615*1000)=F/61.5 Unit:kPa
> LoadPos=O.bodies[NumLoad].state.pos[2]
> if plateF < 100:
> O.bodies[NumLoad].state.vel=(0,0,-0.005) #100mm/s
> O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
> else:
> O.bodies[NumLoad].state.vel=(0,0,0)
> O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
> StepNum=StepNum+1
> #O.pause(),first I got to the 200kPa axial stress, then
> keep loading axial stress
> elif StepNum == 4:
> loadkeep100kPa.dead=False
> StepNum=StepNum+1
> #O.pause()
> elif StepNum == 5:
> LoadPos=O.bodies[NumLoad].state.pos[2]
> A1=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].
> state.pos[0])
> A2=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].
> state.pos[0])
> A3=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].
> state.pos[1])
> A4=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].
> state.pos[1])
> for i in range(1,15):
> Wall1Stressx=Wall1Stressx+abs(O.forces.f(i)[0])
> Wall1Stressy=Wall1Stressy+abs(O.forces.f(i)[1])
> Wall1S=Wall1Stressy/A1
> Wall1Stressx=0
> Wall1Stressy=0
> for i in range(15,29):
> Wall2Stressy=Wall2Stressy+abs(O.forces.f(i)[1])
> Wall2S=Wall2Stressy/A2
> Wall2Stressy=0
> for i in range(29,43):
> Wall3Stressx=Wall3Stressx+abs(O.forces.f(i)[0])
> Wall3S=Wall3Stressx/A3
> Wall3Stressx=0
> #Wall3Stressy=0
> for i in range(43,57):
> Wall4Stressx=Wall4Stressx+abs(O.forces.f(i)[0])
> Wall4S=Wall4Stressx/A4
> Wall4Stressx=0
> ##########################
> WallStress=(Wall1S+Wall2S+Wall3S+Wall4S)/4000 ###
> Unit(kPa)
> ConfDevi=(abs(WallStress-ConfStress))/ConfStress ###
> Unit/kPa # parameter!!!
> for i in range(1,57):
> NumContact=NumContact+len(O.
> bodies[NumLoad].intrs())
> WallContact=NumContact/4+1
> NumContact=4
> MoveVel=1*A1*(WallStress-ConfStress)/((1.0e9)*WallContact*(8.0e-6))
>
> ################check the parameter
> # print "Ini-conf-stress= %.5f, Vel= %.8f, WallContact=
> %.1f, NumContact= %.1f, MoveVel= %.8f, Area= %.5f"%(WallStress,
> O.bodies[1].state.vel[1], WallContact, NumContact, MoveVel, A1)
> ################
> #MoveVel=0.000005*(WallStress-ConfStress)
> if abs(MoveVel) > 0.0001:
> MoveVel=0.000001*(WallStress-ConfStress)
> else:
> print "MoveVel is OK"
> for i in range(1,57):
> xratio=(abs(O.bodies[i].state.
> pos[0]))/(O.bodies[i].state.pos[0]+0.00001)
> yratio=(abs(O.bodies[i].state.
> pos[1]))/(O.bodies[i].state.pos[1]+0.00001)
> O.bodies[i].state.vel=(
> MoveVel*xratio,MoveVel*yratio,0)
> for i in range(1,57):
> O.bodies[i].state.blockedDOFs='xyzXYZ'
> print "Ini-Cal-conf-stress= %.5f, Vel= %.8f"%(WallStress,
> O.bodies[1].state.vel[1])
> if ConfDevi > 0.05: return
> for i in range(1,57):
> O.bodies[i].state.blockedDOFs='xyzXYZ'
> print "Ini-Conf-stress= %.5f"%(WallStress)
> StepNum=StepNum+1
> elif StepNum == 6:
> loadkeep100kPa.dead=True
> O.engines=O.engines+[PyRunner(command='Confining()',
> iterPeriod=1)]
> StepNum=StepNum+1
> O.pause()
> elif StepNum == 7:
> AreaPlate=(O.bodies[50].state.
> pos[0]-O.bodies[36].state.pos[0])*(O.bodies[50].state.pos[0]
> -O.bodies[36].state.pos[0])
> plateF=(O.forces.f(NumLoad)[2])/(1000*AreaPlate)
> LoadPos=O.bodies[NumLoad].state.pos[2]
> AxiDevi=abs((plateF-forceA))/forceA
> forceA=200 #10Hz Omega=2*pi*f//
> 120+80*sin(Omega*t)//40-200kPa
> zratio=1*AreaPlate/((2.0e9)*(
> len(O.bodies[NumLoad].intrs())+1)*(8.0e-6)) #alpha=50==>100
> MoveAxial=1*zratio*(plateF-forceA)
> if abs(MoveAxial) > 0.0001:
> MoveAxial=0.000001*(plateF-forceA)
> else:
> print "MoveAxial is OK"
> print "force= %.5f, ForceA= %.5f, Vel= %.8f"%(plateF,
> forceA, MoveAxial)#Unit (kPa)//0.0615 is Area of loadingplate
> O.bodies[NumLoad].state.vel=(0,0,MoveAxial)
> O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
> if AxiDevi > 0.05: return
> print "Loadplate force= %.5f, ForceA= %.5f"%(plateF,
> forceA)
> O.bodies[NumLoad].state.blockedDOFs='xyXYZ'
> IniLoadPos=LoadPos
> IniTime=O.time
> StepNum=StepNum+1
> #O.pause()
> elif StepNum == 8:
> print "8-force= %.5f"%(plateF)
> O.engines=O.engines+[PyRunner(command='AxialLoading()',
> iterPeriod=1)]+[PyRunner(command='addPlotData()',iterPeriod=1)]
> StepNum=StepNum+1
> O.pause()
> else:
> print "Well Done"
> #O.pause()
>
> ##
> def Confining():
> global Wall1Stressx,Wall2Stressx,Wall3Stressx,Wall4Stressx
> global Wall1Stressy,Wall2Stressy,Wall3Stressy,Wall4Stressy
> global Wall1S,Wall2S,Wall3S,Wall4S
> global WallStress,ConfStress,ConfDevi,MoveVel #stress control
> global A1,A2,A3,A4
> global LoadPos,NumLoad,NumEndBall
> global xratio,yratio,NumContact,WallContact #control the velocity
> of confining walls
> LoadPos=O.bodies[NumLoad].state.pos[2]
> A1=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
> A2=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
> A3=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].state.pos[1])
> A4=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].state.pos[1])
> for i in range(1,15):
> Wall1Stressx=Wall1Stressx+abs(O.forces.f(i)[0])
> Wall1Stressy=Wall1Stressy+abs(O.forces.f(i)[1])
> Wall1S=Wall1Stressy/A1
> Wall1Stressx=0
> Wall1Stressy=0
> for i in range(15,29):
> Wall2Stressy=Wall2Stressy+abs(O.forces.f(i)[1])
> Wall2S=Wall2Stressy/A2
> Wall2Stressy=0
> for i in range(29,43):
> Wall3Stressx=Wall3Stressx+abs(O.forces.f(i)[0])
> Wall3S=Wall3Stressx/A3
> Wall3Stressx=0
> #Wall3Stressy=0
> for i in range(43,57):
> Wall4Stressx=Wall4Stressx+abs(O.forces.f(i)[0])
> Wall4S=Wall4Stressx/A4
> Wall4Stressx=0
> ##########################
> WallStress=(Wall1S+Wall2S+Wall3S+Wall4S)/4000 ### Unit(kPa)
> ConfDevi=(abs(WallStress-ConfStress))/ConfStress ### Unit/kPa #
> parameter!!!
> if ConfDevi > 0.05:
> for i in range(1,57):
> NumContact=NumContact+len(O.
> bodies[NumLoad].intrs())
> WallContact=NumContact/4+1
> NumContact=4
> MoveVel=1*A1*(WallStress-ConfStress)/((1.0e9)*
> WallContact*(8.0e-6))
> if abs(MoveVel) > 0.0001:
> MoveVel=0.000001*(WallStress-ConfStress)
> else:
> print "MoveVel is OK"
> for i in range(1,57):
> xratio=(abs(O.bodies[i].state.
> pos[0]))/(O.bodies[i].state.pos[0]+0.00001)
> yratio=(abs(O.bodies[i].state.
> pos[1]))/(O.bodies[i].state.pos[1]+0.00001)
> O.bodies[i].state.vel=(
> MoveVel*xratio,MoveVel*yratio,0)
> for i in range(1,57):
> O.bodies[i].state.blockedDOFs='xyzXYZ'
> print "Keep-Cal-conf-stress= %.5f, Vel= %.8f"%(WallStress,
> MoveVel)
> else:
> for i in range(1,57):
> O.bodies[i].state.blockedDOFs='xyzXYZ'
> print "Keep-Conf= %.5f"%(WallStress)
> ############## keep confining pressure
>
> def LoadAxial100kPa():
> global forceA,plateF,LoadPos,AxiDevi,IniTime,NumLoad,MoveAxial,
> AreaPlate,zratio
> AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[
> 0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
> plateF=(O.forces.f(NumLoad)[2])/(1000*AreaPlate)
> LoadPos=O.bodies[NumLoad].state.pos[2]
> AxiDevi=abs((plateF-forceA))/forceA
> forceA=100 #10Hz Omega=2*pi*f// 120+80*sin(Omega*t)//40-200kPa
> zratio=5*AreaPlate/((2.0e9)*(len(O.bodies[NumLoad].intrs())
> +1)*(8.0e-6))
> MoveAxial=1*zratio*(plateF-forceA)
> if abs(MoveAxial) > 0.0001:
> MoveAxial=0.000001*(plateF-forceA)
> else:
> print "MoveAxial is OK"
> if AxiDevi > 0.05:
> print "keep-force= %.5f, ForceA= %.5f, Vel= %.8f, CNum=
> %.1f"%(plateF, forceA, MoveAxial, len(O.bodies[NumLoad].intrs()))#Unit
> (kPa)//0.0615 is Area of loadingplate
> O.bodies[NumLoad].state.vel=(0,0,MoveAxial)
> O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
> else:
> print "Done-keep-force= %.5f, ForceA= %.5f, Vel= %.8f,
> CNum= %.1f"%(plateF, forceA, MoveAxial, len(O.bodies[NumLoad].intrs()))
>
> O.bodies[NumLoad].state.blockedDOFs='xyXYZ'
>
> ##AxialLoading
> def AxialLoading():
> global forceA,plateF,LoadPos,AxiDevi,IniTime,NumLoad,MoveAxial,
> AreaPlate,zratio
> AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[
> 0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
> plateF=(O.forces.f(NumLoad)[2])/(1000*AreaPlate)
> LoadPos=O.bodies[NumLoad].state.pos[2]
> AxiDevi=abs((plateF-forceA))/forceA
> forceA=200+50*sin((20*pi)*(O.time-IniTime)) #10Hz Omega=2*pi*f//
> 120+80*sin(Omega*t)//40-200kPa
> zratio=0.5*AreaPlate/((2.0e9)*(len(O.bodies[NumLoad].intrs()
> )+1)*(8.0e-6))
> MoveAxial=1*zratio*(plateF-forceA)
> if abs(MoveAxial) > 0.0001:
> MoveAxial=0.000001*(plateF-forceA)
> else:
> print "MoveAxial is OK"
> if AxiDevi > 0.05:
> print "final-force= %.5f, ForceA= %.5f, Vel=
> %.8f"%(plateF, forceA, MoveAxial)#Unit (kPa)//0.0615 is Area of loadingplate
> O.bodies[NumLoad].state.vel=(0,0,MoveAxial)
> O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
> else:
> print "final-Loadplate force= %.5f, ForceA= %.5f"%(plateF,
> forceA)
> O.bodies[NumLoad].state.blockedDOFs='xyXYZ'
>
> ##Record
> def addPlotData():
> global LoadPos,IniLoadPos,NumLoad,forceA,plateF
> global theta,thega,WallStress,Vol,AreaPlate
> theta=forceA
> theta2=plateF
> LoadPos=O.bodies[NumLoad].state.pos[2]
> thega=((IniLoadPos-LoadPos)/IniLoadPos)*100
> AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[
> 0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
> Vol=IniLoadPos*AreaPlate-LoadPos*(O.bodies[1].state.
> pos[1])*(O.bodies[1].state.pos[1])
> plot.addData(Thega=thega,Theta=theta,Thega2=thega,
> Theta2=theta2,T=O.time,Conf=WallStress,TimeLast=O.time,Volume=Vol)
>
> ##ConfiningWall=>wall(1-5)//wall(2-6)//wall(3-7)//wall(4-8)
>
> plot.plots={'Thega':('Theta',),'Thega2':('Theta2',),'T':('
> Conf',),'TimeLast':('Volume',)}
> plot.plot()
>
> qt.Controller()
> V = qt.View()
> V.screenSize = (550,450)
> V.sceneRadius = 1
> V.eyePosition = (0.7,0.5,0.1)
> V.upVector = (0,0,1)
> V.lookAt = (0.15,0.15,0.1)
> #########################################################
>
> #####----2rd is the using of 'label.dead=True'----#################
> #This simulation for triaxial experiment of ballast which size betweeen
> 30cm~45cm
> #Friction angle for 48 degree
> from yade import pack,qt,plot,utils,polyhedra_utils,ymport,export,pack,
> timing
> from yade import *
> import numpy
> from pprint import pprint
> import random
> from random import uniform
> #from random import randint
> import math
> from math import *
>
> ##################################
> #material:ballast and loadingplate
> global gravel
> global steel
>
> gravel = FrictMat()
> gravel.density = 2600 #kg/m^3
> gravel.young = 2e9
> gravel.poisson = 0.21 # real 0.21
> gravel.frictionAngle = 0.83 #rad radians(48) // change for rad
> math.radians(31)
>
>
> steel = FrictMat()
> steel.density = 7850 #kg/m^3
> steel.young = 10*gravel.young #inital steel was 10*gravel.young
> steel.poisson = 0.3
> steel.frictionAngle = 0.55 #rad radians(31)
> ##next
> #################################################################
>
> ##### make circle dormetory
> ### bottom wall
> bottom_wall=utils.wall(0.00,axis=2,sense=1,material=steel)
> O.bodies.append(bottom_wall)
> bottom_wall.state.blockedDOFs='xyzXYZ'
>
>
> ###Number for 7 walls
> for i in range(1,8):
> O.bodies.append(utils.facet(vertices=((-0.3,-0.15,0.1*(i-
> 1)),(0.3,-0.15,0.1*(i-1)),(0.3,-0.15,0.1*i)),dynamic=None,
> fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
> O.bodies.append(utils.facet(vertices=((-0.3,-0.15,0.1*(i-
> 1)),(-0.3,-0.15,0.1*i),(0.3,-0.15,0.1*i)),dynamic=None,
> fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
>
> for i in range(1,8):
> O.bodies.append(utils.facet(vertices=((-0.3,0.15,0.1*(i-1)
> ),(0.3,0.15,0.1*(i-1)),(0.3,0.15,0.1*i)),dynamic=None,fixed=
> True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
> O.bodies.append(utils.facet(vertices=((-0.3,0.15,0.1*(i-1)
> ),(-0.3,0.15,0.1*i),(0.3,0.15,0.1*i)),dynamic=None,fixed=
> True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
>
> for i in range(1,8):
> O.bodies.append(utils.facet(vertices=((-0.15,-0.3,0.1*(i-
> 1)),(-0.15,0.3,0.1*(i-1)),(-0.15,0.3,0.1*i)),dynamic=None,
> fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
> O.bodies.append(utils.facet(vertices=((-0.15,-0.3,0.1*(i-
> 1)),(-0.15,-0.3,0.1*i),(-0.15,0.3,0.1*i)),dynamic=None,
> fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
>
> for i in range(1,8):
> O.bodies.append(utils.facet(vertices=((0.15,-0.3,0.1*(i-1)
> ),(0.15,0.3,0.1*(i-1)),(0.15,0.3,0.1*i)),dynamic=None,
> fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
> O.bodies.append(utils.facet(vertices=((0.15,-0.3,0.1*(i-1)
> ),(0.15,-0.3,0.1*i),(0.15,0.3,0.1*i)),dynamic=None,fixed=
> True,wire=True,color=(0.35,0.35,0.35),highlight=False,
> noBound=False,material=steel,mask=1,chain=-1))
>
> global Wall1Stressx,Wall2Stressx,Wall3Stressx,Wall4Stressx
> global Wall1Stressy,Wall2Stressy,Wall3Stressy,Wall4Stressy
> global Wall1S,Wall2S,Wall3S,Wall4S
> Wall1Stressx=0
> Wall2Stressx=0
> Wall3Stressx=0
> Wall4Stressx=0
> Wall1Stressy=0
> Wall2Stressy=0
> Wall3Stressy=0
> Wall4Stressy=0
> Wall1S=0
> Wall2S=0
> Wall3S=0
> Wall4S=0
>
> global WallStress,ConfStress,ConfDevi,MoveVel,AxiDevi,MoveAxial #stress
> control
> # Area of the confining Wall
> global A1,A2,A3,A4
> global LoadPos,IniLoadPos,plateF,IniTime,forceA
> global WallStress,ConfStress,ConfDevi,MoveVel,AxiDevi,MoveAxial,AreaPlate
> #unit:m^2
> IniTime=0
> plateF=0 #Unit:kPa
> LoadPos=0.6
> IniLoadPos=LoadPos # (link to Area of
> Walls)
> forceA=200 # Unit:kPa,P=N/A;N=P*0.0615*1000;A=0.0615
> A1=LoadPos*0.3
> A2=LoadPos*0.3
> A3=LoadPos*0.3
> A4=LoadPos*0.3
>
> WallStress=0 # Unit:kPa
> ConfStress=100 # Unit:kPa
> ConfDevi=0
> AxiDevi=0
> MoveVel=0
> MoveAxial=0
>
> AreaPlate=0.09
> # Id of different substances
> global NumLoad,NumEndBall,StepNum,NumEnd
> NumLoad=1
> NumEndBall=1
> StepNum=1
> NumEnd=1
> global xratio,yratio,zratio,NumContact,WallContact
> xratio=1
> yratio=1
> zratio=1
> NumContact=4
> WallContact=1
> # Position and other parameters record
> ######################parameters
>
> sp=pack.SpherePack()
> sp.makeCloud((-0.15,-0.15,0),(0.15,0.15,0.4),rMean=0.016,rRelFuzz=0.25)
> sp.toSimulation(material=gravel)
> NumEndBall=O.bodies[-1].id#Mark Sphere
> global iternum
> iternum=0
> #O.dt=1.0e-6 #Check it!
> O.dt=8.0e-6 #Check it!
>
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_
> Aabb(),Bo1_Facet_Aabb(),Bo1_Sphere_Aabb()]),
> InteractionLoop(
> [Ig2_Wall_Polyhedra_PolyhedraGeom(),Ig2_Polyhedra_
> Polyhedra_PolyhedraGeom(),Ig2_Facet_Polyhedra_PolyhedraGeom(
> ),Ig2_Sphere_Sphere_ScGeom(),Ig2_Sphere_Polyhedra_ScGeom(),
> Ig2_Facet_Sphere_L3Geom(),Ig2_Wall_Sphere_ScGeom()],
> [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys(),
> Ip2_FrictMat_PolyhedraMat_FrictPhys(),Ip2_FrictMat_FrictMat_FrictPhys()],
> [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric(),
> Law2_L3Geom_FrictPhys_ElPerfPl(),Law2_ScGeom_FrictPhys_CundallStrack(),
> Law2_L6Geom_FrictPhys_Linear()],
> ),
> NewtonIntegrator(damping=0.6,gravity=(0,0,-9.81),label='newton'),
> PyRunner(command='TraiStep()',iterPeriod=1,label='checker'),
> PyRunner(command='LoadAxial100kPa()',iterPeriod=
> 1,label='loadkeep100kPa'),
> PyRunner(command='AxialLoading()',iterPeriod=1,label='axialload'),
> PyRunner(command='addPlotData()',iterPeriod=1,label='plotdata'),
> PyRunner(command='Confining()',iterPeriod=1,label='keepconf'),
> #PyRunner(command='LoadAxial200kPa()',iterPeriod=
> 1,label='loadkeep200kPa'),
> ]
>
> ##Fullfill the box
> def TraiStep():
> global Wall1Stressx,Wall2Stressx,Wall3Stressx,Wall4Stressx
> global Wall1Stressy,Wall2Stressy,Wall3Stressy,Wall4Stressy
> global Wall1S,Wall2S,Wall3S,Wall4S
> global WallStress,ConfStress,ConfDevi,MoveVel,AxiDevi,MoveAxial
> #stress control,WallStress,ConfStress,ConfDevi,MoveVel
> global A1,A2,A3,A4
> global LoadPos,NumLoad,NumEndBall,IniLoadPos,plateF,IniTime,
> forceA,StepNum,NumEnd,iternum,AreaPlate
> global xratio,yratio,ztario,NumContact,WallContact
> ######
> #Step1=> add the loadingplate
> #Step2=> apply the initial axial force and confing force
> #Step3=> apply the loading force and confining stress
> if StepNum == 1:
> loadkeep100kPa.dead=True
> axialload.dead=True
> plotdata.dead=True
> keepconf.dead=True
> StepNum=StepNum+1
> elif StepNum == 2:
> #PyRunner(command='WallStressGet()',iterPeriod=1)
> #checker.command='WallStressGet()' #get the wall stress
> print "2-unbalanced forces = %.5f"%(utils.unbalancedForce()
> )
> if O.iter < 30000: return
> if utils.unbalancedForce() > 0.05: return
> #O.bodies.append(utils.wall(O.
> bodies[NumEndBall].state.pos[2]+0.04,axis=2,sense=0,material=steel))
> iternum=O.iter
> m=max([b.state.pos[2]+b.shape.radius for b in O.bodies if
> isinstance(b.shape,Sphere)])
> O.bodies.append(utils.wall(m,
> axis=2,sense=0,material=steel))
> NumLoad=O.bodies[-1].id
> NumEnd=O.bodies[-1].id
> LoadPos=O.bodies[NumLoad].state.pos[2]
> StepNum=StepNum+1
> elif StepNum == 3:
> LoadPos=O.bodies[NumLoad].state.pos[2]
> #PyRunner(command='WallStressGet()',iterPeriod=1)
> print "3-Loadplate force= %.5f"%(plateF)
> AreaPlate=(O.bodies[50].state.
> pos[0]-O.bodies[36].state.pos[0])*(O.bodies[50].state.pos[0]
> -O.bodies[36].state.pos[0])
> plateF=(O.forces.f(NumLoad)[2])/(AreaPlate*1000)
> #P=F/A=F/(0.0615*1000)=F/61.5 Unit:kPa
> LoadPos=O.bodies[NumLoad].state.pos[2]
> if plateF < 50:
> O.bodies[NumLoad].state.vel=(0,0,-0.005) #100mm/s
> O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
> else:
> O.bodies[NumLoad].state.vel=(0,0,0)
> O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
> StepNum=StepNum+1
> #O.pause(),first I got to the 200kPa axial stress, then
> keep loading axial stress
> elif StepNum == 4:
> loadkeep100kPa.dead=False
> StepNum=StepNum+1
> #O.pause()
> elif StepNum == 5:
> LoadPos=O.bodies[NumLoad].state.pos[2]
> A1=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].
> state.pos[0])
> A2=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].
> state.pos[0])
> A3=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].
> state.pos[1])
> A4=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].
> state.pos[1])
> for i in range(1,15):
> Wall1Stressx=Wall1Stressx+abs(O.forces.f(i)[0])
> Wall1Stressy=Wall1Stressy+abs(O.forces.f(i)[1])
> Wall1S=Wall1Stressy/A1
> Wall1Stressx=0
> Wall1Stressy=0
> for i in range(15,29):
> Wall2Stressy=Wall2Stressy+abs(O.forces.f(i)[1])
> Wall2S=Wall2Stressy/A2
> Wall2Stressy=0
> for i in range(29,43):
> Wall3Stressx=Wall3Stressx+abs(O.forces.f(i)[0])
> Wall3S=Wall3Stressx/A3
> Wall3Stressx=0
> #Wall3Stressy=0
> for i in range(43,57):
> Wall4Stressx=Wall4Stressx+abs(O.forces.f(i)[0])
> Wall4S=Wall4Stressx/A4
> Wall4Stressx=0
> ##########################
> WallStress=(Wall1S+Wall2S+Wall3S+Wall4S)/4000 ###
> Unit(kPa)
> ConfDevi=(abs(WallStress-ConfStress))/ConfStress ###
> Unit/kPa # parameter!!!
> for i in range(1,57):
> NumContact=NumContact+len(O.
> bodies[NumLoad].intrs())
> WallContact=NumContact/4+1
> NumContact=4
> MoveVel=1*A1*(WallStress-ConfStress)/((1.0e9)*WallContact*(8.0e-6))
>
> ################check the parameter
> if abs(MoveVel) > 0.0001:
> MoveVel=0.000001*(WallStress-ConfStress)
> else:
> print "MoveVel is OK"
> for i in range(1,57):
> xratio=(abs(O.bodies[i].state.
> pos[0]))/(O.bodies[i].state.pos[0]+0.00001)
> yratio=(abs(O.bodies[i].state.
> pos[1]))/(O.bodies[i].state.pos[1]+0.00001)
> O.bodies[i].state.vel=(
> MoveVel*xratio,MoveVel*yratio,0)
> for i in range(1,57):
> O.bodies[i].state.blockedDOFs='xyzXYZ'
> print "Ini-Cal-conf-stress= %.5f, Vel= %.8f"%(WallStress,
> O.bodies[1].state.vel[1])
> if ConfDevi > 0.05: return
> for i in range(1,57):
> O.bodies[i].state.blockedDOFs='xyzXYZ'
> print "Ini-Conf-stress= %.5f"%(WallStress)
> StepNum=StepNum+1
> elif StepNum == 6:
> loadkeep100kPa.dead=True
> #O.engines=O.engines+[PyRunner(command='Confining()'
> ,iterPeriod=1)]
> keepconf.dead=False
> StepNum=StepNum+1
> #O.pause()
> elif StepNum == 7:
> AreaPlate=(O.bodies[50].state.
> pos[0]-O.bodies[36].state.pos[0])*(O.bodies[50].state.pos[0]
> -O.bodies[36].state.pos[0])
> plateF=(O.forces.f(NumLoad)[2])/(1000*AreaPlate)
> LoadPos=O.bodies[NumLoad].state.pos[2]
> AxiDevi=abs((plateF-forceA))/forceA
> forceA=200 #10Hz Omega=2*pi*f//
> 120+80*sin(Omega*t)//40-200kPa
> zratio=1*AreaPlate/((2.0e9)*(
> len(O.bodies[NumLoad].intrs())+1)*(8.0e-6)) #alpha=50==>100
> MoveAxial=1*zratio*(plateF-forceA)
> if abs(MoveAxial) > 0.0001:
> MoveAxial=0.000001*(plateF-forceA)
> else:
> print "MoveAxial is OK"
> print "force= %.5f, ForceA= %.5f, Vel= %.8f"%(plateF,
> forceA, MoveAxial)#Unit (kPa)//0.0615 is Area of loadingplate
> O.bodies[NumLoad].state.vel=(0,0,MoveAxial)
> O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
> if AxiDevi > 0.05: return
> print "Loadplate force= %.5f, ForceA= %.5f"%(plateF,
> forceA)
> O.bodies[NumLoad].state.blockedDOFs='xyXYZ'
> IniLoadPos=LoadPos
> IniTime=O.time
> StepNum=StepNum+1
> #O.pause()
> elif StepNum == 8:
> print "8-force= %.5f"%(plateF)
> #O.engines=O.engines+[PyRunner(command='
> AxialLoading()',iterPeriod=1)]+[PyRunner(command='
> addPlotData()',iterPeriod=1)]
> axialload.dead=False
> plotdata.dead=False
> StepNum=StepNum+1
> O.pause()
> else:
> print "Well Done"
> #O.pause()
>
> ##
> def Confining():
> global Wall1Stressx,Wall2Stressx,Wall3Stressx,Wall4Stressx
> global Wall1Stressy,Wall2Stressy,Wall3Stressy,Wall4Stressy
> global Wall1S,Wall2S,Wall3S,Wall4S
> global WallStress,ConfStress,ConfDevi,MoveVel #stress control
> global A1,A2,A3,A4
> global LoadPos,NumLoad,NumEndBall
> global xratio,yratio,NumContact,WallContact #control the velocity
> of confining walls
> LoadPos=O.bodies[NumLoad].state.pos[2]
> A1=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
> A2=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
> A3=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].state.pos[1])
> A4=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].state.pos[1])
> for i in range(1,15):
> Wall1Stressx=Wall1Stressx+abs(O.forces.f(i)[0])
> Wall1Stressy=Wall1Stressy+abs(O.forces.f(i)[1])
> Wall1S=Wall1Stressy/A1
> Wall1Stressx=0
> Wall1Stressy=0
> for i in range(15,29):
> Wall2Stressy=Wall2Stressy+abs(O.forces.f(i)[1])
> Wall2S=Wall2Stressy/A2
> Wall2Stressy=0
> for i in range(29,43):
> Wall3Stressx=Wall3Stressx+abs(O.forces.f(i)[0])
> Wall3S=Wall3Stressx/A3
> Wall3Stressx=0
> #Wall3Stressy=0
> for i in range(43,57):
> Wall4Stressx=Wall4Stressx+abs(O.forces.f(i)[0])
> Wall4S=Wall4Stressx/A4
> Wall4Stressx=0
> ##########################
> WallStress=(Wall1S+Wall2S+Wall3S+Wall4S)/4000 ### Unit(kPa)
> ConfDevi=(abs(WallStress-ConfStress))/ConfStress ### Unit/kPa #
> parameter!!!
> if ConfDevi > 0.05:
> for i in range(1,57):
> NumContact=NumContact+len(O.
> bodies[NumLoad].intrs())
> WallContact=NumContact/4+1
> NumContact=4
> MoveVel=1*A1*(WallStress-ConfStress)/((1.0e9)*
> WallContact*(8.0e-6))
> if abs(MoveVel) > 0.0001:
> MoveVel=0.000001*(WallStress-ConfStress)
> else:
> print "MoveVel is OK"
> for i in range(1,57):
> #zratio=0.5*AreaPlate/((2.0e9)
> *(len(O.bodies[NumLoad].intrs())+1)*(8.0e-6))
> xratio=(abs(O.bodies[i].state.
> pos[0]))/(O.bodies[i].state.pos[0]+0.00001)
> yratio=(abs(O.bodies[i].state.
> pos[1]))/(O.bodies[i].state.pos[1]+0.00001)
> O.bodies[i].state.vel=(
> MoveVel*xratio,MoveVel*yratio,0)
> for i in range(1,57):
> O.bodies[i].state.blockedDOFs='xyzXYZ'
> print "Keep-Cal-conf-stress= %.5f, Vel= %.8f"%(WallStress,
> MoveVel)
> else:
> for i in range(1,57):
> O.bodies[i].state.blockedDOFs='xyzXYZ'
> print "Keep-Conf= %.5f"%(WallStress)
> ############## keep confining pressure
>
> def LoadAxial100kPa():
> global forceA,plateF,LoadPos,AxiDevi,IniTime,NumLoad,MoveAxial,
> AreaPlate,zratio
> AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[
> 0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
> plateF=(O.forces.f(NumLoad)[2])/(1000*AreaPlate)
> LoadPos=O.bodies[NumLoad].state.pos[2]
> AxiDevi=abs((plateF-forceA))/forceA
> forceA=100 #10Hz Omega=2*pi*f// 120+80*sin(Omega*t)//40-200kPa
> zratio=5*AreaPlate/((2.0e9)*(len(O.bodies[NumLoad].intrs())
> +1)*(8.0e-6))
> MoveAxial=1*zratio*(plateF-forceA)
> if abs(MoveAxial) > 0.0001:
> MoveAxial=0.000001*(plateF-forceA)
> else:
> print "MoveAxial is OK"
> if AxiDevi > 0.05:
> print "keep-force= %.5f, ForceA= %.5f, Vel= %.8f, CNum=
> %.1f"%(plateF, forceA, MoveAxial, len(O.bodies[NumLoad].intrs()))#Unit
> (kPa)//0.0615 is Area of loadingplate
> O.bodies[NumLoad].state.vel=(0,0,MoveAxial)
> O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
> else:
> print "Done-keep-force= %.5f, ForceA= %.5f, Vel= %.8f,
> CNum= %.1f"%(plateF, forceA, MoveAxial, len(O.bodies[NumLoad].intrs()))
>
> O.bodies[NumLoad].state.blockedDOFs='xyXYZ'
>
> def LoadAxial200kPa():
> global forceA,plateF,LoadPos,AxiDevi,IniTime,NumLoad,MoveAxial,
> AreaPlate,zratio
> AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[
> 0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
> plateF=(O.forces.f(NumLoad)[2])/(1000*AreaPlate)
> LoadPos=O.bodies[NumLoad].state.pos[2]
> AxiDevi=abs((plateF-forceA))/forceA
> forceA=200 #10Hz Omega=2*pi*f// 120+80*sin(Omega*t)//40-200kPa
> zratio=5*AreaPlate/((2.0e9)*(len(O.bodies[NumLoad].intrs())
> +1)*(8.0e-6))
> MoveAxial=1*zratio*(plateF-forceA)
> if abs(MoveAxial) > 0.0001:
> MoveAxial=0.000001*(plateF-forceA)
> else:
> print "MoveAxial is OK"
> if AxiDevi > 0.05:
> print "keep-force= %.5f, ForceA= %.5f, Vel= %.8f, CNum=
> %.1f"%(plateF, forceA, MoveAxial, len(O.bodies[NumLoad].intrs()))#Unit
> (kPa)//0.0615 is Area of loadingplate
> O.bodies[NumLoad].state.vel=(0,0,MoveAxial)
> O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
> else:
> print "Done-keep-force= %.5f, ForceA= %.5f, Vel= %.8f,
> CNum= %.1f"%(plateF, forceA, MoveAxial, len(O.bodies[NumLoad].intrs()))
>
> O.bodies[NumLoad].state.blockedDOFs='xyXYZ'
> ##AxialLoading
> def AxialLoading():
> global forceA,plateF,LoadPos,AxiDevi,IniTime,NumLoad,MoveAxial,
> AreaPlate,zratio
> AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[
> 0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
> plateF=(O.forces.f(NumLoad)[2])/(1000*AreaPlate)
> LoadPos=O.bodies[NumLoad].state.pos[2]
> AxiDevi=abs((plateF-forceA))/forceA
> forceA=200+50*sin((20*pi)*(O.time-IniTime)) #10Hz Omega=2*pi*f//
> 120+80*sin(Omega*t)//40-200kPa
> zratio=0.5*AreaPlate/((2.0e9)*(len(O.bodies[NumLoad].intrs()
> )+1)*(8.0e-6))
> MoveAxial=1*zratio*(plateF-forceA)
> if abs(MoveAxial) > 0.0001:
> MoveAxial=0.000001*(plateF-forceA)
> else:
> print "MoveAxial is OK"
> if AxiDevi > 0.05:
> print "final-force= %.5f, ForceA= %.5f, Vel=
> %.8f"%(plateF, forceA, MoveAxial)#Unit (kPa)//0.0615 is Area of loadingplate
> O.bodies[NumLoad].state.vel=(0,0,MoveAxial)
> O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
> else:
> print "final-Loadplate force= %.5f, ForceA= %.5f"%(plateF,
> forceA)
> O.bodies[NumLoad].state.blockedDOFs='xyXYZ'
>
>
>
>
> ##Record
> def addPlotData():
> global LoadPos,IniLoadPos,NumLoad,forceA,plateF
> global theta,thega,WallStress,Vol,AreaPlate
> theta=forceA
> theta2=plateF
> LoadPos=O.bodies[NumLoad].state.pos[2]
> thega=((IniLoadPos-LoadPos)/IniLoadPos)*100
> AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[
> 0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
> Vol=IniLoadPos*AreaPlate-LoadPos*(O.bodies[1].state.
> pos[1])*(O.bodies[1].state.pos[1])
> plot.addData(Thega=thega,Theta=theta,Thega2=thega,
> Theta2=theta2,T=O.time,Conf=WallStress,TimeLast=O.time,Volume=Vol)
>
> ##ConfiningWall=>wall(1-5)//wall(2-6)//wall(3-7)//wall(4-8)
>
> plot.plots={'Thega':('Theta',),'Thega2':('Theta2',),'T':('
> Conf',),'TimeLast':('Volume',)}
> plot.plot()
>
> qt.Controller()
> V = qt.View()
> V.screenSize = (550,450)
> V.sceneRadius = 1
> V.eyePosition = (0.7,0.5,0.1)
> V.upVector = (0,0,1)
> V.lookAt = (0.15,0.15,0.1)
>
>
>
>
> --
> 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 : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Dear Klaus.thoeni,
Thanks for your reply! First I wanted to do a non-sphere assembly triaxial test actually, so I tried a sphere test in advance. As my script is too long to check, I simplified script as following for you to help me.
####-----####
from yade import *...
....# Definition materials parameters...
...# Generation of six boundary walls and sphere particles...
...O.engines=[
          ForceResetter(),...
          InsertionSortCollider()...
          IneractionLoop()...
          PyRunner(command=‘TraiStep()’,iterPeriod=1,label='checker'),
          PyRunner(command=‘LoadAxial()’,iterPeriod=1,label='loadkeep',dead=True),
          PyRunner(command=‘AxialLoading()’,iterPeriod=1,label='axialload',dead=True),
          PyRunner(command=‘addPlotData()’,iterPeriod=1,label='plotdata',dead=True),
          PyRunner(command=‘Confining()’,iterPeriod=1,label='keepconf',dead=True),
]
###
def TraiStep():
      StepNum == 1:
#all last four PyRunner.dead=True
      StepNum=StepNum+1
     StepNum == 2:
#when the unbalanced Force got to below 0.01, add a loading plate to compact the sphere sample
    StepNum=StepNum+1
     StepNum == 3:
#when loading plate stress get to 100kPa, skip to next StepNum
    StepNum=StepNum+1
     StepNum == 4:
    loadkeep.dead=False
# Open the PyRunner'LoadAxial()’ to keep the loading plate Stress at 100kPa
    StepNum=StepNum+1
   StepNum == 5:
# To make a servo to keep the confining stress to 100kPa
    StepNum=StepNum+1
.....

def LoadAxial(): # to keep the loading plate stress to 100kPa
def AxialLoading(): # To load after sample generation
def Confining(): # To keep confining pressure to 100kPa
#############################################################
The problem is when at the StepNum == 4, using the 'loadkeep.dead=False' to open the PyRunner LoadAxial(), and at the StepNum==5, the loading plate and confining walls disappeared, I think it might be caused the walls and loading plate were applied large velocity, when I cancel the PyRunner LoadAxial(), the following procedure was normal, and confining walls keep moving to maintain the confining pressure.
I have been confused several days for that PyRunner(), how to solve this discontiunous situation?

Dear @bruno-chareyre and @Klaus.thoeni,
       Thank you for your reply, I checked the script again attentively and I found that I didn't fix the loading wall after compaction of assembly which caused that spheres and walls flew away. The PyRunner[command=''] worked regularly. I will check scripts more carefully in the future and follow the questions rule!
Thank you again!
Best regards,
De Zhang

Thanks Bruno Chareyre, that solved my question.