Ellipsoid¶
Elllipsoids
are defined by their local coordinates and three semi-axes
.
Ellipsoid-Ellipsoid contact¶
Thsi algorithm is implemented in by woo.dem.Cg2_Ellipsoid_Ellipsoid_L6Geom
.
Following [PRPraestgaardL96], we note semi-axes
\(a_i\) and \(b_i\) as \(i\in\{1,2,3\}\); local axes are given as sets of orthonormal vectors \(\vec{u}_i\), \(\vec{v}_i\). The intercenter vector reads \(\vec{R}=\vec{r}_b-\vec{r}_a\), with \(\vec{r}_a\) and \(\vec{r}_b\) being positions of ellipsoid centroids. We define matrices
which are both invertible for non-vanishing semi-axes; then we define the function
where \(\lambda\in[0,1]\) is an unknown parameter. The function \(S(\lambda)\) is concave and non-zero for \(\lambda\in[0,1]\) and zero at both end-points. It therefore has a sigle maximum. Using the Brent’s method, we obtain the maximum value \(\Lambda\) for which \(S(\Lambda)\) is maximal.
Todo
This computation can be (perhaps substantially) sped up by using other iterative methods, such as Newton-Raphson, finding root of \(S'(\lambda)\), and re-using the value from the previous step as the initial guess. There are also papers suggesting better algorithms such as [ZIPM09].
We define the Perram-Wertham potential (first introduced in [PW85]) as
for which ([PW85])
\(F_{AB}<1\) if both ellipsoids overlap,
\(F_{AB}=1\) if they are externaly tangent,
\(F_{AB}>1\) if the ellipsoids do not overlap.
[PW85] gives a geometrical interpretation of \(F_{AB}\), where
where \(\mu\) is scaling factor which must be applied to both ellipses to be tangent.
Following [DTS05] eq (18), we can compute the contact normal and the contact point of two ellipsoids as
with \(\mat{G}\) defined in (1); note that \(\vec{n}_c\) is not normalized.
The penetration depth (overlap distance) can be reasoned out as follows. \(\mu\) scaled ellipsoid sizes while keeping their distance, so that they become externally tangent. Therefore \(1/\mu\) scales ellipsoid distance while keeping their sizes. With \(d=|\mat{R}|\) being the current inter-center distance, we obtain
This is the displacement that msut be performed along \(\vec{R}\) while the contact normal may be oriented differently; we therefore project \(u_n'\) along \(\vec{R}\) onto (normalized) \(\vec{n}_c\) obtaining
The \(u_n\), \(\normalized{\vec{n}}\), \(\vec{r}_c\) can be fed to Generic contact routine for further computation.