#-----------------------------------------------------------------------------------
#Simple shearbox test
#like J. K. Morgan 2004 but in 3D 
#
#Compute on own LIGGGHTS
#
#Alexandre Sac--Morane
#-----------------------------------------------------------------------------------

#-----------------------------------------------------------------------------------
#Initialization
#-----------------------------------------------------------------------------------

units         si
atom_style    granular
atom_modify   map array
boundary      p f p
newton        off

communicate   single vel yes

#region of the study domain
region        domain block 0.0 0.02 0.0 0.02 0.0 0.001 units box
create_box    1 domain

#regions of insertion
region       ri block 0.0 0.02 0.0 0.02 0.0 0.001 units box

neighbor      0.0008 bin
neigh_modify  delay 0

timestep    0.000001

fix gravi all gravity 9.81 vector 0.0 -1.0 0.0

#-----------------------------------------------------------------------------------
#Material properties required for new pair styles
#-----------------------------------------------------------------------------------

fix  m1 all property/global youngsModulus peratomtype 340.e6
fix  m2 all property/global poissonsRatio peratomtype 0.2
fix  m3 all property/global coefficientRestitution peratomtypepair 1 1
fix  m4 all property/global coefficientFriction peratomtypepair 1 0.3

#New pair style
pair_style  gran model hertz/own tangential history
pair_coeff  * *

#-----------------------------------------------------------------------------------
#Definition of the bottom plate
#-----------------------------------------------------------------------------------

fix cad1 all mesh/surface/stress file meshes/Face.stl type 1 &
	scale 1.0 move 0. 0. 0. 
fix yWallsLow all wall/gran model hertz tangential history &
	 mesh n_meshes 1 meshes cad1

#-----------------------------------------------------------------------------------
#Definition of the particles insertion
#-----------------------------------------------------------------------------------

#Particles
fix pts1 all particletemplate/sphere 32452843 atom_type 1 &
	density constant 1000 radius constant 0.0002	
fix pts2 all particletemplate/sphere 49979687 atom_type 1 &
	density constant 1000 radius constant 0.00015
fix pts3 all particletemplate/sphere 67867979 atom_type 1 &
	density constant 1000 radius constant 0.0001
	
fix pdd all particledistribution/discrete/massbased 15485863 3 pts1 0.42 pts2 0.36 pts3 0.22

fix ins all insert/pack seed 86028121 distributiontemplate pdd &
	vel constant 0. 0. 0.0 insert_every once overlapcheck yes all_in yes&
	particles_in_region 8750 region ri
	
# Integrator
fix integrate all nve/sphere

#-----------------------------------------------------------------------------------
#Definition of the ouput
#-----------------------------------------------------------------------------------

thermo_style  custom step ke
thermo        1000
thermo_modify lost ignore norm no

#-----------------------------------------------------------------------------------
#Particles insertion
#-----------------------------------------------------------------------------------

#insert the first particles so that dump is not empty

run    1 upto
unfix  ins

#-----------------------------------------------------------------------------------
#Particles settling
#-----------------------------------------------------------------------------------

run    15000 upto

#-----------------------------------------------------------------------------------
#Definition of the top plate
#-----------------------------------------------------------------------------------

fix    cad2 all mesh/surface/stress/servo file meshes/Face.stl type 1&
	scale 1.0 move 0. 0.009 0. com 0.01 0. 0.0005 ctrlPV force &
	axis 0. 1. 0. target_val -100.0 vel_max 60

#Updating the wall/gran contact
unfix  yWallsLow
fix    yWalls all wall/gran model hertz/own tangential history &
	 mesh n_meshes 2 meshes cad1 cad2 

#-----------------------------------------------------------------------------------
#Definition of the dump
#-----------------------------------------------------------------------------------

dump   dmp all custom/vtk 200 post/jkmorgan_*.vtk id &
	x y z vx vy vz omegax omegay omegaz radius  
dump   dumpstl all mesh/stl 200 post/dump*.stl

#-----------------------------------------------------------------------------------
#Update of the ouput
#-----------------------------------------------------------------------------------

thermo_style  custom time ke f_cad1[2] f_cad2[2] f_cad2[8]

#-----------------------------------------------------------------------------------
#Applying normal force by moving down the top plate 
#-----------------------------------------------------------------------------------

run 4000

#-----------------------------------------------------------------------------------
#Applying shear force by moving the bottom plate 
#-----------------------------------------------------------------------------------

fix movecad all move/mesh mesh cad1 linear 1.0e-7 0.0 0.0

#The plate shearing
run 10000