;version: PFC2D 3.10-234
;
;
;E-mail: sheng0619@163.com
;
;2017/06/08 LI ChangSheng @ NanJing University 
;           add viscosity
;           see more 
;           [1]LI ChangSheng, YIN HongWei, JIA Dong. Validation tests for discrete element codes using single-contact systems
;           T1 viscosity=0.0,   T2 viscosity=3e8
;2015/09/20 LI ChangSheng @ NanJing University 
;2014/11/07 LI ChangSheng @ NanJing University 
;----------------------------------------------
new
call fishcall.FIS ; load in the FISHCALL macros
set log on
set disk on
def ini_set
	;change den to the model in our manuscript, 
	;so that we used the same mass in the simulation
	;PFC2D     , mass=den*pi*rad^2
	;ref.
	;[1]Itasca (2004). "PFC2D Users�� Manual (version3.10)." Minnesota: Itasca Consulting Group Inc.
	;our model , mass=den*2*sqrt(3)*rad^2    
	;ref.
	;[1]Place, D., Lombard, F., Mora, P., and Abe, S. (2002). "Simulation of the micro-physics of rocks using LSMearth." Earthquake Processes: Physical Modelling, Numerical Simulation and Data Analysis Part I, Springer, 1911-1932.
	;[2]Liu, C., Pollard, D. D., and Shi, B. (2013). "Analytical solutions and ;numerical tests of elastic and failure behaviors of close-packed lattice for brittle rocks and crystals." Journal of Geophysical Research: Solid Earth, ;118(1), 71-82.
	den=2500
	realden=den*(2*sqrt(3))/pi
	viscosity=0e8 ; T1 0.0
	;viscosity=3e8 ; T2 3e8
	mystep=0;
end
ini_set

def apply_vis
	;
	; ----- If initial ball (bp1) falls below the 5-m level, then
	; split this ball into two balls, each of half the original radius.
	;
	bp = ball_head
	loop while bp # null
		bid = b_id(bp)
		b_xfvis  = - viscosity * b_xvel(bp)
		b_yfvis  = - viscosity * b_yvel(bp) 
		command
			prop xforce b_xfvis yforce b_yfvis range id bid
		end_command
		bp = b_next(bp)
	end_loop
end

;----------------------------------------------------------------------------
ball rad 50 id 1 x 50.0 y -100.0
ball rad 50 id 2 x 50.0 y    0.0

prop density realden kn 11e9 ks 11e9 fric 0.0 damping 0.0 ; contact stiffness = 11e9 * 11e9 / (11e9 + 11e9) = 5.5e9
prop xvel=0.0 yvel=-0.100  spin=0.0 range id 2
fix spin
fix x y range id 1 

set fishcall FC_CYC_MOT apply_vis ; call [apply_vis] just before motion calculation

set dt max 1e-3
set grav 0.0 0.0

hist id=1 nstep=1 ball ypos id=2 

plot add axes black
plot add ball yellow id=on angle on
plot add con blue
plot add hist 1
plot show

cyc 400
;print particle2's y-position
hist write 1  file T1particle2_y_position.dat
;hist write 1  file T2particle2_y_position.dat