yade-users team mailing list archive

[Václav Šmilauer <eudoxos@xxxxxxxx>]

> I want to run a simple test of two balls in displacement control. Is
> there any engine for this purpose? I see in ForceContainer a function
> called move.. Is that one I should use for? Can I run it say every
> time step?
Use TranslationEngine or JumpChangeSe3. You can use move as well, but
you have to set it at every step; besides, it will not set particle's
velocity.

I am attaching script that I use(d) for plotting 1-1 graphs, take it as
inspiration, it is rather hairy. I wanted to add something like this to
some yade standard module (since every time you develop new contact law,
you want to test it on 1-1 configuration first), but I did not (and will
not probably) have time for that. It uses the JumpChangeSe3 engine.

It runs 3 scenarios, the one you might need is distProx perhaps.

HTH, Vaclav

import matplotlib
#matplotlib.rc('font',size=14)
#matplotlib.rc('figure',figsize=(6,4))
matplotlib.rc('text',usetex=True)
matplotlib.rc('text.latex',preamble=r'\usepackage{euler}')

from yade import utils, eudoxos, log
#log.setLevel('JumpChangeSe3',log.INFO)
log.setLevel('CpmPhys',log.TRACE)
log.setLevel('Law2_Dem3DofGeom_CpmPhys_Cpm',log.TRACE)
from math import *
import pylab
import re



# material
young=30e9
poisson=.2
G_over_E=.2
sigmaT=3e6
frictionAngle=atan(0.8)
epsCrackOnset=1e-4
relDuctility=5.
intRadius=2

# geometry
#r1,r2=1.0,1.0 # radii
r1,r2=1e-3,1e-3
dist=r1+r2
p1,p2=[0,0,0],[r1+r2,0,0] # center positions

O.materials.append(CpmMat(young=young,poisson=poisson,frictionAngle=frictionAngle,sigmaT=sigmaT,relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,G_over_E=G_over_E,neverDamage=True,plTau=-1,plRateExp=0,dmgTau=-1,dmgRateExp=0))

O.bodies.append(utils.sphere(p1,r1))
O.bodies.append(utils.sphere(p2,r2))

for b in O.bodies: b.dynamic=False
O.dt=1e-3 #1e-6 #.5*utils.PWaveTimeStep()

O.engines=[
	JumpChangeSe3(subscribedBodies=[1],label='jumper'),
	ForceResetter(),
	BoundDispatcher([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius)]),
	InsertionSortCollider(),
	InteractionDispatchers(
		[Ig2_Sphere_Sphere_Dem3DofGeom(distFactor=intRadius),],
		[Ip2_CpmMat_CpmMat_CpmPhys(label='cpmPhysFunctor')],
		[Law2_Dem3DofGeom_CpmPhys_Cpm(omegaThreshold=1.,label='cpmLaw')]
	),
	NewtonIntegrator(damping=0),
	PeriodicPythonRunner(iterPeriod=1,command='myAddPlotData()',label='plotDataCollector')
]

import yade.plot as yp

yp.labels=dict(Fn=r'$F_N$',Fs=r'$|F_S|$',epsN=r'$\varepsilon_N$',epsT=r'$|\varepsilon_T|$',sigmaN=r'$\sigma_N$',sigmaT=r'$|\sigma_T|$',epsNPl=r'$\eps_N^{\rm pl}$',omega=r'$\omega$')

def myAddPlotData():
	st=O.bodies[1].state
	iPhys=O.interactions[0,1].phys
	zRotArcLen=st.ori.ToAxisAngle()[1]*r2
	yade.plot.addData(dy1=st.pos[1]-p2[1],Fn=iPhys.Fn,Fs=iPhys.Fs.Length(),epsN=iPhys.epsN,epsT=iPhys.epsT.Length(),sigmaN=iPhys.sigmaN,sigmaT=iPhys.sigmaT.Length(),epsNPl=iPhys.epsNPl,
		dx1=st.pos[0]-p2[0],
		dz1=st.pos[2]-p2[2],
		zRotArcLen=zRotArcLen,
		fx0=O.forces.f(0)[0],fy0=O.forces.f(0)[1],fz0=O.forces.f(0)[2],mx0=O.forces.t(0)[0],my0=O.forces.t(0)[1],mz0=O.forces.t(0)[2],
		fx1=O.forces.f(1)[0],fy1=O.forces.f(1)[1],fz1=O.forces.f(1)[2],mx1=O.forces.t(1)[0],my1=O.forces.t(1)[1],mz1=O.forces.t(1)[2],
		omega=iPhys.omega
	)
	
def printInteraction(i,mode=''):
	print '=========== mode %s, iteration %d; bodies %d %d ============'%(mode,O.iter,i.id1,i.id2)
	print "Bodies at ",O.bodies[i.id1].state.pos,O.bodies[i.id2].state.pos
	for attr in ['epsN','epsT','Fn','Fs','E','G','omega','crossSection']:
		print attr,'=',getattr(i.phys,attr)
	print 'contactPoint =',i.geom.contactPoint

if 0:
	try:
		from yade import qt
		r=qt.Renderer()
		r.dispScale=1000,1000,1000
		r.intrPhys=True
		#O.miscParams=[Generic("GLDrawBrefcomContact",{'dmgLabel':False,'dmgPlane':True,'colorStrain':False,'epsNLabel':False,'epsT':False,'epsTAxes':False,'normal':False,'contactLine':True})]
		utils.wireAll()
		qt.View()
	except ImportError: pass

O.saveTmp('init')

import time
for mode in ['viscDist','cycleProxDist','distProx']: #'viscShear','viscDist','distProx','yShear','zRot','distProx'][-1:]:# [0:1]: #[-1:]:
	O.loadTmp('init')
	jumper.setVelocities=True
	mat=O.materials[0]
	yp.reset()
	#O.step()
	#O.tags['id']=mode+'-ROT' #'one-one_'+O.tags['isoTime']+'_'+mode
	O.tags['id']=mode
	#oof=open(O.tags['id']+'.oofem','w').write('\n'.join(eudoxos.oofemTextExport()));
	#oofemDirect(O.tags['id'],negIds=[0],posIds=[1])
	if mode=='viscShear':
		nSteps=10
		yp.plots={'dz1/dist':[]}
		deltaZ=50*epsCrackOnset*dist/2.
		print deltaZ
		tau=1.
		maxima=[]
		for strainRate in [.01,.025,.05,.1,.25,.5,1,2.5,5,10,25,50,100]:
			O.loadTmp('init')
			speed=strainRate*dist
			O.dt=(deltaZ/nSteps)/speed
			key=r'$|F_s|$, $|\dot\varepsilon_T|/\tau_s$=%g'%(strainRate/tau)
			def plotDta():
				iPhys=O.interactions[0,1].phys
				data={'dz1_div_dist':(O.bodies[1].state.pos[2]-p2[2])/dist,key:iPhys.Fs.Length()}
				yp.addData()
			yp.plots['dz1/dist'].append(key)
			plotDataCollector.command='plotDta()'
			jumper.setVelocities=True
			jumper.deltaSe3=(Vector3(0,0,deltaZ/nSteps),Quaternion().IDENTITY)
			mat.plTau=tau; mat.plRateExp=.1; mat.neverDamage=False;
			O.run(nSteps,True)
			#try:
			#	print O.interactions[0,1].phys['cummBetaIter']*1./O.interactions[0,1].phys['cummBetaCount']
			#except ZeroDivisionError: pass
			yp.splitData()
			maxima.append((key,max([x if x==x else -float('inf') for x in yp.data[key]])))
		O.tags['id']=mode
		O.tags['title']='Shear motion, (plTau=%g, plRateExp=%g, dist=%g)'%(tau,mat.plRateExp,dist)
		for a,b in maxima: print a,b
		#eudoxos.oofemDirectExport('/tmp/one-one-shear')
	if mode=='viscDist':
		nSteps=100
		yp.plots={'dx1/dist':[]}
		yp.labels.update({'dx1/dist':r'$\varepsilon_N$'})
		deltaX=10*epsCrackOnset*dist/2.
		tau=1.
		#for speed in [.01,.025,.05,.1,.25,.5,1,2.5,5,10,25,50,100]:
		maxima=[]
		for strainRate in [1e2,1e1,1e0,1e-1,1e-3,1e-6,1e-12]:#[3:6]:
			O.loadTmp('init')
			mat=O.materials[0]
			print 'strainRate =',strainRate
			speed=strainRate*dist
			O.dt=(deltaX/nSteps)/speed
			key=r'$F_N(\dot\varepsilon_N/\tau_d=%g$)'%(strainRate/tau)
			key=re.sub(r'(1e)(-|)0*([0-9]+)',r'10^{\2\3}',key)
			#key='Fn, speed/tau=%g'%(speed/tau)
			def plotDta():
				iPhys=O.interactions[0,1].phys
				yp.addData(**{'dx1/dist':(O.bodies[1].state.pos[0]-p2[0])/dist,key:iPhys.Fn,})
			yp.plots['dx1/dist'].append(key)
			plotDataCollector.command='plotDta()'
			plotDataCollector.initRun=True
			jumper.setVelocities=True
			jumper.deltaSe3=Vector3(deltaX/nSteps,0,0),Quaternion().IDENTITY
			mat.dmgTau=tau; mat.dmgRateExp=.1; mat.neverDamage=False
			O.run(nSteps,True)
			try:
				print O.interactions[0,1].phys.cummBetaIter*1./O.interactions[0,1].phys.cummBetaCount
			except ZeroDivisionError: pass
			yp.splitData()
			maxima.append((key,max([x if x==x else -float('inf') for x in yp.data[key]])))
		fig=yp.plot(noShow=True)
		#fig=yp.plot()
		pylab.ylabel(r'$F_N(\varepsilon,\dot\varepsilon_N/\tau_d)$')
		fig.savefig('cpm-viscDist.pdf')
		O.tags['id']=mode
		O.tags['title']='Distant motion, (dmgTau=%g, dmgExpRate=%g, dist=%g)'%(tau,mat.dmgRateExp,dist)
		for a,b in maxima: print a,b
		#eudoxos.oofemDirectExport('/tmp/one-one')
	if mode=='logPressure':
		yp.plots={'dx1':('Fn',)}
		nSteps=200
		jumper.deltaSe3=Vector3(-dist/nSteps,0,0),Quaternion.IDENTITY
		O.run(int(.8*nSteps),True)
		print yp.data.keys()
		O.tags['title']='Proximal motion till linear strain=%g'%(-float(O.iter)/nSteps)
	elif mode=='distProx':
		mat.neverDamage=False
		mat.relDuctility=1.
		#print mat.dict()
		#yp.plots={'dx1':('Fn',)}
		yp.plots={'epsN':('sigmaN',)}
		nSteps=100
		deltaX=6*epsCrackOnset*dist/2.
		jumper.deltaSe3=Vector3(deltaX/nSteps,0,0),Quaternion().IDENTITY
		O.run(nSteps,True);
		#print O.interactions[0,1].phys.omega
		jumper.deltaSe3=Vector3(-deltaX/nSteps,0,0),Quaternion().IDENTITY
		O.run(int(1.2*nSteps),True);
		fig=yp.plot(noShow=True)
		pylab.annotate(r'$\varepsilon_0$',(mat.epsCrackOnset*1.01,-2e5))
		pylab.annotate(r'$\varepsilon_f$',xytext=(mat.epsCrackOnset*(1+mat.relDuctility),-2e5),xy=(mat.epsCrackOnset,mat.epsCrackOnset*mat.young),arrowprops=dict(arrowstyle="-",linestyle='dashed'))
		#fig.set_dpi(40.)
		#fig.set_size_inches((3,2))
		fig.savefig('cpm-distProx.pdf')
		O.tags['title']='Distant and proximal motion' #Distant motion to eps=10*epsCrackOnset (elongation %g), then proximal motion by 20*epsCrackOnset (until eps=-10*epsCrackOnset)'%(deltaX)
	elif mode=='cycleProxDist':
		mat.neverDamage=False
		#yp.plots={'dx1':('Fn',)}
		yp.plots={'epsN':('sigmaN',),}
		nSteps=100
		cpmLaw.epsSoft=-2*epsCrackOnset
		cpmLaw.relKnSoft=.5
		plotDataCollector.initRun=True
		print mat.dict()
		#O.engines=O.engines+[PeriodicPythonRunner(iterPeriod=1,command='printInteraction(O.interactions[0,1].phys.neverDamage)')]
		O.step()
		for i in range(5):
			deltaX=12*epsCrackOnset*dist/2.
			#if i==4: deltaX*=1.5
			if i%2==1: deltaX*=-1 # proximal motion in odd cycles
			jumper.deltaSe3=Vector3(deltaX/nSteps,0,0),Quaternion().IDENTITY
			O.run(int(nSteps+nSteps*1.0*i),True)
		fig=yp.plot(noShow=True)
		pylab.annotate(r'$\varepsilon_s$',xytext=(1.1*cpmLaw.epsSoft,0),xy=(cpmLaw.epsSoft,cpmLaw.epsSoft*mat.young),arrowprops=dict(arrowstyle='-',linestyle='dashed'),ha='left')
		fig.savefig('cpm-cycleProxDist.pdf',title='Cyclic normal loading of single CPM interaction.')
		O.tags['title']='cycleProxDist'
	elif mode=='yShear':
		#yp.plots={'dy1':('Fs','Fn')}
		yp.plots={'dy1':('epsT','epsN')}
		nSteps=400
		deltaY=50*epsCrackOnset*dist/2.
		jumper.deltaSe3=Vector3(0,deltaY/nSteps,0),Quaternion().IDENTITY
		O.run(nSteps,True)
		#O.pause()
		O.tags['title']='Movement perpendicular to the contact in the +y direction such that eps=200*epsCrackOnset (displacement %g)'
	elif mode=='zRot':
		yp.plots={'zRotArcLen':('Fs','Fn')}
		arcLen=2*epsCrackOnset*dist
		arcAngle=2*pi*r2*(arcLen/r2)
		assert(arcAngle/nSteps>1e-7)
		nSteps=100
		jumper.deltaSe3=Vector3(0,0,0),Quaternion(Vector3(0,0,1),arcAngle/nSteps)
		O0run(nSteps,True)
		#yp.plot()
		O.tags['title']='Rotation by %g, i.e. mutual slip by %g'%(arcAngle,arcLen)
	yp.saveGnuplot(O.tags['id'],title=O.tags['title'])
	print "GRAPH: gnuplot %s.gnuplot"%(O.tags['id'])
quit()