# =========================================================================
# Cantilever beam under static uniform load
# 
# Replace NBONDS with the desired value before running.
# =========================================================================

variable    N       equal NBONDS

# Beam geometry (dimensionless)
variable    L       equal 100
variable    diam    equal 1
variable    l0      equal $(v_L/v_N)

# Cylindrical cross-section
variable    A       equal (PI/4)*v_diam^2
variable    I       equal (PI/64)*v_diam^4
variable    J       equal (PI/32)*v_diam^4

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

# Load chosen such that max tip deflection = 10% of L
variable    dymax   equal "-0.1*v_L"
variable    w       equal 8*v_E*v_I*v_dymax/v_L^4

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

units           lj
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

change_box      all x final -500 500 y final -500 500 z final -500 500

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

region          left block INF 0.1 INF INF INF INF
region          right block 99.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 &
                1e-2 1e-2 1e-2 1e-2

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

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

timestep        0.01

compute         delx right reduce ave x
compute         dely right reduce ave y
variable        dx equal "c_delx/v_L-1"
variable        dy equal "c_dely/v_L"

# Distribute total load to free atoms
variable        Fy equal v_w*v_L/count(free)

# Theoretical deflection curve (Euler-Bernoulli)
variable        ytheo atom "v_w*x^2/(24*v_E*v_I)*(6*v_L^2-4*v_L*x+x^2)"
variable        ydiff atom "y-v_ytheo"
compute         maxdiffy all reduce maxabs v_ydiff

# Energy
variable        Wexta atom v_w*v_l0*y
compute         Wext all reduce sum v_Wexta
variable        Uf_theo equal $(v_L^5*v_w^2/(40*v_E*v_I))

compute         mytemp all temp/sphere
compute         erot all erotate/sphere
variable        U equal "pe"
variable        Uerr equal abs(v_U-v_Uf_theo)/v_Uf_theo

thermo          1000000
thermo_style    custom step temp ke c_erot pe v_U v_Uf_theo v_Uerr c_Wext v_dy c_maxdiffy
thermo_modify   temp mytemp norm no

dump            1 all custom 1000000 ./bend-uniform-N${N}.dump &
                id x y z

# =========================================================================
# Stage 1: Ramp load
# =========================================================================

fix             nve free nve/bpm/sphere

variable        fy equal v_Fy*ramp(0,1)
fix             adf free addforce 0 v_fy 0
fix_modify      adf energy no

run             50000000

# =========================================================================
# Stage 2: Relax under constant load
# =========================================================================

fix             lan free langevin 0 0 1000 83715463 zero yes omega yes

variable        fy equal v_Fy
fix             adf free addforce 0 v_fy 0
fix_modify      adf energy no

run             50000000

# =========================================================================
# Post-processing
# =========================================================================

variable        Wexta_theo atom v_w*v_l0*v_ytheo
compute         Wextf_theo all reduce sum v_Wexta_theo

variable        fddy atom "v_w/(v_E*v_I)*(0.5*v_L^2-v_L*x+0.5*x^2)"
variable        Uf_atom atom "0.5*v_E*v_I*(v_fddy)^2*v_l0"
compute         Uf_discrete all reduce sum v_Uf_atom
run             0

print "======================================================================"
print "External work:       W = $(c_Wext)"
print "Theoretical:         W = $(c_Wextf_theo)"
print "======================================================================"
print "Elastic energy:      U = ${Uf_theo}"
print "Discrete approx:     U = $(c_Uf_discrete)"
print "======================================================================"

dump            final all custom 1 ./bend-uniform-N${N}-final.dump &
                id x y z
run             0