```#impactLoadCell.py:
#=================
#
# An ESyS-Particle simulation script to implement breakage
#  of a cylindrical specimen of DEM rock impacted with
#  a spherical indenter a la the laboratory Short Impact Load Cell.
#
# D.Weatherley, UQ, 2018.
#
#Import the relevant modules:
from esys.lsm import *
from esys.lsm.util import *
from math import sqrt
from ball_mover import SphereBodyMoverRunnable

#Fundamental units of measure:
#	length = millimetres (mm)
#	time = milliseconds (ms)
#	force = Newtons (N)

#Derived units of measure:
#	velocity = metres per second (m/s)
#	stress = megaPascals (MPa)
#	mass = grams (g)
#	density = 1.0e-6 kg/m^3 (g/mm^3)
#	acceleration = 1.0e-3 m/^2 (mm/s^2)

#Model Parameters:
#================
cylRadius = 5.0 	#radius of the cylindrical specimen of rock
cylHeight = 10.0	#length of the cylindrical specimen of rock

Rmin = 0.2		#minimum DEM particle radius
Rmax = 0.5		#maximum DEM particle radius

bondModulus = 1.0e5		#elastic modulus of beam interactions
bondPoissonRatio = 0.25		#Poisson's ratio of beam interactions
bondCohesion = 100.		#cohesive strength of beam interactions
bondTanAngle = 1.0		#tan(friction angle) of beam interactions
frictionCoeff = 0.6		#friction coefficient for frictional contacts

density = 2.6e-3		#density of individual DEM particles
gravity = Vec3(0,-9.81e-3,0)	#acceleration due to gravity

dropHeight = 200.0		#initial height of the spherical indenter
ballRadius = 25.0		#radius of the spherical indenter
ballDensity = 7.5e-3		#density of the spherical indenter (i.e. steel)

impactSpeed = sqrt(2.*9.81e-3*200.0)	#calculated impact speed

wallInteractionType = "glued"	#specifies type of wall interaction
#(valid types: "elastic" or "glued")

numT = 100000	#total number of simulation timesteps

#calculated stable timestep increment:
dt = 0.1*sqrt(4.*density*Rmin**3./3./bondModulus/Rmax)

#Construct a simulation container:
sim = LsmMpi (
numWorkerProcesses = 1,
mpiDimList = [1,1,1]
)

#Initialise the neighbour search algorithm (using RotSphere):
sim.initNeighbourSearch (
particleType = "RotSphere",
gridSpacing = 2.5*Rmax,
verletDist = 0.2*Rmin
)

#Specify the number of simulation timesteps and timestep increment:
sim.setNumTimeSteps (numT)
sim.setTimeStepSize (dt)

#Read the initial DEM particle geometry:

#Set the density of the DEM particles:
sim.setParticleDensity (tag = 1, mask = -1, Density = density)
sim.setParticleDensity (tag = 2, mask = -1, Density = density)

#Add a "SphereBody" to represent the spherical indenter:
sim.createSphereBody (
name = "steel_ball",
)

#Add a planar wall beneath the DEM rock specimen to record forces:
sim.createWall (
name = "floor",
posn = Vec3(0,0,0),
normal = Vec3(0,1,0)
)

#Initialise elastic interactions between DEM particles and the indenter:
sim.createInteractionGroup (
NRotElasticSphereBodyPrms (
name = "pb_repel",
sphereName = "steel_ball",
normalK = bondModulus
)
)

#Initialise elastic interactions between DEM particles and the wall:
if (wallInteractionType=="elastic"):
sim.createInteractionGroup (
NRotElasticWallPrms (
name = "pw_floor",
wallName = "floor",
normalK = bondModulus
)
)
elif (wallInteractionType == "glued"):
sim.createInteractionGroup (
NRotSoftBondedWallPrms (
name = "pw_floor",
wallName = "floor",
normalK = bondModulus,
shearK = 0.0,
particleTag = 2,
scaling = True
)
)

#Specify brittle-elastic beam interactions between bonded particles:
sim.createInteractionGroup (
BrittleBeamPrms (
name = "pp_beams",
youngsModulus = bondModulus,
poissonsRatio = bondPoissonRatio,
cohesion = bondCohesion,
tanAngle = bondTanAngle,
tag = 1
)
)

#Specify rotational frictional interactions between unbonded particles:
sim.createInteractionGroup (
FrictionPrms (
name = "pp_friction",
youngsModulus = bondModulus,
poissonsRatio = bondPoissonRatio,
dynamicMu = frictionCoeff,
staticMu = frictionCoeff
)
)

#Ensure each particle pair only undergoes beam or frictional interactions:
sim.createExclusion ("pp_beams","pp_friction")

#Specify the acceleration due to gravity:
sim.createInteractionGroup (
GravityPrms (
name = "gravity",
acceleration = gravity
)
)

#Attach a Runnable to move the indenter at constant speed:
ball_mover = SphereBodyMoverRunnable (
sim,
"steel_ball",
speed = impactSpeed,
dirn = Vec3 (0,-1,0)
)

#Add a checkpointer to record simulation data for visualisation:
sim.createCheckPointer (
CheckPointPrms (
fileNamePrefix = "snapshot",
beginTimeStep = 0,
endTimeStep = numT,
timeStepIncr = 1000
)
)

#Add a FieldSaver to record wall forces:
sim.createFieldSaver (
WallVectorFieldSaverPrms (
wallName = ["floor"],
fieldName = "Force",
fileName = "wall_force.csv",
fileFormat = "RAW_SERIES",
beginTimeStep = 0,
endTimeStep = numT,
timeStepIncr = 1
)
)

#Add a FieldSaver to record total particle kinetic energy each timestep:
sim.createFieldSaver (
ParticleScalarFieldSaverPrms (
fieldName="e_kin",
fileName="Ekin.csv",
fileFormat="SUM",
beginTimeStep = 0,
endTimeStep = numT,
timeStepIncr = 1
)
)

#Add a FieldSaver to record the number of bonded interactions each timestep:
sim.createFieldSaver (
InteractionScalarFieldSaverPrms (
interactionName="pp_beams",
fieldName="count",
fileName="Nbonds.csv",
fileFormat="SUM",
beginTimeStep = 0,
endTimeStep = numT,
timeStepIncr = 1
)
)

#Add a FieldSaver to record the total stored elastic strain energy:
sim.createFieldSaver (
InteractionScalarFieldSaverPrms (
interactionName="pp_beams",
fieldName="potential_energy",
fileName="Epot.csv",
fileFormat="SUM",
beginTimeStep = 0,
endTimeStep = numT,
timeStepIncr = 1
)
)

#execute the simulation:
sim.run()
```
GlossyBlue theme adapted by David Gilbert 