Tutorial

In this tutorial, we consider two ellipsoids, and check wether or not they overlap.

Ellipsoid E₁ is an oblate spheroid centered at point x₁ = (-0.5, 0.4, -0.7), with equatorial radius a₁ = 10, polar radius c₁ = 0.1 and polar axis (0, 0, 1).

Ellipsoid E₂ is a prolate spheroid centered at point (0.2, -0.3, 0.4), with equatorial radius a₁ = 0.5, polar radius c₁ = 5 and polar axis (1, 0, 0).

To carry out the overlap check, we must first create the representation of ellipsoids Eᵢ as quadratic forms Qᵢ (see Mathematical representation of ellipsoids). Convenience functions are provided to compute the matrix representation of a spheroid.

Note

In principle, the contact function implemented in PW85 applies to any ellipsoids (with unequal axes). However, at the time of writing this tutorial (2019-01-01), convenience functions to compute the matrix representation of a general ellipsoid is not yet implemented. Users must compute the matrices themselves.

We first check for the overlap of E₁ and E₂ using the Python wrapper of pw85. We will then illustrate the C API.

Python tutorial

The Python module relies on NumPy for passing arrays to the underlying C library. We therefore import both modules:

>>> import numpy as np
>>> import pw85

and define the parameters of the simulation:

>>> x1 = np.array([-0.5, 0.4, -0.7])
>>> n1 = np.array([0., 0., 1.])
>>> a1, c1 = 10, 0.1
>>> x2 = np.array([0.2, -0.3, 0.4])
>>> n2 = np.array([1., 0., 0.])
>>> a2, c2 = 0.5, 5.
>>> r12 = x2-x1

where r₁₂ is the vector that joins the center of the first ellipsoid, x₁, to the center of the second ellipsoid, x₂.

We use the function pw85.spheroid() to create the matrix representations q₁ and q₂ of the two ellipsoids. Note that these arrays must be preallocated:

>>> q1 = np.empty((6,), dtype=np.float64)
>>> pw85.spheroid(a1, c1, n1, q1)
>>> q1
array([ 1.e+02, -0.e+00, -0.e+00,  1.e+02, -0.e+00,  1.e-02])
>>> q2 = np.empty_like(q1)
>>> pw85.spheroid(a2, c2, n2, q2)
>>> q2
array([25.  ,  0.  ,  0.  ,  0.25,  0.  ,  0.25])

We can now compute the value of the contact function — see the documentation of pw85.contact_function():

>>> out = np.empty((2,), dtype=np.float64)
>>> pw85.contact_function(r12, q1, q2, out)
>>> mu2, lambda_ = out
>>> print('μ² = {}'.format(mu2))
>>> print('λ = {}'.format(lambda_))
μ² = 3.362706040638343
λ = 0.1668589553405904

We find that μ² > 1, hence μ > 1. In other words, both ellipsoids must be swollen in order to bring them in contact: the ellipsoids do not overlap!

Checking the output

The output of this simulation can readily be checked. First, we can check that q₁ and q₂ indeed represent the ellipsoids E₁ and E₂. To do so, we first construct the symmetric matrices Q₁ and Q₂ from their upper triangular part

>>> Q1 = np.zeros((3, 3), dtype=np.float64)
>>> i, j = np.triu_indices_from(Q1)
>>> Q1[i, j] = q1
>>> Q1[j, i] = q1
>>> Q1
array([[ 1.e+02, -0.e+00, -0.e+00],
       [-0.e+00,  1.e+02, -0.e+00],
       [-0.e+00, -0.e+00,  1.e-02]])

>>> Q2 = np.zeros_like(Q1)
>>> Q2[i, j] = q2
>>> Q2[j, i] = q2
>>> Q2
array([[25.  ,  0.  ,  0.  ],
       [ 0.  ,  0.25,  0.  ],
       [ 0.  ,  0.  ,  0.25]])

We can now check these matrices for some remarkable points, first for ellipsoid E₁

>>> Q1_inv = np.linalg.inv(Q1)
>>> f1 = lambda x: Q1_inv.dot(x).dot(x)
>>> f1((a1, 0., 0.))
1.0
>>> f1((-a1, 0., 0.))
1.0
>>> f1((0., a1, 0.))
1.0
>>> f1((0., -a1, 0.))
1.0
>>> f1((0., 0., c1))
0.9999999999994884
>>> f1((0., 0., -c1))
0.9999999999994884

then for ellipsoid E₂

>>> Q2_inv = np.linalg.inv(Q2)
>>> f2 = lambda x: Q2_inv.dot(x).dot(x)
>>> f2((c2, 0., 0.))
1.0
>>> f2((-c2, 0., 0.))
1.0
>>> f2((0., a2, 0.))
1.0
>>> f2((0., -a2, 0.))
1.0
>>> f2((0., 0., a2))
1.0
>>> f2((0., 0., -a2))
1.0

Note that in the above tests, we have omitted the centers of ellipsoids E₁ et E₂ (both ellipsoids were translated to the origin).

We will now verify the corectness of the value found for the scaling factor μ. To do so, we will find the coordinates of the contact point of the two scaled ellipsoids, and check that the normals to the two ellipsoids at that point coincide.

Although we use formulæ from the Theory section to find the coordinates of the contact point, x₀, it is not essential for the time being to fully understand this derivation. What really matters is to check that the resulting point x₀ is indeed the contact point of the two scaled ellipsoids; how the coordinates of this point were found is irrelevant.

From the value of λ returned by the function pw85.contact_function(), we compute Q defined by Eq. (10) in section Theory, as well as x = Q⁻¹⋅r₁₂

>>> Q = (1-lambda_)*Q1+lambda_*Q2
>>> x = np.linalg.solve(Q, r12)

From which we find x₀, using either Eq. (9a) or Eq. (9b) (and we can check that both give the same result)

>>> x0a = x1+(1-lambda_)*np.dot(Q1, x)
>>> x0a
array([ 0.16662271, -0.29964969, -0.51687799])
>>> x0b = x2-lambda_*np.dot(Q2, x)
>>> x0b
array([ 0.16662271, -0.29964969, -0.51687799])

We can now check that x₀ belongs to the two scaled ellipsoids, that we first define, overriding the matrices of the unscaled ellipsoids, that are no longer needed. We observe that if ellispoid Eᵢ is scaled by μ, then its matrix representation Qᵢ is scaled by μ², and its inverse Qᵢ⁻¹ is scaled by μ⁻².

>>> x0 = x0a
>>> Q1 *= mu2
>>> Q2 *= mu2
>>> Q1_inv /= mu2
>>> Q2_inv /= mu2

>>> x = x0-x1
>>> Q1_inv.dot(x).dot(x)
1.0000000000058238

>>> x = x0-x2
>>> Q2_inv.dot(x).dot(x)
0.9999999999988334

Therefore x₀ indeed belongs to both ellipsoids. We now compute the normal nᵢ to ellipsoid Eᵢ at point x₀. Since ellipsoid Eᵢ is defined by the level-set: (x-xᵢ)ᵀ⋅Qᵢ⁻¹⋅(x-xᵢ) = 1, the normal to Eᵢ is given by Qᵢ⁻¹⋅(x-xᵢ) (which is then suitably normalized)

>>> n1 = Q1_inv.dot(x0-x1)
>>> n1 /= np.linalg.norm(n1)
>>> n1
array([ 3.64031943e-04, -3.82067448e-04,  9.99999861e-01])

>>> n2 = Q2_inv.dot(x0-x2)
>>> n2 /= np.linalg.norm(n2)
>>> n2
array([-3.64031943e-04,  3.82067448e-04, -9.99999861e-01])

We find that n₁ = -n₂. Therefore, E₁ and E₂ are in external contact. QED

Follow this link to download the above Python script.

C++ tutorial

The Python interface to PW85 has been kept close to the undelying C++ API. The following C++ program (download source file) defines the two ellipsoids, then computes μ² and λ:

#include <iostream>
#include <array>

#include "pw85/pw85.hpp"

using Vec = std::array<double, pw85::dim>;
using SymMat = std::array<double, pw85::sym>;

int main() {
  Vec x1{-0.5, 0.4, -0.7};
  Vec n1{0., 0., 1.};
  double a1 = 10.;
  double c1 = 0.1;

  Vec x2{0.2, -0.3, 0.4};
  Vec n2{1., 0., 0.};
  double a2 = 0.5;
  double c2 = 5.;

  Vec r12;
  for (int i = 0; i < pw85::dim; i++) r12[i] = x2[i] - x1[i];

  SymMat q1, q2;
  pw85::spheroid(a1, c1, n1.data(), q1.data());
  pw85::spheroid(a2, c2, n2.data(), q2.data());

  std::array<double, 2> out;
  pw85::contact_function(r12.data(), q1.data(), q2.data(), out.data());
  std::cout << "mu^2 = " <<  out[0] << std::endl;
  std::cout << "lambda = " << out[1] << std::endl;
}

A CMakeLists.txt file is provided for the compilation of the tutorial using CMake. You can reuse it in one of your own projects (download):

cmake_minimum_required(VERSION 3.13)

project("tutorial" LANGUAGES CXX)
set(CMAKE_CXX_STANDARD 17)

add_executable(${PROJECT_NAME} ${PROJECT_NAME}.cpp)

find_library(MATH_LIBRARY m)
if (MATH_LIBRARY)
    target_link_libraries(${PROJECT_NAME} INTERFACE ${MATH_LIBRARY})
endif()

find_package(Boost REQUIRED)
target_link_libraries(${PROJECT_NAME} PRIVATE Boost::headers)

target_include_directories(${PROJECT_NAME} PRIVATE "../../include")

cd into the cpp_tutorial subdirectory. The provided example program should be compiled and linked against pw85:

$ mkdir build
$ cd build
$ cmake ..
$ cmake --build . --config Release

An executable called tutorial should be present in the build/Release subdirectory. On execution, it prints the following lines to stdout:

mu^2 = 3.36271
lambda = 0.166859