# =========================================================================
# Twisting a pre-buckled straight rod
# =========================================================================

variable    N       equal 300

# Beam geometry (mm)
variable    L       equal 300
variable    l0      equal $(v_L/v_N)
variable    d       equal 3.1

# Cross-section properties
variable    A       equal (PI/4)*v_d^2
variable    I       equal (PI/64)*v_d^4
variable    J       equal (PI/32)*v_d^4

# Elastic moduli (kPa)
variable    E       equal 1300
variable    G       equal 433

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

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

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

neighbor        $(5*v_l0) bin
comm_modify     mode single cutoff $(5*v_l0) vel yes

read_data       rod-N${N}.lam

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

group           l1 id 1
group           l2 id 2
group           r2 id $(v_N)
group           r1 id $(v_N+1)

group           left  union l1 l2
group           right union r1 r2
group           ends  union l1 r1
group           free  subtract all ends
group           mobile union free right

# =========================================================================
# 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 overlay/pair no
bond_coeff      1 ${Kr} ${Ks} ${Kt} ${Kb} &
                -1 -1 -1 -1 &
                1e-4 1e-4 1e-4 1e-4

# =========================================================================
# Granular contact (Hertz-Mindlin)
# =========================================================================

# Made up values
variable        YoungsModulus    equal 10*v_E
variable        PoissonsRatio    equal 0.4
variable        FrictionCoeff    equal 0.1
variable        CoeffRestitution equal 0.7

pair_style      granular limit_damping yes
pair_coeff      * * hertz/material $(v_YoungsModulus) $(v_CoeffRestitution) $(v_PoissonsRatio) &
                tangential mindlin/force NULL 0.5 $(v_FrictionCoeff) &
                damping coeff_restitution

# =========================================================================
# Computes and output
# =========================================================================

compute         mytemp all temp/sphere
compute         kerot  all erotate/sphere

# Note: variable torque and compute torque/chunk only works for group of point particles
compute         ft l1 property/atom fx fy fz tqx tqy tqz
compute         f  l1 reduce sum c_ft[1] c_ft[2] c_ft[3]
compute         t  l1 reduce sum c_ft[4] c_ft[5] c_ft[6]

variable        fx equal c_f[1]
variable        tx equal c_t[1]

compute         xmax    all reduce max x
compute         xmin    all reduce min x
variable        dx      equal "c_xmax-c_xmin"
compute         quat    all property/atom quatw quati quatj quatk

compute         ymaxabs all reduce maxabs y
compute         zmaxabs all reduce maxabs z

variable        cycles equal 0

thermo          10000
thermo_style    custom step ke c_kerot pe v_cycles v_dx c_ymaxabs c_zmaxabs v_fx v_tx
thermo_modify   temp mytemp norm no

# =========================================================================
# Stage 1: Axial compression (80 mm)
# =========================================================================

dump            1 all custom 100000 ./compress-N${N}.dump &
                id diameter x y z mux muy muz fx fy fz tqx tqy tqz c_quat[*]

variable        timestep equal 1e-3
timestep        ${timestep}

variable        nsteps equal 5e6
variable        period equal v_timestep*v_nsteps
variable        vx     equal -80/(v_timestep*v_nsteps)

fix             1 free nve/bpm/sphere update dipole
fix             2 r1 move linear ${vx} 0 0
fix             3 free gravity 9.81e-3 vector 0 -1 0
fix             4 free langevin 0 0 1e3 12345 omega yes

run             $(v_nsteps)

# =========================================================================
# Stage 2: Relaxation
# =========================================================================

unfix           2
fix             4 free langevin 0 0 1e3 12345 omega yes

variable        nsteps equal 1e7
run             $(v_nsteps)

# write_restart   compressed-beam.restart

undump          1
unfix           1
unfix           3
unfix           4

# =========================================================================
# Stage 3: Twist (9 full rotations)
# =========================================================================

variable        timestep equal 1e-5
timestep        ${timestep}
reset_timestep  0

balance         1.2 shift x 10 1.1

dump            1 all custom 10000000 ./twist-N${N}.dump &
                id diameter x y z mux muy muz fx fy fz tqx tqy tqz c_quat[*]

variable        nturns       equal 9
variable        stepsperturn equal 1e9
variable        period       equal v_timestep*v_stepsperturn
variable        cycles       equal step/v_stepsperturn
variable        theta        equal 2*PI*v_cycles

bond_coeff      1 ${Kr} ${Ks} ${Kt} ${Kb} &
                -1 -1 -1 -1 &
                1e-2 1e-2 1e-2 1e-2

fix             1 free nve/bpm/sphere update dipole
fix             2 r1 move rotate 0 0 0 -1 0 0 $(v_period) update dipole
fix             3 free gravity 9.81e-3 vector 0 -1 0
fix             4 free langevin 0 0 1e2 12345 omega yes

fix             5 all print 1000000 &
                "${theta} $(c_ymaxabs) $(c_zmaxabs) ${tx}" &
                file twist-N${N}-theta-yzmaxabs-Tx.txt screen no &
                title "theta [rad], y_maxabs [mm], z_maxabs [mm], Tx [uJ]"

variable        i loop ${nturns}
label           loopstart
run             $(v_stepsperturn)
next            i
jump            SELF loopstart