Discrete element model for general polyhedra

We present a version of the Discrete Element Method considering the particles as rigid polyhedra. The Principle of Virtual Work is employed as basis for a multibody dynamics model. Each particle surface is split into sub-regions, which are tracked for contact with other sub-regions of neighboring particles. Contact interactions are modeled pointwise, considering vertex-face, edge-edge, vertex-edge and vertex-vertex interactions. General polyhedra with triangular faces are considered as particles, permitting multiple pointwise interactions which are automatically detected along the model evolution. We propose a combined interface law composed of a penalty and a barrier approach, to fulfill the contact constraints. Numerical examples demonstrate that the model can handle normal and frictional contact effects in a robust manner. These include simulations of convex and non-convex particles, showing the potential of applicability to materials with complex shaped particles such as sand and railway ballast.


Introduction
Modeling of granular materials is challenging, particularly when trying to represent continuum mechanics behaviour. Complex constitutive models appear on this context, usually involving many fitting parameters.
The Discrete Element Method (DEM) provides unique possibilities to handle micromechanical behavior of granular media (see [1] for the method origins, [2] for a useful review of spherical DEM modeling among other topics and [3,4] as textbooks in molecular dynamics, bringing many similarities and common strategies of DEM).
Numerous implementations of DEM exist, considering spherical particles due to their geometrical simplicity. Usually the computational bottleneck of a DEM solver is associated with spatial searching/treatment of contact between particles. Spherical particles are welcome in this context, since evaluating their overlap/proximity is simple and straightfor-B Peter Wriggers wriggers@ikm.uni-hannover.de Alfredo Gay Neto alfredo.gay@usp.br 1 University of São Paulo, São Paulo, Brazil 2 Leibniz Universität-Hannover, Hannover, Germany ward. When the particle shape deviates from a sphere, usually contact computation time is strongly increased.
Unfortunately spherical particle modeling is not sufficient for many applications. There are challenging practical problems of interest that motivate the development of DEM considering arbitrary-shaped particles, such as sandy soil and railway ballast mechanical interactions, and hopper discharge problems, just to mention a few examples that are heavily dependent on the particle shape. Indeed, mechanical macroscale behavior of such materials is evidenced in packing density, compressibility, critical state friction angle and other mechanical properties. According to [5], the macroscale properties result from particle interactions, which are affected by the particle shape.
When there is a need to consider non-spherical particles in DEM models, several alternatives can be utilized. Clustering spheres can be employed, see, e.g. [6]. In [7] normal impact of quasi-spherical particles is explored. Forming clusters of tetrahedra with distinct shapes is introduced in the work of [8]. Clusters of spheres are employed in the context of distinct integration techniques in [9] and, finally, in [10] a scheme is proposed that represents surfaces as spheres from a triangular surface mesh. Clusters of spheres have the advantage of the simplicity of local treatment of spherical particles, but always lack continuity in particle curvature and definition of tangents at the boundary surface. According to [11] this may lead to clumping ofparticles.
Another possibility is the usage of polyhedra to describe particles. In this case, surface singularities (vertices, edges) are present. This approach may be interesting when real particulate material consists of such polyhedra.
Modeling of DEM with polyhedra can be found in [12] and [13], which explore convex and non-convex (by convex decomposition) polyhedra simulation. In [14] convex polyhedra are addressed. In [15] and [16] efforts are presented to modeling convex and non-convex DEM-based polyhedra. A complementary alternative to DEM to model polyhedral particles (among other applications) is the usage of the Contact Dynamics Method (see, e.g. [17,18]). In this case, one solves for the velocity function, and imposes the balance of momentum at each time-step. Implicit schemes for time-integration are available, which permits employing large time-steps. A review of such method, including a discussion on comparisons with DEM and applications for polygons is found in [19].
When the particle shape is convex, a convenient and elegant description given by super-ellipsoids was proposed in [20] and applied to DEM simulations in [21]. Superellipsoids were also addressed in [22]. An initiative for dealing with the geometry of the particles in a more detailed way was presented in [23], which describes the called "granular element method", as a variant of DEM. It considers each grain's boundary as described by a non-uniform rational B-spline (NURBS) surface. This allows DEM to consider realistic and complex granular shapes. In [11] the non-convexity of particle shapes is addressed by a particular contact algorithm (knot-surface, in the context of non-uniform rational B-splines (NURBS) surfaces), which permits important advances in correlation with experimental results of sand specimen tests. Another effort to describe arbitrary-shaped particles is given in [24], which considers level set discrete elements.
When establishing a DEM model for applications in which particle shape plays a role, one can also approximate the physical behaviour by spherical models. However, in this case there is the main drawback that rolling resistance models are needed to achieve the expected macro mechanical properties, which may lead to problem-dependent calibration. For more complex particle shapes, particularly involving non-convex geometries, there is a need to handle contact detection and in a different, more time consuming way. It is highly desired, however, to minimize the need for model calibration parameters. This challenge was addressed in [25], which shows the importance of particle shapes when considering shear band evaluation is sand specimen by numerical simulations. In the context of railway ballast simulation, in [26] model parameters are discussed for shear tests, as also addressed in [27] and [28]. In [29] computational strategies are proposed in the context of railway ballast simulation to handle numerous contact interactions. Sand material simu-lation employing clusters of spheres was addressed in [30], which led to the effort in representing arbitrary shapes, even composed of spheres. Particle shape is important also for hopper discharge simulations. In [31] one can find comparison of clusters of spheres and polyhedra for hopper experiments, representing an interesting discussion on the influence of particle shape and its modeling issues.

Proposed DEM model
Motivated by the need for a realistic shape representation of particles in the aforementioned problems, a new methodology will be developed which includes general polyhedra with convex and non-convex shapes in a discrete element formulation. The model incorporates numerous pointwise contact interactions between particles, according to their shape and orientation. This is achieved by a strategy of division of the external surface of each particle into sub-regions which are described by faces, edges and vertices, as natural geometric components of general polyhedra. All of these entities are considered during the overall contact search.
Each pointwise contact interaction is addressed by degenerations of a specific surface-to-surface treatment, see [32,33], which is a novel approach in the DEM context. We provide a novel systematic usage of the master-to-master degenerated technique, particularly in the DEM context. This is done by the combination of several kinds of contact degenerations to handle the possibly numerous interactions between general polyhedra.
Particles are considered as rigid bodies. A novel treatment for the contact interface law is proposed. At the interface a hybrid normal direction law is formulated which is composed of a physically or numerically (penalty) ruled part and a Barrier-based part. The Barrier-based approach prevents penetration between particles and creates a thin layer of contact activation, borrowing ideas from molecular dynamics and collision detection in the context of computer graphics. The proposed method is tested with basic and complex examples to show its robustness and generality.

Nomenclature
In the context of a 3D Euclidean space, the present work uses the following nomenclature: scalar variables are non-bold (e.g.: v), vectors are lower-case bold (e.g.: v) and secondorder tensors are upper-case bold (e.g.: V ). Column matrices (termed as "vectors") are also lower-case bold and matrices with more than one column are also upper-case bold. Zero column-matrices are denoted by o s , where s in the number of rows.
The derivative of a quantity a with respect to a quantity b is denoted by a ,b . The variation of a quantity a is indicated by δa.

Model description
In present work, discrete elements are named "particles" or "bodies", with the same meaning. The time-evolution behavior of a DEM system is obtained by solving the equations of motion. For rigid bodies the equations of motion are written using the Newton-Euler description, see, e.g.: [34] and [3]. The dynamical behaviour of each rigid body is ruled by a set of six differential equations in time. In some applications rotational part of the motion can be neglected, especially for very small particles. In such case, only three differential equations of translational motion have to be used for each body. In present work translation and rotation are considered for the description of the particle motion while assuming rigid body behavior.
In general an analytically integration of the differential equations of motion is not possible due to complexity and frequent changes of contributions (forces/moments) acting on each body. Therefore, one usually has to solve the strongly coupled set of nonlinear differential equations by adopting a numerical scheme.
In the present work a weak form is employed to describe the equations of motion which is equivalent to Principle of Virtual Work (PVW). This is a common way to describe a multibody system in nonlinear solid mechanicsparticularly when employing numerical methods to consider the flexibility of bodies-such as using the Finite Element Method (FEM)-see, e.g.: [35] and [36].
The general weak form which describes a multibody system composed of rigid or flexible bodies is given by: where δW i is the total virtual work of internal forces, δW e is the total virtual work of external loads and the term δT describes the total inertial contribution to the model. All mechanical contact interactions are included in δW c . Additional terms can be considered, such as fluid loads. Each term in (1) contains contributions stemming from all bodies of the system. Let N be the number of particles in the system. With that, the contributions in (1) can be written as a sum where the superscript B denotes the virtual work contribution of body B. The term δW AB c is related to the virtual work of the contact interactions for a pair of bodies A and B.
Next, we will detail the methodology employed for the evaluation of each of the contributions in (2). Based on the rigid body assumption δW B i is equal to zero. Section 2.1 will provide details on the evaluation of δT B , Sect. 2.3 discusses the external loads and their contributions δW B e and Sect. 2.4 describes the methodology employed to address contact between bodies, leading to the contributions in δW AB c .

Dynamics of a single rigid body
Details of the contribution δT B related to each rigid body B are presented in this section. Rigid body motion and its causes are governed by Newton-Euler equations. As rotations have to be described in a 3D Euclidean space by distinct techniques (such as rotation vectors and quaternions). We will here present some details on these aspects, for completeness of the present work.

Kinematics
We adopt an updated Lagrangian description for the movement of material points. Let P be a generic material point in a body that is tracked along three successive configurations: the reference configuration r , the current configuration i and the next configuration i + 1. The position of the material point P is denoted by x r , x i and x i+1 for configurations r , i and i + 1, respectively. A general vector quantity a r is associated with the material point P at the reference configuration. This quantity experiences rigid body rotations when the configuration changes to i, leading to a i and to i + 1, leading to a i+1 .
The translatory motion is given by the kinematic relation where u Δ is the incremental displacement of the material point P from time t i to t i+1 . The evolution of the vector a r is described by the rotation tensor as The multiplicative decomposition of the rotation tensor holds permitting the description of partial rotations. The link between the rotation tensor operator and kinematic quantities to describe a rotation is provided by a vector parameterization, here given by the called Rodrigues rotation vector α. This description has been already employed in the context of beams and shell structural models, see, e.g.: [37][38][39][40][41][42]) and for rigid bodies in [43]. It was also employed in the context of spherical particles in [44] and [45].
Let a given rotation with magnitude θ about an axis be represented by the unit vector e described by the Euler rotation vector θ = θ e. For a Rodrigues rotation vector α = α e, where α = 2 tan (θ/2), the corresponding rotation tensor is given by where I is the identity tensor, α = α and A = skew(α). Equation (6) can applied to write rotation tensors Q i , Q Δ and Q i+1 . A convenient formula is available for updating the rotation vector directly, without usage of (5), The time-derivative of the vector a i+1 can be computed by considering the previous configuration vector a i as constant. This follows from the updated Lagrangian description where a fixed current configuration i can be assumed. The angular velocity of a i+1 is described by Ω = skew(ω). With some algebraic work a useful relation between the angular velocity and the time-derivative of the corresponding Rodrigues rotation vector can be found, see, e.g. [46], where with A Δ = skew(α Δ ).

Rigid body description
The herein utilized formulation was presented in a completed form in [43] and is here presented, briefly, for completeness. The reader interested in more details may refer to the original paper. We describe the movement of a rigid body B by a choice of a material point P, named "pole". This point does not need to be in the physical domain of the body, but it has to follow the rigid body constraints, that is, it experiences a rigid body movement as if it was a part of the body B. The center of mass of the body is defined as G. Based on the above notation the position of the center of mass is given by Assuming that the vector b r changes only its orientation but keeps its magnitude along time (rigid body assumption) yields where b r is a quantity evaluated at the reference configuration. It denotes the center of mass position with respect to the pole P. A generic material point within the body is given by x i+1 , such that: As all material points experience the same rigid body rotation, the update of s i+1 is given by Inertia tensor at configuration i + 1 is defined by where ρ is the body material specific mass and S i+1 = skew(s i+1 ). The inertia tensor J i+1 at configuration i + 1 is related to the inertia tensor at the reference configuration J r by Hence it can be pre-evaluated for any rigid body just using s r instead of s i+1 in (14). Note that we have to consider the pole P when evaluating J r and J i+1 . Thus, a particular (and convenient) choice for the pole is the center of mass, which simplifies the forthcoming equations. The kinetic energy T B of the rigid body B is By inserting (10)-(15) into (16) the first variation can be evaluated for the kinetic energy. Using furthermore the angular velocity relation presented (8) leads after some algebra to where Here m is the mass of the body B, ω is the body angular velocity,ω is its angular acceleration andü is the acceleration of the pole P. For later use, we define the inertial pseudo- The particular choice of the pole P as the center of mass simplifies (18) yielding which is a more standard representation of the equations of motion for a rigid body.
Note that the only kinematic variables employed in the evaluation of δT B are the incremental displacements of the pole u Δ P and the incremental Rodrigues rotation parameters α Δ P which are the chosen degrees of freedom (DOF). The time-derivatives of these quantities have to be computed, see (18). Moreover, updating of the center of mass position b i+1 and of the inertia tensor J i+1 is also necessary. Both are obtained by the applying the reference configuration values and the rotation tensor Q i+1 , as expressed in (11) and (15). Alternatively, one may save intermediate values at configuration i and then update employing the partial rotation Q Δ .
The inertia tensor is constant for a frame accompanying the rigid body motion, as employed in classical rigid-body derivations. However, the primary kinematic variable associated with rotation herein employed is the Rodrigues rotation vector, which is easily defined with respect to a fixed frame (differently from alternative rotation schemes that are naturally described following the body orientation). This choice makes more natural to work with the rotation tensor and the angular velocity vector written using the same basis associated with the fixed frame. This is the reason for what one must update the rotation tensor, accordingly.

Time-integration scheme
In present work, we use the Newmark integration scheme, see e.g. [47]. This is based on approximations for time-variable quantities as a function of the adopted time-step Δt. Here we assume that a time-step Δt governs the motion from the current configuration i and the next configuration i + 1. Within the Newmark method approximation formulae are where the parameters β and γ determine the behavior of the integration scheme. Depending on the choice numerical damping can be introduced or it is possible to generate an implicit/explicit scheme. Our choice is β = 0.3 and γ = 0.5, leading to an implicit method with small numerical damping. Equation (20) can be re-written such that velocity and acceleration of the unknown next configuration are given as a function of quantities from the previous configuration and the unknown current displacement u Δ P which follows from the equations of motion, see, e.g. [47] where the following constants are used: Δt. The proposed Newmark integration method can generally be applied for nonlinear problems, but lacks conservation of angular momentum. This drawback can be circumvented, see [48] and [49], where the Newmark method was modified for rotations in the context of simulation of beam-like structures and mechanisms involving joints together with flexible bodies, using Euler parameters. In [39] the modified Newmark method was reformulated for Rodrigues parameters which yieldṡ Equations (21) and (22) are convenient since they can be employed directly in the expression (17) for the rigid body contribution to the model weak form. As a result, (17) can be written at configuration i + 1 as a function of the incremental displacements u Δ P and incremental rotations α Δ P . The time integration is then already embedded. All other quantities involved, such as b i+1 and J i+1 are functions of the incremental displacements and rotation parameters. The corresponding values from the previous configuration are considered to be constant and known.

External loads
The contribution of the external loads to the virtual work δW e in (1) stems from each rigid body B, given by δW B e .
This term includes the virtual work of forces and moments where f B e is the external force applied at pole P, m B e is the external moment applied at pole P. When using the same variables as in (18) δW B e is given by the relation which introduces the external pseudo-moment μ B e = Ξ T m B e that is energetically-conjugated with δα Δ P . The pseudomoment depends on the amount of rotation experienced by the pole P, which is embedded in Ξ according to (9).
In the context of DEM, the weight of the particles is usually an important contribution of the external load. This is included in our model according to [43]. Weight in a particle can be represented by a single pointwise load applied at the center of mass. Since the pole P is arbitrary, one has to consider an equivalent system composed of a force and a moment applied at P. Considering a body with gravitational mass m (the same as inertial mass) and the vector g e describing the gravitational field yields

Contact between bodies
The treatment of the interaction of many particles is essential for a successful DEM implementation. Complex strategies have to be developed, especially for the handling of contact between general polyhedra that are possibly non-convex. The main idea is not to assume a single pointwise contact interaction for each pair of bodies AB (terms δW AB c in (2)) but instead to allow all pointwise contributions composing the term δW AB c to interact at the same time.
In the sequence, a master-master contact formulation is presented to handle single contact pointwise action-reaction. Next, a strategy to split the external surface of each particle in sub-regions is proposed, this allows to capture multiple contacts.

Master-master contact formulation
The master-master contact formulation is presented here for the particular case of contact between two rigid polyhedra. The reader interested in more details in the basics of this method may refer to [50,51], in the context of beam-to-beam contact.
For a triangular face of a rigid polyhedron at configuration i + 1 the following surface parameterization is proposed where the quantities x i+1 P and Q i+1 relate to the rigid body kinematics, described in Sect. 2.1. The vertices of the triangular face are denoted by the position vectors x r 1 , x r 2 and x r 3 at reference configuration where the origin is the pole P. Arbitrary material points are described on the triangular surface by the functions N 1 , N 2 and N 3 where the parameters ζ and θ are chosen from the range (ζ, θ ) ∈ (−1, +1). They define a mapping from the parametric space ( Fig. 1) to an actual material point in the three-dimensional Euclidean space (Fig. 2). Employing this parameterization for both contacting bodies, their external surfaces (faces) are represented locally. The polyhedra faces are given by Γ A (ζ A , θ A ) and Γ B (ζ B , θ B ). The parameters (which can be interpreted as convective coordinates) can be organized in vectors c A = ζ A θ A and c B = ζ B θ B . The surfaces depend on the general DOF organized in vectors, here named as d A and d B which leads to can be introduced. We assume that a single pointwise contact interaction occurs locally between the surfaces Γ A and Γ B .
In the original master-master scheme, there is no preselection of a material point where contact takes place within the surfaces. Instead a Local Contact Problem (LCP) is defined to find the material points in both surfaces where contact occurs. Initially it is considered that both surfaces are locally smooth. The LCP is defined for a fixed set of DOFd. Then, one has to find the set of convective coordinatesc = ζ AθAζBθB T such that the conditions are tangent vectors of the surfaces. The conditions in (30) are interpreted as orthogonality relations, meaning that the gap vector is orthogonal to both surfaces, see Fig. 3. The solution of the LCP yields the point of contact interaction between the surfaces. For the definition of interface forces at the contact point it is useful to define the contact normal n as a quantity which depends on the gap vector g So far the strategy is restricted to smooth surfaces because (30) uses information from the tangent directions of the surfaces. In case of singularities (such as the edges and vertices of polyhedra) one cannot define tangent directions uniquely. Instead, it is necessary to employ the sub derivative concept.
A strategy to handle singularities using the LCP degeneration was proposed in [32] and [33]. The basic idea is to alleviate selected orthogonality relations from (30). In a systematic way [32] defined a degenerative operator P s where c s is the vector with degenerated parameters, which dimension is s ∈ N|1 ≤ s ≤ 4. With adequate choices of the degenerative operator, the requirement of orthogonality in selected directions can be alleviated, thus creating an adequate geometrical treatment for contacts involving surfaces with singularities, such as the faces of a polyhedron. The LCP can be re-written after its degeneration as, see details in [32],

LCP degenerations
Discrete element models with polyhedra have complex contact interactions. These are vertex-to-face, edge-to-face, faceto-face, edge-to-edge, vertex-to-edge or vertex-to-vertex interactions. The special case of a face-to-face contact interaction deserves more discussion. As the faces are planar by definition, the only possibility for a face-to-face interaction needs parallel faces with anti parallel external normal directions n A and n B , see Fig. 4a. This would yield a continuous contact in the intersecting areas, see Fig. 4b. The new idea is, to approximate this continuous contact by a set of pointwise contact interactions involving the singularities of the faces, as shown in Fig. 4c for the case of edge-to-edge and vertex-to-face interactions. In the following we will use the representation depicted in Fig. 4c. A similar strategy is also applied for the edge-to-face interaction, not illustrated here, but with an equivalent pointwise contact force representation, involving vertex-to-face or edge-to-edge interactions (see example 1 in Sect. 3.1). This special treatment reflects the fact that a perfect parallel approaching of faces is not only improbable, but any disturbance on parallelism actually would lead to a description of pointwise contact(s) involving neighboring singularities (edges or vertices), as shown in Fig. 4c. Therefore, in the DEM context this seems to be the best choice.
Next the different kinds of degeneration are discussed. The analytical solution of each LCP degeneration can be found in the "Appendix".

Vertex-to-face
In this case, the material point assumed as contact-candidate (a vertex) is known in one surface. As an example: if the vertex is x r 1 in Γ A in the reference configuration, the parameters ζ A = −1 and θ A = −1 are already known and fixed. Contrary, in surface Γ B , no a priori values of parameters are known. This solution can be systematically treated by the master-master formulation with the special degenerative operator The geometric interpretation of this degenerated LCP is the solution of the orthogonal projection of a vertex onto a face. In this case, the contact normal direction is parallel to the external normal direction of the face.

Edge-to-edge
In this case no material points are assumed as contact candidate a priori. However, the edges are associated with particular values of the parameters. For example, in a surface Γ A the edge connecting vertices x r 1 and x r 2 can be considered in reference configuration. With that θ A = −1 is known a priori. In surface Γ B , considering also the edge connecting the corresponding vertices x r 1 and x r 2 leads to the already known value of θ B = −1. Based on that knowledge the LCP can be solved by employing the degenerative operator The geometric interpretation refers to the solution of the minimum distance between two lines in space. The contact normal direction is not related to any face external normal in this case. Instead, it is determined by the LCP solution, disregarding some orthogonality relations from the original LCP in (30) by using the proper degenerative operator defined in (35).

Vertex-to-edge
In this case one vertex is a material point candidate for contact in a surface, while in the edge it is already chosen. As an example, if the vertex under analysis is in Γ A and its coordinate is given by x r 1 in the reference configuration, the parameters ζ A = −1 and θ A = −1 are known and fixed. If the edge of interest is the one connecting x r 1 and x r 2 in the reference configuration on surface Γ B this leads to the known value of θ B = −1. Therefore, the LCP is a singlevariable problem, in which the only parameter still not chosen is ζ B . The treatment by master-master contact formulation is systematized in this case using the following degenerative operator: The geometric interpretation of this case is the orthogonal projection of a point (the vertex) onto a line in space (containing the edge). As in previous case, the contact normal direction is not related to the external normals of the faces. It is determined by the LCP solution using the proper degenerative operator defined in (36).

Vertex-to-vertex
In this case both vertices are already taken as the material points candidate to contact. Hence the vertex x r 1 in Γ A relates to the parameters ζ A = −1 and θ A = −1 which are known. Analogously, for the vertex in Γ B denoted by x r 1 in the reference configuration, the parameters ζ B = −1 and θ B = −1 are known. Hence no LCP has to be solved for this point-topoint contact (similar to node-node FEM descriptions). The contact normal direction follows directly from the direction of the vector connecting both vertices.

Remark
The discussed examples of degenerations are given for particular cases of vertices (e.g.: x r 1 ) and edges of surfaces (e.g.: connecting the points x r 1 and x r 2 ). However, all vertices and edges may be described by such particular choices, just changing the numbering sequence for the selected vertices within a given parameterization.

Surface splitting in sub-regions
To make practical use of the results in the Sect. 2.4.2, a strategy can be established in which the external surface of each body is split into sub-regions. Each of these sub-regions will be associated with a local parameterization Γ A in body A and Γ B in body B. After this split, pairs of Γ A and Γ B are investigated, in order to address possible contacts. The sub-regions are denoted by A i and are part of surface of the body A, with The idea is illustrated in Fig. 5, where two polygons are divided in sub-regions. The left one has sub-regions A 1 to A 5 associated with edges and A 6 to A 10 associated with vertices. The right one has sub-regions B 1 to B 4 associated with edges and B 5 to B 8 associated with vertices.
In the context of DEM with polyhedric bodies, the subregions are: (i) faces; (ii) edges and (iii) vertices. These have to be defined for each polyhedron. The contact check involves investigation between all sub-regions of a polyhedron against all sub-regions of another one. This strategy is quite general and works for arbitrarily shaped polyhedra, possibly non-convex. However, the associated computational cost is high because many sub-problems have to be solved. Figure 6 depicts a non-convex polyhedron, with sub-regions of faces, edges and vertices. Even with its non-convexity, it can be addressed by the proposed strategy.
A sub-region does not only consist of polyhedron faces, but also of edges and vertices. When evaluating the LCP for a pair of sub-regions the associated degeneration scheme has to be selected. Note however, not all pairs of sub-regions have a LCP to be solved. As already pointed out, pairs of face-to-face and edge-to-face are not considered.

Valid solutions of LCP
After obtaining the solution of a given LCP (see "Appendix"), one has to check if the set of parametersc lies within the valid range of the parameterization of both surfaces Γ A and Γ B (see the valid parametric space in Fig. 1). For parameters outside the valid range a contact interaction is not further considered for this sub-region.
As an example, consider the approach of a vertex in Γ B , named a sub-region B 1 towards the body A. Body A has sub-regions defined on its external surface. Some of them are depicted in Fig. 7. In this case, depending on the position of B 1 , its orthogonal projection onto the faces A 1 and A 2 may result in a solution which is outside the valid range. However, the projection of B 1 onto the edge A 3 leads to a valid contact-candidate, with its orthogonal projection falling in a valid range. This interaction describes a vertex-to-edge case.
Another example relates to the approach of two pyramids A and B shown in Fig. 8. In this particular case, the approach induces a contact interaction of their apexes. Depending on the trajectories of sub-regions A 1 and B 1 , the vertex-to-vertex case may lead to the only valid solution.
Testing the four possibilities of degenerations as presented in Sect. 2.4.2 has two drawbacks. The first one is related to the computational cost, because one needs to define sub-regions for all faces, edges and vertices of all polyhedra considered in the DEM model. This leads to a high number of geometric entities to be tested. The second drawback relates to the B1 A1 possibility of simultaneous validity of more than one LCP considered for a single pointwise solution. A simple example for such scenario is shown in Fig. 9 which exhibits a vertex B 1 , as a sub-region of a body B approaching a concave region of a body A. The sub-regions of body A in this case are the faces A 1 and A 2 and the edge A 3 . The LCP, interpreted as the orthogonal projection of the vertex B 1 onto the other geometric entity on A in this case yields valid parameters for all three tests. Therefore, in case of contact occurrence, more than one pointwise action-reaction check is needed for B 1 . This avoids penetration properly.
Remark A contact hierarchy can be introduced to circumvent such difficulties. Strategy rules can be created to favor some contact interactions against others, depending on the expected geometric case. A simple (and preliminary) discussion of such idea is employed in the sixth numerical example, but the topic is worth of further discussions.
Even with the mentioned drawbacks we see our strategy as promising because it is able to detect any kind of contact interaction between convex and non-convex polyhedra, including a natural description of multiple pointwise contacts in one polyhedron, which is essential for DEM to handle multiple interactions due to local concavities. In order to make this approach practical and feasible, a robust strategy for detecting possible contact pairs of sub-regions on bodies A and B are necessary. These ideas are expanded in Sect. 2.5.

Contact detection
Once the LCP is solved, one needs to check for contact and to enforce the associated contact constraint. Many approaches are available for the latter, such as the penalty method, Lagrange multiplier approaches, the augmented Lagrangian method, among others, see, e.g.: [36]. Usually one needs to detect the penetration between bodies in order to activate a contact constraint in the model. A scalar gap quantity with a sign is used for that, as in classical master-slave schemes.
In the master-master contact formulation, the kinematic quantity available to detect contact is the gap g. Its norm can be evaluated and defines a scalar g n = g . This quantity is always positive and cannot solely be used to detect penetration between bodies, but defines only their proximity level. However a simple dot product test of the contact normal at configurations i and i + 1 can be used, as shown in [32]. This leads to the evaluation of n i · n i+1 . If the result is negative, an inversion in the direction of the gap has taken place between configurations i and i + 1 which indicates a penetration.
This strategy has been used in [32] for contact detection when employing the master-master formulation. However, in the present context it is not sufficient. A contact pair could slide from one sub-region pair to a neighbor sub-region. In this case, the normals at time t i and t i+1 belong to different sub-regions and thus the above evaluation does not make any sense. Similar problems occur for thin bodies or vertices at cone tips and other cases. An enormous complexity is related to a robust contact search in the context of general non-convex polyhedra. Especially in DEM where thousands of different cases have to be resolved at every time step. To overcome these difficulties the barrier method is applied, which never permits penetration, but fulfils the contact constraints in a weak sense.
We present next the contact contributions to be included in weak-form. The interface laws described in Sect. 2.4.7, provide details on the strategy for an implementation of the barrier method.

Contact contributions to the weak-form
After establishing sub-regions of each body as introduced in Sect.

one can write
where δW AB i j are the contributions from each pair of subregions i j of a pair of bodies AB to the weak form.
Next, the terms δW AB i j are addressed based on the mastermaster contact formulation with its degenerations, including normal and frictional contributions. The presentation is short, more details can be found in [32] and [33].
The contact contribution δW AB i j can be split into are discussed. They refer to the contact between bodies A and B with sub-regions i and j. For the ease of notation, these indices are omitted in the following and only used when necessary.
The elastic term related to the normal direction yields where the elastic contact force is given by f n = − f n n i+1 and the gap in normal direction by δg n = δ g i+1 = δ g i+1 · n i+1 . The negative sign reflects the way of contact treatment with the barrier method. Note, the smaller g n , the larger is f n . Details are discussed in Sect. 2.4.7.
Damping in normal direction is given by where the damping force is given by f d at the contact point. At this point it is convenient to introduce the relative velocity in normal direction for the contact pair at configuration i + 1, which is derived in [33] and leads tȯ This kinematical quantity has to be employed for the evaluation of f d in Sect. 2.4.7. It is also possible to define a vector quantity for the relative velocitẏ Friction contributions are related to tangential kinematics of the master-master contact which was developed in [33]. The quantity g i+1 t , resembles the tangential gap at configuration i + 1. This term quantifies the sliding amount in tangential direction. It is updated when contact takes place between configurations i and i + 1. The increment of the tangential gap is given by With this increment the accumulated tangential gap at configuration i + 1 can be computed where Q Δ c is a rotation tensor to account for rigid body rotations experienced by the contact normal/tangent from configuration i to i + 1.
Note that in present contact formulation the interpretation of the tangential direction of a contact pair needs no definition of a tangential direction of contacting surfaces. Instead, it is defined by a projection to the tangential direction using the normal n i+1 in (43). Therefore, it can handle singularities such as edge-to-edge contact or other particular cases in a natural and straightforward way, which is fundamental for the robust algorithmic treatment within DEM.
With that, the friction contribution to the weak form yields where f t is the friction force. The variation of the tangential gap at configuration i + 1 can be evaluated in a similar way as the time-derivative of such quantity at configuration i + 1.
Following [33] we obtaiṅ The operator D can be found in [32] and denotes the relation between δc and δd, which depends on the solved LCP. Therefore, contact degeneration plays a role within evaluation of D, see [32], The degenerative operator P s has to be selected according to the degeneration, see Sect. 2.4.2, and the derivatives r 4,c and r 4,d are obtained differentiating Eq. (32).

Remark
The consistent linearization of the weak form involves relevant geometric terms. Indeed, one can evaluate δ g i+1 consistently employing the operator D as follows: With that, the LCP constraint information is embedded in the solution of the global model. This is essential for achieving quadratic convergence in the Newton-Raphson Method, as discussed in Sect. 2.6.

Interface law in normal direction
Contact constraints can be exactly enforced by techniques such as the Lagrange multiplier method, with the drawback of extra unknowns in the model, the Lagrangian parameters.
As an alternative, a popular choice is the penalty method, which permits penetration in a controlled manner, ruled by a penalty stiffness parameter. In this case, no extra unknowns are introduced. Instead, numerical problems like bad conditioning of equations systems can occur, which are related to the contact stiffness. Both techniques need a computation of the penetration within each contact region. As discussed in Sect. 2.4.5, the complexity to detect penetration is very high for polyhedral particles. The main difficulty stems from the need to check for pointwise contact occurrence between sub-regions. This local search needs to be performed for many pairs of sub-regions, leading generally to multiple interactions, possibly involving concave contact regions when evaluating the gap g. A different solution to treat contact problems relies on a technique, the barrier method, that does not permit penetration. Instead, it detects contact based on proximity between bodies. This method is frequently employed by Computer Graphics community, see, e.g.: [52] and [53].
The basic idea of the barrier method is to enforce the contact constraint by introducing a gap-dependent barrier function in the model potential. This function assumes very large values when two bodies (proximity) approach each other closely. The first variation of the barrier potential leads to a term in the weak-form, seen in (39). It yields an expression for elastic force f n = − f n n i+1 in normal direction.
The magnitude of f n follows from an evaluation of the barrier function, which is always active, or active only when two bodies come close. Similar treatments exist from the physical point of view when one is interested in very small geometric scales. In this context, in molecular dynamics forces stem from potentials, such as the Lennard-Johnes, potential, see, e.g.: [54]. The main advantage is that no penetration occurrence is needed to activate a contact constraint but only proximity of two or more bodies.
Based on such ideas we propose a compliance law for the contact behavior in normal direction. The law can translate a local physical behavior of contact (such as e.g. Hertz contact [55]) or can be interpreted as a numerical way to enforce contact constraints approximately-such as in a penalty-based approach.
The interface law will provide nonzero contact forces only within the proximity level g n <ḡ n . With that, one avoids the possibility of spurious contact forces exerted between bodies located arbitrarily distant. The valueḡ n is an activation threshold. Moreover, a very large force stems from the barrier method when g n → 0. But we need to keep a physical/numerical control when the interface law is activated for not too close proximity levels. This can be achieved by a hybrid interface law composed of two parts: (i) a physically or numerically ruled part for not too close proximity values (from the activation threshold value until a certain fraction of it) and (ii) a barrier-based part when the bodies come very close, to avoid penetration where 1 and n 1 are parameters to establish the part (i) of the compliance law,ḡ n is a proximity threshold to change from part (i) to part (ii) of the compliance law. We definē g n = fḡ n , where f ∈ R|0 < f < 1. The parameter n 2 is associated with the intensity of the barrier function. The remaining parameters 2 and c 2 are chosen to ensure the matching conditions resulting in Within the law in (50) we achieve a smooth transition between parts (i) and (ii) of the interface law, which yields good numerical behavior. One can choose parameters for the part (i) to recover a Hertzian compliance law, which is the basis of many DEM implementations for spheres see, e.g. [44,56,57], and for super ellipsoids see, e.g. [20,21]. One can also recover in part (i) a simple penalty-based linear law. Examples in Fig. 10 illustrate such interface laws.
When employing the proposed interface law, a threshold valueḡ n for activation of the contact parameter has to be selected. This can be visualized as a thickness of a contact skin to be superimposed on all surfaces of the bodies. For smaller value ofḡ n less alteration is introduced in the geometry for contact detection/evaluation. The value ofḡ n is also linked to the desired contact stiffness (ruled by 1 ). The choice of smallerḡ n is usually combined with larger 1 . Otherwise, only the barrier part of the interface law will be experienced by the bodies, which does not take advantage of the hybrid aspect of the interface law. Thus, the choice of parameters has to be a compromise between the desired (usually small)ḡ n and the contact stiffness. This is similar to the relation between contact stiffness and allowed pene- trations, when referring to contact formulations that permit penetrations (for example the penalty method). Figure 10(b) shows an example of interface law that is adopted in example 6, Sect. 3.6. In order to estimate the maximum time-step to that is allowed in an explicit integration scheme, one evaluates the Courant criterion. By representing the contact pair of two particles as an equivalent mass-spring oscillator, the tangent stiffness can be obtained as k c = − d f n dg n and an equivalent mass is estimated by m c . A rough estimate for the critical time-step is then given by Δt c = 2 m c k c . Figure 11 shows that this critical time-step achieves small time increments when a close gap is experienced during con-tact. In present work we employ an implicit time-integration scheme, thus no critical time-step has to be respected and the computation can be run with much larger time steps, as demonstrated in the numerical examples.
When addressing contact with the barrier method, one cannot permit penetrations, which would not make physical sense. However, when solving the nonlinear system (1) by the Newton-Raphson Method, it is possible to find (wrong) solutions that include penetration between bodies. Then, one needs an efficient way to verify occurrence of penetration. This can be done at the end of each converged time-step, by executing the simple test n i · n i+1 . The result has always to be positive, for all contact pairs that are present in the model. In case of a negative value, one has to discard such solution and to use a smaller time-step to find a valid one.

Damping model for normal direction
Damping in normal direction is classically related to the coefficient of restitution, which is a measure of the energy dissipated when contact-impact has taken place. Here we follow the alternative way, which describes damping as viscous dissipation. This is common in many DEM applications, see, e.g.: [56-58] and [21]. Here we use the model from [21], but do not assume Hertzian contact. We create a damping ratio input for the model, leading to a instantaneous viscous damping coefficient (dashpot-model). Damping is evaluated in a consistent way with our proposition of interface law the elastic contact force in (49).
The viscous damping force is evaluated using where ζ is the desired damping ratio and d f n dg n depends on the interface law, see (49). The desired level of damping is a choice to match a known experimental data, as in the case of the coefficient of restitution, when modeling contactimpact by distinct approaches. The negative sign inside the square root is inserted because the required derivative of the normal elastic force with respect to the proximity g n is always negative. The parameters m A and m B are the masses of the contacting bodies.
The damping force f d in (40) is given by which ensures that the total normal contact force ( f n + f d ) is always compressive, as in [21]. The damping force in (52) represents a simple way to introduce a desired damping ratio ζ for an equivalent linear spring-dashpot oscillator at each instant.
Note that a given input of ζ does not necessarily represent the classical behavior of a mass-spring-dashpot model (for example ζ = 1 leading to a critical damping ratio). Indeed, the viscous damping can change when expressed by (53). Moreover, in the present context we may have multiple simultaneous contact forces between bodies. In this case, the employed equivalent mass evaluation could be improved, thus imagining that each contact pair would have its own damping related to the masses of the pair. We decided not to propose such a refinement in our damping model, to keep it simpler. However, this topic is worth an investigation, but out of the scope of present work. Thus, ζ may be seen as an equivalent input parameter to reach a desired level of damping, which is defined similarly to a classical mass-spring-dashpot model, but having distinct particularities.

Interface law in tangential direction
The classical Coulomb law is introduced to model friction in each contact point. The friction force follows from a rheological model for the tangential direction based on elastic and damping contributions. This comprises a model, already applied for sphere-to-sphere contact in [57]. We employ enhancements as proposed in [33].
For a given contact pair the tangential gap and the relative velocity can be evaluated and used to compute a trial friction force composed of an elastic part f e t and a viscous part where t is a tangential penalty stiffness and c t is a tangential viscous damping coefficient. The trial friction force is given by Both, f e t and f d t lie in the tangent plane of contact (I − n ⊗ n ), defined by normal direction n. However, forces f e t and f d t are not necessarily parallel. The magnitude of the friction force is limited by Coulomb's law. When the current contact status is sticking, no sliding occurs in next configuration if f tr t ≤ μ s f n , where μ s is the static friction coefficient. In this case, contact keeps in sticking status. Otherwise, a sticking to sliding transition takes place.
When the current contact status is sliding, sticking occurs in next configuration if f tr t ≤ μ d f n , where μ d is the dynamic friction coefficient. Then, a transition sliding to sticking takes place. Otherwise, sliding keeps active.
Next to the Coulomb inequality test the tangential force is updated. For a sticking contact status the update is f t = f tr t .
When sliding occurs, the friction force is given by is the sliding direction. Finally, one has to update the tangential gap g i+1 t for the next step when sliding occurs. This has to be done in such a way that this quantity holds the "recoverable part" of sliding tendency. Here we adopt the rheological model proposed in [33], which considers only the elastic part of friction for a possible update of g i+1 t . One has to check if f e t > μ f n , where μ = μ s if the contact status is sticking and μ = μ d , otherwise. If this inequality holds, sliding in the elastic part of friction is given by: Then, the following updating formula follows: If f e t ≤ μ f n no sliding takes place in the elastic part of friction model.
Note that the status of sliding or sticking for the contact is ruled solely by the inequality test employing f tr t . Thus the composed friction model with elastic and damping parts has to be used. The proposed update for sliding in (57) considering only f e t exhibits some numerical advantages discussed in [33].

Global contact-search
The computational bottleneck of numerical solutions using DEM is related to contact search and contact evaluation. The heavy computational cost stems from (2), which shows the double summation associated with all possible contact contributions. In our formulation the scenario is even more complex, because for each contribution δW AB c one needs to seek possible interactions between sub-regions of bodies, as expressed in (37). Therefore, the number of contact candidates is very high in practical problems involving general polyhedra.
To make the model computational feasible, the global contact-search has to be enhanced, in order to avoid evaluations of contact candidates that are far away from each other. Different global search strategies were developed in the last years, involving, e.g. bounding volume (BV) overlapping search and combined Verlet and Linked-Cell schemes.

Bounding volumes overlapping
A first level of BV is defined for bodies, considering only spheres as shown in Fig. 12a. This provides a quick search, for overlapping/proximity between such spheres and permits to eliminate non-strong candidates. Our strategy here is to always construct inflated BV's, ruled by an inflation factor. Then, in case of overlapping, proximity is detected. Otherwise, no contact is considered between the tested bodies.
Once overlapping is detected in the first-level BV search, a second level is entered, in which a BV is defined for each sub-region of the that overlapped in the first BV search. In this context, we consider oriented and inflated BVs for each geometric entity: (i) faces have a triangular prismatic BV, (ii) edges have a cylindrical BV and vertices have spherical BV. An inflation factor is introduced for each defined BV. Again a search for overlap is started. In the positive case, a pair of sub-regions becomes a strong (probable) contact-candidate. Otherwise, no contact is considered. Examples of BVs for geometric entities are shown in Fig. 12b-d.
Only the strong contact-candidates are checked for contact using the methods in Sect. 2.4 which results in contributions to the weak form (38).

Combined Verlet and linked-cell schemes
Both levels of BV overlapping search are computationally costly . Particularly, the second-level (sub-regions) includes elaborated BV. Due to numerous expected tests in large scale practical problems, one has to improve the contact-search even more.
Instead of checking all combinations of BV overlapping at each time increment, a list of probable contact-candidates can be established by a simple neighboring search. Only the candidates from such lists are worth of a closer investigation employing the inflated BV overlapping search. These lists are based on a Verlet scheme [54]. In the first level, they are generated considering the distance between pairs of body poles, as reference points, together with the radius of the inflated BV. In the second level, the centroid of each sub-region BV is taken as a reference, together with a pre-computed cut-off size. Proximity between sub-region BV is checked, therefore establishing a Verlet list also for the second search level. This process saves substantial computational time when compared to an all-to-all search.
A second scheme to save computational time while pursuing a global contact search is related to the establishment of a Linked-Cell scheme. The basic idea stems from molecular dynamics applications (see, e.g.: [4]). Spatial cells are defined as regions in space. Each body is associated with a cell, according to its current location (in our case, a reference point of the body is considered). Then, an enhanced search is performed considering only neighboring cells as candidates for a body located in a given cell. This leads to another list for enhancing spatial search, not based in each particle neighborhood, but as a spatial search to organize particles in groups, based on their current locations and associating them to cells.
Both Verlet and Linked-Cell schemes were implemented independently, and also in combined ways. The more efficient scheme is problem dependent. The reader may refer to [59] and [60] for some comparisons, as well as for some details on how to enhance the numerical implementation.

Global solution
A strategy for the overall model integration is derived on the basis of (1). We consider a system with N bodies. Each of the bodies has reference point (pole) P, with which six DOF are associated (incremental displacements u Δ P and incremental rotations α Δ P . To begin the time-integration, we have to assume an initial status for all possible contact interactions. No overlaps between bodies are considered in the initial (reference) configuration. This is essential for the usage of the contact normal interface law and the barrier method which cannot handle overlaps but proximity between contact surfaces. We start the simulation with an initial global contact search as discussed in Sect. 2.5. Considering only the strong contact-candidates, the LCP and contact contributions are evaluated. When g n <ḡ n , according to the normal interface law expressed in (49), an active contact pair is established. It has nonzero contributions to the weak form. All the nonzero contributions of (38) are inserted into (37) and thus all contact contributions are included in (1).

Set of equations for the global model
To perform the time-integration, we have to write (1) for a given time instant. Following the scheme of section 2.2 the still unknown (next) configuration i + 1 has to be considered. Thus, incremental DOFs appear together with current (known) values at configuration i. A 6N -dimension global vector of DOFs a is introduced where u Δ B and α Δ B (B = 1, ..., N ) are the incremental DOFs describing the motion of each rigid body. Analogously, it is possible to write a similar vector for the virtual quantities where δu Δ B and δα Δ B are virtual quantities associated with each DOF. Now the weak form (1) can be written as where the 6N -dimension vector of residuals r was introduced where each term r B f and r B μ (B = 1, ..., N ) represents a contribution given by The inertial contributions f B t and μ B t stem from (18) and the external loads contributions f B e and μ B e follow from (24). The terms f B c and μ B c represent contact contributions, stemming from each contact pair between sub-regions related to (38).
The weak form (61) is valid for arbitrary δa. This leads to a nonlinear system of 6N equations, with 6N unknowns. The solution is obtained by applying the Newton-Raphson Method. For that, one needs to obtain the consistent linearization of the residual (62), with respect to the unknowns. The solution of the nonlinear system of equations completes a single time-step evaluation within the Newmark method in Sect. 2.2.

Consistent Linearization
This process usually leads to cumbersome algebraic work. Alternatively (as we did here), one can employ automatic differentiation techniques, leading to automatic code generation for such complex expressions. Here the tool AceGen is used, see [61]. For the automatic differentiation, one needs to use a symbolic computational framework to define all the expressions necessary to evaluate the residual (62). All expressions must be defined as a function of the model unknowns (in our case, the displacements and rotations). With that, the AceGen tool is able to automatically evaluate the partial derivatives of all quantities with respect to the model unknowns (DOF), and perform their consistent linearization. The total consistent linearization is composed of contributions stemming from each particle, which are independent of each other, and of contact terms, which couple pairs of particles. Our strategy was to separate such effects both for residual evaluation and also for its consistent linearization. For each particle, one has to write the contributions f B t , f B e , μ B t , and μ B e . In such quantities, the time-derivative quantities have to be written following the Newmark adopted scheme, yielding only expressions dependent on the model DOF. Programming these expressions in the AceGen tool, one can evaluate the necessary linearizations. The contribution of each particle is included on a global residual and global tangent by a procedure similar to done in the FEM, when including local influences of an element into the global system of equations.
The consistent linearizations of contact terms follow a similar procedure. One has to write f B c and μ B c as a function of DOF and use the same tool to perform the consistent linearization. After, the contact contributions are included on a global residual and global tangent. As each contact contribution involves DOF of two particles, such contributions are the responsible for the overall system coupling between particles. There are efficient ways of programming the AceGen tool to provide directly the residual and its consistent linearization, using potentials (when available) or pseudo-potentials, instead of defining directly the residual and claiming its linearization. The reader can found more details of such techniques in [61].

Remark
The obtained system of nonlinear equations couples DOFs of different bodies due to contact contributions. As contact occurrences obey a physical neighborhood, the consistent linearization of r, represented by a tangent stiffness matrix, has a sparse structure. This permits saving memory and solution time while solving each Newton-Raphson iteration, within a time-step. Hence solution techniques employed when solving nonlinear FEM models can be applied. Hence it was possible to implement the DEM in the Giraffe platform that was developed for FEM, see [62], also aiming for future coupling of FEM and DEM models.

Night owl contacts
In the contact detection we speculate that non-overlapping BV at configuration i would not be responsible for active contact pairs at configuration i + 1. This has to be checked, however, by a second global search at the end of the timestep, after a converged solution has been achieved. All new overlapping BV found are tested to investigate if they represent active contacts, when g n <ḡ n , according to the normal interface law expressed in (49). If no active new contact is found, the speculation has been correct and one can advance to the next time-step, already using the just-obtained list of strong contact-candidates. Otherwise, one has to recompute the current time-step solution, considering additionally the new contacts in the model (we call these as "night owl contacts").
Obviously, night owl contacts are undesirable, since they require time consuming re-evaluation of the time-step. However, the speculation may be well assertive when considering large enough inflation factors for BV, which are related to the expected body velocity and assumed time-steps. Too large inflation factors are not efficient, since many strong contactcandidates are elected, leading to high computational cost for their evaluations with numerous zero contributions to the model. With that, one may expect an optimal computationalefficiency by choosing the smallest inflation factor, which leads to rare occurrence of night owl contacts. This may be estimated prior to simulation considering expected velocities to be experienced by the bodies which lead to distances that a body can travel within the next time step and employed to calculate the inflation of each BV. Similar ideas of speculation while looking for contact between bodies are found in [63] in the context of graphics computing.

Numerical examples
The discrete element model was implemented in the Giraffe platform, an in-house C++ code, initially conceived for FEM models, but extended to encompass DEM implementations [62]. A new modulus for automatic contact search with the above discussed strategies was implemented. Specific codes for evaluating contributions to the global problem as in Sect. 2.6 were implemented in Mathematica TM . AceGen [61] provided automatic differentiation, optimization and auto-matic code generation in C++ and then implemented in the Giraffe platform.
The considered examples rely on definitions of bodies and boundaries. The geometry description for the boundaries is similar to the bodies. Therefore, they are composed of triangular regions. The interaction body-to-boundary is treated the same way as in the case body-to-body.

Block falling on the ground
In this example a single cube is falling on the ground. No damping is considered in the impact occurrences, in order to prof conservation of mechanical energy.
The cube, has external normals of faces given in the local directions ±e 1 , ±e 3 and ±e 1 × e 3 . Each cube face is composed of two triangular regions, for contact detection/treatment purposes, employing the surface parameterization proposed in (28). The ground is flat. All numerical data considered are shown in Table 1.
Distinct initial orientations are proposed for the block, leading to the cases: (a) The block is aligned such that its center of mass is positioned exactly above a given vertex. Successive impacts between this vertex and the ground occur and no rotation is induced to the block. Results are shown in Fig. 13. (b) The block is aligned such that its center of mass is positioned exactly above the center of a given edge. Successive impacts involving pointwise forces on the vertices of such edge occur, representing the line-toface interaction by two pointwise forces. No rotation is induced to the block. Results are shown in Fig. 14. (c) The block is aligned such that its center of mass is positioned exactly above the center of a cube face. Successive impacts involving pointwise forces on the vertices of that face occur, representing the face-to-face interaction in a pointwise manner. No rotation is induced to the block. Results are shown in Fig. 15. (d) The block is aligned such that the contact force changes its angular momentum when it collides with the ground. Successive (bouncing) impacts occur involving distinct vertices and the ground, inducing a composition of translational and rotational movements to the block. Results are shown in Fig. 16.
For all cases one can observe conservation of mechanical energy, due to the absence of dissipative effects. To achieve that in simulations, an adequate time-step guideline (maximum) was chosen, enough to integrate the contact-impact interactions with reasonable precision. The contact parameterḡ n defines the so-called contact skin as an offset in the block. In this case, the chosen value represents 1% of the cube edge size.
Initial orientation (b) Initial orientation (c) Initial orientation (d)

Pile of blocks
In this example a set of blocks is analyzed, which fall in vertical direction under gravity load, forming a pile. A set of ten blocks is considered, each one modeled as in the previous example. The blocks are located such that they will touch on top of each other, after forming the pile. Their initial center of    Table 2. Significant energy dissipation is introduced to obtain a final state that is static. The contact properties are depicted in Table 2. A proper time-step was considered, in order to handle correctly the contact-impact scenarios. Figure 17 shows the final configuration of the formed pile. The red arrows demonstrate contact normal forces acting on the blocks. Their length is proportional to the force magnitude. Note that the face-to-face interactions are handled by the model of equivalent pointwise forces on the edges

Sliding of a block
This example proposes a scenario to test the friction model. A block (the same cube as considered in example 1) rests initially on a flat surface. This configuration is shown in 19a, which exhibits the normal force components evenly distributed on the four vertices touching the surface. Table 3 summarizes the model characteristics, including the contact data with a high friction coefficient. An external force in direction i is applied at the block, leading to a motion. The force magnitude is linearly ramped Normal force reaction F z between pile and ground while the pile of blocks is being formed up in time (from 0 to 0.8 s), leading to force values from 0 to 16 N. We consider three distinct application points for the force which are described by an extra point (node) in the model, defined as rigidly-connected to the block. This additional node increases the number of DOFs, but is constraint to the body motion. The associated equations are solved by the Lagrange multiplier method, see [39] and [43]. With this procedure it is possible to apply forces/moments at any point of the body and not only at the defined pole P. Based on this formulation the following cases are simulated: (a) Force applied at the bottom face centroid of the block: this case does not introduce a toppling tendency to the block. As long as the external force increases, the friction forces react, also increasing their magnitude. Friction is evenly distributed on the four vertices, see Fig. 19(b). This leads to the expected overall result for the total friction force applied on the block, as shown in Fig. 20. When the maximum tangential force (Coulomb limit) is reached, the friction force is constant and the block starts moving (dynamic friction takes place). (b) When the external force is applied at block center of mass toppling can occur. As long as the external force increases, the tangential forces due to stick increase in magnitude. Differently from the previous case, now the normal forces differ at the vertices, to fulfil static moment equilibrium. We define a leading edge of the block which is the base edge related to the positive direction of i and a trailing edge which is the opposite one. The leading edge vertices endure larger normal forces, while at the trailing edge the forces are smaller, see Fig. 19(c). Hence, the Coulomb limit for friction leads to different tangential reactions for vertices located on leading and trailing edges. When the maximum available friction is reached at the trailing edge, there is still stick friction at the leading edge, leading to a sudden friction force redistribution. However, as the block is still in stick mode, the trailing edge recovers its friction in next instants. This process leads to small fluctuations over time in the friction force. When both, trailing and leading edge vertices reach the static Coulomb friction limit, see Fig. 20, motion starts and the fluctuations no longer exist because the friction coefficient keeps the value of μ d for all pointwise interactions. Even with the toppling possibility, the available friction was not enough in this loading case to topple the block. (c) In the last case, the external force is applied at the block top face centroid. This case has a similar behavior as case (b), having again a tendency of toppling. However, when the force reaches a certain level toppling occurs, when the block is still sticking to the surface. Figure 19d shows the block configuration just when toppling starts (note the absence of contact forces at the trailing edge). During toppling, friction forces invert their direction and assume even zero values because the block looses contact with the ground for a while due to inertia effects at the end of the simulation.
The results in Fig. 20 demonstrate consistency between the simulated cases and the expected physical results. The repeatedly sudden drops of the friction level, followed by its recoveries in cases (b) and (c) are a result of the transition between stick and slip states. We did not assume a smoothing in the transition from μ s to μ d . Hence when the friction force reaches the Coulomb limit, there is a non-smooth hang in the friction coefficient dropping suddenly, independently of the occurrence of a motion. This behaviour can be improved in future works by introducing a smooth transition from μ s to μ d based on the relative velocity at the contact interface which is common in the context of FEM contact models. When solving such example considering the same value for the static and dynamic friction coefficients, such oscillations in friction force no longer appear, as depicted in Fig. 21.

Tetrapod body experiencing a twist
A body with tetrapod shape is an example of a non-convex polyhedron. The data describing the body are summarized in Table 4. Furthermore, its stereolithography (stl) CAD file is provided as a supplementary material. We consider that the tetrapod initially rests on the ground under gravity loading, as depicted in Fig. 22. Due to the body-ground interaction vertical reaction forces are present in each touching point, having the same magnitude. These forces balance self-weight of the tetrapod. As an additional external load, we consider a moment applied at the body center of mass, trying to introduce a rotation in direction +k (twist). This external moment is aligned in direction +k, and its magnitude is linearly ramped during the simulation duration of 2 s, in the range from 0 to 750 N.m.
Contact between the tetrapod body and the ground is distributed on 3 vertices, as shown in Fig. 22. When applying the external moment load, the tetrapod body keeps in static equilibrium, due to friction occurrence, as indicated in Fig. 22 by red arrows. As the external moment increases its magnitude, friction also increases, up to the Coulomb limit, when finally the system exhibits a transition from statics to dynamics, and the tetrapod body starts to twist. Fig. 23 shows the reactive moment induced by friction forces, which balances the external moment up to the instant approximately 1.6 s, when the Coulomb limit is achieved. From this instant on, the induced reactive moment keeps constant. Note that in this case we did not consider a difference between static and dynamic friction coefficients, which would result in a drop of such reactive friction moment, when twisting starts.

Box of blocks
The interaction between several bodies are discussed in this example. First, we establish an initial positioning for a set of 564 identical particles. Each one has the same geometry as the cube proposed in the first example. The cubes are initially positioned without overlapping, having arbitrary orientations, as shown in Fig. 24a. We define four planar bounding surfaces, depicted also in Fig. 24 which form a   box. The bottom of the box is a square with side length 1.0 m. Table 5 summarizes environment and solution data, such as contact data considered for body-body and body-boundary interactions.
From the initial state the particles move inside the box under gravity, forming the pack shown in Fig. 24b.
A second phase of the simulation follows by removing the lateral walls of the box inducing the pack to collapse. The particles are spread on the ground. Fig. 25 depicts a sequence of selected configurations assumed by the model  as long as the time advances, until the final configuration has been achieved (static). Note that some blocks close to the middle of the pack still form a small pile at the end of the simulation which is related to friction.
We employed here an automatic time-stepping, which adapts the time-step value according to difficulties for convergence encountered during the solution. This adaptive scheme is ruled by a guideline of time-step range.
Active contact pairs were monitored along the evolution of the simulation in time. Fig. 26 demonstrates for the second phase (pack collapse) of the simulation the number of active body-to-body and body-to-boundary contacts. All contacts are handled automatically by the implemented solver, considering the contact degenerations proposed in Sect. 2.4.2.

Box of tetrapod-shape particles
This example considers the collapse of a pack of tetrapodshape particles (the same as presented in example 4). Handling particles of this shape within a DEM simulation is quite complex, since the tetrapod has thin tips and concave regions. Thin tips are difficult to deal with in contact algorithms that are based on the computation of penetrations, because spurious contact detection can occur. Our scheme takes advantage of the barrier-based normal interface law, making this treatment much simpler and leads to a proper scheme for handling contact of thin tips. Moreover, concavities are challenging due to difficulties to deal with multiple pointwise contacts, as discussed in the scheme shown in Fig. 9. Therefore, we see the tetrapod-shape particles as an interesting test for the new approach to the treatment of contact discussed in this paper.
To establish a initial pack of tetrapod bodies, the following steps are performed (a) establishment of an initial configuration with a set of 1005 bodies with no overlapping and assuming arbitrary orientation, similar to example 5, see Fig. 24a; (b) dropping of all bodies in the box, also similar to the procedure shown in Fig 24b. This process leads to a first pack of bodies. However, in this case we would like to obtain a more compact packing, thus additional steps are followed; (c) moving of lateral walls of the box towards the center of the pack, thus forcing the bodies to adjust their position and increasing the height of the pack. The final shape of the box bottom is a square with side length 7.0 m; (d) turning off the friction forces, as an artificial feature for compaction. This forces the particles to reorganize, filling many voids, due to absence of friction; (e) turning on the friction force again and letting the pack come to rest as shown in Fig. 27a. This represents the initial configuration for the simulation of the pack collapse.
For the evaluation of the pack collapse, we removed suddenly the lateral walls of the filled box shown in Fig. 27(a). We performed a simulation of 7.0 s duration considering only gravitational field as source of external loads (in direction −k) with magnitude 9.81 m/s 2 . The boundary is defined only by the ground, considered as a flat surface. The contact model for all interactions has the same properties considered in example 4, as shown in Table 4. Figure 27b depicts an intermediate configuration after 3.0 s of the simulation and Fig. 27c shows the final configuration obtained, after 7.0 s. The same final configuration may be seen in a lateral view in Fig. 28, which shows the interesting characteristic of interaction of tetrapods. Many interlockings occur due to their particular shape. Figure 29 depicts the number of active contact pairs along time, both considering body-to-body and body-to-boundary interactions.
In examples 5 and 6, we used an automatic time-stepping solver to guide the solution. The solver adjusts automatically  the time-step size according to success or not while trying to find Newton-Raphson convergence. When contact-impact takes place, convergence difficulties are faced and the solver has to decrease the time-step for successful integration. In systems with high dissipation level such as examples 5 and 6, one can expect an appropriate solution to be obtained with this time-stepping strategy, even without a strict control of energy dissipation in each single impact occurrence (as we did in example 1). This approach is only possible when solving the dynamics equations with an implicit solver, which has no upper-bound for the time-step. To set the solver time-step guidelines, however, one needs to use information based on the physics of the problem, as a reference for the time-step range. Figure 30 shows time evolution of the kinetic energy of the system for example 6, demonstrating the expected result during the pack collapse: an initial increase of the system kinetic energy (during the pack dismounting), followed by the system kinetic energy drop to zero, as long as the particles find their new static position on the ground or above other particles, forming a pile. Figure 31 shows the time-steps employed during the implicit solution. One can see clearly that smaller time-steps were required only between 1.0s and 2.0 s, which coincides with the peak in the kinetic energy. During this simulation contact interactions experienced a very small gap, achieving 0.001ḡ n or even smaller values, leading to a very high contact stiffness in (49). An explicit integration scheme would have needed a time step smaller than Δt c = 10 −5 s, as can be seen in Fig. 11. Therefore, using data from Fig. 31, the average time-step of Δt = 8.2 · 10 −4 s is 82 times higher than Δt c . With that, even with the higher computational cost of each time-step integration of an implicit scheme, the practical used time-step value leads to a brake even in computing time and the implicit scheme has advantages of more control on energy and yields very large time steps when the system is close to a static response. As a reference of computational solving time, the step (e) of this example (pack collapse) takes about 4 hours to complete in a CPU Intel Xeon W-2135 3.7 GHz with 6-cores.
Due to the complex shape of the tetrapods, one may face much more difficulties in this example for achieving convergence within each time-step. The main challenge was related to the description of the contacts, particularly in concave regions. Scenarios such as the one shown in Fig. 9 are frequently found and can lead to co-existence of multiple pointwise interactions in concave regions. Moreover, the search for equilibrium in Newton-Raphson iterations may lead to alternating patterns of switching on/off some contact pairs, due to entering/exiting their range of validity of projection (LCP solution). We succeeded in circumventing such problem in the presented example by specific contact hierarchy rules, constructed to foresee some problematic scenarios.   As an example, one may inactivate a vertex-to-edge contact in particular cases, when facing concave edges. When contact has already been established between a given vertex and both neighboring faces of a concave edge, one may inactivate a possible vertex-to-edge contact, such as in the example of Fig. 9. Similar ideas were also implemented for vertexto-vertex interactions. One may inactivate a vertex-to-vertex contact when some of the attached faces or edges to a vertex have already presented an active contact region with the another vertex.
This kind of geometric-based observations, thus transformed into a set of rules to switch off automatically some detected contacts, improves substantially the convergence of the Newton-Raphson iteration within each time-step, because it avoids alternating patterns of switching on/off in contact regions. Here we see need for further research and improvements of the method, especially an in-depth investigation of the above mentioned rules, which is a pure geometrical matter.

Conclusions
We propose in present work a discrete element method considering particles as general rigid polyhedra. Our main novelty is related to the contact treatment, where constraints are set up for multiple pointwise interactions between particles. We presented a strategy of splitting the external surface of each particle into a collection of sub-regions associated with geometric entities: vertices, edges and faces. Combinations of sub-regions form contact pairs. Each one is treated considering a distinct degeneration of a basic surface-tosurface formulation. The strategy is able to handle contact involving bodies with complex and general shapes, including thin tips and concave regions, as demonstrated in the numerical examples.
We propose a novel interface law for the normal contact treatment by a composition of a penalty-based and a barrier-based law, with the advantage of avoiding penetration between particles, but instead allowing a controlled level of proximity. For frictional cases, we employed a Coulombbased treatment, considering distinct static/dynamic friction coefficients.
The proposed formulation was implemented considering an automatic contact detection scheme based on inflated bounding volumes overlapping, associated with each considered geometric entity. Each particle has a bounding volume and a set of sub-bounding-volumes, which are assumed to be contact-candidates. Only strong-candidates are investigated by the introduced Local Contact Problem. This leads to a computationally feasible scheme, even when numerous necessary contact and proximity tests are involved.
We expect high potential for application of the present method to practical engineering problems, especially, in which the particle shape is important.
Note that these parameterizations are compatible with such edges as part of surfaces described by (28). The associated set of orthogonality relations are which permit the straightforward evaluation of the solution of the LCP given byζ A andζ B . One can employ this solution for any edge of face A and any edge on face B, just adapting a list with the sequence of vertices.

Vertex-to-edge
Here we consider a vertex in face A named x A , projected onto an edge in face B. One may find this analytical solution by proposing for face B the alternative parameterization for the edge γ B , as already pointed out in (67). In this case the only orthogonality relation to be fulfilled is: The same ideas are applied when projecting a vertex from face B onto a an edge on face A, just adapting the presented equations, accordingly. Moreover, one may employ this solution for any edge of a face, just adapting a list with the sequence of vertices.