Looking for Direct Shear test script

Hello,

according to the amount of answers, the question provided is not motivating enough to answer.
Please provide more information, for example provide (not only!):
- material
   - real - soil, concrete, ...
   - how you want to simulate it
- 2D / 3D
- loading
   - lateral stress?
   - ...
- what results you are interested in
   - forces
   - failure mode
   - ...
- do you have some ideas / have you tried something already?
   - e.g. boundary conditions - using walls or directly prescribed motion of boundary spheres or ...
   - ...
- what have you searched (googling 'yade "direct shear"' or 'dem "direct shear"' gives plenty of results) and why it did not help?
- ...
- ...

cheers
Jan

Hi,

The following lines may help you build a direct shear test simulation (the script needs other scripts to work properly, i.e., one to generate the sample and one to do the identification of the contact along the joint) . Let us know if that's what you are after.

###

# -*- coding: utf-8 -*-
O=Omega()

from yade import ymport, utils, plot

################## parameters definition
PACKING='global_spheres.xyz' # imported from a text file

intR=1.4
DENS=3000
YOUNG=12e9
FRICT=18
POISSON=0.15
TENS=11e6
COH=11e7
SMOOTH=True
jointNSTIFF=12e9
jointSSTIFF=12e8
jointFRICT=30
jointTENS=0.
jointCOH=0.
jointDIL=0.

DAMPING=0.5
dtCOEFF=0.5
iterMAX=100000
saveVTK=10
recCRACKS=True
Output=PACKING+'_1MPa'

normalSTRESS=1.e6
normalVEL=normalSTRESS/1e8 # optimized for normalVEL=normalSTRESS/1e8?
shearVEL=2*normalVEL # optimized for shearVEL=normalVEL?

###################### Import of the sphere assembly
### material definition
def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson=POISSON,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH,jointNormalStiffness=jointNSTIFF,jointShearStiffness=jointSSTIFF,jointTensileStrength=jointTENS,jointCohesion=jointCOH,jointFrictionAngle=radians(jointFRICT),jointDilationAngle=radians(jointDIL))
def wallMat(): return JCFpmMat(type=0,density=DENS,young=YOUNG,poisson=POISSON,frictionAngle=radians(0))

## copy spheres from the packing into the scene
O.bodies.append(ymport.text(PACKING,scale=1,shift=Vector3(0,0,0),material=sphereMat))

## preprocessing to get dimensions
dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf
## initial surface
S0=X*Z

## spheres factory
R=0
Rmax=0
numSpheres=0.
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   o.shape.color=(0.7,0.5,0.3)
   numSpheres+=1
   R+=o.shape.radius
   if o.shape.radius>Rmax:
     Rmax=o.shape.radius
Rmean=R/numSpheres

#### create the fracture directly here
import gts
ptA = gts.Vertex( xinf-X/4., yinf-Y/4., zinf+Z/2.)
ptB = gts.Vertex( xinf-X/4., ysup+Y/4., zinf+Z/2.)
ptC = gts.Vertex( xsup+X/4., yinf-Y/4., zinf+Z/2.)
ptD = gts.Vertex( xsup+X/4., ysup+Y/4., zinf+Z/2.)
e1 = gts.Edge(ptA,ptB)
e2 = gts.Edge(ptA,ptC)
e3 = gts.Edge(ptC,ptB)
f1 = gts.Face(e1,e2,e3)
e4 = gts.Edge(ptC,ptD)
e5 = gts.Edge(ptD,ptB)
f2 = gts.Face(e4,e5,e3)
s1 = gts.Surface()
s1.add(f1)
s1.add(f2)
facet = gtsSurface2Facets(s1,wire = False,material=wallMat)
O.bodies.append(facet)

execfile('identifBis.py') #you need the script in the same folder (this one is available in the JCFPM examples)

###################### creation of shear box
thickness=Z/100.
oversizeFactor=1.3

### loading platens
O.bodies.append(utils.box(center=(xinf+X/2.,yinf+Y/2.,zinf-thickness+Rmean/5.),extents=(oversizeFactor*X/2,oversizeFactor*Y/2,thickness),material=wallMat,fixed=True))
bottomPlate=O.bodies[-1]
O.bodies.append(utils.box(center=(xinf+X/2.,yinf+Y/2.,zsup+thickness-Rmean/5.),extents=(oversizeFactor*X/2,oversizeFactor*Y/2,thickness),material=wallMat,fixed=True))
topPlate=O.bodies[-1]

###################### Engines
O.engines=[

 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='aabb'),Bo1_Box_Aabb()]),

 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='ss2d3dg'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=SMOOTH,recordCracks=True,label='interactionLaw')]
 ),

 GlobalStiffnessTimeStepper(defaultDt=0.1*PWaveTimeStep(),timestepSafetyCoefficient=dtCOEFF),
 TranslationEngine(ids=[topPlate.id],translationAxis=(0.,0.,-1.),velocity=0.,label='zTranslation'),
 PyRunner(iterPeriod=1,initRun=True,command='servoController()',label='servo'),
 NewtonIntegrator(damping=DAMPING,gravity=(0,0,0),label='damper'),
 PyRunner(iterPeriod=100,initRun=True,command='dataCollector()'),
 VTKRecorder(iterPeriod=int(iterMAX/saveVTK),initRun=True,fileName=Output+'-',recorders=['spheres','stress','velocity','jcfpm','cracks','boxes'])

]

###################### Engines definition ( servoController() and dataCollector() )
shearing=False
sigmaN=0
tau=0
Fs1=0
Fs2=0
Xdispl=0
px0=0
Zdispl=0
pz0=topPlate.state.pos[2]
prevTranslation=0
n=0

def servoController():
    global px0, pz0, sigmaN, n, Fn1, Fn2, shearing, butee, piston1, piston2
    Fn1=abs(O.forces.f(topPlate.id)[2])
    Fn2=abs(O.forces.f(bottomPlate.id)[2])
    sigmaN=Fn1/(S0-2*Xdispl*Y) # necessary?
    #print yTranslation.velocity, sigmaN
    if shearing==False:
      if zTranslation.velocity<normalVEL:
 zTranslation.velocity+=normalVEL/100.
      if sigmaN>(0.95*normalSTRESS):
 zTranslation.velocity=10*normalVEL*((normalSTRESS-sigmaN)/normalSTRESS) # not good enough because not general (coeff needs to be adjusted)
    if shearing==False and abs((normalSTRESS-sigmaN)/normalSTRESS)<0.001 and unbalancedForce()<0.005:
      zTranslation.velocity=0
      print 'stress on joint plane = ', utils.forcesOnPlane((X/2,Y/2,Z/2),(0,0,1))/S0
      ### top box
      O.bodies.append(utils.box(center=(xinf-thickness+Rmean/5.,(yinf+0.5*Y),(zinf+0.75*Z)+3*Rmean)
    ,extents=(thickness,oversizeFactor*Y/2.,Z/4.),material=wallMat,fixed=True))
      butee=O.bodies[-1]
      O.bodies.append(utils.box(center=(xsup+thickness-Rmean/5.,(yinf+0.5*Y),(zinf+0.75*Z)+3*Rmean),
    extents=(thickness,oversizeFactor*Z/2.,Z/4.),material=wallMat,fixed=True))
      O.bodies.append(utils.box(center=((xinf+0.5*X),yinf-thickness+Rmean/5.,(zinf+0.75*Z)+3*Rmean),
    extents=(oversizeFactor*X/2.,thickness,Z/4.),material=wallMat,fixed=True,wire=True))
      O.bodies.append(utils.box(center=((xinf+0.5*X),ysup+thickness-Rmean/5.,(zinf+0.75*Z)+3*Rmean),
    extents=(oversizeFactor*X/2.,thickness,Z/4.),material=wallMat,fixed=True,wire=True))
      ### bottom box
      O.bodies.append(utils.box(center=(xsup+thickness-Rmean/5.,(yinf+0.5*Y),(zinf+0.25*Z)-3*Rmean),
    extents=(thickness,oversizeFactor*Y/2.,Z/4.),material=wallMat,fixed=True))
      piston1=O.bodies[-1]
      O.bodies.append(utils.box(center=(xinf-thickness+Rmean/5.,(yinf+0.5*Y),(zinf+0.25*Z)-3*Rmean),
    extents=(thickness,oversizeFactor*Y/2.,Z/4.),material=wallMat,fixed=True))
      piston2=O.bodies[-1]
      O.bodies.append(utils.box(center=((xinf+0.5*X),yinf-thickness+Rmean/5.,(zinf+0.25*Z)-3*Rmean),
    extents=(oversizeFactor*X/2.,thickness,Z/4),material=wallMat,fixed=True,wire=True))
      O.bodies.append(utils.box(center=((xinf+0.5*X),ysup+thickness-Rmean/5.,(zinf+0.25*Z)-3*Rmean),
    extents=(oversizeFactor*X/2.,thickness,Z/4.),material=wallMat,fixed=True,wire=True))
      ### initialisation
      O.engines=O.engines[:5]+ [TranslationEngine(ids=[piston1.id,piston2.id],translationAxis=(-1.,0.,0.),velocity=0.)]+O.engines[5:]
      px0=piston1.state.pos[0]
      pz0=topPlate.state.pos[2]
      shearing=True
      print 'shearing now! || iteration=', O.iter
    if shearing==True:
      zTranslation.velocity=0.
      ## Removes collider (lattice like) -> not really helpful in this particular case
      #O.engines=O.engines[:1]+O.engines[2:]
      zTranslation.velocity=50*normalVEL*((normalSTRESS-sigmaN)/normalSTRESS) # not good enough because not general (coeff needs to be adjusted)
      if O.engines[5].velocity<shearVEL:
 O.engines[5].velocity+=(shearVEL/1000)

def dataCollector():
    global Xdispl, Zdispl, tau, Fs1, Fs2
    if shearing:
      Fs1=abs(O.forces.f(butee.id)[0])
      Fs2=abs(O.forces.f(piston1.id)[0])
      Xdispl=abs(piston1.state.pos[0]-px0)
    tau=Fs1/(S0-2*Xdispl*Y) # necessary?
    Zdispl=abs(topPlate.state.pos[2]-pz0)
    yade.plot.addData({'t':O.time,'i':O.iter,'Xdispl':Xdispl,'sigmaN':sigmaN,'tau':tau,'Zdispl':Zdispl,'unbF':utils.unbalancedForce(),'tc':interactionLaw.nbTensCracks,'sc':interactionLaw.nbShearCracks})
    plot.saveDataTxt(Output)
    #plot.saveGnuplot(Output)

################# to manage interaction detection factor during the first timestep
O.step();
################# initializes the interaction detection factor to its default value (new contacts will be created between strictly contacting particles only)
ss2d3dg.interactionDetectionFactor=-1.
aabb.aabbEnlargeFactor=-1.

################# simulation starts here
O.run(int(iterMAX))

Hello,
Sorry for the lack of infos. I’m new to the yade community.

Luc, thanks for your support. I’ve used a script that you posted on this forum:

###

from yade import pack,plot,export
import math

sp=pack.SpherePack()
O.periodic=True

# dimensions of sample (fixed by particle size such as L/D~15)
RADIUS=0.05
length=15*(2*RADIUS)
height=length
width=length
thickness=RADIUS

# friction angles
compFRIC=10. # during compaction (controls porosity)
FRIC=30. # during shear

# boundary conditions
PI=1.e5 # sample preparation: pressure applied for the isotropic compaction
SN=5.e6 # normal stress
RATE=0.1 # shear velocity (top plate)

# simulation control
DAMPSHEAR=0.
ITER=2e5
VTK=20
OUT='TEST'

####

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)

O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=radians(0.),normalCohesion=1e100,shearCohesion=1e100,label='boxMat'))
O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=radians(compFRIC),normalCohesion=0,shearCohesion=0,label='sphereMat'))

upBox = utils.box(center=(0,2*height+thickness/2.0,0),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,2*width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(0,height-thickness/2.0,0),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,2*width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
O.bodies.append([upBox,lowBox])

sp.makeCloud((0,height+1.5*RADIUS,0),(length,2*height-1.5*RADIUS,width),rMean=RADIUS,rRelFuzz=0.2,periodic=True)
O.bodies.append([utils.sphere(s[0],s[1],color=(0,0,1),material='sphereMat') for s in sp])

effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol

#print 'volRatio=',volRatio

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep())
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(damping=0.3,label='newton')
 ,PyRunner(command='dataRecorder()',iterPeriod=10,label='recData',dead=True)
 ,VTKRecorder(fileName=OUT+'.',iterPeriod=1,skipNondynamic=1,recorders=['spheres','colors','velocity','bstresses'],label='saveSolid',dead=True)
]

def dataRecorder():
 h=vol=vol_s=nb_s=0.
 h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
 vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
 contactStress=getStress(vol)
 for o in O.bodies:
   if isinstance(o.shape,Sphere) and o.shape.color[0]!=1:
    nb_s+=1
    vol_s += 4.*pi/3.*(o.shape.radius)**3
 n = 1-vol_s/vol
 nbFrictCont=0.
 for i in O.interactions:
  if i.isReal and i.phys.cohesionBroken:
   nbFrictCont+=1
 plot.addData(
  iter=O.iter
  ,stress_upWall0=abs(O.forces.f(0)[0]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])),stress_upWall1=abs(O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])),stress_upWall2=abs(O.forces.f(0)[2]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]))
  ,contactStress00=(contactStress[0,0]),contactStress01=(contactStress[0,1]),contactStress02=(contactStress[0,2]),contactStress10=(contactStress[1,0]),contactStress11=(contactStress[1,1]),contactStress12=(contactStress[1,2]),contactStress20=(contactStress[2,0]),contactStress21=(contactStress[2,1]),contactStress22=(contactStress[2,2])
  ,xW=O.bodies[0].state.pos[0]
  ,height=h
  ,volume=vol
  ,porosity=n
  ,k=2.0*nbFrictCont/nb_s
  ,unbF=unbalancedForce()
 )

def triaxDone():
 global phase
 volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])
 h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
 vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
 contactStress=getStress(vol)
 vol_s=Rmean=Rmax=nbSph=0
 Rmin=1e6
 for o in O.bodies:
  if isinstance(o.shape,Sphere):
   nbSph+=1
   Rmean+=o.shape.radius
   if o.shape.radius>Rmax: Rmax=o.shape.radius
   if o.shape.radius<Rmin: Rmin=o.shape.radius
   vol_s += 4.*pi/3.*(o.shape.radius)**3
 Rmean=Rmean/nbSph
 n = 1-vol_s/vol
 print 'DONE! iter=',O.iter,'| sample generated: nb spheres',nbSph,', Rmean=',Rmean,', Rratio=',Rmax/Rmin,', porosity=',n
 print 'Changing contact properties now'
 tt=TriaxialCompressionEngine()
 tt.setContactProperties(FRIC)
 triax.dead=True
 O.pause()

#### Initialization
print 'SAMPLE PREPARATION!'

#recData.dead=False # uncomment to record what is happening during stress initialization
O.run(1000000,1)
saveSolid.dead=False
#saveSolid.fileName=OUT+'_isoConfined.'
O.step()
saveSolid.dead=True
O.save(OUT+'_initialState.yade')

print 'Normal stress (platen) = ',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
print 'Normal stress (contacts) = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])[1,1]

#### Applying normal stress
print 'NORMAL LOADING! iter=',O.iter

stage=0
stiff=fnPlaten=currentSN=0.
def servo():
 global stage,stiff,fnPlaten,currentSN
 if stage==0:
  currentSN=O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
  unbF=unbalancedForce()
  #print 'SN=',SN,'| current SN = ',currentSN,' | unbF=',unbF
  boundaryVel=copysign(min(0.1,abs(0.5*(currentSN-SN))),currentSN-SN)
  O.bodies[0].state.vel[1]=boundaryVel
  if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.001 ):
   stage+=1
   fnPlaten=O.forces.f(0)[1]
   print 'Normal stress =',currentSN,' | unbF=',unbF
   ## the following computes the stiffness of the plate (used for stress control of the top plate)
   for i in O.interactions.withBody(O.bodies[0].id):
    stiff+=i.phys.kn
   print 'DONE! iter=',O.iter
   O.pause()
 if stage==1:
  fnDesired=SN*(O.cell.hSize[0,0]*O.cell.hSize[2,2])
  #boundaryVel=copysign(abs(0.5*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt),O.forces.f(0)[1]-fnDesired)
  boundaryVel=copysign(min(RATE,abs(0.333*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt)),O.forces.f(0)[1]-fnDesired)
  O.bodies[0].state.vel[1]=boundaryVel

O.engines = O.engines[:4]+[PyRunner(command='servo()',iterPeriod=1,label='servo')]+O.engines[4:]

O.run(1000000,1)
print 'STABILIZING! iter=',O.iter
O.run(1000,1)
# coloring particles to see vertical stripes in material
dxi=4*(2*RADIUS) # can be a function of cell length dxi=O.cell.hSize[0,0]/5.
n=int(length/dxi)
xmin=1e6
for i in range(0,n):
 for o in O.bodies:
  if o.id>1:
   if o.state.pos[0]<xmin: xmin=o.state.pos[0]
   if (o.state.pos[0]>=i*dxi) and (o.state.pos[0]<((i+0.5)*dxi)):
    o.shape.color[1]=1
for o in O.bodies:
 if (o.state.pos[0]>(xmin+0.5*O.cell.hSize[0,0]-RADIUS)) and (o.state.pos[0]<(xmin+0.5*O.cell.hSize[0,0]+3*RADIUS)):
  o.shape.color[2]=0

saveSolid.dead=False
#saveSolid.fileName=OUT+'_normalLoaded.'
O.step()
saveSolid.dead=True
O.save(OUT+'_normallyLoaded.yade')
if recData.dead==False: plot.saveDataTxt(OUT)

print 'Normal stress (platen) = ',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
print 'Normal stress (contacts) = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])[1,1]

#### preparing for shearing
Gl1_Sphere.stripes=1
print 'Gluing spheres to boundary walls'
gluedSpheres=[]

## gluing particles in contact with the walls
for i in O.interactions:
  if i.isReal:
   if isinstance(O.bodies[i.id1].shape,Box):
    O.bodies[i.id2].shape.color[0]=1
    gluedSpheres += [O.bodies[i.id2]]
   if isinstance(O.bodies[i.id2].shape,Box):
    O.bodies[i.id1].shape.color[0]=1
    gluedSpheres += [O.bodies[i.id1]]

for i in O.interactions:
  if i.isReal and ( O.bodies[i.id1].shape.color[0]==1 and O.bodies[i.id2].shape.color[0]==1 ):
   O.bodies[i.id1].mat.normalCohesion=O.bodies[i.id1].mat.normalCohesion
   O.bodies[i.id2].mat.normalCohesion=O.bodies[i.id1].mat.normalCohesion
   O.bodies[i.id1].mat.shearCohesion=O.bodies[i.id1].mat.shearCohesion
   O.bodies[i.id2].mat.shearCohesion=O.bodies[i.id1].mat.shearCohesion
   i.phys.initCohesion=True

print 'nb of glued spheres=',len(gluedSpheres)

#### shearing
print 'SHEARING! iter=',O.iter

saveSolid.dead=False
O.step()
saveSolid.dead=True
O.save(OUT+'_shearInit.yade')

recData.dead=False
newton.damping=DAMPSHEAR
saveSolid.dead=False
saveSolid.iterPeriod=int(ITER/VTK)
shearVel=0
iterShear=O.iter
while 1:
 O.run(int(10),1)
 if shearVel<RATE:
  shearVel+=(RATE/100.)
 # only top wall moves
 O.bodies[0].state.vel[0]=shearVel
 ## top and bottom walls move
 #O.bodies[0].state.vel[0]=0.5*shearVel
 #O.bodies[1].state.vel[0]=-0.5*shearVel
 if ( O.iter >= (int(iterShear+ITER)) ):
  print 'iter=',O.iter,' -> FINISHED!'
  plot.saveDataTxt(OUT)
  O.save(OUT+'_sheared.yade')
  sys.exit(0)

Is it suitable for a model of a direct shear test?
What about units of measurement? What’s the unit of RADIUS for example.

Thanks, and sorry for my poor English.

Amedeo

> What about units of measurement? What’s the unit of RADIUS for example.

Yade has no units, just numbers [1,2,3,4].
If you use basic SI units for all inputs, you get all output in basic SI units, too.
Using only basic SI units is a good practice.

cheers
Jan

[1] https://answers.launchpad.net/yade/+question/691104
[2] https://answers.launchpad.net/yade/+question/690022
[3] https://answers.launchpad.net/yade/+question/687871
[4] https://answers.launchpad.net/yade/+question/672270