.. _theory-geom-ellipsoid: Ellipsoid ----------- :obj:Elllipsoids  are defined by their local coordinates and three :obj:semi-axes . Ellipsoid-Ellipsoid contact ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Thsi algorithm is implemented in by :obj:woo.dem.Cg2_Ellipsoid_Ellipsoid_L6Geom. Following :cite:Perram1996, we note :obj:semi-axes  :math:a_i and :math:b_i as :math:i\in\{1,2,3\}; local axes are given as sets of orthonormal vectors :math:\vec{u}_i, :math:\vec{v}_i. The intercenter vector reads :math:\vec{R}=\vec{r}_b-\vec{r}_a, with :math:\vec{r}_a and :math:\vec{r}_b being positions of ellipsoid centroids. We define matrices .. math:: :nowrap: \begin{align*} \mat{A}&=\sum_{i=1}^{3} a_i^{-1}\vec{u}_i\vec{u}_i^T \\ \mat{B}&=\sum_{i=1}^{3} b_i^{-1}\vec{v}_i\vec{v}_i^T \end{align*} which are both invertible for non-vanishing semi-axes; then we define the function .. math:: :label: perram-wertheim-S-lambda S(\lambda)=\lambda(1-\lambda)\vec{R}^T\big[\underbrace{(1-\lambda)\mat{A}^{-1}+\lambda\mat{B}^{-1}}_{\mat{G}}\big]^{-1}\vec{R} where :math:\lambda\in[0,1] is an unknown parameter. The function :math:S(\lambda) is concave and non-zero for :math:\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 :math:\Lambda for which :math: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 :math:S'(\lambda), and re-using the value from the previous step as the initial guess. There are also papers suggesting better algorithms such as :cite:Zheng2009. We define the Perram-Wertham potential (first introduced in :cite:Perram1985) as .. math:: F_{AB}=\{\max S(\lambda)|\lambda\in[0,1]\}=S(\Lambda) for which (:cite:Perram1985) * :math:F_{AB}<1 if both ellipsoids overlap, * :math:F_{AB}=1 if they are externaly tangent, * :math:F_{AB}>1 if the ellipsoids do not overlap. :cite:Perram1985 gives a geometrical interpretation of :math:F_{AB}, where .. math:: F_{AB}=\mu^2 where :math:\mu is scaling factor which must be applied to both ellipses to be tangent. Following :cite:Donev2005 eq (18), we can compute the contact normal and the contact point of two ellipsoids as .. math:: :nowrap: \begin{align*} \vec{n}_c&=\mat{G}^{-1}\mat{R} \\ \vec{r}_c&=\vec{r}_a+(1-\Lambda)\mat{A}^{-1}\cdot \end{align*} with :math:\mat{G} defined in :eq:perram-wertheim-S-lambda; note that :math:\vec{n}_c is not normalized. The penetration depth (overlap distance) can be reasoned out as follows. :math:\mu scaled ellipsoid sizes while keeping their distance, so that they become externally tangent. Therefore :math:1/\mu scales ellipsoid distance while keeping their sizes. With :math:d=|\mat{R}| being the current inter-center distance, we obtain .. math:: u_n'=d-d_0=d-\frac{1}{\mu}d=d\left(1-\frac{1}{\mu}\right). This is the displacement that msut be performed along :math:\vec{R} while the contact normal may be oriented differently; we therefore project :math:u_n' along :math:\vec{R} onto (normalized) :math:\vec{n}_c obtaining .. math:: u_n=d\left(1-\frac{1}{\mu}\right)\normalized{\vec{R}}\cdot\normalized{\vec{n}_c}. The :math:u_n, :math:\normalized{\vec{n}}, :math:\vec{r}_c can be fed to :ref:contact-geometry-l6gom-generic for further computation. Ellipsoid-Wall intersection ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. todo:: Write.