def initTest():
	global mode
	print "init"
	if O.iter>0:
		O.wait();
		O.loadTmp('initial')
		print "Reversing plot data"; plot.reverseData()
	else: plot.plot(subPlots=False)
	strainer.strainRate=abs(strainRateTension) if mode=='tension' else -abs(strainRateCompression)
	try:
		from yade import qt
		renderer=qt.Renderer()
		renderer.dispScale=(1000,1000,1000) if mode=='tension' else (100,100,100)
	except ImportError: pass
	print "init done, will now run."
	O.step(); # to create initial contacts
	# now reset the interaction radius and go ahead
	ss2sc.interactionDetectionFactor=1.
	is2aabb.aabbEnlargeFactor=1.

	O.run()

	def plot(self):
		try:
			from yade import plot
			if len(plot.plots)==0: return None
			fig=plot.plot(subPlots=True,noShow=True)
			img=O.tmpFilename()+'.'+plotImgFormat
			sqrtFigs=math.sqrt(len(plot.plots))
			fig.set_size_inches(5*sqrtFigs,7*sqrtFigs)
			fig.savefig(img)
			f=open(img,'rb'); data=f.read(); f.close(); os.remove(img)
			#print 'returning '+plotImgFormat
			return xmlrpclib.Binary(data)
		except:
			print 'Error updating plots:'
			import traceback
			traceback.print_exc()
			return None

  TranslationEngine(translationAxis=[0,0,1],velocity=-2.0,ids=idTop,dead=False,label='translat'),
  
  CombinedKinematicEngine(ids=idTop,label='combEngine',dead=True) + 
    ServoPIDController(axis=[0,0,1],maxVelocity=2.0,iterPeriod=1000,ids=idTop,target=1.0e7,kP=1.0,kI=1.0,kD=1.0) + 
    RotationEngine(rotationAxis=(0,0,1), angularVelocity=10.0, rotateAroundZero=True, zeroPoint=(0,0,0)),
  PyRunner(command='addPlotData()',iterPeriod=1000, label='graph'),
  PyRunner(command='switchTranslationEngine()',iterPeriod=45000, nDo = 2, label='switchEng'),
]
O.step()
from yade import qt
qt.View()
r=qt.Renderer()
r.bgColor=1,1,1  

def addPlotData():
  fMove = Vector3(0,0,0)
  
  for i in idTop:
    fMove += O.forces.f(i)
  
  plot.addData(z=O.iter, pMove=fMove[2], pFest=fMove[2])

def switchTranslationEngine():
  print "Switch from TranslationEngine engine to ServoPIDController"
  translat.dead = True
  combEngine.dead = False
  


plot.plots={'z':('pMove','pFest')}; plot.plot()

	jumper.rot=Quaternion(Vector3.UnitZ,1/(1e4*nSteps))

O.run(nSteps,True)
if mode=='mov':
	plot.plots={'zShift-DD':(
		#('epsT-DD','g-'),('epsT-Sc','r^'),
		('dist-DD','g-'),('dist-Sc','r^'),
		None,
		('sigmaT-DD','b-'),('sigmaT-Sc','m^')
		#('relResStr-DD','b-'),('relResStr-Sc','m^')
		#('epsN-DD','b-'),('epsN-Sc','m^')
	)}
elif mode=='rot':
	plot.plots={'zRot-DD':(
		('epsT-DD','g|'),('epsT-Sc','r-'),
		None,
		('sigmaT-DD','b-'),('sigmaT-Sc','mv')
	)}



if 1:
	f='/tmp/cpm-geom-'+mode+'.pdf'
	plot.plot(noShow=True,subPlots=False).savefig(f)
	print 'Plot saved to '+f
	quit()
else:
	plot.plot(subPlots=False)

]

#### dataCollector
from yade import plot
plot.plots={'iterations':'v','x':'z'}

def dataCollector():
	R=O.bodies[refPoint]
	plot.addData(v=R.state.vel[2],z=R.state.pos[2],x=R.state.pos[0],iterations=O.iter,t=O.realtime)

#### joint strength degradation
stableIter=1000
stableVel=0.001
degrade=True
def jointStrengthDegradation():
	global degrade
	if degrade and O.iter>=stableIter and abs(O.bodies[refPoint].state.vel[2])<stableVel :
		print('Equilibrium reached \nJoint cohesion canceled now !', ' | iteration=', O.iter)
		degrade=False
	for i in O.interactions:
		if i.phys.isOnJoint : 
			if i.phys.isCohesive:
				i.phys.isCohesive=False
				i.phys.FnMax=0.
				i.phys.FsMax=0.

print('Seeking after an initial equilibrium state')
print('')
O.run(10000)
plot.plot()# note the straight trajectory (z(x) plot)during sliding step (before free fall) despite the discretization of joint plane with spheres

  f3 = [0,0,0]
  f4 = [0,0,0]
  
  try:
    f1=O.forces.f(id12)
    f2=O.forces.f(id22)
    f3=O.forces.f(id32)
    f4=O.forces.f(id42)
  except:
    f1 = [0,0,0]
    f2 = [0,0,0]
    f3 = [0,0,0]
    f4 = [0,0,0]
  
  s1 = (O.bodies[id12].state.pos[2]-O.bodies[id11].state.pos[2])-(r1 + r2)
  s2 = (O.bodies[id22].state.pos[2]-O.bodies[id21].state.pos[2])-(r1 + r2)
  s3 = (O.bodies[id32].state.pos[2]-O.bodies[id31].state.pos[2])-(r1 + r2)
  s4 = (O.bodies[id42].state.pos[2]-O.bodies[id41].state.pos[2])-(r1 + r2)
  
    
  plot.addData(fc1=f1[2], sc1=-s1, fc2=f2[2], sc2=-s2, fc3=f3[2], sc3=-s3, fc4=f4[2], sc4=-s4)

plot.plots={'sc1':('fc1'), 'sc2':('fc2'), 'sc3':('fc3'), 'sc4':('fc4')}; plot.plot()

from yade import qt
qt.View()

O.run(320000)
O.wait()
plot.saveGnuplot('sim-data_LudigPM')

        return
    Fz = O.forces.f(plate.id)[2]
    plot.addData(
        Fz=Fz,
        w=plate.state.pos[2] - (-4*Diameter),
        unbalanced=unbalancedForce(),
        i=O.iter
    )


def defVisualizer():
   with open("data.dat","a") as f:
       for b in O.bodies:
           if isinstance(b.shape, Sphere):
               rData = "{x},{y},{z},{r},{w}\t".format(x = b.state.pos[0],
                                                y = b.state.pos[1],
                                                z = b.state.pos[2],
                                                r = b.shape.radius + b.state.dR,
                                                w = plate.state.pos[2]
                                               )
               f.write(rData)
       f.write("\n")



O.timingEnabled=True
O.run(1, True)
plot.plots={'w':('Fz', None)}
plot.plot()

from yade import qt
qt.View()
qt.Controller()

############################################
##### now the part pertaining to plots #####
############################################

from yade import plot
## make one plot: step as function of fn
plot.plots={'un':('fn')}

## this function is called by plotDataCollector
## it should add data with the labels that we will plot
## if a datum is not specified (but exists), it will be NaN and will not be plotted

def myAddPlotData():
	if O.interactions[0,1].isReal:
		i=O.interactions[0,1]
		## store some numbers under some labels
		plot.addData(fn=i.phys.normalForce[0],step=O.iter,un=2*s0.shape.radius-s1.state.pos[0]+s0.state.pos[0],kn=i.phys.kn)	

O.run(100,True); plot.plot(subPlots=False)

## We will have:
## 1) data in graphs (if you call plot.plot())
## 2) data in file (if you call plot.saveGnuplot('/tmp/a')
## 3) data in memory as plot.data['step'], plot.data['fn'], plot.data['un'], etc. under the labels they were saved

  ),
  NewtonIntegrator(damping=0,gravity=[0,0,-9.81]),
  PyRunner(command='addPlotData()',iterPeriod=100),
]

vel=-0.50
O.bodies[id12].state.vel=[0,0,vel]
O.bodies[id22].state.vel=[0,0,vel]
O.bodies[id32].state.vel=[0,0,vel]

def addPlotData():
  s1 = (O.bodies[id12].state.pos[2]-O.bodies[id11].state.pos[2])-(r1+r2)
  s2 = (O.bodies[id22].state.pos[2]-O.bodies[id11].state.pos[2])-(r1+r2)
  s3 = (O.bodies[id32].state.pos[2]-O.bodies[id11].state.pos[2])-(r1+r2)
  plot.addData(mat1mat2=s1,mat1mat3=s2,mat2mat3=s3,it=O.iter)
  
  

plot.plots={'it':('mat1mat2','mat1mat3','mat2mat3')}; plot.plot()

O.step()
from yade import qt
qt.View()

print "Friction coefficient for id11 and id12 is %g"%(math.atan(O.interactions[id11,id12].phys.tangensOfFrictionAngle))
print "Friction coefficient for id21 and id22 is %g"%(math.atan(O.interactions[id21,id22].phys.tangensOfFrictionAngle))
print "Friction coefficient for id31 and id32 is %g"%(math.atan(O.interactions[id31,id32].phys.tangensOfFrictionAngle))

O.run(100000, True)
#plot.saveGnuplot('sim-data_')

id1 = O.bodies.append(sphere(center=[0,0,2*r],radius=r,material=mat1))
id2 = O.bodies.append(geom.facetBox(center=(0,-16.0*r,-2*r),orientation=oriBody,extents=(r,17.0*r,0), material=mat1,color=(0,0,1)))

id3 = O.bodies.append(sphere(center=[10*r,0,2*r],radius=r,material=mat2))
id4 = O.bodies.append(geom.facetBox(center=(10*r,-16.0*r,-2*r),orientation=oriBody,extents=(r,17.0*r,0), material=mat2,color=(0,0,1)))

o.engines = [
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
    [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
    [Law2_ScGeom_ViscElPhys_Basic()],
  ),
  NewtonIntegrator(damping=0,gravity=[0,0,-9.81]),
  PyRunner(command='addPlotData()',iterPeriod=10000, dead = False, label='graph'),
]

def addPlotData(): 
  f1 = [0,0,0]
  s1 = O.bodies[id1].state.pos[1]
  s2 = O.bodies[id3].state.pos[1]
  plot.addData(sc=O.time, fc1=s1, fc2=s2)
  
  

plot.plots={'sc':('fc1','fc2')}; plot.plot()

from yade import qt
qt.View()

#### dataCollector
from yade import plot
plot.plots={'iterations':('v',)}
#plot.plots={'x':('z',)}
def dataCollector():
	R=O.bodies[refPoint]
	plot.addData(v=R.state.vel[2],z=R.state.pos[2],x=R.state.pos[0],iterations=O.iter,t=O.realtime)

#### joint strength degradation
stableIter=2000
stableVel=0.001
degrade=True
def jointStrengthDegradation():
    global degrade
    if degrade and O.iter>=stableIter and abs(O.bodies[refPoint].state.vel[2])<stableVel :
	print '!joint cohesion total degradation!', ' | iteration=', O.iter
	degrade=False
	for i in O.interactions:
	    if i.phys.isOnJoint : 
		if i.phys.isCohesive:
		  i.phys.isCohesive=False
		  i.phys.FnMax=0.
		  i.phys.FsMax=0.

O.step()

print 'Seeking after an initial equilibrium state'
print ''
O.run(10000)
plot.plot()# note the straight line (during sliding step, before free fall) despite the discretization of joint plane with spheres

O.bodies[id2].state.vel=[0,0,vel]
O.bodies[id4].state.vel=[0,0,vel]
O.bodies[id6].state.vel=[0,0,vel]
O.bodies[id8].state.vel=[0,0,vel]
O.bodies[id10].state.vel=[0,0,vel]
O.bodies[id12].state.vel=[0,0,vel]

def addPlotData(): 
  f1=O.forces.f(id2)
  f2=O.forces.f(id4)
  f3=O.forces.f(id6)
  f4=O.forces.f(id8)
  f5=O.forces.f(id10)
  f6=O.forces.f(id12)
  
  s1 = (O.bodies[id2].state.pos[2]-O.bodies[id1].state.pos[2])-2*r
  sc=s1

  plot.addData(Willett_numeric=-f1[2], Willett_analytic=-f2[2], Rabinovich=-f3[2], Lambert=-f4[2], Weigert=-f5[2], Soulie=-f6[2], sc=sc)
  
  

plot.plots={'sc':('Willett_numeric','Willett_analytic','Rabinovich','Lambert','Weigert','Soulie')}; plot.plot()

O.step()
from yade import qt
qt.View()

O.run(250000, True)
#plot.saveGnuplot('sim-data_'+CapillarType1)

  global iterN
  if (O.bodies[id2].state.vel[2]<0.0): 
    O.bodies[id2].state.vel*=-1.0
    cDir.iterPeriod = int(iterN/10000.0)
  elif (O.bodies[id2].state.pos[2]-O.bodies[id1].state.pos[2])-(r1 + r2) > 0.0:
    iterN = int(iterN*1.2)
    O.bodies[id2].state.vel[2]*=-1.0
    cDir.iterPeriod = iterN
    
def addPlotData(): 
  f = [0,0,0]
  sc = 0
  try:
    f=O.forces.f(id2)
  except:
    f = [0,0,0]
  
  s1 = (O.bodies[id2].state.pos[2]-O.bodies[id1].state.pos[2])-(r1 + r2)
  
  fc1 = f[2]
  sc1 = -s1/r1
    
  plot.addData(fc1=fc1, sc=sc1)

plot.plots={'sc':('fc1')}; plot.plot()

from yade import qt
qt.View()

plot.saveGnuplot('sim-data_LudigPM')

integratoreng.abs_err=1e-3;

# We use only the integrator engine

# Time step determines the exiting period of the integrator since the integrator performs one step from current_time to current_time+dt;

O.dt=1e-3;

O.engines=[
#		integratoreng,
		ForceResetter(),
#		Apply internal force to the deformable elements and internal force of the interaction element
		FEInternalForceEngine([If2_Lin4NodeTetra_LinIsoRayleighDampElast(),If2_2xLin4NodeTetra_LinCohesiveStiffPropDampElastMat()]),
		PyRunner(iterPeriod=1,command='applyforcetoelements()'),
		NewtonIntegrator(damping=0,gravity=[0,0,0]),
#		Plotting data: adds plots after one step of the integrator engine
		PyRunner(iterPeriod=1,command='addplot()')
	]

from yade import plot

plot.plots={'t':'vel','time':'pos','tm':'force','tt':'postail'}
plot.plot(subPlots=True)

try:
	from yade import qt
	qt.View()
	qt.Controller()
except ImportError: pass

	jumper.rot=Quaternion(Vector3.UnitZ,1/(1e4*nSteps))

O.run(nSteps,True)
if mode=='mov':
	plot.plots={'zShift-DD':(
		#('epsT-DD','g-'),('epsT-Sc','r^'),
		('dist-DD','g-'),('dist-Sc','r^'),
		None,
		('sigmaT-DD','b-'),('sigmaT-Sc','m^')
		#('relResStr-DD','b-'),('relResStr-Sc','m^')
		#('epsN-DD','b-'),('epsN-Sc','m^')
	)}
elif mode=='rot':
	plot.plots={'zRot-DD':(
		('epsT-DD','g|'),('epsT-Sc','r-'),
		None,
		('sigmaT-DD','b-'),('sigmaT-Sc','mv')
	)}



if 1:
	f='/tmp/cpm-geom-'+mode+'.pdf'
	plot.plot(noShow=True).savefig(f)
	print 'Plot saved to '+f
	quit()
else:
	plot.plot()

  SPHEngine(mask=1, k=k, rho0 = rho, KernFunctionDensity=1),
  PyRunner(command='addPlotData()',iterPeriod=1,dead=False),
]

print("Time\tX\tRho\tP\tFpr ")

# Function to add data to plot
def addPlotData(): 
  #print "%.2f\t%.5f\t%.5f\t%.5f\t%.5f" % (O.time+O.dt, O.bodies[id2].state.pos[2], O.bodies[id2].rho, O.bodies[id2].press, O.forces.f(id2)[2])
  s1 = (O.bodies[id12].state.pos[2]-O.bodies[id11].state.pos[2])-(Rad*2.0)
  s2 = (O.bodies[id22].state.pos[2]-O.bodies[id21].state.pos[2])-(Rad*2.0)
  s3 = (O.bodies[id32].state.pos[2]-O.bodies[id31].state.pos[2])-(Rad*2.0)
  s4 = (O.bodies[id42].state.pos[2]-O.bodies[id41].state.pos[2])-(Rad*2.0)
  s5 = (O.bodies[id52].state.pos[2]-O.bodies[id51].state.pos[2])-(Rad*2.0)
  
  f1 = O.forces.f(id12)[2]
  f2 = O.forces.f(id22)[2]
  f3 = O.forces.f(id32)[2]
  f4 = O.forces.f(id42)[2]
  f5 = O.forces.f(id52)[2]
  
  plot.addData(sc=O.iter, s1 = s1, s2 = s2, s3 = s3, s4 = s4, s5 = s5,
               fc=O.iter, f1 = f1, f2 = f2, f3 = f3, f4 = f4, f5 = f5)

plot.plots={'sc':('s1', 's2', 's3', 's4', 's5'), 'fc':('f1', 'f2', 'f3', 'f4', 'f5')}; plot.plot()

O.run(1, True)

qt.View()
O.run(1500)

############################################

from yade import plot
## we will have 2 plots:
## 1. t as function of i (joke test function)
## 2. i as function of t on left y-axis ('|||' makes the separation) and z_sph, v_sph (as green circles connected with line) and z_sph_half again as function of t
plot.plots={'t':('error')}

def myAddPlotData():
	symId=0
	numId=1
	O.bodies[symId].state.update()
	psiDiff=((O.bodies[symId].state)-(O.bodies[numId].state))	
	plot.addData(t=O.time,error=(psiDiff|psiDiff).real)
plot.liveInterval=.2
plot.plot(subPlots=False)


try:
	from yade import qt
	Gl1_QMGeometry().analyticUsesStepOfDiscrete=False
	qt.View()
	qt.Controller()
	qt.controller.setWindowTitle("1D eigenwavefunction in harmonic potential")
	qt.controller.setViewAxes(dir=(0,1,0),up=(0,0,1))
	qt.Renderer().blinkHighlight=False
except ImportError:
	pass

#O.run(20000)

o.engines = [
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(label='is2aabb'),Bo1_Facet_Aabb(label='is3aabb')]),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(label='ss2sc'),Ig2_Facet_Sphere_ScGeom()],
    [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
    [Law2_ScGeom_ViscElPhys_Basic()],
  ),
  NewtonIntegrator(damping=0.05,gravity=[0,-9.81,0]),
  SPHEngine(mask=3, k=k, rho0 = rho, h = h, KernFunctionDensity= 1),
  VTKRecorder(iterPeriod=100,fileName='./cpt/spheres-', recorders=['spheres','velocity','colors','intr','ids','mask','materialId','stress']),
  VTKRecorder(iterPeriod=100,fileName='./cpt/facet-',   recorders=['facets'],label='VTK_box2'),
  PyRunner(command='addPlotData()',iterPeriod=50,dead=False),
]

def addPlotData(): 
  plot.addData(t=O.time, Ekin=utils.kineticEnergy())

plot.plots={'t':('Ekin')}; plot.plot()

enlargeF = h/Rad*1.1
print("enlargeF = %g"%enlargeF)
is2aabb.aabbEnlargeFactor = enlargeF
ss2sc.interactionDetectionFactor = enlargeF

O.run(1, True)

qt.View()
#O.run(10000, True)
#plot.saveGnuplot('sim-data_0.05')

#### initializes now the interaction detection factor
aabb.aabbEnlargeFactor=-1.
Ig2ssGeom.interactionDetectionFactor=-1.


#### define engines for simulation with UniaxialStrainer
O.engines = O.engines[:3] + [ UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=1,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=True,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
NewtonIntegrator(damping=0.5),
PyRunner(initRun=True,iterPeriod=1000,command='addPlotData()'),
]


#### plot some results
plot.plots={'un':('Fn',)}
plot.plot(noShow=False, subPlots=False)

def addPlotData():
	Fn = 0.
	for i in posIds:
		try:
			inter=O.interactions.withBody(i)[0]
			F = abs(inter.phys.normalForce[1])
		except:
			F = 0
		Fn += F
	un = O.bodies[O.bodies[posIds[0]].id].state.pos[1] - O.bodies[O.bodies[posIds[0]].id].state.refPos[1]
	if un > 0.10:
		O.pause()
	plot.addData( un=un*1000, Fn=Fn/1000 )

# Add engines
o.engines = [
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb()]),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom()],
    [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
    [Law2_ScGeom_ViscElPhys_Basic()],
  ),
  NewtonIntegrator(damping=0,gravity=[0,0,-9.81]),
  PyRunner(command='addPlotData()',iterPeriod=100),
]

# Function to add data to plot
def addPlotData(): 
  try:
    delta = (O.bodies[id2].state.pos[2]-O.bodies[id1].state.pos[2])-(2*Rad)
    plot.addData(delta=delta, time1=O.time, time2=O.time, time3=O.time, time4=O.time,
              Fn = O.interactions[0,1].phys.Fn,
              Fv = O.interactions[0,1].phys.Fv,
              deltaDot = O.bodies[id2].state.vel[2] - O.bodies[id1].state.vel[2])
  except:
    pass
  
plot.plots={'time1':('delta'), 'time2':('deltaDot'), 'time3':('Fn'), 'time4':('Fv')}; plot.plot()

O.run(1)
qt.View()

#O.wait() ; plot.saveGnuplot('sim-data_Sphere')