Does jcfpm material have strain softening stage?

Hello,
My question is similar to that mentioned in [1], but we use different materials ([1] is CPM material but mine is JCFPM material).

What I am doing is simulating the triaxial compression test of rock, and using the interaction enlarge factor of 1.5 to consider simulating brittle rock.In my understanding, brittle rocks should change from brittleness to ductility under high confining pressure (i.e. strain softening stage after peak).

I try to set the confining pressure from 20MPa to 100MPa, however, the stress-strain curve is similar to that described in [1], i.e have a linear elastic state and then followed by a sharp drop (the tangent of the stress decrease is almost 90 degree)..

I didn't get the desired strain softening stage. Is it due to material problems?

Thanks for help

----------------------code----------------------
from yade import pack, ymport, plot, utils, export, timing
import numpy as np
import sys

Wy=-20e6
rate=-0.5
damp=0.4
stabilityThreshold=0.001
key='_triax_base_'
young=3000e9
name='JCFPM_triax'
#targetPorosity=0.43
compFricDegree=30
poisson=0.4
OUT=str(Wy)+'_JCFPM_triax'

mn,mx=Vector3(0,0,0),Vector3(0.05,0.05,0.05)
O.materials.append(JCFpmMat(type=1,density=2640,young=young,poisson=poisson,tensileStrength=160e6,cohesion=1600e6,frictionAngle=radians(10),label='sphere'))
O.materials.append(JCFpmMat(type=1,frictionAngle=0,density=0,label='wall'))

walls=aabbWalls([mn,mx],thickness=0,material='wall')
wallIds=O.bodies.append(walls)

O.bodies.append(ymport.text('packing-JCFPM.spheres',scale=1,shift=Vector3(0,0,0),material='sphere',color=(0.6,0.5,0.15)))

triax=TriaxialStressController(
 maxMultiplier=1.+4e8/young,
 finalMaxMultiplier=1.+4e7/young,
 thickness = 0,
 stressMask = 7,
 internalCompaction=True,
)
newton=NewtonIntegrator(damping=damp)

def recorder():
 yade.plot.addData(
 i=O.iter,
 e11=-triax.strain[0],e22=-triax.strain[1],e33=-triax.strain[2],
 s11=-triax.stress(triax.wall_right_id)[0],#0+边界id为right
 s22=-triax.stress(triax.wall_top_id)[1],#1+边界id为top
 s33=-triax.stress(triax.wall_front_id)[2],#2+边界id为front
 numberTc=interactionLaw.nbTensCracks,
 numberSc=interactionLaw.nbShearCracks,
 unb=unbalancedForce()
)
 plot.saveDataTxt(OUT)

def stop_condition():
 extremum=max(abs(s) for s in plot.data['s33'])
 s=abs(plot.data['s33'][-1])
 e=abs(plot.data['e33'][-1])
 #if abs(s) > .5*abs(extremum) or e < 0.005:
 if e < 0.001 :
  return
 if abs(s)/abs(extremum) < 0.9:
  print('Max stress and strain:',extremum,e)
  presentcohesive_count = 0
  for i in O.interactions:
          if hasattr(i.phys, 'isCohesive'):
               if i.phys.isCohesive == True:
                   presentcohesive_count+=1
  print('the number of cohesive bond now is:',presentcohesive_count)
  O.wait()

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.5,label='is2aabb'),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.5,label='ss2sc'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT+'_Crack',label='interactionLaw')]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+key),
 PyRunner(iterPeriod=int(1000),initRun=True,command='recorder()',label='data',dead=0),
 PyRunner(iterPeriod=1000,command='stop_condition()',dead=0),
 VTKRecorder(iterPeriod=500,initRun=True,fileName='triax/JFFPM-',recorders=['spheres','jcfpm','cracks'],Key=OUT+'_Crack',label='vtk',dead=1),
 newton,
]

O.step()
ss2sc.interactionDetectionFactor=-1
is2aabb.aabbEnlargeFactor=-1
cohesiveCount = 0
for i in O.interactions:
 if hasattr(i.phys, 'isCohesive'):
  if i.phys.isCohesive == True:
   cohesiveCount+=1
print('the origin total number of cohesive bond is:',cohesiveCount)

triax.goal1=triax.goal2=triax.goal3=Wy
while 1:
 #global Wy
 O.run(500,1)
 unb=unbalancedForce()
 print('unbalanced force:',unb,'mean stres:',triax.meanStress)
 if unb<stabilityThreshold and abs(Wy-triax.meanStress)/abs(Wy)<0.001:
  break

cohesiveCount = 0
for i in O.interactions:
 if hasattr(i.phys, 'isCohesive'):
  if i.phys.isCohesive == True:
   cohesiveCount+=1
print('the first total number of cohesive bond is:',cohesiveCount)

triax.internalCompaction=False
triax.stressMask=3
triax.goal1=Wy
triax.goal2=Wy
triax.goal3=rate

plot.plots={'e33':('s33',None,'unb'),'i':('numberTc','numberSc',None,'s33')}
plot.plot()
O.run()
--------------------------------------------------------
[1]https://answers.launchpad.net/yade/+question/695734