# Potential Particles and Potential Blocks¶

The origins of scientific development regarding the algorithms described in this section are traced back to: [Boon2012] (PotentialBlocks code), [Boon2013b] (PotentialParticles code) and [Boon2015] (Block Generation code).

## Introduction¶

This section discusses two codes to simulate (i) non-spherical particles using the concept of the Potential Particles [Houlsby2009], with the solution procedures in [Boon2013] for 3-D and (ii) polyhedral blocks using planar linear inequalities, based on linear programming concepts [Boon2012]. These codes define two shape classes in YADE, named PotentialParticle and PotentialBlock, respectively. Besides some similarities in syntax, they have distinct differences, concerning morphological characteristics of the particles and the methods used to facilitate contact detection.

The Potential Particles code (abbreviated herein as PP) is detailed in [Boon2013], where non-spherical particles are assembled as a combination of 2nd degree polynomial functions and a fraction of a sphere, while their edges are rounded with a user-defined radius of curvature.

The Potential Blocks code (abbreviated herein as PB) is used to simulate polyhedral particles with flat surfaces, based on the work of [Boon2012], where a smooth, inner potential particle is used to calculate the contact normal vector. This code is compatible with the Block Generation algorithm defined in [Boon2015], in which potential blocks can be generated by intersections of original, intact blocks with discontinuity planes.

## Enabling the PP & PB codes during installation¶

The PP and PB codes are not currently included in the binary, pre-packaged versions of yade or yadedaily (i.e. they cannot be accessed through the apt-get environment). They can be easily used though, if YADE is built and installed from source, following the installation instructions and the additional steps below. These two codes are independent, in the sense that either one of them can be compiled/used separately, without enabling the other, while they do not interact with each other (i.e. we cannot establish contact between a PB and a PP). Enabling the PB code causes an automatic compilation of the Block Generation algorithm.

To use the PB code, the CLP library must be installed before compilation, using the apt-get environment:

sudo apt-get install coinor-clp coinor-libclp-dev


Then, the PP and PB codes are enabled during configuration, using:

cmake -DENABLE_POTENTIAL_PARTICLES=ON -DENABLE_POTENTIAL_BLOCKS=ON -DCMAKE_INSTALL_PREFIX=../install ../trunk


Alternatively, these codes can be enabled before configuration manually, by changing the options ENABLE_POTENTIAL_PARTICLES and ENABLE_POTENTIAL_BLOCKS from OFF to ON inside the CMakeLists.txt file.

## Potential Particles code (PP)¶

The concept of Potential Particles was introduced and developed by [Houlsby2009]. The problem of contact detection between a pair of potential particles was cast as a constrained optimization problem, where the equations are solved using the Newton-Raphson method in 2-D. In [Boon2013] it was extended to 3-D and more robust solutions were proposed. Many numerical optimization solvers generally cannot cope with discontinuities, ill-conditioned gradients (Jacobians) or curvatures (Hessians), and these obstacles were overcome in [Boon2013], by re-formulating the problem and solving the equations using conic optimization solvers. Users can either make use of MOSEK (using its academic licence) or the simplified code written by [Boon2013], which is used by default in the source code. A potential particle is defined as in (1) [Houlsby2009]:

(1)$\begin{split}f=(1-k) \Bigg(\sum_{i=1}^{N} \langle a_{i}x+b_{i}y+c_{i}z-d_i\rangle ^2-r^2 \Bigg) +k (x^2+y^2+z^2-R^2)\\\end{split}$

where $$(a_i, b_i, c_i)$$ is the normal vector of the $$i^{th}$$ plane, defined with respect to the particle’s local coordinate system and $$d_i$$ is the distance of the plane to the local origin. $$\langle \;\rangle$$ are Macaulay brackets, i.e., $$〈x〉 = x$$ for $$x > 0$$; $$\langle x \rangle = 0$$ for $$x \leq 0$$. The planes are assembled such that their normal vectors point outwards. They are summed quadratically and expanded by a distance r, which is also related to the radius of the curvature at the corners. Furthermore, a “shadow” spherical particle is added; R is the radius of the sphere, with $$0 < k \; \leq \; 1$$, denoting the fraction of sphericity of the particle. The geometry of some cuboidal potential particles is displayed in Fig. fig-pp, for different values of the parameter $$k$$.

Construction of potential particles (a) constituent planes are squared and expanded by a constant r. A fraction of sphere is added. Particles with the spherical term are visible in (b) k=0.9, (c) k=0.7, and (d) k=0.4 (after [Boon2013]).

The potential function is normalized for computational reasons in the form (2) [Houlsby2009]:

(2)$\begin{split}f=(1-k) \Bigg(\sum_{i=1}^{N} \frac{ \langle a_{i}x+b_{i}y+c_{i}z-d_i\rangle^2 }{ r^2 } -1 \Bigg) +k \Bigg( \frac{ x^2+y^2+z^2 }{ R^2 }-1 \Bigg)\\\end{split}$

This potential function takes values:

• $$f=0$$: on the particle surface
• $$f<0$$: inside the particle
• $$f>0$$: outside the particle

To ensure numerical stability, it is not advised to use the extreme values k=0 and k=1 or values approaching these extremes. In particular, k=0 cannot be used from a theoretical standpoint, since the Potential Particles were formulated for strictly convex shapes (curved faces).

## Potential Blocks code (PB)¶

The Potential Blocks code was developed during the D.Phil. thesis of CW Boon [Boon2013b] and discussed in [Boon2012]. It was developed originally for rock engineering applications, to model polygonal and polyhedral blocks with flat surfaces. The blocks are defined with linear inequalities only and unlike the PotentialParticle shape class, no spherical term is considered (so, practically k=0). Although k and R are input parameters of the PotentialBlock shape class, their existence during computation is null. In particular, R is used within the source code, denoting a characteristic dimension of the blocks, but does not reflect the radius of a “shadow particle”, like it does for the Potential Particles. This value of R is used in the Potential Blocks code to calculate the initial bi-section step size for line search, to obtain a point on the particle, which in turn is used to calculate the overlap distance during contact.

For a convex particle defined by N planes, the space that it occupies can be defined using the following inequalities (3):

(3)$\begin{split}a_{i}x + b_{i}y + c_{i}z \; \leq \; d_{i}, i=1:N\\\end{split}$

where $$(a_i, b_i, c_i)$$ is the unit normal vector of the $$i^{th}$$ plane, defined with respect to the particle’s local coordinate system, and $$d_i$$ is the distance of the plane to the local origin. According to [Boon2012], an inner, smooth potential particle is used to calculate the contact normal, formulated as in (4):

(4)$\begin{split}f=\sum_{i=1}^{N} \langle a_{i}x + b_{i}y + c_{i}z - d_i + r\rangle^2\\\end{split}$

This potential particle is defined inner by a distance r inside the actual particle, with edges rounded by a radius or curvature r, as well (see Fig. fig-pbInner).

A potential particle is defined inside the actual particle. The normal vector of the particle at any point can be calculated from the first derivative of the potential particle. (after [Boon2012]).

In YADE, the Potential Blocks have a slightly different mathematical expression, since their shape is generated as an assembly of planes as in (5):

(5)$\begin{split}a_{i}x + b_{i}y + c_{i}z - d_{i} - r = 0, i=1:N\\\end{split}$

while the inner Potential Particle used to calculate the contact normal is defined as in (6):

(6)$\begin{split}f=\sum_{i=1}^{N} \langle a_{i}x + b_{i}y + c_{i}z - d_i\rangle^2.\\\end{split}$

Now, the Potential Block surface is at a distance of $$(d_{i}+r)$$ from the local particle center, while the inner potential particle is at a distance $$d$$ from the local particle center.

It is worth to emphasize on the fact that the shape of a Potential Block is defined using an assembly of planes and not a single, implicit potential function, like we have for the Potential Particles code. The inner potential particle in the Potential Blocks code is only used to calculate the contact normal.

The problem of establishing intersection between a pair of blocks is cast as a standard linear programming problem of finding a feasible region which satisfies all the linear inequalities defining both blocks. The contact point is calculated as the analytic centre of the feasible region, a well-known concept of interior-point methods in convex optimization calculations. The contact normal is obtained from the gradient of a smooth “potential particle” defined inside the block. The overlap distance is calculated through bi-section searching along the contact normal, within the overlap region.

A potential block. The normal vectors of the faces point outwards (after [Boon2013b]).

The linear programming solver for Potential Blocks was originally CPLEX, but has been updated to CLP, developed by COIN-OR, since the latter can be downloaded from Ubuntu or Debian’s distributions without requiring an academic licence.

## Engines¶

The PP and PB codes use their own classes to handle bounding volumes, contact geometry & physics and recording of outputs in vtk format, while they derive the interparticle friction angle from the frictional material class FrictMat. The syntax used to define these classes is similar, unless if specified otherwise.

Shape PotentialParticle PotentialBlock
Material FrictMat FrictMat
Bound functor PotentialParticle2AABB PotentialBlock2AABB
IGeom functor Ig2_PP_PP_ScGeom Ig2_PB_PB_ScGeom
IPhys functor Ip2_FrictMat_FrictMat_KnKsPhys Ip2_FrictMat_FrictMat_KnKsPBPhys
Law2 functor Law2_SCG_KnKsPhys_KnKsLaw Law2_SCG_KnKsPBPhys_KnKsPBLaw
VTK Recorder PotentialParticleVTKRecorder PotentialBlockVTKRecorder

A simple simulation loop using the Potential Blocks reads as:

O.engines=[
ForceResetter(),
InsertionSortCollider([PotentialBlock2AABB()], verletDist=0.01),
InteractionLoop(
[Ig2_PB_PB_ScGeom()],
[Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=1e8, ks_i=1e8, Knormal=1e8, Kshear=1e8, twoDimension=True, unitWidth2D=1.0, calJointLength=True, viscousDamping=0.2)],
[Law2_SCG_KnKsPBPhys_KnKsPBLaw(label='law', neverErase=False)]
),
NewtonIntegrator(damping=0.2, exactAsphericalRot=True, gravity=[0,0,-9.81]),
PotentialBlockVTKRecorder(fileName='./vtk/file_prefix', iterPeriod=1000, twoDimension=True, sampleX=30, sampleY=30, sampleZ=30, maxDimension=0.2, label='vtkRecorder')
]


Attention should be given to the twoDimension parameter, which defines whether a contact should be handled as 2-D or 3-D. If twoDimension=True, then the contactArea parameter is calculated as:

if(phys->twoDimension) { phys->contactArea = phys->unitWidth2D*phys->jointLength;}


The unitWidth2D is given by the user (usually equal to 1.0), while the jointLength parameter is controlled by the boolean parameter calJointLength. If calJointLength=True, the jointLength parameter is calculated for that specific contact. If caljointLength=False, a default jointLength (=1.0) is used for all the contacts that occur in the system.

The calculation of an automatically defined, contact-specific jointLength is currently available only for the Potential Blocks code. Also, for the Potential Blocks code only, if twoDimension=False, the contactArea parameter is calculated for the 3-D case as well, as described in [Boon2013b] for a pair of 3-D blocks.

For the Potential Particles code, the boolean parameter calJointLength should be kept at a value calJointLength=False and the parameter jointLength should be kept equal to 1.0, until a calculation of the contactArea is developed.

## Geometric definition of a PP and a PB¶

A strong merit of the Potential Particles and the Potential Blocks codes lies in the fact that the geometric definition of the particle shape and the contact detection problem are resolved using only the equations of the faces of the particles. In this way, using a single data structure, there is no need to store information about the vertices or their connectivity to establish contact, a feature that makes them computationally affordable, while all contacts are handled in the same way (there is no need to distinguish among face-face, face-edge, face-vertex, edge-edge, edge-vertex or vertex-vertex contacts). Due to this, the geometry of a particle is defined in the shape class using the values of the normal vectors of the faces and the distances of the faces from the local origin.

For example, to define a cuboid (6 faces) with rounded edges, an edge length of D, centred to its local centroid and aligned to its principal axes, using the Potential Particles code, we type:

r=D/100.
k=0.3
R=D/2.
b=Body()
b.shape=PotentialParticle( r=r, k=k, R=R,
a=[   1.0,    -1.0,     0.0,     0.0,     0.0,     0.0],
b=[   0.0,     0.0,     1.0,    -1.0,     0.0,     0.0],
c=[   0.0,     0.0,     0.0,     0.0,     1.0,    -1.0],
d=[D/2.-r,  D/2.-r,  D/2.-r,  D/2.-r,  D/2.-r,  D/2.-r], ...)


The first element of the vector parameters $$a, b, c, d$$ refers to the normal vector of the first plane, the second element to the second plane, and so on.

Using the Potential Particles code, this is not a perfect cube, since the particle geometry is defined by a potential function as in (2). It is reminded that within this potential function, these planes are summed quadratically, the particle edges are rounded by a radius of curvature r and then the particle faces are curved by the addition of a “shadow” spherical particle with a radius R, to a percentage defined by the parameter k. A value r is deducted from each element of the vector parameter d, to compensate for expanding the potential particle by r.

The parameters $$a_{i}, b_{i}, c_{i}, d_{i}$$ stated above correspond to the planes used in (5):

$\begin{split}1.0 x + 0.0 y + 0.0 z = D/2 \Leftrightarrow +x=D/2\\ -1.0 x + 0.0 y + 0.0 z = D/2 \Leftrightarrow -x=D/2\\ 0.0 x + 1.0 y + 0.0 z = D/2 \Leftrightarrow +y=D/2\\ 0.0 x - 1.0 y + 0.0 z = D/2 \Leftrightarrow -y=D/2\\ 0.0 x + 0.0 y + 1.0 z = D/2 \Leftrightarrow +z=D/2\\ 0.0 x + 0.0 y - 1.0 z = D/2 \Leftrightarrow -z=D/2\\\end{split}$

To model a cube with an edge of D, using the Potential Blocks code, we type:

r=D/100.
R=D/2.
b=Body()
b.shape=PotentialBlock( r=r, R=R,
a=[   1.0,    -1.0,     0.0,     0.0,     0.0,     0.0],
b=[   0.0,     0.0,     1.0,    -1.0,     0.0,     0.0],
c=[   0.0,     0.0,     0.0,     0.0,     1.0,    -1.0],
d=[D/2.-r,  D/2.-r,  D/2.-r,  D/2.-r,  D/2.-r,  D/2.-r], ...)


Using the Potential Blocks code, this particle will have sharp edges and flat faces in what regards its geometry (i.e. the space it occupies), defined by the given planes, while for the calculation of the contact normal, an inner potential particle with rounded edges is used, formulated as in (6), located fully inside the actual particle. The distances of the planes from the local origin, stored in the vector parameter d, are reduced by r to achieve an exact edge length of D, using the (5).

To ensure numerical stability, it is advised to normalize the normal vector of each plane, so that $${a_{i}}^2 + {b_{i}}^2 + {c_{i}}^2 = 1$$. There is no limit to the number of the particle faces that can be used, a feature that allows the modelling of a variety of convex particle shapes.

In practice, it is usual for the geometry of a particle to be given in terms of vertices & their connectivity (e.g. in the form of a surface mesh, like in .stl files). In such cases, the user can calculate the normal vector of each face, which will give the coefficients $$a_{i}, b_{i}, c_{i}$$ and using a vertex of each face, then calculate the coefficients $$d_{i}$$. A python routine to perform this without any additional effort by the user is currently being developed.

## Body definition of a PP or a PB¶

To define a body using the PotentialParticle or PotentialBlock shape classes, it has to be assembled using the _commonBodySetup function, which can be found in the file py/utils.py. For example, to define a PotentialBlock:

O.materials.append(FrictMat(young=150e6,poisson=.30,frictionAngle=radians(0.0),density=2650,label='frictionless'))

b=Body()
b.shape=PotentialBlock(...)
b.aspherical=True # To be used in conjunction with exactAsphericalRot=True in the NewtonIntegrator
# V: Volume
# I11, I22, I33: Principal inertias
utils._commonBodySetup(b,V,Vector3(I11,I22,I33), material='frictionless', pos=(0,0,0), ori=Quaternion((1,0,0),0), fixed=False)
b.state.pos=Vector3(xPos,yPos,zPos)
b.state.ori=Quaternion((random.random(),random.random(),random.random()),random.random())
b.shape.volume=V;
O.bodies.append(b)


The particle must be initially defined, so that the local axes coincide with its principal axes, for which the inertia tensor is diagonal.

It should be noted that the principal inertia values used here should be divided with the density of the considered material, since the I11, I22, I33 values given by the user are multiplied with the density inside the _commonBodySetup function. The mass of the particle is calculated within the same function as well, so we do not need to set manually b.mass=V*density.

## Boundary Particles¶

The PP & PB codes support the definition of boundary particles, which interact only with non-boundary ones. These particles can have a variety of uses, e.g. to model loading plates acting on a granular sample, while different uses can emerge for different applications. A particle can be set as a boundary one in both codes, using the boolean parameter isBoundary inside the shape class.

In the PP code, all particles interact with the same normal and shear contact stiffness Knormal and Kshear, defined in the Ip2_FrictMat_FrictMat_KnKsPhys functor.

The PB code supports the definition of different contact stiffness values for interactions between boundary and non-boundary or non-boundary and non-boundary particles. When isBoundary=False, the PotentialBlock in question is handled to interact with normal and shear stiffnesses of Knormal and Kshear, respectively, with other non-boundary particles. When isBoundary=True, the PotentialBlock in question is handled to interact with normal and shear stiffnesses of kn_i and ks_i, respectively, with non-boundary particles.

## Visualization¶

Visualization of the PotentialParticle and PotentialBlock shape classes is offered using the qt environment (OpenGL). Additionally, the PotentialParticleVTKRecorder and PotentialBlockVTKRecorder classes can be used to export geometrical and interaction information of the analyses in vtk format (visualized in Paraview). It should be noted that currently the PotentialBlockVTKRecorder records the inner, rounded potential particle, rather than the actual particle with sharp edges and flat faces.

In the qt environment, the PotentialParticle shape class is visualized using the Marching Cubes algorithm, and the level of display accuracy can be determined by the user. This is controlled by the parameters:

# Potential Particles
Gl1_PotentialParticle.sizeX=20
Gl1_PotentialParticle.sizeY=20
Gl1_PotentialParticle.sizeZ=20


A similar choice exists for output in vtk format, using the PotentialParticleVTKRecorder or PotentialBlockVTKRecorder, syntaxed as:

# Potential Particles
PotentialParticleVTKRecorder(sampleX=30, sampleY=30, sampleZ=30, maxDimension=20)

# Potential Blocks
PotentialBlockVTKRecorder(sampleX=30, sampleY=30, sampleZ=30, maxDimension=20)


The parameters sizeX,Y,Z (for OpenGL visualization) and sampleX,Y,Z (for output in vtk format) represent the number of subdivisions of the Aabb of the particle to a grid, which will be used to draw its geometry, in respect to the global axes X, Y, Z. Larger values will result to a more accurate display of the particles’ shape, but will slow down the visualization speed in qt and writing speed of the .vtk files and increase the size of the .vtk files. For output in vtk format, users can also define the parameter maxDimension, which overrides the selected sampleX,Y,Z values if they are too small, as described below:

$\begin{split}if \; \mid xmax-xmin \mid /sampleX > maxDimension \Rightarrow sampleX = \mid xmax-xmin \mid /maxDimension\\ if \; \mid ymax-ymin \mid /sampleY > maxDimension \Rightarrow sampleY = \mid ymax-ymin \mid /maxDimension\\ if \; \mid zmax-zmin \mid /sampleZ > maxDimension \Rightarrow sampleZ = \mid zmax-zmin \mid /maxDimension \; \\\end{split}$

The PotentialParticleVTKRecorder and PotentialBlockVTKRecorder also support optionally the recording of the particles’ velocities (linear and angular), interaction information (contact point and forces), colors and ids, using:

# Potential Particles
PotentialParticleVTKRecorder(..., REC_VELOCITY=True, REC_INTERACTION=True, REC_COLORS=True, REC_ID=True)

# Potential Blocks
PotentialBlockVTKRecorder(..., REC_VELOCITY=True, REC_INTERACTION=True, REC_COLORS=True, REC_ID=True)


Force chains and other visual outputs are available in qt by default, while they can be extracted in vtk format using the classic VTKRecorder or the export.VTKExporter class.

A boolean parameter twoDimension exists to specify whether the particles will be rendered as 2-D or 3-D in the vtk output:

# Potential Particles
PotentialParticleVTKRecorder(..., twoDimension=False)

# Potential Blocks
PotentialBlockVTKRecorder(..., twoDimension=False)


This parameter should not be mixed up with the Ip2_FrictMat_FrictMat_KnKsPBPhys.twoDimension parameter, which is used to define how the contact forces are calculated, as described in the Engines section.

## Axis-Aligned Bounding Box¶

The PP & PB codes use their own BoundFunctors, called PotentialParticle2AABB and PotentialBlock2AABB, respectively, to define the Axis-Aligned Bounding Box of each particle. In both bound functors, a boolean parameter AabbMinMax exists, allowing the user to choose between an approximate cubic Aabb or a more accurate one.

In particular, if AabbMinMax=False, a cubic Aabb is considered with dimensions 1.0*R. This is implemented for both the PP and PB codes, even though the Potential Blocks do not have a spherical term. In this case, the radius R is used as a reference length, denoting half the diagonal of the cubic Aabb. Usage of this approximate cubic Aabb is not advised, since it can increase the number of empty contacts, increasing thus the time needed to facilitate the approximate contact detection, while it relies on the radius R, the value of which should enclose the whole particle if this option is activated.

If AabbMinMax=True, a more accurate Aabb can be defined. Currently, the initial Aabb of a PotentialParticle has to be defined manually by the user, in the particle local coordinate system and for the initial orientation of the particle. To do so, the user has to manually specify the two extreme points of the Aabb: minAabbRotated, maxAabbRotated inside the shape class. The Aabb for a PotentialBlock, on the other hand, is calculated and updated automatically from the vertices of the particle, if the boolean parameter AabbMinMax =True.

As discussed in the subsection Visualization, the dimensions of the Aabb are used as a drawing space in the code implementing rendering of the particles in the qt environment (for the PP code) and for the creation of the output files in vtk format (for both codes). This is achieved by using two auxiliary parameters: minAabb and maxAabb. For the particles to be properly rendered as closed surfaces in both qt and vtk outputs using the available codes, we need to define a drawing space slightly larger than the actual one. Here, this drawing space is represented by the Aabb of the particles, and thus the differentiation between the minAabb, maxAabb and minAabbRotated, maxAabbRotated stems from the need to satisfy two conditions: 1. The Aabb used for primary contact detection must be as tight as possible, in order to have the least number of empty contacts and 2. The Aabb used as a rendering space must be slightly larger, in order to have proper rendering. If a dimension of the Aabb used for visualization purposes is defined smaller than the actual one, the faces on that side of the particle are rendered as hollow and only the edges are visualised, a functionality that can be used to e.g. see through boundaries, like demonstrated in the vtk output of the examples/PotentialParticles/cubePPscaled.py example.

To recap, in the Potential Particles code, the minAabbRotated and maxAabbRotated parameters are used to define the Aabb used to facilitate primary contact detection, while the minAabb and maxAabb parameters are used for visualization of the particles in qt and vtk outputs. In the Potential Blocks code, the Aabb used to facilitate primary contact detection is calculated automatically from the particles’ vertices, which are also used for visualization in qt, while the parameters minAabb and maxAabb are used for visualization in vtk outputs.

Two brief examples demonstrating the syntax of these features can be found below.

For the Potential Particles code:

b=Body()
b.shape=PotentialParticle(AabbMinMax=True,
minAabbRotated=Vector3(xmin,ymin,zmin),
maxAabbRotated=Vector3(xmax,ymax,zmax),
minAabb=Vector3(xmin,ymin,zmin),
maxAabb=Vector3(xmax,ymax,zmax), ...)


For the Potential Blocks code:

b=Body()
b.shape=PotentialBlock(AabbMinMax=True,
minAabb=Vector3(xmin,ymin,zmin),
maxAabb=Vector3(xmax,ymax,zmax), ...)


## Block Generation algorithm¶

The Potential Blocks code is compatible with the Block Generation algorithm introduced in [Boon2015], which can split particles by their intersection with discontinuity planes, initially developed for the study of rock-masses. This code is hardcoded in YADE in the form of a Preprocessor. Using a single data structure for the definition of the particle shape and the definition of the discontinuities, as well, allows the generation of a large number of particles at a reasonable computational cost. The sequential subdivision concept is used along with a linear programming framework. Non-persistent joints can be modelled by introducing more constraints.

An example to demonstrate the usage of this code exists in examples/PotentialBlocks/WedgeYADE.py The discontinuity planes used in this script are included in a csv format in examples/PotentialBlocks/joints/jointC.csv.

The documentation on how to use this code is currently being written.

## Examples¶

Examples are being included in the folders examples/PotentialParticles and examples/PotentialBlocks/, where the syntax of the codes is demonstrated.

## Disclaimer¶

These codes were developed for academic purposes. Some variables are no longer in use, as the PhD thesis of the original developer spanned over many years, with numerous trials and errors. As this piece of code has many dependencies within the YADE ecosystem, user discretion is advised.

## References¶

To acknowledge our scientific contribution, please cite the following:

$$\underline{\textrm{Potential Blocks}}$$

• Boon CW (2013) Distinct Element Modelling of Jointed Rock Masses: Algorithms and Their Verification. D.Phil. Thesis, University of Oxford
• Boon CW, Houlsby GT, Utili S (2012) A new algorithm for contact detection between convex polygonal and polyhedral particles in the discrete element method. Computers and Geotechnics, 44: 73-82

$$\underline{\textrm{Potential Particles}}$$

• Houlsby GT (2009) Potential particles: a method for modelling non-circular particles in DEM. Computers and Geotechnics, 36(6):953-959
• Boon CW, Houlsby GT, Utili S (2013) A new contact detection algorithm for three dimensional non-spherical particles. Powder Technology, S.I. on DEM, 248: 94-102

$$\underline{\textrm{Block Generation}}$$

• Boon CW, Houlsby GT, Utili S (2015) A new rock slicing method based on linear programming. Computers and Geotechnics, 65:12-29