Perform UCS-TS test

Hi to all,

currently i am trying to perform simulation of UCS and TS as in the paper <<A DEM model for soft and hard rocks: Role of grain interlocking on strength >> of Luc Scholtes and Frédéric-Victor Donzé.

I took the file from trunk/examples/concrete/uniax.py and change the model from CpmMat to JCFpmMat.
Also i took the microparameters from table 2 in the above paper and put it in my script.

Some Remarks i saw in my simulation:
As i put intRadius=1 as in table 2, my Coordination number is 3.37 byt in paper is 6.1
Number of elements: 17412
And Fc = 1.36e06 .
Am i doing something wrong?
 Please have a look in my code :

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division

from yade import plot,pack,timing
import time, sys, os, copy

#import matplotlib
#matplotlib.rc('text',usetex=True)

# default parameters or from table
readParamsFromTable(noTableOk=True, # unknownOk=True,
    density=2640,
 young=230e9,
 poisson=2.5,
 frictionAngle=25,
 tensileStrength=22e6,
 cohesion=220e6, #10*tensileStrength,

 sphereRadius=1.25e-3,
 rRelFuzz=0.34,

 intRadius=1,
 dtSafety=.8,
 damping=0.6,
 strainRateTension=.05,
 strainRateCompression=.5,
 setSpeeds=True,
 # 1=tension, 2=compression (ANDed; 3=both)
 doModes=2,

 # isotropic confinement (should be negative)
 isoPrestress=0,
)

from yade.params.table import *

if 'description' in O.tags.keys(): O.tags['id']=O.tags['id']+O.tags['description']

sample = O.materials.append(JCFpmMat(young=young, poisson=poisson, frictionAngle=radians(frictionAngle), cohesion=cohesion, tensileStrength=tensileStrength,density = density, label='spheres'))

sps=SpherePack()

sp=pack.randomDensePack(pack.inAlignedBox((-25e-3,-25e-3,-62.5e-3),(25e-3,25e-3,62.5e-3)),spheresInCell=2000,radius=sphereRadius,rRelFuzz=rRelFuzz,memoizeDb=False,returnSpherePack=True) #'/tmp/triaxPackCache.sqlite'

sp.toSimulation(color=(1.09,1.0,0.0),material=sample) #color=(1.09,1.0,0.0)

bb=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']

O.dt=dtSafety*PWaveTimeStep()
print 'Timestep',O.dt

mm,mx=[pt[axis] for pt in aabbExtrema()]
coord_25,coord_50,coord_75=mm+.25*(mx-mm),mm+.5*(mx-mm),mm+.75*(mx-mm)
area_25,area_50,area_75=approxSectionArea(coord_25,axis),approxSectionArea(coord_50,axis),approxSectionArea(coord_75,axis)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb')],verletDist=.05*sphereRadius),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc')],
        [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
        [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(Key='cracks',recordCracks=True,cracksFileExist=True,label='interactionLaw')],
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
    NewtonIntegrator(damping=damping,label='damper'),
        VTKRecorder(fileName='3d-vtk-',recorders=['intr','cracks','jcfpm','facets','spheres','colors'],Key='cracks', realPeriod=50),
 PyRunner(virtPeriod=1e-6/strainRateTension,realPeriod=1,command='addPlotData()',label='plotDataCollector',initRun=True),
 PyRunner(realPeriod=4,command='stopIfDamaged()',label='damageChecker'),
]

plot.plots={'eps':('sigma',)}

O.saveTmp('initial');

O.timingEnabled=False

global mode
mode='tension' if doModes & 1 else 'compression'

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 stopIfDamaged():
 global mode
 if O.iter<2 or not plot.data.has_key('sigma'): return # do nothing at the very beginning
 sigma,eps=plot.data['sigma'],plot.data['eps']
 extremum=max(sigma) if (strainer.strainRate>0) else min(sigma)
 minMaxRatio=0.5 if mode=='tension' else 0.5
 if extremum==0: return
 import sys; sys.stdout.flush()
 if abs(sigma[-1]/extremum)<minMaxRatio or abs(strainer.strain)>(5e-3 if isoPrestress==0 else 5e-2):
  if mode=='tension' and doModes & 2: # only if compression is enabled
   mode='compression'
   O.save('uniax-tension.yade.gz')
   print "Saved /tmp/uniax-tension.yade.gz (for use with interaction-histogram.py and uniax-post.py)"
   print "Damaged, switching to compression... "; O.pause()
   # important! initTest must be launched in a separate thread;
   # otherwise O.load would wait for the iteration to finish,
   # but it would wait for initTest to return and deadlock would result
   import thread; thread.start_new_thread(initTest,())
   return
  else:
   print "Damaged, stopping."
   ft,fc=max(sigma),min(sigma)
   print 'Strengths fc=%g, ft=%g, |fc/ft|=%g'%(fc,ft,abs(fc/ft))
   title=O.tags['description'] if 'description' in O.tags.keys() else O.tags['params']
   print 'gnuplot',plot.saveGnuplot(O.tags['id'],title=title)
   print 'Bye.'
   O.pause()
   #sys.exit(0) # results in some threading exception

def addPlotData():
 yade.plot.addData({'t':O.time,'i':O.iter,'eps':strainer.strain,'sigma':strainer.avgStress+isoPrestress,
  'sigma.25':forcesOnCoordPlane(coord_25,axis)[axis]/area_25+isoPrestress,
  'sigma.50':forcesOnCoordPlane(coord_50,axis)[axis]/area_50+isoPrestress,
  'sigma.75':forcesOnCoordPlane(coord_75,axis)[axis]/area_75+isoPrestress,
  })
plot.plot(subPlots=False)
#O.run()
initTest()

Thank You!!