Biswajit Banerjee

Generating periodic RVEs with polydisperse ellipsoids

Part 1: Finding intersections


Many micromechanics studies use periodic representative volume elements to simulate the effective response of a material. The periodicity requirement comes from requirements imposed by the asymptotic theory of homogenization. In that theory, we assume that field quantities in a heterogeneous material can be expressed as functions of the form where is periodic in . Such functions are identified with microscopic values and macroscopic quantities are found by averaging these values.

The above approach is not the only possible way forward. The H-measure theory developed by Luc Tartar and Francois Murat does not have a periodicity requirement but is more difficult to simulate and has, therefore, not attracted the attention it deserves. We will also avoid the H-convergence route in this article. For more information on possible approaches, see Approximation of H-measures.

Periodic microstructures

polydisperse spheres
Figure 1. Polydisperse sphere microstructure generated using a Poisson process (c. 2002).

A common technique used to generate non-periodic microstructures uses the following steps:

  • Create a regular grid in three dimensions and bound the grid with an axis-aligned box
  • Place particles at each grid point based on an input particle size distribution
  • Initialize each particle with a random velocity
  • Run a discrete element or molecular dynamics simulation with rigid particles while moving the boundaries inward until a given packing fraction is achieved.

Once a non-periodic distribution is generated, a periodic one can be created from it by copying particles from one face of the RVE to the opposite face. This process is straightforward if particles are spherical, but more work is needed for ellipsoidal particles.

Ellipsoid-plane intersection

original ellipsoid 2d original ellipsoid 3d
Figure 2. An ellipsoid intersecting an oriented box.

The first step in the periodic particle generation is to identify which particles are boundary particles and need to be copied. An easy approach is to choose a bounding box around each face and select particles whose centers fall inside those boxes. However, there is the possibility of missing some particles if we use that approach.

If we start with an ellipsoidal particle distribution, another way of identifying which particles are at or near a boundary face is to determine whether the ellipsoid intersects the plane face. The advantage of this approach is that we can use it even if the bounding box is not axis-aligned.

We can determine whether an ellipsoid intersects a bounded plane, by using a transformation that takes the ellipsoid to a sphere. This transformation has the advantage that now we can can use sphere-edge intersection algorithms to determine whether the ellipsoid intersects a plane.

The transformation

The transformation is straightforward. We just compute the matrix , where

$$ \mathbf{T} = \begin{bmatrix} e^a_x/r_a & e^a_y/r_a & e^a_z/r_a \\ e^b_x/r_b & e^b_y/r_b & e^b_z/r_b \\ e^c_x/r_b & e^c_y/r_b & e^c_z/r_b \end{bmatrix} $$
original ellipsoid 2d
Figure 3. The transformed ellipsoid.

If the center of the ellipsoid is at , and the bounded plane has vertices , we get the transformed plane by computing the transformed vertices

$$ \mathbf{v}^{(t)}_i = \mathbf{T} (\mathbf{v}_i - \mathbf{c}) $$

where the radii of the three axes of the ellipsoid are and the unit vectors along the axes are . Each unit vector has components .

Notice that the rectangular box has transformed into a parallelepiped.

Finding face-sphere intersections

To find whether the transformed face intersects the sphere, I decided to compute not only the distance from the sphere to the face but also the location of the point of intersection of a line from the center of the sphere to each edge. Of course, such a line would have to be normal to the edge for that information to be useful.

First I computed a normal to the transformed face, assuming a counter-clockwise layout of the vertices around the outward normal to the face. The normal is given by

$$ \mathbf{n} = (\mathbf{v}_2 - \mathbf{v}_1) \times (\mathbf{v}_3 - \mathbf{v}_1) ~,~~ \hat{\mathbf{n}} = \mathbf{n}/\left\lVert\mathbf{n}\right\rVert $$

We, of course, use the unit normal, , in our calculations. The distance from the center of the sphere (which is at the origin after the transformation) to the plane is given by

$$ d = \hat{\mathbf{n}} \cdot (\mathbf{c} - \mathbf{v}_1) $$

and therefore, the projection of the sphere center on to the face is

$$ \mathbf{c}_p = \mathbf{c} - d \hat{\mathbf{n}} $$

Now that we have the closest point on the face (from the sphere), we can find the distances to the four edges of the bounded plane. If we parameterize an edge with the parameter such that any point on the edge is given by , the value of for the closest point on the edge from is

$$ t_c = -\frac{(\mathbf{v}_1 - \mathbf{c}_p)\cdot(\mathbf{v}_2 - \mathbf{v}_1)}{\left\lVert\mathbf{v}_2 - \mathbf{v}_1\right\rVert^2} $$

and the distance to that point is

$$ d_c = \left\lVert \mathbf{c}_p - (1 - t_c) \mathbf{v}_1 - t_c \mathbf{v}_2\right\rVert $$

We could alternatively compute using

$$ d_c = \frac{\left\lVert (\mathbf{c}_p - \mathbf{v}_2) \times (\mathbf{c}_p - \mathbf{v}_1) \right\rVert}{\left\lVert\mathbf{v}_2 - \mathbf{v}_1\right\rVert} $$

With this information, we can find whether the sphere intersects the transformed face if the distances are less than 1 (note that the transformation matrix takes the ellipsoid into a sphere of radius 1).

Let us now look at my C++ implementation of these ideas.


I did not want to lose any of the information computed during the distance computations. So I decided to return a pair containing the values of and for each edge. I created an Edge struct as shown below.

struct Edge
  std::array<Vertex, 2> vertex;
  // Distance computation
  std::pair<REAL, REAL> distance(const Point& point) const
    Vertex v1v0 = vertex[0] - point;
    Vertex v2v0 = vertex[1] - point;
    Vertex v2v1 = vertex[1] - vertex[0];
    REAL lengthSq = v2v1.lengthSq();
    REAL t = -dot(v1v0, v2v1)/lengthSq;
    REAL dist = cross(-v1v0, -v2v0).length()/std::sqrt(lengthSq);
    return std::make_pair(t, dist);

The distance function returns a pair such that first contains while second contains .

I then created a Face struct that had the following form. Once again, instead of returning only the distance to the face, I returned a pair containing the and values for each edge in addition to the distance to the face.

struct Face
  enum class Location
    NONE = 0, VERTEX = 1, EDGE = 2, FACE = 3
  std::array<Vertex, 4> vertex;
  // Normal computation
  Vec normal() const
    Vec normal = cross(vertex[1] - vertex[0], vertex[2] - vertex[0]);
    return normal;
  // Distance computation
  std::pair<std::vector<std::pair<REAL, REAL>>, REAL>
    distance(const Point& point) const
    Vec normal = this->normal();
    REAL dist = dot(normal, (point - vertex[0]));
    Point pointOnFace = point - normal*dist;
    Edge edge1(vertex[0], vertex[1]);
    std::pair<REAL, REAL> t1d1  = edge1.distance(pointOnFace);
    Edge edge2(vertex[1], vertex[2]);
    std::pair<REAL, REAL> t2d2  = edge2.distance(pointOnFace);
    Edge edge3(vertex[2], vertex[3]);
    std::pair<REAL, REAL> t3d3  = edge3.distance(pointOnFace);
    Edge edge4(vertex[3], vertex[0]);
    std::pair<REAL, REAL> t4d4  = edge4.distance(pointOnFace);
    return std::make_pair(
      std::vector<std::pair<REAL, REAL>>({t1d1, t2d2, t3d3, t4d4}), dist);

Now we are ready to find the intersection of the face with the ellipsoid. For the purposes of periodic particle generation, it is important that we know whether the ellipsoid intersects a vertex or edge of the face, or whether the region of intersection is fully contained in the face. We return that information from the intersection function in addition to a boolean that indicates whether there is any intersection at all. The code is listed below.

std::pair<bool, std::pair<Face::Location, int> >
Ellipsoid::intersects(const Face& face) const
  // Compute transformation matrix (ellipsoid to sphere)
  Matrix3 N = toUnitSphereTransformationMatrix();
  // Create transformed face
  Face transformedFace(N * (face.v1() - d_center), N * (face.v2() - d_center), N * (face.v3() - d_center), N * (face.v4() - d_center));
  // Compute distance to face from (0, 0, 0)
  auto face_td = transformedFace.distance(Point(0, 0, 0));
  // If the distance is greater than 1 there is no intersection
  if (std::abs(face_td.second) > 1.0) {
    return std::make_pair(false, std::make_pair(Face::Location::NONE, 0));
  // Check sphere vertex intersections
  std::array<Vec, 4> vertices = { transformedFace.v1(), transformedFace.v2(), transformedFace.v3(), transformedFace.v4() };
  int index = 0;
  for (const auto& vertex : vertices) {
    if (vertex.lengthSq() < 1.0) {
      return std::make_pair(true, std::make_pair(Face::Location::VERTEX, index));
  // Check sphere edge intersections
  std::array<Edge, 4> edges = { Edge(transformedFace.v1(), transformedFace.v2()), Edge(transformedFace.v2(), transformedFace.v3()), Edge(transformedFace.v3(), transformedFace.v4()), Edge(transformedFace.v4(), transformedFace.v1()) };
  index = 0;
  for (const auto& edge : edges) {
    auto edge_td = edge.distance(Point(0, 0, 0));
    if (edge_td.second < 1.0) {
      if (!(edge_td.first < 0 || edge_td.first > 1)) {
        return std::make_pair(true, std::make_pair(Face::Location::EDGE, index));
  // If none of the edges are intersected
  if (!(face_td.first[0].first < 0 || face_td.first[0].first > 1 ||
        face_td.first[1].first < 0 || face_td.first[1].first > 1 ||
        face_td.first[2].first < 0 || face_td.first[2].first > 1 ||
        face_td.first[3].first < 0 || face_td.first[3].first > 1)) {
    return std::make_pair(true, std::make_pair(Face::Location::FACE, 0));
  return std::make_pair(false, std::make_pair(Face::Location::NONE, 0));


With these basic infrastructure objects in place, we are ready to proceed to the actual generation of periodic particles. The details will be provided in part 2 of this article.

If you have questions/comments/corrections, please contact banerjee at parresianz dot com dot zen (without the dot zen).