yade-users team mailing list archive

Thread
Date
[Question #669404]: A bug for using of Pyrunner[]

To: yade-users@xxxxxxxxxxxxxxxxxxx
From: De zhang <question669404@xxxxxxxxxxxxxxxxxxxxx>
Date: Sun, 20 May 2018 03:50:58 -0000
Reply-to: question669404@xxxxxxxxxxxxxxxxxxxxx
Sender: bounces@xxxxxxxxxxxxx
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.