triaxial test on soft rock

I want to model a triaxial test on soft rock using Yade. I generate a sample using Hertz-Mindlin law, and the sample contains coarse and fine grains (specified PSD). Then in triaxial shearing stage, I change the material and contact law into JCF material and contact law. But failed. This way is impossible in Yade? Could you give me some advice? Thanks!

from __future__ import print_function
from yade import pack,plot,os,timing
from yade import export, ymport
import matplotlib;matplotlib.rc('axes',grid=True)
import pylab
import numpy as np

O.load('isocompact200kpa2.yade.bz2')

#O.reset()

for b in O.bodies:
       if isinstance(b.shape, Sphere):
            O.bodies.erase(b.id) #delete grains firstly and they will be generate later

density=2600
frictionangle=0
poisson=0.15
young=30e9
damp=0.25
number=10000

walldensity=0
wallfrictionangle=0
wallpoisson=0.5
wallyoung=30000e9

mn=Vector3(0.,0.,0.)
mx=Vector3(2,2,2)

conpress=-200000
rate=0.1

#mat=O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(frictionangle),density=density,label='spheres'))
wallmat=O.materials.append(JCFpmMat(type=1,frictionAngle=0,density=0,label='walls'))
idRockTest = grainmat=O.materials.append(JCFpmMat(type=1,young=30e9,density=2500.0,poisson=0.1,frictionAngle=np.radians(18.0),tensileStrength=1e6,cohesion=1e6,label='Rock'))

O.bodies[0].material=O.materials[wallmat]
O.bodies[1].material=O.materials[wallmat]
O.bodies[2].material=O.materials[wallmat]
O.bodies[3].material=O.materials[wallmat]
O.bodies[4].material=O.materials[wallmat]
O.bodies[5].material=O.materials[wallmat]
#wallIds=O.bodies.append(aabbWalls([mn,mx],thickness=0.001,material='walls'))
#Material constants
sp3 = ymport.text('/home/hezihao/desk/bond/cloud.txt')
for s in range(0,len(sp3)):
        #sp2[s].shape.color = (1.0,1.0,0)
        sp3[s].state.material=grainmat
        O.bodies.append(sp3[s])

for b in O.bodies:
    if isinstance(b.shape,Sphere):
       r=b.shape.radius
       if r>0.03 and r<0.50:
          b.shape.color=(1.0,0/255.,51/255.)
          #b.mat=O.materials[grainmat]
       if r>0.0 and r<0.03:
          b.shape.color=(1.0,1.0,0)
          #b.mat=O.materials[grainmat]
#for i in range(5,number+1):
    #O.bodies[i].material=O.materials[grainmat]
#FrictionAngle = 35

Damp = 0.25

#Confining variables
ConfPress = -2.0e5

#Loading control
LoadRate = -0.01

#folders for post-processing
#try:
    #os.mkdir('post')
#except:
    #pass

#try:
    #os.mkdir('pdata')
#except:
    #pass

#try:
    #os.mkdir('cdata')
#except:
    #pass

#restart

########################################
#import necessary packages
from yade import pack,plot,os,timing

#def Output():
    #e22=-triax.strain[1]
    #particle info
    #f=open('./pdata/pInfo'+'{:.5f}'.format(e22), 'w')
    #f.write('Box position: \n')
    #f.write('leftwall:%.5f %.5f %.5f\n'%(O.bodies[0].state.pos[0],O.bodies[0].state.pos[1],O.bodies[0].state.pos[2]))
    #f.write('rightwall:%.5f %.5f %.5f\n'%(O.bodies[1].state.pos[0],O.bodies[1].state.pos[1],O.bodies[1].state.pos[2]))
    #f.write('topwall:%.5f %.5f %.5f\n'%(O.bodies[3].state.pos[0],O.bodies[3].state.pos[1],O.bodies[3].state.pos[2]))
    #f.write('bottomwall:%.5f %.5f %.5f\n'%(O.bodies[2].state.pos[0],O.bodies[2].state.pos[1],O.bodies[2].state.pos[2]))
    #f.write('frontwall:%.5f %.5f %.5f\n'%(O.bodies[5].state.pos[0],O.bodies[5].state.pos[1],O.bodies[5].state.pos[2]))
    #f.write('backwall:%.5f %.5f %.5f\n'%(O.bodies[4].state.pos[0],O.bodies[4].state.pos[1],O.bodies[4].state.pos[2]))

    #f.write('ID x y z radius disx disy disz rotx roty rotz \n')
    #for b in O.bodies:
       #if isinstance(b.shape, Sphere):
          #pos=b.state.pos
          #rad=b.shape.radius
          #displ=b.state.displ()
          #rot=b.state.rot()
          #f.write('%-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g\n'%(b.id,pos[0],pos[1],pos[2],rad,displ[0],displ[1],displ[2],rot[0],rot[1],rot[2]))

    #f.close()

    #contact info
    #g = open('./cdata/cInfo'+'{:.5f}'.format(e22), 'w')
    #print ('Contact information at Iter %d' % O.iter)
    #g.write('ctype id1 id2 nfx nfy nfz tfx tfy tfz \n')

    #for k in O.interactions:
       #if isinstance(O.bodies[k.id1].shape,Sphere) and isinstance(O.bodies[k.id2].shape,Box):
           #g.write('01:%.5f %.5f %.5f %.5f %.5f %.5f %.5f %.5f\n'%(k.id1,k.id2,k.phys.normalForce[0], k.phys.normalForce[1],k.phys.normalForce[2],k.phys.shearForce[0],k.phys.shearForce[1],k.phys.shearForce[2]))
       #elif isinstance(O.bodies[k.id1].shape,Box) and isinstance(O.bodies[k.id2].shape,Sphere):
           #g.write ('10:%.5f %.5f %.5f %.5f %.5f %.5f %.5f %.5f\n'%(k.id1,k.id2,k.phys.normalForce[0], k.phys.normalForce[1],k.phys.normalForce[2],k.phys.shearForce[0],k.phys.shearForce[1],k.phys.shearForce[2]))
       #elif isinstance(O.bodies[k.id1].shape,Sphere) and isinstance(O.bodies[k.id2].shape,Sphere):
           #g.write ('00:%.5f %.5f %.5f %.5f %.5f %.5f %.5f %.5f\n'%(k.id1,k.id2,k.phys.normalForce[0], k.phys.normalForce[1],k.phys.normalForce[2],k.phys.shearForce[0],k.phys.shearForce[1],k.phys.shearForce[2]))

    #g.close()

#Output()
#################################
#####Defining triaxil engines####
#################################

triax=TriaxialStressController(
    #define wall ids
    #wall_bottom_id = wallIds[2],
    #wall_top_id = wallIds[3],
    #wall_left_id = wallIds[0],
    #wall_right_id = wallIds[1],
    #wall_back_id = wallIds[4],
    #wall_front_id = wallIds[5],
    thickness = 0.001,
    #maxMultiplier = 1.002,
    internalCompaction = False, # If true the confining pressure is generated by growing particles
    max_vel = 1.0,
    stressMask = 5,
    goal1 = ConfPress,
    goal2 = LoadRate,
    goal3 = ConfPress,

)

ini_e22a = -triax.strain[1]
ini_e22b = -triax.strain[1]

O.usesTimeStepper=True

#O.dt=2e-6

newton=NewtonIntegrator(damping=Damp)

#setContactFriction(radians(FrictionAngle))

O.trackEnergy=True

O.timingEnabled=True

###engine
O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.3),Bo1_Box_Aabb()]),
    InteractionLoop(
     [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.3),Ig2_Box_Sphere_ScGeom()],
     [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
     [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(label='interactionLaw')]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    triax,
    newton,
    PyRunner(realPeriod=2,command='endCheck()',label='check'),
    PyRunner(command='addPlotData()',iterPeriod=100,label='record'),
    #VTKRecorder(iterPeriod=100,recorders=['all'],fileName='./post/p1-'),
    #TriaxialStateRecorder(iterPeriod=10000,file='WallStresses3')
]
# Simulation stop conditions defination
def endCheck():
    unb=unbalancedForce()
    e22=-triax.strain[1]
    if abs(e22-0.40)<0.001:
       O.pause()
       O.save('e0.05.yade.bz2')
# collect history of data
def addPlotData():

    global ini_e22a
    global ini_e22b

    e22=-triax.strain[1]
    unb = unbalancedForce()
    mStress = -(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3.0
    volume = (O.bodies[1].state.pos[0]-O.bodies[0].state.pos[0])*(O.bodies[3].state.pos[1]-O.bodies[2].state.pos[1])*(O.bodies[5].state.pos[2]-O.bodies[4].state.pos[2])
    Porosity=1-sum(4*pi*b.shape.radius**3/3 for b in O.bodies if isinstance(b.shape,Sphere))/volume

    plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
     ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
   s11=-triax.stress(triax.wall_right_id)[0],
   s22=-triax.stress(triax.wall_top_id)[1],
   s33=-triax.stress(triax.wall_front_id)[2],
                 ub=unbalancedForce(),
                 dstress=(-triax.stress(triax.wall_top_id)[1])-(-triax.stress(triax.wall_right_id)[0]),
                 disx=triax.width,
                 disy=triax.height,
                 disz=triax.depth,
   i=O.iter,
                 porosity=Porosity,
                 cnumber=avgNumInteractions(),
                 ke=utils.kineticEnergy(),
                 totale=O.energy.total()
                 )

    #if (e22 - ini_e22a) > 0.001 and e22 < 0.05:
    # O.save('save'+'{:.4f}'.format(e22)+'yade.bz2')
    # ini_e22a = e22

    if (e22 - ini_e22b) > 0.00005:
        print ('unbalanced force: %f, mean stress: %f, s11: %f, s22: %f, s33: %f, coordination number: %f, porosity: %f' %(unb,mStress,-triax.stress(triax.wall_right_id)[0],-triax.stress(triax.wall_top_id)[1],-triax.stress(triax.wall_front_id)[2],avgNumInteractions(),Porosity))
        #Output()
        ini_e22b = e22
        plot.saveDataTxt('loadinglog3.txt.bz2')
    #plot.saveGnuplot('plotScript')

##################################################
######Defining output for postprocessing##########
##################################################

###display
#yade.qt.Controller(),yade.qt.View()

### declare what is to plot. "None" is for separating y and y2 axis
#plot.plots={'i':('ub',),'i ':('e11','e22','e33',),' i':('s11','s22','s33',),'e22':'dstress',' e22':'ev'}
#plot.plots={'i':('ub',),'i ':('s11','s22','s33'),' i':('e11','e22','e33')}

### the traditional triaxial curves would be more like this:
plot.plots={'e22':('s22',None,'ev')}

plot.plot()
#O.run()#; O.wait()