atom_style  sphere	
atom_modify	map array
dimension	2
boundary	p p p
newton		off
comm_modify	vel yes
units		si


# How many grain/atom types
create_box	6 reg

# Size of bins for more efficiently searching for grain contacts.
neighbor	0.001 bin

# Reconstruct the neighbor list without any delay, every time-step
neigh_modify	delay 0

# Glass marbles, tangential force
pair_style gran/hertz/history 36630036630.0 0.0 0.2 NULL 0.0 0

pair_coeff * *

timestep        1e-8

fix             gravi all gravity 0.0 vector 0.0 -1.0 0.0

# Particle insertion in regions

run 1

set     group nve_group density/disc 2.5 
fix		integr nve_group nve/sphere disc
fix		makeit2d all enforce2d

thermo		1
thermo_modify	lost ignore norm no

# Unfix particle insertion

run		100 upto

# Set gouge layer grain density
set             group all density/disc 2.5
   
# Apply NVE integration to all particles
fix		integr all nve/sphere disc

# Output settings
compute		1 all erotate/sphere
compute		2 all contact/atom
compute         3 all ke

variable        Sxx equal pxx
variable        Syy equal pyy
variable        TotalPressure equal (v_Sxx+v_Syy)/2.0

thermo_style    custom step atoms ke pxx pyy v_TotalPressure ly
thermo		50000
thermo_modify	lost ignore norm no

set             group all density/disc 2.5

# Stop the confining pressure once the pressure is superior to the desired pressure
fix             condition all halt 1 v_Syy > ${confinement_pressure} error continue
fix             def all deform 1 x erate -50 y erate -50 z erate 0

run             1500000
unfix           def
unfix           condition

run 10000

# Check if the pressure does not decrease
label       loopa
variable    a loop 1000
  label     loopb
  variable P equal pyy
  variable  b loop 1000
    if "${P} > ${confinement_pressure}" then "jump SELF break"
    fix def all deform 1 x erate -10 y erate -10 z erate 0
    run 1000
    unfix def
    run 10000
  next      b
  jump      SELF loopb
  label       break
  variable    b delete
if "${P} > ${confinement_pressure}" then "jump SELF break2"
next        a
jump        SELF loopa
label break2
print "Over" 

restart 1000000 confined.restart confined2.restart

run 1000000

write_data box_confined.data