chr1shr / voro

Voro++: a three-dimensional Voronoi cell library in C++
Other
164 stars 47 forks source link

Voronoi cell from partice or id #35

Closed vlccek closed 1 year ago

vlccek commented 1 year ago

Hello, I need some way to get from the given particle the Voronoi cell. I can use the find_voronoi_cell(...) but this method returns only ID and practice and I can't find how I can get from this the voronoicell object or even voronoicell_neighbors. Does have voro++ some way of doing it?

I need this for passage through neighbors to minimize the number of past Voronoi cells.

chr1shr commented 1 year ago

Here are two ways to do this.

The first method is to pre-compute all of the Voronoi cells and associate them with the ID. Usually if you look at the code examples, a single voronoicell object is re-used to compute all of the Voronoi cells. This is often highly efficient, because it means you can work with a tiny amount of memory: you compute a cell, store what you need, throw it away, and move on. But it's also possible to create an array of voronoicell objects, and keep them all in memory. The array can be indexed by ID, which makes it easy to look anything you need later. This might be a good approach if you may need to re-use the Voronoi cells.

The second method is to make use of the particle_order class. Internally, Voro++ sorts the particles into a grid of blocks for computational efficiency. Every particle is uniquely associated with two integers: its block ID, and its index within the block. When particles are inserted into the container, that information can optionally be stored in the particle_order class. The primary goal of the particle_order class is to later do a full computation in the order that the particles were added. But the particle_order class is essentially just a large array of two-integer index pairs that make it possible to look up where any particular particle went. The compute_cell function that computes a Voronoi cell can take in those two integers to do a computation of a given particle.

The example code below is adapted from the random_points.cc example. It demonstrates the two methods to compute the volume of the Voronoi cell for the particle with ID 5 (or any other particle).

#include "voro++.hh"
using namespace voro;

// Set up the number of blocks that the container is divided into
const int n_x=6,n_y=6,n_z=6;

// Set the number of particles that are going to be randomly introduced
const int particles=20;

// This function returns a random double between 0 and 1
double rnd() {return double(rand())/RAND_MAX;}

int main() {
    int i;
    double x,y,z;

    // Create a container with the geometry given above, and make it
    // non-periodic in each of the three coordinates. Allocate space for
    // eight particles within each computational block
    container con(0,1,0,1,0,1,n_x,n_y,n_z,false,false,false,8);

    // Randomly add particles into the container
    particle_order po;
    for(i=0;i<particles;i++) {
        x=rnd();
        y=rnd();
        z=rnd();
        con.put(po,i,x,y,z);
    }

    // Method 1: pre-compute all of the Voronoi cells in an array
    voronoicell c[particles];
    bool success[particles];
    c_loop_all cl(con);
    if(cl.start()) do {
        int id=cl.pid();
        success[id]=con.compute_cell(c[id],cl);
    } while(cl.inc());

    // Directly look up a Voronoi cell and do some computation
    const int k=5;
    if(success[k]) printf("Volume of cell %d is %g\n",k,c[k].volume());
    else printf("No cell for particle %d\n",k);

    // Method 2: use the particle_order class to look up the position
    // in the container where the particle went
    voronoicell c_single;
    int *ref=po.o+2*k;
    if(con.compute_cell(c_single,*ref,ref[1])) {
        printf("Volume of cell %d is %g\n",k,c_single.volume());
    } else printf("No cell for particle %d\n",k);
}

Note that there is always a chance that a Voronoi cell computation will fail. For the radical Voronoi tessellation, some Voronoi cells may be completely occluded. Thus the example code above catches those cases, by recording when the compute_cell function passed back false, indicating failure.