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

region reg block 0.0 0.1 0 0.1 -0.0005 0.0005 units box
variable confinement_pressure equal 80000

# 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

region region_gouge_1 block 0.00 0.1 0.0000 0.0025 -0.0005 0.0005 units box
region region_gouge_2 block 0.00 0.1 0.0025 0.0050 -0.0005 0.0005 units box
region region_gouge_3 block 0.00 0.1 0.0050 0.0075 -0.0005 0.0005 units box
region region_gouge_4 block 0.00 0.1 0.0075 0.0100 -0.0005 0.0005 units box
region region_gouge_5 block 0.00 0.1 0.0100 0.0125 -0.0005 0.0005 units box
region region_gouge_6 block 0.00 0.1 0.0125 0.0150 -0.0005 0.0005 units box
region region_gouge_7 block 0.00 0.1 0.0150 0.0175 -0.0005 0.0005 units box
region region_gouge_8 block 0.00 0.1 0.0175 0.0200 -0.0005 0.0005 units box
region region_gouge_9 block 0.00 0.1 0.0200 0.0225 -0.0005 0.0005 units box
region region_gouge_10 block 0.00 0.1 0.0225 0.0250 -0.0005 0.0005 units box
region region_gouge_11 block 0.00 0.1 0.0250 0.0275 -0.0005 0.0005 units box
region region_gouge_12 block 0.00 0.1 0.0275 0.0300 -0.0005 0.0005 units box
region region_gouge_13 block 0.00 0.1 0.0300 0.0325 -0.0005 0.0005 units box
region region_gouge_14 block 0.00 0.1 0.0325 0.0350 -0.0005 0.0005 units box
region region_gouge_15 block 0.00 0.1 0.0350 0.0375 -0.0005 0.0005 units box
region region_gouge_16 block 0.00 0.1 0.0375 0.0400 -0.0005 0.0005 units box
region region_gouge_17 block 0.00 0.1 0.0400 0.0425 -0.0005 0.0005 units box
region region_gouge_18 block 0.00 0.1 0.0425 0.0450 -0.0005 0.0005 units box
region region_gouge_19 block 0.00 0.1 0.0450 0.0475 -0.0005 0.0005 units box
region region_gouge_20 block 0.00 0.1 0.0475 0.0500 -0.0005 0.0005 units box
region region_gouge_21 block 0.00 0.1 0.0500 0.0525 -0.0005 0.0005 units box
region region_gouge_22 block 0.00 0.1 0.0525 0.0550 -0.0005 0.0005 units box
region region_gouge_23 block 0.00 0.1 0.0550 0.0575 -0.0005 0.0005 units box
region region_gouge_24 block 0.00 0.1 0.0575 0.0600 -0.0005 0.0005 units box
region region_gouge_25 block 0.00 0.1 0.0600 0.0625 -0.0005 0.0005 units box
region region_gouge_26 block 0.00 0.1 0.0625 0.0650 -0.0005 0.0005 units box
region region_gouge_27 block 0.00 0.1 0.0650 0.0675 -0.0005 0.0005 units box
region region_gouge_28 block 0.00 0.1 0.0675 0.0700 -0.0005 0.0005 units box
region region_gouge_29 block 0.00 0.1 0.0700 0.0725 -0.0005 0.0005 units box
region region_gouge_30 block 0.00 0.1 0.0725 0.0750 -0.0005 0.0005 units box
region region_gouge_31 block 0.00 0.1 0.0750 0.0775 -0.0005 0.0005 units box
region region_gouge_32 block 0.00 0.1 0.0775 0.0800 -0.0005 0.0005 units box
region region_gouge_33 block 0.00 0.1 0.0800 0.0825 -0.0005 0.0005 units box
region region_gouge_34 block 0.00 0.1 0.0825 0.0850 -0.0005 0.0005 units box
region region_gouge_35 block 0.00 0.1 0.0850 0.0875 -0.0005 0.0005 units box
region region_gouge_36 block 0.00 0.1 0.0875 0.0900 -0.0005 0.0005 units box
region region_gouge_37 block 0.00 0.1 0.0900 0.0925 -0.0005 0.0005 units box
region region_gouge_38 block 0.00 0.1 0.0925 0.0950 -0.0005 0.0005 units box
region region_gouge_39 block 0.00 0.1 0.0950 0.0975 -0.0005 0.0005 units box
region region_gouge_40 block 0.00 0.1 0.0975 0.1000 -0.0005 0.0005 units box

group nve_group region region_gouge_1
group nve_group region region_gouge_2
group nve_group region region_gouge_3
group nve_group region region_gouge_4
group nve_group region region_gouge_5
group nve_group region region_gouge_6
group nve_group region region_gouge_7
group nve_group region region_gouge_8
group nve_group region region_gouge_9
group nve_group region region_gouge_10
group nve_group region region_gouge_11
group nve_group region region_gouge_12
group nve_group region region_gouge_13
group nve_group region region_gouge_14
group nve_group region region_gouge_15
group nve_group region region_gouge_16
group nve_group region region_gouge_17
group nve_group region region_gouge_18
group nve_group region region_gouge_19
group nve_group region region_gouge_20
group nve_group region region_gouge_21
group nve_group region region_gouge_22
group nve_group region region_gouge_23
group nve_group region region_gouge_24
group nve_group region region_gouge_25
group nve_group region region_gouge_26
group nve_group region region_gouge_27
group nve_group region region_gouge_28
group nve_group region region_gouge_29
group nve_group region region_gouge_30
group nve_group region region_gouge_31
group nve_group region region_gouge_32
group nve_group region region_gouge_33
group nve_group region region_gouge_34
group nve_group region region_gouge_35
group nve_group region region_gouge_36
group nve_group region region_gouge_37
group nve_group region region_gouge_38
group nve_group region region_gouge_39
group nve_group region region_gouge_40

fix ins_1 nve_group pour 50000 1 952 region region_gouge_1 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_2 nve_group pour 50000 1 81 region region_gouge_2 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_3 nve_group pour 50000 1 699 region region_gouge_3 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_4 nve_group pour 50000 1 113 region region_gouge_4 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_5 nve_group pour 50000 1 2 region region_gouge_5 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_6 nve_group pour 50000 1 642 region region_gouge_6 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_7 nve_group pour 50000 1 220 region region_gouge_7 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_8 nve_group pour 50000 1 566 region region_gouge_8 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_9 nve_group pour 50000 1 855 region region_gouge_9 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_10 nve_group pour 50000 1 997 region region_gouge_10 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_11 nve_group pour 50000 1 736 region region_gouge_11 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_12 nve_group pour 50000 1 225 region region_gouge_12 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_13 nve_group pour 50000 1 385 region region_gouge_13 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_14 nve_group pour 50000 1 403 region region_gouge_14 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_15 nve_group pour 50000 1 638 region region_gouge_15 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_16 nve_group pour 50000 1 130 region region_gouge_16 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_17 nve_group pour 50000 1 53 region region_gouge_17 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_18 nve_group pour 50000 1 684 region region_gouge_18 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_19 nve_group pour 50000 1 730 region region_gouge_19 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_20 nve_group pour 50000 1 672 region region_gouge_20 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_21 nve_group pour 50000 1 710 region region_gouge_21 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_22 nve_group pour 50000 1 416 region region_gouge_22 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_23 nve_group pour 50000 1 247 region region_gouge_23 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_24 nve_group pour 50000 1 836 region region_gouge_24 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_25 nve_group pour 50000 1 439 region region_gouge_25 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_26 nve_group pour 50000 1 203 region region_gouge_26 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_27 nve_group pour 50000 1 184 region region_gouge_27 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_28 nve_group pour 50000 1 123 region region_gouge_28 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_29 nve_group pour 50000 1 401 region region_gouge_29 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_30 nve_group pour 50000 1 767 region region_gouge_30 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_31 nve_group pour 50000 1 294 region region_gouge_31 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_32 nve_group pour 50000 1 280 region region_gouge_32 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_33 nve_group pour 50000 1 837 region region_gouge_33 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_34 nve_group pour 50000 1 884 region region_gouge_34 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_35 nve_group pour 50000 1 610 region region_gouge_35 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_36 nve_group pour 50000 1 198 region region_gouge_36 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_37 nve_group pour 50000 1 982 region region_gouge_37 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_38 nve_group pour 50000 1 907 region region_gouge_38 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_39 nve_group pour 50000 1 511 region region_gouge_39 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5
fix ins_40 nve_group pour 50000 1 752 region region_gouge_40 vol 0.2 200 diam poly 2 0.0025 0.5 0.0035 0.5

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

unfix ins_1
unfix ins_2
unfix ins_3
unfix ins_4
unfix ins_5
unfix ins_6
unfix ins_7
unfix ins_8
unfix ins_9
unfix ins_10
unfix ins_11
unfix ins_12
unfix ins_13
unfix ins_14
unfix ins_15
unfix ins_16
unfix ins_17
unfix ins_18
unfix ins_19
unfix ins_20
unfix ins_21
unfix ins_22
unfix ins_23
unfix ins_24
unfix ins_25
unfix ins_26
unfix ins_27
unfix ins_28
unfix ins_29
unfix ins_30
unfix ins_31
unfix ins_32
unfix ins_33
unfix ins_34
unfix ins_35
unfix ins_36
unfix ins_37
unfix ins_38
unfix ins_39
unfix ins_40

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