Question #669404 on Yade changed: https://answers.launchpad.net/yade/+question/669404
Status: Open => Answered Klaus Thoeni proposed the following answer: 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 < question669...@answers.launchpad.net> 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 : yade-users@lists.launchpad.net > Unsubscribe : https://launchpad.net/~yade-users > More help : https://help.launchpad.net/ListHelp > -- 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 : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp