How to exclude periodic images of the same face from iteration

[PavelStishenko]

Imagine a Voronoi tesselation of a periodic box with a single particle. Let it be a box with boundaries at 0 and 10 over all dimensions. And let particle be placed in the middle at the point (5,5,5).
The code below iterates over the cells (there is only one cell in the considered case) and prints out all faces. It woks well.

My problem with this code is that it prints periodic images of the same faces twice. Two face for each dimension: up and down, left and right, front and back. 6 facets. But I need only one facet for each direction. For example only up, left, and front.

Of course it is possible to find pairs of parallel faces and filter out one of each pair. But I dont know a simple way to do it - the only thing that comes in mind is to find a subspace spanned by vertices of each facet and compare them. It would take SVD or eigendecomposition and seems laborious.

Is it possible to make Voro++ to not produce periodic images of the same faces at all?

#include <iostream>
#include <vector>
#include <voro++/voro++.hh>
using std::vector;
using std::cout;
using std::endl;

using voro::container;
using voro::c_loop_all;
using voro::voronoicell_neighbor;

int main()
{
  container con(0,10,0,10,0,10, 3, 3, 3, true, true, true, 8);
  con.put(1, 5, 5, 5);

  c_loop_all vl(con);
  voronoicell_neighbor c;
  vector<double> vertices;
  vector<int> neighbors, face_vertices, face_orders;

  if(vl.start()) do if(con.compute_cell(c,vl)) {
    c.vertices(5, 5, 5, vertices);
    c.neighbors(neighbors);
    c.face_vertices(face_vertices);
    c.face_orders(face_orders);
    int vidx=0;
    for (int i=0; i<c.number_of_faces(); ++i)
    {
      const int nv = face_orders.at(i); // number of vertices in i-th face
      cout << "Face # " << i << endl;
      for (int j=0; j<nv; ++j)
      {
        int cvi = face_vertices.at(vidx+1+j); // vertex index in all cell vertices vector
        cout << "   " << j << " : " << vertices.at(cvi*3) << " " << vertices.at(cvi*3+1) << " " << vertices.at(cvi*3+2) << endl;
      }
      vidx += face_vertices.at(vidx) + 1;
    }
  }while(vl.inc());
}

[PavelStishenko]

Knowing face normal vector would help.

[chr1shr]

There are some practical ways to handle this, although there are some limitations (described below) that would prevent a completely robust solution.

First, there is indeed a function called normals that returns all of the normal vectors n=(nx,ny,nz) the faces. For the pairs of faces that you describe above, the normals should point in opposite directions. If you picked a random fixed direction vector p=(px,py,pz) then you could compute p.n, and then only keep the face if p.n > 0. In practice that would likely work very well. The only issue with this is that there could be floating point errors in some marginal cases, where you might end up picking both faces or none at all.

Another approach, in general, is to use the neighbor information. If look at the polygons.cc example you will see that it picks out a unique face from a pair by only keeping those where the neighbor ID is larger than the particle's own ID.

There are unfortunately some issues with this in the periodic case. In your example above, when a single particle meets its periodic images, those get tagged -1,...,-6 in the neighbor ID list. But if there are more particles, then it is possible for particle A to have faces from multiple periodic images of particle B. In this case there are repeated entries of B in the neighbor list and currently no way to differentiate them. At some point I plan to add functionality so that the neighbor information also contains the periodic image displacements, but there are some issues with the code structure that mean it's not super straightforward to do.

Hopefully the normals function will work well enough for you case. A hybrid of normal vectors and neighbor information might be highly reliable in practice.

[PavelStishenko]

voronoicell_neighbor::normals() worked like a charm! Thank you! But I am interested, where are those negative tags for periodic images "tagged -1,...,-6 in the neighbor ID list"?

Also, if particle A have faces from multiple periodic images of particle B, these faces are intrinsically different, because a particle can be asymmetric by it's nature. Therefore, for example, one face "sees" left side of particle A and right side of particle B, and the other face "sees" right side of particle A and left side of particle "B". Genuinely same faces can be filtered-out by the condition "the neighbor ID is larger than the particle's own ID".

Thank you very much!

[chr1shr]

The -1, -2, ..., -6 neighbor IDs correspond to faces that are made in the orthogonal directions by the images of the current particle. This is mentioned in the custom output documentation. This provides a way to disambiguate faces due to the periodic image of the current particle. But there's still the issue with disambiguating faces due to periodic images of any other particle, which appear identical in the neighbor information.