# yade-users team mailing list archive

## Re: [Question #680144]: plot the deviatoric part of TriaxialStressController simulation

```Question #680144 on Yade changed:

Alireza Sadeghi is still having a problem:
Hello Jerome,

Thank you for your response and help. I will use  https://gitlab.com/yade-dev/trunk. Thank you.
about my question, I used plot.resetData() but it did not solve my problem. The simulation of triaxial test has two parts. The particles compact with isometric pressure and in the second step we have deviatoric compaction. The problem is, strain in the first step is added to the strain in the second step. but I want to draw strain-Stress diagram just for the second step of triaxial test simulation. my code is in the below of the mail. Thank you very much for your help.

Best Regards

Alireza

#Triaxial compresion for finding a proper size for RVE (stress control
boundary condition)

from datetime import datetime

qtr=qt.Renderer()
qtr.bgColor=(1,1,1)

#==============================================================
#================= define the materials =======================
#==============================================================

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion=
1e20, isCohesive= True, young=6.81e8, density=1377.5, poisson=0.3,
frictionAngle= 0.31, fragile=False, label='Coke'))

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= False, young=4e9,
density=1523.6, poisson=0.3, frictionAngle= 0.0, fragile=False, label='wall'))

#===============================================================
#=================== define walls ============================
#===============================================================

walls=aabbWalls([(0,0,0),(2e-1,2e-1,2e-1)],thickness=0.0003,oversizeFactor=1.0,material='wall')
wallIds=O.bodies.append(walls)

#===============================================================
#=================== define packing ============================
#===============================================================

nums=['t']

mats=['Coke']

coke=(5e-3,2000)

nums=pack.SpherePack()

nums.makeCloud((0,0,0),(2e-1,2e-1,2e-1),rMean=coke[0],rRelFuzz=1e-4,num=coke[1])

O.bodies.append([utils.sphere(c,r,material=mats[0],color=(0,0,1)) for
c,r in nums])

#===============================================================
#=================== define Engine =============================
#===============================================================

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()]),
InteractionLoop([Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
),
#	TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
#        VTKRecorder(fileName='3d-vtk-',recorders=['all'],iterPeriod=100),
#        qt.SnapshotEngine(fileBase='3d-',iterPeriod=200,label='snapshot'),

NewtonIntegrator(damping=0.4,gravity=[0,0,0])
]
O.dt=2e-6

#=================================================================
#=================  APPLYING CONFINING PRESSURE   ================
#=================================================================
stabilityThreshold=0.2
sigmaIso=-1e5

triax=TriaxialStressController(

maxMultiplier=1.000,
finalMaxMultiplier=1.000,
thickness = 0,
internalCompaction=False,
)

O.engines=O.engines+[triax]

triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=sigmaIso
triax.wall_back_activated=True

while 1:
O.run(100, True)
unb=unbalancedForce()
meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
print 'mean Stress:',triax.meanStress,'porosity:', triax.porosity,'meanS:',unb
if triax.porosity < 0.497 and ((sigmaIso-triax.meanStress)/sigmaIso) < 0.001:
print 'Isotropic strain1:',-triax.strain[0], 'Isotropic strain 2:',-triax.strain[1], 'Isotropic strain 3:',-triax.strain[2]
break

print "###      Isotropic state saved      ###"

#====================================================================
#====================================================================
plot.reset()
triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=-10
#dataOld = plot.data
#plot.saveDataTxt('isotropic_part_deformableRVE.txt')
plot.resetData()
#plot.plot()

O.saveTmp()

#=================================================================
#======================= data collecting =========================
#=================================================================

O.engines=O.engines+[PyRunner(iterPeriod=16000,command='finishIt()',label='calmEngine')]

def finishIt():
#        if abs(-triax.strain[2]) > 0.2:
#           makeVideo(snapshot.snapshots,'3d.mpeg',fps=10,bps=10000)
print O.iter
print 'time is ', O.time
print 'porosity:',triax.porosity
print 'strain 1:',-triax.strain[0], 'strain 2:',-triax.strain[1], 'strain 3:',-triax.strain[2]

def history():
ev5=(triax.strain[0]+triax.strain[1]+triax.strain[2]),
sxx=-triax.stress(triax.wall_right_id)[0],
syy=-triax.stress(triax.wall_top_id)[1],
szz=-triax.stress(triax.wall_front_id)[2],
q=-(triax.stress(triax.wall_front_id)[2]-triax.stress(triax.wall_right_id)[0]),
p=triax.meanStress,
R=3*(triax.stress(triax.wall_front_id)[2]-triax.stress(triax.wall_right_id)[0])/triax.meanStress,
por=triax.porosity,
i=O.iter),

if 1:
O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]

else:
O.engines[4]=PyRunner(iterPeriod=20,command='history()',label='recorder')

O.run(50000,True)

plot.plots={('ezz'):('ev5')}
plot.plot()

--