Biswajit Banerjee

Generating periodic RVEs with polydisperse ellipsoids

Part 3: An implementation


vertex numbering
Figure 1. Initial ellipsoid distribution.

In the Part 2, we discussed the various cases we have to consider to generate periodically distributed ellipsoids inside an axis-aligned bounding box.

A typical starting point for our algorithm can be seen in Figure 1. Our task is to generate periodic particles as discussed in Part 2.

In this article, we will discuss a few more details and a possible way of implementing the algorithm in C++.

Shrinking the bounding box

vertex numbering
Figure 2. A periodic ellipsoid distribution.

If we start with the initial bounding box we may run into situations where the faces of the bounding box touch or intersect only a few ellipsoids. That can lead to a sparse distribution of ellipsoids in the extended bounding box. We can see that in Figure 2 where the expanded bounding box in blue appears to have a lot of empty space. Excessive empty space implies that periodic boundary conditions on the boundary particles will not lead to the desired homogeneous deformation of the exterior of the RVE without a lot of extra checks on the locations of particles.

vertex numbering
Figure 3. A periodic ellipsoid distribution with a shrunk initial bounding box.

One solution to this problem is to move shrink the initial bounding box by a small amount so that each of the three faces () intersects a larger number of particles. There is a possibility that some of the original boundary particles will be missed in this procedure. But keep in mind that all we are interested in is periodicity, and a few missed small particles should not affect the particle distribution in the RVE significantly given that we are already changing the distribution by adding periodic particles. Of course, we do have to check that the particle size and mass distribution remains within reasonable bounds after the particle generation procedure.

Figure 3 shows the distribution of the periodic particles after the initial bounding box has been shrunk by a small factor. There are still gaps, but fewer particles will need special treatment at the boundaries with the new distribution.

vertex numbering
Figure 4. Orthographic projection of periodic ellipsoid distribution.

You can try to convince yourself that the new particle distribution is actually periodic by looking at the image in Figure 4 (orthographic projection). However, a computational check of overlaps is essential if the margin parameter (discussed in Part 2) is smaller than 2.


We will assume that we have arrays of particle positions, radii, and axis orientations and also a bounding box. We will also assume that we have a marginFactor to determine the locations where the periodic particles will be placed, and a faceShiftFactor that will determine how much the initial bounding box will have to be shrunk.

Max-min particle radii

The first step in the particle generation process is the determination of the minimum and maximum particle radii. The maximum radius will determine the minimum size of the margin while the minimum radius will determine how much bounding box shrinkage is acceptable. We will use the binary search based std::minmax_element for this purpose because of its complexity.

  auto minmaxIter = std::minmax_element(particles.begin(), particles.end(),
   [](const DEMParticleP& p1, const DEMParticleP& p2){
     auto max_p1 = std::max({p1->radiusA(), p1->radiusB(), p1->radiusC()});
     auto max_p2 = std::max({p2->radiusA(), p2->radiusB(), p2->radiusC()});
     return max_p1 < max_p2;
  auto minRadius = std::min({(*minmaxIter.first)->radiusA(),
  auto maxRadius = std::max({(*minmaxIter.second)->radiusA(),
  auto boundaryMargin = marginFactor*maxRadius;
Setting up the shrunk bounding box

We use the face shift factor and the minimum particle radius to compute the shrunk domain and get a vector of vertices of the domain.

  auto faceShift = faceShiftFactor*minRadius;
  Box shrunkDomain(
    spatialDomain.minCorner() + Vec(faceShift, faceShift, faceShift),
    spatialDomain.maxCorner() - Vec(faceShift, faceShift, faceShift));
  OrientedBox box(shrunkDomain);
  std::vector<Vec> vertices = box.vertices();
  auto widthX = shrunkDomain.dimX();
  auto widthY = shrunkDomain.dimY();
  auto widthZ = shrunkDomain.dimZ();
Setting up the three faces for intersections

Next we set up a set of face indices (see Part 2 for a discussion of the orientations of faces) and create three Face objects for convenience.

constexpr std::array<std::array<int, 4>, 3> faceIndices = {
      {0, 4, 7, 3} // x-
    , {0, 1, 5, 4} // y-
    , {0, 3, 2, 1} // z-
std::vector<Face> faces;
for (const auto& indices : faceIndices) {
  int i0 = indices[0]; int i1 = indices[1];
  int i2 = indices[2]; int i3 = indices[3];
  Vec  v0 = vertices[i0]; Vec  v1 = vertices[i1];
  Vec  v2 = vertices[i2]; Vec  v3 = vertices[i3];
  Face face(v0, v1, v2, v3);
Check intersections and add periodic particles

We just loop through the particles and check intersections with the three faces. Of course, new particle will need to be assigned IDs. We use the size of the list of initial particles to determine the starting particle ID for the periodic particles. A safer alternative is to find the maximum particle ID in the existing set and then increment that value as new particles are added.

  DEMParticlePArray extraParticles;
  for (const auto& particle : particles) {
    auto position = particle->currentPosition();
    auto axis_a = vcos(particle->currentAnglesAxisA());
    auto axis_b = vcos(particle->currentAnglesAxisB());
    auto axis_c = vcos(particle->currentAnglesAxisC());
    auto radius_a = particle->radiusA();
    auto radius_b = particle->radiusB();
    auto radius_c = particle->radiusC();
    // Create an ellipsoid object for ellipsoid-face intersection tests
    Ellipsoid ellipsoid(position, axis_a, axis_b, axis_c, radius_a, radius_b, radius_c);
       int faceID = 1;
    for (const auto& face : faces) {
      auto status = ellipsoid.intersects(face);
      if (status.first) {
        std::vector<Vec> translations;

        switch (static_cast<Boundary::BoundaryID>(faceID)) {
          case Boundary::BoundaryID::NONE:
          case Boundary::BoundaryID::XMINUS:
            Vec shift(widthX + boundaryMargin, 0, 0);
            addExtraTranslations(shift, boundaryMargin, Vec(0, 1, 1), widthX, widthY, widthZ, face, status, translations);
          case Boundary::BoundaryID::XPLUS:
          case Boundary::BoundaryID::YMINUS:
            Vec shift(0, widthY + boundaryMargin, 0);
            addExtraTranslations(shift, boundaryMargin, Vec(1, 0, 1), widthX, widthY, widthZ, face, status, translations);
          case Boundary::BoundaryID::YPLUS:
          case Boundary::BoundaryID::ZMINUS:
            Vec shift(0, 0, widthZ + boundaryMargin);
            addExtraTranslations(shift, boundaryMargin, Vec(1, 1, 0), widthX, widthY, widthZ, face, status, translations);
          case Boundary::BoundaryID::ZPLUS:
        // Create copies
        for (const auto& translation : translations) {
          DEMParticleP newParticle = std::make_shared<DEMParticle>(*particle);
          newParticle->setCurrentPosition(particle->currentPosition() + translation);
          newParticle->setPreviousPosition(particle->previousPosition() + translation);
      faceID += 2;
Treatment for special cases
vertex numbering
Figure 5. Face orientations assumed for the treatment of special cases.

For special cases where more than one periodic particle has to be inserted for a given particle we use the addExtraTranslations function. Note the special treatment of shared vertices and shared edges.

addExtraTranslations(const Vec& shift, REAL boundaryMargin, Vec inPlaneDiag, REAL widthX, REAL widthY, REAL widthZ, const Face& face, const IntersectionStatus status, std::vector<Vec>& translations)
  // ......

Recall that when we do the ellipsoid-face intersections (see Part 1), we return a pair containing a bool indicating whether there is an intersection and also containing a pair that identifies whether the intersection is at a vertex or an edge along with the vertex or edge index..

Vertex intersection

If the intersection is at a vertex, we check whether the vertex has index 0 (i.e., it is shared by three faces). In that case we create three particles, two parallel to the coordinate axes and the third along a diagonal.

  if (status.second.first == Face::Location::VERTEX) {
    int vertIndex = status.second.second;
    if (vertIndex == 0) {
      for (int ii = 1; ii < 4; ++ii) {
        Vec inPlane = face.vertex[ii] - face.vertex[vertIndex];
        auto length = inPlane.length();
        if (ii == 2) {
          inPlane *= length;
          inPlane += (inPlaneDiag*boundaryMargin);
        } else {
          inPlane *= (length + boundaryMargin);
        Vec outOfPlane = inPlane + shift;

If the vertex has index 1 or 3, it is shared by two faces and a copy along a diagonal is necessary.

    } else if (vertIndex == 1 || vertIndex == 3) {
      Vec normal = Vec(1,1,1) - inPlaneDiag;
      Vec vec1 = normal * Vec(widthX + boundaryMargin,
        widthY + boundaryMargin, widthZ + boundaryMargin);
      Vec vec2 = face.vertex[2] - face.vertex[vertIndex];
        auto length = vec2.length();
      vec2 *= (length + boundaryMargin);
      translations.push_back(vec1 + vec2);
Edge intersection

If the intersection is at an edge, we find an in-plane and an out-of-plane shift.

  else if (status.second.first == Face::Location::EDGE) {
    int edgeIndex = status.second.second;
    if (edgeIndex == 0 || edgeIndex == 3) {
      int oppIndex = (edgeIndex+3) % 4;
      Vec inPlane = face.vertex[oppIndex] - face.vertex[edgeIndex];
      auto length = inPlane.length() + boundaryMargin;
      inPlane *= length;
      Vec outOfPlane = inPlane + shift;
Removing duplicates

At the end of the process, we will have a few duplicate particles which will have to be cleaned up. We store a set of tags in the seen vector and then remove the elements that satisfy the seen condition.

removeDuplicates(DEMParticlePArray& input)
  std::vector<Vec> seen;
  auto newEnd = std::remove_if(input.begin(), input.end(),
    [&seen](const DEMParticleP& particle)
      Vec pos = particle->currentPosition();
      if (std::find_if(seen.begin(), seen.end(), [&pos](const Vec& seen_pos) {return pos == seen_pos;})
           != seen.end()) {
        return true;
      return false;
  input.erase(newEnd, input.end());


The approach discussed here is remarkably fast for RVEs containing a reasonably large number of particles.

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