From c3696f22c404a88e8bbbe7c3d91458cdcb25a8e2 Mon Sep 17 00:00:00 2001
From: Vasileios Angelidakis
Date: Thu, 30 Apr 2020 16:27:34 +0000
Subject: [PATCH] Small docs fix

doc/sphinx/potentialparticles.rst  852 +++++++++++++++
1 file changed, 426 insertions(+), 426 deletions()
diff git a/doc/sphinx/potentialparticles.rst b/doc/sphinx/potentialparticles.rst
index 6677965c5..4c3d9d5ed 100644
 a/doc/sphinx/potentialparticles.rst
+++ b/doc/sphinx/potentialparticles.rst
@@ 1,426 +1,426 @@
.. _PotentialParticles:

########################################
Potential Particles and Potential Blocks
########################################
The origins of scientific development regarding the algorithms described in this section are traced back to:
[Boon2012]_ (*Potential Blocks* code), [Boon2013b]_ (*Potential Particles* code) and [Boon2015]_ (*Block Generation* code).

Introduction
^^^^^^^^^^^^
This section discusses two codes to simulate (i) nonspherical particles using the concept of the Potential Particles [Houlsby2009]_, with the solution procedures in [Boon2013]_ for 3D and (ii) polyhedral blocks using planar linear inequalities, based on linear programming concepts [Boon2012]_. These codes define two shape classes in YADE, namely :yref:`PotentialParticle` and :yref:`PotentialBlock`. 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 nonspherical particles are assembled as a combination of 2\ :sup:`nd` degree polynomial functions and a fraction of a sphere, while their edges are rounded with a userdefined 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 :yref:`Block Generation` algorithm defined in [Boon2015]_, in which Potential Blocks can be generated by intersections of original, intact blocks with discontinuity planes.

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 PP and a PB). Enabling the PB code causes an automatic compilation of the :yref:`Block Generation` algorithm.

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 NewtonRaphson method in 2D. In [Boon2013]_ it was extended to 3D and more robust solutions were proposed. Many numerical optimization solvers generally cannot cope with discontinuities, illconditioned gradients (Jacobians) or curvatures (Hessians), and these obstacles were overcome in [Boon2013]_, by reformulating 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 :eq:`ppFormula` [Houlsby2009]_:

 .. math:: f=(1k) \Bigg(\sum_{i=1}^{N} \langle a_{i}x+b_{i}y+c_{i}zd_i\rangle ^2r^2 \Bigg) +k (x^2+y^2+z^2R^2)\\
 :label: ppFormula

where :math:`(a_i, b_i, c_i)` is the normal vector of the :math:`i^{th}` plane, defined with respect to the particle's local coordinate system and :math:`d_i` is the distance of the plane to the local origin. :math:`\langle \;\rangle` are Macaulay brackets, i.e., :math:`〈x〉 = x` for :math:`x > 0`; :math:`\langle x \rangle = 0` for :math:`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 :math:`0 < k \; \leq \; 1`, denoting the fraction of sphericity of the particle. The geometry of some cuboidal potential particles is displayed in Fig. `figpp`_, for different values of the parameter :math:`k`.

.. _figpp:
.. figure:: /fig/potentialParticle.*
 :width: 14cm
 :align: center

 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 :eq:`ppFormulaNormalized` [Houlsby2009]_:

 .. math:: f=(1k) \Bigg(\sum_{i=1}^{N} \frac{ \langle a_{i}x+b_{i}y+c_{i}zd_i\rangle^2 }{ r^2 } 1 \Bigg) +k \Bigg( \frac{ x^2+y^2+z^2 }{ R^2 }1 \Bigg)\\
 :label: ppFormulaNormalized

This potential function takes values:

* :math:`f=0`: on the particle surface
* :math:`f<0`: inside the particle
* :math:`f>0`: outside the particle

To ensure numerical stability, it is not advised to use values approaching *k=0*. In particular, the extreme value *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 :yref:`PotentialParticle` shape class, no spherical term is considered (so, practically k=0). Although *k* and *R* are input parameters of the :yref:`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 bisection step size for line search, to obtain a point on the particle, which in turn is used to calculate the overlap distance during contact.

.. Its value is suggested to be set to an order or magnitude near the sieve size of the particle, since a small value for *R* (relatively to the particle dimensions) will create warnings during contact detection.

For a convex particle defined by *N* planes, the space that it occupies can be defined using the following inequalities :eq:`pbOriginalPlaneEquations`:

 .. math:: a_{i}x + b_{i}y + c_{i}z \; \leq \; d_{i}, i=1:N\\
 :label: pbOriginalPlaneEquations

where :math:`(a_i, b_i, c_i)` is the unit normal vector of the :math:`i^{th}` plane, defined with respect to the particle's local coordinate system, and :math:`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 :eq:`pbOriginalFormula`:

 .. math:: f=\sum_{i=1}^{N} \langle a_{i}x + b_{i}y + c_{i}z  d_i + r\rangle^2\\
 :label: pbOriginalFormula

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. `figpbInner`_).

.. _figpbInner:
.. figure:: /fig/potentialBlockInner.*
 :width: 8.5cm
 :align: center

 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 :eq:`pbYADEPlaneEquations`:

 .. math:: a_{i}x + b_{i}y + c_{i}z  d_{i}  r = 0, i=1:N\\
 :label: pbYADEPlaneEquations

while the inner *Potential Particle* used to calculate the contact normal is defined as in :eq:`pbYADEFormula`:

 .. math:: f=\sum_{i=1}^{N} \langle a_{i}x + b_{i}y + c_{i}z  d_i\rangle^2.\\
 :label: pbYADEFormula

Now, the *Potential Block* surface is at a distance of :math:`(d_{i}+r)` from the local particle center, while the inner potential particle is at a distance :math:`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 wellknown concept of interiorpoint 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 bisection searching along the contact normal, within the overlap region.

.. _figpb:
.. figure:: /fig/potentialBlock.*
 :width: 15cm
 :align: center

 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 COINOR, since the latter can be downloaded from Ubuntu or Debian’s distributions without requiring an academic licence.

.. More topics for future development of the documentation, derived from the thesis of CW Boon (2013):

.. Contact detection (formulation of the optimisation problem)
.. Contact point (analytic centre of the inequalities defining the particles)
.. Contact normal (partial derivatives) \& average normal vector
.. Contact forces (after Hart et al, 1988)

.. _engines:

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 :yref:`FrictMat`. The syntax used to invoke these classes is similar, unless if specified otherwise.

================== ======================================== ========================================
Shape :yref:`PotentialParticle` :yref:`PotentialBlock`
================== ======================================== ========================================
Material :yref:`FrictMat` :yref:`FrictMat`
BoundFunctor :yref:`PotentialParticle2AABB` :yref:`PotentialBlock2AABB`
IGeom :yref:`ScGeom` :yref:`ScGeom`
IGeomFunctor :yref:`Ig2_PP_PP_ScGeom` :yref:`Ig2_PB_PB_ScGeom`
IPhys :yref:`KnKsPhys` :yref:`KnKsPBPhys`
IPhysFunctor :yref:`Ip2_FrictMat_FrictMat_KnKsPhys` :yref:`Ip2_FrictMat_FrictMat_KnKsPBPhys`
LawFunctor :yref:`Law2_SCG_KnKsPhys_KnKsLaw` :yref:`Law2_SCG_KnKsPBPhys_KnKsPBLaw`
VTK Recorder :yref:`PotentialParticleVTKRecorder` :yref:`PotentialBlockVTKRecorder`
================== ======================================== ========================================

A simple *simulation loop* using the *Potential Blocks* reads as:

.. codeblock:: python

 O.engines=[
 ForceResetter(),
 InsertionSortCollider([PotentialBlock2AABB()], verletDist=0.01),
 InteractionLoop(
 [Ig2_PB_PB_ScGeom(twoDimension=True, unitWidth2D=1.0, calContactArea=True)],
 [Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=1e8, ks_i=1e7, Knormal=1e8, Kshear=1e7, viscousDamping=0.2)],
 [Law2_SCG_KnKsPBPhys_KnKsPBLaw(label='law', neverErase=False, allowViscousAttraction=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 :yref:`twoDimension` parameter, which defines whether a contact should be handled as 2D or 3D.

Contact Law
^^^^^^^^^^^
In both codes, the normal force is calculated as:

 .. math:: \mathbf{F_{n}}=Knormal \cdot A_{c} \cdot u_{n} \cdot \mathbf{n}\\
 :label: normalForce

where :math:`Knormal` the normal stiffness coefficient :math:`[kN/m^{3}]`; :math:`A_{c}` the contact area :math:`[m^{2}]` and :math:`u_{n}` the overlap distance.
The normal stiffness of each contact :math:`[kN/m]` is thus :math:`k_{n} = Knormal \cdot A_{c}`, where :math:`A_{c}` is updated in every timestep.

The shear force is calculated incrementally, using a similar logic. The increment of the shear force vector before slipping of the contact is calculated as:

 .. math:: \mathbf{\Delta{}F_{s}}=Kshear \cdot A_{c} \cdot \mathbf{\Delta{}u_{s}}\\
 :label: shearForce

where :math:`Kshear` the shear stiffness coefficient :math:`[kN/m^{3}]` and :math:`\Delta{}u_{s}` the current relative shear displacement.

Contact Area

The contact area is calculated using a heuristic algorithm to detect points on the surface of the overlap volume, searching along the contact shear direction. In essence, it is calculated as the area of a 2D slice of the overlap volume along the shear direction, passing from the contact point. If ``twoDimension=True``, the :yref:`contactArea` parameter is calculated as::

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

The :yref:`unitWidth2D` parameter is defined by the user (usually equal to 1.0), denoting the outofplane width in 2D simulations. The :yref:`contactArea` and :yref:`jointLength` parameters are calculated if :yref:`calContactArea` ``=True``. In the opposite case, they are considered equal to 1.0 and the contact law is degenerated to a linear law with constant stiffness. A minimum value is considered for the :yref:`contactArea`, to represent cases where the overlap volume is practically a point.

Overlap distance

The overlap distance :math:`u_{n}` is calculated using a bracketed bisection search algorithm along the contact normal direction, to find two opposite points on the surface of the overlap region, starting from the contact point. It is stored in the parameter :yref:`penetrationDepth`, as the distance between these two opposite points.

Shape 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 faceface, faceedge, facevertex, edgeedge, edgevertex or vertexvertex 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 set:

.. codeblock:: python

 r=D/10.
 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 :math:`a, b, c, d` refers to the normal vector of the first plane and its distance from the local origin, 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 :eq:`ppFormulaNormalized`.
It is reminded that within this potential function, these planes are summed quadratically, the particle edges are rounded by a radius of curvature :yref:`r` and then the particle faces are curved by the addition of a "shadow" spherical particle with a radius :yref:`R`, to a percentage defined by the parameter :yref:`k`.
A value :yref:`r` is deducted from each element of the vector parameter :yref:`d`, to compensate for expanding the potential particle by :yref:`r in :eq:`ppFormulaNormalized`.

The parameters :math:`a_{i}, b_{i}, c_{i}, d_{i}` stated above correspond to the planes used in :eq:`pbYADEPlaneEquations`:

 .. math:: 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\\

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

.. codeblock:: python

 r=D/10.
 R=D/2.*sqrt(3)
 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 :eq:`pbYADEFormula`, located fully inside the actual particle.
The distances of the planes from the local origin, stored in the vector parameter :yref:`d`, are reduced by :yref:`r` to achieve an exact edge length of *D*, using :eq:`pbYADEPlaneEquations`. The value of :yref:`r` must be sufficiently small, so that :math:`d_{min}r>0`, while it should be sufficiently large, to allow for a proper calculation of the gradient of the inner Potential Particle at the contact point. A recommended value is :math:`r \approx 0.5*d_{min}`.

To ensure numerical stability, it is advised to normalize the normal vector of each plane, so that :math:`{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 :math:`a_{i}, b_{i}, c_{i}` and using a vertex of each face, then calculate the coefficients :math:`d_{i}`. A python routine to perform this without any additional effort by the user is currently being developed.



Body definition of a PP and a PB
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
To define a body using the :yref:`PotentialParticle` or :yref:`PotentialBlock` shape classes, it has to be assembled using the ``_commonBodySetup`` function, which can be found in the file :ysrc:`py/utils.py`. For example, to define a :yref:`PotentialParticle`:

.. codeblock:: python

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

 b=Body()
 b.shape=PotentialParticle(...)
 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 :yref:`PotentialParticle` must be initially defined, so that the local axes coincide with its principal axes, for which the inertia tensor is diagonal. More specifically, the plane coefficients :math:`(a_i, b_i, c_i)` defining the plane normals must be rotated, so that when the orientation of the particle is zero, the :yref:`PotentialParticle` is oriented to its principal axes.

It should be noted that the principal inertia values ``I11, I22, I33`` mentioned here are divided with the density of the considered material, since they 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``.

For the *Potential Particles*, the volume and inertia must be calculated manually and assigned to the body as demonstrated above. For the *Potential Blocks*, an automatic calculation has been implemented for the volume and inertia tensor, the user does not have to define the particle to its principal axes, since this is handled automatically within the source code, while if no value is given for the parameter :yref:`R`, it is calculated as half the distance of the farthest vertices.

For example, to define a :yref:`PotentialBlock`:

.. codeblock:: python

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

 b=Body()
 b.shape=PotentialBlock(R=0.0, ...) #here we set R=0.0 to trigger automatic calculation of R
 b.aspherical=True # To be used in conjunction with exactAsphericalRot=True
 utils._commonBodySetup(b,b.shape.volume,b.shape.inertia, material='frictionless', pos=Vector3(xPos,yPos,zPos), ori=Quaternion((1,0,0),0), fixed=False)
 b.state.ori=b.shape.orientation # this will rotate the particle to its initial random system. If b.state.ori=Quaternion.Identity, the PB is oriented to its principal axes
 O.bodies.append(b)


Boundary Particles
^^^^^^^^^^^^^^^^^^
The PP \& PB codes support the definition of *boundary* particles, which interact only with *nonboundary* 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 :yref:`isBoundary` inside the shape class.

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

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



.. _visualization:

Visualization
^^^^^^^^^^^^^
Visualization of the :yref:`PotentialParticle` and :yref:`PotentialBlock` shape classes is offered using the qt environment (OpenGL). Additionally, the :yref:`yade.export.VTKExporter.exportPotentialBlocks` function and :yref:`PotentialParticleVTKRecorder` and :yref:`PotentialBlockVTKRecorder` engines 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 :yref:`PotentialBlockVTKRecorder` records the inner, rounded potential particle, rather than the actual particle with sharp edges and flat faces.

In the qt environment, the :yref:`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:

.. codeblock:: python

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

.. # Potential Blocks
.. Gl1_PotentialBlock.sizeX=20
.. Gl1_PotentialBlock.sizeY=20
.. Gl1_PotentialBlock.sizeZ=20

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

.. codeblock:: python

 # 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 :yref:`maxDimension`, which overrides the selected sampleX,Y,Z values if they are too small, as described below:

.. math::
 if \; \mid xmaxxmin \mid /sampleX > maxDimension \Rightarrow sampleX = \mid xmaxxmin \mid /maxDimension\\
 if \; \mid ymaxymin \mid /sampleY > maxDimension \Rightarrow sampleY = \mid ymaxymin \mid /maxDimension\\
 if \; \mid zmaxzmin \mid /sampleZ > maxDimension \Rightarrow sampleZ = \mid zmaxzmin \mid /maxDimension \; \\

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

.. codeblock:: python

 # 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 :yref:`VTKRecorder` or the :yref:`yade.export.VTKExporter` class.

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

.. codeblock:: python

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

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

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



AxisAligned Bounding Box
^^^^^^^^^^^^^^^^^^^^^^^^^
The PP \& PB codes use their own BoundFunctors, called :yref:`PotentialParticle2AABB` and :yref:`PotentialBlock2AABB`, respectively, to define the AxisAligned Bounding Box of each particle. In both bound functors, a boolean parameter :yref:`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 :yref:`R` is used as a reference length, denoting half the diagonal of the cubic Aabb. Usage of this approximate cubic Aabb is not advised in general, since it can increase the number of empty contacts, adding thus to the time needed to facilitate the approximate contact detection, while it relies on the radius :yref:`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 :yref:`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: :yref:`minAabbRotated`, :yref:`maxAabbRotated` inside the shape class.
The Aabb for a :yref:`PotentialBlock`, on the other hand, is calculated and updated automatically from the vertices of the particle, if the boolean parameter :yref:`AabbMinMax` ``=True``.

As discussed in the subsection :ref:`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: :yref:`minAabb` and :yref:`maxAabb`. For the *Potential Blocks* code only, if these parameters are left unassigned, the drawing space is configured automatically inside the :yref:`PotentialBlockVTKRecorder` using the Aabb of the particle. 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 :ysrc:`examples/PotentialParticles/cubePPscaled.py` example.

To recap, in the *Potential Particles* code, the :yref:`minAabbRotated` and :yref:`maxAabbRotated` parameters define the initial Aabb used to facilitate primary contact detection, while the :yref:`minAabb` and :yref:`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 :yref:`minAabb` and :yref:`maxAabb` are used for visualization in vtk outputs and can be left unassigned, to trigger an automatic configuration of the drawing space of the particle in the :yref:`PotentialBlockVTKRecorder`.

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

For the *Potential Particles* code:

.. codeblock:: python

 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:

.. codeblock:: python

 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 :yref:`Block Generation` algorithm introduced in [Boon2015]_, which can split particles by their intersection with discontinuity planes, initially developed for the study of rockmasses. 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.
Nonpersistent joints can be modelled by introducing more constraints.

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

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



Examples
^^^^^^^^^^^^
Examples can be found in the folders :ysrc:`examples/PotentialParticles` and :ysrc:`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:

:math:`\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: 7382

:math:`\underline{\textrm{Potential Particles}}`

 Houlsby GT (2009) Potential particles: a method for modelling noncircular particles in DEM. Computers and Geotechnics, 36(6):953959
 Boon CW, Houlsby GT, Utili S (2013) A new contact detection algorithm for three dimensional nonspherical particles. Powder Technology, S.I. on DEM, 248: 94102

:math:`\underline{\textrm{Block Generation}}`

 Boon CW, Houlsby GT, Utili S (2015) A new rock slicing method based on linear programming. Computers and Geotechnics, 65:1229
+.. _PotentialParticles:
+
+########################################
+Potential Particles and Potential Blocks
+########################################
+The origins of scientific development regarding the algorithms described in this section are traced back to:
+[Boon2012]_ (*Potential Blocks* code), [Boon2013b]_ (*Potential Particles* code) and [Boon2015]_ (*Block Generation* code).
+
+Introduction
+^^^^^^^^^^^^
+This section discusses two codes to simulate (i) nonspherical particles using the concept of the Potential Particles [Houlsby2009]_, with the solution procedures in [Boon2013]_ for 3D and (ii) polyhedral blocks using planar linear inequalities, based on linear programming concepts [Boon2012]_. These codes define two shape classes in YADE, namely :yref:`PotentialParticle` and :yref:`PotentialBlock`. 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 nonspherical particles are assembled as a combination of 2\ :sup:`nd` degree polynomial functions and a fraction of a sphere, while their edges are rounded with a userdefined 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 :yref:`Block Generation` algorithm defined in [Boon2015]_, in which Potential Blocks can be generated by intersections of original, intact blocks with discontinuity planes.
+
+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 PP and a PB). Enabling the PB code causes an automatic compilation of the :yref:`Block Generation` algorithm.
+
+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 NewtonRaphson method in 2D. In [Boon2013]_ it was extended to 3D and more robust solutions were proposed. Many numerical optimization solvers generally cannot cope with discontinuities, illconditioned gradients (Jacobians) or curvatures (Hessians), and these obstacles were overcome in [Boon2013]_, by reformulating the problem and solving the equations using conic optimization solvers. Previous versions made use of MOSEK (using its academic licence), while currently an inhouse code written by [Boon2013]_ is used to solve the conic optimization problem.
+A potential particle is defined as in :eq:`ppFormula` [Houlsby2009]_:
+
+ .. math:: f=(1k) \Bigg(\sum_{i=1}^{N} \langle a_{i}x+b_{i}y+c_{i}zd_i\rangle ^2r^2 \Bigg) +k (x^2+y^2+z^2R^2)\\
+ :label: ppFormula
+
+where :math:`(a_i, b_i, c_i)` is the normal vector of the :math:`i^{th}` plane, defined with respect to the particle's local coordinate system and :math:`d_i` is the distance of the plane to the local origin. :math:`\langle \;\rangle` are Macaulay brackets, i.e., :math:`〈x〉 = x` for :math:`x > 0`; :math:`\langle x \rangle = 0` for :math:`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 :math:`0 < k \; \leq \; 1`, denoting the fraction of sphericity of the particle. The geometry of some cuboidal potential particles is displayed in Fig. `figpp`_, for different values of the parameter :math:`k`.
+
+.. _figpp:
+.. figure:: /fig/potentialParticle.*
+ :width: 14cm
+ :align: center
+
+ 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 :eq:`ppFormulaNormalized` [Houlsby2009]_:
+
+ .. math:: f=(1k) \Bigg(\sum_{i=1}^{N} \frac{ \langle a_{i}x+b_{i}y+c_{i}zd_i\rangle^2 }{ r^2 } 1 \Bigg) +k \Bigg( \frac{ x^2+y^2+z^2 }{ R^2 }1 \Bigg)\\
+ :label: ppFormulaNormalized
+
+This potential function takes values:
+
+* :math:`f=0`: on the particle surface
+* :math:`f<0`: inside the particle
+* :math:`f>0`: outside the particle
+
+To ensure numerical stability, it is not advised to use values approaching *k=0*. In particular, the extreme value *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 :yref:`PotentialParticle` shape class, no spherical term is considered (so, practically k=0). Although *k* and *R* are input parameters of the :yref:`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 bisection step size for line search, to obtain a point on the particle, which in turn is used to calculate the overlap distance during contact.
+
+.. Its value is suggested to be set to an order or magnitude near the sieve size of the particle, since a small value for *R* (relatively to the particle dimensions) will create warnings during contact detection.
+
+For a convex particle defined by *N* planes, the space that it occupies can be defined using the following inequalities :eq:`pbOriginalPlaneEquations`:
+
+ .. math:: a_{i}x + b_{i}y + c_{i}z \; \leq \; d_{i}, i=1:N\\
+ :label: pbOriginalPlaneEquations
+
+where :math:`(a_i, b_i, c_i)` is the unit normal vector of the :math:`i^{th}` plane, defined with respect to the particle's local coordinate system, and :math:`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 :eq:`pbOriginalFormula`:
+
+ .. math:: f=\sum_{i=1}^{N} \langle a_{i}x + b_{i}y + c_{i}z  d_i + r\rangle^2\\
+ :label: pbOriginalFormula
+
+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. `figpbInner`_).
+
+.. _figpbInner:
+.. figure:: /fig/potentialBlockInner.*
+ :width: 8.5cm
+ :align: center
+
+ 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 :eq:`pbYADEPlaneEquations`:
+
+ .. math:: a_{i}x + b_{i}y + c_{i}z  d_{i}  r = 0, i=1:N\\
+ :label: pbYADEPlaneEquations
+
+while the inner *Potential Particle* used to calculate the contact normal is defined as in :eq:`pbYADEFormula`:
+
+ .. math:: f=\sum_{i=1}^{N} \langle a_{i}x + b_{i}y + c_{i}z  d_i\rangle^2.\\
+ :label: pbYADEFormula
+
+Now, the *Potential Block* surface is at a distance of :math:`(d_{i}+r)` from the local particle center, while the inner potential particle is at a distance :math:`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 wellknown concept of interiorpoint 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 bisection searching along the contact normal, within the overlap region.
+
+.. _figpb:
+.. figure:: /fig/potentialBlock.*
+ :width: 15cm
+ :align: center
+
+ 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 COINOR, since the latter can be downloaded from Ubuntu or Debian’s distributions without requiring an academic licence.
+
+.. More topics for future development of the documentation, derived from the thesis of CW Boon (2013):
+
+.. Contact detection (formulation of the optimisation problem)
+.. Contact point (analytic centre of the inequalities defining the particles)
+.. Contact normal (partial derivatives) \& average normal vector
+.. Contact forces (after Hart et al, 1988)
+
+.. _engines:
+
+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 :yref:`FrictMat`. The syntax used to invoke these classes is similar, unless if specified otherwise.
+
+================== ======================================== ========================================
+Shape :yref:`PotentialParticle` :yref:`PotentialBlock`
+================== ======================================== ========================================
+Material :yref:`FrictMat` :yref:`FrictMat`
+BoundFunctor :yref:`PotentialParticle2AABB` :yref:`PotentialBlock2AABB`
+IGeom :yref:`ScGeom` :yref:`ScGeom`
+IGeomFunctor :yref:`Ig2_PP_PP_ScGeom` :yref:`Ig2_PB_PB_ScGeom`
+IPhys :yref:`KnKsPhys` :yref:`KnKsPBPhys`
+IPhysFunctor :yref:`Ip2_FrictMat_FrictMat_KnKsPhys` :yref:`Ip2_FrictMat_FrictMat_KnKsPBPhys`
+LawFunctor :yref:`Law2_SCG_KnKsPhys_KnKsLaw` :yref:`Law2_SCG_KnKsPBPhys_KnKsPBLaw`
+VTK Recorder :yref:`PotentialParticleVTKRecorder` :yref:`PotentialBlockVTKRecorder`
+================== ======================================== ========================================
+
+A simple *simulation loop* using the *Potential Blocks* reads as:
+
+.. codeblock:: python
+
+ O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([PotentialBlock2AABB()], verletDist=0.01),
+ InteractionLoop(
+ [Ig2_PB_PB_ScGeom(twoDimension=True, unitWidth2D=1.0, calContactArea=True)],
+ [Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=1e8, ks_i=1e7, Knormal=1e8, Kshear=1e7, viscousDamping=0.2)],
+ [Law2_SCG_KnKsPBPhys_KnKsPBLaw(label='law', neverErase=False, allowViscousAttraction=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 :yref:`twoDimension` parameter, which defines whether a contact should be handled as 2D or 3D.
+
+Contact Law
+^^^^^^^^^^^
+In both codes, the normal force is calculated as:
+
+ .. math:: \mathbf{F_{n}}=Knormal \cdot A_{c} \cdot u_{n} \cdot \mathbf{n}\\
+ :label: normalForce
+
+where :math:`Knormal` the normal stiffness coefficient :math:`[kN/m^{3}]`; :math:`A_{c}` the contact area :math:`[m^{2}]` and :math:`u_{n}` the overlap distance.
+The normal stiffness of each contact :math:`[kN/m]` is thus :math:`k_{n} = Knormal \cdot A_{c}`, where :math:`A_{c}` is updated in every timestep.
+
+The shear force is calculated incrementally, using a similar logic. The increment of the shear force vector before slipping of the contact is calculated as:
+
+ .. math:: \mathbf{\Delta{}F_{s}}=Kshear \cdot A_{c} \cdot \mathbf{\Delta{}u_{s}}\\
+ :label: shearForce
+
+where :math:`Kshear` the shear stiffness coefficient :math:`[kN/m^{3}]` and :math:`\Delta{}u_{s}` the current relative shear displacement.
+
+Contact Area
+
+The contact area is calculated using a heuristic algorithm to detect points on the surface of the overlap volume, searching along the contact shear direction. In essence, it is calculated as the area of a 2D slice of the overlap volume along the shear direction, passing from the contact point. If ``twoDimension=True``, the :yref:`contactArea` parameter is calculated as::
+
+ if(twoDimension) { phys>contactArea = phys>jointLength*unitWidth2D;}
+
+The :yref:`unitWidth2D` parameter is defined by the user (usually equal to 1.0), denoting the outofplane width in 2D simulations. The :yref:`contactArea` and :yref:`jointLength` parameters are calculated if :yref:`calContactArea` ``=True``. In the opposite case, they are considered equal to 1.0 and the contact law is degenerated to a linear law with constant stiffness. A minimum value is considered for the :yref:`contactArea`, to represent cases where the overlap volume is practically a point.
+
+Overlap distance
+
+The overlap distance :math:`u_{n}` is calculated using a bracketed bisection search algorithm along the contact normal direction, to find two opposite points on the surface of the overlap region, starting from the contact point. It is stored in the parameter :yref:`penetrationDepth`, as the distance between these two opposite points.
+
+Shape 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 faceface, faceedge, facevertex, edgeedge, edgevertex or vertexvertex 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 set:
+
+.. codeblock:: python
+
+ r=D/10.
+ 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 :math:`a, b, c, d` refers to the normal vector of the first plane and its distance from the local origin, 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 :eq:`ppFormulaNormalized`.
+It is reminded that within this potential function, these planes are summed quadratically, the particle edges are rounded by a radius of curvature :yref:`r` and then the particle faces are curved by the addition of a "shadow" spherical particle with a radius :yref:`R`, to a percentage defined by the parameter :yref:`k`.
+A value :yref:`r` is deducted from each element of the vector parameter :yref:`d`, to compensate for expanding the potential particle by :yref:`r in :eq:`ppFormulaNormalized`.
+
+The parameters :math:`a_{i}, b_{i}, c_{i}, d_{i}` stated above correspond to the planes used in :eq:`pbYADEPlaneEquations`:
+
+ .. math:: 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\\
+
+To model a cube with an edge of *D*, using the *Potential Blocks* code, we define:
+
+.. codeblock:: python
+
+ r=D/10.
+ R=D/2.*sqrt(3)
+ 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 :eq:`pbYADEFormula`, located fully inside the actual particle.
+The distances of the planes from the local origin, stored in the vector parameter :yref:`d`, are reduced by :yref:`r` to achieve an exact edge length of *D*, using :eq:`pbYADEPlaneEquations`. The value of :yref:`r` must be sufficiently small, so that :math:`d_{min}r>0`, while it should be sufficiently large, to allow for a proper calculation of the gradient of the inner Potential Particle at the contact point. A recommended value is :math:`r \approx 0.5*d_{min}`.
+
+To ensure numerical stability, it is advised to normalize the normal vector of each plane, so that :math:`{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 :math:`a_{i}, b_{i}, c_{i}` and using a vertex of each face, then calculate the coefficients :math:`d_{i}`. A python routine to perform this without any additional effort by the user is currently being developed.
+
+
+
+Body definition of a PP and a PB
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+To define a body using the :yref:`PotentialParticle` or :yref:`PotentialBlock` shape classes, it has to be assembled using the ``_commonBodySetup`` function, which can be found in the file :ysrc:`py/utils.py`. For example, to define a :yref:`PotentialParticle`:
+
+.. codeblock:: python
+
+ O.materials.append(FrictMat(young=1,poisson=1,frictionAngle=radians(0.0),density=2650,label='frictionless'))
+
+ b=Body()
+ b.shape=PotentialParticle(...)
+ 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), 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 :yref:`PotentialParticle` must be initially defined, so that the local axes coincide with its principal axes, for which the inertia tensor is diagonal. More specifically, the plane coefficients :math:`(a_i, b_i, c_i)` defining the plane normals must be rotated, so that when the orientation of the particle is zero, the :yref:`PotentialParticle` is oriented to its principal axes.
+
+It should be noted that the principal inertia values ``I11, I22, I33`` mentioned here are divided with the density of the considered material, since they 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``.
+
+For the *Potential Particles*, the volume and inertia must be calculated manually and assigned to the body as demonstrated above. For the *Potential Blocks*, an automatic calculation has been implemented for the volume and inertia tensor, the user does not have to define the particle to its principal axes, since this is handled automatically within the source code, while if no value is given for the parameter :yref:`R`, it is calculated as half the distance of the farthest vertices.
+
+For example, to define a :yref:`PotentialBlock`:
+
+.. codeblock:: python
+
+ O.materials.append(FrictMat(young=1,poisson=1,frictionAngle=radians(0.0),density=2650,label='frictionless'))
+
+ b=Body()
+ b.shape=PotentialBlock(R=0.0, ...) #here we set R=0.0 to trigger automatic calculation of R
+ b.aspherical=True # To be used in conjunction with exactAsphericalRot=True
+ utils._commonBodySetup(b,b.shape.volume,b.shape.inertia, material='frictionless', pos=Vector3(xPos,yPos,zPos), fixed=False)
+ b.state.ori=b.shape.orientation # this will rotate the particle to its initial random system. If b.state.ori=Quaternion.Identity, the PB is oriented to its principal axes
+ O.bodies.append(b)
+
+
+Boundary Particles
+^^^^^^^^^^^^^^^^^^
+The PP \& PB codes support the definition of *boundary* particles, which interact only with *nonboundary* 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 :yref:`isBoundary` inside the shape class.
+
+In the PP code, all particles interact with the same normal and shear contact stiffness :yref:`Knormal` and :yref:`Kshear`, defined in the :yref:`Ip2_FrictMat_FrictMat_KnKsPhys` functor.
+
+The PB code supports the definition of different contact stiffness values for interactions between *boundary* and *nonboundary* or *nonboundary* and *nonboundary* particles.
+When ``isBoundary=False``, the :yref:`PotentialBlock` in question is handled to interact with normal and shear stiffness coefficients :yref:`Knormal` and :yref:`Kshear`, respectively, with other nonboundary particles.
+When ``isBoundary=True``, the :yref:`PotentialBlock` in question is handled to interact with normal and shear stiffness coefficients :yref:`kn_i` and :yref:`ks_i`, respectively, with nonboundary particles.
+
+
+
+.. _visualization:
+
+Visualization
+^^^^^^^^^^^^^
+Visualization of the :yref:`PotentialParticle` and :yref:`PotentialBlock` shape classes is offered using the qt environment (OpenGL). Additionally, the :yref:`yade.export.VTKExporter.exportPotentialBlocks` function and :yref:`PotentialParticleVTKRecorder` and :yref:`PotentialBlockVTKRecorder` engines 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 :yref:`PotentialBlockVTKRecorder` records a rounded approximation of the particle, rather than the actual particle with sharp corners and edges.
+
+In the qt environment, the :yref:`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:
+
+.. codeblock:: python
+
+ # Potential Particles
+ Gl1_PotentialParticle.sizeX=20
+ Gl1_PotentialParticle.sizeY=20
+ Gl1_PotentialParticle.sizeZ=20
+
+.. # Potential Blocks
+.. Gl1_PotentialBlock.sizeX=20
+.. Gl1_PotentialBlock.sizeY=20
+.. Gl1_PotentialBlock.sizeZ=20
+
+A similar choice exists for output in vtk format, using the :yref:`PotentialParticleVTKRecorder` or :yref:`PotentialBlockVTKRecorder`, syntaxed as:
+
+.. codeblock:: python
+
+ # 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 :yref:`maxDimension`, which overrides the selected sampleX,Y,Z values if they are too small, as described below:
+
+.. math::
+ if \; \mid xmaxxmin \mid /sampleX > maxDimension \Rightarrow sampleX = \mid xmaxxmin \mid /maxDimension\\
+ if \; \mid ymaxymin \mid /sampleY > maxDimension \Rightarrow sampleY = \mid ymaxymin \mid /maxDimension\\
+ if \; \mid zmaxzmin \mid /sampleZ > maxDimension \Rightarrow sampleZ = \mid zmaxzmin \mid /maxDimension \; \\
+
+The :yref:`PotentialParticleVTKRecorder` and :yref:`PotentialBlockVTKRecorder` also support optionally the recording of the particles' velocities (linear and angular), interaction information (contact point and forces), colors and ids, using:
+
+.. codeblock:: python
+
+ # 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 :yref:`VTKRecorder` or the :yref:`yade.export.VTKExporter` class.
+
+A boolean parameter :yref:`twoDimension` exists to specify whether the particles will be rendered as 2D or 3D in the vtk output:
+
+.. codeblock:: python
+
+ # Potential Particles
+ PotentialParticleVTKRecorder(..., twoDimension=False)
+
+ # Potential Blocks
+ PotentialBlockVTKRecorder(..., twoDimension=False)
+
+This parameter should not be mixed up with the :yref:`Ip2_FrictMat_FrictMat_KnKsPBPhys.twoDimension` parameter, which is used to define how the contact forces are calculated, as described in the :ref:`engines` section.
+
+
+
+AxisAligned Bounding Box
+^^^^^^^^^^^^^^^^^^^^^^^^^
+The PP \& PB codes use their own BoundFunctors, called :yref:`PotentialParticle2AABB` and :yref:`PotentialBlock2AABB`, respectively, to define the AxisAligned Bounding Box of each particle. In both bound functors, a boolean parameter :yref:`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 :yref:`R` is used as a reference length, denoting half the diagonal of the cubic Aabb. Usage of this approximate cubic Aabb is not advised in general, since it can increase the number of empty contacts, adding thus to the time needed to facilitate the approximate contact detection, while it relies on the radius :yref:`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 :yref:`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: :yref:`minAabbRotated`, :yref:`maxAabbRotated` inside the shape class.
+The Aabb for a :yref:`PotentialBlock`, on the other hand, is calculated and updated automatically from the vertices of the particle, if the boolean parameter :yref:`AabbMinMax` ``=True``.
+
+As discussed in the subsection :ref:`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: :yref:`minAabb` and :yref:`maxAabb`. For the *Potential Blocks* code only, if these parameters are left unassigned, the drawing space is configured automatically inside the :yref:`PotentialBlockVTKRecorder` using the Aabb of the particle. 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 :ysrc:`examples/PotentialParticles/cubePPscaled.py` example.
+
+To recap, in the *Potential Particles* code, the :yref:`minAabbRotated` and :yref:`maxAabbRotated` parameters define the initial Aabb used to facilitate primary contact detection, while the :yref:`minAabb` and :yref:`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 :yref:`minAabb` and :yref:`maxAabb` are used for visualization in vtk outputs and can be left unassigned, to trigger an automatic configuration of the drawing space of the particle in the :yref:`PotentialBlockVTKRecorder`.
+
+Two brief examples demonstrating the syntax of these features can be found below.
+
+For the *Potential Particles* code:
+
+.. codeblock:: python
+
+ 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:
+
+.. codeblock:: python
+
+ 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 :yref:`Block Generation` algorithm introduced in [Boon2015]_, which can split particles by their intersection with discontinuity planes, initially developed for the study of rockmasses. 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.
+Nonpersistent joints can be modelled by introducing more constraints.
+
+An example to demonstrate the usage of this code exists in :ysrc:`examples/PotentialBlocks/WedgeYADE.py`
+The discontinuity planes used in this script are included in a csv format in :ysrc:`examples/PotentialBlocks/joints/jointC.csv`.
+
+The documentation on how to use this code is currently being written.
+
+
+
+Examples
+^^^^^^^^^^^^
+Examples can be found in the folders :ysrc:`examples/PotentialParticles` and :ysrc:`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:
+
+:math:`\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: 7382
+
+:math:`\underline{\textrm{Potential Particles}}`
+
+ Houlsby GT (2009) Potential particles: a method for modelling noncircular particles in DEM. Computers and Geotechnics, 36(6):953959
+ Boon CW, Houlsby GT, Utili S (2013) A new contact detection algorithm for three dimensional nonspherical particles. Powder Technology, S.I. on DEM, 248: 94102
+
+:math:`\underline{\textrm{Block Generation}}`
+
+ Boon CW, Houlsby GT, Utili S (2015) A new rock slicing method based on linear programming. Computers and Geotechnics, 65:1229

GitLab