#############################	Vertical compression	#############################
#Made by GONCALVES DE OLIVEIRA CHUEIRE Joao Augusto as part of PhD at INSA de Lyon

#
##
### Description : insert a set of particles into a rectangular box allowing them 
### to settle under gravity. After, add controlled movement to 5 pistons, one in z-,
### two in x and two in y. A constant stress will be applied untill equilibrium. 
### Then the stress on the z piston will be increased wa certain vertical strain value. 
### We will then observe the way the grains organize themselfs
##
#

#First part	=> Grains Insertion, save grains values & restart simulation
#Second part	=> Load grains, restart simulation properties & pistZ calibration
#Third part	=> Consolidation
#Fourth part	=> Compression

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#############################		Variables		#############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

########## Exterior parameterss
variable	ExpWidth 	equal 0.1	#experiment widith, x in meters
variable	ExpLength 	equal 0.35	#experiment length, y in meters
variable	ExpHeight	equal 0.45	#experiment height, z in meters
variable	MaxRadius	equal 0.008	#maximal radius of the simulation
variable	MaxE		equal 0.3	#strain target to stop simulation 0.2=20%

########## Material parameters
variable	Eg 		equal 1e8	#Young modulus grains 
variable	Ew 		equal 1e8	#Young modulus walls
variable	V 		equal 0.25	#Possoin coef
variable	u 		equal 0.57	#friction coefficient
variable	e 		equal 0.5	#restitution coefficient
variable	rf 		equal 0.3	#rolling friction coefficient

########## Simulation parameters
	#global values
variable	NeighSize	equal 0.01	#size of neighbor bins
variable	dt 	  	equal 2e-5	#timestep value in seconds
	#insertion values
variable	Ngrains   	equal 16e3	#number of grains
variable	InsrNb		equal 10e3 	#Number of particles inserted per second
variable	InsrInt		equal 2.5e3	#Interval between insertions
	#saving intervals
variable	Intv 	  	equal 1250	#Interval between savefiles
	#running times
variable	run1 	  	equal 100e3	#nb of steps insertion (first part)
variable	run1b	  	equal 25e3	#nb of steps calibration (second part)
variable	run2 	  	equal 100e3	#nb of steps consolidation (third part)
variable	run3 	  	equal 700e3	#nb of steps compression (fourth part)
variable	run3b 	  	equal 25e3	#interval between pistonZ remakes (also fourth part)
	#piston parameters
variable	CpStress 	equal 1.5e6	#maximal vertical stress, plan x.y in Pa
variable	CsStress 	equal 1e5	#maximal horizontal stress, plan y.z in Pa
variable	CsSpeed  	equal 0.02	#speed of pistons in Consolidation part
variable	CpSpeedH 	equal 0.5	#speed of horizontal pistons in Compression part
variable	CpSpeedV 	equal 0.02	#speed of vertical pistons in Compression part
variable	KP 	  	equal 1		#value to control piston advancing speed on manual mode	
variable	KD		equal 1e-6

########## File details
variable	part 		string post1
variable	part2 		string post2
variable	partg 		string generalvariables

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#############################		Prepare Simulations		#############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

########## Small Calculations
	#geometric calculations
variable	ExpWidth2 	equal $(v_ExpWidth/2) 			#half box witdh, value used on the grain insertion
	#vertical force
variable 	Fz1		equal $(v_CsStress*v_ExpLength*v_ExpWidth) #used for consolidation
variable 	Fz2		equal $(v_CpStress*v_ExpLength*v_ExpWidth) #used for compression
variable	dFz		equal $((v_Fz2-v_Fz1)/5)
	#horizontal force
variable 	Fy		equal $(v_CsStress*v_ExpWidth*v_ExpHeight)
variable 	Fx		equal $(v_CsStress*v_ExpLength*v_ExpHeight)
	#horizontal consolidation speed - take into account the diferent surfaces
variable 	CsSpeedY	equal $(v_CsSpeed*v_ExpLength/v_ExpHeight)
variable 	CsSpeedX	equal $(v_CsSpeed*v_ExpWidth/v_ExpHeight)
	#step calculations
variable	Intvmesh	equal $(20*v_Intv)		#value between each savefiles for non-moving mesh files
variable	consoStep 	equal $(v_run1b)
variable	compStep 	equal $(v_consoStep+v_run2)
variable	endStep 	equal $(v_compStep+v_run3)

########## Loop variables
variable 	cntr 		equal 0 #How many times run3b was executed
variable	eps		equal 0	#Calculate lagrangian deformation
variable	dh		equal 0	#Height storing for deformation calculation

########## Print For Matlab
#The following lines will be printed in the top of the 'log.liggghts' file and will be read with matlab
	#geometry
print		""
print		"ExpWidth=${ExpWidth}"
print		"ExpLength=${ExpLength}"
print		"ExpHeight=${ExpHeight}"
	#step
print		"Timestep=${dt}"
print		"StartCons=${consoStep}"
print		"EndCons=${compStep}"
print		"EndComp=${endStep}"
print		"Intv=${Intv}"
print		"PistZchInt=${run3b}"
	#folder
print		"vtkFolder=${part2}"
print		"genFolder=${partg}"
print		""

########## Execute STL files
variable	stlModes index  1 21 22 31 32 4 5 -1	#radius for stress/strain calculation
variable	stlFname index  Boxbase PistonY1 PistonY2 PistonX1 PistonX2 PistonZ InsHelp	
#command to create a new stl using openscad file contained in folder meshes
label doStl
shell cd meshes
shell openscad -Dmode=${stlModes} -Dw=$(v_ExpWidth) -Dl=$(v_ExpLength) -Dh=$(v_ExpHeight) &
	 -o ${stlFname}.stl BOX.scad
next stlModes
next stlFname
shell cd ..
jump Box.liggghts doStl

### LOAD DATA ###
#if file containing grains data was found, the insertion of grains will not be done
if '$(v_tfI==1)' then 'jump Box.liggghts loadFile'
if '$(v_tfI>1)' then 'read_restart comp*.resData' &
	'jump Box.liggghts loadFile2'  

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#############################		First Part -Insertion		#############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## First part of the simulation. Start the simulation parameters, then create a volume for insertion 
## and start inserting the grains. A great amount of grains will be inserted, excess grains will then 
## be removed and the remaining will be saved in a 'part1Data.txt' file.

##DONE BY OTHER FILE

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#############################		Second Part - Calibration		#############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Grains file will be Loaded. As the simulation was cleared, all the initialization parameters must be remade.
## Then the pisntonZ will be slowly advanced to touch the grains. The max stress it can aplly is really small
## so it should not impact the simulation.

label loadFile 

##############Prepare simulation
#### Initialization ####

# Preliminaries
units		si
atom_style	sphere
atom_modify	map array
boundary	f f f	# describes boundaries of the domain as fixed(f) or periodic (p)
label 		loadFile2 	#load from restart file
newton		off      #for when shortrange operations
communicate	single vel yes
if '$(v_Eg>1e9)' then 'hard_particles yes'

if '$(v_tfI>1)'	then 'jump Box.liggghts loadFile2a' 
read_data	part1Data.txt #Read data from the insertion file
label		loadFile2a

#Neighbor Listing  
neighbor	${NeighSize} bin
neigh_modify	delay 0 

#Material and interaction properties
fix		m1 all property/global youngsModulus peratomtype ${Eg} ${Ew}
fix		m2 all property/global poissonsRatio peratomtype ${V} ${V}
fix		m3 all property/global coefficientRestitution peratomtypepair 2 ${e} 0.1 0.1 0.1
fix		m4 all property/global coefficientFriction peratomtypepair 2 ${u} 0. 0. 0. 
fix		m5 all property/global coefficientRollingFriction peratomtypepair 2 ${rf} 0. 0. 0.

#### Detailed Settings ####

#Define the particle physics
pair_style	gran model hertz tangential history rolling_friction epsd2
pair_coeff	* *

#Integrator
fix		integrate all nve/sphere 

#Gravity
fix		grav all gravity 9.81 vector 0. 0. -1.0

#Time step
timestep	${dt}
    #Typically try to keep <= 20% of  the Rayleigh time

#Thermodynamic output settings
thermo_style	custom step atoms 	#quantity to be printed in the logfile
thermo		$(2*v_Intv)		#how often to write the thermo quantities
thermo_modify	norm no lost ignore 	#Set options. 'norm' normalizes output on a per

##############Prepare for calibration
#Import Mesh from Cad
	#all cad files were drawn in cm, so they have to be scaled by a 0.01 factor
fix		mBse all mesh/surface file meshes/Boxbase.stl 	type 2 scale 1. #box base

	#pistons as servo
fix 		sPz  all mesh/surface/stress/servo file meshes/PistonZ.stl type 2 scale 1. &
	com 0. 0. 0. ctrlPV force axis 0. 0. -1. target_val $(v_Fz1*0.01) vel_max ${CsSpeed} kp ${KP} kd ${KD}	
fix 		sPy1 all mesh/surface/stress/servo file meshes/PistonY1.stl type 2 scale 1. &
	com 0. 0. 0. ctrlPV force axis 0. 1. 0. target_val $(v_Fy*0.01) vel_max ${CsSpeedY} kp ${KP} kd ${KD}
fix 		sPy2 all mesh/surface/stress/servo file meshes/PistonY2.stl type 2 scale 1. &
	com 0. 0. 0. ctrlPV force axis 0. -1. 0. target_val $(v_Fy*0.01) vel_max ${CsSpeedY} kp ${KP} kd ${KD}
fix		sPx1 all mesh/surface/stress/servo file meshes/PistonX1.stl type 2 scale 1. &
	com 0. 0. 0. ctrlPV force axis 1. 0. 0. target_val $(v_Fx*0.01) vel_max ${CsSpeedX} kp ${KP} kd ${KD}
fix		sPx2 all mesh/surface/stress/servo file meshes/PistonX2.stl type 2 scale 1. & 
	com 0. 0. 0. ctrlPV force axis -1. 0. 0. target_val $(v_Fx*0.01) vel_max ${CsSpeedX} kp ${KP} kd ${KD} 

#block the movement of all pistons but the pistonZ, since it is the only one that does not touch the grains
fix_modify	sPy1 		servo/integrate stop 
fix_modify	sPy2 		servo/integrate stop 
fix_modify	sPx1 		servo/integrate stop 
fix_modify	sPx2 		servo/integrate stop

#Change thermo so we cna see what is happening in the terminal
thermo_style    custom step atoms f_sPz[9] f_sPy1[8] f_sPy2[8] f_sPx1[7] f_sPx2[7] 

#Remaking the walls
fix		walls1 all wall/gran model hertz tangential history  rolling_friction epsd2 &
	mesh n_meshes 6 meshes mBse sPy1 sPy2 sPx1 sPx2 sPz 

#Adding new computing elements, to calculate important information
compute 	ctc all 	contact/atom	#contacts nb 
compute 	dpl all 	displace/atom	#each particle displacement
compute		pgl all 	pair/gran/local pos id force contactPoint delta #grain-grain contact
compute		sts all 	stress/atom 	#particle stress through Love-Webber
#compute		ken all ke/atom		#atom knectic energy
#compute		pen all pe/atom		#atom potential energy
	
#adding new dump files for drawing data -- paraview (.vtk)
dump		dmpSph2 all custom/vtk ${Intv} ${part2}/partcles_*.vtk x y z radius #c_ken
dump		dmpSph3 all custom/vtk ${Intv} ${part2}/displacement_*.vtk c_dpl[1] &
	c_dpl[2] c_dpl[3] c_dpl[4]
dump		frchain all local/gran/vtk ${Intv} ${part2}/forcechain*.vtk pgl

#adding new dump files for the meshes -- paraview (.vtk)
dump		dmpBas all mesh/gran/VTK ${Intvmesh} ${part2}/mesh_base_*.vtk stress mBse
dump		dmpsPz all mesh/gran/VTK ${Intv} ${part2}/mesh_sPz_*.vtk stress sPz
dump		dmpPy1 all mesh/gran/VTK ${Intv} ${part2}/mesh_sPy1_*.vtk stress sPy1
dump		dmpPy2 all mesh/gran/VTK ${Intv} ${part2}/mesh_sPy2_*.vtk stress sPy2
dump		dmpPx1 all mesh/gran/VTK ${Intv} ${part2}/mesh_sPx1_*.vtk stress sPx1
dump		dmpPx2 all mesh/gran/VTK ${Intv} ${part2}/mesh_sPx2_*.vtk stress sPx2

if '$(v_tfI>1)' then 'jump Box.liggghts loadFile2b' 
run ${run1b}
label loadFile2b 

variable	VertDpl1 equal $(f_sPz[9])
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#############################		Third Part - Consolidation		#############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## In this part all the pistons will apply the same stress to the grains. The function 'servo' is actually con-
## trolled by a Force, this value is calculated at the start of the file.

#restart pistons movement
fix_modify	sPy1 servo/integrate start 
fix_modify	sPy2 servo/integrate start 
fix_modify	sPx1 servo/integrate start 
fix_modify	sPx2 servo/integrate start

#modify the target force of the all pistons
fix_modify	sPz  servo/target_val ${Fz1} 
fix_modify	sPy1 servo/target_val ${Fy} 
fix_modify	sPy2 servo/target_val ${Fy} 
fix_modify	sPx1 servo/target_val ${Fx} 
fix_modify	sPx2 servo/target_val ${Fx} 

#adding new dump files for general data -- matlab (.txt)
fix		servostress all ave/time 1 1 ${Intv} f_sPz[3] f_sPz[9] f_sPy1[2] & 
	f_sPy1[8] f_sPy2[2] f_sPy2[8] f_sPx1[1] f_sPx1[7] f_sPx2[1] f_sPx2[7] &
	file ${partg}/servoForce.txt
dump            grains all custom ${Intv} ${partg}/grains*.txt id x y z radius c_ctc &
	c_dpl[1] c_dpl[2] c_dpl[3] c_sts[1] c_sts[2] c_sts[3] c_sts[4] c_sts[5] c_sts[6] &
	omegax omegay omegaz # c_ken c_pen
dump            contactForces all local ${Intv} ${partg}/contForce*.txt c_pgl[7] &
	c_pgl[8] c_pgl[10] c_pgl[11] c_pgl[12] c_pgl[13] c_pgl[14] c_pgl[15] c_pgl[16]
	#cont forces ID1 ID2 Fx Fy Fz Cpx Cpy Cpz Delta

if '$(v_tfI>1)' then 'jump Box.liggghts loadFile2c' 
run ${run2}
label loadFile2c 

variable	VertDpl2 equal $(f_sPz[9])

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#############################		Fourth Part - Compression		#############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## All pistons will mantain the same stress, while the pistonZ will have it's value augmented with the objective
## of causing the failure of the created granular material.

#Change the piston speeds to the right value
fix_modify	sPz  servo/vel_max ${CpSpeedV}
fix_modify	sPy1 servo/vel_max ${CpSpeedH}
fix_modify	sPy2 servo/vel_max ${CpSpeedH}
fix_modify	sPx1 servo/vel_max ${CpSpeedH}
fix_modify	sPx2 servo/vel_max ${CpSpeedH}

variable 	dh equal $(f_sPz[9])

#Save restart file
restart ${run3b} reFolder/comp*.resData

### Compression Loop
label doComp 

#the target value of the Piston Z will be increased gradually in the start so it does not pack a punch right into the beggining
if "$(v_cntr+1) <= 5" then "fix_modify sPz servo/target_val $(v_Fz1 +(v_cntr+1)*v_dFz )"

#Check the deformation
variable 	eps 	equal $(v_eps-(f_sPz[9]-v_dh)/(v_ExpHeight+v_dh))
if "$(v_eps) >= $(v_MaxE)" then "jump  Box.liggghts goEnd"
variable 	dh 	equal $(f_sPz[9])	#update dh for next loop

#Recalculate piston Y and X target value to adapt to changing Pz
variable 	nFy	equal $(v_CsStress*(v_ExpWidth-f_sPx1[7]+f_sPx2[7])*(v_ExpHeight+f_sPz[9]))
variable 	nFx	equal $(v_CsStress*(v_ExpLength-f_sPy1[8]+f_sPy2[8])*(v_ExpHeight+f_sPz[9]))

fix_modify	sPy1 	servo/target_val ${nFy} 
fix_modify	sPy2 	servo/target_val ${nFy} 
fix_modify	sPx1 	servo/target_val ${nFx} 
fix_modify	sPx2 	servo/target_val ${nFx} 

variable 	cntr 	equal $(v_cntr+1) #count another run3b
run ${run3b} 

#if we run enough run3b to equal or surpass run3 then we can end the simulation, else go back to newPiston
if " $(v_run3b * v_cntr) < ${run3}" then "jump  Box.liggghts doComp"

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#############################		End Print Var		#############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
label goEnd

variable	endStep equal $(v_compStep+v_run3b*v_cntr)
variable	consDpl equal $(v_VertDpl2-v_VertDpl1) #consolidation vertical displacement
variable	compDpl equal $(f_sPz[9]-v_VertDpl2)#compression vertical displacement

#print values
print		"     "
print		"Vertical consolidation force : ${Fz1} N"
print		"Vertical compression force : ${Fz2} N"
print		"Horizontal starting force Y: ${Fy} N"
print		"Horizontal ending force Y: ${nFy} N"
print		"Horizontal starting force X: ${Fx} N"
print		"Horizontal ending force X: ${nFx} N"
print		"     "
print		"Consolidation Displacement :${consDpl} m"
print		"Mean Consolidation Speed :$(v_consDpl/((v_compStep-v_consoStep)*v_dt)) m/s"
print		"Consolidation Displacement :${compDpl} m"
print		"Mean Consolidation Speed :$(v_compDpl/((v_endStep-v_compStep)*v_dt)) m/s"


########## Print For Matlab
#The following lines will be printed in the top of the 'LiggghtstoMatlab.txt' file and will be read with matlab
log LiggghtstoMatlab.txt
	#geometry
print		"The following values give the Matlab app information about the simulaiton"
print		""
print		"ExpWidth=${ExpWidth}"
print		"ExpLength=${ExpLength}"
print		"ExpHeight=${ExpHeight}"
	#step
print		"Timestep=${dt}"
print		"StartCons=${consoStep}"
print		"EndCons=${compStep}"
print		"EndComp=${endStep}"
print		"Intv=${Intv}"
print		"PistZchInt=${run3b}"
	#folder
print		"vtkFolder=${part2}"
print		"genFolder=${partg}"
print		""