Fast contact detection for ellipsoids using optimization approaches

The paper is concerned with fast and robust contact detection methods for arbitrary ellipsoids. An iterative procedure, namely the Levenberg‐Marquardt method, based on the common normal concept for parametric ellipsoids, is employed together with an implementation of the widely used GJK algorithm for comparison. The performance and accuracy of both are analysed and compared to each other on the basis of two test sets, each containing a total of 70 000 pairs of prolates or oblates. Emphasis is placed on the specific error measure relating the iterative solution to the exact one, which was chosen to be the maximum angle between the normal vector and the distance vector between two ellipsoids. The results indicate increased performance when using the Levenberg‐Marquardt method over the GJK algorithm with no loss of accuracy.


INTRODUCTION
The transport of particles in fluid flow plays an extensive role in both industrial [1] and environmental [2] processes.Spatially resolved simulations of particle-laden flows, considering the interactions between the particles and the fluid as well as among the particles themselves, allow the physics underlying these processes to be studied.Recent simulations of this kind considering sediment transport with ellipsoidal particles indicate a significant influence of the particle shape on the characteristics of the flow [3] as well as on the motion of the particles themselves [4].These findings were obtained using direct numerical simulations of a channel flow with different ellipsoidal shapes, with the phase-resolved particles being coupled to the fluid using the immersed boundary method (IBM) [5].While the IBM was employed for fluid-particle coupling, the inter-particle forces were computed using an impulse-based collision model [6] that relies on the knowledge of the points of contact, which were determined using an iterative procedure.It accounts for both, forces created by direct contact between the particle surfaces, as well as lubrication forces.In fact, the latter are highly relevant in such a situation and require the distance between the particle surfaces to be known, even if these do not touch, as this is the input for the lubrication model.A similar approach shall be used in future simulations as well, with a substantially larger number of particles.Hence, an efficient and robust contact detection method for ellipsoids has to be available.Contact detection methods for ellipsoids can be classified by the definition of the contact points and by the their surface descriptions.Regarding the contact point, either the common normal concept or the geometric potential concept can be applied [7].The former aims to minimise the Euclidean distance between the two ellipsoids, while the latter minimises the geometric potential.These approaches can be combined with either an implicit or a parametric surface description to generate root-finding or optimisation problems of varying dimension and complexity [8].It is understood that such a contact detection is employed in cases where particle surfaces are in direct contact, as well as for the situation that the surfaces are still apart but close.In this case the method has to determine the shortest distance between the surfaces, as well as two points, one on each surface, constituting the anchor of the segment realising the closest distance.These points are called contact points here.Ellipsoids are strictly convex, so that a unique set of contact points exists.
It has to be recognised that the multiphase-flow community of researchers has long employed relatively simple algorithms for contact detection, while in the computer graphics community other methods are being used.One of the algorithms widespread in the latter is the Gilbert-Johnson-Keerthi (GJK) algorithm [9].While generally devised for piecewise linear surfaces it has recently been applied to analytically defined ellipsoids with promising results [10].For this reason it has been implemented during the present study as well in the same computational framework, so as to allow a comparison of execution times.
In the present contribution, the common normal concept together with a parametric surface description is discussed, with the GJK algorithm as a reference.The methods are tested and compared on large sets of prolate and oblate ellipsoids, providing essential insights into the performance and accuracy of both.

Statement of contact problem
To determine the contact points for a large number of ellipsoids, individual pairs are considered separately.Such a pair may lead to three possible scenarios: Either a real distance occurs between their surfaces (gap), or they touch at exactly one point (perfect contact), or they penetrate each other (intersection).For the purposes of this study, only the first case is investigated.This, however, entails a certain difficulty.If the contact pair has no actual contact point, a virtual contact point must be specified.The virtual contact points are, thus, defined as the points that would most likely touch each other, as described in the introduction.The problem of finding the points of contact is, therefore, embedded into the problem of finding the points  * 1 and  * 2 on the surfaces of two ellipsoids  1 and  2 the distance  * of which is minimal.This requirement can be formulated in terms of the constrained optimization problem as Here, ( 1 ,  2 ) =  2 −  1 is the vector connecting the two points and  = 1 2 |( 1 ,  2 )| 2 is the objective function to be minimized.The two constraints ensure that the points are located on the surfaces of the respective ellipsoids.To fulfill these conditions, the parametric surface description for ellipsoids   =   (  ,   ) is introduced, with  = 1, 2 referring to one of the two ellipsoids.This allows transforming the constrained optimization problem into an unconstrained one, that is where the four angles  1 ,  1 ,  2 ,  2 are summarized by the vector  .A necessary condition for an optimal solution to (2) is a vanishing gradient of the objective function.This yields where the gradient ∇ comprises the derivatives of the two points  1 and  2 , and the scalar product of each of them with the distance vector  must vanish.This states that the tangent vectors F I G U R E 1 Two neighboring ellipsoids  1 and  2 together with the exact contact points in the context of the common normal concept.Both contact points have to share the same normal vector, apart from orientation, resulting in the requirement that each normal vector has to be co-linear to the distance vector.
are perpendicular to the distance vector at the exact points of contact.This fact identifies the present approach as a common normal method, meaning that the global solution to ( 2) is chosen based on the direction of the two surface normal vectors This is illustrated in Figure 1 for the two ellipsoids  1 and  2 .
The requirement that the distance ||, indeed, assumes a minimum and not a maximum can be rephrased into restrictions on the direction of the vector  with respect to the normal vectors  1 and  2 at the two contact points that is

Solution procedure
Due to its non-linear nature, problem (2) has to be solved by an iterative approach.Starting at an initial guess ( S ), successive positions (  ) are determined, using the step Δ  to decrease .For the present study, this was achieved by employing the Levenberg-Marquardt algorithm [11], resulting in which requires to solve a system of equations.This can be interpreted as a combination of the steepest descent and the Gauss-Newton method based on the blending parameter   .The latter is updated in every step based on the model quality, which relates the actual change in the objective function  +1 −   to the predicted change.The initial condition of the iterations is obtained from connecting the centers of the two ellipsoids and taking the intersection points of the two ellipsoids and the vector connecting their centers.

Convergence criterion
An important issue to be addressed here is the choice of the convergence criterion.It is used to terminate the iterative solution procedure if in some step  the computed distance ||  between the two ellipsoids is sufficiently close to the exact distance || * , such that ||  − || * ≤   .Since  * is not known for an arbitrary arrangement of two ellipsoids, the difference must be estimated.This estimation is realized using quantities which are readily given in each iteration of the solution procedure.Considering condition (6), the maximum angle between one of the two normal vectors and the distance vector (  , ) was chosen as the convergence criterion, such that the solution procedure is terminated if max =1,2 (  , ) ≤   (8) with 0 <   ≪ 1, and the final iteration called .

F I G U R E 2
Test setup to evaluate the error in the calculation of the angle  between two unit vectors  and  with the length of the former multiplied by .

F I G U R E 3
Relative error for decreasing exact angle θ, where θsin and θcos refer to the vector product and the scalar product computation, respectively.
Here, two possibilities for computing the angle   are considered.One is based on the scalar product The other is based on the vector product The criterion is, thus, closely related to the gradient of the objective function (6), where both the normal vectors are required to be parallel to the distance vector at the exact contact points, meaning that (  , ) = 0.The two possible implementations were evaluated in terms of their accuracy.For this reason, a test setup was constructed, as shown in Figure 2, and the relative error for the two implementations analysed for decreasing angles .To make the test more difficult, the parameter  was introduced making one of the vectors very short in all tests reported, with  = 10 −10 .It turned out, however, that the results did not change when  = 1 was employed.
The diagram in Figure 3 shows a linear increase in the relative error for decreasing angle when using (9).With (10), quadratic behaviour is observed.For very small angles, the relative error might be considered too large for being used in a complex multiphase computation.It is therefore necessary to assess whether this effect has an actual impact on the overall performance and accuracy of the contact detection method.

SETUP
To evaluate the accuracy of the method in a large-scale simulation, a framework inspired by [12] was employed that allows the construction of pairs of ellipsoids with exactly known contact points.This was used to create seven test groups, referring to seven different dimensionless distances  * ∕ eq,m , each containing a total of 10 4 pairs of ellipsoids.The sphericity Ψ = 3 √ ∕ 2 according to Krumbein [13] of the ellipsoids was fixed at Ψ = 2∕3 to agree with the simulations of Jain et al. [3].This leads to two cases, a case of oblate ellipsoids and a case of prolate ellipsoids.For each individual pair, the mean equivalent diameter  eq,m and the ratio of equivalent diameters  eq,a were drawn from a logarithmic uniform distribution.Then the three semi-axes , ,  of the two ellipsoids were constructed based on Ψ,  eq,m ,  eq,a , resulting in ellipsoids of different size and shape.In addition to the computed distance   and the contact points   1 and   2 after the final iteration, different error measures were determined to evaluate the accuracy of the procedure.Since the contact points are exactly known by construction, it is straightforward to compare these exact quantities with the solutions of the iterative procedure.The following measures were determined: which is the relative error in distance, the relative error in the position of the contact points on the two ellipsoids in their local coordinate systems, and the ratio of the two median execution times, respectively.For the latter, the own implementation of the GJK algorithm adapted to ellipsoids by Girault et al. [10] was used as a reference.Implementing both algorithms in the very same framework provides optimal conditions for assessment.

Prolate-prolate contact
The results for the medians of the measures (11) for the pairs of prolates are shown in Figure 4, 5 and 6 together with an evaluation of the robustness in Figure 7.In some of the figures only a single curve can be discerned, since continuous and broken lines strictly overlap.The findings justify the utilization of the angle max =1,2 (  , ) as the error measure for the iterative procedure, which can be observed in Figure 4 and 5.In the calculations performed, the medians of the maximum angles correspond to both the error in the distance    and the error in the contact points    .Moreover, med   decreases monotonically as the distance  *  eq,m increases, while med   remains approximately constant.This is due to the definition of   , given in (11).In addition, Figure 6 reveals that the Levenberg-Marquardt algorithm converges significantly faster than the GJK algorithm for small distances  * ∕ eq,m ≤ 10 −1 , while the opposite behavior is observed for distances  * ∕ eq,m = 10 0 .Such cases, however, are not of interest, since the distances are too large to be needed in a true simulation with ellipsoids of this sphericity.The robustness of the Levenberg-Marquardt algorithm is illustrated in Figure 7.It is in accordance to the evaluation of the computation time since it reveals a guaranteed convergence for distances  * ∕ eq,m ≤ 10 −1 .Finally, the algorithm failed to converge in 9 out of 10000 cases for which  * ∕ eq,m = 10 0 , corresponding to a percentage of 0.09%.If this should occur in an application appropriate measures can be taken, such as an under-relaxation of the increment, etc.

Oblate-oblate contact
The results for the medians of the measures (11) for the oblates are presented in the same way as the results of the prolates.This is shown in Figure 8, 9 and 10 together with the assessment of the robustness in Figure 11.
A very similar pattern to the above case can be observed as well, with the Levenberg-Marquardt algorithm failing to converge to the contact points in only 4 cases for  * ∕ eq,m ≤ 10 −1 .These cases can be addressed as indicated above.Hence, it can be concluded that the medians of the relative performance and accuracy of the algorithms are independent of the specific configuration of the ellipsoids with the given shapes, for the given sphericity.