# =========================================================================
# Cantilever beam under sinusoidal point load on free end
# 
# Replace NBONDS with the desired value before running.
# =========================================================================

variable    N       equal NBONDS

# Beam geometry
variable    L       equal 10.
variable    l0      equal $(v_L/v_N)
variable    w       equal 0.5
variable    h       equal 0.25

# Rectangular cross-section
variable    A       equal v_w*v_h
variable    I       equal v_w*v_h^3/12
variable    J       equal v_w*v_h/12*(v_w^2+v_h^2)

# Material properties
variable    E       equal 210e9
variable    nu      equal 0.3
variable    G       equal "v_E/(2*(1+v_nu))"

# Loading
variable    P0      equal 10e6
variable    omega   equal 50

# =========================================================================
# Initialization
# =========================================================================

units           si
dimension       3
boundary        f f f
processors      * 1 1

atom_style      bpm/sphere
special_bonds   lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
newton          on off

neighbor        0.5 bin
neigh_modify    once yes
comm_modify     mode single cutoff $(2*v_l0) vel yes

read_data       beam-${N}.lam

# =========================================================================
# Groups
# =========================================================================

region          left block INF 0.1 INF INF INF INF
region          right block 9.9 INF INF INF INF INF
group           left region left
group           right region right
group           free subtract all left

# =========================================================================
# BPM bond style
# =========================================================================

variable        Kr equal v_E*v_A/v_l0
variable        Ks equal 12*v_E*v_I/v_l0^3
variable        Kt equal v_G*v_J/v_l0
variable        Kb equal v_E*v_I/v_l0

bond_style      bpm/rotational smooth no break no
bond_coeff      1 ${Kr} ${Ks} ${Kt} ${Kb} &
                -1 -1 -1 -1 &
                0 0 0 0

pair_style      zero $(2*v_l0)
pair_coeff      * * $(2*v_l0)

# =========================================================================
# Simulation time and output
# =========================================================================

variable        timestep equal 1e-8
variable        runtime  equal 1.5
variable        nsteps   equal v_runtime/v_timestep
variable        time     equal step*v_timestep

timestep        ${timestep}

compute         tx right reduce ave x
compute         ty right reduce ave y
variable        dx equal "c_tx/v_L-1"
variable        dy equal "c_ty/v_L"

compute         mytemp all temp/sphere
compute         erot all erotate/sphere

thermo          $(v_nsteps/1000)
thermo_style    custom step temp v_time v_dx v_dy
thermo_modify   temp mytemp

dump            1 all custom $(v_nsteps/100) ./bend-point-N${N}.dump &
                id x y z fx fy fz

fix             prnt all print $(v_nsteps/1000) &
                "${time} ${dx} ${dy}" &
                file deflection-N${N}.txt screen no &
                title "t dx/L dy/L"

# =========================================================================
# Dynamics
# =========================================================================

fix             nve free nve/bpm/sphere

variable        p equal v_P0*sin(v_omega*step*v_timestep)
fix             adf right addforce 0 v_p 0
fix_modify      adf energy no

run             ${nsteps}