Rigid Blocks

Introduction

The rigid block module introduces the capability of modeling closed, convex and manifold {polygons in 2D; polyhedra in 3D} to PFC. The facets of the rigid blocks are {linear segments in 2D; triangles in 3D}. As with clumps, the full inertia tensor is used to update the rotational equations of motion (see the description here). The inertial properties of each rigid block are computed as presented here. Facet connectivity information is retained to assist in contact detection and resolution. Rigid blocks exist within the model domain and can interact with balls, clumps, walls, zones (via the wall-zone coupling) and structural elements in 3D (via the wall-structure coupling and also directly). They can be created from vertices, imported from geometry objects, cut, and they also support edge/vertex rounding.

In PFC parlance, a rigid block is a piece. Pieces are objects between which a single contact can exist. Bodies, as described here, are composed of one or more pieces. In other words, only one contact can exist between a rigid block and other pieces, including balls, clump pebbles, wall facets and other rigid blocks. This is in contrast to clumps where each pebble of a clump is considered a piece, often resulting in multiple contacts between a multi-pebble clump and another piece. Such a reduction in the number of contacts can be beneficial for performance as runtime profiling suggests that L3 cache misses due to looping through contacts quickly dominates runtime degradation as the number of CPU cores increases. On the other hand, the existence of only one contact between a rigid block and another piece means that robust and perhaps computationally intensive contact detection and resolution strategies must be employed. The sections below outline how contacts between rigid blocks and other pieces are identified and resolved.

Contact Detection and Resolution

Contact detection is the process by which one can determine if pieces collide. Contact resolution, on the other hand, is used to define the contact normal, gap and position. This section focuses on describing the broad- and narrow-phases of contact detection where at least one contacting piece is a rigid block.

Broad-Phase Contact Detection

As discussed here, rigid blocks are entered into a cell space for rapid broad-phase contact detection. Separate cell spaces exist for each piece type. This construct allows one to efficiently cull the number of narrow-phase operations, making DEM simulations less computationally intensive. In order to expedite the process, each rigid block is abstracted as an axis aligned bounding box, termed an extent, and entered into cells (i.e., spatial bins) in the cell space based on an enlarged extent (termed the cell extent). Rigid blocks are candidates for collision with other pieces provided they map into the cells covering the same spatial regions in the cell spaces of other pieces. Once the cell extents of two pieces overlap, the pair of pieces is considered a contacting candidate; the narrow-phase contact detection algorithm is applied to such candidates. By limiting the amount peices can displace during each timestep, one can ensure that all contacts are detected prior to overlap.

Narrow-Phase Contact Detection

Once a pair of pieces is identified as potentially colliding, a contact is created. Contact creation triggers a narrow-phase contact detection algorithm depending on the piece types in question. As long as the contact exists, the narrow-phase contact detection algorithm runs each timestep prior to the force-displacement law. If the pieces are deemed to be in contact, an additional contact-resolution algorithm may be required to determine the contact gap, normal and position.

Rigid Block Contacts with Balls/Pebbles

For the case of contacts with balls or pebbles, an approach inspired by [Nezami2006] is used for narrow-phase contact detection. This approach seeks to determine the closest point on the rigid block surface to the ball or pebble centroid. Provided that the centroid falls outside of the rigid block, there is only one such closest point. The closest point is sought by using the facet connectivity information to traverse the rigid block surface, finding the surface point with global minimum distance from the ball or pebble centroid. During each timestep, all adjacent facets must be checked to see if the closest point on the surface has moved, and the surface can be locally traversed yet again to find the closest point if necessary. If the closest point on the surface is less than the ball radius, assuming the interaction distance is set to 0, the contact is deemed active.

Rigid Block Contacts with Rigid Blocks

Efficient narrow-phase contact detection between convex {polygons in 2D; polyhedra in 3D} can be very computationally intensive. The discussion below applies in both 2D and 3D. The Gilbert–Johnson–Keerthi (GJK) algorithm ([Gilbert1988]) is used to determine the overlap state of rigid blocks. This iterative scheme uses the concept of the Minkowski difference to efficiently detect the overlap status of convex shapes in both 2D and 3D. The full Minkowski difference of two closed, convex objects can be constructed by taking the entire set of vertices on the surface of object A and subtracting them from the entire set of vertices on the other object, object B. Thus, if object A has N vertices on its exterior surface and object B has M vertices on its exterior surface, the convex hull of the N * M vertices forms the full Minkowski difference. The Minkowski difference has the interesting property that, if the origin is inside the Minkowski difference, the objects overlap. In addition, the vector from the origin to the closest point on Minkowski difference is the minimum vector by which the objects must be displaced from overlapping so that they come into touching contact; this is termed the penetration vector. Note that, in general, the penetration vector may not be unique. If the convex objects are not overlapping, on the other hand, the vector from the origin to the closest point on the Minkowski difference denotes the separation vector, or the vector by which the objects must be displaced from separation to come into touching contact.

Instead of determining the full Minkowski difference, the GJK algorithm attempts to build a simplex (i.e., a {triangle in 2D; tetrahrdron in 3D}) composed of vertices from the Minkowski difference that encloses the origin iteratively. If such a simplex exists then the objects used to create the Minkowski difference overlap. The GJK algorithm relies on obtaining points on the surface of the Minkowski difference without computing the entire surface of the Minkowski difference. This can be accomplished via the use of convex support points. A convex support point is the point on the surface of a convex object that is farthest in a specified direction; this point may not be unique. By querying support points on each of the convex objects, one can find a point on the Minkowski surface that is farthest in a specified direction. Note that the farthest point in a specified direction on the Minkowski difference may also not be unique. For rigid blocks, the determination of a support point requires finding the dot product of each vertex of the rigid block and the specified direction. As a result, the computational time of this algorithm depends directly on the number of vertices of each block.

The GJK algorithm begins by choosing an initial direction (usually the difference of the object centroids) and querying the point on the Minkowski difference in that direction. An incomplete simplex is created composed of that point. During each subsequent step the direction from the closest point on the incomplete simplex to the origin is used to determine a new point on the Minkowski difference. If any incomplete simplex contains the origin the objects overlap. When a complete simplex is created, the algorithm checks whether or not the origin is inside the simplex. If it is not, the algorithm reduces to the largest incomplete simplex closest to the origin. If no other support points are closer to the origin than those already in the simplex the objects do not collide.

The GJK algorithm is sufficient to determine the separation vector. When objects overlap, though, an additional algorithm is required to compute the penetration vector. When the origin is inside the Minkowski difference, the Expanding Polytope Algorithm (EPA) ([Bergen1999]) allows one to determine the closest point on the surface to the origin. This algorithm is fundamentally similar to the Quickhull ([Barber1996]) algorithm to compute the convex hull of a set of points in space. An additional point on the surface of the Minkowski difference is computed based on the facet outward facing normal of the simplex closest to the origin. The simplex is expanded into a hull to include the new point by computing the faces for which the point is located on the opposite side of the origin. These faces are removed and new faces are created to close the hull. The hull is incrementally expanded in the fashion until the newly selected point to add is already in the hull. This incremental procedure allows one to effectively build a local representation of the Minkowski difference to determine the penetration vector.

In order to speed computation of the narrow-phase contact detection with rigid blocks, rigid block rounding has been introduced in an attempt to eliminate the EPA algorithm (see the rblock erode and rblock dilate commands). If rounding is used it is generally advisable to round the templates used to generate a packing of rigid blocks or to round before cycling has occurred. When rounding is active and the overlap is less than the sum of the rigid block rounding distances, the EPA algorithm is unnecessary as the separation vector is sufficient to determine the penetration vector with rounding.

Rigid Block Contacts with Wall Facets

Walls are composed of {linear facets in 2D; triangular facets in 3D}. The current implementation adds a point at the center of the facet to the set of points that are used in the GJK and EPA algorithms as those algorithms are generic and can be applied to planar geometries as well. Note that currently the wall resolution (see the wall resolution command) is not applied to contacts between wall facets and rigid blocks.

Contact Resolution

Once a contact exists and the narrow-phase contact detection determines that the pieces overlap (i.e., a penetration vector is computed and the gap is negative), the remaining information required to compute forces is the contact location. For contacts with rigid blocks, there is only one contact location where forces are computed. This is in contrast to the logic used in UDEC and 3DEC where multiple subcontacts exist between contacting convex blocks. The reason for this difference is to avoid the potentially significant memory and computational overhead associated with numerous contact locations between each pair of contacting objects. For contacts with rigid blocks there are two modes of contact resolution, namely reduced and full (see the rblock contact-resolution command). Both modes will be discussed in detail for each contact type.

In addition to these modes of contact resolution, one must also deal with updating the contact contact locations bonding is active. When bonding is active the contact location at the time of bonding is used as an initial location that is incrementally updated in response to the relative velocities at the contact location. This location may differ significantly from the contact location determined by the methods outlined below. If one conceives of bonds as providing a mechanism to glue pieces together, such incremental contact location makes sense as the material of the bond conceptually adheres to the surfaces of the pieces at the time of bonding. This conception requires that rigid blocks overlap for bonds to be installed. One might use a combination of the rblock erode and rblock scale commands to prepare a zero-porosity packing of rigid blocks for bonding.

Rigid Block Contacts with Balls/Pebbles

Narrow-phase contact detection produces the closest point on the rigid block surface to the centroid of the ball or pebble. The contact location can be readily computed as the average of the the point on the ball surface in the direction of the closest-point position and the closest-point position. This definition is consistent with the ball-ball, ball-pebble, and ball-facet contact locations.

Rigid Block Contacts with Rigid Blocks

Determination of a singular contact location for rigid block contacts is a difficult task. As noted in [Feng2012], for rigid blocks to be at rest with a single contact the contact location must be placed at the centroid of the overlap {area in 2D; volume in 3D}.

Computing the overlap for 2D convex polygons can be accomplished via the plane sweep as outlined in [Shamos1976]. This marching algorithm requires the points be ordered with increasing direction in one of the global axes directions. By partitioning the vertices into slabs based on the ordered vertices of both polygons, one can compute the overlap polygon.

Computing the overlap for 3D polyhedra, on the other hand, is a more involved computation. One method for doing this is outlined by [Preparata1979]. An affine transformation is applied to the facet normals to produce a new set of points in the transformed space. The transformation requires that a point entirely inside the overlap volume be specified. The previous contact location incrementally updated by the relative velocity is often located within the overlap volume. If this point fails to be entirely inside the overlap volume, a new point must be computed. In this situation the rigid blocks are sliced by a plane using the previous contact location; the [Shamos1976] method is used to find the center of the overlap area which is likely a point entirely inside the intersection volume. With a viable set of points in the transformation space defined, the convex hull is computed via the Quickhull algorithm (Barber et al. 1996). Note that it is often the case, for very small penetration, that the transformed points are nearly colinear, highlighting the need for a robust convex hull implementation. PFC uses an internally developed Quickhull implementation for this purpose. Once a convex hull in transformation space is defined, the reverse of the affine transformation is applied to the facet normals of the convex hull to produce the points on the exterior of the volume of intersection. The transformed hull facet connectivity can be used to deduce the connectivity of the convex hull of the intersection volume in real space. This computation depends directly on the number of block vertices. For systems with a fair amount of spatial coherence, one can speed this computation up by caching information about the affine transformation hull, checking to see whether or not the previous hull is still valid with the new block positions. This is quite an effective strategy to reduce the contact resolution time.

The previous discussion outlines the full method of contact resolution. When in reduced mode, the exact centroid of the overlap region is not computed. Instead the points of deepest penetration are determined from the penetration vector and the plane halfway between these points is deemed the contact plane. Both blocks are cut by this plane and the contact location is set as the centroid of the {line of intersection in 2D; area of intersection in 3D} of these cuts. This computation can also be sped up by caching the facet where the cut starts between successive timesteps. By traversing the facet connectivity information the cut can be determined quickly once an initial cut has been identified.

Either the full or reduced method can also be used when rounding is active. In this case the core shapes are displaced in equal amounts so that they overlap by the computed gap and the schemes outlined above are applied.

Rigid Block Contacts with Wall Facets

In order for the preceding logic to be applied to contacts with wall facets, one must create a convex shape out of the facet. This is accomplished by extruding the facet in the penetration direction away from the rigid block centroid. Such an extrusion allows for contact locations with facet edges to act as expected and is also consistent with using rigid blocks in the place of walls.

References

[Barber1996]Barber, C.B., Dobkin, D.P., and H. Huhdanpaa, “The quickhull algorithm for convex hulls,” ACM Transactions on Mathematical Software, 22(4), 469-483, 1996.
[Bergen1999]van den Bergen, G. “A fast and robust GJK implementation for collision detection of convex objects,” Journal of Graphics Tools, 4(2), 7-25, 1999.
[Feng2012]Feng, Y.T., Han, K., and D.R.J. Owen, “Energy-conserving contact interaction models for arbitrarily shaped discrete elements,” Computer Methods in Applied Mechanics and Engineering, 205–208, 169-177, 2012.
[Gilbert1988]Gilbert, E.G., D.W. Johnson and S.S. Keerthi, “A fast procedure for computing the distance between convex objects in three-dimensional space,” IEEE Journal of Robotics and Automation, 4(2), 193-203, 1988.
[Nezami2006]Nezami G., Hashash Y., Zhao D., and J. Ghaboussi, “Shortest link method for contact detection in discrete element method,” International Journal for Numerical and Analytical Methods in Geomechanics, 30(8), 783–801, 2006.
[Preparata1979]Preparata, F.P. and D.E. Muller, “Finding the intersection of n half-spaces in time O(nlogn),” Theoretical Computer Science, 8(1), 45-55, 1979.
[Shamos1976](1, 2) Shamos, M.I. and D. Hoey, “Geometric intersection problems,” SFCS ‘76 Proceedings of the 17th Annual Symposium on Foundations of Computer Science, October 25 - 27, 208-215, 1976.