∎

11institutetext: Andrew D. Davis 22institutetext: Cold Regions Research and
Engineering Laboratory, 72 Lyme Rd, Hanover, NH 03755
NYU—Courant Institute of Mathematical Sciences, 251 Mercer Street, New York,
NY 10012
22email<EMAIL_ADDRESS>33institutetext: Brendan A. West
44institutetext: Cold Regions Research and Engineering Laboratory, 72 Lyme Rd,
Hanover, NH 03755 55institutetext: Nathanael J. Frisch 66institutetext: Cold
Regions Research and Engineering Laboratory, 72 Lyme Rd, Hanover, NH 03755
77institutetext: Devin T. O’Connor 88institutetext: Cold Regions Research and
Engineering Laboratory, 72 Lyme Rd, Hanover, NH 03755 99institutetext: Matthew
D. Parno 1010institutetext: Cold Regions Research and Engineering Laboratory,
72 Lyme Rd, Hanover, NH 03755

# ParticLS: Object-oriented software for discrete element methods and
peridynamics

Andrew D. Davis  Brendan A. West  Nathanael J. Frisch  Devin T. O’Connor
Matthew D. Parno

(Received: date / Accepted: date)

###### Abstract

ParticLS (_Partic_ le _L_ evel _S_ ets) is a software library that implements
the discrete element method (DEM) and meshfree methods. ParticLS tracks the
interaction between individual particles whose geometries are defined by level
sets capable of capturing complex shapes. These particles either represent
rigid bodies or material points within a continuum. Particle-particle
interactions using various contact laws numerically approximate solutions to
energy and mass conservation equations, simulating rigid body dynamics or
deformation/fracture. By leveraging multiple contact laws, ParticLS can
simulate interacting bodies that deform, fracture, and are composed of many
particles. In the continuum setting, we numerically solve the peridynamic
equations—integro-differential equations capable of modeling objects with
discontinuous displacement fields and complex fracture dynamics. We show that
the discretized peridynamic equations can be solved using the same software
infrastructure that implements the DEM. Therefore, we design a unique software
library where users can easily add particles with arbitrary geometries and new
contact laws that model either rigid-body interaction or peridynamic
constitutive relationships. We demonstrate ParticLS’ versatility on test
problems meant to showcase features applicable to a broad selection of fields
such as tectonics, granular media, multiscale simulations, glacier calving,
and sea ice.

###### Keywords:

Discrete element method Peridynamics Level set method Optimization Contact
mechanics Fracture dynamics Meshfree methods Geophysical modeling

††journal: ArXiv

##  1 Introduction

Small-scale processes, such as the interaction of sand grains or fracture
dynamics, impact the large-scale behavior of many materials. We devise a new
software library that numerically simulates such phenomena using meshless
methods or the discrete element method (DEM). These techniques solve equations
derived from conservation of mass and momentum and constitutive models
relating material deformation and stress. Meshless methods, like peridynamics
Sillingetal2007 , smoothed-particle hydrodynamics Monaghan1992 , or
reproducing kernel particle methods Bessaetal2014 , use partial differential
equations or integro-differential equations to describe the material behavior.
Alternatively, the DEM explicitly models small-scale behavior through the
motion and interaction of discrete particles (e.g., sand grains)
Andradeetal2012GranularElement ; CundallStrack1979 ; GuoCurtis2015 ;
Hopkins1991 ; Hopkins2004 . Although DEM and meshless methods have
fundamentally different material representations, we use abstract
computational structures to implement both methods in a single software
library, called ParticLS (_Partic_ le _L_ evel _S_ ets). We will refer to both
the DEM and peridynamics as “meshfree” techniques to highlight their similar
computational structure. ParticLS is uniquely designed as an object-oriented
software library that leverages software concepts such as abstraction and
inheritance to organize and structure the code. Therefore, ParticLS provides
core functionality that users can use to simulate complex geophysical systems
such as sea ice or soils.

The macroscale properties of granular media require simulating computationally
expensive interactions between irregularly shaped particles Choetal2006 ;
GuoMorgan2004 . Methods in the literature for representing complex shapes
include clumping spherical particles Garciaetal2009 , 3-dimensional polyhedra
Lathametal2001 ; LathamMunjiza2004 , tablet shaped particles Songetal2006 , or
axisymmetrical particles Favieretal1999 . However, computing particle-particle
contact information (e.g., the location where two particles collide) becomes
computationally expensive given many particles or complex geometries
Johnsonetal2004 ; Songetal2006 . Alternatively, irregular particle boundaries
can be represented using the zero level set of a signed distance function
(SDF) Kawamotoetal2016 ; Kawamotoetal2018 .

ParticLS uses this level set approach, as formulated by Kawamotoetal2016 ;
Kawamotoetal2018 , which allows us to formulate a novel contact detection
scheme as a constrained optimization problem. The contact point between a home
particle and its neighbor is a point on the home particle’s boundary that
minimizes the neighbor particle’s SDF. This generalized approach allows
ParticLS to represent many unique particle geometries within one simulation
without having to specify different contact models for every potential
particle geometry combination. Users only need to define an SDF to add a new
particle geometry that works with this contact detection algorithm.

ParticLS can also model deformable bodies using peridynamics, which can
simulate fracturing and discontinuous materials Fosterestal2010 ;
MadenciOterkus2016 . Peridynamics replaces the divergence operator on the
Cauchy stress with an integral operator; these equations are valid even along
crack tips and discontinuities where the spatial derivatives used to compute
stress in traditional continuum formulations generate singularities
Askarietal2008 ; Rogula1982 ; Silling2000 ; Sillingetal2007 ;
SillingLehoucq2010 ; SillingLehoucq2008 . Peridynamic implementations share
many features with the DEM and we leverage this similarity to provide a
software framework for both approaches.

ParticLS’s object-oriented framework mimics mathematical formulations,
allowing users to easily add functionality specific to their application. In
this paper, we highlight four major contributions of ParticLS:

  1. 1.

We address computational limitations when computing contact between particles
that have complex geometries by deriving an optimization-based contact
detection algorithm between neighboring particles.

  2. 2.

Generic meshfree algorithms (e.g., contact detection) are easily adapted to
specific applications using abstract base classes that provide default
implementations for important functionality. Users can easily add to the
library by inheriting the functionality of these classes and prescribing their
own definitions.

  3. 3.

Users can couple DEM and peridynamic methods with other time-dependent
simulations or data. ParticLS uses state variables that are defined by an
ordinary differential equation, such as temperature-dependent material
properties or contact shear forces.

  4. 4.

Deformable bodies can interact by implementing DEM and peridynamic models with
the same code.

Section 2 describes the model formulation for both the DEM and peridynamics,
specifically highlighting how the discretized peridynamic equations are
analogous to a DEM. Section 3 describes the optimization-based contact
detection algorithm. Section 4 describes ParticLS’ object-oriented framework
and Section 5 shows selected examples to demonstrate ParticLS’ unique
capabilities.

##  2 DEM and peridynamic background

###  2.1 Discrete element method (DEM)

The equations for the DEM are derived from force balance and conservation of
angular momentum CundallStrack1979 . ParticLS’s DEM implementation evolves the
global position and orientation of a specific particle relative to a reference
position. We use lowercase symbols and letters to denote quantities in the
global frame and uppercase symbols to represent quantities in the local
reference frame. Let $X$ denote the $d\in\\{2,3\\}$ dimensional coordinate in
a reference frame such that the particle’s center of mass is at $X=0$. In this
frame, a particle’s geometry is defined by a SDF
$\varphi:\mathbb{R}^{d}\mapsto\mathbb{R}$ Kawamotoetal2016 ; Kawamotoetal2018
such that

$\varphi(X)\begin{cases}&<0\mbox{ if $X$ is inside the particle}\\\ &=0\mbox{
if $X$ is on the particle's boundary}\\\ &>0\mbox{ if $X$ is outside the
particle.}\end{cases}$ (1)

Furthermore, the magnitude $|\varphi(X)|$ defines the minimum distance between
$X$ and a point on the boundary. For the purpose of this paper, we assume that
the particles are convex. Let $\bar{y}_{i}(t)$ be the center of mass of
particle $i\in\\{1,2,\ldots,n\\}$ in a global coordinate frame at time $t$ and
$u_{i}(t)$ be the displacement at time $t$ such that $u_{i}(0)=0$ and
$\bar{y}_{i}(t)=\bar{x}_{i}+u_{i}(t)$, where $\bar{x}_{i}$ is the center of
mass at $t=0$. We employ quaternions to represent the orientation of the
particles, which can be defined in terms of a unit orientation vector
$\chi_{i}(t)$ ($\|\chi_{i}(t)\|_{2}=1$) and angle $\theta_{i}(t)$

$\begin{split}\mathpzc{q}_{i}(t)=&\cos{\left(\frac{\theta_{i}(t)}{2}\right)}+\\\
&\left(\chi_{i}^{(0)}(t)\hat{\imath}+\chi_{i}^{(1)}(t)\hat{\jmath}+\chi_{i}^{(2)}(t)\hat{k}\right)\sin{\left(\frac{\theta_{i}(t)}{2}\right)},\end{split}$
(2)

where $\chi_{i}^{(k)}(t)$ is the $k^{\text{th}}$ component of $\chi_{i}(t)$,
such that rotating the particle about $\chi_{i}(t)$ by $\theta_{i}(t)$ radians
aligns the global and reference coordinate frame. The transformation

$X(x,\bar{y}_{i}(t),\mathpzc{q}_{i}(t))=\mathpzc{q}_{i}(t)\circ(x-\bar{y}_{i}(t))$
(3)

maps the global coordinate $x$ into the $i^{\text{th}}$ particle’s local
reference frame. The $\circ$ operator is the action of rotating the vector $x$
by the quaternion $\mathpzc{q}$.

The particles are translated and rotated in the global frame according to
Newton’s second law and conservation of angular momentum. Given the
translational velocity $v_{i}(t)=\dot{u}_{i}(t)$ where the dot denotes time
derivatives,

$\displaystyle\dot{u}_{i}(t)$ $\displaystyle=$ $\displaystyle v_{i}(t)$ (4a)
$\displaystyle m_{i}\ddot{u}_{i}(t)=m_{i}\dot{v}_{i}(t)$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{n}f_{ij}(t)+b_{i}(t),$ (4b)

where $m_{i}$ is the mass of particle $i$, $f_{ij}(t)$ is the force acting on
particle $i$ from particle $j$, and $b_{i}(t)$ are body forces acting on
particle $i$ (e.g., gravity). For particle $i$, let the vector $\omega_{i}(t)$
be the angular velocity, which defines the pure quaternion
$\mathpzc{w}_{i}(t)=(0,\omega_{i}(t))$. The angular velocity describes the
change in orientation

$\dot{\mathpzc{q}}_{i}(t)=\frac{1}{2}\mathpzc{w}_{i}(t)*\mathpzc{q}_{i}(t),$
(5a) and evolves via conservation of angular momentum
$I_{i}(t)\dot{\omega}_{i}(t)=\sum_{j=1}^{n}\tau_{ij}(t)-\dot{I}_{i}(t)\omega_{i}(t),$
(5b)

where $I(t)$ and $\dot{I}_{i}(t)$ is the moment of inertia for particle $i$ in
the global frame and its time derivative and $\tau_{ij}(t)$ is the torque
acting on particle $i$ from particle $j$. In two dimensions, the rotational
equations (5b) significantly simplify because the particle is restricted to
the $x$-$y$ plane, and therefore, we only need the orientation angle
$\theta_{i}(t)$. The angular velocity is now also a scalar quantity such that
(5a) is replaced by $\dot{\theta}_{i}(t)=\omega_{i}(t)$.

###  2.2 Peridynamic formulation

In the peridynamic framework, non-local forces acting across finite distances
allow discontinuities and cracks to form within continuous materials
Silling2000 ; Sillingetal2007 . Let $\mathcal{B}_{\delta}(x)$ be a ball with
radius $\delta$ centered at $x$. Forces are modeled as bonds between a point
$x$ and points in its horizon $\mathcal{B}_{\delta}(x)$, as illustrated in
Figure 1.

Figure 1: Illustration of the non-local forces acting on the middle particle
in a peridynamic setting. Each circle represents the material point
$\mathcal{B}_{r}(x^{\prime})$. The gray region corresponds to the horizon
$\mathcal{B}_{\delta}(x)$ of the middle particle. Each arrow in the right
image corresponds to the pair-wise force per volume squared
$\eta_{i}(x,x_{i},t)$ contributing to the force on the middle particle after
discretization.

Given a material with density $\rho(x)$, the peridynamic equation is

$\rho(x)\ddot{u}(x,t)=\int_{\mathcal{B}_{\delta}(x)}\eta(x,x^{\prime},t)\,dx^{\prime}+\beta(x,t),$
(6)

where $u(x,t)$ is the unknown displacement, $\eta(x,x^{\prime},t)$ is the
pairwise force per volume squared, and $\beta(x,t)$ is the body force per
volume. Mechanical restrictions require that

$\displaystyle\eta(x,x^{\prime},t)=-\eta(x^{\prime},x,t),$ (7a)
$\displaystyle\eta(x,x^{\prime},t)=0\mbox{ when }\|x-x^{\prime}\|_{2}>\delta,$
(7b)

and

$\int_{\Omega}y(x,t)\times\int_{\mathcal{B}_{\delta}(x)}\eta(x,x^{\prime},t)\,dx^{\prime}\,dx=0,$
(8)

where $y(x,t)=x+u(x,t)$ is the deformation Silling2000 ; Sillingetal2007 .
These restrictions are sufficient conditions to conserve linear and angular
momentum. To discretize, we first integrate over a ball centered at material
point $x_{i}$ with radius $r\ll\delta$

$\begin{split}\int_{\mathcal{B}_{r}(x_{i})}\rho(x)&\ddot{u}(x,t)\,dx=\\\
&\int_{\mathcal{B}_{r}(x_{i})}\int_{\mathcal{B}_{\delta}(x)}\eta(x,x^{\prime},t)\,dx^{\prime}+\beta(x,t)\,dx.\end{split}$
(9)

Let the density $\rho_{i}=\rho(x_{i})$ and assume that it is approximately
constant in $\mathcal{B}_{r}(x_{i})$. We also assume that the discrete force
per volume squared $\eta_{ij}(t)=\eta(x_{i},x_{j},t)$ is approximately
constant with respect to the position $x$ (or $x^{\prime}$) in the ball
$\mathcal{B}_{r}(x)$ (or $\mathcal{B}_{r}(x^{\prime})$) and that
$\eta_{ij}(t)=0$ if $\|x_{i}-x_{j}\|>\delta$. Similarly, we assume that the
acceleration $\ddot{u}_{i}(t)=\ddot{u}(x_{i},t)$ and body force density
$\beta_{i}(t)=\beta(x_{i},t)$ are constant in $\mathcal{B}_{r}(x_{i})$. The
discrete peridynamic equation is, therefore,

$\rho_{i}V_{r}\ddot{u}_{i}(t)\approx\sum_{j=1}^{n}\eta_{ij}(t)V_{r}^{2}+\beta_{i}(t)V_{r},$
(10)

where $V_{r}$ is the volume of a $r$-radius ball. Defining
$m_{i}=\rho_{i}V_{r}$, $f_{ij}(t)=\eta_{ij}(t)V_{r}^{2}$, and
$b_{i}(t)=\beta_{i}(t)V_{r}$ derives the peridynamic analog to (4). The
restriction on $\eta(x,x^{\prime},t)$ requiring (8) ensures that the system
conserves angular momentum and, therefore, we do not need to explicitly evolve
the orientation of each particle Sillingetal2007 .

####  2.2.1 State-based peridynamics

ParticLS implements state-based peridynamics, which is an extension of the
formulation above but a detailed description is beyond the scope of this
paper—see Sillingetal2007  for details. Briefly, state-based peridynamic
theory defines the pairwise force per squared volume $\eta(x,x^{\prime},t)$
via a force vector state $F[x,x^{\prime},t]\langle x^{\prime}-x\rangle$, which
is a potentially nonlinear and discontinuous operator $F[x,x^{\prime},t]$ that
acts on the bond $\xi=x^{\prime}-x$ Sillingetal2007 . State-based peridynamics
is the special case of (6) with

$\eta(x,x^{\prime},t)=F[x,x^{\prime},t]\langle
x^{\prime}-x\rangle-F[x^{\prime},x,t]\langle x-x^{\prime}\rangle.$ (11)

We assume forces only act across bonds, therefore, define

$Y[x,x^{\prime},t]\langle
x^{\prime}-x\rangle=\frac{y(x^{\prime},t)-y(x,t)}{\|y(x^{\prime},t)-y(x,t)\|_{2}}$
(12)

to be the unit vector state that points across the deformed bond
$y(x,t)-y(x^{\prime},t)=(x^{\prime}+u(x^{\prime},t))-(x+u(x,t))$. The force
state of ordinary materials is, therefore, written in terms of a scalar state
$f(x,x^{\prime},t)$

$F[x,x^{\prime},t]\langle x^{\prime}-x\rangle=f(x,x^{\prime},t)Y[x,t]\langle
x^{\prime}-x\rangle.$ (13)

Given $f(x,x^{\prime},t)$, (13) defines the force state that is substituted
into (11).

##  3 Optimization-based contact detection

As introduced in Section 1, ParticLS implements a novel contact detection
routine that is formulated as an optimization problem based on the particles’
SDFs. The LS-DEM, developed by Kawamotoetal2016 , stores discrete points along
a particle’s boundary (zero level set surface). Contact is detected by looping
over these boundary points and determining if a neighbor’s SDF is negative for
each boundary point. This approach has two potential drawbacks: (i) it is
computationally expensive to iterate through each set of boundary points for
simulations with many particles and (ii) the contact force between the two
particles can be sensitive to the number of points that are placed on the
particle boundary li2019capturing . Instead, ParticLS formulates an
optimization problem that finds the point on the “home” particle that
penetrates deepest into its neighbor. This contact routine is defined
abstractly in terms of the SDF so that new particle geometries are easily
added by defining the corresponding SDF.

Consider two particles defined by their centers of mass ($\bar{x}$ and
$\bar{x}^{\prime}$), orientations ($\mathpzc{q}$ and $\mathpzc{q}^{\prime}$),
and SDFs ($\varphi(X(x,\bar{x},\mathpzc{q})$ and
$\varphi^{\prime}(X(x,\bar{x}^{\prime},\mathpzc{q}^{\prime})$); denote these
particles $p=(\bar{x},\mathpzc{q},\varphi)$ and
$p^{\prime}=(\bar{x}^{\prime},\mathpzc{q}^{\prime},\varphi^{\prime})$. We
refer to a ‘point of contact’ as the location on one particle’s boundary that
is closest to the neighboring particle, or the deepest point of penetration if
the particles are penetrating each other. A potential point of contact is
located at a point $x_{1}$ on the boundary of particle $p$ that is closest to
particle $p^{\prime}$.111Although we have assumed that particles are convex,
this method is compatible for non-convex particles with minimal alterations.
In the non-convex case, multiple contact points are possible and each
potential contact point corresponds to a local minimum of Equation (14). We
can identify these local solutions by repeatedly solving the optimization
problem with different initial guesses. The details of the non-convex case are
beyond the scope of this paper and we, therefore, leave further discussion to
future work. There is a corresponding point $x_{2}$ on the boundary of
particle $p^{\prime}$ that minimizes the distance to particle $p$.
Mathematically, these contact points solve the optimization problem

$\displaystyle\underset{x_{1},x_{2}\in\mathbb{R}^{d}}{\text{minimize}}$
$\displaystyle\varphi(X(x_{2},\bar{x},\mathpzc{q}))+\varphi^{\prime}(X(x_{1},\bar{x}^{\prime},\mathpzc{q}^{\prime}))$
(14) subject to $\displaystyle\varphi(X(x_{1},\bar{x},\mathpzc{q}))=0$
$\displaystyle\varphi^{\prime}(X(x_{2},\bar{x}^{\prime},\mathpzc{q}^{\prime}))=0.$

In principle, we can use any standard optimization algorithm to solve this
problem. In some cases, such as sphere-sphere contact, this problem has a
known analytical solution. In others, we can leverage the specific form of the
signed distance function to devise a specialized algorithm; for example, the
Gauss-Seidel can efficiently solve this problem for polygon-polygon contact.
More generally, when we do not have an analytical solution or specialized
algorithm we solve this optimization problem using the alternating direction
method of multipliers (ADMM)—also called Douglas-Rachford splitting
Benningetal2015 ; Goldsteinetal2014 ; ParikhBoyd2014 . Importantly, ADMM does
not require derivative information of the SDF, instead relying on the proximal
operator, which allows particle geometry to have sharp corners. Since the SDF
is continuous by definition, this algorithm does not impose any restrictions
on the particle geometry. A detailed description of these algorithms is beyond
the scope of this paper.

Previous “brute-force” contact detection algorithms, which stores $n$ points
along the zero level set, are often significantly more computationally
expensive than solving the optimization problem posed in Equation (14).
Specifically, in the brute force algorithm, the “home” particle needs to
evaluate its neighbor’s SDF at all $n$ points on its boundary. In comparison,
ADMM typically converges very quickly. We also use the solution from the
previous timestep as an initial condition, which often leads to convergence in
1-2 iterations. Additionally, ADMM finds the deepest penetration point up to a
user-specified tolerance, rather than approximating this point using a pre-
defined set along the particle boundary. Finally, storing the required $n$
boundary points for every particle in the simulation can require a huge amount
of computer memory for simulations with many thousands of particles.

For convex particles, the points $x_{1}$ and $x_{2}$ are unique ParikhBoyd2014
and represent potential contact locations. The particles are penetrating if
both distance functions are negative.222Note that the SDFs will always have
the same sign. When this occurs, it is useful to look at the bond $\xi$, which
is the vector pointing from the first particle to the second, as shown in
Figure 2. Mathematically, the bond is defined by

$\xi=x_{2}-x_{1}$ (15a) and the bond length is $\|\xi\|_{2}$. When the
particles are penetrating, the bond length is the penetration distance. Given
the translational velocities $v(t)$ and $v^{\prime}(t)$, the relative velocity
is
$v_{r}(t)=(v(t)+\omega\times(x-x_{1}))-(v^{\prime}(t)+\omega^{\prime}\times(x^{\prime}-x_{2}))$
(15b) and the tangential velocity, as defined by Kawamotoetal2016 , is
$v_{a}(t)=v_{r}(t)-(v_{r}(t)\cdot\xi)\xi/\|\xi\|_{2}.$ (15c)

The torque and force between the two particles depends on the contact summary
as defined by Equation (15c). For non-convex particles, we divide their
geometry into convex sub-particles and then find the potential contact points
between each sub-particle and the neighbor. The resultant set of points
describes potential contact points for the non-convex particles, and allows
for multiple contact points, which we expect to happen in certain non-convex
contact scenarios.

$x_{1}$$x_{2}$$x_{1}$$x_{2}$Particle $p=(\bar{x},\mathpzc{q},\varphi)$Particle
$p^{\prime}=(\bar{x}^{\prime},\mathpzc{q}^{\prime},\varphi^{\prime})$ Figure
2: The red and blue curves correspond to
$\varphi(X(x_{1},\bar{x},\mathpzc{q}))=0$ and
$\varphi^{\prime}(X(x_{1},\bar{x}^{\prime},\mathpzc{q}^{\prime}))=0$ level
sets of two particles $p=(\bar{x},\mathpzc{q},\varphi)$ and
$p^{\prime}=(\bar{x}^{\prime},\mathpzc{q}^{\prime},\varphi^{\prime})$. The
arrows represent the bond $\xi=x_{1}-x_{2}$.

###  3.1 Example: elastic contact

Elastic contact highlights an example where the force and torque are
parameterized by the contact summary (15c) CundallStrack1979 . Given two
particles $p(t)=(\bar{x}(t),\bar{\mathpzc{q}}(t),\varphi)$ and
$p^{\prime}(t)=(\bar{x}^{\prime}(t),\bar{\mathpzc{q}}^{\prime}(t),\varphi^{\prime})$,
the normal force is

$F_{n,e}(t)=-\kappa_{n}\llbracket\varphi^{\prime}(X(x_{1},\bar{x}^{\prime}(t),\mathpzc{q}^{\prime}(t)))<0\rrbracket\xi,$
(16)

where $\llbracket\cdot\rrbracket$ is the indicator function. The shear force
evolves according to

$\dot{F}_{s,e}(t)=-\kappa_{s}\llbracket\varphi^{\prime}(X(x_{1},\bar{x}^{\prime}(t),\mathpzc{q}^{\prime}(t)))<0\rrbracket
v_{a}(t)$ (17)

with $F_{s,e}(0)=0$. If the particles are not in contact then we reset
$F_{s,e}=0$—this is important if the particles ever deflect away from each
other then come back into contact. A Coulomb friction criterion is applied to
the shear force according to

$F^{c}_{s,e}(t)=\left\\{\begin{array}[]{ll}F_{s,e}(t)&\text{if
}\|F_{s,e}(t)\|\leq\mu\|F_{n,e}(t)\|\\\ \mu
F_{n,e}(t)&\text{otherwise}\end{array}\right.$ (18)

where $\mu$ is the friction coefficient. The total force acting on particle
$p$ due to particle $p^{\prime}$ is therefore $F(t)=F_{n,e}(t)+F^{c}_{s,e}(t)$
and the torque is $\tau(t)=((x_{1}+x_{2})/2-\bar{x})\times F(t)$. In this
example, we define the moment arm $(x_{1}+x_{2})/2$ to be halfway between the
points that solve (14). In two dimensions, this simplifies by restricting all
translation and rotation to the $x$-$y$ plane.

###  3.2 Example: viscoelastic contact

We also highlight linear viscoelastic contact, which can be conceptualized as
a linear spring-dashpot system Shafer1996 . This model is similar to the
elastic contact model but with additional viscous damping components that
remove energy from particle interactions. The normal force due to damping is

$F_{n,v}(t)=-\gamma_{n}v_{r}(t),$ (19)

where $\gamma_{n}$ is the normal viscous damping coefficient. Therefore the
total normal force is $F_{n}(t)=F_{n,e}(t)+F_{n,v}(t)$. The same evolution
equation (17) models the shear force, however, viscoelastic contact adds an
additional shear force component that accounts for damping:

$F_{s,v}(t)=-\gamma_{s}v_{a}(t),$ (20)

where $\gamma_{s}$ is the tangential viscous damping coefficient. Therefore,
the total normal force is $F_{s}(t)=F_{s,e}(t)+F_{s,v}(t)$. The total shear
force applies the same Coulomb friction criterion in (18). The torque remains
$\tau(t)=((x_{1}+x_{2})/2-\bar{x})\times F(t)$.

##  4 ParticLS object-oriented implementation

ParticLS defines abstract objects that dictate an interface to important
functionality for meshfree methods. Through inheritance and abstraction users
can easily customize base functionality to implement new particle geometries,
body forces, contact forces, etc.. This makes ParticLS a useful tool for a
wide variety of applications, e.g. soil mechanics, fracturing bodies, and/or
sea ice. We refer to objects/data structures as “classes” and denote them with
the typewritter font.

ParticLS evolves time-dependent quantities of interest by evolving a
collection of abstract objects called StateObjects, which are stored in the
Simulation. The simulation-level state is the combined state of its
StateObjects. ParticLS compartmentalizes Simulation into (i) a collection of
time-dependent StateObjects, (ii) numerical algorithms required to evolve the
states (finding particle nearest neighbors and time-integration), and (iii)
observers that extract quantities of interest at specified times.

###  4.1 Collections of time-dependent state objects

Unlike other particle methods, this allows ParticLS to additionally evolve
non-particle states—for example, time-dependent geometries. Each StateObject
stores a time-dependent parameter, such as a rigid body’s position and
velocity or the shear force defined by (17) and (20). Rather than explicitly
integrating the shear force, many algorithms use first-order finite difference
to approximate an explicit expression for $F_{s}$. However, if we use a
higher-order method to integrate particle position, the time-discretization of
the shear force must also be updated—otherwise we lose higher-order accuracy
in time. Rather than tie expressions for time-dependent quantities to their
discretization, we incorporate them into the time-dependent state of the
Simulation object via the StateObject structure.

Rigid particles (material points in the peridynamic setting) are implemented
as RigidBodys, which inherit from StateObject. RigidBody stores the

  1. 1.

center of mass (a $d$ dimensional vector $x_{i}(t)$),

  2. 2.

reference particle geometry (a pointer to a Geometry object),

  3. 3.

translational velocity (a $d$ dimensional vector $v_{i}(t)$),

  4. 4.

orientation (an angle $\theta_{i}(t)$ in two dimensions and a quaternion
$\mathpzc{q}_{i}(t)$ in three dimensions), and

  5. 5.

angular velocity (a scalar in two dimensions and a three dimensional vector in
three dimensions $\omega_{i}(t)$).

The global coordinate $x$ is transformed into the particle’s reference
coordinate $X$ via (3) using the center of mass and orientation. Given the
reference coordinate $X$, the RigidBody evaluates its Geometry’s SDF (1) to
determine its shape. The Geometry object also computes the particle’s mass and
moment of inertia, given density $\rho(X)$. The RigidBody object also contains
vectors of BodyForces and ContactForces, which represent forces acting on the
particle. Given the net force, torque, translational velocity, and angular
velocity, RigidBody computes the time-derivative of its state using (4) and
(5b). The Simulation object storing the RigidBodys uses their time derivatives
to evolve in time. Simulation also evolves other time-dependent StateObjects,
such as shear forces stored in ContactForces.

###  4.2 Algorithms for evolving simulations

####  4.2.1 Finding nearest neighbors with $k$-d trees

The Searcher object finds the nearest neighbors for each particle in the
simulation. We define “nearest neighbors” with a distance metric between
particles $d(p_{i},p_{j})$ and say $p_{i}$ and $p_{j}$ are neighbors if
$d(p_{i},p_{j})<h$. For reference, the worst case algorithm is to first loop
over every particle, and for each of these ‘outer’ particles we again loop
over all of the particles and decide if they are neighbors. The overall
complexity of the algorithm is therefore $O(n^{2})$. Instead, we organize the
particles in a $k$-d tree to significantly reduce the complexity of finding
the nearest neighbors. Assuming
$d(p_{i},p_{j})=\|\bar{x}_{i}-\bar{x}_{j}\|_{2}$, a $k$-d tree stores each
particle as a node in a binary tree Bentley1975 . ParticLS uses the open
source software package nanoflann to construct $k$-d trees nanoflann . On
average, the complexity of constructing $k$-d trees is $O(n\log{(n)})$ and
finding the nearest neighbors is $O(\log{(n)})$ per particle. Therefore, the
overall complexity of constructing the tree and finding each particle’s
nearest neighbors is $O(2n\log{(n)})$. We note that this complexity estimate
assumes that particles are relatively evenly distributed throughout the domain
or that the search finds the $k$ nearest neighbors (rather than the neighbors
within radius $\Delta$). It is also important to note that the $k$-d tree
approach is limited by the Euclidean distance metric between points.

####  4.2.2 Integrating simulation state in time

ParticLS implements second order Runge-Kutta for all StateObjects except
RigidBodys, which use the second order velocity Verlet. Given timestep $\Delta
t$ and any StateObject that is not a RigidBody, the time discretization is

$\displaystyle s_{t+\frac{1}{2}}$ $\displaystyle=$ $\displaystyle
s_{t}+\frac{\Delta t}{2}f(s_{t},t\Delta t)$ (21a) $\displaystyle s_{t+1}$
$\displaystyle=$ $\displaystyle s_{t+\frac{1}{2}}+\frac{\Delta
t}{2}f(s_{t+\frac{1}{2}},t\Delta t+\Delta t/2),$ (21b)

where $f(s_{t},\tau)$ is the result of computing the StateObject’s right hand
side at state $s_{t}$ and time $\tau$. The RigidBody’s right hand side
computes the time derivatives of $v_{t}$ and $\omega_{t}$, $\dot{v}_{t}$ and
$\dot{\omega}_{t}$, respectively. This defines the velocity and acceleration
of the RigidBody and, therefore,

$x_{t+1}=x_{t}+v_{t}\Delta t+\dot{v}_{t}\Delta t^{2}/2$ (22)

is a second-order discretization. Similarly, we use $\omega_{t}$ and
$\dot{\omega}_{t}$ to define a second order scheme to update the orientation
quaternion lim2014granular

$\displaystyle\mathpzc{q}_{t+1}=\frac{\Delta
t}{2}(\mathpzc{w}_{t}+\frac{\Delta t}{2}\dot{\omega}_{t})*\mathpzc{q}_{t},$
(23)

where $\omega_{t}$ is the angular velocity and $\mathpzc{w}_{t}$ is the
corresponding pure quaternion.

###  4.3 Observers and monitoring

ParticLS provides functionality for tracking different quantities of interest
throughout a simulation. This capability makes it easy for users to measure
time-series of a specific particle’s properties, or quantities related to the
whole particle collection. The base class for this functionality is the
Observer class. The ParticleTracker, for example, tracks the state variables
of a RigidBody, which allows the simulation to extract the position, velocity,
orientation, and/or angular velocity of a particular particle-of-interest.
Similarly, the CollectionTracker computes quantities that depend on the entire
collection. For example, in the DEM setting, we can use the CollectionTracker
observer to measure the homogenized Cauchy stress tensor of the collection:

$\sigma=\frac{1}{V}\sum_{i=1}^{n}\sum_{j=1}^{n}d_{ij}\otimes F_{ij},$ (24)

where $V$ is the volume of the smallest axis-aligned bounding box containing
all of the particles, $n$ is the number of particles in the collection,
$d_{ij}=\bar{x}_{j}-\bar{x}_{i}$ is the vector pointing from particle $i$’s
center of mass to particle $j$’s center of mass, and $F_{ij}$ is the force
vector acting between particles $i$ and $j$ GuoZhao2014 . The operator
$\otimes$ denotes the outer product. We note that $F_{ij}=0$ for most pairs
and that (24) can be computed much more efficiently by replacing the second
summation with the nearest neighbors of particle $i$. In practice, this is how
we compute the Cauchy stress tensor. The Observer class and its children
(e.g., ParticleTracker and CollectionTracker) lets ParticLS monitor observable
quantities that can be directly compared to observations.

##  5 Example applications

###  5.1 Non-spherical particle geometries

Our first example highlights various particle geometries, which are
implemented by creating children of the Geometry base class. Each new Geometry
must compute its SDF in reference coordinates $X$ and, optionally, can
implement computationally efficient functions to calculate its mass, moment of
inertia, and other properties. Given this information, ParticLS can
incorporate the new geometry into DEM or peridynamic simulations. We define
three example shapes: (i) spheres, (ii) boxes, and (iii) polyhedra—see Figure
3. The SDF for a sphere with radius $R$ is

$\varphi(X)=\|X\|_{2}-R.$ (25)

The SDF for a box is defined by the $d\in\\{2,3\\}$ dimensional vector $Y$
such that $Y_{i}$ is the length of the box in the $i^{\mbox{th}}$ dimension.
Considering the $d\in\\{2,3\\}$ dimensional vector $D$ whose components are
$D_{i}=\max{(|X_{i}|-Y_{i}/2),0)}$; the SDF is

$\varphi(X)=\|D\|_{2}-\min_{i\in\\{1,\ldots,d\\}}{(D_{i})}.$ (26)

Figure 3: Top row: SDFs for the sphere and the box classes—(25) and (26),
respectively. Bottom row: SDFs for two different polygons.

ParticLS can also simulate more complicated geometries. For example, convex
polyhedra, which generalizes the box geometry (see Figure 2), are defined by a
set of faces outlining the shape’s perimeter—lines in two dimensions and
polygons in there dimensions. The signed distance for a polyhedra with $m$
faces is determined by calculating the minimum distance $S_{i}$ between the
point $X$ and the $i^{\mbox{th}}$ face. The SDF for a polygon is therefore

$\varphi(X)=\min_{i\in\\{1,\ldots,m\\}}{(S_{i})}.$ (27)

Figures 4 and 5 show irregularly shaped polygons in two dimensions, and
polyhedra in three dimensions, respectively.

(a) $t=0$

(b) $t=0.5$

(c) $t=1$

Figure 4: We define polygon geometries with density $\delta=1$ that interact
via the viscoelastic contact (section 3.2) with $\kappa_{n}=10^{5}$ $N/m$,
$\kappa_{s}=10^{4}$ $N/m$, $\mu=1$, $\gamma_{n}=1$ $Ns/m$, and
$\gamma_{s}=\sqrt{0.1}$ $Ns/m$. We prescribe an initial velocity of
$v_{0}=(0.25,0.24)$ $m/s$ to the red particle.

(a) $t=0$

(b) $t=0.5$

(c) $t=1$

Figure 5: We define polyhedra geometries by extruding the polygons in Figure 4
in the $z$ direction—the height of each polygon is randomly choosen from
$\\{0.1,0.25,0.5\\}$. We use the same viscoelastic contact law and prescribe
an initial velocity of $v_{0}=(0.25,0.25,0)$ $m/s$ to the red particle.

###  5.2 Time-dependent particle geometries

ParticLS allows users to create objects with time-dependent properties by
implementing children of the StateObject class. In this example, we simulate a
dilating sphere. We implement a class SphereRadius, which is a one dimensional
child of StateObject whose right hand side defines how quickly the sphere’s
radius changes as a function of time. In this example, the radius simply
increases linearly with respect to time. However, users can implement more
complex models that change the particle’s radius according to a function
$dR(t)/dt=g(R(t),t;\theta)$, where $g$ is a generic function, and $\theta$
parameterizes a model that could simulate thermal expansion or a change in
particle shape due to other material properties.

We demonstrate this time-dependent functionality with a simple example in
which normal sphere particles are displaced by one that expands. We begin with
$25$ equally-sized spheres that are configured in a hexagonal packing
configuration with $5$ rows of $5$ spheres (Figure 6). The radius of each
sphere is defined by a SphereRadius object where $dR(t)/dt=0$ $m/s$, except
the middle sphere where $dR(t)/dt=1$ $m/s$. In addition to adding each
RigidBody to the Collection object, we also add the SphereRadius with non-zero
right hand side. The VelocityVerlet time-integrator evolves both the particles
and the radius in time. The contact law between objects is the elastic contact
described in section 3.1 with with parameters $\kappa_{n}=20$ $N/m$,
$\kappa_{s}=\kappa_{n}/2$, and $\mu=0.3$. The results in Figure 6 show how
this time-dependent capability gives users one approach to infuse additional
physics into their ParticLS simulation beyond the rigid body mechanics
provided with DEM.

(a) t=1

(b) t=2

(c) t=3

Figure 6: Example of a particle with a time-dependent geometry. Initially each
of the hex packed particles are the same size but the radius of the middle
particle grows according to $dR(t)/dt=1$ $m/s$.

###  5.3 Multiscale granular media mechanics

ParticLS tracks micro-scale particles (e.g. soil grains), whose aggregate
properties determine the macro-scale response of the granular material (e.g.
soil strength). In general, the DEM is a useful tool for modeling the
mechanical response of granular media and studies have been able to match
modeled results with laboratory experiments Kawamotoetal2016 ; Mcdowell2006 ;
Shmulevich2007 ; Asaf2007 . However, large domains with millions or billions
of particles can be prohibitively expensive. One approach for modeling these
large systems is with multiscale methods, of which there are a few different
options. One method, called the hierarchical multiscale method (HMM),
leverages the finite element method (FEM) over the full domain and uses small-
scale DEM simulations as local representative volume elements (RVEs) to
capture micro-scale dynamics where necessary GuoZhao2014 . In this example, we
present the general principles of the HMM approach and show how the results of
a ParticLS DEM simulation can fit into the framework.

A classic problem for illustrating the benefit of the HMM is simulating the
response of granular media under different loading conditions. In this case we
are concerned with relating the applied loads to the material’s resultant
deformation. In a normal FEM approach to this problem, the governing equation
can be described as:

$\sigma=D\varepsilon$ (28)

where $\sigma$ is the stress, $D$ is the tangent operator, and $\varepsilon$
is the strain GuoZhao2015 . Unfortunately, this requires a constitutive law to
define $D$ beforehand, which introduces assumptions and biases regarding the
material’s response. However, the HMM approach avoids this issue by explicitly
solving for the stress in the granular material by simulating a DEM particle
collection with the given strain boundary conditions. In the HMM, information
is transferred between the two models in an iterative manner such that the FEM
solution provides strain boundary conditions ($\varepsilon$) to each DEM
simulation, and the resolved DEM solution returns stress ($\sigma$) and
tangent operator ($D$) to the FEM GuoZhao2015 . Therefore, we can solve the
equation 28 without making any assumptions about the material’s constitutive
response.

The stress in the DEM particle collection is calculated with the average
Cauchy stress tensor (24), and the tangent operator is calculated with

$D=\frac{1}{V}\sum_{N_{c}}(k_{n}\boldsymbol{n^{c}}\otimes\boldsymbol{d^{c}}\otimes\boldsymbol{n^{c}}\otimes\boldsymbol{d^{c}}+k_{t}\boldsymbol{t^{c}}\otimes\boldsymbol{d^{c}}\otimes\boldsymbol{t^{c}}\otimes\boldsymbol{d^{c}})$
(29)

where $k_{n}$ and $k_{t}$ are the normal and tangential stiffness,
$\boldsymbol{n^{c}}$ and $\boldsymbol{t^{c}}$ are the normal and tangential
directions of each particle contact, and $\boldsymbol{d^{c}}$ is the vector
between each pair of particles in contact. This feedback loop between DEM and
FEM solvers is repeated every time the FEM model requires fine scale
information.

In this example, we show how ParticLS’ DEM capabilities can be used to
simulate an RVE in uniaxial compression—in HMM this strain would come from the
FEM solution. We track the average Cauchy stress tensor of a confined particle
assembly as it is compressed in the $z$ direction. The RVE consists of $16500$
hex-packed spheres that start in a $0.85\times 0.85\times 1.00$ box. The sides
and bottom of the box are planar particles whose location and orientation are
fixed in time. The top of the box is another planar particle whose velocity is
fixed at $v=(0,0,-0.35)$ $m/s$. We run the simulation for $t=1$ units of time
with a viscoelastic contact model with parameters $\kappa_{n}=20$ $N/m$,
$\kappa_{s}=\kappa_{n}/2$, $\gamma_{n}=0.08$ $Ns/m$,
$\gamma_{s}=\gamma_{n}/2$, and $\mu=0.3$; Figure 7 shows the initial and final
configurations. A CollectionTracker was used to track the average Cauchy
stress tensor as a function of time.

Figure 7: Example of uniaxial compression of a sphere particle assembly with
16500 particles and initial dimensions of $0.85\times 0.85\times 1.00$.

This example illustrates ParticLS’ ability to track macro-scale properties
that are important when modeling granular media—specifically, how ParticLS
fits into a multiscale framework. The CollectionTracker measures the
collection’s average Cauchy stress and in this case a ParticleTracker
estimates the strain as the change in top plane’s position. The collection’s
stress-strain response throughout the experiment is illustrated in Figure 8,
which plots the $\sigma_{zz}$ component of the stress against the strain. The
stress-strain response oscillates at the beginning as the top plate meets the
top row of particles and they spread out to fill the box domain. After initial
contact, the stress increases as the collection is compressed. The example RVE
simulation (Figures 7 and 8) shows how ParticLS could be used to model the
RVE’s used in an HMM multi-scale framework for different geomechanic systems
such as soil, snow, or sea ice GuoZhao2014 ; GuoZhao2016 ; Bobillier12018 ;
Gaume2014 ; Hagenmuller2015 .

$0.00$$0.05$$0.10$$0.15$$0.20$$0.25$$0.30$$0.35$$-40.0$$-20.0$$0.0$Strain
$\varepsilon$Stress $\sigma_{zz}$ Figure 8: Stress-strain plot for uniaxial
compression experiment in Figure 7, where the plotted stress value is the
$\sigma_{zz}$ component of (24) and the strain is estimated as the change in
position of the top plane.

###  5.4 Peridynamic model of ordinary material

We demonstrate ParticLS’ ability to model peridynamic solids by stretching a
plate that is initially defined by the material domain $\Omega=1\times 1\times
0.01$ plate. We denote the components of $x\in\Omega$ by
$x=(x_{0},x_{1},x_{2})$. The horizon is determined by a ball with radius
$\delta$,
$\bar{\mathcal{B}}_{\delta}(x)=\mathcal{B}_{\delta}(x)\cap\Omega$—each
material point $x\in\Omega$ interacts with the points $x^{\prime}$ such that
$\|x-x^{\prime}\|<\delta$ and $x^{\prime}\in\Omega$. We consider a linear
peridynamic solid Sillingetal2007  with a hole cut through it, that fractures
given a driving force on its right boundary.

We implement a linear peridynamic solid—refer to sections 13-15 of
Sillingetal2007  for additional details. Given $x\in\Omega$ and
$x^{\prime}\in\bar{\mathcal{B}}_{\delta}(x)$, define the extension

$e(x,x^{\prime},t)=\|(x^{\prime}+u(x^{\prime},t))-(x+u(x,t))\|-\|x^{\prime}-x\|,$
(30)

which is a proxy for the strain in the bond between $x$ and $x^{\prime}$.
Additionally, define the weighted volume

$m(x)=\int_{\bar{\mathcal{B}}_{\delta}(x)}\omega(\|x^{\prime}-x\|)\|x^{\prime}-x\|^{2}\,dx^{\prime}$
(31)

and the dilation parameter

$\theta(x,t)=\int_{\bar{\mathcal{B}}_{\delta}(x)}\omega(\|x^{\prime}-x\|)\|x^{\prime}-x\|e(x,x^{\prime},t)\,dx^{\prime},$
(32)

where $\omega:\mathbb{R}^{+}\mapsto\mathbb{R}^{+}$ is the influence function.
We choose

$\omega(\|\xi\|)=\exp{(-\|\xi\|^{2}/\delta^{2})}.$ (33)

The isotropic and deviatoric parts of $e(x,x^{\prime},t)$ are

$\displaystyle e^{i}(x,x^{\prime},t)$ $\displaystyle=$
$\displaystyle\theta(x,t)\|x^{\prime}-x\|/3$ (34a) $\displaystyle
e^{d}(x,x^{\prime},t)$ $\displaystyle=$ $\displaystyle
e(x,x^{\prime},t)-e^{i}(x,x^{\prime},t),$ (34b)

respectively. A linear peridynamic solid is defined by the force state (13)

$\begin{split}f(x,x^{\prime},t)=&\omega(\|x^{\prime}-x\|)(3k\theta(x,t)\|x^{\prime}-x\|\\\
&-15\mu e^{d}(x,x^{\prime},t))/m(x),\end{split}$ (35)

where $k>0$ and $\mu\geq 0$ are the bulk and shear moduli. The velocities at
$x=0$ and $x=1$ are zero in the $y$ and $z$ directions, but we apply a
velocity of $0.1\text{/sec}$ in the $x$ direction.

We discretize the domain $\Omega$ using $1/N_{x}\times 1/N_{y}\times 1/N_{z}$
cubes whose SDF is defined by (26). ParticLS integrates the discrete
peridynamic equation (10) with respect to time. We cut a hole through the
plate by removing all particles whose initial center of mass $\bar{x}$
satisfies

$(\bar{x}_{1}-0.5)^{2}/r_{1}^{2}+(\bar{x}_{2}-0.5)^{2}/r_{2}^{2}<1,$ (36)

where $\bar{x}_{1}$ and $\bar{x}_{2}$ are the first and second components of
$\bar{x}$ and $r_{1}$ and $r_{2}$ are radii of an ellipse in the $x$ and $y$
directions, respectively. Figure 9 shows the initial location of the particles
with $N_{x}=N_{y}=125$ and $N_{z}=1$. We model internal forces using the
linear peridynamic solid force state defined by (35). However, we drive the
simulation by prescribing a body force

$\beta(x,t)=\begin{cases}3(1-\exp{(-10t)})&\mbox{if }x>0.95\\\
-3(1-\exp{(-10t)})&\mbox{if }x<0.05\\\ 0&\mbox{else}\end{cases}$ (37)

on the portion of the plate such that $\bar{x}_{1}>0.05$. In Figure 9, we show
results for our plate-with-hole experiment. The plate deforms most at the top
and bottom boundaries of the hole, as we would expect, and fractures nucleate
in these regions due to the stress concentrations. Traditional continuum
mechanic models—which define a differential constitutive law—generate a
singularity at crack tip, which makes modeling with these methods changeling,
if not impossible. The peridynamic framework in ParticLS allows us to model
deforming and fracturing bodies. After fracture, the resulting bodies interact
using DEM contact laws, which allows us to model materials that break into
several discrete pieces within one simulation. In Figure 9, we visualize the
material points in the plate using cubes. In the peridynamic setting, this
geometric choice is irrelevant: the peridynamic law only requires the distance
between the center-of-masses. However, after fracture occurs, we replace the
contact law with elastic contact allowing the resulting halves of the plate to
interact as discrete, deforming bodies. Therefore, the chosen particle
geometry only matters _after_ fracture.

(a) $t=0s$

(b) $t=0.25s$

(c) $t=1.0s$

Figure 9: Results of a linear peridynamic plate with a hole that is stretched
to the right under a driving force. The left and right columns of (b) and (c)
are the $x$ and $y$ components of the force density, respectively.

###  5.5 Peridynamic model of a cantilevered beam

To illustrate ParticLS’ ability to match analytical solutions for continuum
mechanic problems, we simulate the deflection of a cantilevered beam in two
and three dimensions. We discretize the beam using material points—similar to
Section 5.4—and connect nearby points (within radius $r$) with a beam that
models internal forces. As the particles move, the beam deforms (compresses,
stretches, and bends), exerting a material force on the connected particles
following Euler-Bernoulli beam theor—Andreetal2012CohesiveBeam  derive the
corresponding expression for the force and torque that this exerts on the
material points.

We propose that this beam model is akin to a peridynamic model. As the number
of particles $n$ discretizing the beam increases, each particle responds to
deformation to the material within a radius $r$. We let $r\rightarrow 0$ as
$n\rightarrow\infty$—SillingLehoucq2008  to show that in this limit the
peridynamic model converges to classical elastic theory. We, therefore,
compare the deflection of the discretized beam with bonded adjacent particles
to the classical analytic solution for a cantilevered beam bending due to an
external acceleration. There are two key differences between this simulation
and how we typically define a peridynamic model. First, we have not explicitly
defined the force state $F[x,x^{\prime},t]$ that is used in Equation (11).
Second, the beams impose a torque on the material point, which requires us to
evolve the material point rotation and angular velocity in order to conserve
angular momentum. Traditional peridynamics—even state-based
peridynamics—assumes angular momentum is conserved by imposing requirements on
the force state. We propose this as a peridynamic-like model but exploring how
to formally include angular momentum in the peridynamic framework is beyond
the scope of the current paper. For our purposes, we numerically validate the
model and our code by matching an analytic solution and showing expected
dynamic behavior.

####  5.5.1 Two dimensional case

The analytical solution for the maximum deflection, $d$, of a cantilevered
beam acting under an external acceleration, $a$, is as follows:

$d=\frac{MaL^{4}}{8EI}.$ (38)

where $M$ is the beam mass, $L$ is the beam length, and $I$ is the beam moment
of inertia. Figure 10 shows the comparison between $n=11$ material points
along the beam and the expected analytical result. It is obvious from this
result that ParticLS is able to accurately approximate the deflection of a
continuous beam.

$0$$100$$200$$300$$400$$500$$600$$700$$800$$900$$1{,}000$$-4$$-3$$-2$$-1$$0$$\cdot
10^{-2}$Young’s Modulus $E$Deflection $d$Analytical solutionParticLS solution
Figure 10: Comparison of numerical beam deflection to analytical solution for
a constant external acceleration of $a=-2E-6$ $m/s^{2}$. We set the physical
parameters $M=1$ $kg$, $L=1$ $m$, and $I=7.85E-9$ $m^{4}$.

####  5.5.2 Three dimensional case

The left column of Figure 11 shows a three dimensional cantilevered beam where
the material points are only bonded with their nearest neighbors. As the
number of points increases, the solution neighborhood containing the connected
points decreases and the numerical model recovers the analytical solution. The
right column fixes a constant radius—as the number of points increases a given
point connects with more nearest neighbors. The latter simulation reflects the
non-local forces used by peridynamic models. Physically, this could represent
internal forces that evolve on near-instantaneous time scales. Mathematically,
this example is more complicated than a typical peridynamic model. In
particular, we do not have an explicit expression for the force state and we
also need to explicitly evolve the angular momentum conservation equation.
However, the peridynamic-like model with non-local forces has the intuitive
effect of stiffening the beam—as shown in Figure 11(c). We leave exploring the
mathematics more rigorously to future work.

(a) $t=0$

(b) Material point bonds

(c) $t=10$

Figure 11: A three dimensional cantilevered beam: (a) the initial condition,
(b) the bonds between material points, and (c) the state given a constant
external acceleration $a=-10^{-1}$ $m/s^{2}$ in the $z$ direction at time
$t=10$ $s$. We set the physical parameters $M=1$ $kg$ and $E=10^{3}$ $Pa$. The
left column shows bonds only between adjacent material points and the right
column shows bonds within a fixed radius $r=0.06$ $m$.

##  6 Conclusions and future directions

ParticLS is a software package for discrete element and peridynamic methods
designed with geoscience applications in mind. Our novel optimization-based
contact detection routine allows users to define complex particle geometries
via an SDF. Leveraging abstraction and inheritance allows users to easily add
new geometries without re-implementing contact-detection algorithms. Also, the
ability to evolve general time-dependent objects allows users to infuse
additional physics into their simulations beyond typical particle mechanics.
ParticLS can also couple rigid body dynamics and peridynamic interactions in
the same simulation, which allows users to model large bodies that fracture
into smaller fragments that then interact with each other in a discrete
manner. These interactions can be fine-tuned for a particular material through
the definitions of new contact models or peridynamic states that capture the
appropriate mechanics.

ParticLS focuses on creating computational objects that mimic mathematical
counterparts (e.g. abstract classes that implement signed distance functions
to represent particle geometries). Without this general framework the software
must either assume a particular geometry (typically spheres) or implement
different pair-wise contact laws for each possible pair of particle
geometries.

In summary, ParticLS provides a platform with an intuitive implementation that
will assist further study of the constitutive laws governing materials prone
to fracture and plastic deformation, specifically rock and ice, as well as
empower the simulation of systems that are relevant to various geoscience
fields. In particular, ParticLS’ ability to simultaneously approximate the
peridynamic equations and evolve time-dependent parameters, such as
temperature, could allow users to model temperature-dependent interaction
laws. For example, materials may become brittle as temperatures decrease or
ductile as temperatures increase. This capability is particularly useful when
modeling the mechanical response of land or sea ice for simulating glacier
calving or sea ice pressure ridging. The coupling of DEM and peridynamics is
particularly useful for modeling ice ridging, where individual ice fragments
coalesce to form sails and keels that catch atmospheric and oceanic currents
Timco1997 ; Hopkins1991 . These examples of future model developments and
applications illustrate how ParticLS provides a useful set of base
functionality that can be extended and used to study a variety of geoscience
problems with complicated physics.

###### Acknowledgements.

The authors gratefully acknowledge funding provided by the Engineering
Research and Development Center (ERDC) Future Innovation Funds (FIF) program.
Since writing this paper, author Andrew D. Davis has moved from the Cold
Regions Research and Engineering Labs to the Courant Institute at New York
University, where he is funded by the Multidisciplinary University Research
Initiative—ONR N00014-19-1-2421.

##### Conflict of Interest

On behalf of all authors, the corresponding author states that there is no
conflict of interest.

## References

  * (1) Andrade, J.E., Lim, K.W., Avila, C.F., Vlahinić, I.: Granular element method for computational particle mechanics.  Computer Methods in Applied Mechanics and Engineering 241, 262–274 (2012)
  * (2) André, D., Iordanoff, I., Charles, J.l., Néauport, J.: Discrete element method to simulate continuous material by using the cohesive beam model.  Computer methods in applied mechanics and engineering 213, 113–125 (2012)
  * (3) Asaf, Z., Rubinstein, D., Shmulevich, I.: Determination of discrete element model parameters required for soil tillage.  Soil and Tillage Research 92(1-2), 227–242 (2007)
  * (4) Askari, E., Bobaru, F., Lehoucq, R.B., Parks, M.L., Silling, S.A., Weckner, O.: Peridynamics for multiscale materials modeling.  Journal of Physics: Conference Series 125, 012078 (2008).  DOI 10.1088/1742-6596/125/1/012078
  * (5) Benning, M., Knoll, F., Schönlieb, C.B., Valkonen, T.: Preconditioned admm with nonlinear operator constraint.  In: IFIP Conference on System Modeling and Optimization, pp. 117–126. Springer (2015)
  * (6) Bentley, J.L.: Multidimensional binary search trees used for associative searching.  Communications of the ACM 18(9), 509–517 (1975)
  * (7) Bessa, M., Foster, J., Belytschko, T., Liu, W.K.: A meshfree unification: reproducing kernel peridynamics.  Computational Mechanics 53(6), 1251–1264 (2014)
  * (8) Blanco, J.L., Rai, P.K.: nanoflann: a C++ header-only fork of FLANN, a library for nearest neighbor (NN) with kd-trees.  https://github.com/jlblancoc/nanoflann (2014)
  * (9) Bobillier1, G., Gaume, J., Van Herwijnen, A., Schweizer, J.: Modeling crack propagation for snow slab avalanche release with discrete elements.  In: Proceedings of the European conference on computational mechanics (2018)
  * (10) Cho, G.C., Dodds, J., Santamarina, J.C.: Particle shape effects on packing density, stiffness, and strength: natural and crushed sands.  Journal of geotechnical and geoenvironmental engineering 132(5), 591–602 (2006)
  * (11) Cundall, P.A., Strack, O.D.: A discrete numerical model for granular assemblies.  Geotechnique 29(1), 47–65 (1979)
  * (12) Favier, J., Abbaspour-Fard, M., Kremmer, M., Raji, A.: Shape representation of axi-symmetrical, non-spherical particles in discrete element simulation using multi-element model particles.  Engineering computations 16(4), 467–480 (1999)
  * (13) Foster, J.T., Silling, S.A., Chen, W.W.: Viscoplasticity using peridynamics.  International journal for numerical methods in engineering 81(10), 1242–1258 (2010)
  * (14) Garcia, X., Latham, J.P., Xiang, J.s., Harrison, J.: A clustered overlapping sphere algorithm to represent real particles in discrete element modelling.  Geotechnique 59(9), 779–784 (2009)
  * (15) Gaume, J., van Herwijnen, A., Schweizer, J., Chambon, G., Birkeland, K.: Discrete element modeling of crack propagation in weak snowpack layers.  In: Proceedings ISSW (2014)
  * (16) Goldstein, T., O’Donoghue, B., Setzer, S., Baraniuk, R.: Fast alternating direction optimization methods.  SIAM Journal on Imaging Sciences 7(3), 1588–1623 (2014)
  * (17) Guo, N., Zhao, J.: A coupled fem/dem approach for hierarchical multiscale modelling of granular media.  International Journal for Numerical Methods in Engineering 99(11), 789–818 (2014)
  * (18) Guo, N., Zhao, J.: A multiscale investigation of strain localization in cohesionless sand.  In: K.T. Chau, J. Zhao (eds.) Bifurcation and Degradation of Geomaterials in the New Millennium, pp. 121–126. Springer International Publishing, Cham (2015)
  * (19) Guo, N., Zhao, J.: Multiscale insights into classical geomechanics problems.  International Journal for Numerical and Analytical Methods in Geomechanics 40(3), 367–390 (2016).  DOI 10.1002/nag.2406.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nag.2406
  * (20) Guo, Y., Curtis, J.S.: Discrete element method simulations for complex granular flows.  Annual Review of Fluid Mechanics 47, 21–46 (2015)
  * (21) Guo, Y., Morgan, J.K.: Influence of normal stress and grain shape on granular friction: Results of discrete element simulations.  Journal of Geophysical Research: Solid Earth 109(B12) (2004)
  * (22) Hagenmuller, P., Chambon, G., Naaim, M.: Microstructure-based modeling of snow mechanics: a discrete element approach.  The Cryosphere 9(5), 1969–1982 (2015).  DOI 10.5194/tc-9-1969-2015.  URL https://www.the-cryosphere.net/9/1969/2015/
  * (23) Hopkins, M., Hibler, W.: On the ridging of a thin sheet of lead ice.  Annals of Glaciology 15, 81–86 (1991)
  * (24) Hopkins, M.A.: A discrete element lagrangian sea ice model.  Engineering Computations 21(2/3/4), 409–421 (2004)
  * (25) Johnson, S., Williams, J.R., Cook, B.: Contact resolution algorithm for an ellipsoid approximation for discrete element modeling.  Engineering Computations 21(2/3/4), 215–234 (2004)
  * (26) Kawamoto, R., Andò, E., Viggiani, G., Andrade, J.E.: Level set discrete element method for three-dimensional computations with triaxial case study.  Journal of the Mechanics and Physics of Solids 91, 1–13 (2016)
  * (27) Kawamoto, R., Andò, E., Viggiani, G., Andrade, J.E.: All you need is shape: Predicting shear banding in sand with ls-dem.  Journal of the Mechanics and Physics of Solids 111, 375–392 (2018)
  * (28) Latham, J., Lu, Y., Munjiza, A.: A random method for simulating loose packs of angular particles using tetrahedra.  Geotechnique 51(10), 871–879 (2001)
  * (29) Latham, J., Munjiza, A.: The modelling of particle systems with real shapes.  Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 362(1822), 1953–1972 (2004)
  * (30) Li, L., Marteau, E., Andrade, J.E.: Capturing the inter-particle force distribution in granular material using ls-dem.  Granular Matter 21(3), 43 (2019)
  * (31) Lim, K.W., Andrade, J.E.: Granular element method for three-dimensional discrete element calculations.  International Journal for Numerical and Analytical Methods in Geomechanics 38(2), 167–188 (2014)
  * (32) Madenci, E., Oterkus, S.: Ordinary state-based peridynamics for plastic deformation according to von mises yield criteria with isotropic hardening.  Journal of the Mechanics and Physics of Solids 86, 192–219 (2016)
  * (33) McDowell, G., Harireche, O., Konietzky, H., Brown, S., Thom, N.: Discrete element modelling of geogrid-reinforced aggregates.  Proceedings of the Institution of Civil Engineers-Geotechnical Engineering 159(1), 35–48 (2006)
  * (34) Monaghan, J.J.: Smoothed particle hydrodynamics.  Annual review of astronomy and astrophysics 30(1), 543–574 (1992)
  * (35) Parikh, N., Boyd, S.: Proximal algorithms.  Foundations and Trends® in Optimization 1(3), 127–239 (2014)
  * (36) Rogula, D.: Introduction to nonlocal theory of material media.  In: Nonlocal theory of material media, pp. 123–222. Springer (1982)
  * (37) Shäfer, J., Dippel, S., Wolf, D.: Force schemes in simulations of granular materials.  Journal de physique I 6(1), 5–20 (1996)
  * (38) Shmulevich, I., Asaf, Z., Rubinstein, D.: Interaction between soil and a wide cutting blade using the discrete element method.  Soil and Tillage Research 97(1), 37–50 (2007)
  * (39) Silling, S.A.: Reformulation of elasticity theory for discontinuities and long-range forces.  Journal of the Mechanics and Physics of Solids 48(1), 175–209 (2000)
  * (40) Silling, S.A., Epton, M., Weckner, O., Xu, J., Askari, E.: Peridynamic states and constitutive modeling.  Journal of Elasticity 88(2), 151–184 (2007)
  * (41) Silling, S.A., Lehoucq, R.: Peridynamic theory of solid mechanics.  In: Advances in applied mechanics, vol. 44, pp. 73–168. Elsevier (2010)
  * (42) Silling, S.A., Lehoucq, R.B.: Convergence of peridynamics to classical elasticity theory.  Journal of Elasticity 93(1), 13 (2008)
  * (43) Song, Y., Turton, R., Kayihan, F.: Contact detection algorithms for dem simulations of tablet-shaped particles.  Powder Technology 161(1), 32–40 (2006)
  * (44) Timco, G., Burden, R.: An analysis of the shapes of sea ice ridges.  Cold Regions Science and Technology 25(1), 65 – 77 (1997).  DOI https://doi.org/10.1016/S0165-232X(96)00017-1.  URL http://www.sciencedirect.com/science/article/pii/S0165232X96000171


Towards Effective and Robust Neural Trojan Defenses via Input Filtering
Kien Do^1The corresponding author, Haripriya Harikumar^1,
Hung Le^1, Dung Nguyen^1,
Truyen Tran^1,

Santu Rana^1, Dang Nguyen^1,
Willy Susilo^2, Svetha Venkatesh^1

^1 Applied Artificial Intelligence Institute (A2I2),
Deakin University, Australia

^2 University of Wollongong, Australia

^1 {k.do, h.harikumar, thai.le,
dung.nguyen, truyen.tran,

santu.rana, d.nguyen<EMAIL_ADDRESS>
^2<EMAIL_ADDRESS>
§ INTRODUCTION

§ STANDARD TROJAN ATTACKS

§ DIFFICULTY IN FINDING INPUT-SPECIFIC TRIGGERS

§ PROPOSED TROJAN DEFENSES

§ EXPERIMENTS

§ RELATED WORK

§ CONCLUSION

§ LIMITATIONS

[1]
Bryant Chen, Wilka Carvalho, Nathalie Baracaldo, Heiko Ludwig, Benjamin
Edwards, Taesung Lee, Ian Molloy, and Biplav Srivastava.
Detecting backdoor attacks on deep neural networks by activation
arXiv preprint arXiv:1811.03728, 2018.

[2]
Huili Chen, Cheng Fu, Jishen Zhao, and Farinaz Koushanfar.
DeepInspect: A Black-box Trojan Detection and Mitigation
Framework for Deep Neural Networks.
In Proceedings of the 28th International Joint Conference on
Artificial Intelligence. AAAI Press, pages 4658–4664, 2019.

[3]
Xinyun Chen, Chang Liu, Bo Li, Kimberly Lu, and Dawn Song.
Targeted backdoor attacks on deep learning systems using data
arXiv preprint arXiv:1712.05526, 2017.

[4]
Bao Gia Doan, Ehsan Abbasnejad, and Damith C Ranasinghe.
Februus: Input purification defense against trojan attacks on deep
neural network systems.
In Annual Computer Security Applications Conference, pages
897–912, 2020.

[5]
Yinpeng Dong, Xiao Yang, Zhijie Deng, Tianyu Pang, Zihao Xiao, Hang Su, and Jun
Black-box detection of backdoor attacks with limited information and
arXiv preprint arXiv:2103.13127, 2021.

[6]
Alhussein Fawzi, Hamza Fawzi, and Omar Fawzi.
Adversarial vulnerability for any classifier.
arXiv preprint arXiv:1802.08686, 2018.

[7]
Yansong Gao, Change Xu, Derui Wang, Shiping Chen, Damith C Ranasinghe, and
Surya Nepal.
Strip: A Defence Against Trojan Attacks on Deep
Neural Networks.
In Proceedings of the 35th Annual Computer Security Applications
Conference, pages 113–125, 2019.

[8]
Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley,
Sherjil Ozair, Aaron Courville, and Yoshua Bengio.
Generative adversarial nets.
Advances in neural information processing systems, 27, 2014.

[9]
Ian J Goodfellow, Jonathon Shlens, and Christian Szegedy.
Explaining and harnessing adversarial examples.
arXiv preprint arXiv:1412.6572, 2014.

[10]
Tianyu Gu, Brendan Dolan-Gavitt, and Siddharth Garg.
Badnets: Identifying vulnerabilities in the machine learning model
supply chain.
arXiv preprint arXiv:1708.06733, 2017.

[11]
Wenbo Guo, Lun Wang, Xinyu Xing, Min Du, and Dawn Song.
Tabor: A highly accurate approach to inspecting and restoring trojan
backdoors in ai systems.
arXiv preprint arXiv:1908.01763, 2019.

[12]
H Harikumar, Vuong Le, Santu Rana, S Bhattacharya, Sunil Gupta, and Svetha
Scalable backdoor detection in neural networks.
In Joint European Conference on Machine Learning and Knowledge
Discovery in Databases (ECML PKDD), pages 289–304. Springer, 2020.

[13]
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.
Deep residual learning for image recognition.
In Proceedings of the IEEE conference on computer vision and
pattern recognition, pages 770–778, 2016.

[14]
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.
Identity mappings in deep residual networks.
In European conference on computer vision, pages 630–645.
Springer, 2016.

[15]
Yujie Ji, Xinyang Zhang, and Ting Wang.
Backdoor attacks against learning systems.
In IEEE Conference on Communications and Network Security,
pages 1–9. IEEE, 2017.

[16]
Diederik P Kingma and Jimmy Ba.
Adam: A Method for Stochastic Optimization.
arXiv preprint arXiv:1412.6980, 2014.

[17]
Diederik P Kingma and Max Welling.
Auto-encoding variational bayes.
arXiv preprint arXiv:1312.6114, 2013.

[18]
Alex Krizhevsky.
Learning Multiple Layers of Features from Tiny Images.
Technical report, 2009.

[19]
Yuezun Li, Yiming Li, Baoyuan Wu, Longkang Li, Ran He, and Siwei Lyu.
Backdoor attack with sample-specific triggers.
arXiv preprint arXiv:2012.03816, 2020.

[20]
Yige Li, Xixiang Lyu, Nodens Koren, Lingjuan Lyu, Bo Li, and Xingjun Ma.
Neural attention distillation: Erasing backdoor triggers from deep
neural networks.
arXiv preprint arXiv:2101.05930, 2021.

[21]
Kang Liu, Brendan Dolan-Gavitt, and Siddharth Garg.
Fine-pruning: Defending against backdooring attacks on deep neural
In International Symposium on Research in Attacks, Intrusions,
and Defenses, pages 273–294. Springer, 2018.

[22]
Yingqi Liu, Wen-Chuan Lee, Guanhong Tao, Shiqing Ma, Yousra Aafer, and Xiangyu
ABS: Scanning Neural Networks for Back-doors by Artificial
Brain Stimulation.
In Proceedings of the ACM SIGSAC Conference on Computer and
Communications Security, pages 1265–1282, 2019.

[23]
Yunfei Liu, Xingjun Ma, James Bailey, and Feng Lu.
Reflection backdoor: A natural backdoor attack on deep neural
In European Conference on Computer Vision, pages 182–199.
Springer, 2020.

[24]
Yuntao Liu, Yang Xie, and Ankur Srivastava.
Neural trojans.
In 2017 IEEE International Conference on Computer Design
(ICCD), pages 45–48. IEEE, 2017.

[25]
Ziwei Liu, Ping Luo, Xiaogang Wang, and Xiaoou Tang.
Deep learning face attributes in the wild.
In Proceedings of the IEEE international conference on computer
vision, pages 3730–3738, 2015.

[26]
Christos Louizos, Karen Ullrich, and Max Welling.
Bayesian compression for deep learning.
arXiv preprint arXiv:1705.08665, 2017.

[27]
Aleksander Madry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and
Adrian Vladu.
Towards deep learning models resistant to adversarial attacks.
arXiv preprint arXiv:1706.06083, 2017.

[28]
Pim Moeskops, Mitko Veta, Maxime W Lafarge, Koen AJ Eppenhof, and Josien PW
Adversarial training and dilated convolutions for brain mri
In Deep learning in medical image analysis and multimodal
learning for clinical decision support, pages 56–64. Springer, 2017.

[29]
Dmitry Molchanov, Arsenii Ashukha, and Dmitry Vetrov.
Variational dropout sparsifies deep neural networks.
In International Conference on Machine Learning, pages
2498–2507. PMLR, 2017.

[30]
Luis Muñoz-González, Bjarne Pfitzner, Matteo Russo, Javier
Carnerero-Cano, and Emil C Lupu.
Poisoning attacks with generative adversarial nets.
arXiv preprint arXiv:1906.07773, 2019.

[31]
Anh Nguyen and Anh Tran.
Input-aware dynamic backdoor attack.
arXiv preprint arXiv:2010.08138, 2020.

[32]
Anh Nguyen and Anh Tran.
Wanet–imperceptible warping-based backdoor attack.
International Conference on Learning Representations, 2021.

[33]
Omkar M Parkhi, Andrea Vedaldi, and Andrew Zisserman.
Deep face recognition.

[34]
Ximing Qiao, Yukun Yang, and Hai Li.
Defending neural backdoors via generative distribution modeling.
arXiv preprint arXiv:1910.04749, 2019.

[35]
Joseph Redmon, Santosh Divvala, Ross Girshick, and Ali Farhadi.
You only look once: Unified, real-time object detection.
In Proceedings of the IEEE conference on computer vision and
pattern recognition, pages 779–788, 2016.

[36]
Aniruddha Saha, Akshayvarun Subramanya, and Hamed Pirsiavash.
Hidden trigger backdoor attacks.
In Proceedings of the AAAI Conference on Artificial
Intelligence, volume 34, pages 11957–11965, 2020.

[37]
Ahmed Salem, Rui Wen, Michael Backes, Shiqing Ma, and Yang Zhang.
Dynamic backdoor attacks against machine learning models.
arXiv preprint arXiv:2003.03675, 2020.

[38]
Florian Schroff, Dmitry Kalenichenko, and James Philbin.
Facenet: A unified embedding for face recognition and clustering.
In Proceedings of the IEEE conference on computer vision and
pattern recognition, pages 815–823, 2015.

[39]
Ramprasaath R Selvaraju, Michael Cogswell, Abhishek Das, Ramakrishna Vedantam,
Devi Parikh, and Dhruv Batra.
Grad-cam: Visual explanations from deep networks via gradient-based
In Proceedings of the IEEE international conference on computer
vision, pages 618–626, 2017.

[40]
Ali Shafahi, W Ronny Huang, Mahyar Najibi, Octavian Suciu, Christoph Studer,
Tudor Dumitras, and Tom Goldstein.
Poison frogs! targeted clean-label poisoning attacks on neural
arXiv preprint arXiv:1804.00792, 2018.

[41]
Johannes Stallkamp, Marc Schlipsing, Jan Salmen, and Christian Igel.
Man vs. Computer: Benchmarking Machine Learning Algorithms
for Traffic Sign Recognition.
Neural Networks, pages 323–332, 2012.

[42]
Waqas Sultani, Chen Chen, and Mubarak Shah.
Real-world anomaly detection in surveillance videos.
In Proceedings of the IEEE conference on computer vision and
pattern recognition, pages 6479–6488, 2018.

[43]
Simen Thys, Wiebe Van Ranst, and Toon Goedemé.
Fooling automated surveillance cameras: adversarial patches to attack
person detection.
In Proceedings of the IEEE/CVF Conference on Computer Vision and
Pattern Recognition Workshops, pages 0–0, 2019.

[44]
Brandon Tran, Jerry Li, and Aleksander Madry.
Spectral Signatures in Backdoor Attacks.
In Advances in Neural Information Processing Systems, pages
8000–8010, 2018.

[45]
Mart van Baalen, Christos Louizos, Markus Nagel, Rana Ali Amjad, Ying Wang,
Tijmen Blankevoort, and Max Welling.
Bayesian bits: Unifying quantization and pruning.
arXiv preprint arXiv:2005.07093, 2020.

[46]
Pascal Vincent, Hugo Larochelle, Yoshua Bengio, and Pierre-Antoine Manzagol.
Extracting and composing robust features with denoising autoencoders.
In Proceedings of the 25th international conference on Machine
learning, pages 1096–1103, 2008.

[47]
Bolun Wang, Yuanshun Yao, Shawn Shan, Huiying Li, Bimal Viswanath, Haitao
Zheng, and Ben Y Zhao.
Neural Cleanse: Identifying and Mitigating Backdoor Attacks
in Neural Networks.
In IEEE Symposium on Security and Privacy, pages 707–723.
IEEE, 2019.

[48]
Ren Wang, Gaoyuan Zhang, Sijia Liu, Pin-Yu Chen, Jinjun Xiong, and Meng Wang.
Practical detection of trojan neural networks: Data-limited and
data-free cases.
In Computer Vision–ECCV 2020: 16th European Conference,
Glasgow, UK, August 23–28, 2020, Proceedings, Part XXIII 16, pages
222–238. Springer, 2020.

[49]
Emily Wenger, Josephine Passananti, Arjun Nitin Bhagoji, Yuanshun Yao, Haitao
Zheng, and Ben Y Zhao.
Backdoor attacks against deep learning systems in the physical world.
In Proceedings of the IEEE/CVF Conference on Computer Vision and
Pattern Recognition, pages 6206–6215, 2021.

[50]
Dong Yang, Daguang Xu, S Kevin Zhou, Bogdan Georgescu, Mingqing Chen, Sasa
Grbic, Dimitris Metaxas, and Dorin Comaniciu.
Automatic liver segmentation using an adversarial image-to-image
In International Conference on Medical Image Computing and
Computer-Assisted Intervention, pages 507–515. Springer, 2017.

[51]
Pu Zhao, Pin-Yu Chen, Payel Das, Karthikeyan Natesan Ramamurthy, and Xue Lin.
Bridging mode connectivity in loss landscapes and adversarial
arXiv preprint arXiv:2005.00060, 2020.

[52]
Chen Zhu, W Ronny Huang, Hengduo Li, Gavin Taylor, Christoph Studer, and Tom
Transferable clean-label poisoning attacks on deep neural nets.
In International Conference on Machine Learning, pages
7614–7623. PMLR, 2019.

§ APPENDIX


# Exploring The Relationship Between Road Infrastructure and Crimes in
Memphis, Tennessee

Alexandre Signorel

(December 1, 2022)

###### Abstract

Memphis, Tennessee is one of the cities with highest crime rate in the United
States. In this work, we explore the relationship between road infrastructure,
especially potholes, and crimes. The pothole and crime data are collected from
[1] between 2020 and 2022. The crime data report various crimes in the Memphis
area, which contain the location, time, and type of the crime. The pothole
data is part of the Open 311 data [1], which contains information of different
infrastructure projects, including the location of the project, and the
starting and ending dates of the project. We focus on infrastructure projects
regarding pothole repairs.

##  1 Introduction

Memphis, Tennessee is one of the cities with highest crime rate in the United
States. In this task, we are interested in studying the potential relationship
between the road infrastructure, potholes in particular, and crimes in
Memphis. The data is collected from the Memphis Data Hub [1], There are three
main data categories that are currently accessible: _Public Safety_ which
contains crime data; _Open 311_ which contains data about current and past
infrastructure projects; and _Civic Assets_ which contains data about public
services.

For our project, we focus mainly on Public Safety and Open 311 where in Open
311 we specifically focus on data regarding potholes. Our goal is to explore
potential correlations between crimes that happened in an area and the road
condition of that area, evaluated by the amount of potholes reported.

##  2 Working with Data

The data are downloaded in the CSV format. We process the data and extract the
data from 2020 to 2022. This is the time period when both pothole data and
crime data are available.

###  2.1 Crime Data

To process the crime data, we convert crime names into numerical values that
range from 100 to 1700. This indicates the overall 17 types of crimes. The
last two digits of the numerical values indicate the sub-categories of crimes.
Next, we remove information that does not concern the crime’s location and
time. Lastly, we convert the date to the Date-Time format, which displays the
year, month, and day of a crime. An example of the processed crime data is
shown in Figure 1.

Figure 1: An example of the processed crime data.

###  2.2 Pothole Data

We extract the pothole data from the Open 311 data. We convert the extracted
date to the Date-Time format as of the pothole data.

Figure 2: An example of the processed pothole data.

###  2.3 Visualizing Data

In order to qualitatively study the data, we visualize both types of data by
dividing them into windows of three months. We then use Folium [2], a python
package that uses OpenStreetMap [3], to generate a view of the Memphis area.
Next, we use Folium.plugins to generate a heat map which is overlaid onto the
map of Memphis. The heat map of crimes is shown in Figure 4 and the heat map
of potholes is shown in Figure 4.

Figure 3: The heat map of crimes in Memphis between April 2020 and June 2020.

Figure 4: The heat map of potholes reported in Memphis between April 2020 and
June 2020.

###  2.4 Analyzing Data

In order to quantitatively study the data, we convert both types of data into
windows of a week between January 2020 and October 2022. The scatter plots of
crimes and potholes are shown in Figure 6 and Figure 6, respectively.

Figure 5: The scatter plot of crimes in Memphis between January 2020 and
October 2022.

Figure 6: The scatter plot of potholes between January 2020 and October 2022.

##  3 Results

After examining the data, we find that both types of data exhibit wave
patterns. The crime data first go up and then back down with an upward trend
within the last wave. The pothole data exhibit the wave pattern more
explicitly with the lowest being around the week 40 of each year. In all three
years, we can see that the number of crimes reaches a peak while the amount of
potholes reaches a bottom. However, with the limited data and fluctuating
patterns, it is difficult to claim any correlation between the two.

##  4 Conclusion and Future Work

In this task, we explore the relationship between potholes and crimes in
Memphis, Tennessee. We find that both types of data exhibit wave patterns.
However, it is difficult to claim that they are strictly correlated.

In the future, we would like to expand our study to cover more infrastructure
such as buildings and public services within an area. In addition, instead of
conducting a correlation study of the entire Memphis area, we can compare
higher-crime-rate areas with lower-crime-rate areas to explore the
relationships between various infrastructure factors and the crimes.

## References

  * [1] Memphis  “Memphis Data Hub”, https://data.memphistn.gov/, 2022
  * [2] Folium  “Folium visualization package”, https://python-visualization.github.io/folium/, 2022
  * [3] OSM  “OpenStreetMap”, https://www.openstreetmap.org/, 2022