1st state

# coding: utf-8
# loose specimen
from yade import pack
import os
#from utils import *

nRead=utils.readParamsFromTable(
    num_spheres=10000,# number of spheres
    compFricDegree = 32, # contact friction during the confining phase the final porosity
    unknownOk=True,
    isoForce=20000,
    conStress=100000
)
from yade.params import table

num_spheres=table.num_spheres
mn,mx=Vector3(-0.1,-0.1,-0.1),Vector3(0.1,0.1,0.1)
thick = 0.1
compFricDegree = table.compFricDegree
finalFricDegree=35
#targetPorosity=0.382
rate=0.05
damp=0.06
stabilityThreshold=0.001 # leave it too low can have a super stable state, but will eat onion for waiting so long
#stopStep=160000 # pas de calcul on choisit pour faire calculer travail du second-ordre
key='_probing_envelope_' # just used for adding name to the output, kind of boring

## create material #0, which will be used as default
O.materials.append(FrictMat(young=356e6,poisson=0.42,frictionAngle=radians(compFricDegree),density=3000,label='spheres'))
O.materials.append(FrictMat(young=356e6,poisson=0.5,frictionAngle=0,density=0,label='walls'))

## create walls around the packing
walls=utils.aabbWalls([mn,mx],thickness=thick,oversizeFactor=5,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
#sp.makeCloud(mn,mx,-1,0,num_spheres,False, 0.95,psdSizes,psdCumm,False,seed=1)
sp.makeCloud(mn,mx,0.001,0.65,num_spheres,False,0.8,seed=1)
#sp.makeCloud(mn,mx,-1,0,-1,False, 0.95,psdSizes,psdCumm,False,seed=1)
sp.toSimulation(material='spheres')

volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
mean_rad = pow(0.09*volume/num_spheres,0.3333)

#clumps=False
#if clumps:
#   c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])
#   sp.makeClumpCloud((-0.24,0.24,-0.24),(0.24,0.24,0.24),[c1],periodic=False)
#   O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])
#   standalone,clumps=sp.getClumps()
#   for clump in clumps:
#       O.bodies.clump(clump)
#       for i in clump[1:]: O.bodies[i].shape.color=O.bodies[clump[0]].shape.color
#   #sp.toSimulation()
#else:
#   O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])

O.dt=.5*utils.PWaveTimeStep() # initial timestep, to not explode right away
O.usesTimeStepper=True

triax=TriaxialStressController(
    maxMultiplier=1.001,
    finalMaxMultiplier=1.0001,
    thickness = thick,
    radiusControlInterval=10,
    stressMask = 7,
    max_vel=0.01,
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()],verletDist=-mean_rad*0.06),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack(label='law')]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
    triax,
    TriaxialStateRecorder(iterPeriod=50,file='WallStresses'+key,addIterNum=False,initRun=False), # 20 is too small, will make big file size
    newton
]

#O.save('initial'+key+'.xml')
#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=1
#if nRead==0: yade.qt.Controller(), yade.qt.View()
print 'Number of elements: ', len(O.bodies)
print 'Box Volume: ',  triax.boxVolume
print 'Box Volume calculated: ', volume

# Applying confining force


#########################################
# Phase 1: ISOTROPIC GENERATOR OF 20kPa #
#########################################
while 1:
  triax.internalCompaction=True # Growing the particles is what we use in this step
  setContactFriction=radians(compFricDegree)
  triax.stressmask=7
  triax.goal1=table.isoForce
  triax.goal2=table.isoForce
  triax.goal3=table.isoForce
  O.run(1000,True)
  #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  unb=unbalancedForce()
  #average stress
  #note: triax.stress(k) returns a stress vector, so we need to keep only the normal component
  meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
  print 'unbalanced force:',unb,' mean stress: ',meanS
#  print 'void ratio=',triax.porosity/(1-triax.porosity), 'porosity=', triax.porosity
  print 'mean stress engine', triax.meanStress, 'kineticE=', utils.kineticEnergy()
  print '-----State_01: Isotropic compression 20kPa--(^_^)---'
  if unb<stabilityThreshold and abs(meanS-table.isoForce)/table.isoForce<0.001:
    break

O.save('confinedState'+key+'.xml') # remember this impotant part, will load is later
print "##   Isotropic state saved (20 kPa)  ##"
print 'current porosity=',triax.porosity
print 'current void ratio=',triax.porosity/(1-triax.porosity)

###### porosity reaching #####
#import sys #this is only for the flush() below
#while triax.porosity>targetPorosity:
#   # we decrease friction value and apply it to all the bodies and contacts
#   compFricDegree = 0.95*compFricDegree
#   setContactFriction(radians(compFricDegree))
#   print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
#   sys.stdout.flush()
#   # while we run steps, triax will tend to grow particles as the packing
#   # keeps shrinking as a consequence of decreasing friction. Consequently
#   # porosity will decrease
#   O.run(500,1)

#O.save('compactedState'+key+'.yade.gz')
#print "###    Compacted state saved      ###"

################################################################################
# 2nd phase: Confining the specimen to desired value of stress of confinement,
# reduce the cinetic energy
# will implement in this state the code for generating the granulometric curve
################################################################################

while 1:
  triax.internalCompaction=False
  triax.stressMask=7
  setContactFriction=radians(compFricDegree)
  triax.goal1=table.conStress
  triax.goal2=table.conStress
  triax.goal3=table.conStress
  O.run(1000, True)
  #the global unbalanced force on dynamic bodies, thus excluding boundaries,
  #which are not at equilibrium
  unb=unbalancedForce()
  #average stress
  #note: triax.stress(k) returns a stress vector, so we need to keep only the normal component
  meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
  print 'unbalanced force:',unb,' mean stress: ',meanS
  print 'void ratio=',triax.porosity/(1-triax.porosity), 'porosity=', triax.porosity
  print 'mean stress engine', triax.meanStress, 'kineticE=', utils.kineticEnergy()
  print '-----Current State_02: Compress 100kPa--v(-___-)v--'
  if unb<stabilityThreshold and abs(meanS-table.conStress)/table.conStress<0.001:
    break

O.save('compressState'+key+'.xml') # important, put to meromy this file to start the deviatoric loading later

print 'step that starts the deviatoric loading ', O.iter
e22Check=triax.strain[1]
# We will save a check point, to know that from this point the deviatoric loading process begins, so we know where to start the ploting shit
from yade import plot
from pprint import pprint
plot.reset()
plot.addData(CheckpointStep=O.iter,Porosity=triax.porosity,e11=triax.strain[0],e22=triax.strain[1],e33=triax.strain[2],elasticE=law.elasticEnergy(),kineticE=utils.kineticEnergy())
plot.saveDataTxt('checkpoint.txt',vars=('CheckpointStep','Porosity','e11','e22','e33','elasticE','kineticE'))

###############################################
# 3rd phase: DEVIATORIC LOADING (BIG SHOW!!!) #
###############################################

#We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False

# Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction=radians(35)
#TriaxialStateRecorder(iterPeriod=20,file='WallStresses_load_deviatoric'+key,addIterNum=True,initRun=True)
while 1:
    triax.stressMask=5
    triax.goal2=-rate
    triax.goal1=triax.goal3=table.conStress
    O.run(100, True)
    #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
    unb=unbalancedForce()
    #note: triax.stress(k) returns a stress vector, so we need to keep only the normal component
    axialS=triax.stress(triax.wall_top_id)[1]
    print 'step=', O.iter, 'unbalanced force:',unb,' sigma2: ',axialS, 'q=', axialS-table.conStress
    print 'axial deformation (%)', (triax.strain[1]-e22Check)*100
    if abs((triax.strain[1]-e22Check)-0.01)<=0.001:
        O.save('firstpoint.xml')
    if abs((triax.strain[1]-e22Check)-0.02)<=0.001:
        O.save('secondpoint.xml')
    if abs((triax.strain[1]-e22Check)-0.03)<=0.001:
        O.save('thirdpoint.xml')
    if abs((triax.strain[1]-e22Check)-0.05)<=0.001:
        O.save('fourthpoint.xml')
    if triax.strain[1]-e22Check>=0.3 : #triax.sigma2
        break

O.save('final.xml')