The C++ API

Note

functions whose name is prefixed with an underscore should be considered as “private”: these functions are exposed for testing purposes. They should not be used, since they are susceptible of incompatible changes (or even removal) in future versions.

Representation of vectors and matrices

An ellipsoid is defined from its center c (a 3×1, column-vector) and quadratic form Q (a 3×3, symmetric, positive definite matrix) as the set of points m such that:

(m-c)ᵀ⋅Q⁻¹⋅(m-c) ≤ 1.

In this module, objects referred to as “vectors” are double[3] arrays of coordinates. In other words, the representation of the vector x is the double[3] array x such that:

    ⎡ x[0] ⎤
x = ⎢ x[1] ⎥.
    ⎣ x[2] ⎦

Objects referred to as “symmetric matrices” (or “quadratic forms”) are of type double[6]. Such arrays list in row-major order the coefficients of the triangular upper part. In other words, the representation of a the symmetric matrix A is the double[6] array a such that:

    ⎡ a[0] a[1] a[2] ⎤
A = ⎢      a[3] a[4] ⎥.
    ⎣ sym.      a[5] ⎦

API

namespace pw85

Functions

void _cholesky_decomp(double const *a, double *l)

Compute the Cholesky decomposition of a symmetric, positive matrix.

Let A be a symmetric, positive matrix, defined by the double[6] array a. This function computes the lower-triangular matrix L, defined by the double[6] array l, such that Lᵀ⋅L = A.

The array l must be pre-allocated; it is modified by this function. Note that storage of the coefficients of L is as follows

    ⎡ l[0]    0    0 ⎤
L = ⎢ l[1] l[3]    0 ⎥.
    ⎣ l[2] l[4] l[5] ⎦

void _cholesky_solve(const double *l, const double *b, double *x)

Compute the solution to a previously Cholesky decomposed linear system.

Let L be a lower-triangular matrix, defined by the double[6] array l (see pw85::_cholesky_decomp() for ordering of the coefficients). This function solves (by substitution) the linear system Lᵀ⋅L⋅x = b, where the vectors x and b are specified through their double[3] array of coordinates; x is modified by this function.

void spheroid(double a, double c, const double *n, double *q)

Compute the quadratic form associated to a spheroid.

The spheroid is defined by its equatorial radius a, its polar radius c and the direction of its axis of revolution, n (unit-vector, double[3] array).

q is the representation of a symmetric matrix as a double[6] array. It is modified in-place.

double f_neg(double lambda, const double *r12, const double *q1, const double *q2)

Return the value of the opposite of the function f defined as (see Theory).

\[f(\lambda)=\lambda\bigl(1-\lambda\bigr)r_{12}^{\mathsf{T}} \cdot Q^{-1}\cdot r_{12},\]

with

\[Q = \bigl(1-\lambda\bigr)Q_1 + \lambda Q_2,\]

where ellipsoids 1 and 2 are defined as the sets of points \(m\) (column-vector) such that

\[\bigl(m-c_i\bigr)\cdot Q_i^{-1}\cdot\bigl(m-c_i\bigr)\leq1.\]

In the above inequality, \(c_i\) is the center; \(r_{12}=c_2-c_1\) is the center-to-center radius-vector, represented by the double[3] array r12. The symmetric, positive-definite matrices \(Q_1\) and \(Q_2\) are specified through the double[6] arrays q1 and q2.

The value of \(\lambda\) is specified through the parameter lambda.

This function returns the value of \(−f(\lambda)\) (the “minus” sign comes from the fact that we seek the maximum of \(f\), or the minimum of \(−f\)).

This implementation uses Cholesky decompositions*.

void _residual(double lambda, const double *r12, const double *q1, const double *q2, double *out)

Compute the residual \(g(\lambda)=\mu_2^2-\mu_1^2\).

See Optimization of the function f for the definition of \(g\). The value of \(\lambda\) is specified through the parameter lambda. See contact_function() for the definition of the parameters r12, q1 and q2.

The preallocated double[3] array out is updated as follows: out[0]= \(f(\lambda)\), out[1]= \(g(\lambda)\) and out[2]= \(g'(\lambda)\).

This function is used in function pw85::contact_function() for the final Newton–Raphson refinement step.

int contact_function(const double *r12, const double *q1, const double *q2, double *out)

Compute the value of the contact function of two ellipsoids.

The center-to-center radius-vector \(r_{12}\) is specified by the double[3] array r12. The symmetric, positive-definite matrices \(Q_1\) and \(Q_2\) that define the two ellipsoids are specified through the double[6] arrays q1 and q2.

This function computes the value of \(\mu^2\), defined as

\[\mu^2=\max_{0\leq\lambda\leq 1}\bigl\{\lambda\bigl(1-\lambda\bigr) r_{12}^{\mathsf{T}}\cdot\bigl[\bigl(1-\lambda\bigr)Q_1+\lambda Q_2\bigr]^{-1} \cdot r_{12}\bigr\},\]

and the maximizer \(\lambda\), see Theory . Both values are stored in the preallocated double[2] array out: out[0] = \(\mu^2\) and out[1] = \(\lambda\).

\(\mu\) is the common factor by which the two ellipsoids must be scaled (their centers being fixed) in order to be tangentially in contact.

This function returns 0.

Todo

This function should return an error code.

Variables

constexpr size_t dim = 3

The dimension of the physical space (3).

constexpr size_t sym = 6

The dimension of the space of symmetric matrices (6).

constexpr double lambda_atol = 1e-6

The absolute tolerance for the stopping criterion of Brent’s method (in function contact_function()).

constexpr size_t max_iter = 25

The maximum number of iterations of Brent’s method (in function contact_function()).

constexpr size_t nr_iter = 3

The total number of iterations of the Newton–Raphson refinement phase (in function contact_function()).

namespace metadata

Variables

constexpr std::string_view author = {"S. Brisard"}
constexpr std::string_view description = {"Implementation of the \"contact function\" defined by Perram and Wertheim (J. Comp. Phys. 58(3), 409-416, DOI:10.1016/0021-9991(85)90171-8) for two ellipsoids."}
constexpr std::string_view author_email = {"sebastien.brisard@univ-eiffel.fr"}
constexpr std::string_view license = {"BSD 3-Clause License"}
constexpr std::string_view name = {"pw85"}
constexpr std::string_view url = {"https://github.com/sbrisard/pw85"}
constexpr std::string_view version = {"2.0"}
constexpr std::string_view year = {"2021"}